povray

How to compute normal vector to least square plane in PovRay only


Is it possible to calculate normal vector to a plane defined by set of points using PovRay only (proper set has more than 3 points)? At the moment I'm using external program that computes via Jacobi eigenvalues of a least square plane.

Still it would be nice not to have to switch for this step to different program but just to use internal procedures of PovRay.

Kris


Solution

  • Below Jacobi algorithm in povray format. I hope, it would be useful to somebody. As the example a ring is being drawn in the least-square plane defined by 5 points Pts[0..4]. The coordinates of points need to be initialized in the 3rd line. Of course, if the count is other than 5, you need to change variable NoPts in 2nd line.

    #declare ooo=<0,0,0>;
    #declare NoPts=5; 
    #declare Pts=array[NoPts]{p1,p2,p3,p4,p5}
    
    // compute centroid              
    #declare Centroid=ooo;
    #declare i = 0;
    #while(i < NoPts )
      #declare Centroid = Centroid+Pts[i];
      #declare i = i + 1;
    #end
    #declare Centroid = Centroid / NoPts;
    
    // Move points to centroid
    #declare nPts=array[NoPts]{ooo,ooo,ooo,ooo,ooo}
    #declare i = 0;
    #while(i < NoPts )
      #declare nPts[i] = Pts[i] -  Centroid;
      #declare i = i + 1;
    #end                   
    
    #declare SM=array[3][3]{ {0,0,0}, {0,0,0}, {0,0,0} } // 2nd momentum matrix for coordinates moved to the centroid
    #declare EV=array[3][3]{ {1,0,0}, {0,1,0}, {0,0,1} } // eigen vectors 
    
    #declare i = 0;
    #while(i < NoPts )
      #declare SM[0][0] = nPts[i].x  * nPts[i].x  + SM[0][0];
      #declare SM[0][1] = nPts[i].x  * nPts[i].y  + SM[0][1];
      #declare SM[0][2] = nPts[i].x  * nPts[i].z  + SM[0][2];
      #declare SM[1][1] = nPts[i].y  * nPts[i].y  + SM[1][1];
      #declare SM[1][2] = nPts[i].y  * nPts[i].z  + SM[1][2];
      #declare SM[2][2] = nPts[i].z  * nPts[i].z  + SM[2][2];
      #declare i = i + 1;
    #end                   
    #declare SM[1][0] = SM[0][1];
    #declare SM[2][0] = SM[0][2];
    #declare SM[2][1] = SM[1][2];
    
    // =============== JACOBI rotation                               
    
    #declare tol = 0.00000001; // tolerance
    #declare iterMax = 8;
    #declare summ = 0; // control sum
    
    
    #declare l = 0;            //  ----------- sanity check 
    #while(l < 3 )
      #declare m = 0;
      #while(m < 3 )
          #declare summ = summ + abs(SM[l][m]);   
          #declare m = m + 1;     
      #end    
      #declare l = l + 1;     
    #end
    
    #if (summ < 0.00005)
      #debug concat("=============== STH WRONG \n")       
    #end 
    
    #declare j=0;
    #while(j < iterMax)
       #declare ssum = 0;
       #declare amax = 0;
       #declare row=1;
       #while(row < 3)     // ------- loop over rows
         #declare col=0;
         #while(col < row)     // ------- loop over columns
                  #declare atmp = abs(SM[row][col]);    
                  #if(atmp > amax)
                    #declare amax = atmp;
                  #end
                  #declare ssum = ssum + atmp;
                  #if(atmp < amax * 0.1)
                     #declare col = 5; 
                  //#end 
                  #else
                     #if(abs(ssum / summ) > tol )       //   --- only for "huge" elements ;-)
                        #declare th = atan(2 * SM[row][col] / (SM[col][col] - SM[row][row])) / 2;
                        #declare sin_th = sin(th);
                        #declare cos_th = cos(th); 
                        #declare k=0;
                        #while(k < 3)
                          #declare tmp2 = SM[k][col];                        
                          #declare SM[k][col] =  cos_th * tmp2 + sin_th * SM[k][row];
                          #declare SM[k][row] = -sin_th * tmp2 + cos_th * SM[k][row];
                          #declare tmp2 = EV[k][col];
                          #declare EV[k][col] =  cos_th * tmp2 + sin_th * EV[k][row];
                          #declare EV[k][row] = -sin_th * tmp2 + cos_th * EV[k][row];
                          #declare k = k + 1;
                        #end   
                        #declare SM[col][col] =  cos_th * SM[col][col] + sin_th * SM[row][col];
                        #declare SM[row][row] = -sin_th * SM[col][row] + cos_th * SM[row][row];
                        #declare SM[col][row] = 0;
                        #declare k=0;
                        #while(k < 3)
                          #declare SM[col][k] = SM[k][col];
                          #declare SM[row][k] = SM[k][row];
                          #declare k = k + 1;
                        #end  
                     #end  // end of loop for huge elements
                  #end
            #declare col = col + 1; 
         #end // end col
         #declare row = row + 1; 
      #end  //end row
      #declare j = j + 1;  // ' ------- next iteration if not bigger than iterMAX   
    #end                                     
    // =============== end JACOBI
    
    // find EV with the smallest eigenvalue
    #declare EVmin = 10000; 
    #declare k=0;
    #while (k < 3)
      #if(SM[k][k] < EVmin)
        #declare EVmin = SM[k][k];
        #declare EVindex = k;
      #end
      #declare k=k+1;
    #end   
    
    
    // TEST - draw the ring   
    
    #declare ringRadius=1;
    #declare    w = < EV[0][EVindex], EV[1][EVindex], EV[2][EVindex] >;      
    object{ 
          torus { ringRadius*.75, Dash_Radius texture { Bond_Texture  } }
          Reorient_Trans(<0, 1, 0>,w)
          translate Centroid
          }