Inverse Bilinear Interpolation?
https://stackoverflow.com/questions/808441
-
03-07-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianPergunta
I têm quatro pontos 2d, p0 = (x0, y0), P1 = (x1, y1), etc, que formam um quadrilátero. No meu caso, o quad não é retangular, mas deve pelo menos ser convexo.
p2 --- p3
| |
t | p |
| |
p0 --- p1
s
Eu estou usando interpolação bilinear. S e T são, dentro [0..1] e o ponto interpolado é dada por:
bilerp(s,t) = t*(s*p3+(1-s)*p2) + (1-t)*(s*p1+(1-s)*p0)
Aqui está o problema .. Eu tenho um ponto p 2d que eu sei que está dentro do quad. Eu quero encontrar os s, t que vai me dar esse ponto quando usando interpolação bilinear.
Existe uma fórmula simples para reverter a interpolação bilinear?
Obrigado pelas soluções. Eu postei minha implementação da solução da Naaff como um wiki.
Solução
Eu acho que é mais fácil pensar em seu problema como um problema de interseção: o que é o local parâmetro (s, t), onde o ponto P intercepta a 2D arbitrária bilinear superfície definida por p0, p1, p2 e p3
A abordagem que eu vou tomar para resolver este problema é semelhante à sugestão de tspauld.
Comece com duas equações em termos de X e 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 )
Resolvendo 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) )
Agora podemos definir essas duas equações iguais uns aos outros para eliminar t. Mover tudo para o lado esquerdo e simplificando obtemos uma equação da forma:
A*(1-s)^2 + B*2s(1-s) + C*s^2 = 0
Onde:
A = (p0-p) X (p0-p2)
B = ( (p0-p) X (p1-p3) + (p1-p) X (p0-p2) ) / 2
C = (p1-p) X (p1-p3)
Note que eu usei o operador X para designar o produto cruz cruz 2D produto (por exemplo, X p0 p1 = x0 * y1 - y0 * x1). Eu já formatado esta equação quadrática como um Bernstein polinomial , pois torna as coisas mais elegantes e é mais numericamente estável. As soluções para s são as raízes dessa equação. Podemos encontrar as raízes utilizando a fórmula quadrática para polinômios Bernstein:
s = ( (A-B) +- sqrt(B^2 - A*C) ) / ( A - 2*B + C )
A fórmula quadrática dá duas respostas devido à + -. Se você está interessado apenas em soluções onde P encontra-se dentro da superfície bilinear então você pode descartar qualquer resposta, onde s não é entre 0 e 1. Para encontrar t, simplesmente volta substituto s em uma das duas equações acima onde nós resolvida para t em termos de s.
Gostaria de salientar um importante caso especial. Se o denominador A - 2 * B + C = 0, então o seu polinomial quadrática é realmente linear. Neste caso, você deve usar uma equação muito mais simples para encontrar s:
s = A / (A-C)
Isto lhe dará exatamente uma solução, a menos que AC = 0. Se A = C, então você tem dois casos: A = C = 0 significa todas valores para s contêm p, caso contrário, não valores para s contêm p.
Outras dicas
Pode-se usar método de Newton para iteratively resolver o seguinte sistema não linear de equações:
p = p0*(1-s)*(1-t) + p1*s*(1-t) + p2*s*t + p3*(1-s)*t.
Note-se que existem duas equações (igualdade em ambas as componentes x e y da equação), e duas incógnitas (S e T).
Para aplicar o método de Newton, precisamos do residual r, que é:
r = p - (p0*(1-s)*(1-t) + p1*s*(1-t) + p2*s*t + p3*(1-s)*t),
e o Jacobiana matriz J, que é encontrado pela diferenciação da residual. Para o nosso problema o Jacobiano é:
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 utilizar o método de Newton, começa-se com uma estimativa inicial (S, T), e, em seguida, executa o seguinte iteração um punhado de vezes:
(s,t) = (s,t) - J^-1 r,
recalcular J e r cada iteração com os novos valores de s e t. Em cada iteração o custo é dominante na aplicação do inverso da Jacobiana para o residual (J ^ -1 r), através da resolução de um sistema linear de 2x2 com J como a matriz de coeficientes e r como do lado direito.
A intuição para o método:
Intuitivamente, se o quadrilátero fosse um paralelogramo , seria muito mais fácil de resolver o problema . O método de Newton é resolver o problema quadrilátero com sucessivas aproximações paralelogramo. A cada iteração nós
-
Uso de informação derivado local no ponto (C, T) para aproximar o quadrilátero com um paralelogramo.
-
Localizar os valores correctos (s, t) sob a aproximação paralelograma, através da resolução de um sistema linear.
-
Ir para este novo ponto e repita.
Vantagens do método:
Como é esperado para os métodos de Newton-tipo, a convergência é extremamente rápido. Como as iterações prosseguir, não só o método chegar mais perto e mais perto do verdadeiro ponto, mas as aproximações paralelogramo locais tornam-se mais preciso, bem como, de modo que a própria taxa de convergência aumenta (no métodos jargão iterativo, dizemos que o método de Newton é quadraticamente convergente ). Na prática isto significa que o número de dígitos corretos aproximadamente dobra a cada iteração.
Aqui está uma tabela representante de iterações vs. erro em um julgamento aleatório eu fiz (veja abaixo para 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)
Depois de duas iterações o erro é pequeno o suficiente para ser bastante eficaz unnoticable e bom para a maioria dos fins práticos, e após 5 iterações o resultado é precisa até os limites de precisão da máquina .
Se você fixar o número de iterações (digamos, 3 iterações para a maioria das aplicações práticas e 8 iterações Se você precisar de muito alta precisão), em seguida, o algoritmo tem uma lógica muito simples e direto a ele, com uma estrutura que se presta bem a computação de alto desempenho. Não há necessidade de verificação de todos os tipos de casos de borda especiais *, e usando lógica diferente, dependendo dos resultados. É apenas um loop que contém algumas fórmulas simples. Abaixo destaco várias vantagens desta abordagem sobre a tradicional fórmula baseada abordagens encontradas em outras respostas aqui e ao redor da Internet:
-
Fácil de código . Basta fazer um loop e digite algumas fórmulas.
-
Não condicionais ou ramificação (se / então), que normalmente permite muito melhor pipelining eficiência .
-
Não há raízes quadradas, e apenas 1 divisão necessária por iteration (se bem escrito). Todas as outras operações são adições simples, subtrações e multiplicações. raízes quadradas e divisões são tipicamente várias vezes mais lento do que adições / multiplicações / multiplicações, e pode estragar cache de eficiência em algumas arquiteturas (mais notavelmente em determinados sistemas embarcados). Na verdade, se você olhar sob o capô em como raízes quadradas e divisões são realmente calculados pelo linguagens de programação modernas, ambos usam variantes do método de Newton, às vezes em hardware e às vezes em software, dependendo do arquitetura.
-
pode ser facilmente vectorized para o trabalho em matrizes com grande número de quadriláteros de uma só vez. Veja meu código vetorizado abaixo um exemplo de como fazer isso.
-
Estende a arbitrária dimensões . O algoritmo estende-se de uma maneira simples para inverter interpolação multilinear em qualquer número de dimensões (2D, 3D, 4D, ...). Eu incluí uma versão 3D abaixo, e pode-se imaginar escrevendo uma versão recursiva simples, ou usando bibliotecas de diferenciação automática para ir para n-dimensões. Método de taxas de convergência dimensão independente tipicamente exposições de Newton, por isso, em princípio, o método deve ser escalável para alguns milhares de dimensões (!) sobre o hardware atual (após o qual ponto de construir e resolver do n-por-n matriz J será, provavelmente, a limitação fator).
É claro que a maioria destas coisas também depende do hardware, compilador, e uma série de outros fatores, assim que sua milhagem pode variar.
Código:
De qualquer forma, aqui está o meu código Matlab: (I liberar tudo aqui para o domínio público)
versão 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
Exemplo de utilização:
% 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)])
Exemplo de saída:
>> test_bilinearInverse
Error after 5 iterations: 1.5701e-16
rápido vetorizado versão 2D:
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 a velocidade este código implicitamente usa a seguinte fórmula para a inversão de uma matriz 2x2:
[a,b;c,d]^-1 = (1/(ad-bc))[d, -b; -c, a]
Exemplo de utilização:
% 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)])
Exemplo de saída:
>> 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
versão 3D:
Inclui algum código para exibir o andamento de convergência.
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
Exemplo de utilização:
%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]))])
Exemplo de saída:
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
Um tem que ter cuidado com a ordenação dos pontos de entrada, uma vez que a interpolação multilinear inverso só está bem definido se a forma tem volume positivo, e em 3D é muito mais fácil de escolher pontos que fazem a forma sua vez dentro para fora de si mesmo.
Aqui está minha implementação da solução da Naaff, como um wiki comunidade. Obrigado mais uma vez.
Esta é uma implementação C, mas deve funcionar em c ++. Ele inclui uma função de teste 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;
}
Uma vez que você está trabalhando em 2D, a sua função bilerp é realmente 2 equações, 1 para x e 1 para y. Eles podem ser reescrita na 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
Onde:
A = p3 - p2 - p1 + p0
B = p2 - p0
C = p1 - p0
D = p0
Reescreva a primeira equação para obter t em termos de s, substituto para o segundo, e resolver para s.
Esta é minha implementação ... Eu acho que é mais rápido do 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;
*/
}
e derivação
// 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 tudo que você tem é um único valor de p, tal que p está entre os valores mínimos e máximos para os quatro cantos da praça, então não, não é possível, em geral, para encontrar uma solução única (s, t ) tal que a interpolação bilinear lhe dará esse valor.
De um modo geral, haverá um número infinito de soluções de (S, T) dentro do quadrado. Eles encontram-se ao longo de uma (hiperbólico) percurso curvo através da praça.
Se a sua função é um vetor valorizado um, então você tem dois valores conhecidos em algum ponto desconhecido na praça? Dado valores conhecidos de dois parâmetros em cada um dos cantos do quadrado, em seguida, uma solução pode existir, mas não há garantia de que. Lembre-se que nós podemos ver isso como dois problemas distintos e independentes. A solução para cada um deles vai estar ao longo de uma linha de contorno hiperbólica através do quadrado. Se o par de contornos cruzam no interior do quadrado, então existe uma solução. Se eles não se cruza, então não existe uma solução.
Você também perguntar se uma fórmula simples existe para resolver o problema. Desculpe, mas não é realmente o que eu vejo. Como eu disse, as curvas são hiperbólica.
Uma solução é mudar para um método diferente de interpolação. Então, ao invés de bilinear, quebrar o quadrado em um par de triângulos. Dentro de cada triângulo, podemos agora usar interpolação verdadeiramente linear. Então agora podemos resolver o sistema linear de 2 equações em 2 incógnitas dentro de cada triângulo. Pode haver uma solução em cada triângulo, excepto para um caso degenerado raras onde os seccionalmente linear contorno linhas correspondentes acontecer para ser co-incidente.
Algumas das respostas ter interpretado mal um pouco a sua pergunta. ie. eles estão supondo que você está dado o valor de uma função interpolada desconhecido, ao invés de uma posição interpolada p (x, y) dentro da quad que você deseja encontrar o (s, t) coordenadas. Este é um problema mais simples e não é garantido para ser uma solução que é a intersecção de duas linhas retas através do quad.
Uma das linhas de corte completamente a p0p1 segmentos e p2p3, o outro corte vai através p0p2 e p1p3, semelhante ao caso alinhado ao eixo. Essas linhas são definidas de forma única pela posição de P (x, y), e, obviamente, intersectam-se neste ponto.
Considerando apenas a linha que corta p0p1 e p2p3, podemos imaginar uma família de tais linhas, para cada s-valor diferente que nós escolhemos, cada corte do quad com uma largura diferente. Se fixarmos uma s-valor, podemos encontrar os dois pontos finais, definindo t = 0 e t = 1.
Então, primeiro assumir a linha tem a forma: y = a0 * x + B0
Então sabemos dois pontos finais desta linha, se fixar um valor dado s. Eles são:
(1-s) P0 + (s) p1
(1-s) p2 + (s) P3
Tendo em conta estes dois pontos finais, podemos determinar a família de linhas, ligando-los na equação para a linha, e resolvendo para a0 e b0 como funções de s . A definição de uma s-valor dá a fórmula para uma linha específica. Tudo o que precisamos agora é descobrir qual linha nesta família atinge nosso ponto p (x, y). Simplesmente conectando as coordenadas de p (x, y) em nossa fórmula linha, podemos resolver para o valor alvo de s.
A abordagem correspondente pode ser feito para encontrar t bem.
Bem, se p é um ponto 2D, sim, você pode obtê-lo facilmente. Nesse caso, o componente S é fraccionada da largura total do quadrilátero em T, T é do mesmo modo o componente fraccionada da altura total do quadrilátero em S.
Se, no entanto, p é um escalar, não é necessariamente possível, porque a função bilinear interpolação não é necessariamente monolítica.
