Обратная билинейная интерполяция?
https://stackoverflow.com/questions/808441
-
03-07-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianВопрос
У меня есть четыре 2d точки, p0 = (x0, y0), p1 = (x1, y1) и т.д.которые образуют четырехугольник.В моем случае четырехугольник не прямоугольный, но он должен быть, по крайней мере, выпуклым.
p2 --- p3
| |
t | p |
| |
p0 --- p1
s
Я использую билинейную интерполяцию.S и T находятся в пределах [0..1], и интерполированная точка задается:
bilerp(s,t) = t*(s*p3+(1-s)*p2) + (1-t)*(s*p1+(1-s)*p0)
Вот в чем проблема..У меня есть двумерная точка p, которая, как я знаю, находится внутри квадрата.Я хочу найти s, t, которые дадут мне эту точку при использовании билинейной интерполяции.
Существует ли простая формула для обратной билинейной интерполяции?
Спасибо за решения.Я опубликовал свою реализацию решения Naaff в виде wiki.
Решение
Я думаю, проще всего рассматривать вашу проблему как проблему пересечения:каково местоположение параметра (s, t), где точка p пересекает произвольную двумерную билинейную поверхность, определяемую p0, p1, p2 и p3.
Подход, который я применю к решению этой проблемы, аналогичен предложению тсполда.
Начните с двух уравнений в терминах 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)
Это даст вам ровно одно решение, если только A-C = 0.Если A = C, то у вас есть два случая:A=C=0 означает ВСЕ значения для s содержат p, в противном случае НЕТ значения для s содержат p.
Другие советы
Можно было бы использовать Метод Ньютона итеративно решить следующую нелинейную систему уравнений:
p = p0*(1-s)*(1-t) + p1*s*(1-t) + p2*s*t + p3*(1-s)*t.
Обратите внимание, что существует два уравнения (равенство как в x, так и в y компонентах уравнения) и два неизвестных (s и t).
Чтобы применить метод Ньютона, нам нужен остаточный 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 в качестве правой части.
Интуиция для этого метода:
Интуитивно, если бы четырехугольник был параллелограмм, было бы гораздо проще решить эту проблему.Метод Ньютона заключается в решении четырехугольной задачи последовательными параллелограммными приближениями.На каждой итерации мы
Используйте информацию о локальной производной в точке (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 и введите несколько формул.
Никаких условностей или ветвлений (if / then), что обычно позволяет гораздо лучше эффективность конвейеризации.
Никаких квадратных корней и только 1 деление требуется для каждой итерации (если написано хорошо).Все остальные операции - это простые сложения, вычитания и умножения.Квадратные корни и деления обычно выполняются в несколько раз медленнее, чем сложения / умножения / умножения, и могут привести к сбоям Кэш эффективность на определенных архитектурах (в первую очередь на определенных встроенных системах).Действительно, если вы заглянете под капот на то, как квадратные корни и подразделения фактически рассчитываются современными языками программирования, они оба используют варианты метода Ньютона, иногда аппаратно, а иногда и программно, в зависимости от архитектуры.
Может быть легко векторизован работать с массивами с огромным количеством четырехугольников одновременно.Смотрите мой векторизованный код ниже для примера того, как это сделать.
Распространяется на произвольные размеры.Алгоритм простым способом распространяется на обратную многолинейную интерполяцию в любом количестве измерений (2d, 3d, 4d, ...).Я включил 3D-версию ниже, и можно представить себе написание простой рекурсивной версии или использование библиотек автоматического дифференцирования для перехода к n-измерениям.Метод Ньютона обычно демонстрирует не зависящие от размера скорости сходимости, поэтому в принципе метод должен быть масштабируемым до нескольких тысяч измерений (!) на текущем оборудовании (после чего построение и решение матрицы J n-by-n, вероятно, будет ограничивающим фактором).
Конечно, большинство из этих факторов также зависит от аппаратного обеспечения, компилятора и множества других факторов, поэтому ваш пробег может варьироваться.
Код:
В любом случае, вот мой код 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'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
Пример использования:
% 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'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
Для повышения скорости этот код неявно использует следующую формулу для инверсии матрицы 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 гораздо проще выбрать точки, которые заставляют фигуру выворачиваться наизнанку.
Вот моя реализация решения Naaff в качестве вики сообщества.Еще раз спасибо.
Это реализация на C, но должна работать на c ++.Он включает в себя функцию тестирования нечеткости.
#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;
}
Поскольку вы работаете в 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.
это моя реализация ...Я думаю, это быстрее, чем в 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
Если все, что у вас есть, это единственное значение p, такое, что p находится между минимальным и максимальным значениями в четырех углах квадрата, то нет, вообще невозможно найти ЕДИНСТВЕННОЕ решение (s, t), такое, что билинейный интерполятор даст вам это значение.
В общем случае внутри квадрата будет бесконечное число решений (s, t).Они будут лежать вдоль изогнутой (гиперболической) траектории через квадрат.
Если ваша функция имеет векторное значение, значит, у вас есть два известных значения в какой-то неизвестной точке квадрата?Учитывая известные значения двух параметров в каждом углу квадрата, решение МОЖЕТ существовать, но уверенности в этом нет.Помните, что мы можем рассматривать это как две отдельные, независимые проблемы.Решение каждого из них будет лежать вдоль гиперболической контурной линии, проходящей через квадрат.Если пара контуров пересекается внутри квадрата, то решение существует.Если они не пересекаются, то никакого решения не существует.
Вы также спрашиваете, существует ли простая формула для решения проблемы.Извините, но на самом деле я ничего такого не вижу.Как я уже сказал, кривые являются гиперболическими.
Одним из решений является переход на другой метод интерполяции.Поэтому вместо билинейности разбейте квадрат на пару треугольников.Внутри каждого треугольника теперь мы можем использовать действительно линейную интерполяцию.Итак, теперь мы можем решить линейную систему из 2 уравнений с 2 неизвестными внутри каждого треугольника.В каждом треугольнике может быть одно решение, за исключением редкого вырожденного случая, когда соответствующие кусочно-линейные контурные линии совпадают.
Некоторые из ответов немного неверно истолковали ваш вопрос.т.е.они предполагают, что вам даны значение неизвестной интерполированной функции, а не интерполированной позиции 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 .
Ну, если p - двумерная точка, то да, вы можете легко ее получить.В этом случае S является дробной составляющей общей ширины четырехугольника в точке T, T также является дробной составляющей общей высоты четырехугольника в точке S.
Если, однако, p является скаляром, это не обязательно возможно, поскольку функция билинейной интерполяции не обязательно монолитна.
