劉天弼, 馮 瑞
(復旦大學 計算機科學技術學院, 上海 201203)
基因芯片與第二代DNA測序是兩種重要的高通量基因組學研究技術[1]. 要想準確快速地讀取高通量基因芯片的數據, 需要準確地拍攝到芯片圖像. 第二代DNA基因測序是通過對基因芯片拍攝圖片, 再對圖像進行掃描讀取堿基序列. 因此, 測序儀對基因芯片拍攝出高質量的圖片是測序過程中堿基識別的重要前提. 在拍攝圖像的分辨率已經確定的情況下, 單個芯片單元所占的像素數越少, 則整張圖像包含的單元越多, 測序的通量就越高. 然而, 要準確地讀取基因芯片上的信息,需要拍攝的芯片單元處于非常準確的位置, 對于僅占少量像素的芯片單元而言, 少量的位置偏差就會對信息的讀取造成很大的影響. 因此, 將拍攝的視場與芯片位置精準對齊, 對提高第二代基因測序技術的準確度和通量都具有十分重要的意義.
要對基因芯片拍攝出理想的圖片, 需要每個芯片單元都有固定的位置坐標. 攝像儀拍攝的圖片,反映的是當前視場的位置信息;而計算當前視場與芯片位置之間的偏差, 可視為當前拍攝圖片與理想情況下拍攝的圖片之間的圖像配準問題.
圖像配準技術種類繁多, 通常每一種配準技術都是針對某一具體場景而設計的, 不同的應用環境要綜合各方面的因素來選取相應的圖像配準技術.
傳統的圖像配準算法主要分為基于灰度的圖像配準和基于特征的圖像配準兩大類[2?3]. 基于灰度的圖像配準需要遍歷計算圖像中的灰度、梯度[4]、互信息[5?8]等, 從而確定兩幅圖像間的空間變換關系,因此運行速度較慢, 計算量較大. 自21世紀以來, 圖像配準主要采用基于特征的方法, 利用圖像的線特征、面特征、塊特征、邊緣特征等確定空間變換關系. 從圖像中提取關鍵點的技術歷史悠久且種類繁多, 如SIFT (Scale Invariant Feature Transform)[9]、SURF (Speeded Up Robust Features)[10]、SUSAN(Smallest Univalue Segment Assimilating Nucleus)角點[11]、MIC (Minimum Intensity Change)角點[12]、BRISK (Binary Robust Invariant Scalable Keypoints)檢測[13]等, 近些年來廣泛使用的特征點均有良好的尺度不變特征和抗噪聲性能. 此外, 可使用圖像中的直線段、邊緣曲線等線特征[14?16], 或閉合區域[17?18]等面特征進行配準運算. 基于特征的配準算法自由度較高, 但是快速的特征提取一般僅利用圖像中的小部分信息, 會導致精度降低; 具有良好區分度的特征描述往往需要復雜的計算, 反而導致效率降低.
目前,大多數關于圖像配準的研究都涉及深度學習. 在過去的幾年中, 深度學習使計算機視覺任務具有先進的性能, 如圖像分類、目標檢測和語義分割等. 使用深度學習做圖像配準需要準備合適的數據集進行訓練, 以期得到優秀的配準算法模型. 從功能上看, 基于深度學習的圖像配準可以分為兩大類[19]:①利用深度卷積神經網絡估計兩幅圖像的相似性度量; ②直接利用深度回歸網絡預測轉換參數. 前者進行相似性度量, 仍然需要傳統數字圖像處理進行迭代做配準優化, 運算耗時較久; 后者是通過大量學習獲得配準映射能力, 但只能用于非剛性配準, 無法達到高精度要求.
根據深度學習的種類劃分, 此類配準方法也可以分為基于監督學習的配準與基于非監督學習的配準兩大類. 基于監督學習的配準, 是通過兩張作為配準對的圖片輸入得到映射向量, 即真實變形場.其關鍵在于構建數據集時需提供與配準對相對應的真實變形場作為Ground Truth, 作為數據集標簽[20?22]. 基于非監督學習的配準方法在訓練時只需要提供配準對, 不需要標簽[23?25], 完成訓練后可通過一張圖像計算獲得配準圖像.
然而, 基于深度學習本身的特點, 其用于高精度配準的場合尚存在諸多問題: ①深度學習依賴于龐大的訓練數據集; ②算法訓練成果是基于機器學習經驗, 無法滿足高精準的配準; ③單一場景下的圖像配準對泛化能力要求不高, 而深度神經網絡計算量龐大, 面對高分辨率圖像的復雜計算反而成為負擔.
本文通過在基因芯片上的特定位置設置標記, 并確定在理想狀態下拍攝的圖片中所有標記的位置坐標, 從而達到圖像配準的目標;在得到實時拍攝的圖像后, 通過對標記的位置進行捕捉可初步將視場中心與標準圖像的中心對齊; 然后對標記上的關鍵點進行捕捉, 得到一系列特征點的坐標; 再對關鍵點的映射關系進行擬合, 計算出高精度的坐標和角度偏差結果.
本文設計的用于準確捕捉圖像特征的新穎卷積方法, 具有良好的抗干擾能力, 其算法簡捷, 且易于移植使用硬件加速; 同時使用圖像全范圍多處采樣進行擬合的計算方法, 保證了圖像配準達到較高的精確度. 經實踐驗證和仿真實驗分析, 本文算法能取得精確的偏差計算結果, 并且在應用場景下具有很好的魯棒性.
基因芯片上布滿了攜帶DNA信息的探針對, 在布滿探針對的區域內, 開辟一塊不放置探針對的獨立區域, 即可看作芯片內的標記. 標記的作用在拍攝的芯片圖像中體現出來, 每個探針均攜帶熒光信號, 在圖像上可看作是1個帶有灰度信息的單元(cell), 而標記處不放置探針, 所以始終呈現黑色狀態.
在芯片的水平或豎直方向上, 設置1條直線型的標記, 稱為track標記; 1條水平track標記與1條豎直track標記相交于一點, 稱為cross標記.
先驗標記的作用就是人為地在圖像中設置關鍵特征, 作為圖像配準過程中實現映射的特征采樣點. 幾種先驗track標記設計示例如圖1所示.

