У меня есть билинейная поверхность (поверхность, определяемая линиями между 4 вершинами) в 3D.Я хочу отобразить эту поверхность на единицу площади.Я знаю, как сделать это отображение, и я знаю, как рассчитать координаты x, y, z точки с учетом параметрических координат единичного квадрата (скажем, r, s).Однако мне нужно пойти другим путем.У меня есть точки x, y, z (которые, как я знаю, лежит на поверхности), и я хочу определить параметрическую координату этой точки.
Я подумал, что могу использовать griddata, чтобы сделать это, так как это всего лишьлинейное отображение координат.В некоторых случаях это работает отлично, и я получаю правильные значения r, которые я ищу.Но в других случаях я получаю сообщение об ошибке от Qhull.Я хочу знать, есть ли другой способ получить координаты r, s, заданные x, y, z.
Вот рабочий пример (работает так, как я хочу):
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 могут быть любыми, включая все нули ...
Так же, как к вашему сведению ... Я заметил, что это работает правильно, если я беру последний столбец (z = 0) из P1вне, а также координата z от Point_On_Surface.Но мне интересно делать это без особых случаев, если я игнорирую некоторые столбцы и т. Д.
В общем, это простая функция для преобразования параметрических координат билинейной поверхности в соответствующие им 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.
РЕДАКТИРОВАТЬ:
После того, как я принял совет Пола и П.В., я создал следующий код, который, кажется, подходит.Извините за длинную длину.Решение с наименьшими квадратами НАМНОГО, НАМНОГО быстрее, чем использование 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.