Odwrotna Dwuliniowa Interpolacja?

Można by użyć Metoda Newtona aby iteracyjnie rozwiązać następujący nieliniowy układ równań:

p = p0*(1-s)*(1-t) + p1*s*(1-t) + p2*s*t + p3*(1-s)*t.

Zauważ, że istnieją dwa równania (równość w składowych x i y równania) i dwie niewiadome (s I t).

Aby zastosować metodę Newtona, potrzebujemy resztek r, czyli:

r = p - (p0*(1-s)*(1-t) + p1*s*(1-t) + p2*s*t + p3*(1-s)*t),

I macierz Jakobiańska J, którą można znaleźć poprzez różnicowanie szczątków. Dla naszego problemu Jakobian jest:

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).

(s,t) = (s,t) - J^-1 r,

Przeliczanie J i R każdej iteracji z nowymi wartościami s I t. przy każdej iteracji dominującym kosztem jest zastosowanie odwrotności Jakobiańskiej do resztkowej (j^-1 r), rozwiązując układ liniowy 2x2 z J jako macierzą współczynnika i R jako prawą stroną.

Intuicja dla metoda:

Zgodnie z oczekiwaniami dla metod typu Newtona, zbieżność jest niezwykle szybka. W miarę postępu iteracji, nie tylko metoda zbliża się coraz bliżej do prawdziwego punktu, ale także lokalne przybliżenia równoległoboku stają się bardziej dokładne, więc sama szybkość konwergencji wzrasta (w żargonie metod iteracyjnych mówimy, że Metoda Newtona na zbieżność kwadratowa). W praktyce oznacza to, że liczba poprawnych cyfr w przybliżeniu podwaja każdą iterację.

Iteration  Error
1          0.0610
2          9.8914e-04
3          2.6872e-07
4          1.9810e-14
5          5.5511e-17 (machine epsilon)

Po dwóch iteracjach błąd jest na tyle mały, że jest skutecznie niezauważalny i wystarczająco dobry dla najbardziej praktycznych celów, a po 5 iteracjach wynik jest dokładny do granic maszyny precyzja .

Jeśli ustalisz liczbę iteracji (powiedzmy 3 iteracji dla większości praktycznych zastosowań i 8 iteracji, jeśli potrzebujesz bardzo wysokiej dokładności), algorytm ma bardzo prostą i prostą logikę, ze strukturą, która dobrze nadaje się do wysokowydajnych obliczeń. Nie ma potrzeby sprawdzania różnego rodzaju specjalnych przypadków krawędzi * i korzystania z innej logiki w zależności od wyników. Jest to tylko pętla for zawierająca kilka prostych formuł. Poniżej I wyróżnij kilka zalet tego podejścia w stosunku do tradycyjnych metod opartych na formułach znalezionych w innych odpowiedziach tutaj i wokół Internetu: {]}

• łatwy do kodowania . Wystarczy wykonać pętlę for i wpisać kilka formuł.

• brak uwarunkowań lub rozgałęzień (if / then), co zazwyczaj pozwala na znacznie lepszą wydajność pipelining.

• brak pierwiastków kwadratowych i tylko 1 Podział wymagany na iterację (jeśli dobrze napisane). Wszystkie pozostałe operacje to proste dodawanie, odejmowanie i mnożenie. Pierwiastki kwadratowe i podziały są zazwyczaj kilka razy wolniejsze niż dodawanie / mnożenie/mnożenie i mogą zepsuć wydajność pamięci podręcznej {19]} na niektórych architekturach (zwłaszcza na niektórych systemach wbudowanych). W rzeczywistości, jeśli spojrzysz pod maskę, jak pierwiastki kwadratowe i dywizje są faktycznie obliczane przez współczesne języki programowania, oba używają wariantów Metoda Newtona, czasami w sprzęcie, a czasami w oprogramowaniu w zależności od architektury.

• może być łatwo wektoryzowany do pracy na tablicach z ogromną liczbą czworokątów na raz. Zobacz mój wektoryzowany kod poniżej, aby zobaczyć przykład, jak to zrobić.

• rozszerza do dowolnych wymiarów . Algorytm rozszerza się w prosty sposób na odwrotną interpolację wieloliniową w dowolnej liczbie wymiarów (2d, 3d, 4d, ...). Włączyłem wersja 3D poniżej i można sobie wyobrazić napisanie prostej rekurencyjnej wersji, lub użycie automatycznych bibliotek różnicujących, aby przejść do n-wymiarów. Metoda Newtona zazwyczaj wykazuje współczynniki zbieżności niezależne od wymiarów, więc w zasadzie metoda powinna być skalowalna do kilku tysięcy wymiarów (!) na aktualnym sprzęcie (po którym to punkcie konstruowanie i rozwiązywanie N-PO-N macierzy J będzie prawdopodobnie czynnikiem ograniczającym).

Oczywiście większość z tych rzeczy również zależy na sprzęcie, kompilatorze i wielu innych czynnikach, więc przebieg może się różnić.

Kod:

Podstawowa wersja 2D:

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

Przykładowe użycie:

% 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)])

Przykładowe wyjście:

>> test_bilinearInverse
Error after 5 iterations: 1.5701e-16

Szybka wektorowa wersja 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

Dla prędkości ten kod domyślnie używa następującego wzoru na odwrotność 2x2 macierz:

[a,b;c,d]^-1 = (1/(ad-bc))[d, -b; -c, a]

Przykładowe użycie:

% 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)])

Przykładowe wyjście:

>> 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

Wersja 3D:

