如何仅在 PovRay 中计算法向量到最小二乘平面答案

【问题标题】：How to compute normal vector to least square plane in PovRay only如何仅在 PovRay 中计算法向量到最小二乘平面
【发布时间】：2012-05-28 08:36:35
【问题描述】：

是否可以仅使用 PovRay 计算由点集定义的平面的法线向量（正确的点集超过 3 个点）？目前我正在使用通过最小二乘平面的雅可比特征值计算的外部程序。

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

克里斯

【问题讨论】：

标签： povray

【解决方案1】：

在这里，我假设平面的位置是已知的，并且您要计算法线向量。我不回答首先如何计算这架飞机的问题。

（3D）空间中由三个点（例如 A、B、C）定义的平面的法线向量可以用两个向量 a=A-C 和 b=B-C 的叉积计算（在下图中，点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，你就得到了单位长度的法线向量。

法向量的方向取决于三个点在平面上的排列方式。逆时针表示法向量从平面向上（或向外），反之亦然。我用它来检查一个表面是否可见（即它的法向量指向相机所在的点或远离相机的点）。

【讨论】：

很抱歉这么长时间没有评论（最近工作太多）。谢谢你花时间写所有这些东西。不幸的是，你假设我知道这架飞机是错误的。正如我所写的，我想首先计算最小二乘平面，并且我只有一组点（如上所述，超过 3 个点；其中只有 3 个点是精确定义的，我们根本不需要任何回归方法）。所以让我换一种说法 - 如何仅在 povray 中对角化 [3x3] 矩阵？

【解决方案2】：

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

【讨论】：