;
;[ Infomation ]
; Name : ドロネー三角形分割
; SubName :
; Version :
; copyright :
;
;[ Update history ]
;yyyy/mm/dd : ver : comment
;
;#ifndef __MODULE_NAME__
;#define global __MODULE_NAME__
#module
;---------------------------------------------------------------------------------------------------
;///////////////////////////////////////////////////////////////////////////////////////////////////
;
; 定数の定義
;
;#const global
;#enum global
;///////////////////////////////////////////////////////////////////////////////////////////////////
;
; 外接円の中心座標と半径を取得
;
;[ Infomation ]
; GetCircumcircle double ptx1, double pty1, double ptx2, double pty2, double ptx3, double pty3, var resPtx, var resPty, var resR
; double ptx1, double pty1 : [IN] 1点目のx,y座標
; double ptx2, double pty2 : [IN] 2点目のx,y座標
; double ptx3, double pty3 : [IN] 3点目のx,y座標
; var resPtx, var resPty : [OUT]外接円の中心座標
; var resR : [OUT]外接円の半径
;
; return : 成功判定
;
;[ comment ]
; ３点を通る外接円の中心座標と半径を算出します。
;
; statの値で外接円の作成に成功したかどうかの判定ができます。
;
; stat = 0 :
; 外接円の作成に成功しました。
;
; stat = 1 :
; 外接円の作成に失敗しました。
; ３点が直線上に並んでいたり、２点以上が同じ場所にある場合、外接円は作れません。
; 半径 resR には0.0が返ります。
;
#deffunc GetCircumcircle double ptx1, double pty1, double ptx2, double pty2, double ptx3, double pty3, var resPtx, var resPty, var resR
;resPt 中心座標
;resR 半径
;
;式１：2 * x * (x2-x1) + 2 * y * (y2-y1) = (x2*x2 - x1*x1 + y2*y2 - y1*y1)
;式２：2 * x * (x3-x1) + 2 * y * (y3-y1) = (x3*x3 - x1*x1 + y3*y3 - y1*y1)
; 次のように置き換える
;式１：2*x * a + 2*y * b = c
;式２：2*x * d + 2*y * e = f
;
a = ptx2 - ptx1
d = ptx3 - ptx1
b = pty2 - pty1
e = pty3 - pty1
c = ptx2*ptx2 - ptx1*ptx1 + pty2*pty2 - pty1*pty1
f = ptx3*ptx3 - ptx1*ptx1 + pty3*pty3 - pty1*pty1
;
; 3点が同一直線上にあるかチェック
;
n = 2.0 * ( b*d - a*e)
if n = 0.0 {
; ３点が同一直線上に並ぶ場合は点１をそのまま返す
resPty = ptx1
resPtx = pty1
resR = 0.0
return 1
}
;
; 中心座標
;
resPty = (c*d - a*f) / n
if a ! 0.0 {
resPtx = (0.5*c - b*resPty) / a
} else {
; 直線には並んでいないので、aもdも0.0はありえない。
resPtx = (0.5*f - e*resPty) / d
}
;
; 半径
;
x = resPtx - ptx1
y = resPty - pty1
resR = sqrt( x * x + y * y )
return 0
;///////////////////////////////////////////////////////////////////////////////////////////////////
;
; 指定範囲を包む三角形の座標を取得
;
;[ Infomation ]
; GetHugeTriangle double ltx, double lty, double rux, double ruy, var tx0, var ty0, var tx1, var ty1, var tx2, var ty2
; double ltx : 左上座標X
; double lty : 左上座標Y
; double rux : 右下座標X
; double ruy : 右下座標Y
; tx0, ty0 : 三角形頂点座標 0
; tx1, ty1 : 三角形頂点座標 1
; tx2, ty2 : 三角形頂点座標 2
;
; return : 0
;
;[ comment ]
;指定した矩形領域(ltx, lty)-(rux, ruy)の外接円を内接円とした正三角形の３点の座標を取得します。
;ドロネー三角形を求める際の外部三角形の作成に利用できます。（無駄が多いけど気にしない）
;
#deffunc GetHugeTriangle double ltx, double lty, double rux, double ruy, var tx0, var ty0, var tx1, var ty1, var tx2, var ty2
hw = (rux - ltx)/2.0
hh = (ruy - lty)/2.0
r = sqrt( hw * hw + hh * hh )
r3 = r * sqrt(3)
tx0 = ltx + hw - r3
ty0 = lty + hh - r
tx1 = ltx + hw
ty1 = lty + hh + r*2
tx2 = rux - hw + r3
ty2 = ty0
return
;-------------------------------------------------------------------------------
;---------------------------------------------------------------------------------------------------
#global
;#endif ;__MODULE_NAME__
;#################################################
; モジュール型を作成
;三角形：triangle
;外接円：circumcircle
;
#module modTriangle tri0, tri1, tri2, ccm_x, ccm_y, ccm_r
;------------------------
; コンストラクタ
; モジュール変数に値を登録
#modinit int t0, int t1, int t2, array pts
tri0 = t0
tri1 = t1
tri2 = t2
; 外接円取得
GetCircumcircle pts( 0, t0), pts( 1, t0), pts( 0, t1), pts( 1, t1), pts( 0, t2), pts( 1, t2), xx, yy, rr
ccm_x = xx
ccm_y = yy
ccm_r = rr
return
;------------------------
; デストラクタ
; モジュール変数を開放したときの処理
#modterm
; mes "削除しました。"
return
;------------------------
; 外接円の内側にあるか調べる
;0:外側
;1:内側
#modfunc modTri_IsInsideCircumcircle double ptx, double pty
x = ptx - ccm_x
y = pty - ccm_y
if (x*x + y*y) < ccm_r*ccm_r : return 1 ;内側
return 0 ;外側
;------------------------
; 三角形の頂点を取得
#modfunc modTri_GetPoints var tt0, var tt1, var tt2
; 三角形の頂点番号を返す
tt0 = tri0
tt1 = tri1
tt2 = tri2
return
;------------------------
; 三角形の一致判定
;引数３点と三角形の頂点三点が一致するか知らべる。
;一致する場合１を返す。
#modfunc modTri_IsEquals int t0, int t1, int t2
f = (tri0 = t0) & (tri1 = t1) & (tri2 = t2)
f |= (tri0 = t0) & (tri1 = t2) & (tri2 = t1)
f |= (tri0 = t1) & (tri1 = t0) & (tri2 = t2)
f |= (tri0 = t2) & (tri1 = t0) & (tri2 = t1)
f |= (tri0 = t1) & (tri1 = t2) & (tri2 = t0)
f |= (tri0 = t2) & (tri1 = t1) & (tri2 = t0)
return f
#global
;#################################################
; モジュール
;
#module
;------------------------
; モジュール初期化
#deffunc DelaunayIni
ddim points, 2, 1000
pointsCount = 0 ;有効な点の数
; 点の追加可能領域
;とりあえず固定領域とした
pointAreaX = 0
pointAreaY = 0
pointAreaW = ginfo_winx
pointAreaH = ginfo_winy
; 外部三角形取得
GetHugeTriangle pointAreaX, pointAreaY, pointAreaX + pointAreaW - 1, pointAreaY + pointAreaH - 1, tx0, ty0, tx1, ty1, tx2, ty2
points( 0, 0) = double(tx0), double(ty0)
points( 0, 1) = double(tx1), double(ty1)
points( 0, 2) = double(tx2), double(ty2)
pointsCount = 3
;三角形登録
newmod Triangle, modTriangle, 0,1,2, points
return
;------------------------
; 点の追加
#deffunc DelaunayAddPoint double tx, double ty
; 点リストに追加
points( 0, pointsCount) = tx, ty
pc = pointsCount
pointsCount++
;
; すべての三角形リストから点が外接円に入っているか調べる
; 内側にある場合は、三角形を再分割
;
foreach Triangle
modTri_IsInsideCircumcircle Triangle(cnt), tx, ty
if stat {
; 点が三角形の外接円に入っている場合
; 三角形削除
modTri_GetPoints Triangle(cnt), t0, t1, t2 ; 再分割のため頂点取得
delmod Triangle(cnt)
; 点を中心に三角形に再分割
; 無限ループ回避、重複チェックのため仮領域に一旦登録
newmod TriangleSub, modTriangle, pc, t0, t1, points
newmod TriangleSub, modTriangle, pc, t1, t2, points
newmod TriangleSub, modTriangle, pc, t2, t0, points
}
loop
; 再分割三角形で重複するものは削除
; 比較して一致する場合は両方削除
foreach TriangleSub
modTri_GetPoints TriangleSub(cnt), t0, t1, t2
i = cnt
foreach TriangleSub
; 自分以外と比較する
if i ! cnt {
modTri_IsEquals TriangleSub(cnt), t0, t1, t2
if stat {
delmod TriangleSub(cnt)
delmod TriangleSub(i)
;mes "("+cnt+", "+i+")"+t0+", "+t1+", "+t2
}
}
loop
loop
; 再分割三角形を三角形リストに登録
; 本登録
foreach TriangleSub
modTri_GetPoints TriangleSub(cnt), t0, t1, t2
newmod Triangle, modTriangle, t0, t1, t2, points
; 仮領域削除
delmod TriangleSub(cnt)
loop
return pointsCount
;------------------------
; 描画
#deffunc DelaunayDraw
dim gx, 4
dim gy, 4
dim gc, 4
gc = $FF0000, $00FF00, $0000FF, $0000FF
; 点を描画
r = 2
repeat pointsCount
x = points( 0, cnt)
y = points( 1, cnt)
circle x-r, y-r, x+r, y+r, 1
pos x, y
;mes ""+cnt
loop
; 三角形描画
;境界線は重複してもお構いなしで描画します。
i = 0
foreach Triangle
modTri_GetPoints Triangle(cnt), t0, t1, t2
; 外部三角形につながる三角形は採用しない
if (t0 = 0)|(t0 = 1)|(t0 = 2) | (t1 = 0)|(t1 = 1)|(t1 = 2) | (t2 = 0)|(t2 = 1)|(t2 = 2) : continue
; 面に色付けしてみる
gx = int(points( 0, t0)), int(points( 0, t1)), int(points( 0, t2)), int(points( 0, t2))
gy = int(points( 1, t0)), int(points( 1, t1)), int(points( 1, t2)), int(points( 1, t2))
;gsquare -257, gx, gy, gc
; 辺を描画
line points( 0, t0), points( 1, t0), points( 0, t1), points( 1, t1)
line points( 0, t1), points( 1, t1), points( 0, t2), points( 1, t2)
line points( 0, t2), points( 1, t2), points( 0, t0), points( 1, t0)
i++
loop
pos 0,0
mes "三角形個数："+i
return
#global
;#################################################
;-------------------------------------------------------------------------------
;
; 仮実行スクリプト（デバッグ作業用）
;
;#####################################################################
;ここにはデバッグ作業用のスクリプトを記述します。
;ここを有効にするとこのファイル単独での実行が可能になります。
;
;0 ：リリースモード 本体側から#includeで連結して動作させる場合です。
;1 ：デバッグモード このファイル単品で動作確認が出来ます。
#if 1
; ここにデバッグ用コードを記述する
DelaunayIni
#define POINT_MAX 50
;
; 座標データを登録しながら分割
;
randomize
repeat POINT_MAX
; x = double( rnd( ginfo_winx/2 )+ginfo_winx/4 )
; y = double( rnd( ginfo_winy/2 )+ginfo_winy/4 )
x = double( rnd( ginfo_winx/2 ) )*2.0
y = double( rnd( ginfo_winy/2 ) )*2.0
DelaunayAddPoint x, y
loop
; 作成したドロネー三角形を描画
DelaunayDraw
#endif
;#####################################################################