圖 1 3種track標記設計示意圖Fig. 1 Schematic diagram of three track mark designs
圖1中黑色線段即track標記, 而cross標記是依賴于track標記的相交呈現的. 設計track標記遵循以下準則.
(1) 所有track標記寬度統一.
(2) 單條track標記必須為水平方向或者豎直方向.
(3) track標記在視場范圍內的分布應呈左右、上下軸對稱.
(4) cross標記的分布能夠體現出對圖像的采樣比較均衡.
在芯片上設定好先驗目標之后, 根據芯片成像的分辨率確定全部標記的位置、寬度及長度的標準數據, 作為標準位置圖像的特征. 這樣就將視場對芯片位置做校準的問題, 轉化為當前視場拍攝到的圖片與標準圖像的配準問題.
2.2.1 track標記的捕捉
視場中, 帶有亮度信息的cell隨機分布, 而track標記則全部呈現為黑色, 可通過提取其灰度特征進行捕捉.
以豎直的track標記為例, 設其成像寬度為w像素, 長度為L像素. 若拍攝圖像的高為H像素, 寬為W像素, 則對于任意一種track標記設計方式, 都有L≤H. 如圖2所示, 在豎直方向上截取一段含有track標記的圖像, 其中track標記的數量為N個.

圖 2 于豎直方向截取圖像做水平卷積示意圖Fig. 2 Schematic diagram of the image in the vertical direction for horizontal convolution
取卷積核大小為l×w, 其中,l為豎直方向長度,w為水平方向寬度, 且l<L. 在水平方向上對當前截取的l×W范圍上做卷積操作. 卷積核有以下兩種設計方式.
(1)均型: 卷積核所有元素值等于1, 相當于卷積范圍內覆蓋的像素響應均等.
(2)谷型: 卷積核所有元素值為正, 均值等于1. 卷積核行向量的左右邊緣元素值最大, 中央最小,列方向上是行向量的復制,相當于卷積覆蓋的像素以中間列的響應最小.
均型卷積適用于對抗拍攝場景中噪聲的情況; 谷型卷積更適合對抗cell間光線存在串擾的情況.關于卷積核設計對算法魯棒性的影響, 本文將在實驗部分做進一步分析.
設卷積核為k, 截取l×W范圍的圖像為I′, 執行卷積運算

即

