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定义的任意2D双线性曲面相交。
我将采取的解决此问题的方法与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 cross产品（例如，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替换回上面我们解决的两个方程之一就s而言。
我应该指出一个重要的特例。如果分母A - 2 * B + C = 0，则二次多项式实际上是线性的。在这种情况下，您必须使用更简单的公式来查找s：
s = A / (A-C)
这将为您提供一个解决方案，除非AC = 0.如果A = C，那么您有两种情况：A = C = 0意味着所有值包含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,
每次迭代时使用 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 循环并输入一些公式即可。
没有条件或分支 （如果/那么），这通常可以带来更好的结果 流水线效率.
没有平方根，只有 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 中，选择使形状本身从内到外翻转的点要容易得多。
这是我对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
根据 s 重写第一个等式得到 t ，替换为第二个，并求解 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）。它们将沿着弯曲的（双曲线）路径穿过广场。
如果你的函数是一个值为1的向量，那么你在广场的某个未知点有两个已知值？给定广场每个角落的两个参数的已知值，则可能存在解决方案，但无法保证。请记住，我们可以将其视为两个独立的独立问题。它们中的每一个的解决方案将沿着通过正方形的双曲线轮廓线。如果这对轮廓在正方形内交叉，则存在解决方案。如果他们没有交叉，那么就没有解决方案。
您还会问是否存在一个简单的公式来解决问题。对不起，但不是我看到的。正如我所说，曲线是双曲线的。
一种解决方案是切换到不同的插值方法。因此，不是双线性，而是将正方形分成一对三角形。在每个三角形内，我们现在可以使用真正的线性插值。所以现在我们可以解决每个三角形内2个未知数的2个方程的线性系统。每个三角形可能有一个解决方案，除了罕见的退化情况，其中相应的分段线性轮廓线碰巧是共同发生的。
有些回复略微误解了您的问题。即。他们假设给你一个未知插值函数的值，而不是你想要找到（s，t）坐标的四元组内的插值位置p（x，y）。这是一个更简单的问题，并且保证有一个解决方案是通过四边形的两条直线的交叉点。
其中一条线将穿过区段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是2D点，是的，你可以很容易地得到它。在这种情况下，S是T处四边形总宽度的分数分量，T同样是S处四边形总高度的分数分量。
但是，如果p是标量，则不一定可能，因为双线性插值函数不一定是单片的。
