可以使用牛顿法迭代求解以下非线性方程组:
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,
使用新的 s 和 t 值在每次迭代中重新计算 J 和 r。在每次迭代中,主要成本是将雅可比矩阵的逆应用于残差 (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 维。牛顿方法通常表现出与维度无关的收敛速度,因此原则上该方法应该在当前硬件上可扩展到几千维(!)(在此之后,构建和求解 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 中选择使形状从内向外翻转的点要容易得多。