Python/Scipy 3D 双线性表面映射到单位正方形

Question

我有一个双线性表面（就是由四个顶点之间的线定义的表面），它是在三维空间中的。我想把这个表面映射到一个单位正方形上。我知道怎么进行这种映射，也知道如何根据单位正方形的参数坐标（比如 r 和 s）来计算出一个点的 x、y、z 坐标。不过，我现在需要反过来做。我有一个 x、y、z 的点（我知道它在表面上），我想找出这个点的参数坐标。

我想我可以用 griddata 来实现这个，因为这只是坐标的线性映射。在某些情况下，这个方法非常有效，我得到了我想要的 r 和 s 值。但是在其他情况下，我却收到了 Qhull 的错误信息。我想知道有没有其他方法可以根据 x、y、z 来获取 r 和 s 坐标。

这里有一个有效的示例（按我想要的方式工作）：

from numpy import array
from scipy.interpolate import griddata

P1 = array([[-1.0, -1.0,  1.0],
            [ 1.0, -1.0,  0.0],
            [ 1.0,  1.0,  0.0],
            [-1.0,  1.0,  0.0]])

Mapped_Corners = array([array([0.0, 0.0, 0.0]),
                        array([1.0, 0.0, 0.0]),
                        array([1.0, 1.0, 0.0]),
                        array([0.0, 1.0, 0.0])])

Point_On_Surface = array([[0.5, 0.5, 0.25]])
print griddata(P1, Mapped_Corners, Point_On_Surface, method='linear')

现在，把 P1 和 Point_On_Surface 改成下面的内容，我就会收到 Qhull 的长错误信息：

P1 = array([[-1.0, -1.0,  0.0],
            [ 1.0, -1.0,  0.0],
            [ 1.0,  1.0,  0.0],
            [-1.0,  1.0,  0.0]])
Point_On_Surface = array([[0.5, 0.5, 0.0]])

有没有其他方法可以计算双线性表面上一个点的参数坐标？我注意到在失败的情况下，所有的 'z' 值都是零。我觉得这可能是失败的原因之一，但我想要一个通用的算法，可以把任何表面片（包括平面或所有 z 值为零的表面）映射到单位正方形（z 值为零）。在我的程序中，这些 z 值可以是任何值，包括全为零的情况……

顺便说一下……我确实注意到，如果我把 P1 的最后一列（z=0）和 Point_On_Surface 的 z 坐标去掉，它就能正常工作。但我在想有没有办法在不忽略某些列等特殊情况下来实现这一点。

一般来说，把双线性表面的参数坐标转换为对应的 x、y、z 坐标是一个简单的函数（确保 P1-P3 的方向正确）：

xyz = P0*(1-r)*(1-s) + P1*(r)*(1-s) + P2*(r)*(s) + P3*(1-r)*(s)

其中 P0、P1、P2 和 P3 就是 P1 中的点。

编辑：

在听取了 Paul 和 pv 的建议后，我创建了以下代码，似乎这是个不错的方向。抱歉内容有点长。使用最小二乘法的解决方案比用 fmin 快得多。而且，我通过自己指定雅可比矩阵而不是让 leastsq 来估计它，节省了不少时间。

当按原样运行时，最小化的结果应该等于输入的 'p'（因为指定的 SurfPts 就等于单位正方形本身）。然后你可以取消注释第二个 SurfPts 的定义，以及 'p' 列表中的其他值，以查看其他结果（[-20,15,4] 这个点投影到表面上就是下一个点，-19..... --> 因此 r 和 s 的输出应该是相同的，但偏移的 z 应该不同）。

输出的 [r,s,z] 将是 r,s 系统中的坐标，其中 z 是点相对于表面的法线偏移。这非常方便，因为这个输出提供了很多有用的信息，因为这个函数本质上是“外推”表面（因为它在每个方向上都是线性的）。所以如果 r 和 s 超出了 [0,1] 的范围，那么这个点就会投影到双线性表面之外（所以我会忽略这个情况）。然后你还会得到点相对于表面的法线偏移（z 值）。我会忽略 z，因为我只关心点投影到哪里，但这对某些应用可能会有用。我已经用几个表面和点（包括偏移和不偏移，以及在表面内外的点）测试了这个代码，似乎对我尝试过的所有情况都有效。

代码：

from numpy          import array, cross, set_printoptions
from numpy.linalg   import norm
from scipy.optimize import fmin, leastsq
set_printoptions(precision=6, suppress=True)


#-------------------------------------------------------------------------------
#---[ Shape Functions ]---------------------------------------------------------
#-------------------------------------------------------------------------------
# These are the 4 shape functions for an isoparametric unit square with vertices
# at (0,0), (1,0), (1,1), (0,1)
def N1(r, s): return (1-r)*(1-s)
def N2(r, s): return (  r)*(1-s)
def N3(r, s): return (  r)*(  s)
def N4(r, s): return (1-r)*(  s)


