如何在PovRay中计算最小二乘平面的法向量

Question

是否可以仅使用PovRay计算由一组点定义的平面的法向量(真集合具有3个以上的点)？目前我使用的是外部程序，它通过最小二乘平面的Jacobi特征值来计算。

不过，如果这一步不用切换到不同的程序，而只需使用PovRay的内部过程，那就更好了。

克里斯

Answer 1

以下是povray格式的Jacobi算法。我希望，这会对某些人有用。例如，在由5个点Pts0..4定义的最小二乘平面上绘制一个环。点的坐标需要在第三行初始化。当然，如果计数不是5，则需要更改第二行中的变量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
      }

Answer 2

在这里，我假设平面的位置是已知的，并且您想要计算法向量。我没有回答这个问题，这个平面是如何计算的。

(3D)空间中由三个点(例如A，B，C)定义的平面的法向量可以用两个向量a=A-C和b=B-C的叉积来计算(在下图中，点A和B将分别位于箭头a和b的顶端。点C位于两个向量的起始点)。我还假设这些点不在一条线上。

得到的向量垂直于a和b，因此它垂直于给定的平面。

要得到长度为1的向量，必须将其除以其长度。现在来看问题:在POVRay中，我假设你在一些变量中有坐标。然后(省略任何#declare语句)计算法向量( 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，你就得到了单位长度的法向量。

法向量的方向取决于这三个点在平面上的排列方式。逆时针方向意味着法向量从平面向上(或向外)方向，反之亦然。我用它来检查曲面是否可见(例如，它的法向量指向相机所在的点或远离相机)。