Zawiera kod wyświetlający postęp konwergencji.

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

Przykładowe użycie:

%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]))])

Przykładowe wyjście:

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

Należy uważać na kolejność punktów wejściowych, ponieważ odwrotna interpolacja wieloliniowa jest dobrze zdefiniowana tylko wtedy, gdy kształt ma dodatnią objętość, a w 3D o wiele łatwiej jest wybrać punkty to sprawia, że kształt obraca się wewnątrz siebie.

Oto moja implementacja rozwiązania Naaff, jako wiki społeczności. Jeszcze raz dziękuję.

Jest to implementacja C, ale powinna działać na c++. Zawiera funkcję testowania 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;
}

Ponieważ pracujesz w 2D, twoja bilerp Funkcja to naprawdę 2 równania, 1 dla x i 1 dla y. można je przepisać w postaci:

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

Gdzie:

A = p3 - p2 - p1 + p0
B = p2 - p0
C = p1 - p0
D = p0

Przepisz pierwsze równanie, aby uzyskać t w kategoriach s, zamień na drugie i rozwiąż dla s.

To moja realizacja ... Myślę, że jest szybszy niż z 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;
*/ 
}

I wyprowadzenie

// 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

Niektóre odpowiedzi nieco błędnie zinterpretowały twoje pytanie. ie. zakładają, że otrzymujesz wartość nieznanej interpolowanej funkcji, a nie interpolowaną pozycję P (x,y) wewnątrz kwadratu,którego chcesz znaleźć współrzędne (s, t). Jest to prostszy problem i gwarantowane jest rozwiązanie, które jest przecięciem dwóch linii prostych przez kwadrat.

Jedna z linii będzie przecinać segmenty p0p1 i p2p3, druga będzie przecinać przez p0p2 i p1p3, podobnie jak w przypadku osi. Linie te są jednoznacznie określone przez położenie p (x, y)i oczywiście przecinają się w tym punkcie.

Biorąc pod uwagę tylko linię przecinającą p0p1 i p2p3, możemy sobie wyobrazić rodzinę takich linii, dla każdej wybranej przez nas wartości s, każda przecinająca kwadrat na innej szerokości. Jeśli naprawimy wartość s, możemy znaleźć dwa punkty końcowe, ustawiając T = 0 i T = 1.

Więc najpierw Załóżmy, że linia ma postać: y = a0 * x + b0

Wtedy znamy dwa punkty końcowe tej linii, jeśli ustalimy daną wartość s. Są to:

(1-s)p0 + (s)p1

(1-s)p2 + (s)p3

Biorąc pod uwagę te dwa punkty końcowe, możemy wyznaczyć rodzinę linii poprzez włączenie ich do równania dla linii i rozwiązanie dla A0 i b0 jako funkcji s. Ustawienie wartości s daje formułę dla określonej linii. Wszystko, czego teraz potrzebujemy, to dowiedzieć się, która linia w tej rodzinie uderza w nasz punkt p (x, y). Wystarczy podłączyć współrzędne P (x, y) do naszego wzoru liniowego, możemy wtedy rozwiązać dla wartości docelowej s.

Odpowiednie podejście może być wykonane, aby znaleźć t, jak również.

Jeśli masz tylko jedną wartość p, taką, że p znajduje się pomiędzy minimalnymi i maksymalnymi wartościami na czterech rogach kwadratu, to nie, ogólnie nie jest możliwe znalezienie pojedynczego rozwiązania (s,t) takiego, że interpolant dwuliniowy da ci tę wartość.

Ogólnie rzecz biorąc,wewnątrz kwadratu będzie nieskończona liczba rozwiązań (s, t). Będą leżeć wzdłuż zakrzywionej (hiperbolicznej) ścieżki przez kwadrat.

Jeśli twoja funkcja jest wektorem o wartości 1, więc masz dwa znane wartości w jakimś nieznanym punkcie na placu? Biorąc pod uwagę znane wartości dwóch parametrów na każdym rogu kwadratu, może istnieć rozwiązanie, ale nie ma pewności co do tego. Pamiętajmy, że możemy to postrzegać jako dwa oddzielne, niezależne problemy. Rozwiązanie każdego z nich będzie leżeć wzdłuż hiperbolicznej linii konturu przez kwadrat. Jeśli para konturów krzyżuje się wewnątrz kwadratu, to istnieje rozwiązanie. Jeśli się nie krzyżują, to nie ma rozwiązania.

Pytasz też, czy prosty istnieje formuła pozwalająca rozwiązać problem. Przykro mi, ale nie bardzo to widzę. Jak mówiłem, krzywe są hiperboliczne.

Jednym z rozwiązań jest przejście na inną metodę interpolacji. Więc zamiast bilinear, rozbić kwadrat na parę trójkątów. W każdym trójkącie możemy teraz użyć prawdziwie liniowej interpolacji. Więc teraz możemy rozwiązać układ liniowy 2 równań w 2 niewiadomych w każdym trójkącie. W każdym trójkącie może być jedno rozwiązanie, z wyjątkiem rzadkiego przypadku degeneratu, w którym odpowiednie liniowe linie konturu są współwystępujące.

Cóż, jeśli p jest punktem 2D, tak, możesz go łatwo zdobyć. W tym przypadku S jest składnikiem ułamkowym całkowitej szerokości czworokąta w T, T jest również składnikiem ułamkowym całkowitej wysokości czworokąta w s.

Jeśli jednak p jest skalarem, to niekoniecznie jest to możliwe, ponieważ Dwuliniowa funkcja interpolacyjna niekoniecznie jest monolityczna.