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

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