反双线性插值？

Question

我有四个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)

这是问题..我有一个2d点p,我知道它在四边形内.我想在使用双线性插值时找到能给我这一点的s,t.

是否有一个简单的公式来反转双线性插值？

---

谢谢你的解决方案.我将我的Naaff解决方案的实现发布为维基.

Answer 1

我认为将问题视为交叉问题是最简单的:参数位置(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叉积(例如,p0 X p1 = x0*y1 - y0*x1).我将这个方程式格式化为二次Bernstein多项式,因为这样可以使事物更加优雅,并且在数值上更稳定.s的解决方案是这个等式的根源.我们可以使用Bernstein多项式的二次公式找到根:

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意味着s的所有值都包含p,否则s的值不包含p.

Answer 2

可以使用牛顿方法迭代求解以下非线性方程组:

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为系数矩阵且r为右手侧的2x2线性系统,将雅可比行列式的逆应用于残差(J ^ -1 r).

该方法的直觉:

直观地说,如果四边形是平行四边形,那么解决问题会容易得多.牛顿方法用连续的平行四边形近似解决四边形问题.在每次迭代我们

1. 在点(s,t)处使用局部导数信息来使用平行四边形来近似四边形.

2. 通过求解线性系统,在平行四边形近似下找到正确的(s,t)值.

3. 跳到这个新点并重复.

方法的优点:

正如牛顿型方法所预期的那样,收敛速度非常快.随着迭代的进行,方法不仅越来越接近真实点,而且局部平行四边形近似也变得更加准确,因此收敛速度本身也增加了(在迭代方法术语中,我们说牛顿方法是二次收敛).实际上,这意味着每次迭代时正确数字的数量大约翻倍.

这是我所做的随机试验的迭代与错误的代表性表格(见下面的代码):

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中,更容易选择使形状从内部转出的点.

Answer 3

这是我对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;
}

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