Inverse Bilineare Interpolation?
https://stackoverflow.com/questions/808441
-
03-07-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianFrage
Ich habe vier 2d-Punkte p0 = (x0,y0), p1 = (x1,y1), etc.dass die form eines Vierecks.In meinem Fall, das quad ist nicht rechteckig, aber es sollte zumindest konvex sein.
p2 --- p3
| |
t | p |
| |
p0 --- p1
s
Ich bin mit der bilinearen interpolation.S und T sind innerhalb von [0..1] und der interpolierte Punkt ist gegeben durch:
bilerp(s,t) = t*(s*p3+(1-s)*p2) + (1-t)*(s*p1+(1-s)*p0)
Hier ist das problem..Ich habe ein 2d-point-p, ich weiß, ist im inneren des quad.Ich möchte zu finden die s,t, mir diesen Punkt bei Verwendung der bilinearen interpolation.
Gibt es eine einfache Formel zur Umkehrung der bilinearen interpolation?
Vielen Dank für die Lösungen.Ich habe meine Implementierung von Naaff die Lösung, wie ein wiki.
Lösung
Ich denke, es ist am einfachsten des Problems als ein Schnitt Problem zu denken. Was die Parameter Position ist (s, t), wobei der Punkt p die beliebige 2D bilineare Oberfläche von p0, p1, p2 und p3
Der Ansatz, den ich zur Lösung dieses Problems nehmen werde, ist ähnlich wie tspauld Vorschlag.
Starten Sie mit zwei Gleichungen in Bezug auf x und y:
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 )
Solving für 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) )
Wir können nun diese beiden Gleichungen gleich zueinander gesetzt t zu eliminieren. Umzug alles auf die linke Seite und vereinfachen wir eine Gleichung der Form erhalten:
A*(1-s)^2 + B*2s(1-s) + C*s^2 = 0
Dabei gilt:
A = (p0-p) X (p0-p2)
B = ( (p0-p) X (p1-p3) + (p1-p) X (p0-p2) ) / 2
C = (p1-p) X (p1-p3)
Beachten Sie, dass ich den Operator X verwendet haben die 2D-Quer zu bezeichnen Produkt (zB p0 X p1 = x0 * y1 - y0 * x1). Ich habe diese Gleichung als quadratische Bernsteinpolynom wie dies macht Dinge elegant und numerisch stabil ist. Die Lösungen s sind die Wurzeln dieser Gleichung. Wir können die Wurzeln mit der quadratischen Formel für Bernstein-Polynome finden:
s = ( (A-B) +- sqrt(B^2 - A*C) ) / ( A - 2*B + C )
Die quadratische Formel gibt zwei Antworten aufgrund der + -. Wenn Sie daran interessiert sind nur in Lösungen sind, wo p innerhalb der bilinearen Oberfläche liegt, dann können Sie eine Antwort verwerfen, wenn s nicht zwischen 0 und 1 ist t zu finden, einfach ersetzen s zurück in eine der beiden Gleichungen, die oben in dem wir für t gelöst in Bezug auf die s.
Ich soll einen wichtigen Spezialfall hinweisen. Wenn der Nenner A - 2 * B + C = 0 dann quadratisches Polynom tatsächlich linear ist. In diesem Fall müssen Sie eine viel einfachere Gleichung zu finden s verwenden:
s = A / (A-C)
Dies wird Ihnen genau eine Lösung, es sei denn, AC = 0. Ist A = C, dann haben Sie zwei Fälle: A = C = 0 bedeutet alle Werte für s p enthalten, sonst no Werte für s enthalten p.
Andere Tipps
Könnte man verwenden Newton ' s Methode um iterativ lösen Sie die folgenden nichtlinearen Gleichungssystems:
p = p0*(1-s)*(1-t) + p1*s*(1-t) + p2*s*t + p3*(1-s)*t.
Beachten Sie, dass es sind zwei Gleichungen (Gleichheit in der x-und y-Komponenten der Gleichung), und zwei unbekannten (s und t).
Gelten Newtons Methode müssen wir die Restwert r, ist:
r = p - (p0*(1-s)*(1-t) + p1*s*(1-t) + p2*s*t + p3*(1-s)*t),
und die Die Jacobi-matrix J, die sich durch die Differenzierung der Rest.Für unser problem der Jacobi-ist:
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).
Zu Newton ' s Methode startet man mit einer Schätzung (s,t) und führt dann folgende iteration eine Handvoll von Zeiten:
(s,t) = (s,t) - J^-1 r,
Neuberechnung J und r jede iteration mit der neuen Werte von s und t.Bei jeder iteration die dominierenden Kosten in die Anwendung der inversen der Jacobi-um den Rest - (J^-1 r), durch lösen eines 2x2 lineares system mit J als die koeffizientenmatrix und r der rechten Seite.
Intuition für die Methode:
Intuitiv, wenn das viereck der waren Parallelogramm, es wäre viel einfacher, das problem zu lösen.Newton ' s Methode ist die Lösung des Vierecks problem Parallelogramm mit aufeinander folgenden Näherungswerte.Bei jeder iteration werden wir
Verwenden Sie die lokale abgeleitete Informationen auf den Punkt (s,t) zur Angleichung der viereck mit einem Parallelogramm.
Finden Sie die richtige (s,t) - Werte unter der Parallelogramm-approximation, bei der Lösung eines linearen Systems.
Sprung zu diesem neuen Punkt und wiederholen Sie.
Vorteile der Methode:
Wie erwartet, Newton-Typ-Verfahren, die Konvergenz ist extrem schnell.Als die Iterationen gehen, nicht nur, dass die Methode näher und näher an den wahren Punkt, aber die lokalen Parallelogramm Annäherungen werden mehr genaue als gut, so dass die rate der Konvergenz selbst erhöht (in der iterative Methoden jargon, wir sagen, dass Newtons Methode ist quadratisch konvergente).In der Praxis bedeutet dies, dass die Anzahl der korrekten Ziffern in etwa verdoppelt sich mit jeder iteration.
Hier ist eine repräsentative Tabelle der Iterationen vs.Fehler an einem zufälligen Versuch habe ich (siehe unten für den code):
Iteration Error
1 0.0610
2 9.8914e-04
3 2.6872e-07
4 1.9810e-14
5 5.5511e-17 (machine epsilon)
Nach zwei Iterationen der Fehler ist klein genug, um effektiv unnoticable und gut genug für die meisten praktischen Zwecke, und nach 5 Iterationen das Ergebnis ist genau auf die Grenzen der Maschine Präzision.
Wenn Sie Update die Anzahl der Iterationen (sagen wir, 3 Iterationen für die meisten praktischen Anwendungen und 8 Iterationen, wenn Sie brauchen, sehr hohe Genauigkeit), dann ist der Algorithmus hat eine sehr einfache und einfache Logik, mit einer Struktur, die sich besonders gut zu high-performance computation.Es gibt keine Notwendigkeit für Kontrolle, alle Arten von edge-Fällen*, und mit unterschiedlichen Logik-abhängig von den Ergebnissen.Es ist nur eine for-Schleife mit ein paar einfachen Formeln.Im folgenden werde ich markieren mehrere Vorteile, die dieser Ansatz über die traditionellen Rezeptur basierte Ansätze gefunden, die in anderen Antworten hier, und rund um das internet:
Einfach code.Machen Sie einfach eine for-Schleife und geben Sie ein paar Formeln.
Keine Bedingungen oder Verzweigung (wenn/dann), die in der Regel viel besser pipelining die Effizienz.
Keine Wurzeln, und nur 1 Abteilung Aufwand pro iteration (wenn gut geschrieben).Alle anderen Operationen sind einfache Ergänzungen, Subtraktionen und Multiplikationen.Quadratwurzeln und Divisionen sind in der Regel mehrere Male langsamer als Ergänzungen/Multiplikationen/Multiplikationen, und kann Schraube up cache Effizienz auf bestimmten Architekturen (insbesondere auf bestimmte eingebettete Systeme).In der Tat, wenn Sie schauen Sie unter der Haube auf, wie Quadratwurzel und Divisionen tatsächlich sind berechnet durch die modernen Programmiersprachen, die Sie verwenden beide Varianten des Newton-Verfahrens, manchmal in hardware und manchmal in der software abhängig von der Architektur.
Kann leicht vektorisiert arbeiten auf arrays mit riesigen Anzahl von Quadraten auf einmal.Siehe meine vektorisierter code unten ein Beispiel, wie Sie dies tun.
Erstreckt sich auf beliebige Dimensionen.Der Algorithmus, erstreckt sich in einer einfachen Art und Weise zu inverse multilineares interpolation in beliebigen Anzahl von Dimensionen (2d, 3d, 4d, ...).Ich habe eine 3D-version unter, und man kann sich vorstellen, schreiben eine einfache rekursive version, oder verwenden die automatische Differenzierung Bibliotheken zu gehen, um n-Dimensionen.Newton ' s Methode in der Regel Exponate unabhängig von der dimension der Konvergenz-raten, also im Prinzip die Methode sollte skalierbar sein, um ein paar tausend Abmessungen (!) auf aktuelle hardware (nach dem Punkt, den Bau und die Lösung des n-by-n-matrix J wird wahrscheinlich der limitierende Faktor).
Natürlich, die meisten dieser Dinge hängt auch von der hardware, compiler, und eine Vielzahl von anderen Faktoren, so kann Ihre Laufleistung variieren.
Code:
Sowieso, hier ist mein Matlab-code:(Ich lasse hier alles in die public domain)
Grundlegende 2D-version:
function q = bilinearInverse(p,p1,p2,p3,p4,iter)
%Computes the inverse of the bilinear map from [0,1]^2 to the convex
% quadrilateral defined by the ordered points p1 -> p2 -> p3 -> p4 -> p1.
%Uses Newton's method. Inputs must be column vectors.
q = [0.5; 0.5]; %initial guess
for k=1:iter
s = q(1);
t = q(2);
r = p1*(1-s)*(1-t) + p2*s*(1-t) + p3*s*t + p4*(1-s)*t - p;%residual
Js = -p1*(1-t) + p2*(1-t) + p3*t - p4*t; %dr/ds
Jt = -p1*(1-s) - p2*s + p3*s + p4*(1-s); %dr/dt
J = [Js,Jt];
q = q - J\r;
q = max(min(q,1),0);
end
end
Beispiel Verwendung:
% Test_bilinearInverse.m
p1=[0.1;-0.1];
p2=[2.2;-0.9];
p3=[1.75;2.3];
p4=[-1.2;1.1];
q0 = rand(2,1);
s0 = q0(1);
t0 = q0(2);
p = p1*(1-s0)*(1-t0) + p2*s0*(1-t0) + p3*s0*t0 + p4*(1-s0)*t0;
iter=5;
q = bilinearInverse(p,p1,p2,p3,p4,iter);
err = norm(q0-q);
disp(['Error after ',num2str(iter), ' iterations: ', num2str(err)])
Beispiel-Ausgabe:
>> test_bilinearInverse
Error after 5 iterations: 1.5701e-16
Schnell Vektorgrafik-2D-version:
function [ss,tt] = bilinearInverseFast(px,py, p1x,p1y, p2x,p2y, p3x,p3y, p4x,p4y, iter)
%Computes the inverse of the bilinear map from [0,1]^2 to the convex
% quadrilateral defined by the ordered points p1 -> p2 -> p3 -> p4 -> p1,
% where the p1x is the x-coordinate of p1, p1y is the y-coordinate, etc.
% Vectorized: if you have a lot of quadrilaterals and
% points to interpolate, then p1x(k) is the x-coordinate of point p1 on the
% k'th quadrilateral, and so forth.
%Uses Newton's method. Inputs must be column vectors.
ss = 0.5 * ones(length(px),1);
tt = 0.5 * ones(length(py),1);
for k=1:iter
r1 = p1x.*(1-ss).*(1-tt) + p2x.*ss.*(1-tt) + p3x.*ss.*tt + p4x.*(1-ss).*tt - px;%residual
r2 = p1y.*(1-ss).*(1-tt) + p2y.*ss.*(1-tt) + p3y.*ss.*tt + p4y.*(1-ss).*tt - py;%residual
J11 = -p1x.*(1-tt) + p2x.*(1-tt) + p3x.*tt - p4x.*tt; %dr/ds
J21 = -p1y.*(1-tt) + p2y.*(1-tt) + p3y.*tt - p4y.*tt; %dr/ds
J12 = -p1x.*(1-ss) - p2x.*ss + p3x.*ss + p4x.*(1-ss); %dr/dt
J22 = -p1y.*(1-ss) - p2y.*ss + p3y.*ss + p4y.*(1-ss); %dr/dt
inv_detJ = 1./(J11.*J22 - J12.*J21);
ss = ss - inv_detJ.*(J22.*r1 - J12.*r2);
tt = tt - inv_detJ.*(-J21.*r1 + J11.*r2);
ss = min(max(ss, 0),1);
tt = min(max(tt, 0),1);
end
end
Für die Geschwindigkeit dieser code verwendet implizit die folgende Formel für die inverse einer 2x2 matrix:
[a,b;c,d]^-1 = (1/(ad-bc))[d, -b; -c, a]
Beispiel Verwendung:
% test_bilinearInverseFast.m
n_quads = 1e6; % 1 million quads
iter = 8;
% Make random quadrilaterals, ensuring points are ordered convex-ly
n_randpts = 4;
pp_xx = zeros(n_randpts,n_quads);
pp_yy = zeros(n_randpts,n_quads);
disp('Generating convex point ordering (may take some time).')
for k=1:n_quads
while true
p_xx = randn(4,1);
p_yy = randn(4,1);
conv_inds = convhull(p_xx, p_yy);
if length(conv_inds) == 5
break
end
end
pp_xx(1:4,k) = p_xx(conv_inds(1:end-1));
pp_yy(1:4,k) = p_yy(conv_inds(1:end-1));
end
pp1x = pp_xx(1,:);
pp1y = pp_yy(1,:);
pp2x = pp_xx(2,:);
pp2y = pp_yy(2,:);
pp3x = pp_xx(3,:);
pp3y = pp_yy(3,:);
pp4x = pp_xx(4,:);
pp4y = pp_yy(4,:);
% Make random interior points
ss0 = rand(1,n_quads);
tt0 = rand(1,n_quads);
ppx = pp1x.*(1-ss0).*(1-tt0) + pp2x.*ss0.*(1-tt0) + pp3x.*ss0.*tt0 + pp4x.*(1-ss0).*tt0;
ppy = pp1y.*(1-ss0).*(1-tt0) + pp2y.*ss0.*(1-tt0) + pp3y.*ss0.*tt0 + pp4y.*(1-ss0).*tt0;
pp = [ppx; ppy];
% Run fast inverse bilinear interpolation code:
disp('Running inverse bilinear interpolation.')
tic
[ss,tt] = bilinearInverseFast(ppx,ppy, pp1x,pp1y, pp2x,pp2y, pp3x,pp3y, pp4x,pp4y, 10);
time_elapsed = toc;
disp(['Number of quadrilaterals: ', num2str(n_quads)])
disp(['Inverse bilinear interpolation took: ', num2str(time_elapsed), ' seconds'])
err = norm([ss0;tt0] - [ss;tt],'fro')/norm([ss0;tt0],'fro');
disp(['Error: ', num2str(err)])
Beispiel-Ausgabe:
>> test_bilinearInverseFast
Generating convex point ordering (may take some time).
Running inverse bilinear interpolation.
Number of quadrilaterals: 1000000
Inverse bilinear interpolation took: 0.5274 seconds
Error: 8.6881e-16
3D-version:
Enthält einige code, um die Konvergenz Fortschritt.
function ss = trilinearInverse(p, p1,p2,p3,p4,p5,p6,p7,p8, iter)
%Computes the inverse of the trilinear map from [0,1]^3 to the box defined
% by points p1,...,p8, where the points are ordered consistent with
% p1~(0,0,0), p2~(0,0,1), p3~(0,1,0), p4~(1,0,0), p5~(0,1,1),
% p6~(1,0,1), p7~(1,1,0), p8~(1,1,1)
%Uses Gauss-Newton method. Inputs must be column vectors.
tol = 1e-9;
ss = [0.5; 0.5; 0.5]; %initial guess
for k=1:iter
s = ss(1);
t = ss(2);
w = ss(3);
r = p1*(1-s)*(1-t)*(1-w) + p2*s*(1-t)*(1-w) + ...
p3*(1-s)*t*(1-w) + p4*(1-s)*(1-t)*w + ...
p5*s*t*(1-w) + p6*s*(1-t)*w + ...
p7*(1-s)*t*w + p8*s*t*w - p;
disp(['k= ', num2str(k), ...
', residual norm= ', num2str(norm(r)),...
', [s,t,w]= ',num2str([s,t,w])])
if (norm(r) < tol)
break
end
Js = -p1*(1-t)*(1-w) + p2*(1-t)*(1-w) + ...
-p3*t*(1-w) - p4*(1-t)*w + ...
p5*t*(1-w) + p6*(1-t)*w + ...
-p7*t*w + p8*t*w;
Jt = -p1*(1-s)*(1-w) - p2*s*(1-w) + ...
p3*(1-s)*(1-w) - p4*(1-s)*w + ...
p5*s*(1-w) - p6*s*w + ...
p7*(1-s)*w + p8*s*w;
Jw = -p1*(1-s)*(1-t) - p2*s*(1-t) + ...
-p3*(1-s)*t + p4*(1-s)*(1-t) + ...
-p5*s*t + p6*s*(1-t) + ...
p7*(1-s)*t + p8*s*t;
J = [Js,Jt,Jw];
ss = ss - J\r;
end
end
Beispiel Verwendung:
%test_trilinearInverse.m
h = 0.25;
p1 = [0;0;0] + h*randn(3,1);
p2 = [0;0;1] + h*randn(3,1);
p3 = [0;1;0] + h*randn(3,1);
p4 = [1;0;0] + h*randn(3,1);
p5 = [0;1;1] + h*randn(3,1);
p6 = [1;0;1] + h*randn(3,1);
p7 = [1;1;0] + h*randn(3,1);
p8 = [1;1;1] + h*randn(3,1);
s0 = rand;
t0 = rand;
w0 = rand;
p = p1*(1-s0)*(1-t0)*(1-w0) + p2*s0*(1-t0)*(1-w0) + ...
p3*(1-s0)*t0*(1-w0) + p4*(1-s0)*(1-t0)*w0 + ...
p5*s0*t0*(1-w0) + p6*s0*(1-t0)*w0 + ...
p7*(1-s0)*t0*w0 + p8*s0*t0*w0;
ss = trilinearInverse(p, p1,p2,p3,p4,p5,p6,p7,p8);
disp(['error= ', num2str(norm(ss - [s0;t0;w0]))])
Beispiel-Ausgabe:
test_trilinearInverse
k= 1, residual norm= 0.38102, [s,t,w]= 0.5 0.5 0.5
k= 2, residual norm= 0.025324, [s,t,w]= 0.37896 0.59901 0.17658
k= 3, residual norm= 0.00037108, [s,t,w]= 0.40228 0.62124 0.15398
k= 4, residual norm= 9.1441e-08, [s,t,w]= 0.40218 0.62067 0.15437
k= 5, residual norm= 3.3548e-15, [s,t,w]= 0.40218 0.62067 0.15437
error= 4.8759e-15
Man muss vorsichtig sein, über die Bestellung des Eingabe-Punkte, da die inverse multilineares interpolation ist nur dann gut definiert, wenn die Form hat positive Volumen-und in 3D-es ist viel einfacher Punkte, die die Form drehen inside-out von selbst.
Hier ist meine Implementierung von Naaff-Lösung, als Community Wiki. Nochmals vielen Dank.
Dies ist eine C-Implementierung, soll aber auf c ++ arbeiten. Es enthält eine Fuzzing-Funktion.
#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;
}
Da Sie in 2D, Ihre bilerp Funktion arbeiten ist wirklich 2 Gleichungen, 1 für x und y 1 für. Sie können in der Form neu geschrieben werden:
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
Dabei gilt:
A = p3 - p2 - p1 + p0
B = p2 - p0
C = p1 - p0
D = p0
Schreiben Sie die erste Gleichung t in Bezug auf s, Ersatz in die zweite zu bekommen, und lösen für s.
Das ist meine Implementierung ... Ich denke, es ist schneller als die 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;
*/
}
und Ableitung
// 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
Wenn alles, was Sie haben ein einzelner Wert von p ist, so dass p zwischen der Minimal- und Maximalwerten an den vier Ecken des Platzes ist, dann nein, ist es nicht möglich, in der Regel eine einzige Lösung (s zu finden, t ), so dass der bilineare Interpolant Sie diesen Wert geben.
In der Regel gibt es innerhalb des Quadrats eine unendliche Anzahl von Lösungen (s, t) sein. Sie werden liegen entlang einem gekrümmten (hyperbolischen) Pfad durch das Quadrat.
Wenn Ihre Funktion wird ein Vektor eines Wert, so dass man zwei bekannte Werte an einem unbekannten Punkt auf dem Platz hat? Bei bekannten Werten von zwei Parametern an jeder Ecke des Platzes, dann kann eine Lösung gibt, aber es gibt keine Garantie dafür. Denken Sie daran, dass wir dies als zwei getrennte, unabhängige Probleme anzeigen können. Die Lösung wird auf jedem von ihnen wird entlang einer hyperbolischen Konturlinie durch das Quadrat liegen. Wenn das Paar von Konturen innerhalb des Quadrats überquert, dann existiert eine Lösung. Wenn sie nicht überqueren, dann existiert keine Lösung.
Sie auch fragen, ob eine einfache Formel, um das Problem zu lösen, besteht. Sorry, aber nicht wirklich, dass ich sehe. Wie gesagt, sind die Kurven hyperbolische.
Eine Lösung ist, eine andere Methode der Interpolation zu wechseln. Anstatt also bilinear, brechen Sie den Platz in ein Paar Dreiecke. Innerhalb jedes Dreieck können wir nun wirklich eine lineare Interpolation verwendet werden. So, jetzt können wir das lineare System von 2 Gleichungen in zwei Unbekannten innerhalb jedes Dreiecks lösen. Es kann in jedem Dreieck eine Lösung sein, mit Ausnahme eines seltenen entarteten Fall, in dem die entsprechenden abschnittsweise linearen Konturlinien gerade ist koinzident.
Einige der Antworten Ihre Frage etwas falsch interpretiert. dh. sie gehen davon aus Sie die Wert einer unbekannten Funktion interpoliert, anstatt einer interpolierten Position p (x, y) innerhalb des Quad gegeben sind, dass Sie die (s, t) Koordinaten von finden möchten. Dies ist ein einfacheres Problem und es ist garantiert eine Lösung sein, dass der Schnittpunkt zweier Geraden durch den Quad ist.
Eine der Linien der Segmente P0P1 und P2P3 durchtrennt, wird der andere Schnitt durch P0P2 und p1p3, ähnlich zu der Achse ausgerichteten Fall. Diese Leitungen sind eindeutig durch die Position des p definiert (x, y), und wird natürlich an diesem Punkt schneiden.
Unter Berücksichtigung nur die Zeile, die durch P0P1 und P2P3 schneidet, können wir eine Familie solcher Linien vorstellen, für jeden unterschiedlichen s-Wert wählen wir, die jeweils das Quad mit einer anderen Breite zu schneiden. Wenn wir einen s-Wert zu beheben, können wir die beiden Endpunkte finden von t = 0 und t = 1 zu setzen.
So nimmt zuerst die Zeile hat die Form: y = a0 * x + b 0
Dann wissen wir, zwei Endpunkte dieser Linie, wenn wir einen bestimmten Wert s zu beheben. Sie sind:
(1-s) p0 + (e) p1
(1-s) p2 + (e) p3
Vor dem Hintergrund dieser beiden Endpunkte, können wir die Familie von Linien bestimmen, indem sie in die Gleichung für die Linie, und die Lösung für a0 und b0 als Funktionen von s anschließen. einen s-Wert Einstellung gibt die Formel für eine bestimmte Leitung. Alles, was wir jetzt brauchen, ist die Linie in dieser Familie trifft unseren Punkt p (x, y), um herauszufinden. Einfach die Koordinaten von p Aufstecken (x, y) in unsere Linie Formel können wir dann für den Zielwert von s lösen.
Der entsprechende Ansatz kann getan werden, t als auch zu finden.
Nun, wenn p ein 2D-Punkt ist, ja, man kann es leicht zu bekommen. In diesem Fall ist S die Bruchkomponente der Gesamtbreite des Vierecks an T, T ebenfalls die Bruchkomponente der Gesamthöhe des Vierecks an S ist.
Wenn jedoch p ein Skalar ist, ist es nicht unbedingt möglich, weil die bilineare Interpolation Funktion nicht notwendigerweise monolithisch ist.
