3

2つの変数で記述される不規則なメッシュがあります。各面を構成する頂点のインデックスを格納するfaces配列と、各頂点の座標を格納するverts配列です。また、各面で区分的に一定であると想定される関数があり、面ごとの値の配列の形式で格納されます。

fこのデータから関数を構築する方法を探しています。次の線に沿った何か：

faces = [[0,1,2], [1,2,3], [2,3,4] ...]
verts = [[0,0], [0,1], [1,0], [1,1],....]
vals = [0.0, 1.0, 0.5, 3.0,....]

f = interpolate(faces, verts, vals)

f(0.2, 0.2) = 0.0 # point inside face [0,1,2]
f(0.6, 0.6) = 1.0 # point inside face [1,2,3]

手動で評価f(x,y)する方法は、ポイントが存在する対応する面を見つけて、x,yその面に格納されている値を返すことです。scipy（またはmatlab）でこれをすでに実装している関数はありますか？

4

5 に答える 5

1

ブリッジCGAL-pythonモジュールを使用することをお勧めします。私の記憶が正しければ、CGALは三角形検索の機能を提供します。ただし、組み込みの表現を使用して三角形分割を段階的に構築する場合は、おそらく最も効率的に機能します。簡単で汚いものの場合は、ボロノイ図（Matlabでのこの機能は優れていません）を使用するか、単一のクエリの場合はすべてを計算することで、クエリポイントに最も近いメッシュ頂点を見つけることができます。距離を求め、最小値を見つけて、その頂点を持つすべての三角形を検索します。

于 2010-02-20T02:01:02.353 に答える
1

Matlabには組み込み関数inpolygonがあり、三角形の内側にいるかどうかをテストできます。あなたがどちらの顔をしているのかを特定する機能はわかりません。

このような関数を作成する場合は、最初にポイントが最も近い頂点をテストし、次に、一致するものが見つかるまで、頂点を共有するすべての面のinpolygonを評価します。これはかなり速いはずです。

于 2010-02-20T02:01:14.533 に答える
1

これは、ポイントが内部にある三角形の面を見つけるほど、補間のようには聞こえません。各三角形の面をテストする方法については、このサイトをチェックしてください。関数は、内部にある面を判別し、対応する値を返します。もちろん、顔がたくさんある場合、またはこれをたくさん行う場合は、これを最適化する方法を見つける必要があります（x方向とy方向の両方で各三角形の最も遠い+点と-点を保存し、保存しますたとえば、ポイントがこの境界ボックス内にない場合は、三角形の内側にあるかどうかを確認しない方がよいでしょう）。

Matlabまたはscipyに組み込まれているものが、あなたが望むことを実行するために見つかるとは本当に思えませんが、私は間違っている可能性があります。

于 2010-02-19T23:22:19.880 に答える
1

を見てくださいmatplotlib.delaunay.interpolate。これは、Cコードの十分に文書化されたラッパーです。
（ただしclass LinearInterpolator、「現時点では、通常の長方形グリッドのみが補間にサポートされています」と記載されています。）

于 2010-03-07T20:08:49.927 に答える
1

MATLABには、必要な処理を実行する組み込み関数はありません。Jonasが提案したように、関数INPOLYGONを使用して独自のアルゴリズムを構築することもできますが、ポイントがポリゴン内にあるかどうかを確認するための標準アルゴリズムを使用して、より高速な実装を自分で作成できる場合があります。

しばらく前に、線分と三角形のサーフェスのセットとの交点を3Dで見つけるための独自のコードを作成しました。このソフトサーファーのリンクが、アルゴリズムの実装に最も役立つことがわかりました。私の場合はあなたの場合よりも複雑でした。2次元で作業しているため、セグメントが三角形の平面と交差する点を見つけることに関するリンクの最初のセクションは無視してかまいません。

