球體表面隨機均勻采樣方法shell程序(awk)


要在單位球面上隨機選取一個點,從 [-180,180) 中的均勻分布經度和 [-90,90] 中均勻分布的緯度中隨機選擇經緯度是不正確的,因為赤道和兩極的格點面積元並不相同,因此以這種方式選取的點將在極點附近“聚集”。(圖片引自 https://mathworld.wolfram.com/SpherePointPicking.html )

 

 一個簡單的均勻采樣方法是:三維笛卡爾坐標系中生成隨機點,並投影到球面上,這樣生成的點才會均勻。即生成三個高斯隨機變量,x, y,和z,那么以下變量在三維球面上均勻分布的(Muller 1959, Marsaglia 1972):

 

在shell腳本中寫一個awk函數即可實現,以下代碼實現隨機選取兩個點,輸出它們的經緯度和大圓弧距離,輸出五列單位均為度。

 1 seed=10000
 2 sum=$(echo $seed*36/3.1416/3.1416 | bc)
 3 awk 'BEGIN{srand();
 4     for ( i=0; i<='"$sum"'; i++ ){
 5         xa=2*rand()-1;
 6         ya=2*rand()-1;
 7         za=2*rand()-1;
 8         xb=2*rand()-1;
 9         yb=2*rand()-1;
10         zb=2*rand()-1;
11         norma=sqrt(xa*xa+ya*ya+za*za)
12         normb=sqrt(xb*xb+yb*yb+zb*zb)
13         if (norma<0.00001 || normb<0.00001) continue
14         if (xa*xa+ya*ya+za*za<1.0 && xb*xb+yb*yb+zb*zb<1.0){
15             xa=xa/norma
16             ya=ya/norma
17             za=za/norma
18             xb=xb/normb
19             yb=yb/normb
20             zb=zb/normb
21             dis=sqrt((xb-xa)^2+(yb-ya)^2+(zb-za)^2)
22             lona=atan2(ya,xa)/0.017453
23             lonb=atan2(yb,xb)/0.017453
24             lata=90-atan2(sqrt(xa*xa+ya*ya),za)/0.017453
25             latb=90-atan2(sqrt(xb*xb+yb*yb),zb)/0.017453
26             arc=2*atan2(dis/2,sqrt(1-dis*dis/4))/0.017453
27             printf("%8.2f %7.2f %8.2f %7.2f %8.3f\n",
28                 lona,lata,lonb,latb,arc)
29         }
30     }
31 }' > sphere_random_point.txt

備注:

seed代表最終生成文件大致行數

norma代表A點到球心距離

normb代表B點到球心距離

dis代表AB點之間弦長

arc代表AB點之間大圓弧長(度為單位)

 

References

Marsaglia, G. "Choosing a Point from the Surface of a Sphere." Ann. Math. Stat. 43, 645-646, 1972.

Muller, M. E. "A Note on a Method for Generating Points Uniformly on N-Dimensional Spheres." Comm. Assoc. Comput. Mach. 2, 19-20, Apr. 1959.

 

 
 

,弧

noun: 弧, 弧形, 弧光


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM