空間關(guān)系

變化始于相遇，所以交點(diǎn)是一切的核心。

相交判定

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

如果我們的光學(xué)系統(tǒng)中有大量的光學(xué)元件，那么如果有一種方法可以快速判斷光線(xiàn)是否與光學(xué)元件有交點(diǎn)，將會(huì)顯得更加快捷。那么，如果一個(gè)直線(xiàn)穿過(guò)某個(gè)線(xiàn)段，就意味著這條線(xiàn)段的兩個(gè)端點(diǎn)必然在直線(xiàn)的兩側(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ù)組的點(diǎn)積
    flag2 = abc.dot(list(dots[1])+[1])
    return flag1*flag2

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

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

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

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

在numpy中，+、-、*、/這幾個(gè)運(yùn)算符表示對(duì)應(yīng)位置元素的運(yùn)算。如果想使用點(diǎn)乘等其他運(yùn)算，需要調(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表點(diǎn)積
32

所以， f l a g 1 = [ a , b , c ] ⋅ [ x 0 , y 0 , 1 ] = a x 0 + b y 0 + c 當(dāng)然，我們也可以寫(xiě)成flag1 = abc[0]*dots[0][0]+abc[1]*dots[0][1]+c，只是看上去不太優(yōu)雅。

然后，得到了flag1和flag2的值，如果二者異號(hào)，那么就可以斷定，直線(xiàn)在兩個(gè)點(diǎn)的中間。也就是說(shuō)，只要flag1*flag2<0,即可說(shuō)明直線(xiàn)與線(xiàn)段有焦點(diǎn)。

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

雖然這個(gè)理解沒(méi)什么實(shí)際價(jià)值，但對(duì)于python學(xué)習(xí)來(lái)說(shuō)，卻是非常有意義的一個(gè)例子，代碼如下，看懂了這個(gè)代碼，那么就差不多可以算是一個(gè)初級(jí)的python程序員了。

def isCross(abc,dots):
    abc = np.array(abc)
    flags = [abc.dot(list(p)+[1]) for p in dots]    #for表達(dá)式
    poses,negs,zeros = [0,0,0]
    for flag in flags:                              #循環(huán)語(yǔ)句
        if flag > 0:                                #判斷語(yǔ)句
            poses += 1
        elif flag <0:
            negs += 1
        else:
            zeros += 1
    return poses*negs+zeros

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

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

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

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

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

#這是偽代碼，假設(shè)dots中有n個(gè)變量，表示創(chuàng)建flag1、flag2一直到flagn一共n個(gè)變量。
flag1 = abc.dot(list(dots[0])+[1])
flag2 = abc.dot(list(dots[1])+[1])
...
flagn = abc.dot(list(dots[n])+[1])

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

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

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

射線(xiàn)排序

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

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

剛剛我們?cè)谂卸ㄖ本€(xiàn)與線(xiàn)段的交點(diǎn)時(shí)，提到了直線(xiàn)族的概念。發(fā)現(xiàn)對(duì)于a、b取值相同的一組直線(xiàn)來(lái)說(shuō)，其c值的大小與直線(xiàn)族的順序是密切相關(guān)的。如Fig2-2所示。其 c 1 到 c 4依次遞減。

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

考慮到a、b的值可能為0，所以不適合求倒數(shù)，故采用[b,-a]作為特征直線(xiàn)族，用以評(píng)價(jià)點(diǎn)在射線(xiàn)上的位置，最終代碼如下。

def sortByRay(abc,points):
    ab = np.array([abc[1],-abc[0]])         #特征直線(xiàn)族
    pDict = {ab.dot(point):point for point in points}
    keyList = list(pDict.keys())            #將pDict的兼職轉(zhuǎn)化成列表
    keyList.sort(reverse=True)              #對(duì)鍵列表進(jìn)行排序
    return [pDict[key] for key in keyList]

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

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

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

這個(gè)偽列表就可以由字典來(lái)實(shí)現(xiàn)。這種索引關(guān)系就叫做鍵值對(duì)，我們通過(guò)一個(gè)鍵來(lái)索引一個(gè)值。

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

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

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

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

線(xiàn)弧關(guān)系

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

首先，弧是圓的一部分，所以如果一條直線(xiàn)與弧有交點(diǎn)，那么必然與這段弧所在的圓有交點(diǎn)。而直線(xiàn)與圓的交點(diǎn)判定相對(duì)來(lái)說(shuō)還是非常容易的，只要圓心到直線(xiàn)的距離小于半徑即可。

而直線(xiàn)和圓的交點(diǎn)問(wèn)題則可以歸結(jié)為求解方程組：

# abc為光線(xiàn)參數(shù)；cir為圓參數(shù)
# point為光源位置，當(dāng)其為[]時(shí)表示不考慮
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 []          #如果無(wú)交點(diǎn)，返回空列表[]
    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]]

這段程序雖然短，但信息量還是很大的,而且使用了一個(gè)lambda表達(dá)式。

plusMinus = lambda a,b : np.array([a-b,a+b])定義了一個(gè)名為plusMius的函數(shù)

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

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

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

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

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

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

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

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

然后我們?cè)賮?lái)看這個(gè)算法的邏輯，由于我們求解的是直線(xiàn)和圓的交點(diǎn)，而真實(shí)的光線(xiàn)卻是射線(xiàn)，那么必然要考慮交點(diǎn)和光源的位置關(guān)系。

故代碼

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

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

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

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

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

最后，又出現(xiàn)了一個(gè)似曾相識(shí)的表達(dá)式

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

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

這是對(duì)我們之前所使用的[... for ... in ...]的一種擴(kuò)張，這種寫(xiě)法簡(jiǎn)潔而強(qiáng)大，非常推薦使用。

點(diǎn)弧關(guān)系

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

判定點(diǎn)在劣弧上的方法有很多種，例如弧ACB上任意一點(diǎn)關(guān)于AB的對(duì)稱(chēng)點(diǎn)如果落入圓內(nèi)，則為劣弧；如果落到園外，則為優(yōu)??；如果在圓上，說(shuō)明AB是直徑，弧ACB為半圓。

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

所以，代碼為：

def isOnArc(point,arc):
    arc = np.array(arc)
    dAB = 0.5*np.linalg.norm(arc[0]-arc[1])             #AB/2長(zhǎng)度
    dCrossA = np.linalg.norm(0.5*(arc[0]+arc[1])-point) #ED長(zhǎng)度
    return dAB > dCrossA

因此，當(dāng)滿(mǎn)足劣弧判定時(shí)，線(xiàn)圓交點(diǎn)即為線(xiàn)弧交點(diǎn)。

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

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

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

例如，

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

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

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