執行卷積運算時, 為捕捉track標記的中心位置, 圖像需要在左右方向做像素填充(padding), 其上下方向不做. padding像素的灰度為純白色. padding數量為

卷積結果為W個元素的一維張量, 即

卷積結果r的索引j與截取圖像的x坐標相對應. 將r的元素做升序(Ascend, Asc)排列, 前N個元素對應的索引即為搜索到的track標記的x坐標

再將此N個坐標重新排序, 則截取圖像范圍內從左至右的track標記的x坐標就全部搜索完成, 即

豎直方向上捕捉track標記的位置, 使用上述卷積核的轉置kT, 在豎直方向上執行卷積運算, 即可通過類似的算法得到Ytrack.
2.2.2 cross標記的捕捉
如圖3所示, 在cross標記范圍的橫縱方向上執行track標記捕捉運算, 即可得到cross標記上的交點坐標

其中,X是在cross標記位置上捕捉的x坐標集合,Y是在cross標記位置上捕捉的y坐標集合.

圖 3 卷積運算捕捉cross標記示意圖Fig. 3 Schematic diagram of the convolution operation to capture cross marks
根據當前視場拍攝到的圖片與標準圖像的配準, 即可計算出拍攝的圖片與標準位置之間的坐標偏差和角度偏差. 提高配準精度的基本依據是準確的關鍵點位置坐標, 提高配準精度的基本方法是通過多采樣做最優擬合.
2.3.1 坐標原點對齊
在做位置偏差的精確計算之前, 需要圖像的中心位置與標準位置對齊, 將中心位置視作坐標原點.其原因有如下兩點.
(1)便于分析特征點配準的映射關系.
(2)尋找cross標記既定范圍, 便于捕捉cross標記.
基于2.1節先驗track標記的設計準則(3), 視場范圍內的track標記呈軸對稱分布, 則分別計算相互對稱的2個track標記的坐標偏差, 其均值能夠消除旋轉對坐標偏差的影響, 獲得中心點的坐標偏差.
豎直track標記有u行, 每行v個; 水平track標記有k列, 每列l個. 特別是, 對于貫穿型的track標記, 可將豎直track線合理地劃分為u段, 使用相應的卷積核進行捕捉, 仍看作u×v的track分布;同理, 將水平track線合理地劃分為k段, 仍看作k×l的track分布.

通過坐標變換將當前拍攝圖像的中心與標準位置對齊, 即旋轉中心與標準中心位置對齊. 因位置坐標是以整數表示, 因此中心位置對齊的誤差在1像素左右.
2.3.2 旋轉角度與錯位擬合
根據先驗設計, 視場內存在M×N個cross標記, 通過對cross標記的捕捉, 可得到M×N個交點的位置坐標.
假設當前拍攝圖像與標準位置的角度偏差為φ, 因圖像具有旋轉不變特征, 所以圖像坐標系原點位置可任意指定且不影響旋轉角度的計算. 取圖像中任意一點研究其旋轉的位置特性, 如圖4所示,像素點A經旋轉一定角度之后, 處于A′的位置.

圖 4 單個像素錯位映射Fig. 4 Mapping of single pixel misalignment

前文根據cross標記已捕捉到M×N個交點坐標, 是符合式(15)的解. 對點坐標進行處理, 得

根據式(15)可得


兩幅圖像的配準結果是否優秀, 目前并沒有一個很明確的標準. 但是本文具體討論的配準是一幅圖像向標準位置的配準, 針對位置校準問題的評價標準較為明確, 即坐標與角度偏差的準確度.
實驗中, 使用不同型號的基因測序儀進行拍攝: 一種拍攝的圖像為矩形, 分辨率為(2 560 × 2 160)像素; 另一種拍攝的圖像為正方形, 成像分辨率為(5 012 × 5 012)像素. 使用不同型號的基因芯片, 芯片上的探針數量不同, track標記均使用非均勻貫穿型設計.
基因芯片上的探針經拍攝后形成單元cell, 其對應的分辨率為(a×a)像素. 可認為位置校準的結果誤差符合正態分布, 通過多次校準, 取cell位置處灰度最大值的校準結果為標準位, 則可估算出坐標校準的誤差期望, 并以圖像邊緣位置的多個cell位置誤差推測出角度配準的準確度. 不同芯片在不同儀器上進行拍攝和配準, 位置與角度配準誤差如表1所示, 其中, 儀器型號用分辨率表示. 基于實際拍攝情況, cell周圍的像素會受到光線串擾的影響, 因此在估算時對cell周圍的像素保留了灰度串擾值. 實際使用中, 位置誤差在0.5像素以下即可視為鏡頭與基因芯片精準對齊.

