空間關(guān)系

變化始于相遇，所以交點是一切的核心。

相交判定

首先考察一束光線能否打在某個平面鏡上。光線被抽象成了一個列表[a,b,c]，平面鏡則被抽象成為由兩個點構(gòu)成的線段[(x1,y1),(x2,y2)]。兩條直線的交點問題屬于初等數(shù)學(xué)范疇，需要先將線段轉(zhuǎn)換成直線的形式，然后再求交點。但是兩條直線的交點可能落在線段的外面，從而不具有判定的意義。

如果我們的光學(xué)系統(tǒng)中有大量的光學(xué)元件，那么如果有一種方法可以快速判斷光線是否與光學(xué)元件有交點，將會顯得更加快捷。那么，如果一個直線穿過某個線段，就意味著這條線段的兩個端點必然在直線的兩側(cè)。

import numpy as np
def isCross(abc,dots):
    abc = np.array(abc)                 #將abc的格式轉(zhuǎn)成數(shù)組
    flag1 = abc.dot(list(dots[0])+[1])    #數(shù)組的點積
    flag2 = abc.dot(list(dots[1])+[1])
    return flag1*flag2

我們非常熟悉運算符"+"，不過目前來說，只有數(shù)值和數(shù)組是支持加法運算的。所以，對于list(dots[0])+[1]這種表示著實讓人有些摸不到頭腦。

這個含義其實是符合人類直覺的。列表內(nèi)的元素個數(shù)是可變的，兩個列表相加可以理解為兩個列表銜接在一起。當(dāng)然，元組并不支持這種運算。
例如

>>> [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

通過加法運算符來銜接兩個列表，實際上相當(dāng)于新建了一個變量，需要開辟新的內(nèi)存空間。好在對于初學(xué)者來說這樣不容易出錯。

在numpy中，+、-、*、/這幾個運算符表示對應(yīng)位置元素的運算。如果想使用點乘等其他運算，需要調(diào)用numpy中的其他函數(shù)。

>>> 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 當(dāng)然，我們也可以寫成flag1 = abc[0]*dots[0][0]+abc[1]*dots[0][1]+c，只是看上去不太優(yōu)雅。

然后，得到了flag1和flag2的值，如果二者異號，那么就可以斷定，直線在兩個點的中間。也就是說，只要flag1*flag2<0,即可說明直線與線段有焦點。

當(dāng)然，從實際的角度出發(fā)，我們可以不去考慮光線通過鏡片的端點或者平行鏡片并掠過這種情況，從而只要函數(shù)返回值小于零，就認定二者相交，否則不相交。然而，從數(shù)學(xué)的角度出發(fā)，直線和線段之間可能存在三種關(guān)系：不相交、有一個交點、線段在直線上。

雖然這個理解沒什么實際價值，但對于python學(xué)習(xí)來說，卻是非常有意義的一個例子，代碼如下，看懂了這個代碼，那么就差不多可以算是一個初級的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:                              #循環(huán)語句
        if flag > 0:                                #判斷語句
            poses += 1
        elif flag <0:
            negs += 1
        else:
            zeros += 1
    return poses*negs+zeros

這個短短的程序涵蓋了循環(huán)語句、判斷語句以及for表達式的內(nèi)容，前兩者是最基礎(chǔ)的編程知識，后者是python中非常亮眼的一種功能。

首先來認識一下運算符+=，poses += 1即為poses = poses + 1，即相當(dāng)于將poses+1賦值給poses。賦值前后flag1在內(nèi)存中的位置發(fā)生了變化，也就是說flag1已經(jīng)不是原來的flag1了。在這里，等號也并不能讀成等于，而是讀成被賦值為。即poses被賦值為poses+1。前面略有提過，雙等號==才表示真正的相等。

然后來看判斷語句，對于表達式if a A elif b B else C，我們按照人類的語法去讀即可：如果a成立，則執(zhí)行A，如果b成立，則執(zhí)行B，否則的話執(zhí)行C。在上述代碼中，也可以很方便地讀成：如果flag>0,那么poses被賦值為poses+1；如果flag<0，那么negs變成negs+1；否則的話zeros變成zeros+1。

這幾個變量顧名思義，poses表示正數(shù)的個數(shù)，negs表示負數(shù)的個數(shù)，zeros表示0的個數(shù)。

然后來看[abc.dot(list(p)+[1]) for p in dots]，我們首先使用一種沒有陌生字符的方式書寫：

#這是偽代碼，假設(shè)dots中有n個變量，表示創(chuàng)建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])

這個表達式非常規(guī)律，這n個變量相當(dāng)于是在dots中循環(huán)一遍，然后逐個賦值。for p in dots表示的就是將dots中的元素取出，賦值給p，然后再對p進行操作abc.dot(list(p)+[1])，最后將所有操作得到的值包裹再一個list中。

最后再記一下這個表達形式 [abc.dot(list(p)+[1]) for p in dots]，以后會經(jīng)常用到。

最后來看for flag in flags，即拿出flags中的所有元素，循環(huán)操作其下方的代碼塊。flags中的元素即為兩個點分別帶入 a x + b y + c 之后的值。那么對于這兩個點來說，如果一正一負，則poses*negs=1*1=1,此時代表直線和線段有一交點，否則這個值便為0。當(dāng)poses*negs==0時，則zeros的個數(shù)表示端點與直線相交的個數(shù)，zeros為0，表示無交點，為1，表示有一個端點在直線上，為2表示兩個端點都在直線上。

射線排序

現(xiàn)在，我們可以判斷某一個線段與一條直線是否有交點了，那么如果空間中有多個平面鏡，光線所在的直線又與許多平面鏡有交點，那么應(yīng)該如何找到最近的那個呢？最簡單的方法是分別求取這些點到光源的距離，距離越近相交越早。但這樣會產(chǎn)生一個問題，即難以判定這個最近點是否在光的傳播路徑上，如果這個點在光源的后面，那就比較尷尬了。

所以，比較穩(wěn)妥的方法是，按照射線的方向?qū)λ悬c進行排序，那么光源后面的那個點，就是光線傳播過程中的第一個交點。

