是否可以仅使用 PovRay 计算由点集定义的平面的法线向量(正确的点集超过 3 个点)?目前我正在使用通过最小二乘平面的雅可比特征值计算的外部程序。
不必将此步骤切换到不同的程序,而只需使用 PovRay 的内部程序,仍然会很好。
克里斯
是否可以仅使用 PovRay 计算由点集定义的平面的法线向量(正确的点集超过 3 个点)?目前我正在使用通过最小二乘平面的雅可比特征值计算的外部程序。
不必将此步骤切换到不同的程序,而只需使用 PovRay 的内部程序,仍然会很好。
克里斯
在这里,我假设平面的位置是已知的,并且您想要计算法线向量。我不回答首先如何计算这架飞机的问题。
(3D) 空间中由三个点(例如 A、B、C)定义的平面的法线向量可以用两个向量 a=AC 和 b=BC 的叉积来计算(在下图中,点 A 和 B将分别位于箭头 a 和 b 的尖端。点 C 位于两个向量的起点)。我还假设这些点不在一条线上。
结果向量垂直于 a 和 b,因此它垂直于给定平面。
要获得长度为 1 的向量,您必须将其除以它的长度。现在问题来了:在 POVRay 中,我假设您在某些变量中有坐标。然后(省略任何#declare 语句)计算法线向量(inices 1、2、3 对应于 x、y、z 方向的分量):
n1 = a2*b3 - a3*b2;
n2 = a3*b1 - a1*b3;
n3 = a1*b2 - a2*b1;
向量 n 的长度 L 为 L=sqrt(n1*n1 + n2*n2 + n3*n3)。现在你将每个分量除以 L,你就有了单位长度的法线向量。
法向量的方向取决于三个点在平面上的排列方式。逆时针表示法向量从平面向上(或向外),反之亦然。我用它来检查一个表面是否可见(即它的法向量指向相机所在的点或远离相机的点)。
下面是 povray 格式的 Jacobi 算法。我希望,这对某人有用。例如,在由 5 个点 Pts[0..4] 定义的最小二乘平面中绘制一个环。点的坐标需要在第 3 行初始化。当然,如果计数不是 5,则需要在第 2 行更改变量 NoPts。
#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
}