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