目次ページ 前ページ 次ページ
4. 初等幾何学の作図例題
4.12 ボロノイ図の作図
ボロノイ図の双対グラフをドロネー(Delaunay)図と言います。計算の過程でドロネー図のデータも得られますので、この二つを同時に作図することにしました。ランダムな点の座標は乱数で作成しました。点の数はコンピュータの計算時間を幾何級数的に増やしますので、最初は10点くらいで試して見て、段々と数を増やして結果を楽しむとよいでしょう。
空間的に分布する点群についても、その点の勢力圏は隣接する点との境界を面で区切り、一つの点を囲う多面体を求める問題になります。このプログラミングには未だ挑戦していません。
ボロノイ図は、計算幾何学的な問題としての興味の他に、種々の問題に応用する可能性があります。例えば、FEM(有限要素法)の自動メッシュ切りに応用することも研究されています。
 図4.11 乱数で発生した点群データ  図4.12 ボロノイ図  図4.13 ドロネー図 20 REM ------- Compute Voronoi Polygon ----- 30 REM --- First Version in 1987, rewritten in 2007--- 50 DEF2PT P: DEF2CR C: DEF2LN L: DEF2ED E: DEF2BX B 60 CLG: LET BOX=0,0,300,200 : BOX=BOX 70 NP=50 : REM ---- try the program changing NP values 80 NPV=NP+NP+2 : NED=NPV+NP-1 90 DIM P[NP], PV[NPV], ID[NP], IED[NED,4], LV[NP] 100 FOR I=1 TO NP 110 A1=600*(RND(0)-0.5):A2=400*(RND(0)-0.5) 111 LET P[I]=A1,A2 120 P[I]=P[I] :REM --- plot the random NP points 130 NEXT I 140 DPTEXT -100,180, "Random Points " 150 PAUSE 160 REM -- Decide a pair of closest points 170 RMAX=1E5 : IPA=1 : LET XX,YY=P[IPA] 180 FOR I=2 TO NP 190 LET X,Y=P[I] 200 IF ABS(X-XX)>RMAX THEN GOTO 240 210 IF ABS(Y-YY)>RMAX THEN GOTO 240 220 R=DIS(P[IPA],P[I]): IF R>RMAX THEN GOTO 240 230 RMAX=R : IPB=I 240 NEXT I 250 ID[IPA]=1 : ID[IPB]=2 : NG1=1 : NG2=2 260 IED[1,1]=IPA : IED[1,2]=IPB : NES=1 : NEE=1 270 REM --- Main procedure to find the circum centre 280 ITRY=0 : NE=NES 290 IF IED[NE,3]=0 THEN IDR1=3 : IDR2=4 : GOTO 330 300 IF IED[NE,4]=0 THEN IDR1=4 : IDR2=3 : GOTO 330 310 IF ITRY=0 THEN NES=NE 320 GOTO 550 330 ITRY=1 : IPA=IED[NE,1] : NGN=ID[IPA] 340 IF NGN<>NG1 THEN GOTO 550 350 IPB=IED[NE,2] 351 IF IDR1=4 THEN KSIDE=1 ELSE KSIDE=-1 360 GOSUB 880 : IF IPC=0 THEN GOTO 540 |
