Python/Scipy 3D Bilinear Surface Mapped to Unit Square

I have a bilinear surface (surface defined by lines between the 4 vertices) in 3D. I want to map this surface to the unit square. I know how to do this mapping and I know how to calculate the x,y,z coordinates of a point given the parametric coordinates of the unit square (say r,s). However, I need to go the other way. I have the x,y,z point (that I know lies on the surface) and I want to determine the parametric coordinate of this point.

I thought I could use griddata to do this since it is just a linear mapping of the coordinates. For some cases, this works perfectly and I get the correct r,s values that I'm looking for. But in other cases, I get an error returned from Qhull. I want to know if there is another way to get the r,s coordinates given x,y,z.

Here is a working example (works as I want it to):

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')

Now, change P1 and Point_On_Surface to the following and I get a long error from 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]])

Is there some other way to calculate the parametric coordinates of a point on a bilinear surface? I realize that in the failed case, all 'z' values are zero. I think this is part of the reason why it fails, but I want a generic algorithm that would map any surface patch (including a planar one or one with all z values as zero) to the unit square (with z values of zero). In my program, those z values could be anything, including all zero...

Just as an FYI...I did notice that it works correctly if I take the last column (z=0) from P1 out as well as the z coordinate from Point_On_Surface. But I'm wondering about doing this without special cases were I ignore some column, etc.

In general, it is a simple function to convert the parametric coordinates of a bilinear surface to their corresponding x,y,z coordinates (with the correct orientation of P1-P3):

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

Where P0, P1, P2, and P3 would be the points in P1.

EDIT:

After taking the advice from Paul and pv, I've created the following code which seems to be the way to go. Sorry for the long length. The solution with the least squares is MUCH, MUCH quicker then using fmin. Also, I gain a significant time improvement by specifying my Jacobian myself instead of letting leastsq estimate it.

When run as-is, the result from the minimization should equal the input 'p' for any input (since the SurfPts specified are equal to the the unit square itself). Then you can un-comment the second definition of SurfPts and any of the other values listed for 'p' to see other results (the [-20,15,4] point projected onto the surface is that next point, -19..... --> thus the r,s output should be the same, but the offset z should be different).

The output [r,s,z] will be the coordinate in the r,s system with z being the offset of the point normal to the surface. This is really convenient because the output of this gives loads of useful information since the function essentially "extrapolates" the surface (since it is linear in each direction). So if r and s are outside of the range [0,1], then the point gets projected outside of the bilinear surface (so I'll ignore if for my purposes). Then you also get the offset of the point normal to the surface (the z value). I will ignore z since I just care about where the point projects to, but this could be useful for some applications. I've tested this code with several surfaces and points (offset and not offset as well as outside and inside the surface) and it seems to work for everything I've tried so far.

CODE:

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)

Also, if you don't care about the 'z' offset, you can set the last row of the Jacobian matrix to all zeros instead of ones as it is now. This will speed up the calculations even more and the values for r and s are still correct, you just don't get the z offset value.

1) Griddata cannot do this, since it only does volume interpolation. Surface has zero volume, so this use case is out of scope for griddata, hence the errors. Moreover, since floating point has only finite precision, usually no point lies exactly on the surface.

Your initial attempt with griddata also cannot work correctly in any case: griddata is triangle-based, and cannot produce correct maps for bilinear surfaces.

Do you really need the surface to be bilinear? If yes, there is in general no way to avoid solving a nonlinear minimization problem to find the r, s coordinates, as the coordinate map itself is nonlinear (you multiply r with s in it). Such problems can be solved e.g. with scipy.optimize.leastsq:

# untested code, but gives the idea
def func(c):
    r, s = c
    xyz = P0*(1-r)*(1-s) + P1*(r)*(1-s) + P2*(r)*(s) + P3*(1-r)*(s)
    return xyz - Point_on_Surface

res = scipy.optimize.leastsq(func, [0, 0])
r, s = res[0]

2) If you are happy with a triangular surface where each triangle is flat, then you can solve the problem with linear projections onto the triangles.

So let's assume your surface consists of flat triangles. The problem splits into two subproblems: (A) finding the projection of a 3D point onto a given triangle, and (B) finding the triangle inside which the point projects to.

