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

Question

如標題所示，這是一個PovRay問題。

我想知道，如果有可能使用PovRay（適當的集合超過3點）計算法線向量由一組點定義的平面。目前我正在使用通過最小二乘平面的Jacobi特徵值計算的外部程序。

儘管如此，不必爲了這個步驟切換到不同的程序，而只是使用PovRay的內部程序。

克里斯

Answer 1

下面是以povray格式的雅可比算法。我希望，這對某個人有用。例如，一個環正在由5個點Pts [0..4]定義的最小平方平面中繪製。點的座標需要在第3行中初始化。當然，如果計數不是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 = AC和b = BC的叉積來計算（在下面的圖中，點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，並且具有單位長度的法向量。

法向量的方向取決於三點如何排列在平面上。逆時針意味着法向矢量從平面向上（或向外），反之亦然。我用它來檢查表面是否可見（即其法線矢量指向相機所在的位置或遠離相機的位置）。