Scipy ODR python

我正试图使用非球面镜片配方将一个9点的云调到圆锥形：

z(r) = r /(R*(1+sqrt(1-(1+K)*(r/R))))

其中R是曲率半径,K是圆锥常数,r = sqrt(xy). K保持不变(已知值),R正是我正在寻找的.我从http://wiki.scipy.org/Cookbook/Least_Squares_Circle开始用python写它.
我用于圆锥曲线的隐式形式是r – 2.R.Z(1 K).Z

这就是我写的：

# -*- coding: cp1252 -*-
from scipy import  odr
from numpy import *

# Coordinates of the 3D points
X   =   [   0,1,-1,0.5,-0.5,0   ]
Y   =   [   0,-0.5    ]
Z   =   [   0,0.113696489,0.027933838,0.027933838]

#constantes
Rc = 8
K = -0.8

def calc_r(x,y):
    return (x**2 + y**2)

def calc_z(r,R):
    return r**2 /(R*(1+sqrt(1-(1+K)*(r**2/R**2))))

def f_3(beta,M):

    r = calc_r(M[0],M[1])
    Z = calc_z(r,beta[0])

    return r**2 - 2*beta[0]*Z + (1+K)*Z**2


beta0 = [Rc]

lsc_data  = odr.Data(row_stack([X,Y]),y=1)
lsc_model = odr.Model(f_3,implicit = True)
lsc_odr   = odr.ODR(lsc_data,lsc_model,beta0)
lsc_out   = lsc_odr.run()

点描述曲率半径为4.5且圆锥常数为-0.8的圆锥曲线.我的代码不起作用：通过ODR,代码返回R = 8(初始点),而不是4.5.知道我的代码有什么问题吗？

谢谢你的帮助

def f_3(beta,M[1])
    Z = M[2]

    return r**2 - 2*beta[0]*Z + (1+K)*Z**2

...
lsc_data  = odr.Data(row_stack([X,Y,Z]),y=1)

这次运行的结果给出R = 4.34911251 – 0.30341252,这似乎符合您的预期.

（编辑：岳阳站长网）

【声明】本站内容均来自网络，其相关言论仅代表作者个人观点，不代表本站立场。若无意侵犯到您的权利，请及时与联系站长删除相关内容!