Given corners A, B, C of a triangle (in 3D), any point inside the triangle can be written in the form P = c_1 A + c_2 B + (1 - c_1 - c_2) C. c_k specify the barycentric coordinate system on the triangle. Given a point P in 3D, we can then find the barycentric coordinates of a point closest to it on the triangle via c = np.linalg.lstsq([A-C, B-C], P-C). So, with V[0], V[1], V[2] = A, B, C,

def project_onto_triangle(P, V):
    c, distance, rank, s = np.linalg.lstsq(np.c_[V[0] - V[2], V[1] - V[2]], P - V[2])
    return c, distance

Point_on_Surface = 0.25 * P1[0] + 0.25 * P1[1] + 0.5 * P1[2]
c, dist = project_onto_triangle(Point_on_Surface, P1[0:3])
# >>> c
# array([ 0.25,  0.25])      # barycentric coordinates
# >>> dist
# array([  4.07457566e-33])  # the distance of the point from the triangle plane

# second triangle
c2, dist2 = project_onto_triangle(Point_on_Surface, [P1[0], P1[2], P1[3]])
# >>> c2
# array([ 0.45,  0.75])
# >>> dist2
# array([ 0.05])

If you need to do this repeatedly, note that np.linalg.lstsq(V, P - C) == Q.dot(P - C) if Q = np.linalg.pinv(V) --- so you can precompute the projector matrixes Q. The c coordinates obtained then then tell you the coordinates on the unit square.

You will need to go through all the triangles to find the one in which the point lies. A point is on the triangle if both of the following are true:

• dist < epsilon, where epsilon is some small tolerance (point near the triangle plane)

• 0 <= c_1 <= 1, 0 <= c_2 <= 1 and 0 <= 1 - c_1 - c_2 <= 1 (point projection inside triangle)

In the above example case, we find that the point is in the first triangle. For he second triangle, the distance is small (0.05), but the barycentric coordinates indicate that the projection lies outside the triangle.

There are some optimizations one can make in the above code if many points need to be projected at once etc.

If this is still too slow, the problem gets harder. griddata uses special trick for finding the correct simplex, but these work only for volume grids. I don't know a fast algorithm in this case, although there might be something suitable in the computational geometry literature.

This may answer your question. I described the issues I have with the code in a comment. This describes a method of getting coordinates which is an attempt to answer the text of the question, ignoring the code.

This is just one method to find a point on an arbitrary surface at height z_target. There are much more efficient ways to do this if, for example, you know the surface is piecewise planar, or monotonically increasing or whatever.. This is kind of overkill for a situation where you are intentionally generating the surface

There are, of course, infinitely many points on the surface that meet z_target (Think contour lines). Finding just one arbitrary point may not be of much value.

Anyway, this will always seek the x,y pair of coords that meets a target surface height. There will be more than one coordinate that meets the target, but this will find one, if it exists.

from numpy import array
from scipy.interpolate import LinearNDInterpolator
from scipy.optimize import fmin

# X,Y points (unstructured)
points = array([
    [-1.0, -1.0],
    [ 1.0, -1.0],
    [ 1.0,  1.0],
    [-1.0,  1.0]])

# Z - height of surface at points given above
Z = array([5.,13.,23.,4.0])

# interpolant (function that interpoates)
# this function will return the interpolated z value for any given point
# this is the interpolant being used for griddata(method='linear')
# check numpy.source(griddata) to see for yourself.
interp = LinearNDInterpolator(points, Z)

# this is the desired z value for which we want to know the coordinates
# there is no guarantee that there is only one pair of points that
# returns this value.
z_target = 7.0

# so we select an initial guess for the coordinates.
# we will make a root finding function to find a point near this one
point0 = array((0.5, 0.5))

# this is the cost function.  It tells us how far we are from our target
# ideally we want this function to return zero
def f(arr):
    return (interp(*arr) - z_target)**2

# run the downhill simplex algorithm to root-find the zero of our cost function.
# the result should be the coords of a point who's height is z_target
# there are other solvers besides fmin
print fmin(f, x0)