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
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).
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
}