python光學仿真實現光線追跡之空間關系
空間關系
變化始於相遇，所以交點是一切的核心。
相交判定
首先考察一束光線能否打在某個平面鏡上。光線被抽象成瞭一個列表[a,b,c]，平面鏡則被抽象成為由兩個點構成的線段[(x1,y1),(x2,y2)]。兩條直線的交點問題屬於初等數學范疇，需要先將線段轉換成直線的形式，然後再求交點。但是兩條直線的交點可能落在線段的外面，從而不具有判定的意義。
如果我們的光學系統中有大量的光學元件，那麼如果有一種方法可以快速判斷光線是否與光學元件有交點，將會顯得更加快捷。那麼，如果一個直線穿過某個線段，就意味著這條線段的兩個端點必然在直線的兩側。
import numpy as np
def isCross(abc,dots):
abc = np.array(abc) #將abc的格式轉成數組
flag1 = abc.dot(list(dots[0])+[1]) #數組的點積
flag2 = abc.dot(list(dots[1])+[1])
return flag1*flag2
我們非常熟悉運算符”+”，不過目前來說，隻有數值和數組是支持加法運算的。所以，對於list(dots[0])+[1]這種表示著實讓人有些摸不到頭腦。
這個含義其實是符合人類直覺的。列表內的元素個數是可變的，兩個列表相加可以理解為兩個列表銜接在一起。當然，元組並不支持這種運算。
例如
>>> [1,2,3]+[4] [1, 2, 3, 4] >>> (1,2,3)+(4) Traceback (most recent call last): File "<stdin>", line 1, in <module> TypeError: can only concatenate tuple (not "int") to tuple
通過加法運算符來銜接兩個列表，實際上相當於新建瞭一個變量，需要開辟新的內存空間。好在對於初學者來說這樣不容易出錯。
在numpy中，+、-、*、/這幾個運算符表示對應位置元素的運算。如果想使用點乘等其他運算，需要調用numpy中的其他函數。
>>> np.array([1,2,3])*np.array([4,5,6])
array([ 4, 10, 18])
>>> np.array([1,2,3])+np.array([4,5,6])
array([5, 7, 9])
>>> np.array([1,2,3]).dot([4,5,6]) #.dot表點積
32
所以， f l a g 1 = [ a , b , c ] ⋅ [ x 0 , y 0 , 1 ] = a x 0 + b y 0 + c 當然，我們也可以寫成flag1 = abc[0]*dots[0][0]+abc[1]*dots[0][1]+c，隻是看上去不太優雅。
然後，得到瞭flag1和flag2的值，如果二者異號，那麼就可以斷定，直線在兩個點的中間。也就是說，隻要flag1*flag2<0,即可說明直線與線段有焦點。
當然，從實際的角度出發，我們可以不去考慮光線通過鏡片的端點或者平行鏡片並掠過這種情況，從而隻要函數返回值小於零，就認定二者相交，否則不相交。然而，從數學的角度出發，直線和線段之間可能存在三種關系：不相交、有一個交點、線段在直線上。
雖然這個理解沒什麼實際價值，但對於python學習來說，卻是非常有意義的一個例子，代碼如下，看懂瞭這個代碼，那麼就差不多可以算是一個初級的python程序員瞭。
def isCross(abc,dots):
abc = np.array(abc)
flags = [abc.dot(list(p)+[1]) for p in dots] #for表達式
poses,negs,zeros = [0,0,0]
for flag in flags: #循環語句
if flag > 0: #判斷語句
poses += 1
elif flag <0:
negs += 1
else:
zeros += 1
return poses*negs+zeros
這個短短的程序涵蓋瞭循環語句、判斷語句以及for表達式的內容，前兩者是最基礎的編程知識，後者是python中非常亮眼的一種功能。
首先來認識一下運算符+=，poses += 1即為poses = poses + 1，即相當於將poses+1賦值給poses。賦值前後flag1在內存中的位置發生瞭變化，也就是說flag1已經不是原來的flag1瞭。在這裡，等號也並不能讀成等於，而是讀成被賦值為。即poses被賦值為poses+1。前面略有提過，雙等號==才表示真正的相等。
然後來看判斷語句，對於表達式if a A elif b B else C，我們按照人類的語法去讀即可：如果a成立，則執行A，如果b成立，則執行B，否則的話執行C。在上述代碼中，也可以很方便地讀成：如果flag>0,那麼poses被賦值為poses+1；如果flag<0，那麼negs變成negs+1；否則的話zeros變成zeros+1。
這幾個變量顧名思義，poses表示正數的個數，negs表示負數的個數，zeros表示0的個數。
然後來看[abc.dot(list(p)+[1]) for p in dots]，我們首先使用一種沒有陌生字符的方式書寫：
#這是偽代碼，假設dots中有n個變量，表示創建flag1、flag2一直到flagn一共n個變量。 flag1 = abc.dot(list(dots[0])+[1]) flag2 = abc.dot(list(dots[1])+[1]) ... flagn = abc.dot(list(dots[n])+[1])
這個表達式非常規律，這n個變量相當於是在dots中循環一遍，然後逐個賦值。for p in dots表示的就是將dots中的元素取出，賦值給p，然後再對p進行操作abc.dot(list(p)+[1])，最後將所有操作得到的值包裹再一個list中。
最後再記一下這個表達形式 [abc.dot(list(p)+[1]) for p in dots]，以後會經常用到。
最後來看for flag in flags，即拿出flags中的所有元素，循環操作其下方的代碼塊。flags中的元素即為兩個點分別帶入 a x + b y + c 之後的值。那麼對於這兩個點來說，如果一正一負，則poses*negs=1*1=1,此時代表直線和線段有一交點，否則這個值便為0。當poses*negs==0時，則zeros的個數表示端點與直線相交的個數，zeros為0，表示無交點，為1，表示有一個端點在直線上，為2表示兩個端點都在直線上。
射線排序
現在，我們可以判斷某一個線段與一條直線是否有交點瞭，那麼如果空間中有多個平面鏡，光線所在的直線又與許多平面鏡有交點，那麼應該如何找到最近的那個呢？最簡單的方法是分別求取這些點到光源的距離，距離越近相交越早。但這樣會產生一個問題，即難以判定這個最近點是否在光的傳播路徑上，如果這個點在光源的後面，那就比較尷尬瞭。
所以，比較穩妥的方法是，按照射線的方向對所有點進行排序，那麼光源後面的那個點，就是光線傳播過程中的第一個交點。
剛剛我們在判定直線與線段的交點時，提到瞭直線族的概念。發現對於a、b取值相同的一組直線來說，其c值的大小與直線族的順序是密切相關的。如Fig2-2所示。其 c 1 到 c 4依次遞減。
這啟發我們需要構建出一組和光線想垂直的直線族[a,b,~]，則對於空間中任意一點 ( x , y )其所對應的 a x + b y 的值即能夠對射線上的點進行排序。
考慮到a、b的值可能為0，所以不適合求倒數，故采用[b,-a]作為特征直線族，用以評價點在射線上的位置，最終代碼如下。
def sortByRay(abc,points):
ab = np.array([abc[1],-abc[0]]) #特征直線族
pDict = {ab.dot(point):point for point in points}
keyList = list(pDict.keys()) #將pDict的兼職轉化成列表
keyList.sort(reverse=True) #對鍵列表進行排序
return [pDict[key] for key in keyList]
這裡又涉及到瞭一個新的數據類型，即字典。在理解字典之前，我們可以先回顧一下列表，我們可以把列表想象成一組值和自然序列的一一對應。對於列表test = [a,b,c,d]來說，有如下的對應關系{0：a,1:b,2:c,3:d}，所以我們可以通過test[0]來索引a，test[1]來索引b，以此類推。
那麼，現在我不想用自然數來索引瞭，我想通過一個標記來索引，所以希望能夠創建一個偽列表
dic = {3:5,4:15,12:22}，於是我們可以對此列表進行索引dic[3]==5,dic[4]==15,dic[12]==22。
這個偽列表就可以由字典來實現。這種索引關系就叫做鍵值對，我們通過一個鍵來索引一個值。
對於表達式pDict = {ab.dot(point):point for point in points}
表示通過point對points進行遍歷，即對於每個points中的point都進行ab.dot(point)這樣的點乘操作。於是得到瞭由特征直線族得到的特征值與點之間的一一對應關系。
pDict.keys()即可提取出字典中所有的鍵，pDict.values()可以提取出字典中所有的值。在此我們將所有的鍵提取出來之後，再將其轉化為列表。
然後即可調用列表的排序函數，將所有的值進行排序。即keyList.sort()，其中reverse參數默認為False，此時為降序。我們選擇True，此時表示升序。
線弧關系
目前，我們已經能夠精確地衡量射線與線段之間的關系瞭，接下來需要思考如何確定射線與透鏡的位置關系。這一點當然也要從交點說起。
首先，弧是圓的一部分，所以如果一條直線與弧有交點，那麼必然與這段弧所在的圓有交點。而直線與圓的交點判定相對來說還是非常容易的，隻要圓心到直線的距離小於半徑即可。
而直線和圓的交點問題則可以歸結為求解方程組：
# abc為光線參數；cir為圓參數
# point為光源位置，當其為[]時表示不考慮
def getCrossCircle(abc=[1,1,1],cir=[0,0,2],point=()):
c=np.array(cir[0:2]+[1]).dot(abc)
a2b2 = abc[0]**2+abc[1]**2
delt = a2b2*cir[2]**2-c**2
if delt<0: return [] #如果無交點，返回空列表[]
else: delt=np.sqrt(delt) #否則求解判別式
plusMinus = lambda a,b : np.array(set([a-b,a+b])) #定義函數plusMinus
yCross = plusMinus(-abc[1]*c,abc[0]*delt)/a2b2*[1,1]+cir[1]
xCross = plusMinus(-abc[0]*c,-abc[1]*delt)/a2b2*[1,1]+cir[0]
if point==[]:
return [(xCross[i],yCross[i]) for i in [0,1]]
yFlag = (yCross-point[1])*abc[0] >= 0
xFlag = (point[0]-xCross)*abc[1] >= 0
zFlag = np.abs(xFlag)+np.abs(yFlag) > 0
flag = yFlag*xFlag*zFlag
return [(xCross[i],yCross[i])
for i in range(len(yFlag)) if flag[i]]
這段程序雖然短，但信息量還是很大的,而且使用瞭一個lambda表達式。
plusMinus = lambda a,b : np.array([a-b,a+b])定義瞭一個名為plusMius的函數
這個函數寫成常規形式即為：
def plusMinus(a,b)
return np.array([a-b,a+b])
需要註意的是，lambda表達式的後面隻能有一個表達式，即隻能定義一行的函數。
在這段代碼中，我們還看到瞭一個陌生的運算符set，這也是python的一種數據類型，集合。和我們數學上認識的集合一樣，在集合中，不允許出現相同的值。所以，如果b==0的話，那麼set(a,a)=set(a)，即起到瞭去重的作用。然後再通過np.array將集合轉換成可計算的數組數據。
此外，這裡引入瞭比較運算符。我們目前所提到的運算都是數值型的，但實際上我們所接觸到的運算還有其他的類型。例如，當我們進行判斷的時候，if delt<0:中，<也是一種運算符，代表比較，如果delt的確小於0，那麼將返回一個值True，否則自然返回一個False。其中，True表示真，False表示假，這個是符合上過英語課的小學生的直覺的。
類似的運算符有>,<,>=,<=,==,!=，都可以望文生義地理解，其中!=表示不等於。這些都是比較運算符，其返回值為True和False。True和False是一種不同於數值的數據類型，叫佈爾型。
關系運算符不僅可以作用於數值類型，還可以作用於其他數據類型，一般情況下的使用方法都是符合直覺的。
>>> 1==2
False
>>> [3,3]==[3,3]
True
>>> 3>3
False
>>> 3>=3
True
然後我們再來看這個算法的邏輯，由於我們求解的是直線和圓的交點，而真實的光線卻是射線，那麼必然要考慮交點和光源的位置關系。
故代碼
yFlag = (yCross-point[1])*abc[0] >= 0 xFlag = (point[0]-xCross)*abc[1] >= 0
分別定義瞭這兩個判據xFlag、yFlag。但是當二者同時為0時，說明此時 x = x 0 , y = y 0 x=x_0,y=y_0 x=x0,y=y0，此時交點即光源，故不能算作光線與圓有交點。所以又有判據zFlag。
隻有當這三個判據都滿足時，我們所得到的值才是有效的，故總判據與這三個分判據是’與’的關系，所以寫為flag = yFlag*xFlag*zFlag。
此外，我們並不知道交點的個數，當判別式為0的時候，lambda表達式將隻有一個值傳出，這時的xCross和yCross中分別隻有一個元素；如果判別式大於0，則分別有兩個元素。這裡又涉及到array的一個優良特性，當維度不想等的兩個變量進行計算時，其會自動對低維數據進行合理的擴張，例如
>>> np.array([1,2,3])+4
array([5, 6, 7])
>>> np.array([[1],[2],[3]])+4
array([[5],
[6],
[7]])
最後，又出現瞭一個似曾相識的表達式
return [(xCross[i],yCross[i]) for i in range(len(yFlag)) if flag[i]]
這個表達式可以很容易地讀出來：遍歷flag，如果flag的值為真，則將對應的交點坐標放入列表，並返回有效的交點坐標。
這是對我們之前所使用的[... for ... in ...]的一種擴張，這種寫法簡潔而強大，非常推薦使用。
點弧關系
一般來說，在一個光學系統中很少出現一整個球，大部分情況下是由部分球面組成的各種透鏡。所以，作為二維的光路系統，可能更需要被處理的是光線和圓弧的關系，尤其是和劣弧的關系。
判定點在劣弧上的方法有很多種，例如弧ACB上任意一點關於AB的對稱點如果落入圓內，則為劣弧；如果落到園外，則為優弧；如果在圓上，說明AB是直徑，弧ACB為半圓。
在此，我們選取另一種方式。如圖所示，E為對於劣弧上任意一點，其與AB中點D的連線必然小於AB，否則即在優弧上。
所以，代碼為：
def isOnArc(point,arc):
arc = np.array(arc)
dAB = 0.5*np.linalg.norm(arc[0]-arc[1]) #AB/2長度
dCrossA = np.linalg.norm(0.5*(arc[0]+arc[1])-point) #ED長度
return dAB > dCrossA
因此，當滿足劣弧判定時，線圓交點即為線弧交點。
def getCrossArc(abc=[1,1,1],arc=[[0,1],[0,-1],[1,0]],point=[]):
if point == []:
return []
crossDict = {np.sqrt((p[0]-point[0])**2+(p[1]-point[1])**2):p
for p in getCrossCircle(abc,arc2cir(arc),point)
if (isOnArc(p,arc) and (p!=point))}
if crossDict == {}:
return []
return crossDict[min(crossDict.keys())]
機靈的同學其實很早就註意到，在定義函數的時候，其傳入參數竟然被賦瞭值。這在python中並不值得大驚小怪，此時輸入的值便是默認值。
此外，函數在被調用的時候，我們當然可以通過參數的順序進行傳參，但也可以使用變量名來對參數進行賦值，此時參數的順序便沒有意義瞭。
例如，
test1 = getCrossArc([1,1,1],[[0,1],[0,-1],[1,0]],(0,0)]
test2 = getCrossArc(abc=[1,1,1],point=(0,0),
arc=[[0,1],[0,-1],[1,0]])
上述兩種寫法都能得到正確的結果。
以上就是python光學仿真實現光線追跡之空間關系的詳細內容，更多關於python實現光線追跡的空間關系的資料請關註WalkonNet其它相關文章！
推薦閱讀：
- python光學仿真實現光線追跡折射與反射的實現
- 圖文詳解感知機算法原理及Python實現
- python Graham求凸包問題並畫圖操作
- java中lambda(函數式編程)一行解決foreach循環問題
- python opencv鼠標畫點之cv2.drawMarker()函數