#-------------------------------------------------------------------------------
#---[ Surface Tangent (Derivative) along r and s Parametric Directions ]--------
#-------------------------------------------------------------------------------
# This is the dervative of the unit square with respect to r and s
# The vertices are at (0,0), (1,0), (1,1), (0,1)
def dnT_dr(r, s, Pt):
   return (-Pt[0]*(1-s) + Pt[1]*(1-s) + Pt[2]*(s) - Pt[3]*(s))
def dnT_ds(r, s, Pt):
   return (-Pt[0]*(1-r) - Pt[1]*(r) + Pt[2]*(r) + Pt[3]*(1-r))


#-------------------------------------------------------------------------------
#---[ Jacobian ]----------------------------------------------------------------
#-------------------------------------------------------------------------------
# The Jacobian matrix.  The last row is 1 since the parametric coordinates have
# just r and s, no third dimension for the parametric flat surface
def J(arg, Pt, SurfPoints):
   return array([dnT_dr(arg[0], arg[1], SurfPoints),
                 dnT_ds(arg[0], arg[1], SurfPoints),
                 array([1.,     1.,     1.])])


#-------------------------------------------------------------------------------
#---[ Surface Normal ]----------------------------------------------------------
#-------------------------------------------------------------------------------
# This is the normal vector in x,y,z at any location of r,s
def Norm(r, s, Pt):
   cross_prod = cross(dnT_dr(r, s, Pt), dnT_ds(r, s, Pt))
   return cross_prod / norm(cross_prod)


#-------------------------------------------------------------------------------
#---[ Bilinear Surface Function ]-----------------------------------------------
#-------------------------------------------------------------------------------
# This function converts coordinates in (r, s) to (x, y, z).  When 'normOffset'
# is non-zero, then the point is projected in the direction of the local surface
# normal by that many units (in the x,y,z system)
def nTz(r, s, normOffset, Pt):
   return Pt[0]*N1(r,s) + Pt[1]*N2(r,s) + Pt[2]*N3(r,s) + Pt[3]*N4(r,s) + \
          normOffset*Norm(r, s, Pt)


#-------------------------------------------------------------------------------
#---[ Cost Functions (to minimize) ]--------------------------------------------
#-------------------------------------------------------------------------------
def minFunc(arg, Pt, SurfPoints):
   return norm(nTz(arg[0], arg[1], arg[2], SurfPoints) - Pt)
def lstSqFunc(arg, Pt, SurfPoints):
   return (nTz(arg[0], arg[1], arg[2], SurfPoints) - Pt)


#-------------------------------------------------------------------------------
# ---[ The python script starts here! ]-----------------------------------------
#-------------------------------------------------------------------------------
if __name__ == '__main__':
   # The points defining the surface
   SurfPts = array([array([0.0, 0.0, 0.0]),
                    array([1.0, 0.0, 0.0]),
                    array([1.0, 1.0, 0.0]),
                    array([0.0, 1.0, 0.0])])
   #SurfPts = array([array([-48.62664,        68.346764,  0.3870956   ]),
   #                 array([-56.986549,      -27.516319, -0.70402116  ]),
   #                 array([  0.0080659632,  -32.913471,  0.0068369969]),
   #                 array([ -0.00028359704,  66.750908,  1.6197989   ])])

   # The input point (in x,y,z) where we want the (r,s,z) coordinates of
   p = array([0.1, 1.0, 0.05])
   #p = array([-20., 15.,  4. ])
   #p = array([-19.931894, 15.049718,  0.40244904 ])
   #p = array([20., 15.,  4. ])

   # Make the first guess be the center of the parametric element
   FirstGuess = [0.5, 0.5, 0.0]

   import timeit
   repeats = 100
   def t0():
      return leastsq(lstSqFunc, FirstGuess, args=(p,SurfPts), Dfun=None)[0]
   def t1():
      return leastsq(lstSqFunc, FirstGuess, args=(p,SurfPts), Dfun=J,
                     col_deriv=True)[0]
   def t2():
      return fmin(minFunc, FirstGuess, args=(p,SurfPts), xtol=1.e-6, ftol=1.e-6,
                  disp=False)

   print 'Order:'
   print 'Least Squares, No Jacobian Specified'
   print 'Least Squares, Jacobian Specified'
   print 'Fmin'

   print '\nResults:'
   print t0()
   print t1()
   print t2()

   print '\nTiming for %d runs:' % repeats
   t = timeit.Timer('t0()', 'from __main__ import t0')
   print round(t.timeit(repeats), 6)
   t = timeit.Timer('t1()', 'from __main__ import t1')
   print round(t.timeit(repeats), 6)
   t = timeit.Timer('t2()', 'from __main__ import t2')
   print round(t.timeit(repeats), 6)

另外，如果你不在乎 'z' 的偏移，可以把雅可比矩阵的最后一行设置为全零，而不是现在的全一。这会进一步加快计算速度，r 和 s 的值仍然是正确的，只是你得不到 z 的偏移值。

线性插值雅可比矩阵坐标转换最小二乘法双线性表面三维映射参数坐标 Qhull

Python/Scipy 3D 双线性表面映射到单位正方形

2 个回答

撰写回答