ドロネー図の三角形をリストとして取得する
前提オブジェクト
点のリスト originalPointList
ドロネー図 graph1 = DelaunayTriangulation[originalPointList]
課題
ドロネー図の三角形を、リストオブジェクトとして取得したい。
上図の例では、下図のように、5つの三角形からなるリストを取得したいということである。
取得手順
パスパラメータの割り当てを分析する
上図のように、ドロネー図のパスパラメータの割り当ては、良く分からない。
しかし、一定の法則性を見出すことができる。
法則①：パスパラメータを0から等速で増やしていくと、点が線分上を通過していく挙動、線分以外の良く分からない場所を動いていく挙動、ある点上で停止する挙動が次々と観察される。
これら観察されるそれぞれの挙動を、「フェーズ」と総称することにする。点が線分上を通過していく挙動を「通過フェーズ」、線分以外の良く分からない場所を通過していく挙動を「ジャンプフェーズ」、ある点上で停止する挙動を「停止フェーズ」と呼ぶことにする。
法則②：フェーズの切り替わりは、等間隔に訪れる。すなわち、フェーズの始まりにおけるパスパラメータと、フェーズの終わりにおけるパスパラメータとの差は、常に（どのフェーズでも）一定である。
法則③：どの通過フェーズおよびジャンプフェーズにおいても、パスパラメータを等速で増加させると、点は等速で動く。
法則④：パスパラメータを0から増やしていくと、最初に訪れるフェーズは通過フェーズである。そのフェーズでは、パスパラメータ0の点が始点、パスパラメータ1の点が終点である。
以下では、これらの法則を利用していく。
「パスパラメータ0の点と1の点との中点」のパスパラメータを2倍すると、1フェーズに割り当てられているパスパラメータの量を計算できる（unitParameter）。
unitParameter = 2PathParameter[ClosestPoint[graph1, Midpoint[Point[graph1, 0], Point[graph1, 1]]]]
unitParameterの逆数をとれば、フェーズの数を計算できる。
frequency = 1 / unitParameter
各フェーズの中間における点のリストを取得する
half = unitParameter / 2
halfParamList = Sequence[Point[graph1, half + k unitParameter], k, 0, frequency - 1]
halfParamListは、各フェーズの中間における点のリストである（下図の赤い点）。
各フェーズの1/4段階における点のリストを取得する
quarter = unitParameter / 4
quarterParamList = Sequence[Point[graph1, quarter + k unitParameter], k, 0, frequency - 1]
quarterParamListは、各フェーズの1/4段階における点のリストである（下図の緑色の点）。
各フェーズの1/4段階〜中間を線分で結ぶ
quarterToHalfSegmentListFull = Unique[Zip[Segment[s, t], s, quarterParamList, t, halfParamList]]
quarterToHalfSegmentListFullは、各フェーズの1/4段階〜中間を線分で結んだものである（下図の赤い線分）。
quarterToHalfSegmentListRemoved = Unique[RemoveUndefined[Zip[If[Point[s, 1] ≟ ClosestPoint[graph1, Point[s, 1]], s, ?], s, quarterToHalfSegmentListFull]]]
quarterToHalfSegmentListFullのうち、ドロネー図上にないもの（ジャンプフェーズにおけるもの）を削除する（quarterToHalfSegmentListRemoved、下図）
Dilateで拡大して、ドロネー図の線分と揃える
dilateSegmentListFull = Unique[Zip[Dilate[s, 4, Point[s, 1 / 3]], s, quarterToHalfSegmentListRemoved]]
quarterToHalfSegmentListRemovedをDilateで適切に拡大して、ドロネー図の線分と長さ・位置を揃える（dilateSegmentListFull、下図のオレンジの線分）。
delauneySegmentList = RemoveUndefined[Zip[If[Length[s] ≠ 0, s, ?], s, dilateSegmentListFull]]
dilateSegmentListFullのうち、長さ0のもの（停止フェーズがあると生じる）を削除する（delauneySegmentList、下図の水色の線分）。
delauneySegmentListは、本稿の課題の最終目標ではないが、これ自体も、ドロネー図の各線分をリストで返すものとして、意義があるだろう。
点のリスト originalPointListによって作成しうる全ての三角形をリストで返す
polygonListFull = Unique[Flatten[Sequence[Sequence[Sequence[Polygon[Element[originalPointList, s], Element[originalPointList, t], Element[originalPointList, u]], s, 1, Length[originalPointList]], t, 1, Length[originalPointList]], u, 1, Length[originalPointList]]]]
polygonListRemoved = RemoveUndefined[Zip[If[Area[s] ≠ 0, s, ?], s, polygonListFull]]
必要な三角形のみを抽出する
delauneyPolygonList = RemoveUndefined[Zip[If[Sides[s] ⊂ delauneySegmentList, s, ?], s, polygonListRemoved]]
polygonListRemovedの要素である各三角形のうち、3辺がすべてdelauneySegmentList（ドロネー図の線分のリスト）の要素であるような三角形だけを抽出する（delauneyPolygonList、下図の青い三角形）。これが、目的のリストである。
なお、Sides[<多角形>]は、オリジナルツールである。詳細は、下記記事参照。
多角形オブジェクトを用いて、頂点または辺のリストを返す - うしブログ
参考情報
以下の数式は、点のリストoriginalPointListから、ドロネー図の三角形のリストを返す。
処理落ち回避のため、CopyFreeObjectを入れてある。
Zip[Polygon[τ],τ,CopyFreeObject[Zip[{Vertex[σ]}, σ, RemoveUndefined[Zip[If[Sides[ρ] ⊂ RemoveUndefined[Zip[If[Length[η] ≠ 0, η, ?], η, Unique[Zip[Dilate[ζ, 4, Point[ζ, 1 / 3]], ζ, Unique[RemoveUndefined[Zip[If[Point[ε, 1] ≟ ClosestPoint[DelaunayTriangulation[originalPointList], Point[ε, 1]], ε, ?], ε, Unique[Zip[Segment[γ, δ], γ, Sequence[Point[DelaunayTriangulation[originalPointList], (2PathParameter[ClosestPoint[DelaunayTriangulation[originalPointList], Midpoint[Point[DelaunayTriangulation[originalPointList], 0], Point[DelaunayTriangulation[originalPointList], 1]]]]) / 4 + β (2PathParameter[ClosestPoint[DelaunayTriangulation[originalPointList], Midpoint[Point[DelaunayTriangulation[originalPointList], 0], Point[DelaunayTriangulation[originalPointList], 1]]]])], β, 0, 1 / (2PathParameter[ClosestPoint[DelaunayTriangulation[originalPointList], Midpoint[Point[DelaunayTriangulation[originalPointList], 0], Point[DelaunayTriangulation[originalPointList], 1]]]]) - 1], δ, Sequence[Point[DelaunayTriangulation[originalPointList], (2PathParameter[ClosestPoint[DelaunayTriangulation[originalPointList], Midpoint[Point[DelaunayTriangulation[originalPointList], 0], Point[DelaunayTriangulation[originalPointList], 1]]]]) / 2 + α (2PathParameter[ClosestPoint[DelaunayTriangulation[originalPointList], Midpoint[Point[DelaunayTriangulation[originalPointList], 0], Point[DelaunayTriangulation[originalPointList], 1]]]])], α, 0, 1 / (2PathParameter[ClosestPoint[DelaunayTriangulation[originalPointList], Midpoint[Point[DelaunayTriangulation[originalPointList], 0], Point[DelaunayTriangulation[originalPointList], 1]]]]) - 1]]]]]]]]]], ρ, ?], ρ, RemoveUndefined[Zip[If[Area[ξ] ≠ 0, ξ, ?], ξ, Unique[Flatten[Sequence[Sequence[Sequence[Polygon[Element[originalPointList, κ], Element[originalPointList, λ], Element[originalPointList, μ]], κ, 1, Length[originalPointList]], λ, 1, Length[originalPointList]], μ, 1, Length[originalPointList]]]]]]]]]]]
