Обратная билинейная интерполяция?

Ответ 1

Я думаю, что легче всего думать о вашей проблеме как о проблеме пересечения: что такое местоположение параметра (s, t), где точка p пересекает произвольную двумерную двулинейную поверхность, определяемую p0, p1, p2 и p3.

Подход, который я возьму для решения этой проблемы, похож на предложение tspauld.

Начните с двух уравнений в терминах x и 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 )

Решение для 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) )

Теперь мы можем установить эти два уравнения равными друг другу, чтобы исключить t. Перемещая все в левую часть и упрощая, получим уравнение вида:

A*(1-s)^2 + B*2s(1-s) + C*s^2 = 0

Где:

A = (p0-p) X (p0-p2)
B = ( (p0-p) X (p1-p3) + (p1-p) X (p0-p2) ) / 2
C = (p1-p) X (p1-p3)

Обратите внимание, что я использовал оператор X для обозначения 2D поперечного продукта (например, p0 X p1 = x0 * y1 - y0 * x1). Я отформатировал это уравнение как квадратичное полином Бернштейна, поскольку это делает вещи более элегантными и более численно устойчивыми. Решениями s являются корни этого уравнения. Мы можем найти корни, используя квадратичную формулу для многочленов Бернштейна:

s = ( (A-B) +- sqrt(B^2 - A*C) ) / ( A - 2*B + C )

Квадратичная формула дает два ответа из-за + -. Если вас интересуют только решения, в которых p лежит внутри билинейной поверхности, тогда вы можете отказаться от любого ответа, где s не находится между 0 и 1. Чтобы найти t, просто замените s обратно на одно из двух уравнений выше, где мы решили для t в терминах s.

Я должен указать на один важный частный случай. Если знаменатель A - 2 * B + C = 0, то ваш квадратичный многочлен фактически линейный. В этом случае вы должны использовать гораздо более простое уравнение для поиска s:

s = A / (A-C)

Это даст вам ровно одно решение, если AC = 0. Если A = C, то у вас есть два случая: A = C = 0 означает, что все значения для s содержат p, в противном случае значения s не содержат p.

Ответ 2

Можно использовать метод Ньютона, чтобы итеративно решить следующую нелинейную систему уравнений:

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

Заметим, что существуют два уравнения (равенство как для х, так и для y компонент уравнения) и два неизвестных (s и t).

Чтобы применить метод Ньютона, нам нужен residual r, который:

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

и якобиевая матрица J, которая найдена путем дифференцирования остатка. Для нашей задачи якобиан равен:

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) = (s,t) - J^-1 r,

пересчитывая J и r каждую итерацию с новыми значениями s и t. На каждой итерации доминирующая стоимость заключается в применении обратного к якобиану к остатку (J ^ -1 r) путем решения линейной системы 2x2 с J в качестве матрицы коэффициентов и r в правой части.

Интуиция для метода:

Интуитивно, если четырехугольник был parallelogram, было бы намного легче решить проблему. Метод Ньютона решает проблему четырехугольника с последовательными параллелограммными приближениями. На каждой итерации мы

• Используйте локальную производную информацию в точке (s, t), чтобы приблизить четырехугольник с параллелограммом.

• Найти правильные значения (s, t) в параллелограммном приближении, решая линейную систему.

• Перейти к этой новой точке и повторить.

Преимущества метода:

Как и ожидалось для методов типа Ньютона, сходимость чрезвычайно быстро. По мере продолжения итераций метод не только приближается, но и приближается к истинной точке, но локальные параллелограммные аппроксимации также становятся более точными, поэтому сама скорость сходимости возрастает (в жаргоне итерационных методов мы говорим, что метод Ньютона квадратично конвергентный). На практике это означает, что количество правильных цифр приблизительно удваивается на каждую итерацию.

Здесь представлена ​​репрезентативная таблица итераций против ошибки в случайном пробном произведении (см. ниже для кода):

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

После двух итераций ошибка достаточно мала, чтобы быть практически незаметной и достаточно хорошей для большинства практических целей, а после 5 итераций результат с точностью до пределов точность машины.

Если вы исправите количество итераций (скажем, 3 итерации для большинства практических приложений и 8 итераций, если вам нужна очень высокая точность), то алгоритм имеет очень простую и понятную логику, со структурой, которая хорошо себя зарекомендовала к высокопроизводительному вычислению. Нет необходимости проверять все виды особых случаев края * и использовать различную логику в зависимости от результатов. Это просто цикл for, содержащий несколько простых формул. Ниже я выделил несколько преимуществ этого подхода в отношении традиционных подходов, основанных на формулах, найденных в других ответах здесь и в Интернете:

• Легко кодировать. Просто создайте цикл for и введите несколько формул.

• Нет условностей или ветвлений (если/тогда), что обычно позволяет значительно улучшить эффективность конвейерной обработки .

• Нет квадратных корней, и для каждой итерации требуется только 1 деление (если хорошо написано). Все остальные операции - это простые дополнения, вычитания и умножения. Квадратные корни и деления обычно в несколько раз медленнее, чем дополнения/умножения/умножения, и могут испортить эффективность cache на некоторых архитектурах (в особенности на некоторые встроенные системы). Действительно, если вы посмотрите под капотом на то, как квадратные корни и divisions фактически вычисляются на современных языках программирования, они используют варианты метода Ньютона, иногда в аппаратных средствах и иногда в программном обеспечении в зависимости от архитектуры.

• Может быть легко векторизован для работы с массивами с огромным количеством четырехугольников одновременно. См. Мой векторный код ниже для примера того, как это сделать.

