だから、私は今日一日中、率直に言って腹立たしい問題に苦しんでいます。
平面上の三角形の一連の頂点 (3 つの点、6 つの自由パラメーター) が与えられた場合、この三角形と {0,0} および {1,1} で定義される単位正方形との交点の面積を計算する必要があります。(2D の任意の正方形をこれに変換でき、同じ変換で 3 つの頂点を移動できるため、これを選択します)。
そのため、問題は 6 つのパラメーター、3 つのポイントだけに単純化されました...これは、完全なソリューションをコード化して完全なソリューションを見つけるのに十分なほど短いと思います。
(可能であれば、<0.5 秒ごとに文字通り 200 万を超える三角形を GPU で実行したいと考えています。単純化の必要性/データ構造/ライブラリがないため)
解決策への私の試みに関して、私は...私が思いついた方法のリストを手に入れましたが、どれも速くないように見えたり...素敵なケースに固有のものではありません（一般的すぎます）。
オプション 1: 囲まれた多角形を見つけます。三角形から 6 角形まで何でもかまいません。これは、私が見つけた O(n) 時間アルゴリズムで凸多角形の交点を使用して行います。次に、これらの交点 (新しい頂点、最大 7 つの O(n log n) ) を CW または CCw のいずれかの順序で並べ替え、(グリーン関数に基づいて) ポイントで単純な面積アルゴリズムを実行できるようにします ( O(n) 再び）。これは、別の m-gon と交差する任意の凸 n-gon に対して、私ができる最速の方法です。ただし...私の問題はそれほど複雑ではなく、特殊なケースであるため、より良い解決策が必要です...
オプション 2: 私はそれが三角形と単位正方形であることを知っているので、交点のリストをより強引な方法で見つけることができます (率直に言って、上記のように実装するのが少しイライラするアルゴリズムを使用するのではなく)。
チェックポイントは19点のみ。4 点は、三角形の内側にある正方形の角です。3点は四角の中の三角です。次に、三角形の各線について、それぞれが正方形から 4 つの線と交差します (例: y=0、y=1、x=0、x=1 の線)。つまり、さらに 12 ポイントです。したがって、12+3+4 = 19 のチェック ポイントがあります。この交差を行うポイントが最大で 6 つ、最小で 3 つになると、考えられる 2 つの方法のいずれかでフォローアップできます。
2a: x の値を大きくして並べ替え、単純に形状をサブ三角形 / 4 角形に分解します。それぞれの形状は、上端と下端の境界線に基づいた簡単な式を使用します。エリアを合計します。
または 2b: 交点を循環的にソートし、緑関数に基づいて面積を計算します。
残念ながら、私が知る限り、これは依然として複雑です。交点を見つけるために、すべてのケースをもう少し分解し始めることができます。これは、正方形の 0 と 1 だけを知っているため、数学がいくつかの項を脱落させる..しかし、必ずしも単純ではありません。
オプション 3: さまざまな条件に基づいて問題の分離を開始します。例えば。正方形内の三角形の 0、1、2、または 3 点。次に、それぞれのケースについて、考えられるすべての交差数を実行し、多角形のケースごとに、面積の解を一意に書き留めます。
オプション 4: ヘヴィサイド ステップ関数を使用した数式。これはおそらく私が最も望んでいるものです。少し... 大きくなると思いますが、それが可能であり、公式が得られれば計算上最速の実行時間になると楽観的かもしれません。
---全体として、高レベルのライブラリ（クリッパーなど）を使用して解決できることはわかっています。また、さまざまな種類のデータ構造 (リンクされたリストの後に並べ替え) を使用する場合、一般的なソリューションを記述することはそれほど難しくないことも認識しています。そして、これを数回行うだけでよい場合は、これらすべてのケースで問題ありません。しかし、画像処理ステップとして実行する必要があるため、画像ごとに9 * 1024 * 1024回以上実行する必要があり、画像を撮影しています.. 1 fpsとしましょう（技術的には、この速度を押し上げたいと思います可能な限り速くアップしますが、900 万個のこれらの三角形の交差領域の問題を計算するための下限は 1 秒です)。これはCPUでは不可能かもしれませんが、それは問題ありません。とにかくCudaで実装することになるでしょうが、この問題の速度の限界を押し上げたいと思っています。
編集：したがって、オプション2bを使用することになりました。可能な交点は 19 個だけであり、そのうち最大 6 個が形状を定義するため、最初にそれらの 3 ～ 6 個の頂点を見つけます。次に、循環 (CCW) 順序で並べ替えます。そして、その多角形の面積を計算して面積を求めます。
そのために書いたテスト コードを次に示します (Igor 用ですが、疑似コードとして読めるはずです)。 、より良い並べ替えを書くためのオーバーヘッドはそれほど多くありません）...その並べ替え以外では、これ以上速くすることはできないと思います。ただし、このオプションを選択する際に私が持っていた可能性のある提案や見落としを受け入れる.
function calculateAreaUnitSquare(xPos, yPos)
wave xPos
wave yPos
// First, make array of destination. Only 7 possible results at most for this geometry.
Make/o/N=(7) outputVertexX = NaN
Make/o/N=(7) outputVertexY = NaN
variable pointsfound = 0
// Check 4 corners of square
// Do this by checking each corner against the parameterized plane described by basis vectors p2-p0 and p1-p0.
// (eg. project onto point - p0 onto p2-p0 and onto p1-p0. Using appropriate parameterization scaling (not unit).
// Once we have the parameterizations, then it's possible to check if it is inside the triangle, by checking that u and v are bounded by u>0, v>0 1-u-v > 0
variable denom = yPos[0]*xPos[1]-xPos[0]*yPos[1]-yPos[0]*xPos[2]+yPos[1]*xPos[2]+xPos[0]*yPos[2]-xPos[1]*yPos[2]
//variable u00 = yPos[0]*xPos[1]-xPos[0]*yPos[1]-yPos[0]*Xx+yPos[1]*Xx+xPos[0]*Yx-xPos[1]*Yx
//variable v00 = -yPos[2]*Xx+yPos[0]*(Xx-xPos[2])+xPos[0]*(yPos[2]-Yx)+yPos[2]*Yx
variable u00 = (yPos[0]*xPos[1]-xPos[0]*yPos[1])/denom
variable v00 = (yPos[0]*(-xPos[2])+xPos[0]*(yPos[2]))/denom
variable u01 =(yPos[0]*xPos[1]-xPos[0]*yPos[1]+xPos[0]-xPos[1])/denom
variable v01 =(yPos[0]*(-xPos[2])+xPos[0]*(yPos[2]-1)+xPos[2])/denom
variable u11 = (yPos[0]*xPos[1]-xPos[0]*yPos[1]-yPos[0]+yPos[1]+xPos[0]-xPos[1])/denom
variable v11 = (-yPos[2]+yPos[0]*(1-xPos[2])+xPos[0]*(yPos[2]-1)+xPos[2])/denom
variable u10 = (yPos[0]*xPos[1]-xPos[0]*yPos[1]-yPos[0]+yPos[1])/denom
variable v10 = (-yPos[2]+yPos[0]*(1-xPos[2])+xPos[0]*(yPos[2]))/denom
if(u00 >= 0 && v00 >=0 && (1-u00-v00) >=0)
outputVertexX[pointsfound] = 0
outputVertexY[pointsfound] = 0
pointsfound+=1
endif
if(u01 >= 0 && v01 >=0 && (1-u01-v01) >=0)
outputVertexX[pointsfound] = 0
outputVertexY[pointsfound] = 1
pointsfound+=1
endif
if(u10 >= 0 && v10 >=0 && (1-u10-v10) >=0)
outputVertexX[pointsfound] = 1
outputVertexY[pointsfound] = 0
pointsfound+=1
endif
if(u11 >= 0 && v11 >=0 && (1-u11-v11) >=0)
outputVertexX[pointsfound] = 1
outputVertexY[pointsfound] = 1
pointsfound+=1
endif
// Check 3 points for triangle. This is easy, just see if its bounded in the unit square. if it is, add it.
variable i = 0
for(i=0; i<3; i+=1)
if(xPos[i] >= 0 && xPos[i] <= 1 )
if(yPos[i] >=0 && yPos[i] <=1)
if(!((xPos[i] == 0 || xPos[i] == 1) && (yPos[i] == 0 || yPos[i] == 1) ))
outputVertexX[pointsfound] = xPos[i]
outputVertexY[pointsfound] = yPos[i]
pointsfound+=1
endif
endif
endif
endfor
// Check intersections.
// Procedure is: loop over 3 lines of triangle.
// For each line
// Check if vertical
// If not vertical, find y intercept with x=0 and x=1 lines.
// if y intercept is between 0 and 1, then add the point
// Check if horizontal
// if not horizontal, find x intercept with y=0 and y=1 lines
// if x intercept is between 0 and 1, then add the point
for(i=0; i<3; i+=1)
variable iN = mod(i+1,3)
if(xPos[i] != xPos[iN])
variable tx0 = xPos[i]/(xPos[i] - xPos[iN])
variable tx1 = (xPos[i]-1)/(xPos[i] - xPos[iN])
if(tx0 >0 && tx0 < 1)
variable yInt = (yPos[iN]-yPos[i])*tx0+yPos[i]
if(yInt > 0 && yInt <1)
outputVertexX[pointsfound] = 0
outputVertexY[pointsfound] = yInt
pointsfound+=1
endif
endif
if(tx1 >0 && tx1 < 1)
yInt = (yPos[iN]-yPos[i])*tx1+yPos[i]
if(yInt > 0 && yInt <1)
outputVertexX[pointsfound] = 1
outputVertexY[pointsfound] = yInt
pointsfound+=1
endif
endif
endif
if(yPos[i] != yPos[iN])
variable ty0 = yPos[i]/(yPos[i] - yPos[iN])
variable ty1 = (yPos[i]-1)/(yPos[i] - yPos[iN])
if(ty0 >0 && ty0 < 1)
variable xInt = (xPos[iN]-xPos[i])*ty0+xPos[i]
if(xInt > 0 && xInt <1)
outputVertexX[pointsfound] = xInt
outputVertexY[pointsfound] = 0
pointsfound+=1
endif
endif
if(ty1 >0 && ty1 < 1)
xInt = (xPos[iN]-xPos[i])*ty1+xPos[i]
if(xInt > 0 && xInt <1)
outputVertexX[pointsfound] = xInt
outputVertexY[pointsfound] = 1
pointsfound+=1
endif
endif
endif
endfor
// Now we have all 6 verticies that we need. Next step: find the lowest y point of the verticies
// if there are multiple with same low y point, find lowest X of these.
// swap this vertex to be first vertex.
variable lowY = 1
variable lowX = 1
variable m = 0;
for (i=0; i<pointsfound ; i+=1)
if (outputVertexY[i] < lowY)
m=i
lowY = outputVertexY[i]
lowX = outputVertexX[i]
elseif(outputVertexY[i] == lowY)
if(outputVertexX[i] < lowX)
m=i
lowY = outputVertexY[i]
lowX = outputVertexX[i]
endif
endif
endfor
outputVertexX[m] = outputVertexX[0]
outputVertexY[m] = outputVertexY[0]
outputVertexX[0] = lowX
outputVertexY[0] = lowY
// now we have the bottom left corner point, (bottom prefered).
// calculate the cos(theta) of unit x hat vector to the other verticies
make/o/N=(pointsfound) angles = (p!=0)?( (outputVertexX[p]-lowX) / sqrt( (outputVertexX[p]-lowX)^2+(outputVertexY[p]-lowY)^2) ) : 0
// Now sort the remaining verticies based on this angle offset. This will orient the points for a convex polygon in its maximal size / ccw orientation
// (This sort is crappy, but there will be in theory, at most 25 swaps. Which in the grand sceme of operations, isn't so bad.
variable j
for(i=1; i<pointsfound; i+=1)
for(j=i+1; j<pointsfound; j+=1)
if( angles[j] > angles[i] )
variable tempX = outputVertexX[j]
variable tempY = outputVertexY[j]
outputVertexX[j] = outputVertexX[i]
outputVertexY[j] =outputVertexY[i]
outputVertexX[i] = tempX
outputVertexY[i] = tempY
variable tempA = angles[j]
angles[j] = angles[i]
angles[i] = tempA
endif
endfor
endfor
// Now the list is ordered!
// now calculate the area given a list of CCW oriented points on a convex polygon.
// has a simple and easy math formula : http://www.mathwords.com/a/area_convex_polygon.htm
variable totA = 0
for(i = 0; i<pointsfound; i+=1)
totA += outputVertexX[i]*outputVertexY[mod(i+1,pointsfound)] - outputVertexY[i]*outputVertexX[mod(i+1,pointsfound)]
endfor
totA /= 2
return totA
終わり
