用Python求解射线与曲面交点

用Python求解射线与曲面交点

在光学射线追迹问题中，往往需要计算光线与光学元件表面的交点 如果是球面或者有限阶数的非球面，可以找到明确的解析解，（5阶以下） 但是如果是面对高阶非球面以及一些更加复杂的自由曲面时，就只能通过插值近似求解。 一般情况下，这样求解的话可以通过二分法来求 ————先找到射线位于表面的两边的两个端点，然后设置精度并二分查找，就可以得到交点

定义射线

每一个射线可以用两个点来定义，但也可以用点以及向量来定义这里选择用点与向量——主要时为了后期计算折射反射方便。

定义曲面

如果不考虑后期计算可以直接用二元函数z = f(x,y)，但是由于后期要进行其他计算，因此需要将表达式封装到类中，但是这里直接用了表达式 例：

def face_fun(val): return val[0]**2+val[1]**2 + 20

根据以上代码来说具体运算如下

class GetNodicalLineFaceFreedom: def __init__(self,rp,direct,face_fun,acu): """ rp:参考点，direct：直线方向，face_fun：曲面表达式，acu：精度 """ self.rp = rp self.direct=direct self.face_fun = face_fun self.acu = acu def get_range(self): """ 得到两个端点 """ if self.face_fun(self.rp[:2]): afirm_value = self.face_fun(self.rp[:2]) step = self.direct*(afirm_value-self.rp[2]) point =self.rp+step start = self.rp mun = (self.face_fun(point[:2])-point[2])* (self.face_fun(start[:2])-start[2]) while mun>0:# 寻找第二个端点位于曲面另一边 point += step mun = (self.face_fun(point[:2])-point[2])* (self.face_fun(start[:2])-start[2]) return start,point def get_nodical(self): """ 求取交点 """ start,tail = self.get_range() while self.norm(start,tail) > self.acu:#两个点的距离。 newpoint = (start+tail)/2 if (self.face_fun(newpoint[:2])-newpoint[2])*(self.face_fun(start[:2])-start[2])>0: start = newpoint elif (self.face_fun(newpoint[:2])-newpoint[2])*(self.face_fun(start[:2])-start[2])<0: tail = newpoint else: return newpoint return start def norm(self,point1,point2): """ 两个点的距离 """ dis = point1-point2 return np.sum(np.square(dis))

二分的好处就在于用非常有限的迭代（N = log₂(length/acu)）得到很高的精度

import numpy as np rp = np.array([4,2,13]) direct = np.array([1,2,9]) acu = 0.0001 nodi = GetNodicalLineFaceFreedom(rp,direct,face_fun,acu)

当然如果不是实在得不到解析解，尽量不要使用这种方法，因为光线追迹以及镜头优化有着惊人的计算量。 所以一般都需要将表面进行分类，尽量寻找解析解。