• Допускается до произвольных размеров. Алгоритм распространяется простым способом на инверсную полилинейную интерполяцию в любом числе измерений (2d, 3d, 4d,...). Я включил трехмерную версию ниже, и можно представить себе простую рекурсивную версию или использовать автоматические библиотеки дифференцирования, чтобы перейти к n-измерениям. Метод Ньютона обычно демонстрирует независимые от размерности скорости конвергенции, поэтому в принципе метод должен быть масштабируемым до нескольких тысяч измерений (!) На текущем оборудовании (после чего построение и решение n-на-n-матрицы J, вероятно, будет ограничивающим фактор).

Конечно, большинство из этих вещей также зависит от оборудования, компилятора и множества других факторов, поэтому ваш пробег может меняться.

Код:

В любом случае, здесь мой код Matlab: (я публикую все здесь в общественном достоянии)

Основная 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 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

Пример использования:

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

Пример вывода:

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

Быстрая векторная 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 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

Для скорости этот код неявно использует следующую формулу для инверсии матрицы 2x2:

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

Пример использования:

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

Пример вывода:

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

3D-версия:

Включает некоторый код для отображения прогресса конвергенции.

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

Пример использования:

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

Пример вывода:

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

Нужно быть осторожным при упорядочении входных точек, так как инверсная полилинейная интерполяция только четко определена, если форма имеет положительный объем, а в 3D гораздо проще выбирать точки, которые делают форму вывернутой изнутри само по себе.

Ответ 3

Здесь моя реализация решения Naaff, как сообщество wiki. Еще раз спасибо.

Это реализация C, но она должна работать на С++. Он включает функцию тестирования 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;
}

Ответ 4

Поскольку вы работаете в 2D, ваша функция bilerp - это действительно 2 уравнения, 1 для x и 1 для y. Их можно переписать в форме:

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

Где:

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

Перепишите первое уравнение, получив t через s, замените его на второе и решите для s.

Ответ 5

это моя реализация... Я думаю, что это быстрее, чем 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;
*/ 
}

и вывод

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

Ответ 6

Некоторые ответы слегка неверно истолковали ваш вопрос. то есть. они предполагают, что вам присваивается значение неизвестной интерполированной функции, а не интерполированное положение p (x, y) внутри квада, которое вы хотите найти (s, t). Это более простая проблема, и гарантировано будет решение, являющееся пересечением двух прямых через квад.

Одна из линий будет прорезать сегменты p0p1 и p2p3, другая будет прорезать p0p2 и p1p3, аналогично случаю, выровненному по оси. Эти линии однозначно определяются положением p (x, y) и, очевидно, пересекаются в этой точке.

Рассматривая только линию, которая разрезает через p0p1 и p2p3, мы можем представить семейство таких линий для каждого различного s-значения, которое мы выбираем, каждый разрезающий квадрат с другой шириной. Если мы исправим s-значение, мы можем найти две конечные точки, установив t = 0 и t = 1.

Итак, сначала предположим, что линия имеет вид: y = a0 * x + b0

Тогда мы знаем две конечные точки этой строки, если зафиксировать заданное значение s. Это:

(1-s) p0 + (s) p1

(1-s) p2 + (s) p3

Учитывая эти две конечные точки, мы можем определить семейство линий, подключив их к уравнению для линии и решая для a0 и b0 как функции s. Установка значения s дает формулу для конкретной строки. Теперь нам нужно выяснить, какая линия в этом семействе попадает в нашу точку p (x, y). Просто подключив координаты p (x, y) в нашу формулу линии, мы затем можем решить для целевого значения s.

Соответствующий подход может быть выполнен и для поиска t.

Ответ 7

Если все, что у вас есть, представляет собой одно значение p, такое, что p находится между минимальным и максимальным значениями в четырех углах квадрата, тогда нет, в общем случае невозможно найти одно решение (s, t), так что билинейный интерполятор даст вам это значение.

В общем случае в квадрате будет бесконечное число решений (s, t). Они будут лежать вдоль изогнутого (гиперболического) пути через квадрат.

Если ваша функция является векторнозначной, значит, у вас есть два известных значения в некоторой неизвестной точке квадрата? Учитывая известные значения двух параметров в каждом углу квадрата, тогда решение МОЖЕТ существовать, но об этом нет никакой уверенности. Помните, что мы можем рассматривать это как две отдельные самостоятельные проблемы. Решение каждой из них будет лежать вдоль гиперболической контурной линии через квадрат. Если пара контуров пересекает квадрат, то существует решение. Если они не пересекаются, то решения не существует.

Вы также спрашиваете, существует ли простая формула для решения проблемы. Извините, но на самом деле я этого не вижу. Как я уже сказал, кривые являются гиперболическими.

Одним из решений является переход на другой метод интерполяции. Поэтому вместо билинейного разбивайте квадрат на пару треугольников. Внутри каждого треугольника мы можем теперь использовать по-настоящему линейную интерполяцию. Итак, теперь мы можем решить линейную систему из двух уравнений с двумя неизвестными в каждом треугольнике. В каждом треугольнике может быть одно решение, за исключением редкого вырожденного случая, когда соответствующие кусочно-линейные контурные линии оказываются совпадающими.

Ответ 8

Хорошо, если p - двумерная точка, да, вы можете легко ее получить. В этом случае S является дробной составляющей полной ширины четырехугольника при T, T также является дробной составляющей полной высоты четырехугольника в точке S.

Если, однако, p является скаляром, это необязательно возможно, так как билинейная функция интерполяции не обязательно монолитна.