表 1 位置與角度配準誤差Tab. 1 Position and angle deviation in registration
從實驗結果可以看到, 在位置校準能夠保障精度的情況下, 可以提高研究對象的拍攝分辨率, 還可以減少單個信息單元占用的像素數, 能夠大大提升計算通量.
一般來說, 當出現track標記不在視場范圍內、track標記錯位對齊的現象, 已屬于設備故障問題,本算法不展開故障檢測問題的討論. 對于角度偏差的適用范圍, 與標記捕捉時的卷積核設計有關.
3.2.1 卷積核長度與容錯角度
無論是對track標記還是cross標記的捕捉, 在搜索x坐標或y坐標時, 卷積運算的作用是完全一致的, 即通過灰度特征提取得到track標記位置. 因此, 卷積核應設計成寬度與track線的寬度一致.
視場中cell的邊長等于track標記的寬度, 以均值為最高亮度的25%的灰度隨機填充cell. 以track標記中心位置為圓心逐步增加圖像的旋轉角度, 當通過卷積捕捉到track標記位置坐標偏差時,記當前角度為卷積核的容錯角度.
越大的卷積核其抗干擾能力越強, 所以設計卷積核時應當使用較大的長度. 然而圖5的結果表明,越長的卷積核的容錯角度越小. 在實際應用中, 拍攝系統配準前的角度偏差很容易控制在1°以下, 因此本文討論的卷積運算能夠很好地應用在拍攝系統中.

圖 5 卷積核長度與容錯角度的關系Fig. 5 Relationship between the length of the convolution kernel and the angle of fault tolerance
3.2.2 卷積切面與干擾
圖像配準的干擾主要有兩個: 噪聲和光線串擾. 噪聲在拍攝環境中不可避免, 而拍攝對象中的帶有亮度的信息單元在成像時會對周圍的像素產生光線串擾.
對先驗標記進行捕捉的卷積核, 其寬度等于track標記的寬度. 為方便起見, 取豎直方向track標記對應的卷積核為研究對象, 卷積核內元素是相同的行向量的堆疊. 因此, 取行向量作為卷積核的切面, 研究其形狀與魯棒性的關系. 同樣地, 對應于水平方向的track標記, 將列視為切面.
選取如圖6所示的幾種切面形狀的卷積核, 分別是均型、峰型、谷型. 圖6 中, 縱坐標表示卷積元素的值(切面不同位置的數值), 橫坐標表示卷積切面元素的位置.
在沒有噪聲和串擾的理想情況下, track線上的卷積結果應當為0, 其他位置的卷積期望為一個正值. 這種卷積特性可以很好地區別track標記與業務區域, 并準確定位track的中心位置.
為了研究卷積核設計對track的捕捉能力, 需要研究track附近的卷積結果之間的區分能力. 在卷積運算接近track標記并經過整個標記直到離開, 這段范圍的卷積值應當是先遞減再遞增, 出現一個極小值, 極小值出現的位置就是捕捉到的track標記位置.
圖7為通過卷積捕捉track標注的示意圖. 從圖7中可以看到, 普通場景下的業務區域卷積結果是明顯高于track標記上的, 在track標記位置附近呈現凸函數形狀, 因此track捕捉能力僅需要考察track線及其左右相鄰范圍的卷積值的梯度絕對值, 其值越大表明對位置的區分能力越強.下面通過仿真, 分析理想狀態下的卷積運算的區分能力,并進一步分析不同卷積核面對噪聲和光線串擾的抗干擾能力.

圖 6 不同的卷積核切面Fig. 6 Cross sections of various convolution kernels

