https://stackoverflow.com/questions/808441
-
03-07-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
Russian質問
4つの2dポイント、p0 =（x0、y0）、p1 =（x1、y1）などがあり、四辺形を形成しています。私の場合、四角形は長方形ではありませんが、少なくとも凸である必要があります。
p2 --- p3
| |
t | p |
| |
p0 --- p1
s
バイリニア補間を使用しています。 SとTは[0..1]内にあり、補間点は次のようになります。
bilerp(s,t) = t*(s*p3+(1-s)*p2) + (1-t)*(s*p1+(1-s)*p0)
ここに問題があります。クワッド内にあることがわかっている2dポイントpがあります。バイリニア補間を使用するときにそのポイントを与えるs、tを見つけたいです。
双一次補間を逆にする簡単な式はありますか？
ソリューションに感謝します。 Naaffのソリューションの実装をwikiとして投稿しました。
解決
問題を交差点の問題と考えるのが最も簡単だと思います：ポイントpがp0、p1、p2、p3で定義される任意の2D双線形面と交差するパラメーター位置（s、t）は何ですか
この問題を解決するためのアプローチは、tspauldの提案に似ています。
xとyに関する2つの方程式から始めます：
x = (1-s)*( (1-t)*x0 + t*x2 ) + s*( (1-t)*x1 + t*x3 )
y = (1-s)*( (1-t)*y0 + t*y2 ) + s*( (1-t)*y1 + t*y3 )
tの解決：
t = ( (1-s)*(x0-x) + s*(x1-x) ) / ( (1-s)*(x0-x2) + s*(x1-x3) )
t = ( (1-s)*(y0-y) + s*(y1-y) ) / ( (1-s)*(y0-y2) + s*(y1-y3) )
これら2つの方程式を互いに等しく設定して、tを排除することができます。すべてを左側に移動して単純化すると、次の形式の方程式が得られます。
A*(1-s)^2 + B*2s(1-s) + C*s^2 = 0
場所：
A = (p0-p) X (p0-p2)
B = ( (p0-p) X (p1-p3) + (p1-p) X (p0-p2) ) / 2
C = (p1-p) X (p1-p3)
2Dクロスを示すために演算子Xを使用したことに注意してください。製品（例：p0 X p1 = x0 * y1-y0 * x1）。この式を2次Bernstein多項式としてフォーマットしました物事はよりエレガントで、数値的に安定しています。 sの解はこの方程式の根です。バーンスタイン多項式の2次式を使用して根を見つけることができます。
s = ( (A-B) +- sqrt(B^2 - A*C) ) / ( A - 2*B + C )
2次式は、+-による2つの答えを与えます。 pが双一次曲面内にある解のみに関心がある場合は、sが0から1の間にない回答を破棄できます。tを見つけるには、tを解いた上記の2つの方程式のいずれかにsを代入しますsに関して。
重要な特別なケースを1つ指摘する必要があります。分母A-2 * B + C = 0の場合、2次多項式は実際には線形です。この場合、より簡単な方程式を使用してsを見つける必要があります。
s = A / (A-C)
これは、AC = 0でない限り、正確に1つの解を与えます。A= Cの場合、2つのケースがあります。A= C = 0は、sの all 値がpを含むことを意味します。 にはsの値にはpが含まれません。
他のヒント
Newtonの方法 を繰り返し使用できます次の非線形方程式系を解きます。
p = p0*(1-s)*(1-t) + p1*s*(1-t) + p2*s*t + p3*(1-s)*t.
2つの方程式（方程式のx成分とy成分の両方に等しい）と2つの未知数（sとt）があることに注意してください。
ニュートンの方法を適用するには、残余 rが必要です。 / p>
r = p - (p0*(1-s)*(1-t) + p1*s*(1-t) + p2*s*t + p3*(1-s)*t),
およびヤコビ行列J。残差を微分することで求められます。私たちの問題では、ヤコビアンは次のとおりです。
J(:,1) = dr/ds = -p0*(1-t) + p1*(1-t) + p2*t - p3*t (first column of matrix)
J(:,2) = dr/dt = -p0*(1-s) - p1*s + p2*s + p3*(1-s) (second column).
Newtonの方法を使用するには、最初の推測（s、t）から始めて、次の反復を数回実行します。
(s,t) = (s,t) - J^-1 r,
sとtの新しい値を使用して、各反復でJとrを再計算します。各反復で、Jを係数行列、rを右辺とする2x2線形システムを解くことにより、ヤコビアンの逆数を残差（J ^ -1 r）に適用するのが支配的なコストです。
メソッドの直観：
直感的に、四辺形が平行四辺形である場合、問題を解決するのははるかに簡単です。ニュートンの方法は、連続した平行四辺形近似で四辺形の問題を解いています。各反復で
-
点（s、t）で局所微分情報を使用して、四辺形を平行四辺形で近似します。
-
線形システムを解くことにより、平行四辺形近似の下で正しい（s、t）値を見つけます。
-
この新しいポイントにジャンプして繰り返します。
メソッドの利点：
ニュートン型の方法で予想されるように、収束は非常に高速です。反復が進むと、メソッドが真のポイントに近づくだけでなく、局所的な平行四辺形近似もより正確になるため、収束速度自体が増加します（反復メソッドの用語では、Newtonのメソッドは 2次収束 ）。実際には、これは正しい桁数が反復ごとに約2倍になることを意味します。
これは、私が行ったランダム試行の反復対エラーの代表的な表です（コードについては以下を参照）：
Iteration Error
1 0.0610
2 9.8914e-04
3 2.6872e-07
4 1.9810e-14
5 5.5511e-17 (machine epsilon)
2回の反復後、エラーは事実上気付かないほど小さく、ほとんどの実用的な目的に十分であり、5回の反復後、結果は機械精度。
反復回数を修正する場合（たとえば、ほとんどの実用的なアプリケーションでは3回の反復、非常に高い精度が必要な場合は8回の反復）、アルゴリズムには非常にシンプルで単純なロジックがあり、構造がうまく機能します高性能計算へ。あらゆる種類の特別なエッジケース*をチェックしたり、結果に応じて異なるロジックを使用したりする必要はありません。これは、いくつかの単純な式を含むforループにすぎません。以下に、このアプローチがインターネットやその他の回答で見られる従来のフォーミュラベースのアプローチに勝るいくつかの利点を強調します。
-
コーディングが簡単。 forループを作成し、いくつかの式を入力するだけです。
-
条件または分岐なし（if / then）。これにより、通常、パイプラインの効率。
-
平方根はなく、反復ごとに1除算しか必要ありません（適切に記述されている場合）。他のすべての操作は単純な追加で、subtrac
コミュニティWikiとしてのNaaffのソリューションの実装です。どうもありがとう。
これはC実装ですが、c ++で動作するはずです。ファズテスト機能が含まれています。
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
int equals( double a, double b, double tolerance )
{
return ( a == b ) ||
( ( a <= ( b + tolerance ) ) &&
( a >= ( b - tolerance ) ) );
}
double cross2( double x0, double y0, double x1, double y1 )
{
return x0*y1 - y0*x1;
}
int in_range( double val, double range_min, double range_max, double tol )
{
return ((val+tol) >= range_min) && ((val-tol) <= range_max);
}
/* Returns number of solutions found. If there is one valid solution, it will be put in s and t */
int inverseBilerp( double x0, double y0, double x1, double y1, double x2, double y2, double x3, double y3, double x, double y, double* sout, double* tout, double* s2out, double* t2out )
{
int t_valid, t2_valid;
double a = cross2( x0-x, y0-y, x0-x2, y0-y2 );
double b1 = cross2( x0-x, y0-y, x1-x3, y1-y3 );
double b2 = cross2( x1-x, y1-y, x0-x2, y0-y2 );
double c = cross2( x1-x, y1-y, x1-x3, y1-y3 );
double b = 0.5 * (b1 + b2);
double s, s2, t, t2;
double am2bpc = a-2*b+c;
/* this is how many valid s values we have */
int num_valid_s = 0;
if ( equals( am2bpc, 0, 1e-10 ) )
{
if ( equals( a-c, 0, 1e-10 ) )
{
/* Looks like the input is a line */
/* You could set s=0.5 and solve for t if you wanted to */
return 0;
}
s = a / (a-c);
if ( in_range( s, 0, 1, 1e-10 ) )
num_valid_s = 1;
}
else
{
double sqrtbsqmac = sqrt( b*b - a*c );
s = ((a-b) - sqrtbsqmac) / am2bpc;
s2 = ((a-b) + sqrtbsqmac) / am2bpc;
num_valid_s = 0;
if ( in_range( s, 0, 1, 1e-10 ) )
{
num_valid_s++;
if ( in_range( s2, 0, 1, 1e-10 ) )
num_valid_s++;
}
else
{
if ( in_range( s2, 0, 1, 1e-10 ) )
{
num_valid_s++;
s = s2;
}
}
}
if ( num_valid_s == 0 )
return 0;
t_valid = 0;
if ( num_valid_s >= 1 )
{
double tdenom_x = (1-s)*(x0-x2) + s*(x1-x3);
double tdenom_y = (1-s)*(y0-y2) + s*(y1-y3);
t_valid = 1;
if ( equals( tdenom_x, 0, 1e-10 ) && equals( tdenom_y, 0, 1e-10 ) )
{
t_valid = 0;
}
else
{
/* Choose the more robust denominator */
if ( fabs( tdenom_x ) > fabs( tdenom_y ) )
{
t = ( (1-s)*(x0-x) + s*(x1-x) ) / ( tdenom_x );
}
else
{
t = ( (1-s)*(y0-y) + s*(y1-y) ) / ( tdenom_y );
}
if ( !in_range( t, 0, 1, 1e-10 ) )
t_valid = 0;
}
}
/* Same thing for s2 and t2 */
t2_valid = 0;
if ( num_valid_s == 2 )
{
double tdenom_x = (1-s2)*(x0-x2) + s2*(x1-x3);
double tdenom_y = (1-s2)*(y0-y2) + s2*(y1-y3);
t2_valid = 1;
if ( equals( tdenom_x, 0, 1e-10 ) && equals( tdenom_y, 0, 1e-10 ) )
{
t2_valid = 0;
}
else
{
/* Choose the more robust denominator */
if ( fabs( tdenom_x ) > fabs( tdenom_y ) )
{
t2 = ( (1-s2)*(x0-x) + s2*(x1-x) ) / ( tdenom_x );
}
else
{
t2 = ( (1-s2)*(y0-y) + s2*(y1-y) ) / ( tdenom_y );
}
if ( !in_range( t2, 0, 1, 1e-10 ) )
t2_valid = 0;
}
}
/* Final cleanup */
if ( t2_valid && !t_valid )
{
s = s2;
t = t2;
t_valid = t2_valid;
t2_valid = 0;
}
/* Output */
if ( t_valid )
{
*sout = s;
*tout = t;
}
if ( t2_valid )
{
*s2out = s2;
*t2out = t2;
}
return t_valid + t2_valid;
}
void bilerp( double x0, double y0, double x1, double y1, double x2, double y2, double x3, double y3, double s, double t, double* x, double* y )
{
*x = t*(s*x3+(1-s)*x2) + (1-t)*(s*x1+(1-s)*x0);
*y = t*(s*y3+(1-s)*y2) + (1-t)*(s*y1+(1-s)*y0);
}
double randrange( double range_min, double range_max )
{
double range_width = range_max - range_min;
double rand01 = (rand() / (double)RAND_MAX);
return (rand01 * range_width) + range_min;
}
/* Returns number of failed trials */
int fuzzTestInvBilerp( int num_trials )
{
int num_failed = 0;
double x0, y0, x1, y1, x2, y2, x3, y3, x, y, s, t, s2, t2, orig_s, orig_t;
int num_st;
int itrial;
for ( itrial = 0; itrial < num_trials; itrial++ )
{
int failed = 0;
/* Get random positions for the corners of the quad */
x0 = randrange( -10, 10 );
y0 = randrange( -10, 10 );
x1 = randrange( -10, 10 );
y1 = randrange( -10, 10 );
x2 = randrange( -10, 10 );
y2 = randrange( -10, 10 );
x3 = randrange( -10, 10 );
y3 = randrange( -10, 10 );
/*x0 = 0, y0 = 0, x1 = 1, y1 = 0, x2 = 0, y2 = 1, x3 = 1, y3 = 1;*/
/* Get random s and t */
s = randrange( 0, 1 );
t = randrange( 0, 1 );
orig_s = s;
orig_t = t;
/* bilerp to get x and y */
bilerp( x0, y0, x1, y1, x2, y2, x3, y3, s, t, &x, &y );
/* invert */
num_st = inverseBilerp( x0, y0, x1, y1, x2, y2, x3, y3, x, y, &s, &t, &s2, &t2 );
if ( num_st == 0 )
{
failed = 1;
}
else if ( num_st == 1 )
{
if ( !(equals( orig_s, s, 1e-5 ) && equals( orig_t, t, 1e-5 )) )
failed = 1;
}
else if ( num_st == 2 )
{
if ( !((equals( orig_s, s , 1e-5 ) && equals( orig_t, t , 1e-5 )) ||
(equals( orig_s, s2, 1e-5 ) && equals( orig_t, t2, 1e-5 )) ) )
failed = 1;
}
if ( failed )
{
num_failed++;
printf("Failed trial %d\n", itrial);
}
}
return num_failed;
}
int main( int argc, char** argv )
{
int num_failed;
srand( 0 );
num_failed = fuzzTestInvBilerp( 100000000 );
printf("%d of the tests failed\n", num_failed);
getc(stdin);
return 0;
}
2Dで作業しているので、 bilerp 関数は実際には2つの方程式です。1つはx、1はyです。それらは次の形式で書き直すことができます。
x = t * s * A.x + t * B.x + s * C.x + D.x
y = t * s * A.y + t * B.y + s * C.y + D.y
場所：
A = p3 - p2 - p1 + p0
B = p2 - p0
C = p1 - p0
D = p0
最初の方程式を書き換えて、 s の観点から t を取得し、2番目の方程式に置き換えて、 s を解きます。
これは私の実装です... tfinigaよりも高速だと思います
void invbilerp( float x, float y, float x00, float x01, float x10, float x11, float y00, float y01, float y10, float y11, float [] uv ){
// substition 1 ( see. derivation )
float dx0 = x01 - x00;
float dx1 = x11 - x10;
float dy0 = y01 - y00;
float dy1 = y11 - y10;
// substitution 2 ( see. derivation )
float x00x = x00 - x;
float xd = x10 - x00;
float dxd = dx1 - dx0;
float y00y = y00 - y;
float yd = y10 - y00;
float dyd = dy1 - dy0;
// solution of quadratic equations
float c = x00x*yd - y00y*xd;
float b = dx0*yd + dyd*x00x - dy0*xd - dxd*y00y;
float a = dx0*dyd - dy0*dxd;
float D2 = b*b - 4*a*c;
float D = sqrt( D2 );
float u = (-b - D)/(2*a);
// backsubstitution of "u" to obtain "v"
float v;
float denom_x = xd + u*dxd;
float denom_y = yd + u*dyd;
if( abs(denom_x)>abs(denom_y) ){ v = -( x00x + u*dx0 )/denom_x; }else{ v = -( y00y + u*dy0 )/denom_y; }
uv[0]=u;
uv[1]=v;
/*
// do you really need second solution ?
u = (-b + D)/(2*a);
denom_x = xd + u*dxd;
denom_y = yd + u*dyd;
if( abs(denom_x)>abs(denom_y) ){ v = -( x00x + u*dx0 )/denom_x; }else{ v2 = -( y00y + u*dy0 )/denom_y; }
uv[2]=u;
uv[3]=v;
*/
}
および派生
// starting from bilinear interpolation
(1-v)*( (1-u)*x00 + u*x01 ) + v*( (1-u)*x10 + u*x11 ) - x
(1-v)*( (1-u)*y00 + u*y01 ) + v*( (1-u)*y10 + u*y11 ) - y
substition 1:
dx0 = x01 - x00
dx1 = x11 - x10
dy0 = y01 - y00
dy1 = y11 - y10
we get:
(1-v) * ( x00 + u*dx0 ) + v * ( x10 + u*dx1 ) - x = 0
(1-v) * ( y00 + u*dy0 ) + v * ( y10 + u*dy1 ) - y = 0
we are trying to extract "v" out
x00 + u*dx0 + v*( x10 - x00 + u*( dx1 - dx0 ) ) - x = 0
y00 + u*dy0 + v*( y10 - y00 + u*( dy1 - dy0 ) ) - y = 0
substition 2:
x00x = x00 - x
xd = x10 - x00
dxd = dx1 - dx0
y00y = y00 - y
yd = y10 - y00
dyd = dy1 - dy0
// much nicer
x00x + u*dx0 + v*( xd + u*dxd ) = 0
y00x + u*dy0 + v*( yd + u*dyd ) = 0
// this equations for "v" are used for back substition
v = -( x00x + u*dx0 ) / ( xd + u*dxd )
v = -( y00x + u*dy0 ) / ( yd + u*dyd )
// but for now, we eliminate "v" to get one eqution for "u"
( x00x + u*dx0 ) / ( xd + u*dxd ) = ( y00y + u*dy0 ) / ( yd + u*dyd )
put denominators to other side
( x00x + u*dx0 ) * ( yd + u*dyd ) = ( y00y + u*dy0 ) * ( xd + u*dxd )
x00x*yd + u*( dx0*yd + dyd*x00x ) + u^2* dx0*dyd = y00y*xd + u*( dy0*xd + dxd*y00y ) + u^2* dy0*dxd
// which is quadratic equation with these coefficients
c = x00x*yd - y00y*xd
b = dx0*yd + dyd*x00x - dy0*xd - dxd*y00y
a = dx0*dyd - dy0*dxd
pが単一の値であり、pが正方形の四隅の最小値と最大値の間にある場合、いいえ、単一の解（s、t ）双線形内挿がその値を与えるように。
一般に、正方形の内側には無限の数の解（s、t）があります。それらは、正方形を通る曲線（双曲線）のパスに沿って配置されます。
関数がベクトル値の場合、正方形の未知のポイントに2つの既知の値がありますか？正方形の各コーナーにある2つのパラメーターの既知の値を考えると、解が存在する場合がありますが、それを保証するものではありません。これは、2つの独立した独立した問題として見ることができることを忘れないでください。それぞれの解は、正方形を通る双曲線の等高線に沿って配置されます。輪郭のペアが正方形内で交差する場合、解が存在します。それらが交差しない場合、ソリューションは存在しません。
また、問題を解決するための簡単な式が存在するかどうかを尋ねます。申し訳ありませんが、実際にはそうではありません。私が言ったように、曲線は双曲線です。
1つの解決策は、別の補間方法に切り替えることです。したがって、双線形ではなく、正方形を三角形のペアに分割します。各三角形内で、真の線形補間を使用できるようになりました。これで、各三角形内の2つの未知数における2つの方程式の線形システムを解くことができます。対応する区分的線形等高線が偶然に一致するまれな縮退ケースを除いて、各三角形に1つのソリューションがあります。
一部の回答では、質問が若干誤解されています。すなわち。 （s、t）座標を検索するクワッド内の補間位置p（x、y）ではなく、未知の補間関数の値が与えられていると仮定しています。これはより単純な問題であり、クワッドを通る2つの直線の交点である解決策が保証されています。
軸の場合と同様に、線の1つはセグメントp0p1とp2p3をカットし、もう1つはp0p2とp1p3をカットします。これらの線はp（x、y）の位置によって一意に定義され、この点で明らかに交差します。
p0p1とp2p3をカットするラインだけを考えると、そのようなラインのファミリーを想像できます。選択した異なるs値ごとに、それぞれが異なる幅でクワッドをカットします。 s値を修正する場合、t = 0およびt = 1を設定することで2つのエンドポイントを見つけることができます。
最初に、行の形式が次のとおりであると仮定します。 y = a0 * x + b0
特定のs値を修正すると、この行の2つのエンドポイントがわかります。それらは：
（1-s）p0 +（s）p1
（1-s）p2 +（s）p3
これらの2つのエンドポイントを考えると、線の方程式にそれらを接続し、s0の関数としてa0とb0を解くことにより、線のファミリーを決定できます。 S値を設定すると、特定の行の式が得られます。ここで必要なのは、このファミリのどの行がポイントp（x、y）にヒットするかを把握することです。 p（x、y）の座標を直線式に差し込むだけで、sのターゲット値を解くことができます。
対応するアプローチを実行して、tを見つけることもできます。
まあ、pが2Dポイントであれば、はい、簡単に取得できます。その場合、SはTでの四辺形の全幅の小数成分であり、Tは同様にSでの四辺形の全高の小数成分です。
ただし、pがスカラーの場合、双線形補間関数は必ずしもモノリシックではないため、必ずしも可能ではありません。
