Interpolation bilinéaire inverse?
https://stackoverflow.com/questions/808441
-
03-07-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianQuestion
J'ai quatre 2d points, p0 = (x0, y0), p1 = (x1, y1), etc. qui forment un quadrilatère. Dans mon cas, le quad n’est pas rectangulaire, mais il devrait au moins être convexe.
p2 --- p3
| |
t | p |
| |
p0 --- p1
s
J'utilise une interpolation bilinéaire. S et T sont dans [0..1] et le point interpolé est donné par:
bilerp(s,t) = t*(s*p3+(1-s)*p2) + (1-t)*(s*p1+(1-s)*p0)
Voici le problème .. J'ai un point 2d p que je sais être à l'intérieur du quad. Je veux trouver le s, t qui me donnera ce point lors de l'utilisation de l'interpolation bilinéaire.
Existe-t-il une formule simple pour inverser l'interpolation bilinéaire?
Merci pour les solutions. J'ai posté mon implémentation de la solution de Naaff en tant que wiki.
La solution
Je pense qu'il est plus facile de concevoir votre problème comme un problème d'intersection: quel est l'emplacement du paramètre (s, t) où le point p intersecte la surface bilinéaire bidimensionnelle arbitraire définie par p0, p1, p2 et p3.
L’approche que je vais adopter pour résoudre ce problème est similaire à la suggestion de tspauld.
Commencez par deux équations en termes de x et de 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 )
Résoudre pour 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) )
Nous pouvons maintenant définir ces deux équations égales pour éliminer t. En déplaçant tout vers la gauche et en simplifiant, nous obtenons une équation de la forme:
A*(1-s)^2 + B*2s(1-s) + C*s^2 = 0
Où:
A = (p0-p) X (p0-p2)
B = ( (p0-p) X (p1-p3) + (p1-p) X (p0-p2) ) / 2
C = (p1-p) X (p1-p3)
Notez que j'ai utilisé l'opérateur X pour désigner la croix 2D. produit (par exemple, p0 X p1 = x0 * y1 - y0 * x1). J'ai formulé cette équation en tant que quadratique polynôme de Bernstein , ce qui en fait les choses sont plus élégantes et plus stables numériquement. Les solutions à s sont les racines de cette équation. Nous pouvons trouver les racines en utilisant la formule quadratique pour les polynômes de Bernstein:
s = ( (A-B) +- sqrt(B^2 - A*C) ) / ( A - 2*B + C )
La formule quadratique donne deux réponses en raison du + -. Si vous ne vous intéressez qu'aux solutions où p se trouve dans la surface bilinéaire, vous pouvez ignorer toute réponse où s n'est pas compris entre 0 et 1. Pour trouver t, substituez simplement s dans l'une des deux équations ci-dessus où nous avons résolu en termes de s.
Je devrais signaler un cas spécial important. Si le dénominateur A - 2 * B + C = 0, votre polynôme quadratique est en réalité linéaire. Dans ce cas, vous devez utiliser une équation beaucoup plus simple pour trouver s:
s = A / (A-C)
Cela vous donnera exactement une solution, sauf si AC = 0. Si A = C, vous avez deux cas: A = C = 0 signifie que toutes les valeurs de s contiennent p, sinon aucune valeur pour s ne contient p.
Autres conseils
On pourrait utiliser la la méthode de Newton pour effectuer une itération résolvez le système d'équations non linéaire suivant:
p = p0*(1-s)*(1-t) + p1*s*(1-t) + p2*s*t + p3*(1-s)*t.
Notez qu'il existe deux équations (égalité dans les composantes x et y de l'équation) et deux inconnues (s et t).
Pour appliquer la méthode de Newton, nous avons besoin du reste , qui est:
r = p - (p0*(1-s)*(1-t) + p1*s*(1-t) + p2*s*t + p3*(1-s)*t),
et la matrice jacobienne J, qui est trouvée en différenciant le résidu. Pour notre problème, le jacobien est:
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).
Pour utiliser la méthode de Newton, on commence par une première estimation (s, t), puis on effectue l'itération suivante plusieurs fois:
(s,t) = (s,t) - J^-1 r,
recalculant J et r chaque itération avec les nouvelles valeurs de s et t. À chaque itération, le coût dominant consiste à appliquer l'inverse du jacobien au résiduel (J ^ -1 r) en résolvant un système linéaire 2x2 avec J comme matrice de coefficients et r comme côté droit.
Intuition de la méthode:
Intuitivement, si le quadrilatère était un parallélogramme , il serait beaucoup plus facile de résoudre le problème. . La méthode de Newton consiste à résoudre le problème du quadrilatère avec des approximations successives de parallélogrammes. A chaque itération, nous
-
Utilisez les informations dérivées locales au point (s, t) pour approcher le quadrilatère avec un parallélogramme.
-
Trouvez les valeurs correctes (s, t) sous l'approximation du parallélogramme en résolvant un système linéaire.
-
Aller à ce nouveau point et répéter.
Avantages de la méthode:
Comme prévu pour les méthodes de type Newton, la convergence est extrêmement rapide. Au fur et à mesure que les itérations avancent, non seulement la méthode se rapproche-t-elle du point réel, mais les approximations du parallélogramme local deviennent également plus précises, de sorte que le taux de convergence lui-même augmente (dans le jargon des méthodes itératives, nous disons que la méthode de Newton est de manière quadratique convergente ). En pratique, cela signifie que le nombre de chiffres corrects double environ à chaque itération.
Voici un tableau représentatif des itérations par rapport aux erreurs d'un essai aléatoire que j'ai effectué (voir le code ci-dessous):
Iteration Error
1 0.0610
2 9.8914e-04
3 2.6872e-07
4 1.9810e-14
5 5.5511e-17 (machine epsilon)
Après deux itérations, l'erreur est suffisamment petite pour être efficacement imperceptible et suffisante pour la plupart des tâches pratiques. Après cinq itérations, le résultat est précis dans les limites de précision de la machine .
Si vous corrigez le nombre d'itérations (par exemple, 3 itérations pour la plupart des applications pratiques et 8 itérations si vous avez besoin d'une très grande précision), alors l'algorithme a une logique très simple et directe, avec une structure qui se prête bien. au calcul haute performance. Il n'est pas nécessaire de vérifier toutes sortes de cas particuliers * et d'utiliser une logique différente en fonction des résultats. C'est juste une boucle for contenant quelques formules simples. Ci-dessous, je souligne plusieurs avantages de cette approche par rapport aux approches traditionnelles basées sur des formules trouvées dans d’autres réponses ici et sur Internet:
-
Facile à coder . Il suffit de créer une boucle for et de saisir quelques formules.
-
Aucune condition ni ramification (si / alors), ce qui permet généralement de beaucoup mieux efficacité du traitement en pipeline .
-
Pas de racines carrées, et seulement 1 division requise par itération (si bien écrit). Toutes les autres opérations sont de simples additions, soustractions et multiplications. Les racines carrées et les divisions sont généralement plusieurs fois plus lentes que les additions / multiplications / multiplications, et peuvent foirer cache efficacité sur certaines architectures (plus particulièrement sur certains systèmes embarqués). En effet, si vous regardez sous le capot, comment racines carrées et divisions sont en fait calculées par les langages de programmation modernes, elles utilisent toutes deux des variantes de la méthode de Newton, tantôt dans le matériel, tantôt dans le logiciel, architecture.
-
Peut être facilement vectorisé pour travailler à la fois sur des tableaux comportant un grand nombre de quadrilatères. Voir mon code vectorisé ci-dessous pour un exemple de procédure.
-
S'étend aux dimensions arbitraires . L'algorithme s'étend de manière simple à l'interpolation multilinéaire inverse dans un nombre quelconque de dimensions (2d, 3d, 4d, ...). J'ai inclus une version 3D ci-dessous, et on peut imaginer écrire une version récursive simple ou utiliser des bibliothèques de différenciation automatique pour passer à n dimensions. La méthode de Newton présente généralement des vitesses de convergence indépendantes de la dimension. En principe, elle doit donc être extensible à quelques milliers de dimensions (!) Sur le matériel actuel (après quoi la construction et la résolution de la matrice n-par-n J seront probablement les limites facteur).
Bien sûr, la plupart de ces facteurs dépendent également du matériel, du compilateur et de nombreux autres facteurs, de sorte que votre kilométrage peut varier.
Code:
Quoi qu'il en soit, voici mon code Matlab: (je libère tout ici du domaine public)
Version 2D de 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
Exemple d'utilisation:
% 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)])
Exemple de sortie:
>> test_bilinearInverse
Error after 5 iterations: 1.5701e-16
Version 2D vectorisée rapide:
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
Pour accélérer ce code utilise implicitement la formule suivante pour l'inverse d'une matrice 2x2:
[a,b;c,d]^-1 = (1/(ad-bc))[d, -b; -c, a]
Exemple d'utilisation:
% 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)])
Exemple de sortie:
>> 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
Version 3D:
Inclut du code pour afficher la progression de la convergence.
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
Exemple d'utilisation:
%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]))])
Exemple de sortie:
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
Il faut faire attention à l'ordre des points en entrée, car l'interpolation multilinéaire inverse n'est bien définie que si la forme a un volume positif, et en 3D, il est beaucoup plus facile de choisir des points qui font tourner la forme à l'envers de lui-même.
Voici mon implémentation de la solution de Naaff, en tant que wiki de communauté. Merci encore.
Ceci est une implémentation en C, mais devrait fonctionner en c ++. Il inclut une fonction de 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;
}
Puisque vous travaillez en 2D, votre fonction bilerp est en réalité 2 équations, 1 pour x et 1 pour y. Ils peuvent être réécrits sous la forme:
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
Où:
A = p3 - p2 - p1 + p0
B = p2 - p0
C = p1 - p0
D = p0
Réécrivez la première équation pour obtenir t en termes de s , remplacez-le dans le second et résolvez s .
c'est ma mise en oeuvre ... Je suppose que c'est plus rapide que de 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;
*/
}
et dérivation
// 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
Si tout ce que vous avez est une valeur unique de p, telle que p se situe entre les valeurs minimale et maximale aux quatre coins du carré, alors non, il est en général impossible de trouver une solution UNIQUE (s, t ) de telle sorte que l’interpolant bilinéaire vous donne cette valeur.
En général, il y aura un nombre infini de solutions (s, t) dans le carré. Ils se situeront le long d'un chemin incurvé (hyperbolique) à travers le carré.
Si votre fonction est une valeur vectorielle, vous avez donc deux valeurs connues en un point inconnu du carré? Étant donné les valeurs connues de deux paramètres à chaque coin du carré, une solution PEUT alors exister, mais rien ne le garantit. Rappelez-vous que nous pouvons voir cela comme deux problèmes distincts et indépendants. La solution à chacun d'eux sera située le long d'une ligne de contour hyperbolique à travers le carré. Si la paire de contours se croise à l'intérieur du carré, une solution existe. S'ils ne se croisent pas, alors aucune solution n'existe.
Vous demandez également s'il existe une formule simple pour résoudre le problème. Désolé, mais pas vraiment que je vois. Comme je l'ai dit, les courbes sont hyperboliques.
Une solution consiste à changer de méthode d’interpolation. Donc, au lieu de bilinéaire, divisez le carré en une paire de triangles. Dans chaque triangle, nous pouvons maintenant utiliser une interpolation vraiment linéaire. Nous pouvons maintenant résoudre le système linéaire de 2 équations à 2 inconnues dans chaque triangle. Il peut y avoir une solution dans chaque triangle, sauf dans un cas rare et dégénéré dans lequel les courbes de niveau linéaires par morceaux correspondantes sont coïncidentes.
Certaines réponses ont légèrement mal interprété votre question. c'est à dire. ils supposent que vous recevez la valeur d'une fonction interpolée inconnue, plutôt qu'une position interpolée p (x, y) à l'intérieur du quad dont vous voulez trouver les coordonnées (s, t). C’est un problème plus simple et il est garanti qu’une solution sera l’intersection de deux lignes droites à travers le quad.
L’une des lignes coupe les segments p0p1 et p2p3, l’autre les lignes p0p2 et p1p3, comme dans le cas des alignements d’axes. Ces lignes sont uniquement définies par la position de p (x, y) et se recouperont évidemment à cet endroit.
En considérant uniquement la ligne qui coupe p0p1 et p2p3, nous pouvons imaginer une famille de telles lignes, pour chaque valeur s choisie, coupant le quad à une largeur différente. Si nous fixons une valeur s, nous pouvons trouver les deux extrémités en définissant t = 0 et t = 1.
Supposons d’abord que la ligne a la forme suivante: y = a0 * x + b0
Ensuite, nous connaissons deux extrémités de cette ligne, si nous corrigeons une valeur s donnée. Ils sont:
(1-s) p0 + (s) p1
(1-s) p2 + (s) p3
Étant donné ces deux extrémités, nous pouvons déterminer la famille de lignes en les connectant à l’équation de la ligne et en résolvant pour a0 et b0 en tant que fonctions de s . Définir une valeur s donne la formule pour une ligne spécifique. Tout ce dont nous avons besoin maintenant, c'est de déterminer quelle ligne de cette famille touche notre point p (x, y). En branchant simplement les coordonnées de p (x, y) dans notre formule de ligne, nous pouvons alors résoudre la valeur cible de s.
L’approche correspondante peut également être utilisée pour trouver t.
Eh bien, si p est un point 2D, oui, vous pouvez l’obtenir facilement. Dans ce cas, S est la composante fractionnelle de la largeur totale du quadrilatère en T, T est également la composante fractionnelle de la hauteur totale du quadrilatère en S.
Si, toutefois, p est un scalaire, ce n'est pas forcément possible, car la fonction d'interpolation bilinéaire n'est pas nécessairement monolithique.