370 NV=NV+1 : PV[NV]=PVC : IED[NE,IDR1]=NV 371 NGN=ID[IPC] 380 IF NGN < >0 THEN GOTO 400 390 NG2=NG2+1 : ID[IPC]=NG2 : GOTO 520 400 K=0 : IP1=IPC : IP2=IPA 410 GOSUB 1010 420 IF IEDN=0 THEN K=1 430 IP1=IPB : IP2=IPC 440 GOSUB 1010 450 IF IEDN=0 THEN K=K+2 460 IF K=0 THEN GOTO 550 470 IF K=1 THEN GOTO 500 480 IF K=2 THEN GOTO 510 490 IF K=3 THEN GOTO 520 500 NEE=NEE+1 : IED[NEE,1]=IPA 501 IED[NEE,2]=IPC : IED[NEE,IDR2]=NV: GOTO 550 510 NEE=NEE+1 : IED[NEE,1]=IPB : IED[NEE,2]=IPC 511 IED[NEE,IDR1]=NV: GOTO 550 520 NEE=NEE+1 : IED[NEE,1]=IPA : IED[NEE,2]=IPC 521 IED[NEE,IDR2]=NV 530 NEE=NEE+1 : IED[NEE,1]=IPB : IED[NEE,2]=IPC 531 IED[NEE,IDR1]=NV: GOTO 550 540 NLV=NLV+1 : LV[NLV]=L1 : IED[NE,IDR1]=-NLV 550 ILOOP=ILOOP+1 : IF ILOOP=1 THEN GOTO 290 560 NE=NE+1 : IF NE<=NEE THEN GOTO 290 570 NG1=NG1+1: IF NG1=NP THEN GOTO 590 ELSE GOTO 280 580 REM ---------------- Results Display 590 CLG : BOX=BOX : REM ------ Voronoi Polygon 600 FOR I=1 TO NP : P[I]=P[I] : NEXT I 610 FOR I=1 TO NEE 620 I1=IED[I,3] : I2=IED[I,4] 630 IF I1<0 THEN GOTO 690 640 IF I2>0 THEN GOTO 730 650 GROFF 660 I2=ABS(I2) : L1=LV[I2] : L1=REV(L1) : E1=L1 & BOX 670 GRON 680 E1=PV[I1] : GOTO 760 690 GROFF 700 I1=ABS(I1) : L1=LV[I1] : E1=L1 & BOX 710 GRON 720 E1=PV[I2] : GOTO 760 730 I1=ABS(I1) 740 I2=ABS(I2) 750 E1=PV[I1]@PV[I2] 760 NEXT I: 770 DPTEXT -100,180,"Voronoi Polygon" 780 PAUSE : REM -------- Delaunay Polygon 790 CLG : BOX=BOX 800 FOR I=1 TO NEE 810 I1=IED[I,1] : I2=IED[I,2] 820 E1=P[I1]@P[I2] 830 NEXT I 840 DPTEXT -100,180,"Delaunay Figure" 850 END 860 REM ---------Subroutines ------------ 870 REM -- Find the minimum circumcircle through A,B,C 880 I=2 : L1=LBSEC(P[IPA],P[IPB]) : RMAX=1E5 : IPC=0 890 LAB=P[IPA]@P[IPB] : IF KSIDE<0 THEN LAB=REV(LAB) 900 IF I=IPA THEN GOTO 980 910 IF I=IPB THEN GOTO 980 920 NGN=ID[I] : IF NGN=0 THEN GOTO 940 930 IF NGN<NG1 THEN GOTO 980 940 RD=DIS(P[I],LAB) : IF RD<0 THEN GOTO 980 950 PR=L1 & LBSEC(P[I],P[IPA]): IF FLAG < >0 THEN GOTO 980 960 R=DIS(PR,LAB) : IF R>RMAX THEN GOTO 980 970 RMAX=R : IPC=I : PVC=PR 980 I=I+1 : IF I>NP THEN RETURN 990 GOTO 900 1000 REM ------ Search edge table 1010 IEDN=0 : IE=NES 1020 IQ1=IED[IE,1] : IQ2=IED[IE,2] 1030 ID1=IED[IE,IDR1] : ID2=IED[IE,IDR2] 1040 IF ((ID1< >0) & (ID2< >0)) THEN GOTO 1130 1050 IF IP1=IQ1 THEN GOTO 1070 1060 IF IP2=IQ1 THEN GOTO 1100 ELSE GOTO 1130 1070 IF IP2< >IQ2 THEN GOTO 1130 1080 IF ID1< >0 THEN GOTO 1130 1090 IEDN=IE : IED[IEDN,IDR1]=NV : RETURN 1100 IF IP1< >IQ2 THEN GOTO 1130 1110 IF ID2< >0 THEN GOTO 1130 1120 IEDN=IE : IED[IEDN,IDR2]=NV : RETURN 1130 IE=IE+1 : IF IE>NEE THEN RETURN 1140 GOTO 1020 1150 REM ------ end of Voronoi ----- |