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

Question

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

Answer 1

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
      }

Answer 2

Here, I assume that the position of the plane is known and that you want to calculate the normal vector. I do not answer the question how this plane is calculated in the first place.

The normal vector to a plane defined by three points (e.g. A, B, C) in (3D) space can be calculated with the cross product of two vectors a=A-C and b=B-C (in the graph below, point A and B would be at the tips of the arrows a and b, respectively. Point C is at the starting point of the two vectors). I also assume that the points are not in a line.

The resulting vector is perpendicular to a and b, therefore it is normal to the given plane.

To get a vector of length 1, you would have to divide it by its length. Now to the question: In POVRay, I assume you have the coordinates in some variables. Then (omitting any #declare statements) you calculate the normal vector (the inices 1,2,3 correspond to the components in x,y,z direction):

n1 = a2*b3 - a3*b2;

n2 = a3*b1 - a1*b3;

n3 = a1*b2 - a2*b1;

The length L of vector n is L=sqrt(n1*n1 + n2*n2 + n3*n3). Now you divide each component by L and you have the normal vector of unit length.

The direction of the normal vector depends on how the three points are arranged on the plane. Counter-clockwise means that the normal vector is directed upwards (or outwards) from the plane, and vice versa. I used this to check whether a surface is visible or not (i.e. its normal vector points towards the point where the camera is located or away from the camera).