四角形を形成する 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 として投稿しました。
あなたの問題を交差問題と考えるのが最も簡単だと思います.p0、p1、p2、p3で定義された任意の2D双一次曲面と点pが交差するパラメータ位置(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)
演算子 X を使用して2D 外積を示していることに注意してください(たとえば、p0 X p1 = x0*y1 - y0*x1)。この方程式を二次バーンスタイン多項式としてフォーマットしました。これにより、物事がよりエレガントになり、数値的に安定します。s の解は、この方程式の根です。バーンスタイン多項式の二次公式を使用して根を見つけることができます。
s = ( (A-B) +- sqrt(B^2 - A*C) ) / ( A - 2*B + C )
二次方程式は、+- のために 2 つの答えを与えます。p が双一次曲面内にある解だけに関心がある場合は、s が 0 と 1 の間にない場合の答えを破棄できます。t を見つけるには、t について解いた上記の 2 つの方程式のいずれかに s を代入しますsに関して。
重要な特殊なケースを 1 つ指摘しておきます。分母 A - 2*B + C = 0 の場合、二次多項式は実際には線形です。この場合、より単純な式を使用して s を見つける必要があります。
s = A / (A-C)
これにより、AC = 0 でない限り、厳密に 1 つの解が得られます。A = C の場合、2 つのケースがあります。A=C=0 は、s のすべての値に p が含まれていることを意味し、それ以外の場合、s の値に p が含まれていないことを意味します。
これは、コミュニティ 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;
}
これは私の実装です... 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 つの解がある場合があります。
p が 2 次元の点であれば、簡単に取得できます。その場合、S は T での四角形の全幅の分数成分であり、T は同様に S での四角形の全高の分数成分です。
ただし、p がスカラーの場合、双一次補間関数は必ずしもモノリシックではないため、必ずしも可能ではありません。