圖 7 通過卷積捕捉track標記Fig. 7 Capturing a track mark by convolution
(1) 理想狀態下卷積核的區分能力
取卷積核寬度為9像素, 在沒有任何干擾的理想狀態下, 業務區域以25%的亮度單元隨機分布, 各卷積核在此狀態下的區分度如圖8所示. 關注區域在[–5, 5]上, 此區域之外的梯度會逐步趨向0附近.
理想狀態下的卷積以均值型的最為理想, 在任何位置處都有穩定的區分能力. 峰型卷積核在track標記左右邊緣處區分能力最強. 谷型卷積核在track中心位置處區分能力最強. track捕捉能力只體現在關注區域中, 因此紅色標識的峰型卷積核的魯棒性最差.
(2) 噪聲環境下卷積核的區分能力
對理想環境增加不同程度的白噪聲, 以噪聲均值N0區別噪聲的嚴重程度, 同樣對上文列舉的5種卷積核進行對比. 對比后發現5種卷積核對track的區分能力有所變化, 結果如圖9所示.

圖 8 不同卷積核在理想狀態下的區分度Fig. 8 Differentiation of various convolution kernels under ideal conditions

圖 9 不同卷積核在噪聲環境下的區分度Fig. 9 Differentiation of various convolution kernels in a noisy environment
從圖9中可以看出, 隨著噪聲能量的增加, 所有的卷積核對track標記的區分能力都在下降. 在[–5, 5]之間, 梯度最低值代表該卷積核區分能力最差值. 當最差區分能力低到一定程度時, 應對光線的極端分布情況的能力就會變差. 因此,從圖9中還可以看出, 峰型卷積核并不適用于track標記的捕捉, 而均型卷積核保持著穩定的區分能力.
(3)光線串擾環境下卷積核的區分能力
增加了40%的光線串擾之后, 5種卷積核對track的區分能力有所變化, 結果如圖10所示. 圖10表明, 紅色曲線代表的峰型卷積核區分能力變得更差, 因此峰型卷積核完全不適合在捕捉track標記時使用; 在有光線串擾的情況下, 均值型卷積核的區分能力也有所下降, 而谷型卷積核則表現出很好的區分能力, 并且比均值型卷積核具有更好的穩定性.

圖 10 不同卷積核在串擾環境下的區分度Fig. 10 Differentiation of various convolution kernels in a light crosstalk environment
鑒于谷型卷積核對光線串擾具有優秀的抗干擾能力, 因此本文對谷型形狀與串擾之間的相應關系做了進一步觀察. 如圖11所示, 設計5種谷型深度不同的卷積切面, 其卷積核各元素之和與均值型卷積核保持一致. 圖11中, 縱坐標表示卷積元素的值(切面不同位置的數值), 橫坐標表示卷積切面元素的位置. 選取10%、20%、30%、40% 這4個等級的串擾, 5種卷積核對track標記的區分能力如圖12所示.

圖 11 5種抗串擾卷積核切面Fig. 11 Cross sections of five anti?crosstalk convolution kernels

圖 12 5種抗串擾效果Fig. 12 Anti?crosstalk effects of five convolution kernels
從圖12可看到, 10%串擾時卷積核(2)(圖12( a))的抗串擾性能最好, 20%串擾時卷積核(3)(圖12 (b))的抗串擾性能最好, 30%串擾時卷積核(5)(圖12( c))的抗串擾性能最好, 串擾超過40%之后, 上述5種谷型卷積核性能均變得很差. 因此可以得出結論: 谷型卷積核的“谷深”越深, 應對光線串擾的能力越強; 谷型卷積核的設計需對應具體的串擾強度, 從頻域角度計算其濾波性能.
綜上所述, 從實驗結果可知,卷積操作對捕捉track標記具有很高的準確性和很強的魯棒性, 在白噪聲的環境下, 均型卷積核始終保持了較好的track標記區分能力, 而在光線串擾的環境中, 需要使用谷型卷積核以保證其區分能力.
本文針對基因測序儀需要將鏡頭與基因芯片精準對齊的問題, 提出了一種基于先驗標記的圖像配準方法: 計算當前視場與理想位置的偏差. 該方法通過在基因芯片上人工設計標記, 再分析所拍攝圖像的灰度特征, 利用卷積運算對圖像特征的良好提取能力捕捉標記, 使用擬合的方式計算偏差. 從實際應用和仿真分析可知, 本文算法對位置偏差的計算結果具有良好的精確度, 誤差限制在0.5像素以下,并且有良好的抗干擾能力, 在第二代基因測序技術提高測序準確度和通量的需求下具有很好的應用前景.