math - interpolacion - ¿Interpolación bilineal inversa?
interpolacion lineal 3 variables (8)
Tengo cuatro puntos 2d, p0 = (x0, y0), p1 = (x1, y1), etc. que forman un cuadrilátero. En mi caso, el quad no es rectangular, pero al menos debería ser convexo.
p2 --- p3
| |
t | p |
| |
p0 --- p1
s
Estoy usando interpolación bilineal. S y T están dentro de [0..1] y el punto interpolado viene dado por:
bilerp(s,t) = t*(s*p3+(1-s)*p2) + (1-t)*(s*p1+(1-s)*p0)
Aquí está el problema ... Tengo un punto 2d que sé que está dentro del patio. Quiero encontrar la s, t que me dará ese punto al usar la interpolación bilineal.
¿Existe una fórmula simple para revertir la interpolación bilineal?
Gracias por las soluciones. Publiqué mi implementación de la solución de Naaff como wiki.
Algunas de las respuestas han malinterpretado ligeramente tu pregunta. es decir. asumen que se le asigna el valor de una función interpolada desconocida, en lugar de una posición interpolada p (x, y) dentro del cuadrante en el que desea encontrar las coordenadas (s, t). Este es un problema más simple y se garantiza que existe una solución que es la intersección de dos líneas rectas a través del cuadrante.
Una de las líneas cortará a través de los segmentos p0p1 y p2p3, la otra cortará a través de p0p2 y p1p3, similar a la caja alineada con el eje. Estas líneas están definidas de manera única por la posición de p (x, y), y obviamente se intersecarán en este punto.
Si consideramos solo la línea que corta p0p1 y p2p3, podemos imaginar una familia de esas líneas, por cada valor de s diferente que elijamos, cada uno cortando el cuadrángulo en un ancho diferente. Si fijamos un valor de s, podemos encontrar los dos puntos finales estableciendo t = 0 y t = 1.
Así que primero asuma que la línea tiene la forma: y = a0 * x + b0
Entonces conocemos dos puntos finales de esta línea, si fijamos un valor de s dado. Son:
(1-s) p0 + (s) p1
(1-s) p2 + (s) p3
Dados estos dos puntos finales, podemos determinar la familia de líneas insertándolas en la ecuación de la línea y resolviendo para a0 y b0 como funciones de s . Establecer un valor s da la fórmula para una línea específica. Todo lo que necesitamos ahora es averiguar qué línea de esta familia llega a nuestro punto p (x, y). Simplemente conectando las coordenadas de p (x, y) en nuestra fórmula de línea, podemos resolver el valor objetivo de s.
El enfoque correspondiente se puede hacer para encontrar t también.
Aquí está mi implementación de la solución de Naaff, como un wiki de la comunidad. Gracias de nuevo.
Esta es una implementación de C, pero debería funcionar en c ++. Incluye una función de prueba 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;
}
Bueno, si p es un punto 2D, sí, puedes obtenerlo fácilmente. En ese caso, S es el componente fraccionario del ancho total del cuadrilátero en T, T es también el componente fraccionario de la altura total del cuadrilátero en S
Sin embargo, si p es un escalar, no es necesariamente posible, porque la función de interpolación bilineal no es necesariamente monolítica.
Creo que es más fácil pensar en su problema como un problema de intersección: ¿cuál es la ubicación de los parámetros (s, t) donde el punto p intersecta la superficie bilineal 2D arbitraria definida por p0, p1, p2 y p3?
El enfoque que tomaré para resolver este problema es similar a la sugerencia de tspauld.
Comience con dos ecuaciones en términos de x y 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 )
Resolviendo para 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) )
Ahora podemos establecer estas dos ecuaciones iguales para eliminar t. Moviendo todo hacia el lado izquierdo y simplificando obtenemos una ecuación de la forma:
A*(1-s)^2 + B*2s(1-s) + C*s^2 = 0
Dónde:
A = (p0-p) X (p0-p2)
B = ( (p0-p) X (p1-p3) + (p1-p) X (p0-p2) ) / 2
C = (p1-p) X (p1-p3)
Tenga en cuenta que he usado el operador X para indicar el producto cruzado 2D (por ejemplo, p0 X p1 = x0 * y1 - y0 * x1). He formateado esta ecuación como un polinomio de Bernstein cuadrático ya que hace las cosas más elegantes y numéricamente más estables. Las soluciones a s son las raíces de esta ecuación. Podemos encontrar las raíces usando la fórmula cuadrática para los polinomios de Bernstein:
s = ( (A-B) +- sqrt(B^2 - A*C) ) / ( A - 2*B + C )
La fórmula cuadrática da dos respuestas debido a la + -. Si solo está interesado en soluciones donde p se encuentra dentro de la superficie bilineal, puede descartar cualquier respuesta donde s no esté entre 0 y 1. Para encontrar t, simplemente sustituya s en una de las dos ecuaciones anteriores, donde resolvimos t en términos de s.
Debería señalar un caso especial importante. Si el denominador A - 2 * B + C = 0, entonces su polinomio cuadrático es en realidad lineal. En este caso, debes usar una ecuación mucho más simple para encontrar s:
s = A / (A-C)
Esto le dará exactamente una solución, a menos que AC = 0. Si A = C entonces tiene dos casos: A = C = 0 significa que todos los valores para s contienen p, de lo contrario, ningún valor para s contiene p.
Se podría usar el método de Newton para resolver iterativamente el siguiente sistema de ecuaciones no lineales:
p = p0*(1-s)*(1-t) + p1*s*(1-t) + p2*s*t + p3*(1-s)*t.
Tenga en cuenta que hay dos ecuaciones (igualdad en las componentes x e y de la ecuación) y dos incógnitas (s y t).
Para aplicar el método de Newton, necesitamos la r residual , que es:
r = p - (p0*(1-s)*(1-t) + p1*s*(1-t) + p2*s*t + p3*(1-s)*t),
y la matriz jacobiana J, que se encuentra al diferenciar lo residual. Para nuestro problema el jacobiano es:
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).
Para usar el método de Newton, uno comienza con una estimación inicial (s, t), y luego realiza la siguiente iteración un puñado de veces:
(s,t) = (s,t) - J^-1 r,
recalcular J y r en cada iteración con los nuevos valores de s y t. En cada iteración, el costo dominante es aplicar el inverso del jacobiano al residual (J ^ -1 r), resolviendo un sistema lineal 2x2 con J como la matriz de coeficientes y r como el lado derecho.
Intuición para el método:
Intuitivamente, si el cuadrilátero fuera un parallelogram , sería mucho más fácil resolver el problema. El método de Newton es resolver el problema cuadrilátero con aproximaciones de paralelogramos sucesivos. En cada iteración nosotros
Use información derivada local en el punto (s, t) para aproximar el cuadrilátero con un paralelogramo.
Encuentre los valores correctos (s, t) en la aproximación del paralelogramo, resolviendo un sistema lineal.
Salta a este nuevo punto y repite.
Ventajas del método:
Como se espera para los métodos de tipo Newton, la convergencia es extremadamente rápida. A medida que avanzan las iteraciones, no solo el método se acerca más y más al punto verdadero, sino que las aproximaciones del paralelogramo local también se vuelven más precisas, por lo que la tasa de convergencia aumenta (en la jerga de métodos iterativos, decimos que el método de Newton es cuadráticamente convergente ). En la práctica, esto significa que el número de dígitos correctos duplica aproximadamente cada iteración.
Aquí hay una tabla representativa de iteraciones contra error en un ensayo aleatorio que hice (ver más abajo para el código):
Iteration Error
1 0.0610
2 9.8914e-04
3 2.6872e-07
4 1.9810e-14
5 5.5511e-17 (machine epsilon)
Después de dos iteraciones, el error es lo suficientemente pequeño como para que sea efectivamente imperceptible y lo suficientemente bueno para los propósitos más prácticos, y después de 5 iteraciones, el resultado es preciso hasta los límites de la precisión de la máquina .
Si fija el número de iteraciones (por ejemplo, 3 iteraciones para la mayoría de las aplicaciones prácticas y 8 iteraciones si necesita una precisión muy alta), entonces el algoritmo tiene una lógica muy simple y directa, con una estructura que se presta bien a alto Cálculo del rendimiento. No hay necesidad de verificar todo tipo de casos de bordes especiales *, y usar una lógica diferente según los resultados. Es solo un bucle for que contiene algunas fórmulas simples. A continuación, resalto varias ventajas de este enfoque sobre los enfoques basados en fórmulas tradicionales que se encuentran en otras respuestas aquí y en Internet:
Fácil de codificar . Solo haz un bucle for y escribe algunas fórmulas.
Sin condicionales o ramificaciones (if / then), lo que generalmente permite una eficiencia de canalización mucho mejor.
No hay raíces cuadradas, y solo se requiere 1 división por iteración (si está bien escrito). Todas las demás operaciones son simples sumas, restas y multiplicaciones. Las raíces cuadradas y las divisiones suelen ser varias veces más lentas que las adiciones / multiplicaciones / multiplicaciones, y pueden arruinar la eficiencia de la cache en ciertas arquitecturas (sobre todo en ciertos sistemas integrados). De hecho, si observa cómo las raíces cuadradas y las divisions son calculadas por los lenguajes de programación modernos, ambos usan variantes del método de Newton, algunas veces en hardware y otras en software, dependiendo de la arquitectura.
Se puede vectorizar fácilmente para trabajar en arreglos con un gran número de cuadriláteros a la vez. Vea mi código vectorizado a continuación para ver un ejemplo de cómo hacer esto.
Se extiende a dimensiones arbitrarias . El algoritmo se extiende de forma directa a la interpolación multilineal inversa en cualquier número de dimensiones (2d, 3d, 4d, ...). He incluido una versión en 3D a continuación, y uno puede imaginar escribir una versión recursiva simple, o usar bibliotecas de diferenciación automática para ir a n-dimensiones. El método de Newton normalmente muestra tasas de convergencia independientes de la dimensión, por lo que en principio el método debería ser escalable a unos pocos miles de dimensiones (!) En el hardware actual (después de lo cual la construcción y resolución de la matriz J de n por n probablemente será la limitación factor).
Por supuesto, la mayoría de estas cosas también dependen del hardware, el compilador y muchos otros factores, por lo que su millaje puede variar.
Código:
De todos modos, aquí está mi código de Matlab: (Libero todo aquí al dominio público)
Versión 2D básica:
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
Ejemplo de uso:
% 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)])
Ejemplo de salida:
>> test_bilinearInverse
Error after 5 iterations: 1.5701e-16
Versión 2D vectorizada rápida:
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
Para la velocidad, este código utiliza implícitamente la siguiente fórmula para la inversa de una matriz 2x2:
[a,b;c,d]^-1 = (1/(ad-bc))[d, -b; -c, a]
Ejemplo de uso:
% 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)])
Ejemplo de salida:
>> 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
Versión 3D:
Incluye algún código para mostrar el progreso de convergencia.
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
Ejemplo de uso:
%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]))])
Ejemplo de salida:
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
Hay que tener cuidado con el orden de los puntos de entrada, ya que la interpolación multilineal inversa solo está bien definida si la forma tiene un volumen positivo, y en 3D es mucho más fácil elegir los puntos que hacen que la forma se convierta de adentro hacia afuera.
Si todo lo que tiene es un solo valor de p, tal que p está entre los valores mínimo y máximo en las cuatro esquinas del cuadrado, entonces no, no es posible en general encontrar una SOLA solución (s, t) tal que El interpolante bilineal te dará ese valor.
En general, habrá un número infinito de soluciones (s, t) dentro del cuadrado. Se ubicarán a lo largo de un camino curvo (hiperbólico) a través del cuadrado.
Si su función es un vector de valores, ¿tiene dos valores conocidos en algún punto desconocido del cuadrado? Dados los valores conocidos de dos parámetros en cada esquina del cuadrado, entonces puede existir una solución, pero no hay seguridad de eso. Recuerde que podemos ver esto como dos problemas separados e independientes. La solución para cada uno de ellos será una línea de contorno hiperbólica a través del cuadrado. Si el par de contornos se cruza dentro del cuadrado, entonces existe una solución. Si no se cruzan, entonces no existe solución.
También preguntas si existe una fórmula simple para resolver el problema. Lo siento, pero en realidad no es lo que veo. Como dije, las curvas son hiperbólicas.
Una solución es cambiar a un método diferente de interpolación. Así que en lugar de bilineal, divide el cuadrado en un par de triángulos. Dentro de cada triángulo, ahora podemos usar una interpolación verdaderamente lineal. Así que ahora podemos resolver el sistema lineal de 2 ecuaciones en 2 incógnitas dentro de cada triángulo. Puede haber una solución en cada triángulo, a excepción de un caso degenerado raro en el que las líneas de contorno lineales por partes correspondientes resulten ser coincidentes.
Ya que estás trabajando en 2D, tu función bilerp es realmente 2 ecuaciones, 1 para x y 1 para y. Se pueden reescribir en la forma:
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
Dónde:
A = p3 - p2 - p1 + p0
B = p2 - p0
C = p1 - p0
D = p0
Reescriba la primera ecuación para obtener t en términos de s , sustituya en la segunda y resuelva para s .
esta es mi implementación ... supongo que es más rápida que la 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;
*/
}
y derivación
// 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