Interpolazione bilineare inversa?
https://stackoverflow.com/questions/808441
-
03-07-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianDomanda
Ho quattro punti 2d, p0 = (x0, y0), p1 = (x1, y1), ecc. che formano un quadrilatero. Nel mio caso, il quad non è rettangolare, ma dovrebbe almeno essere convesso.
p2 --- p3
| |
t | p |
| |
p0 --- p1
s
Sto usando l'interpolazione bilineare. S e T sono entro [0..1] e il punto interpolato è dato da:
bilerp(s,t) = t*(s*p3+(1-s)*p2) + (1-t)*(s*p1+(1-s)*p0)
Ecco il problema .. Ho un punto 2d p che so sia all'interno del quad. Voglio trovare la s, t che mi darà quel punto quando uso l'interpolazione bilineare.
Esiste una formula semplice per invertire l'interpolazione bilineare?
Grazie per le soluzioni. Ho pubblicato la mia implementazione della soluzione di Naaff come wiki.
Soluzione
Penso che sia più facile pensare al tuo problema come a un problema di intersezione: qual è la posizione del parametro (s, t) in cui il punto p interseca la superficie bilineare 2D arbitraria definita da p0, p1, p2 e p3.
L'approccio che seguirò per risolvere questo problema è simile al suggerimento di tspauld.
Inizia con due equazioni in termini di xey:
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 )
Risoluzione per 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) )
Ora possiamo impostare queste due equazioni uguali tra loro per eliminare t. Spostando tutto sul lato sinistro e semplificando otteniamo un'equazione del modulo:
A*(1-s)^2 + B*2s(1-s) + C*s^2 = 0
Dove:
A = (p0-p) X (p0-p2)
B = ( (p0-p) X (p1-p3) + (p1-p) X (p0-p2) ) / 2
C = (p1-p) X (p1-p3)
Nota che ho usato l'operatore X per indicare la cross 2D prodotto (ad es. p0 X p1 = x0 * y1 - y0 * x1). Ho formattato questa equazione come Bernstein polynomial quadratica cose più eleganti ed è più numericamente stabile. Le soluzioni a s sono le radici di questa equazione. Possiamo trovare le radici usando la formula quadratica per i polinomi di Bernstein:
s = ( (A-B) +- sqrt(B^2 - A*C) ) / ( A - 2*B + C )
La formula quadratica dà due risposte a causa del + -. Se sei interessato solo a soluzioni in cui p si trova all'interno della superficie bilineare, puoi scartare qualsiasi risposta in cui s non è compreso tra 0 e 1. Per trovare t, sostituisci semplicemente s in una delle due equazioni sopra dove abbiamo risolto t in termini di s.
Dovrei sottolineare un caso speciale importante. Se il denominatore A - 2 * B + C = 0, il polinomio quadratico è effettivamente lineare. In questo caso, devi usare un'equazione molto più semplice per trovare s:
s = A / (A-C)
Questo ti darà esattamente una soluzione, a meno che AC = 0. Se A = C hai due casi: A = C = 0 significa tutti i valori per s contengono p, altrimenti nessun valori per s contiene p.
Altri suggerimenti
Si potrebbe usare Metodo di Newton per iterativamente risolvere il seguente sistema di equazioni non lineari:
p = p0*(1-s)*(1-t) + p1*s*(1-t) + p2*s*t + p3*(1-s)*t.
Nota che ci sono due equazioni (uguaglianza in entrambi i componenti xey dell'equazione) e due incognite (s e t).
Per applicare il metodo di Newton, abbiamo bisogno della residual r, che è:
r = p - (p0*(1-s)*(1-t) + p1*s*(1-t) + p2*s*t + p3*(1-s)*t),
e Matrice giacobina J, che si trova differenziando il residuo. Per il nostro problema il giacobino è:
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).
Per usare il metodo di Newton, si inizia con un'ipotesi iniziale (s, t), quindi si esegue la seguente iterazione una manciata di volte:
(s,t) = (s,t) - J^-1 r,
ricalcola J e r ogni iterazione con i nuovi valori di s e t. Ad ogni iterazione il costo dominante sta nell'applicare l'inverso del giacobino al residuo (J ^ -1 r), risolvendo un sistema lineare 2x2 con J come matrice di coefficiente e r come lato destro.
Intuizione per il metodo:
Intuitivamente, se il quadrilatero fosse un parallelogramma , sarebbe molto più semplice risolvere il problema . Il metodo di Newton sta risolvendo il problema quadrilatero con successive approssimazioni in parallelogramma. Ad ogni iterazione
-
Usa le informazioni derivate locali nel punto (s, t) per approssimare il quadrilatero con un parallelogramma.
-
Trova i valori corretti (s, t) sotto l'approssimazione del parallelogramma, risolvendo un sistema lineare.
-
Vai a questo nuovo punto e ripeti.
Vantaggi del metodo:
Come previsto per i metodi di tipo Newton, la convergenza è estremamente rapida. Man mano che le iterazioni procedono, non solo il metodo si avvicina sempre di più al punto vero, ma anche le approssimazioni del parallelogramma locale diventano più accurate, quindi il tasso di convergenza stesso aumenta (nel gergo dei metodi iterativi, diciamo che il metodo di Newton è convergentemente quadratico ). In pratica ciò significa che il numero di cifre corrette raddoppia approssimativamente ogni iterazione.
Ecco una tabella rappresentativa delle iterazioni vs. errore su una prova casuale che ho fatto (vedi sotto per il codice):
Iteration Error
1 0.0610
2 9.8914e-04
3 2.6872e-07
4 1.9810e-14
5 5.5511e-17 (machine epsilon)
Dopo due iterazioni l'errore è abbastanza piccolo da essere effettivamente inosservabile e abbastanza buono per la maggior parte degli scopi pratici, e dopo 5 iterazioni il risultato è accurato ai limiti di precisione della macchina .
Se correggi il numero di iterazioni (diciamo, 3 iterazioni per la maggior parte delle applicazioni pratiche e 8 iterazioni se hai bisogno di una precisione molto elevata), l'algoritmo ha una logica molto semplice e diretta, con una struttura che si presta bene al calcolo ad alte prestazioni. Non è necessario controllare tutti i tipi di casi limite speciali * e utilizzare una logica diversa a seconda dei risultati. È solo un ciclo for contenente alcune semplici formule. Di seguito metto in evidenza diversi vantaggi di questo approccio rispetto agli approcci basati sulla formula tradizionale che si trovano in altre risposte qui e su Internet:
Facile da codificare . Basta creare un ciclo for e digitare alcune formule.
Nessun condizionale o ramificazione (if / then), che in genere consente di migliorare molto efficienza di pipeline .
Nessuna radice quadrata e solo 1 divisione richiesta per iterazione (se scritta bene). Tutte le altre operazioni sono semplici aggiunte, sottrazioni e moltiplicazioni. Le radici e le divisioni quadrate sono in genere molte volte più lente delle aggiunte / moltiplicazioni / moltiplicazioni e possono rovinare cache efficienza su alcune architetture (in particolare su determinati sistemi integrati). In effetti, se guardi sotto il cofano come radici quadrate e divisioni sono in realtà calcolati dai moderni linguaggi di programmazione, entrambi utilizzano varianti del metodo di Newton, a volte in hardware e talvolta in software a seconda del l'architettura.
Può essere facilmente vettorializzato per lavorare su array con un numero enorme di quadrilateri contemporaneamente. Vedi il mio codice vettoriale qui sotto per un esempio di come farlo.
Si estende a dimensioni arbitrarie . L'algoritmo si estende in modo semplice all'inverso dell'interpolazione multilineare in qualsiasi numero di dimensioni (2d, 3d, 4d, ...). Ho incluso una versione 3D di seguito, e si può immaginare di scrivere una semplice versione ricorsiva o di utilizzare librerie di differenziazione automatiche per passare alle dimensioni n. Il metodo di Newton mostra tipicamente tassi di convergenza indipendenti dalla dimensione, quindi in linea di principio il metodo dovrebbe essere scalabile a qualche migliaio di dimensioni (!) Sull'hardware corrente (dopodiché la costruzione e la risoluzione della matrice n-by-n J sarà probabilmente il limite fattore).
Naturalmente, la maggior parte di queste cose dipende anche dall'hardware, dal compilatore e da una serie di altri fattori, quindi il tuo chilometraggio può variare.
Codice:
Comunque, ecco il mio codice Matlab: (rilascio tutto qui di dominio pubblico)
Versione 2D di base:
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
Esempio di utilizzo:
% 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)])
Esempio di output:
>> test_bilinearInverse
Error after 5 iterations: 1.5701e-16
Versione 2D vettoriale veloce:
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
Per velocità questo codice usa implicitamente la seguente formula per l'inverso di una matrice 2x2:
[a,b;c,d]^-1 = (1/(ad-bc))[d, -b; -c, a]
Esempio di utilizzo:
% 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)])
Esempio di output:
>> 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
Versione 3D:
Include del codice per visualizzare l'avanzamento della convergenza.
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
Esempio di utilizzo:
%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]))])
Esempio di output:
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
Bisogna stare attenti all'ordinamento dei punti di input, poiché l'interpolazione multilineare inversa è ben definita solo se la forma ha volume positivo, e in 3D è molto più facile scegliere i punti che fanno girare la forma al rovescio di se stesso.
Ecco la mia implementazione della soluzione di Naaff, come wiki della community. Grazie ancora.
Questa è un'implementazione in C, ma dovrebbe funzionare su c ++. Include una funzione di test fuzz.
#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;
}
Dato che lavori in 2D, la tua funzione bilerp è in realtà 2 equazioni, 1 per xe 1 per y. Possono essere riscritti nel modulo:
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
Dove:
A = p3 - p2 - p1 + p0
B = p2 - p0
C = p1 - p0
D = p0
Riscrivi la prima equazione per ottenere t in termini di s , sostituisci con la seconda e risolvi s .
questa è la mia implementazione ... immagino sia più veloce di 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;
*/
}
e derivazione
// 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
Se tutto ciò che hai è un singolo valore di p, tale che p sia compreso tra i valori minimo e massimo ai quattro angoli del quadrato, allora no, in generale non è possibile trovare una SOLA soluzione (s, t ) tale che l'interpolante bilineare ti darà quel valore.
In generale, ci sarà un numero infinito di soluzioni (s, t) all'interno del quadrato. Stenderanno lungo un percorso curvo (iperbolico) attraverso il quadrato.
Se la tua funzione ha un valore vettoriale, allora hai due valori noti in un punto sconosciuto nel quadrato? Dati i valori noti di due parametri in ciascun angolo del quadrato, allora può esistere una soluzione, ma non c'è alcuna garanzia. Ricorda che possiamo vederlo come due problemi separati e indipendenti. La soluzione a ciascuno di essi si troverà lungo una linea di contorno iperbolica attraverso il quadrato. Se la coppia di contorni si interseca all'interno del quadrato, allora esiste una soluzione. Se non attraversano, allora non esiste alcuna soluzione.
Chiedete anche se esiste una formula semplice per risolvere il problema. Scusa, ma non proprio quello che vedo. Come ho detto, le curve sono iperboliche.
Una soluzione è passare a un diverso metodo di interpolazione. Quindi, anziché bilineare, spezza il quadrato in una coppia di triangoli. All'interno di ogni triangolo, ora possiamo usare un'interpolazione veramente lineare. Quindi ora possiamo risolvere il sistema lineare di 2 equazioni in 2 incognite all'interno di ciascun triangolo. Può esserci una soluzione in ciascun triangolo, ad eccezione di un raro caso degenerato in cui le corrispondenti linee di contorno lineari a tratti risultano coincidenti.
Alcune delle risposte hanno leggermente frainteso la tua domanda. vale a dire. stanno assumendo che ti venga dato il valore di una funzione interpolata sconosciuta, piuttosto che una posizione interpolata p (x, y) all'interno del quad di cui vuoi trovare le coordinate (s, t). Questo è un problema più semplice e c'è sicuramente una soluzione che è l'intersezione di due rette attraverso il quad.
Una delle linee taglierà attraverso i segmenti p0p1 e p2p3, l'altra taglierà attraverso p0p2 e p1p3, in modo simile al caso allineato agli assi. Queste linee sono definite in modo univoco dalla posizione di p (x, y) e si intersecheranno ovviamente a questo punto.
Considerando solo la linea che taglia p0p1 e p2p3, possiamo immaginare una famiglia di tali linee, per ogni diverso valore s che scegliamo, tagliando il quadrante con una larghezza diversa. Se fissiamo un valore s, possiamo trovare i due endpoint impostando t = 0 et t = 1.
Quindi prima supponiamo che la linea abbia la forma: y = a0 * x + b0
Quindi conosciamo due endpoint di questa linea, se fissiamo un dato valore s. Sono:
(1-s) p0 + (s) p1
(1-s) p2 + (s) p3
Dati questi due punti finali, possiamo determinare la famiglia di linee collegandole all'equazione per la linea e risolvendo per a0 e b0 come funzioni di s . L'impostazione di un valore s fornisce la formula per una riga specifica. Tutto ciò di cui abbiamo bisogno ora è capire quale linea di questa famiglia colpisce il nostro punto p (x, y). Inserendo semplicemente le coordinate di p (x, y) nella nostra formula di linea, possiamo quindi risolvere il valore target di s.
L'approccio corrispondente può essere fatto anche per trovare t.
Bene, se p è un punto 2D, sì, puoi ottenerlo facilmente. In tal caso, S è la componente frazionaria della larghezza totale del quadrilatero a T, T è anche la componente frazionaria dell'altezza totale del quadrilatero a S.
Se, tuttavia, p è uno scalare, non è necessariamente possibile, poiché la funzione di interpolazione bilineare non è necessariamente monolitica.
