python光學仿真實現光線追跡的空間關系是什么

今天就跟大家聊聊有關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光學仿真實現光線追跡的空間關系是什么有進一步的了解嗎？如果還想了解更多知識或者相關內容，請關注億速云行業資訊頻道，感謝大家的支持。