Можно использовать метод Ньютона , чтобы итеративно решить следующую нелинейную систему уравнений:
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 и введите несколько формул.
Нет условий или ветвлений (если / затем), что обычно позволяет значительно улучшить эффективность конвейерной обработки .
Нет квадратных корней, и только 1 деление требуется на итерацию (если написано правильно). Все остальные операции представляют собой простые сложения, вычитания и умножения. Квадратные корни и деления, как правило, в несколько раз медленнее, чем сложения / умножения / умножения, и могут испортить эффективность кеша в определенных архитектурах (особенно в определенных встроенных системах). В самом деле, если вы посмотрите изнутри на то, как квадратные корни и деления фактически рассчитываются современными языками программирования, они оба используют варианты метода Ньютона, иногда в аппаратном, а иногда в программном обеспечении в зависимости от по архитектуре.
Может быть легко векторизовано для работы с массивами с огромным количеством четырехугольников одновременно. См. Мой векторизованный код ниже для примера того, как это сделать.
Распространяется на произвольные размеры . Алгоритм простирается прямым способом для обратной мультилинейной интерполяции в любом количестве измерений (2d, 3d, 4d, ...). Я включил 3D-версию ниже, и можно представить, как написать простую рекурсивную версию или использовать библиотеки автоматического дифференцирования для перехода к 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'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 гораздо проще выбрать точки, которые заставляют форму вывернуться наизнанку о себе.