使用できるように、MATLABコードの簡略化されたバージョンを以下に含めました。この関数interpolateはfaces、、、、verticesおよび行列を入力として受け取り、指定された（x、y）ポイントで評価できるvalues関数ハンドルを返し、境界三角形内の区分的値を取得します。fこのコードのいくつかの機能は次のとおりです。

  • 評価時に実行される処理は、入れ子関数fに含まれています。この関数は内の他の変数にアクセスできるため、三角形内テストに必要ないくつかの変数が事前に計算され、可能な限り高速に実行されます。 evaluate_functioninterpolateevaluate_function
  • 三角形がたくさんある場合、ポイントがすべての三角形の内側にあるかどうかをテストすると、コストがかかる可能性があります。したがって、コードは最初に、ポイントの特定の半径（つまり、三角形の最長の脚の長さ）内にある三角形を見つけます。これらの近くの三角形だけが、ポイントがそれらの中にあるかどうかを確認するためにテストされます。
  • ポイントが三角形の領域内にない場合、NaNの値はによって返されfます。

コードに含まれていないものがいくつかあり、使用目的に応じて追加することができます。

  • 入力チェック：facesコードは現在、がN行3列の行列、verticesM行2列の行列、およびvalues長さNのベクトルであると想定しています。入力がこれらの要件に準拠していることを確認するためにエラーチェックを追加し、それらの1つ以上が正しくない場合にエラーをスローすることをお勧めします。
  • 縮退三角形チェック：facesおよび入力によって定義された1つまたは複数の三角形verticesが縮退している可能性があります（つまり、面積が0である可能性があります）。これは、2つ以上の三角形の頂点が正確に同じ点である場合、または三角形の3つの頂点すべてが直線上にある場合に発生します。評価に関しては、このような三角形を無視するチェックを追加することをお勧めしfます。
  • エッジケースの処理：ポイントが2つ以上の三角形領域のエッジに到達する可能性があります。したがって、ポイントがとる値（つまり、最大の面の値、面の値の平均など）を決定する必要があります。このようなエッジケースの場合、以下のコードは、faces変数内の面のリストの先頭に近い面の値を自動的に選択します。

最後に、コードは次のとおりです。

function f = interpolate(faces,vertices,values)

  %# Precompute some data (helps increase speed):

  triVertex = vertices(faces(:,2),:);              %# Triangles main vertices
  triLegLeft = vertices(faces(:,1),:)-triVertex;   %# Triangles left legs
  triLegRight = vertices(faces(:,3),:)-triVertex;  %# Triangles right legs
  C1 = sum(triLegLeft.*triLegRight,2);  %# Dot product of legs
  C2 = sum(triLegLeft.^2,2);            %# Squared length of left leg
  C3 = sum(triLegRight.^2,2);           %# Squared length of right leg
  triBoundary = max(C2,C3);             %# Squared radius of triangle boundary
  scale = C1.^2-C2.*C3;
  C1 = C1./scale;
  C2 = C2./scale;
  C3 = C3./scale;

  %# Return a function handle to the nested function:

  f = @evaluate_function;

  %# The nested evaluation function:

  function val = evaluate_function(x,y)

    w = [x-triVertex(:,1) y-triVertex(:,2)];
    nearIndex = find(sum(w.^2,2) <= triBoundary);  %# Find nearby triangles
    if isempty(nearIndex)
      val = NaN;           %# Return NaN if no triangles are nearby
      return
    end
    w = w(nearIndex,:);
    wdotu = sum(w.*triLegLeft(nearIndex,:),2);
    wdotv = sum(w.*triLegRight(nearIndex,:),2);
    C = C1(nearIndex);
    s = C.*wdotv-C3(nearIndex).*wdotu;
    t = C.*wdotu-C2(nearIndex).*wdotv;
    inIndex = find((s >= 0) & (t >= 0) & (s+t <= 1),1);
    if isempty(inIndex)
      val = NaN;         %# Return NaN if point is outside all triangles
    else
      val = values(nearIndex(inIndex));
    end

  end

end
于 2010-03-08T16:17:50.493 に答える