gts: 三角形網内の任意の点の高さを求めるには

三角形内のXY座標を指定して高さを得ることは簡単にできます。
これを使って、三角形網内の任意の点の高さを求めることができるのですが、
問題は指定したXY座標に対応する三角形を検索するところです。

空間上のデータ検索ではk-d木がよく使われると思われますが、
面積を持ち、不定形の三角形を検索するのはなかなか難しいようです。
そのためか、gtsライブラリにはどうやらそういった関数は無いようです（たぶん）。

代わりかどうかはわかりませんが、gtsライブラリにはBBTree（BoundingBoxTree）というのがあって、
これを使ってある程度高速に検索はできるようですが、
これは求めるべき唯一つの図形を返すわけではなく、
”確実にこの中にある”というリストを返してきます。
なので、最終的にはこちらで候補リストから対象を検索する必要があります。

使用法

まずはBBTreeを作成します。
幸い、GtsSurfaceオブジェクトからBBTreeを作成する関数があるので、
これを使って以下のようにBBTreeを作成します。

  GNode *tree = gts_bb_tree_surface( surface );

これでBBTreeはできたので、次は指定した点に対応する三角形を
BBTreeから探し出します。

まず、指定した点のGtsPointオブジェクトを作成します。

  GtsPoint *pnt = gts_point_new
    ( gts_point_class(), dX, dY, -1000.0 );

Z値に-1000を指定していますが、これはこの後使います。

次に、BBTreeから”確実にこの中にある”リストを取得します。

  GSList *close = gts_bb_tree_point_closest_bboxes( tree, pnt );

変数closeに候補リストが格納されています。
リストの各要素はGtsBBoxオブジェクトですが、
GtsBBoxオブジェクトは、
そのBoundingBoxの元になっている三角形オブジェクトへのポインタも持っています。
ちなみに、筆者が実行したときの候補リストはこんな感じでした。

候補リスト抽出結果

ここからはこちらで該当する三角形をリストから探し出します。

該当する三角形を探すには、gts_triangle_is_stabbed()
を使います。
これは、指定した点と、その点のZ座標が+無限大という線分が三角形を
貫通したらそのオブジェクトを返すというものです。
そのため、上記のGtsPointオブジェクトを作成したときに
-1000を指定した次第です。

三角形が特定されたら、その三角形から点の高さを求めるには
gts_triangle_interpolate_height()を使います。
この関数は、引数に指定したGtsPointオブジェクトにZ値を入れて返します。

というわけで、一応、以下のようなコードで該当する点の高さを求めることができます。

  do
  {
    GtsBBox *d = GTS_BBOX( close->data );
    GtsTriangle *tr = GTS_TRIANGLE( d->bounded );
    if ( gts_triangle_is_stabbed( tr, pnt, NULL ) )
    {
      gts_triangle_interpolate_height( tr, pnt );
      printf( "height = %.3f", pnt->z );
    }

    close = close->next;
  }
  while ( close );

まだ特訓中なので、もしかしたらもう少しいい方法があるかもしれません。