TWI594208B - 基於二維影像及三維完整八象限定位之內視鏡手術器械追蹤方法 - Google Patents

基於二維影像及三維完整八象限定位之內視鏡手術器械追蹤方法 Download PDF

Info

Publication number
TWI594208B
TWI594208B TW105135385A TW105135385A TWI594208B TW I594208 B TWI594208 B TW I594208B TW 105135385 A TW105135385 A TW 105135385A TW 105135385 A TW105135385 A TW 105135385A TW I594208 B TWI594208 B TW I594208B
Authority
TW
Taiwan
Prior art keywords
ring
image
dimensional
ring body
angle
Prior art date
Application number
TW105135385A
Other languages
English (en)
Other versions
TW201818347A (zh
Inventor
沈岱範
許家銓
Original Assignee
國立雲林科技大學
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 國立雲林科技大學 filed Critical 國立雲林科技大學
Priority to TW105135385A priority Critical patent/TWI594208B/zh
Application granted granted Critical
Publication of TWI594208B publication Critical patent/TWI594208B/zh
Publication of TW201818347A publication Critical patent/TW201818347A/zh

Links

Landscapes

  • Length Measuring Devices By Optical Means (AREA)
  • Image Analysis (AREA)

Description

基於二維影像及三維完整八象限定位之內視鏡手術器械追蹤方法
本發明係關於一種基於二維影像及三維完整八象限定位之內視鏡手術器械追蹤方法,可運用於實際手術當中,以及結合3D擴增實境技術(Augmented Reality,AR)在遠距手術環境之應用。
當前微創手術(MIS)在近幾十年來,已成為越來越流行的手術,也因微創手術大大的減低了病患的疼痛與縮短了恢復時間,使得微創手術廣為病人所接受。傳統上要完成的MIS,至少一個助手控制腹腔鏡攝像頭,和一個外科醫生,控制腹腔鏡器械進行手術是必要的。為了幫助外科醫生專注於手術,助理需要的責任,要確保手術器械的尖端要位於攝影機視圖的中央,隨者助理的幫助下,外科醫生可以直接觀察2D攝影機視圖去操作工具,但對於一個助理很難保持攝影相機一段時間的穩定型,此外由於2D照相機視圖缺乏深度資訊,一個外科醫生應該要經過足夠的訓練來進行手術,因此一個系統可以自動跟蹤器械並重建其方向和位置,在微創手術(MIS)的過程中變得更加重要。
按,目前電腦視覺系統取得二維影像中的深度訊息的傳統方法是利用同一物件在兩台以上的攝影機位置來估測。而本國專利申請案第104143279號「標記式圓桿狀物體二維影像之三維定位方法及其應用於內視鏡器械定位與追蹤之方法」之專利申請案中,曾提出的以均勻圓桿狀物體(如內視鏡之手術器 械)在單一鏡頭相機拍攝下,藉由在桿上兩個位置分別貼上兩個相同形狀,但不同顏色之環形標記,運用現代數位影像處理的相關技術,能有效的偵測出標記物以及圓桿狀物體在影像平面上所有的2D資訊,且透過鏡頭映射在相機感測器成像位置之幾何關係,快速求出圓桿狀物體3D定位資訊。此圓桿狀物體之3D定位參數共有七個。前三個參數為3D座標值,最後是In-plane角和Out-plane角三個角度,該專利申請案也提出在桿狀物體上環之標記,以進制編碼方式估測自軸旋轉角度。此七個參數決定唯一圓桿狀物體之3D姿勢。
然而,前述專利申請案係以均勻圓桿狀物體(如內視鏡之手術器械)在單一鏡頭環境中,可藉由在桿上兩個位置分別貼上兩個相同形狀,但不同顏色之環形標記來快速算出物體三維參數。但其只推導單一姿態的推導公式並不足以估測圓桿狀物任意姿態,並無在真實MIS微創手術的背景下執行,精準度上有改進空間且未做量測誤差的分析。
鑑於上述第104143279號專利申請案的技術問題點,於是本發明人窮盡心思研發出一種基於二維影像及三維完整八象限定位之內視鏡手術器械追蹤方法,故本發明最主要目的在於:提供快速又精確的內視鏡手術器械追蹤方法;本發明次要目的在於:提供在複雜器官背景中可辨識出雙環的完整輪廓。
為達上述目的,本發明運用如下的技術手段:一種運用三維完整八象限定位之內視鏡手術器械追蹤方法,係包含有:一雙環輪廓取出步驟,係將所輸入的投影於平面之手術器械影像擷取其第一環體(A)與第二環體(B)之平面影像;一影像座標轉換步驟,係設定前述第一環體及第二環體平面影像的中心為原點,並進行其座標轉換;一找出雙環重心的步驟,係分別計算出第一環 體及第二環體的重心座標;一取得雙環二維資訊步驟,將該第一環體及第二環的重心座標連線形成軸線段之直線方程式(L)並透過邊緣點以計算出該第一環體的二中心端點座標A1、A2及該第二環體的二中心端點座標B1、B2;一像素轉換毫米單位步驟,係將前述四個端點座標A1、A2、B1、B2轉換成毫米單位;一進行三維姿態估測步驟,係使用上個步驟終的端點座標A1、A2、B1、B2及已知的攝影機焦距(λ),環軸長(L),雙環重心軸長(LAB),雙環兩端投射點(xa1、xa2、xb1、xb2)和(ya1、ya2、yb1、yb2)之參數,並配合使用三維八像限定位系統,進而求出求該手術器械的三維定位參數{XA1,YA1,ZA1,α,β,γ}之參數,進而可快速算該手術器械在空間的任意姿態。
上述該雙環輪廓取出步驟中,係定義該第一環體為(Ra,Ga,Ba)及該第二環體為(Rb,Gb,Bb),運用通道差的概念,將該第一環體、第二環體影像中的每個像素做三通道相減取絕對值,在對相減後的三通道分別設其定門檻值,並定義如下數學式:【數學式1】第一環體:|R-Ra|<TRa,|G-Ga|<TGa,|B-Ba|<TBa第二環體:|R-Rb|<TRb,|G-Gb|<TGb,|B-Bb|<TBb;其中該門檻值設定主要對於每種顏色對於光源的容忍度,假設以上條件成立則屬於該環體的像素。
本發明藉由上述技術手段,可以達成如下功效:
1、本發明係將二維雙環輪廓推導出三維器械姿態六參數的關係式由原來只有三維第一象限,延伸至完整八象限,故相較於其他三維器械姿態定位參數估計演算法,本發明所提出的標記式的演算法相當快速精準且簡便。
2、本發明應用於MIS微創手術環境之雙環顏色選取,係對大量的複雜器官背景影像進行色彩模型分析,決定雙環的最優顏色,分析找出能夠在腹腔鏡環境中易於辨識的顏色,並對於因光照明而導致手術器械反光,挑選了理想RGB門檻值,能於實際手術中完整辨識出手術器械上的雙環輪廓。
a、b、c、d、e、f‧‧‧步驟
圖1:雙環圓桿狀物體與影像平面的關係示意圖。
圖2:參考點AA1為圓桿狀物軸線與環A上邊緣平面之交叉點示意圖。
圖3:α角度為X軸與投射於x-y(X-Y)影像平面線向量的夾角示意圖。
圖4:β角度為X軸與投射於X(x)-Z平面線向量的夾角(β為銳角)示意圖。
圖5:標記式圓桿狀物體在二維的影像X軸投射面(第一象限方向)示意圖。
圖6:標記式圓桿狀物體在二維的影像X軸投射面(第二象限方向)示意圖。
圖7:標記式圓桿狀物體在二維的影像X軸投射面(第三象限方向)示意圖。
圖8:標記式圓桿狀物體在二維的影像X軸投射面(第四象限方向)示意圖。
圖9:所有β角的可能投射於X(x)-Z平面(β角為銳角)示意圖。
圖10:γ角度為Y軸與投射於Y(y)-Z平面線向量的夾角(γ角為銳角)示意圖。
圖11:標記式圓桿狀物體在二維的影像Y軸投射面(第一象限方向)示意圖。
圖12:標記式圓桿狀物體在二維的影像Y軸投射面(第二象限方向示意圖。
圖13:標記式圓桿狀物體在二維的影像Y軸投射面(第三象限方向)示意圖。
圖14:標記式圓桿狀物體在二維的影像Y軸投射面(第四象限方向)示意圖。
圖15:所有γ角的可能投射於Y(y)-Z平面(γ角為銳角)示意圖。
圖16:座標點A1及AA1之間的幾何關係示意圖。
圖17:圓桿狀體一至四象限AA1點座標推導結果示意圖。
圖18:三維八象限姿態概念圖(Z軸正的方向為鏡頭方向)示意圖。
圖19:三維八象限姿態關係表示意圖。
圖20:本發明所提出獲取雙環輪廓演算法之示意圖。
圖21:MIS微創手術環境雙環顏色選擇之Step1至Step5流程圖。
圖22:本發明關於RP(x,y)與RIJ(x,y)做XOR之示意圖。
圖23:本發明雙環可適用顏色排名示意圖。
圖24:本發明挑選RGB門檻值流程圖。
圖25:本發明以單調的全黑為背景,擷取雙環之示意圖。
圖26:將多張的手術器械合成圖利用繪圖軟體PhotoImpactX3調整光源強度之示意圖。
圖27:本發明另一獲取雙環輪廓演算法流程圖之示意圖。
圖28:另一RP(x,y)與RJ(x,y)做邏輯運算XOR示意圖。
圖29:本發明標記式圓桿狀物體基於二維影像之三維八象限定位演算法之流圖。
圖30:本發明另一獲取雙環輪廓演算法流程圖示意圖。
圖31:本發明圖29演算法的座標轉換示意圖。
圖32:本發明圖29演算法的重心座標計算示意圖。
圖33:本發明圖29演算法的AB兩環重心點座標連線之示意圖。
圖34:本發明圖29演算法獲取AB兩環上下沿點圖。
關於標記式圓桿狀物體二維影像之三維定位系統八象限延伸推導:本發明以2015年本實驗室余明駿碩士論文(即第104143279號專利申請 案),提出的「標記式圓桿狀物體二維影像之三維定位系統及內視鏡器械追蹤應用」為基礎再延伸,該專利申請案可藉由在桿狀物體上分別在兩個位置上貼上相同形狀且不同顏色的環形標記,藉由兩個標記透過單台攝影機感測器成像位置之幾何關係,推導出圓桿狀物體的3D姿態,但目前該演算法只有單一姿態推導公式並不足以估測圓桿狀物任意姿態,且精準度尚有改進空間,故,本發明目的係推出圓桿狀物體3D任意姿態公式與精準度的提升。
關於桿狀物體物體3D六參數之推導原理: 如圖1所示者,假設原點的3D座標系統為(X,Y,Z)與相機底片之2D座標系統(x,y)重疊,X和x為垂直軸,Y和y為水平軸,Z軸為縱軸。在桿狀物體兩端之靠前端,塗上或標記相同寬度為L的第一環體(以下稱A環)和第二環體(以下稱B環),兩環預定距離則為LAB,為了精確和方便辨識2D影像中的環A及環B,建議使用能與使用情境容易區別的兩種不同的顏色,分別配給環A及環B。
將該圓桿狀物3D姿態定義成6個參數分別為(XAA1,YAA1,ZAA1and α,β,γ),其中(XAA1,YAA1,ZAA1)3D座標設為參考點AA1,如圖2所示,點AA1(BB1)是桿狀物軸線與環A(環B)上邊緣的平面焦點,因此該桿狀物的軸線可被表示為方向,令他為線向量。然而參考點AA1他是在桿狀物體的軸線上,它無法被看見,並不能從二維影像去中去做偵測。環A表面上的上邊緣點A1(XA1,YA1,ZA1)座標如圖2是用來估測先前推估的參考點AA1(XAA1,YAA1,ZAA1)3D座標,邊緣點A1(XA1,YA1,ZA1)座標可在二維影像中被偵測出。
令線向量是交叉於桿狀體的表面和軸線向量投射於影像平面的軌跡,如圖2所示,由於假設該桿狀物的直線是均勻的,在桿狀物 表面上的線是平行於軸線。要注意的是,點A1是在向量與環A上邊緣的交叉點,如圖2所示。
而圓桿狀物三個角度為α、β、γ,如圖1詳細定義如下:(1)α角度為X軸與投射於x-y(X-Y)影像平面線向量的夾角。α角度圍繞的Z軸旋轉,也被稱為In-Plan角。(2)β角度為X軸與投射於X(x)-Z平面線向量的夾角,β角度圍繞著Y軸旋轉,也被稱為out-plan角。(3)γ角度為y軸和投射於Z-Y(y)平面線向量的夾角,γ角圍繞著X旋轉,也被稱為Out-Plan角。上面描述的六參數(XAA1,YAA1,ZAA1and α,β,γ)足以決定該圓桿狀物的3D姿態。
關於In-planeα角之估測:
2D影像平面中的線向量投射於x-y平面上3D線向量。然而α角度是與X軸之間的夾角。線向量則是為影像中點A(x a ,y a )和B(x b ,y b )所連接起來的。A(x a ,y a )是環A透過相機鏡頭長度為λ的焦距投射於影像平面的區域的重心。同樣的,B(x b ,y b )是環A是透過同一台相機鏡頭長度為γ的焦距投射於影像平面的區域的重心。所以可得到如下公式:【數學式2】α=tan-1((|y a -y b |)/(|x a -x b |))
其中α為銳角。這邊列出了四種可能的α角,如圖3所示:1.假設xa<xb和ya>yb成立,α則為第一象限角度,即α I =α。2.假設xa<xb和ya>yb成立,α則為第二象限角度,即α II =180-α。3.假設xa<xb和ya<yb成立,α則為第三象限角度,即α III =180+α。4.假設xa>xb和ya<yb成立,α則為第四象限角度,即α IV =-α。整理成表格如表1所示(X-Y平面四象限α角):【表1】
關於Out-planeβ角及深度Z值估測: 在2D影像平面向量可以決定投射於3D方向線向量,藉由環A與環B的尋找投射於2D區域的質心,如圖4所示,x a1,x a2,x b1,x b2座標點是2D投影線向量與環A及環B上邊緣和下邊緣的交叉點如圖5所示,同樣的y a1,y a2,y b1,y b2座標點也是如此。
如圖4定義在X-Z平面,桿狀物向量與X軸夾角的角度稱為β角,然後再定義了該圓桿狀物在每個一象限姿態的β角,第一象限為β I =β,第二象限為β II =180-β,第三象限為β III =180+β,第四象限為β IV =-β,之後在分別詳細推導各象限的β角度及深度ZA1
在已知條件有雙環投射點x a1,x a2,x b1,x b2和相機焦距λ,環A及環B寬度L,環A及環B之間的間距L AB ,且使用相機的成像幾何關係及相似三小比例關係來計算出out-planeβ角及深度Z資訊。
關於第一象限out-planeβ及深度ZA1估測: 此小節介紹第一象限圓桿狀物體在上的標記在二維影像投射面差以及三角比例關係計算β角度及深度ZA1(本發明採用虛擬平面座標推導,較符合所取得影像的方位)。
從圖5(環A),將用環A與相機的成像幾何關係做推導公式: X A1-X A2=Lcosβ (I-1)
使用相機的成相幾何關係,得到
替代Eq.(I-2)(I-3)到Eq.(I-1),得到
進一步簡化Eq.(I-4)使用3D幾何關係 Z A2=Z A1-Lsinβ (I-5)
得到 x a1 Z A1-λx a1-x a2 Z A1+λx a2+Lx a2 sinβ=λLcosβ (I-6)
同樣的,從圖5(環B),得到下列公式:x b1 Z B1-λx b1-x b2 Z B1+λx b2+Lx b2 sinβ=λLcosβ (I-7)
為了簡單推導令:v 1=cosβ,v 2=sinβ, (I-8)注意v 1 2+v 2 2=1
簡化Eq.(I-6)(I-7)與Eq.(I-8),得到(x a1-x a2)Z A1-λ(x a1-x a2)+Lv 2 x a2=λLv 1 (I-9)
(x b1-x b2)Z B1-λ(x b1-x b2)+Lv 2 x b2=λLv 1 (I-10)
為了進一步簡化(I-9)(I-10),令A=(x a1-x a2),B=(x b1-x b2) (I-11)
並獲得Eq.(I-12)及(I-13)AZ A1-λA+Lv 2 x a2=Lv 1 λ (I-12)
BZ B1-λB+Lv 2 x b2=Lv 1 λ (I-13)
Eq.(I-12)和(I-13)分別除以A和B,得到:
Eq.(I-14)減去Eq.(I-15),得到
得到Z A1-Z B1+Lv 2 C=DLv 1 λ (I-18)
使用環A和環B在圖5中的相似三角形關係,得到:Z A1-Z B1=L AB sinβ=L AB v 2 (I-19)
Eq.(I-18)可以使用Eq.(I-19)做簡化L AB v 2+Lv 2 C=DLv 1 λ
(L AB +LC)v 2=DLv 1 λ (I-20)
Let
Eq.(I-21)被簡化如下v 2=Fv 1 (I-22)
由於v 1 2+v 2 2=1,使用Eq.(I-22),得到 (由於假設β角為銳角,所以v1取正值) (I-23)
(I-25-1)詳細如下:
β I =cos-1 v 1 (I-26)
第二象限out-planeβ及深度ZA1估測:此小節接介紹第二象限圓桿狀物體在上的標記在二維影像投射面差以及三角比例關係計算β角度及深度ZA1(本發明採用虛擬平面座標推導,較符合所取得影像的方位)。
從圖6(環A),將用環A與相機的成像幾何關係做推導公式:X A2-X A1=Lcosβ (II-1)
使用相機的成相幾何關係,得到
替代Eq.(II-2)(II-3)到Eq.(II-1),得到
進一步簡化Eq.(I-4)使用3D幾何關係Z A2=Z A1-Lsinβ (II-5)
得到x a2 Z A1-λx a2-x a1 Z A1+λx a1-Lx a2 sinβ=λLcosβ (II-6)
同樣的,從圖6(環B),得到下列公式:x b2 Z B1-λx b2-x b1 Z B1+λx b1-Lx b2 sinβ=λLcosβ (II-7)
為了簡單推導令:v 1=cosβ,v 2=sinβ, (II-8)
注意
簡化Eq.(II-6)(II-7)與Eq.(II-8),得到(x a2-x a1)Z A1-λ(x a2-x a1)-Lv 2 x a2=λLv 1 (II-9)
(x b2-x b1)Z B1-λ(x b2-x b1)-Lv 2 x b2=λLv 1 (II-10)
為了進一步簡化(II-9)(II-10),令A=(x a2-x a1),B=(x b2-x b1) (II-11)
並獲得Eq.(II-12)及(II-13)AZ A1-λA-Lv 2 x a2=λLv 1 (II-12)
BZ B1-λB-Lv 2 x b2=λLv 1 (II-13)
Eq.(II-12)和(II-13)分別除以A和B,得到:
Eq.(II-14)減去Eq.(II-15),得到
得到Z A1-Z B1-Lv 2 C=DLv 1 λ (II-18)
使用環A和環B在圖6中的相似三角形關係,得到: Z A1-Z B1=L AB sinβ=L AB v 2 (II-19)
Eq.(II-18)可以使用Eq.(II-19)做簡化L AB v 2-Lv 2 C=DLv 1 λ
(L AB -LC)v 2=DLv 1 λ (II-20)
Let
Eq.(II-21)被簡化如下v 2=Fv 1 (II-22)
由於v 1 2+v 2 2=1,使用Eq.(II-22),得到v 1 2+F 2 v 1 2=1和所以, (由於假設β角為銳角,所以v1取正值) (II-23)
From eq.(II-14)
(II-25-1)詳細如下:
β II =180-cos-1 v 1 (II-26)
第三象限out-planeβ及深度ZA1估測:
此小節介紹第三象限圓桿狀物體在上的標記在二維影像投射面差以及三角比例關係計算β角度及深度ZA1(本發明採用虛擬平面座標推導,較符合所取得影像的方位)。
從圖7(環A),將用環A與相機的成像幾何關係做推導公式:X A2-X A1=Lcosβ (III-1)
使用相機的成相幾何關係,得到
替代Eq.(III-2)(III-3)to Eq.(III-1),得到
進一步簡化Eq.(III-4)使用3D幾何關係 Z A2=Z A1+Lsinβ (III-5)
得到x a2 Z A1-λx a2-x a1 Z A1+λx a1+Lx a2 sinβ=λLcosβ (III-6)
同樣的,從圖7(環B),得到下列公式:x b2 Z B1-λx b2-x b1 Z B1+λx b1+Lx b2 sinβ=λLcosβ (III-7)
為了簡單推導令:v 1=cosβ,v 2=sinβ, (III-8)
注意v 1 2+v 2 2=1
簡化Eq.(III-6)(III-7)與Eq.(III-8),得到(x a2-x a1)Z A1-λ(x a2-x a1)+Lv 2 x a2=λLv 1 (III-9)
(x b2-x b1)Z B1-λ(x b2-x b1)+Lv 2 x b2=λLv 1 (III-10)
為了進一步簡化(III-9)(III-10),令A=(x a2-x a1),B=(x b2-x b1) (III-11)
並獲得AZ A1-λA+Lv 2 x a2=λLv 1 (III-12)
BZ B1-λB+Lv 2 x b2=λLv 1 (III-13)
Eq.(III-12)和(III-13)分別除以A和B,得到:
Eq.(III-14)減去Eq.(III-15),得到
得到(Z A1-Z B1)+Lv 2 C=DLv 1 λ (III-18)
使用環A和環B在圖7中的相似三角形關係,得到:Z A1-Z B1=-L AB sinβ=-L AB v 2 (III-19)
Eq.(III-18)可以使用Eq.(III-19)做簡化-L AB v 2+Lv 2 C=DLv 1 λ
(LC-L AB )v 2=DLv 1 λ (III-20)
Let
Eq.(III-21)被簡化如下v 2=Fv 1 (III-22)
由於v 1 2+v 2 2=1,使用Eq.(III-22),得到 v 1 2+F 2 v 1 2=1 and v 1 2=,thus, (由於假設β角為銳角,所以v1取正值) (III-23)
From eq.(III-14)
(III-25-1)詳細如下:
β III =180+cos-1 v 1 (III-26)
第四象限out-planeβ及深度ZA1估測:此小節介紹第四象限圓桿狀物體在上的標記在二維影像投射面差以及三角比例關係計算β角度及深度ZA1(本發明採用虛擬平面座標推導,較符合所取得影像的方位)。
從圖8(環A),將用環A與相機的成像幾何關係做推導公式:X A1-X A2=Lcosβ (IV-1)
使用相機的成相幾何關係,得到
替代Eq.(I-2)(I-3)到Eq.(I-1),得到
進一步簡化Eq.(IV-4)使用3D幾何關係 Z A2=Z A1+Lsinβ (IV-5)
得到x a1 Z A1-λx a1-x a2 Z A1+λx a2-Lx a2 sinβ=λLcosβ (IV-6)
同樣的,從圖8(環B),得到下列公式:x b1 Z B1-λx b1-x b2 Z B1+λx b2-Lx b2 sinβ=λLcosβ (IV-7)
為了簡單推導令:v 1=cosβ,v 2=sinβ, (IV-8)
注意v 1 2+v 2 2=1
簡化Eq.(IV-6)(IV-7)與Eq.(IV-8),得到(x a1-x a2)Z A1-λ(x a1-x a2)-Lv 2 x a2=λLv 1 (IV-9)
(x b1-x b2)Z B1-λ(x b1-x b2)-Lv 2 x b2=λLv 1 (IV-10)
為了進一步簡化(IV-9)(IV-10),令A=(x a1-x a2),B=(x b1-x b2) (IV-11)
並獲得Eq.(IV-12)及(IV-13)AZ A1-λA-Lv 2 x a2=Lv 1 λ (IV-12)
BZ B1-λB-Lv 2 x b2=Lv 1 λ (IV-13)
Eq.(IV-12)和(IV-13)分別除以A和B,得到:
Eq.(IV-14)減去Eq.(IV-15),得到
得到Z A1-Z B1-Lv 2 C=DLv 1 λ (IV-18)
使用環A和環B在圖8中的相似三角形關係,得到:Z A1-Z Bi=-L AB sinβ=-L AB v 2 (IV-19)
Eq.(IV-18)可以使用Eq.(IV-19),做簡化-L AB v 2-Lv 2 C=DLv 1 λ
-(L AB +LC)v 2=DLv 1 λ (IV-20)
Let
Eq.(IV-21)被簡化如下v 2=Fv 1 (IV-22)
由於v 1 2+v 2 2=1,使用Eq.(IV-22),得到 (由於假設β角為銳角,所以v1取正值)
(IV-25-1)詳細如下:
β IV =-cos-1 v 1 (IV-26)
下面列出了所有β角的可能如圖9所示,SD=size(ringA)-size(ringB),其中size()為環A或環B的總面積(像素點的數目)1.假設SD<0與xa>xb,成立,則圓桿狀物朝向第一象限,A I =(x a1-x a2),B I =(x b1-x b2)
β I =cos-1 v 1 (I-26)
2.假設SD<0與xa<xb,成立,則圓桿狀物朝向第二象限A II =(x a2-x a1),B II =(x b2-x b1)
β II =180-cos-1 v 1 (II-26)
3.假設SD>0與xa<xb,成立,則圓桿狀物朝向第三象限A III =(x a2-x a1),β III =(x b2-x b1)
β III =180+cos-1 v 1 (III-26)
4.假設SD>0與xa>xb,成立,則圓桿狀物朝向第四象限A IV =(x a1-x a2),B IV =(x b1-x b2)
β IV =-cos-1 v 1 (IV-26)
一至四象限β推導結果比較:根據上面在X-Z平面上圓桿狀物姿態的,每個象限的角度和深度ZA1所推導的結果都是不一樣的,如下表2:
Out-planeγ角及深度Z值估測
在2D影像平面向量可以決定投射於3D方向線向量,藉由環A與環B的尋找投射於2D區域的質心,如圖10所示,x a1,x a2,x b1,x b2座標點是2D投影線向量與環A及環B上邊緣和下邊緣的交叉點如圖11所示,同樣的y a1,y a2,y b1,y b2座標點也是如此。
如圖10所示,定義在Y-Z平面,桿狀物向量與X軸夾角的角度稱為γ角,定義了該圓桿狀物在每個象限姿態的γ角,第一象限為γ I =γ,第二象限為γ II =180-γ,第三象限為γ III =180+γ,第四象限為γ IV =-γ,在之後,會分別詳細推導各象限的γ角度及深度ZA1
已知條件有雙環投射點y a1,y a2,y b1,y b2和相機焦距λ,環A及環B寬度L,環A及環B之間的間距L AB ,使用相機的成像幾何關係及相似三小比例關係來計算出out-planeγ角及深度Z資訊。
第一象限out-planeγ及深度ZA1估測
此小節介紹第一象限圓桿狀物體在上的標記在二維影像投射面差以及三角比例關係計算γ角度及深度ZA1(本發明採用虛擬平面座標推導,較符合所取得影像的方位)。
從圖11(環A),將用環A與相機的成像幾何關係做推導公式:Y A1-Y A2=Lcosγ (I-1)
使用相機的成相幾何關係到
替代Eq.(I-2)(I-3)到Eq.(I-1),得到
進一步簡化Eq.(I-4)使用3D幾何關係Z A2=Z A1-Lsinγ (I-5)
得到y a1 Z A1-λy a1-y a2 Z A1+λy a2+Ly a2 sinγ=λLcosγ (I-6)
同樣的,從圖11(環B),得到下列公式:y b1 Z B1-λy b1-y b2 Z B1+λy b2+Ly b2 sinγ=λLcosγ (I-7)
為了簡單推導令:v 1=cosγ,v 2=sinγ, (I-8)
注意v 1 2+v 2 2=1
簡化Eq.(I-6)(I-7)與Eq.(I-8),得到(y a1-y a2)Z A1-λ(y a1-y a2)+Lv 2 y a2=λLv 1 (I-9)
(y b1-y b2)Z B1-λ(y b1-y b2)+Lv 2 y b2=λLv 1 (I-10)
為了進一步簡化(I-9)(I-10),令A=(y a1-y a2),B=(y b1-y b2) (I-11)
並獲得Eq.(I-12)及(I-13)AZ A1-λA+Lv 2 y a2=Lv 1 λ (I-12)
BZ B1-λB+Lv 2 y b2=Lv 1 λ (I-13)
Eq.(I-12)和(I-13)分別除以A和B,得到:
Eq.(I-14)減去Eq.(I-15),得到
得到Z A1-Z B1+Lv 2 C=DLv 1 λ (I-18)
使用環A和環B在圖11中的相似三角形關係,得到:Z A1-Z B1=L AB sinγ=L AB v 2 (I-19)
Eq.(I-18)可以使用Eq.(I-19)做簡化L AB v 2+Lv 2 C=DLv 1 λ
(L AB +LC)v 2=DLv 1 λ (I-20)
Let
Eq.(I-21)被簡化如下v 2=Gv 1 (I-22)
由於v 1 2+v 2 2=1,使用Eq.(I-22),得到v 1 2+G 2 v 1 2=1和v 1 2=,所以,
(由於假設γ角為銳角,所以v1取正值) (I-23)
(I-25-1)詳細如下:
γ I =cos-1 v 1
第二象限out-planeγ及深度ZA1估測: 此小節介紹第二象限圓桿狀物體在上的標記在二維影像投射面差以及三角比例關係計算β角度及深度ZA1(本發明採用虛擬平面座標推導,較符合所取得影像的方位)。
從圖12(環A),將用環A與相機的成像幾何關係做推導公式: Y A2-Y A1=Lcosγ (II-1)
使用相機的成相幾何關係,得到
替代Eq.(II-2)(II-3)到Eq.(II-1),得到
進一步簡化Eq.(II-4)使用3D幾何關係 Z A2=Z A1-Lsinγ (II-5)
得到y a2 Z A1-λy a2-y a1 Z A1+λy a1-Ly a2 sinγ=λLcosγ (II-6)
同樣的,從圖12(環B),得到下列公式:y b2 Z B1-λy b2-y b1 Z B1+λy b1-Ly b2 sinγ=λLcosγ (II-7)
為了簡單推導令v 1=cosγ,v 2=sinγ, (II-8)
注意v 1 2+v 2 2=1
簡化Eq.(II-6)(II-7)與Eq.(II-8),得到(y a2-y a1)Z A1-λ(y a2-y a1)-Lv 2 y a2=λLv 1 (II-9)
(y b2-y b1)Z B1-λ(y b2-y b1)-Lv 2 v b2=λLv 1 (II-10)
為了進一步簡化(II-9)(II-10),令A=(y a2-y a1),B=(y b2-y b1) (II-11)
並獲得Eq.(II-12)及(II-13)AZ A1-λA-Lv 2 y a2=λLv 1 (II-12)
BZ B1-λB-Lv 2 y b2=λLv 1 (II-13)
Eq.(II-12)和(II-13)分別除以A和B,得到:
Eq.(II-14)減去Eq.(II-15),得到
得到Z A1-Z B1-Lv 2 C=DLv 1 λ (II-18)
使用環A和環B在圖12中的相似三角形關係,得到:Z A1-Z B1=L AB sinγ=L AB v 2 (II-19)
Eq.(II-18)可以使用Eq.(II-19)做簡化L AB v 2-Lv 2 C=DLv 1 λ
(L AB -LC)v 2=DLv 1 λ (II-20)
Let
Eq.(II-21)被簡化如下v 2=Gv 1 (II-22)
由於v 1 2+v 2 2=1,使用Eq.(II-22),得到v 1 2+G 2 v 1 2=1和v 1 2=,所以, (由於假設γ角為銳角,所以v1取正值) (II-23)
From eq.(II-14)
(II-25-1)詳細如下:
β II =180-cos-1 v 1 (II-26)
第三象限out-planeγ及深度ZA1估測:此小節介紹第三象限圓桿狀物體在上的標記在二維影像投射面差以及三角比例關係計算角度及深度ZA1(本發明採用虛擬平面座標推導,較符合所取得影像的方位)。
從圖13(環A),將用環A與相機的成像幾何關係做推導公式:Y A2-Y A1=Lcosγ (III-1)
使用相機的成相幾何關係,得到
替代Eq.(III-2)(III-3)到Eq.(III-1),得到
進一步簡化Eq.(III-4)使用3D幾何關係Z A2=Z A1+Lsinγ (III-5)
得到y a2 Z A1-λy a2-y a1 Z A1+λy a1+Ly a2 sinγ=λLcosγ (III-6)
同樣的,從圖13(環B),得到下列公式:y b2 Z B1-λy b2-y b1 Z B1+λy b1+Ly b2 sinγ=λLcosγ (III-7)
為了簡單推導令:v 1=cosβ,v 2=sinβ, (III-8)
注意v 1 2+v 2 2=1
簡化Eq.(III-6)(III-7)與Eq.(III-8),得到(y a2-y a1)Z A1-λ(y a2-y a1)+Lv 2 y a2=λLv 1 (III-9)
(y b2-y b1)Z B1-λ(y b2-y b1)+Lv 2 y b2=λLv 1 (III-10)
為了進一步簡化(III-9)(III-10),令A=(y a2-y a1),B=(y b2-y b1) (III-11)
並獲得AZ A1-λA+Lv 2 y a2=λLv 1 (III-12)
BZ B1-λB+Lv 2 y b2=λLv 1 (III-13)
Eq.(III-12)和(III-13)分別除以A和B,得到:
Eq.(III-14)減去Eq.(III-15),得到
得到(Z A1-Z B1)+Lv 2 C=DLv 1 λ (III-18)
使用環A和環B在圖13中的相似三角形關係,得到:Z A1-Z B1=-L AB sinγ=-L AB v 2 (III-19)
Eq.(III-18)可以使用Eq.(III-19)做簡化-L AB v 2+Lv 2 C=DLv 1 λ
(LC-L AB )v 2=DLv 1 λ (III-20)
Let
Eq.(III-21)被簡化如下v 2=Gv 1 (III-22)
由於v 1 2+v 2 2=1,使用Eq.(III-22),得到 (由於假設γ角為銳角,所以v1取正值) (III-23)
From eq.(III-14)
(III-25-1)詳細如下:
γ III =180+cos-1 v 1
第四象限out-planeγ及深度ZA1估測:此小節介紹第四象限圓桿狀物體在上的標記在二維影像投射面差以及三角比例關係計算β角度及深度ZA1(本發明採用虛擬平面座標推導,較符合所取得影像的方位)。
從圖14(環A),將用環A與相機的成像幾何關係做推導公式:Y A1=Y A2=Lcosβ (IV-1)
使用相機的成相幾何關係,得到
替代Eq.(IV-2)(IV-3)到Eq.(IV-1),得到
進一步簡化Eq.(IV-4)使用3D幾何關係 Z A2=Z A1+Lsinγ (IV-5)
得到y a1 Z A1-λy a1-y a2 Z A1+λy a2-Ly a2 sinγ=λLcosγ (IV-6)
同樣的,從圖14(環B),得到下列公式:y b1 Z B1-λy b1-y b2 Z B1+λy b2-Ly b2 sinγ=λLcosγ (IV-7)
為了簡單推導令: v 1=cosγ,v 2=sinγ, (IV-8)
注意v 1 2+v 2 2=1
簡化Eq.(IV-6)(IV-7)與Eq.(IV-8),得到(y a1-y a2)Z A1-λ(y a1-y a2)-Lv 2 y a2=λLv 1 (IV-9)
(y b1-y b2)Z B1-λ(y b1-y b2)-Lv 2 y b2=λLv 1 (IV-10)
為了進一步簡化(IV-9)(IV-10),令A=(y a1-y a2),B=(y b1-y b2) (IV-11)
並獲得Eq.(IV-12)及(IV-13)AZ A1-λA-Lv 2 y a2=Lv 1 λ (IV-12)
BZ B1-λB-Lv 2 y b2=Lv 1 λ (IV-13)
Eq.(IV-12)和(IV-13)分別除以A和B,得到:
Eq.(IV-14)減去Eq.(IV-15),得到
得到Z A1-Z B1-Lv 2 C=DLv 1 λ (IV-18)
使用環A和環B在圖14中的相似三角形關係,得到:Z A1-Z B1=-L AB sinγ=-L AB v 2 (IV-19)
Eq.(IV-18)可以使用Eq.(IV-19),做簡化-L AB v 2-Lv 2 C=DLv 1 λ
-(L AB +LC)v 2=DLv 1 λ (IV-20)
Let
Eq.(IV-21)被簡化如下v 2=Gv 1 (IV-22)
由於v 1 2+v 2 2=1,使用Eq.(IV-22),得到v 1 2+G 2 v 1 2=1和v 1 2=,所以, (由於假設β角為銳角,所以v1取正值)
(IV-25-1)詳細如下:
γ IV =-cos-1 v 1
下面列出了所有γ角的可能,如圖15所示SD=size(ringA)-size(ringB),其中Size()為環A或環B的總面積:1.假設SD<0與ya>yb,成立,則圓桿狀物朝向第一象限, A I =(y a1-y a2),B I =(y b1-y b2)
γ I =cos-1 v 1 (I-26)
2.假設SD<0與ya<xb,成立,則圓桿狀物朝向第二象限A II =(y a2-y a1),B II =(y b2-y b1)
sinγ=v 2=G II v 1
γ II =180-cos-1 v 1 (II-26)
3.假設SD>0與ya<yb,成立,則圓桿狀物朝向第三象限A III =(y a2-y a1),B III =(y b2-y b1)
γ III =180+cos-1 v 1 (III-26)
4.假設SD>0與ya>yb,成立,則圓桿狀物朝向第四象限A IV =(y a1-y a2),B IV =(y b1-y b2)
γ IV =-cos-1 v 1 (IV-26)
一至四象限γ推導結果比較:根據上面在Y-Z平面上圓桿狀物姿態的γ,每個象限的γ角度和深度ZA1所推導的的結果都是獨立的,如下表3:
估測3D座標參考點AA1:在圖16中,以估測表面的3D座標A1(XA1,YA1,ZA1)為基礎,來找出軸心參考點座標(XAA1,YAA1,ZAA1)。給一個影像座標(x a1,y a1),ZA1λ,如果想要估測桿狀體表面的3D座標A1(XA1,YA1,ZA1)及參考點AA1,步驟如下:Step1:在桿狀體表面估測A1(XA1,YA1,ZA1)3D座標時,由於A1點沒有深度ZA1資訊,將藉由求得β與γ角的過程得到ZA1,可利用影像幾何中的逆透視轉換得到XA1和YA1,公式如下:
Step2:估測AA 1(X AA1,Y AA1,Z AA1)3D座標。
從圖16中,
XAA1=XA1-△x YAA1=YA1-△y ZAA1=ZA1+△z (29) △x=Rsinβ,△z=Rcosβ,△y=Rsinγ,R是桿狀體的半徑。所以可以獲得定義雙環圓桿狀物3D姿態6參數:XAA1,YAA1,ZAA1and α,β,γ
延伸至其他象限的AA1推導結果: 上述推導出了圓桿狀物在第一象限AA1點的推導,本節將對每個象限做推導結果的比較。將該圓桿狀物擺出以鏡心為中心點的四個象限姿態,發現當桿狀體在鏡心左邊時表面看的到的,當桿狀體在鏡心右邊時,表面則是完全擋住,所以依據圖17,將所有不同的AA1點座標的公式做推導整理如表4:
無論對X軸或Y軸,相對於靜心的位置為正的(右邊或上面)則需扣掉△z,反之若為負(左邊或下邊)則須加上△z。
總合圓桿3D八象限三角度姿態與位置推導關係: 可以想像圓桿狀物在3維空間裡面會有八種象限的姿態,如圖18,因三維空間有三個平面,也將這3個平面定義成3三個角度α、β、γ,α角則是以X-Y平面定義、β角則是以Y-Z平面定義、γ角則是以X-Z平面定義,再以圖18去分析出八種象限姿態的關係表,如圖19,即可更明確的知道八種姿態所需條件關係。其中Q1~Q4與鏡頭同向,Q5~Q8指向鏡頭。
總結,原演算法只有單一姿態象限推導公式,並不足以計算出圓桿狀物體的任意姿態,所以將3D空間定義為八象限,推出八種象限圓桿狀物姿態,α角則是以X-Y平面定義、γ角則是以Y-Z平面定義、β角則是以X-Z平面定義,推導結果α、β、γ、ZA1分別列於表1~表3。另,該桿狀物軸心點AA1則會依據以鏡頭為中心的四個象限姿態,無論對X軸或Y軸,相對於靜心的位置為正的(右邊或上面)則需扣掉△z,反之若為負(左邊或下邊)則須加上△z。
關於標記式圓桿狀物體應用於MIS微創手術環境之雙環顏色選取: 本發明將延伸至MIS微創手術的實際環境,且在複雜器官背景找出雙環完整輪廓,進而取得3D定位資訊所需要的2D資訊,且對複雜器官背景進行色彩模型分析,決定雙環的顏色,分析找出能夠在腹腔鏡環境中易於辨識的顏色,以利找出環形標記完整輪廓,在複雜器官背景快速且完整辨識出雙環之輪廓,並取得精準的二維資訊。
本發明的另一目的是找出適合在真實手術環境中容易辨識的兩種標記環顏色,係運用如下步驟: Step1:收集大量具代表性的內視鏡複雜器官影像,並將做拼貼動作及Histogram分析找出使用最少的2種顏色;Step2:再以單調背景的雙環手術刀,找出雙環任一環當作精準答案影像(這邊以A環當精準答案);Step3:利用繪圖軟體PhotoImpactX3,將物件化的手術器械刀與真實環境做結合,並對環上的可能顏色(除第一步驟分析出來的顏色,再加入8個容易分辨的極端顏色作為可能的顏色選項稱為可能顏色),令他為實驗影像;Step4:將實驗影像套至本發明所提出之獲取雙環輪廓演算法,獲取A環輪廓,令他為濾出A環實驗影像; Step5:完美答案R_P(x,y)與濾出A環實驗影像R_j^A(x,y)做XOR找出差異像素點並計算變異數分析(j=1~50);Step6:依據step5分析結果做出雙環顏色選擇。圖21為前述Step1至Step5之流程圖。
獲取雙環輪廓演算法:套至本發明所提出之獲取雙環輪廓演算法,令影像上的顏色定義為(R,G,B),對環設定的顏色定義A環為(Ra,Ga,Ba)及B環為(Rb,Gb,Bb),利用通道差的概念,將影像中的每個Pixel做三通道相減取絕對值,在對相減後的三通道分別設定門檻值,分別定義為:【數學式3】A環:|R-Ra|<TRa,|G-Ga|<TGa,|B-Ba|<TBa
其中門檻值設定主要對於每種顏色對於光源的容忍度。假設以上條件成立則屬於該環Pixel。流程圖如圖20所示者(例如:PixelA環)。
各個顏色性能分析:
1.令精確Marker位置答案為Rp(x,y),令該Step5的實際複雜器官背景Marker位置為RIJ(x,y)做邏輯運算XOR(exclusive or)互斥或閘:D IJ =R P R IJ
2.計算每N個Case,精確Marker位置Rp(x,y)二元影像與實際複雜器官背景Marker位置二元影像RIJ(x,y)做邏輯運算XOR(exclusive or),找出差異像素點,如圖22所示者:
3.在這邊同時進行變異數(Variance)的運算,變異數公式如下:
其中公式中x i 相當於我的實際複雜器官背景Marker位置二元影像RIJ(x,y),μ相當於精確Marker位置Rp(x,y)二元影像,n則是實驗的張數。
4.記錄每一種顏色XOR(exclusive or)運算後的總數與變異數制實驗數據表。
以下是關於每種可能實驗顏色使用50張複雜器官背景所計算出的實驗數據,接下來皆適當的調整門檻值TR、TG、TB觀測每種顏色N IJ 與變異數(Variance)挑出最適當的兩種顏色當作本發明雙環顏色。
從下表5的數據發現(255,255,255)與(0,0,0)這兩總顏色已經受的影響, N IJ 數在門檻值均設為1時變大,變異數(Variance)也隨之變大(255,255,255)全白顏色受的影響是最大的,基本上可以淘汰(255,255,255)與(0,0,0)。
下表6(門檻值均設為10)、表7(門檻值均設為30)的數據觀察除了(255,255,255)與(0,0,0)的N IJ 與變異數(Variance)隨著門檻值越高遞增以外,其餘顏色皆不受影響。
由上述表6到表7可看出門檻值10調整到30的情況下,還看不出來所要選取的兩種顏色,所以決定將門檻值直接調到120觀測數據。下述表8(門檻值均設為10)的數據觀察除了極端值顏色(0,255,255)與(0,255,0)受影響較小,其餘顏色的N IJ 與變異數(Variance)隨門檻值越高遞增,這時可以選出標 記環在真實手術環境最適合的兩種顏色,圖23還做了真實環境標記環理想顏色的排名。
MIS微創手術環境決定標記環RGB門檻值:
因真實手術環境手術器械會因照明而反光,將以繪圖軟體(PhotoimpacX3t)模擬之,依據本發明獲取雙環輪廓演算法,調整並決定適當的RGB門檻值(光源容忍度)以獲取清晰完整雙環輪廓,其方法流程圖,如圖24所示者,其目的是找出適合在真實手術環境中合適的RGB門檻值。
Step1:以單調的全黑為背景,找出A環與B環位置,當作實驗完美(perfect)答案Rp(x,y);Step2:將多張的手術器械合成圖利用繪圖軟PhotoimpactX3打入反光,令影像為測試影像;Step3:將實驗影像套至本發明所提出之獲取雙環輪廓演算法,獲取雙環輪廓,令他為濾出A環實驗影像;Step5:完美答案R P (x,y)與濾出雙環實驗影像R P (x,y)做XOR找出差異像素點並計算變異數分析(j=1~20);及Step6:依據前述step5分析結果做出門檻值選擇。
模擬真實環境方法: 本發明以單調的全黑為背景如圖25,找出A環與B環位置。
將多張的手術器械合成圖利用繪圖軟體PhotoImpactX3調整光源強度,在PhotoImpactX3中的環境光線參數,可以調整光源的強度。將利用此參隊原影像(零級光源),數調整光源的強、中、弱。並定義強、中、弱,如下:環境光線參數為30時,令他為一級光源(弱),如圖26-1;環境光線參數為60時,令他為二級光源(中),如圖26-2;環境光線參數為90時,令他為三級光源(強),如圖26-3。
獲取雙環輪廓擷取演算法:令影像上的顏色定義為(R,G,B),對環設定的顏色定義A環為(Ra,Ga,Ba)及B環為(Rb,Gb,Bb),利用通道差的概念,將影像中的每個Pixel做三通道相減取絕對值,在對相減後的三通道分別設定門檻值,分別定義為:【數學式6】A環:|R-Ra|<TRa,|G-Ga|<TGa,|B-Ba|<TBa B環:|R-Rb|<TRb,|G-Gb|<TGb,|B-Bb|<TBb其中,門檻值設定主要對於每種顏色對於光源的容忍度。假設以上條件成立則屬於該環Pixel。流程圖如圖27(例如:PixelA環)。
各RGB門檻值性能分析:令精確Marker位置答案為RP(x,y),令實際複雜器官背景Marker位置為RJ(x,y)做邏輯運算XOR(exclusive or)互斥或閘:【數學式7】D j =|R P R j |(集合的點數)其示意圖,如圖28。其目的為計算兩者相異的點數D J 。(XOR規則:兩者相同為0,不同則為1)。
將上述RP(x,y)與Rj(X,Y)經過互斥或運算的得到的 變更門檻值並重複上述步驟,並比較各D Th 。其目的:挑選各級光源出最佳RGB門檻值。
零級光源實驗結果如下: 為了要找出零級光源的理想門檻值,首先從大範圍開始尋找,如表9,由1開始之後從等差為10的加至250。可以觀察到門檻值為1時DTH為最小的,這時可縮小範圍到門檻值30左右去尋找理想門檻值,結果如表10。
故,從表9擴大範圍尋找,可以觀察到理想門檻值可能再為30附近,這時再縮小範圍從門檻值30附近去尋找,如表10,發現從門檻值25開始的D Th 是遞增的狀態,這時可以知道零級光源理想門檻值為25。
一級光源實驗結果如下: 為了要找出一級光線的理想門檻值,首先從大範圍開始尋找,如表11,由1開始之後從等差為10的加至250。可以觀察到門檻值為130時DTH為最小的,這時可縮小範圍到門檻值130左右去尋找理想門檻值,結果如表12。
從表11擴大範圍尋找,可以觀察到理想門檻值可能再為130附近,這時再縮小範圍從門檻值130附近去尋找,如表12,發現在130以上或以下的D Th 是遞增的狀態,這時可以知道一級光源理想門檻值為130。
二級光源實驗結果如下: 為了要找出二級光線的理想門檻值,首先從大範圍開始尋找,如表13,由1開始之後從等差為10的加至250。可以觀察到門檻值為60時DTH為最小的,這時可縮小範圍到門檻值60左右去尋找理想門檻值,結果如表14。
從表13擴大範圍尋找,可以觀察到理想門檻值可能再為60附近,這時再縮小範圍從門檻值60附近去尋找,如表14,發現從門檻值62以上開始的D Th 是遞增的狀態,這時可以知道二級光源理想門檻值為62。
三級光源實驗結果如下: 為了要找出三級光線的理想門檻值,首先從大範圍開始尋找,如表15,由1開始之後從等差為10的加至250。可以觀察到門檻值為1時DTH為最小的,這時可縮小範圍到門檻值1左右去尋找理想門檻值,結果如表16。
【表15】
從表15擴大範圍尋找,可以觀察到理想門檻值可能再為1附近,這時再縮小範圍從門檻值1附近去尋找,如表16,發現從門檻值1以上開始的D Th 是遞增的狀態,這時可以知道結果三級光源理想門檻值為1。
經過上述實驗,將光源調整強、中、弱、無,決定每個光源級別(零級、一級、二級、三級)的理想RGB門檻值。1.零級光源(無光源效果)門檻值為25。2.一級光源(弱)門檻值為130。3.二級光源(中)門檻值為62。4.三級光源(強)門檻值為1。
因此,將對大量的複雜器官背景影像進行色彩模型分析,分析找出能夠在腹腔鏡環境中易於辨識的兩種作為雙環的顏色,A環為(0,255,255),B環為(0,255,0)。而對於因光照明而導致手術器械反光,利用繪圖軟體PhotoimpactX3模擬光源,並定義光源強度(無光源、強、中、弱),決定了各級光源的理想RGB門檻值,零級光源(無光源效果)門檻值為25。一級強光(弱)門檻值 為130。二級強光(中)門檻值為62。三級強光(強)門檻值為1。接者擺出器械姿態並觀察光源對六參數的影響。觀察結果光源強度對3D六參數的結果,還是微小的誤差,可能原因:因在獲取雙環輪廓時各級光線與完美答案的所濾出之輪廓的相異位置分布在雙環的邊緣,造成在獲取雙環上沿點與下沿點時會有差異性,導致微小的誤差。
標記式圓桿狀物體三維定位演算法應用於內視鏡手術:
本發明基於前述標記式圓桿狀物體二維影像之三維八象限定位演算法流程,並於圓桿狀物在不同的姿態及距離對三維定位參數之影響。基於前述推導公式,提出標記式圓桿狀物體基於單張二維影像之三維八象限參數估計演算法並應用於內視鏡手術環境,做詳細步驟介紹。
首先,於前述MIS微創手術環境中雙環顏色選擇的實驗結果如,圖23(只列出前十名)決定雙環的顏色,本發明選擇了第一名與第二名的顏色,分別為:(1)A環:(0,255,255);(2)B環:(0,255,0),這兩種顏色在真實手術環境(如圖5.1)中能有效地分辨出圓桿狀物體的雙環標記。
標記式圓桿狀物體基於二維影像之三維八象限定位演算法流程:
本發明的標記式圓桿狀物體單張二維影像作輸入,再透過現代影像處理的技術,取得圓桿狀物的二維資訊,再透過第三章所推導之公式即可快速算出標記式圓桿狀物體在三維空間中的八象限任意姿態參數。演算法流程圖,如圖29。
Step1.取出雙環輪廓: 令影像上的顏色定義為(R,G,B),對環設定的顏色定義A環為(Ra,Ga,Ba)及B環為(Rb,Gb,Bb),利用通道差的概念,將影像中的每個Pixel 做三通道相減取絕對值,在對相減後的三通道分別設定門檻值,分別定義為,如前述的數學式6;其中門檻值設定主要對於每種顏色對於光源的容忍度。假設以上條件成立則屬於該環Pixel,如圖30所示者。例如:PixelA環。
Step2.影像座標轉換: 因影像開發軟體的(0,0)是由圖片左上方開始,但推導公式的過程中則是以影像的中心為(0,0),所以必須將座標轉換將原點移至中心,如圖31,令輸入的是一張MXN大小影像或視訊的影格為f(i,j),要將f(i,j)轉成f(x,y),公式如下: 【數學式9】f(x,y)=(i-M/2,N/2-j)
Step3.找出雙環重心: 由前述的Step1敘述將圓桿狀物體上的A,B環切割出來,之後分別計算其重心座標,如圖32,並透過數學式9可計算出物件重心點座標值(Xc,Yc),其中P為pixel。
Step4.取得雙環二維資訊: 利用前述的Step3,將兩環重心點座標連線形成軸線段之直線方程式L,如圖33,且兩環二值化影像分別透過Canny邊緣檢測或形態學方法,取得兩環所有之邊緣點座標,帶入軸線段之直線方程式L通過邊緣點來計算出環 型標記兩端點座標位置A(Xa1,Ya1)、A(Xa2,Ya2)、B(Xb1,Yb1)、B(Xb2,Yb2),如圖34。
Step5.pixel轉換成mm單位: 將前述的Step4所取得的兩端點座標位置A(Xa1,Ya1)、A(Xb2,Ya2)、B(Xb1,Yb1)、B(Xb2,Yb2)單位從像素點Pixel轉換成mm單位(10mm=1cm),使用影像的解析度以及相機鏡頭大小換算。而本發明係使用Iphone5相機解析度為:3264x2448,相機鏡頭大小為:4.54x3.42mm,轉換係數為dy=4.54/3264,dx=3.42/2448。
Step6.進行三維姿態估測: 透過Step1到Step5方法可以推算出所需要的二維空間中的參數,再利用前述所推出的三維八象限定位系統,即可快速算出圓桿狀物體在3D空間的任意姿態。已知攝影機焦距(λ),環軸長(L),雙環重心軸長(L AB ),雙環兩端投射點(x a1,x a2,x b1,x b2)和(y a1,y a2,y b1,y b2),求桿狀物體三維定位參數{X A1Y A1Z A1αβγ},參數求解步驟如下: α為銳角。這邊列出了四種可能的α角,如表17(X-Y平面四象限α角)、表18(X-Z平面四象限β角與ZA1)、表19(Y-Z平面四象限γ角與ZA1)。
【表18】
帶入數學式11
便獲得XA1與YA1,進而獲得參考點座標點AA1(XAA1,YAA1,ZAA1)、如表20(圓桿狀體一至四象限AA1點座標推導結果):
綜上所述各實施說明,本發明係關於一種「基於二維影像及三維完整八象限定位之內視鏡手術器械追蹤方法」,且其構成結構未曾見於諸書刊 或公開使用,誠符合專利申請要件,懇請 鈞局明鑑,早日准予專利,至為感禱;需陳明者,以上所述乃是本專利申請案之具體實施例及所運用之技術原理,若依本專利申請案之構想所作之改變,其所產生之功能作用仍未超出說明書及圖式所涵蓋之精神時,均應在本專利申請案之範圍內,合予陳明。
a、b、c、d、e、f‧‧‧步驟

Claims (5)

  1. 一種運用三維完整八象限定位之內視鏡手術器械追蹤方法,係包含有:一雙環輪廓取出步驟,係將所輸入的投影於平面之手術器械影像擷取其第一環體(A)與第二環體(B)之平面影像;一影像座標轉換步驟,係設定前述第一環體及第二環體平面影像的中心為原點,並進行其座標轉換;一找出雙環重心的步驟,係分別計算出第一環體及第二環體的重心座標;一取得雙環二維資訊步驟,將該第一環體及第二環的重心座標連線形成軸線段之直線方程式(L)並透過邊緣點以計算出該第一環體的二中心端點座標A1、A2及該第二環體的二中心端點座標B1、B2;一像素轉換毫米單位步驟,係將前述四個端點座標A1、A2、B1、B2轉換成毫米單位;一進行三維姿態估測步驟,係使用上個步驟終的端點座標A1、A2、B1、B2及已知的攝影機焦距(λ),環軸長(L),雙環重心軸長(LAB),雙環兩端投射點(xa1、xa2、xb1、xb2)和(ya1、ya2、yb1、yb2)之參數,並配合使用三維八像限定位系統,進而求出求該手術器械的三維定位參數{XA1,YA1,ZA1,α,β,γ}之參數,進而可快速算該手術器械在空間的任意姿態;另該進行三維姿態估測步驟中,係透過下列步驟的數學式,以求得該桿狀物體三維定位參數{X A1Y A1Z A1αβγ}:步驟(1)先求出四象限的可能的α角:[數學式1] 步驟(2)再求出四象限的可能的β角與ZA1 步驟(3)再求出四象限的可能的γ角與ZA1 步驟(4)再求出XA1與YA1[數學式4] 如此,就可獲得參考點座標點AA1(XAA1,YAA1,ZAA1)
  2. 如請求項1所述運用三維完整八象限定位之內視鏡手術器械追蹤方法,其中該雙環輪廓取出步驟中,係定義該第一環體為(Ra,Ga,Ba)及該第二環體為(Rb,Gb,Bb),運用通道差的概念,將該第一環體、第二環體影像中的每個像素做三通道相減取絕對值,在對相減後的三通道分別設其定門檻值,並定義如下數學式:[數學式6]第一環體:|R-Ra|<TRa,|G-Ga|<TGa,|B-Ba|<TBa第二環體:|R-Rb|<TRb,|G-Gb|<TGb,|B-Bb|<TBb;其中該門檻值設定主要對於每種顏色對於光源的容忍度,假設以上條件成立則屬於該環體的像素。
  3. 如請求項1所述運用三維完整八象限定位之內視鏡手術器械追蹤方法,其中該影像座標轉換步驟中,係設一張MXN大小影像或視訊的影格為f(i,j),運用如下數學式得將f(i,j)轉成f(x,y):[數學式7]f(x,y)=(i-M/2,N/2-j)。
  4. 如請求項1所述運用三維完整八象限定位之內視鏡手術器械追蹤方法,其中該找出雙環重心步驟中,係透過下列數學式計算出其重心點座標值(Xc,Yc),其中P為像數(pixel):
  5. 如請求項1所述運用三維完整八象限定位之內視鏡手術器械追蹤方法,其中該取得雙環二維資訊驟中,係兩環二值化影像分別透過Canny邊緣檢測或形態學方法,取得兩環所有之邊緣點座標,帶入該軸線段之直線方程式(L)通過邊緣點來計算出環型標記兩端點座標位置A1(Xa1,Ya1)、A2(Xa2,Ya2)、B1(Xb1,Yb1)、B2(Xb2,Yb2)。
TW105135385A 2016-11-01 2016-11-01 基於二維影像及三維完整八象限定位之內視鏡手術器械追蹤方法 TWI594208B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
TW105135385A TWI594208B (zh) 2016-11-01 2016-11-01 基於二維影像及三維完整八象限定位之內視鏡手術器械追蹤方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
TW105135385A TWI594208B (zh) 2016-11-01 2016-11-01 基於二維影像及三維完整八象限定位之內視鏡手術器械追蹤方法

Publications (2)

Publication Number Publication Date
TWI594208B true TWI594208B (zh) 2017-08-01
TW201818347A TW201818347A (zh) 2018-05-16

Family

ID=60189218

Family Applications (1)

Application Number Title Priority Date Filing Date
TW105135385A TWI594208B (zh) 2016-11-01 2016-11-01 基於二維影像及三維完整八象限定位之內視鏡手術器械追蹤方法

Country Status (1)

Country Link
TW (1) TWI594208B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
TWI697317B (zh) * 2019-08-30 2020-07-01 國立中央大學 應用於與手術導航整合之混合實境系統之數位影像實境對位套件與方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20060235298A1 (en) * 2005-03-31 2006-10-19 Robert Kotmel Internal biopsy marking
TW200937350A (en) * 2008-02-22 2009-09-01 Univ Nat Cheng Kung Three-dimensional finger motion analysis system and method

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20060235298A1 (en) * 2005-03-31 2006-10-19 Robert Kotmel Internal biopsy marking
TW200937350A (en) * 2008-02-22 2009-09-01 Univ Nat Cheng Kung Three-dimensional finger motion analysis system and method

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
TWI697317B (zh) * 2019-08-30 2020-07-01 國立中央大學 應用於與手術導航整合之混合實境系統之數位影像實境對位套件與方法

Also Published As

Publication number Publication date
TW201818347A (zh) 2018-05-16

Similar Documents

Publication Publication Date Title
Mahmoud et al. Live tracking and dense reconstruction for handheld monocular endoscopy
Lin et al. Video‐based 3D reconstruction, laparoscope localization and deformation recovery for abdominal minimally invasive surgery: a survey
Schmalz et al. An endoscopic 3D scanner based on structured light
US10706610B2 (en) Method for displaying an object
CA2931529C (en) 3d corrected imaging
CN107167093B (zh) 一种激光线扫描与阴影莫尔的复合式测量系统及测量方法
JP7379704B2 (ja) 視覚化カメラ及び光学コヒーレンストモグラフィを統合するシステム及び方法
CN109357633B (zh) 三维扫描方法、装置、存储介质和处理器
JP6580761B1 (ja) 偏光ステレオカメラによる深度取得装置及びその方法
CN112184653B (zh) 一种基于双目内窥镜的病灶三维尺寸测量及显示方法
CN102831601A (zh) 基于联合相似性测度和自适应支持权重的立体匹配方法
WO2017199285A1 (ja) 画像処理装置及び画像処理方法
CN113008195B (zh) 一种基于空间点云的三维曲面距离测量方法及系统
CN110223355A (zh) 一种基于双重极线约束的特征标志点匹配方法
CN106236264A (zh) 基于光学跟踪和图像匹配的胃肠手术导航方法及系统
TWI594208B (zh) 基於二維影像及三維完整八象限定位之內視鏡手術器械追蹤方法
Shintaku et al. Three-dimensional surface models of autopsied human brains constructed from multiple photographs by photogrammetry
CN115311405A (zh) 一种双目内窥镜的三维重建方法
Malian et al. Medphos: A new photogrammetric system for medical measurement
CN108090930A (zh) 基于双目立体相机的障碍物视觉检测系统及方法
CN106303501B (zh) 基于图像稀疏特征匹配的立体图像重构方法及装置
Ettl Introductory review on ‘Flying Triangulation’: a motion-robust optical 3D measurement principle
CN109636903A (zh) 一种基于抖动的双目三维重建方法
Speidel et al. Interventional imaging: vision
US20210052146A1 (en) Systems and methods for selectively varying resolutions

Legal Events

Date Code Title Description
MM4A Annulment or lapse of patent due to non-payment of fees