剛剛我們在判定直線與線段的交點時，提到了直線族的概念。發(fā)現(xiàn)對于a、b取值相同的一組直線來說，其c值的大小與直線族的順序是密切相關(guān)的。如Fig2-2所示。其 c 1 到 c 4依次遞減。

這啟發(fā)我們需要構(gòu)建出一組和光線想垂直的直線族[a,b,~]，則對于空間中任意一點 ( x , y )其所對應(yīng)的 a x + b y 的值即能夠?qū)ι渚€上的點進行排序。

考慮到a、b的值可能為0，所以不適合求倒數(shù)，故采用[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的兼職轉(zhuǎn)化成列表
    keyList.sort(reverse=True)              #對鍵列表進行排序
    return [pDict[key] for key in keyList]

這里又涉及到了一個新的數(shù)據(jù)類型，即字典。在理解字典之前，我們可以先回顧一下列表，我們可以把列表想象成一組值和自然序列的一一對應(yīng)。對于列表test = [a,b,c,d]來說，有如下的對應(yīng)關(guān)系{0：a,1:b,2:c,3:d}，所以我們可以通過test[0]來索引a，test[1]來索引b，以此類推。

那么，現(xiàn)在我不想用自然數(shù)來索引了，我想通過一個標(biāo)記來索引，所以希望能夠創(chuàng)建一個偽列表

dic = {3:5,4:15,12:22}，于是我們可以對此列表進行索引dic[3]==5,dic[4]==15,dic[12]==22。

這個偽列表就可以由字典來實現(xiàn)。這種索引關(guān)系就叫做鍵值對，我們通過一個鍵來索引一個值。

對于表達式pDict = {ab.dot(point):point for point in points}

表示通過point對points進行遍歷，即對于每個points中的point都進行ab.dot(point)這樣的點乘操作。于是得到了由特征直線族得到的特征值與點之間的一一對應(yīng)關(guān)系。

pDict.keys()即可提取出字典中所有的鍵，pDict.values()可以提取出字典中所有的值。在此我們將所有的鍵提取出來之后，再將其轉(zhuǎn)化為列表。

然后即可調(diào)用列表的排序函數(shù)，將所有的值進行排序。即keyList.sort()，其中reverse參數(shù)默認為False，此時為降序。我們選擇True，此時表示升序。

線弧關(guān)系

目前，我們已經(jīng)能夠精確地衡量射線與線段之間的關(guān)系了，接下來需要思考如何確定射線與透鏡的位置關(guān)系。這一點當(dāng)然也要從交點說起。

首先，弧是圓的一部分，所以如果一條直線與弧有交點，那么必然與這段弧所在的圓有交點。而直線與圓的交點判定相對來說還是非常容易的，只要圓心到直線的距離小于半徑即可。

而直線和圓的交點問題則可以歸結(jié)為求解方程組：

# abc為光線參數(shù)；cir為圓參數(shù)
# point為光源位置，當(dāng)其為[]時表示不考慮
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]))  #定義函數(shù)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的函數(shù)

這個函數(shù)寫成常規(guī)形式即為：

def plusMinus(a,b)
    return np.array([a-b,a+b])

需要注意的是，lambda表達式的后面只能有一個表達式，即只能定義一行的函數(shù)。

在這段代碼中，我們還看到了一個陌生的運算符set，這也是python的一種數(shù)據(jù)類型，集合。和我們數(shù)學(xué)上認識的集合一樣，在集合中，不允許出現(xiàn)相同的值。所以，如果b==0的話，那么set(a,a)=set(a)，即起到了去重的作用。然后再通過np.array將集合轉(zhuǎn)換成可計算的數(shù)組數(shù)據(jù)。

此外，這里引入了比較運算符。我們目前所提到的運算都是數(shù)值型的，但實際上我們所接觸到的運算還有其他的類型。例如，當(dāng)我們進行判斷的時候，if delt<0:中，<也是一種運算符，代表比較，如果delt的確小于0，那么將返回一個值True，否則自然返回一個False。其中，True表示真，False表示假，這個是符合上過英語課的小學(xué)生的直覺的。

類似的運算符有>,<,>=,<=,==,!=，都可以望文生義地理解，其中!=表示不等于。這些都是比較運算符，其返回值為True和False。True和False是一種不同于數(shù)值的數(shù)據(jù)類型，叫布爾型。

關(guān)系運算符不僅可以作用于數(shù)值類型，還可以作用于其他數(shù)據(jù)類型，一般情況下的使用方法都是符合直覺的。

>>> 1==2
False
>>> [3,3]==[3,3]
True
>>> 3>3
False
>>> 3>=3
True

然后我們再來看這個算法的邏輯，由于我們求解的是直線和圓的交點，而真實的光線卻是射線，那么必然要考慮交點和光源的位置關(guān)系。

故代碼

yFlag = (yCross-point[1])*abc[0] >= 0
xFlag = (point[0]-xCross)*abc[1] >= 0

分別定義了這兩個判據(jù)xFlag、yFlag。但是當(dāng)二者同時為0時，說明此時 x = x 0 , y = y 0 x=x_0,y=y_0 x=x0​,y=y0​，此時交點即光源，故不能算作光線與圓有交點。所以又有判據(jù)zFlag。

只有當(dāng)這三個判據(jù)都滿足時，我們所得到的值才是有效的，故總判據(jù)與這三個分判據(jù)是'與'的關(guān)系，所以寫為flag = yFlag*xFlag*zFlag。

此外，我們并不知道交點的個數(shù)，當(dāng)判別式為0的時候，lambda表達式將只有一個值傳出，這時的xCross和yCross中分別只有一個元素；如果判別式大于0，則分別有兩個元素。這里又涉及到array的一個優(yōu)良特性，當(dāng)維度不想等的兩個變量進行計算時，其會自動對低維數(shù)據(jù)進行合理的擴張，例如

>>> np.array([1,2,3])+4
array([5, 6, 7])
>>> np.array([[1],[2],[3]])+4
array([[5],
       [6],
       [7]])

最后，又出現(xiàn)了一個似曾相識的表達式

return [(xCross[i],yCross[i]) for i in range(len(yFlag)) if flag[i]]

這個表達式可以很容易地讀出來：遍歷flag，如果flag的值為真，則將對應(yīng)的交點坐標(biāo)放入列表，并返回有效的交點坐標(biāo)。

這是對我們之前所使用的[... for ... in ...]的一種擴張，這種寫法簡潔而強大，非常推薦使用。

點弧關(guān)系

一般來說，在一個光學(xué)系統(tǒng)中很少出現(xiàn)一整個球，大部分情況下是由部分球面組成的各種透鏡。所以，作為二維的光路系統(tǒng)，可能更需要被處理的是光線和圓弧的關(guān)系，尤其是和劣弧的關(guān)系。

判定點在劣弧上的方法有很多種，例如弧ACB上任意一點關(guān)于AB的對稱點如果落入圓內(nèi)，則為劣?。蝗绻涞綀@外，則為優(yōu)弧；如果在圓上，說明AB是直徑，弧ACB為半圓。

在此，我們選取另一種方式。如圖所示，E為對于劣弧上任意一點，其與AB中點D的連線必然小于AB，否則即在優(yōu)弧上。

所以，代碼為：

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

因此，當(dāng)滿足劣弧判定時，線圓交點即為線弧交點。

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())]

機靈的同學(xué)其實很早就注意到，在定義函數(shù)的時候，其傳入?yún)?shù)竟然被賦了值。這在python中并不值得大驚小怪，此時輸入的值便是默認值。

此外，函數(shù)在被調(diào)用的時候，我們當(dāng)然可以通過參數(shù)的順序進行傳參，但也可以使用變量名來對參數(shù)進行賦值，此時參數(shù)的順序便沒有意義了。

例如，

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

上述兩種寫法都能得到正確的結(jié)果。

以上就是python光學(xué)仿真實現(xiàn)光線追跡之空間關(guān)系的詳細內(nèi)容，更多關(guān)于python實現(xiàn)光線追跡的空間關(guān)系的資料請關(guān)注腳本之家其它相關(guān)文章！

最新評論