GPU 加速下的影象視覺
越來越多的移動計算裝置都開始攜帶照相機鏡頭，這對於攝影界來說是一個好事情，不僅如此攜帶鏡頭也為這些裝置提供了更多的可能性。除了最基本的拍攝功能，結合合適的軟體這些更為強大的硬體裝置可以像人腦一樣理解它看到了什麼。
僅僅具備一點點的理解能力就可以催生一些非常強大的應用，比如說條形碼識別，文件識別和成像，手寫文字的轉化，實時影象防抖，增強現實等。隨著處理能力變得更加強大，鏡頭保真程度更高，演算法效率更好，機器視覺 (machine vision) 這個技術將會解決更加重大的問題。
有些人認為機器視覺是個非常複雜的領域，是程式設計師們的日常工作中絕不會遇到的。我認為這種觀點是不正確的。我發起了一個開源專案
GPU 是一種理想的處理圖片和視訊的裝置，因為它是專門為並行處理大量資料而生的，圖片和視訊中的每一幀都包含大量的畫素資料。在某些情況下 GPU 處理圖片的速度可以是 CPU 的成千上百倍。
在我開發 GPUImage 的過程中我學到了一件事情，那就是即使是圖片處理這樣看上去很複雜的工作依然可以分解為一個個更小更簡單的部分。這篇文章裡我想將一些機器視覺中常見的過程分解開來，並且展示如何在現代的 GPU 裝置上讓這些過程執行地更快。
以下的每一步在 GPUImage 中都有完整的實現，你可以下載包含了 OS X 和 iOS 版本的示例工程 FilterShowcase，在其中體驗一下各個功能。此外，這些功能都有基於 CPU (有些使用了 GPU 加速) 的實現，這些實現是基於 OpenCV 庫的，在
索貝爾 (Sobel) 邊界探測
我將要描述的第一種操作事實上在濾鏡方面的應用比機器視覺方面更多，但是從這個操作講起是比較合適的。索貝爾邊界探測用於探測一張圖片中邊界的出現位置，邊界是指由明轉暗的突然變化或者反過來由暗轉明的區域[1]。在被處理的圖片中一個畫素的亮度反映了這個畫素周圍邊界的強度。
下面是一個例子，我們來看看同一張圖片在進行索貝爾邊界探測之前和之後：
正如我上面提到的，這項技術通常用來實現一些視覺效果。如果在上面的圖片中將顏色進行反轉，最明顯的邊界用黑色代表而不是白色，那麼我們就得到了一張類似鉛筆素描效果的圖片。
那麼這些邊界是如何被探測出來的？第一步這張彩色圖片需要減薄成一張亮度 (灰階) 圖。Janie Clayton 在她的文章中解釋了這一步是如何在一個片斷著色器 (fragment shader) 中完成的。簡單地說這個過程就是將每個畫素的紅綠藍部分加權合為一個代表這個畫素亮度的值。
有的視訊裝置和相機提供的是 YUV 格式的圖片，而不是 RGB 格式。YUV 這種色彩格式已經將亮度資訊 (Y) 和色度資訊 (UV) 分離，所以如果原圖片是這種格式，顏色轉換這個步驟就可以省略，直接用其中亮度的部分就可以。
圖片一旦減薄到僅剩亮度資訊，一個畫素周圍的邊界強度就可以由它周圍 3*3 個臨近畫素計算而得。在一組畫素上進行圖片處理的計算過程涉及到一個叫做卷積矩陣 (參考：convolution matrix)) 的東西。卷積矩陣是一個由權重資料組成的矩陣，中心畫素周圍畫素的亮度乘以這些權重然後再相加就能得到中心畫素的轉化後數值。
圖片上的每一個畫素都要與這個矩陣計算出一個數值。在處理的過程中畫素的處理順序是無關緊要的，所以這種計算很容易並行執行。因此，這個計算過程可以通過一個片斷著色器的方式執行在可程式設計的 GPU 上，來極大地提高處理效率。正如在 Janie 的文章中所提到的，片斷著色器是一些 C 語言風格的程式，執行在 GPU 上可以進行一些非常快速的圖片處理。
下面這個是索貝爾運算元的水平處理矩陣：
| −1 | 0 | +1 |
| −2 | 0 | +2 |
| −1 | 0 | +1 |
為了進行某一個畫素的計算，每一個臨近畫素的亮度資訊都要讀取出來。如果要處理的圖片已經被轉化為灰階圖，亮度可以從紅綠藍任意顏色通道中抽樣。臨近畫素的亮度乘以矩陣中對應的權重，然後加到最終值裡去。
在一個方向上尋找邊界的過程是這樣的：轉化之後對比一個畫素左右兩邊畫素的亮度差。如果當前這個畫素左右兩邊的畫素亮度相同也就是說在圖片上是一個柔和的過度，它們的亮度值和正負權重會相互抵消，於是這個區域不會被判定為邊界。如果左邊畫素和右邊畫素的亮度差別很大也就是說是一個邊界，用其中一個亮度減去另一個，這種差異越大這個邊界就越強 (越明顯)。
索貝爾過程有兩個步驟，首先是水平矩陣進行，同時一個垂直矩陣也會進行，這個垂直矩陣中的權重如下
| −1 | −2 | −1 |
| 0 | 0 | 0 |
| +1 | +2 | +1 |
兩個方向轉化後的加權和會被記錄下來，它們的平方和的平方根也會被計算出來。之所以要進行平方是因為計算出來的值可能是正值也可能是負值，但是我們需要的是值的量級而不關心它們的正負。有一個好用內建的 GLSL 函式能夠幫助我們快速完成這個過程。
最終計算出來的這個值會用來作為輸出圖片中畫素的亮度。因為索貝爾運算元會突出顯示兩邊畫素亮度的不同的地方，所以圖片中由明轉暗或者相反的突然轉變會成為結果中明亮的畫素。
索貝爾邊界探測有一些相似的變體，例如普里維特 (Prewitt) 邊界探測[2]。普里維特邊界探測會在橫向豎向矩陣中使用不同的權重，但是它們運作的基本過程是一樣的。
作為索貝爾邊界探測如何用程式碼實現的一個例子，下面是用 OpenGL ES 進行索貝爾邊界探測的片斷著色器：
precision mediump float;
varying vec2 textureCoordinate;
varying vec2 leftTextureCoordinate;
varying vec2 rightTextureCoordinate;
varying vec2 topTextureCoordinate;
varying vec2 topLeftTextureCoordinate;
varying vec2 topRightTextureCoordinate;
varying vec2 bottomTextureCoordinate;
varying vec2 bottomLeftTextureCoordinate;
varying vec2 bottomRightTextureCoordinate;
uniform sampler2D inputImageTexture;
void main()
{
float bottomLeftIntensity = texture2D(inputImageTexture, bottomLeftTextureCoordinate).r;
float topRightIntensity = texture2D(inputImageTexture, topRightTextureCoordinate).r;
float topLeftIntensity = texture2D(inputImageTexture, topLeftTextureCoordinate).r;
float bottomRightIntensity = texture2D(inputImageTexture, bottomRightTextureCoordinate).r;
float leftIntensity = texture2D(inputImageTexture, leftTextureCoordinate).r;
float rightIntensity = texture2D(inputImageTexture, rightTextureCoordinate).r;
float bottomIntensity = texture2D(inputImageTexture, bottomTextureCoordinate).r;
float topIntensity = texture2D(inputImageTexture, topTextureCoordinate).r;
float h = -bottomLeftIntensity - 2.0 * leftIntensity - topLeftIntensity + bottomRightIntensity + 2.0 * rightIntensity + topRightIntensity;
float v = -topLeftIntensity - 2.0 * topIntensity - topRightIntensity + bottomLeftIntensity + 2.0 * bottomIntensity + bottomRightIntensity;
float mag = length(vec2(h, v));
gl_FragColor = vec4(vec3(mag), 1.0);
}
上面這段著色器中中心畫素周圍的畫素都有使用者定義的名稱，是由一個自定義的頂點著色器提供的，這麼做可以優化減少對移動裝置環境的依賴。從 3*3 網格中抽樣出這些命名了的畫素，然後用自定義的程式碼來進行橫向和縱向索貝爾探測。為簡化計算權重為 0 的部分會被忽略。GLSL 函式 length() 計算出水平和垂直矩陣轉化後值的平方和的平方根。然後這個代表量級的值會被拷貝進輸出畫素的紅綠藍通道中，這樣就可以用來代表邊界的明顯程度。
坎尼 (Canny) 邊界探測
索貝爾邊界探測可以給你一張圖片邊界強度的直觀印象，但是並不能明確地說明某一個畫素是否是一個邊界。如果要判斷一個畫素是否是一個邊界，你要設定一個類似閾值的東西，亮度高於這個閾值的畫素會被判定為邊界的一部分。然而這樣並不是最理想的，因為這樣的做法判定出的邊界可能會有好幾個畫素寬，並且不同的圖片適合的閾值不同。
這裡你更需要一種叫做坎尼邊界探測[3]的邊界探測方法。坎尼邊界探測可以在一張圖片中探測出連貫的只有一畫素寬的邊界：
坎尼邊界探測包含了幾個步驟。和索貝爾邊界探測以及其他我們接下來將要討論的方法一樣，在進行邊界探測之前首先圖片需要轉化成亮度圖。一旦轉化為灰階亮度圖緊接著進行一點點的高斯模糊，這麼做是為了降低感測器噪音對邊界探測的影響。
一旦圖片已經準備好了，邊界探測就可以開始進行。這裡的 GPU 加速過程原本是在 Ensor 和 Hall 的文章 "GPU-based Image Analysis on Mobile Devices" [4]中所描述的。
首先，一個給定畫素的邊界強度和邊界梯度要確定下來。邊界梯度是指亮度發生變化最大的方向，也是邊界延伸方向的垂直方向。
為了尋找邊界梯度，我們要用到上一章中的索貝爾矩陣。索貝爾轉化得到的橫豎值加合後就是邊界梯度的強度，這個值會編碼進輸出畫素的紅色通道。然後橫向豎向索貝爾結果值會與八個方向 (對應一箇中心畫素周圍的八個畫素) 中的一個結合起來，一個方向上 X 部分值會作為輸出畫素的綠色通道值，Y 部分則作為藍色通道值。
這個方法使用的著色器和索貝爾邊界探測使用的類似，只是最後一個計算步驟用下面這段程式碼：
vec2 gradientDirection;
gradientDirection.x = -bottomLeftIntensity - 2.0 * leftIntensity - topLeftIntensity + bottomRightIntensity + 2.0 * rightIntensity + topRightIntensity;
gradientDirection.y = -topLeftIntensity - 2.0 * topIntensity - topRightIntensity + bottomLeftIntensity + 2.0 * bottomIntensity + bottomRightIntensity;
float gradientMagnitude = length(gradientDirection);
vec2 normalizedDirection = normalize(gradientDirection);
normalizedDirection = sign(normalizedDirection) * floor(abs(normalizedDirection) + 0.617316); // Offset by 1-sin(pi/8) to set to 0 if near axis, 1 if away
normalizedDirection = (normalizedDirection + 1.0) * 0.5; // Place -1.0 - 1.0 within 0 - 1.0
gl_FragColor = vec4(gradientMagnitude, normalizedDirection.x, normalizedDirection.y, 1.0);
為確保坎尼邊界一畫素寬，只有邊界中強度最高的部分會被保留下來。因此，我們需要在每一個切面邊界梯度的寬度之內尋找最大值。
這就是我們在上一步中算出的梯度方向起作用的地方。對每一個畫素，我們根據梯度值向前和向後取出最近的相鄰畫素，然後對比他們的梯度強度 (邊界明顯程度)。如果當前畫素的梯度強度高於梯度方向前後的畫素我們就保留當前畫素。如果當前畫素的梯度強度低於任何一個臨近畫素，我們就不再考慮這個畫素並且將他變為黑色。
執行這個步驟的著色器如下：
precision mediump float;
varying highp vec2 textureCoordinate;
uniform sampler2D inputImageTexture;
uniform highp float texelWidth;
uniform highp float texelHeight;
uniform mediump float upperThreshold;
uniform mediump float lowerThreshold;
void main()
{
vec3 currentGradientAndDirection = texture2D(inputImageTexture, textureCoordinate).rgb;
vec2 gradientDirection = ((currentGradientAndDirection.gb * 2.0) - 1.0) * vec2(texelWidth, texelHeight);
float firstSampledGradientMagnitude = texture2D(inputImageTexture, textureCoordinate + gradientDirection).r;
float secondSampledGradientMagnitude = texture2D(inputImageTexture, textureCoordinate - gradientDirection).r;
float multiplier = step(firstSampledGradientMagnitude, currentGradientAndDirection.r);
multiplier = multiplier * step(secondSampledGradientMagnitude, currentGradientAndDirection.r);
float thresholdCompliance = smoothstep(lowerThreshold, upperThreshold, currentGradientAndDirection.r);
multiplier = multiplier * thresholdCompliance;
gl_FragColor = vec4(multiplier, multiplier, multiplier, 1.0);
}
其中 texelWidth 和 texelHeight 是要處理的圖片中臨近畫素之間的距離，lowerThreshold 和 upperThreshold 分別設定了我們預期的邊界強度上下限。
在坎尼邊界探測的最後一步，邊界上出現畫素間間隔的地方要被填充，出現間隔是因為有一些點不在閾值範圍之內或者是因為非最大值轉化沒有起作用。這一步會完善邊界使邊界連續起來。
在最後一步中需要考慮一箇中心畫素周圍的所有畫素。如果這個中心畫素是最大值，上一步中非最大值轉化就不會影響它，它依然是白色。如果它不是最大值，就會變成黑色。對於中間的灰色畫素，會考察它周圍畫素的資訊。凡是與超過一個白色畫素挨著的都會變為白色，相反就會變成黑色。這樣就可以將邊界分離的部分接合起來。
正如你所看到的，坎尼邊界探測會比索貝爾邊界探測更復雜一些，但是它會探測出一條物體邊界的乾淨線條。這是線條探測、輪廓探測或者其他圖片分析很好的起點。同時也可以被用來生成一些有趣的美學效果。
哈里斯 (Harris) 邊角探測
雖然利用上一章中的邊界探測技術我們可以獲取關於圖片邊界的資訊，我們會得到一張可以直觀觀察到邊界所在位置的圖片，但是並沒有更高層面有關圖片中所展示內容的資訊。為了得到這些資訊，我們需要一個可以處理場景中的畫素然後返回場景中所展示內容的描述性資訊的演算法。
進行物體探測和匹配時一個常見的出發點是特徵探測。特徵是指一個場景中具有特殊意義的點，這些點可以唯一的區分出一些結構或者物體。由於邊角的出現往往意味著亮度或者顏色的突然變化，所以邊角常常會作為特徵的一種。
在 Harris 和 Stephens 的文章 "A Combined Corner and Edge Detector."[5] 中他們提出一個邊角探測的方法。這個命名為哈里斯邊角探測的方法採用了一個多步驟的方法來探測場景中的邊角。
像我們已經討論過的其他方法一樣，圖片首先需要減薄到只剩亮度資訊。通過索貝爾矩陣，普里維特矩陣或者其他相關的矩陣計算出一個畫素 X 和 Y 方向上的梯度值，計算出的值並不會合併為邊界的量級。而是將 X 梯度傳入紅色部分，Y 梯度傳入綠色部分，X 與 Y 梯度的乘積傳入藍色部分。
然後對上述計算結果進行一個高斯模糊。從模糊後的照片中取出紅綠藍部分中的編碼過的值，並將值帶入一個計算畫素是邊角點可能性的公式：
R = Ix2 × Iy2 − Ixy × Ixy − k × (Ix2 + Iy2)2
其中 Ix 是 X 方向梯度值 (模糊後圖片中紅色部分)，Iy 是 Y 梯度值 (綠色部分)，Ixy 是 XY 值的乘積 (藍色部分)，k 是一個靈敏性常數，R 是計算出來的這個畫素是邊角的確定程度。Shi，Tomasi[6] 和 Noble[7] 提出過這種計算的另一種實現方法但是結果其實是十分接近的。
在公式中你可以會覺得頭兩項會抵消掉。但這就是前面高斯模糊那一步起作用的地方。通過在一些畫素上分別模糊 X、Y 和 XY 的乘積，在邊角附近就會出現可以被探測到的差異。
我們從 Stack Exchange 訊號處理分站中的一個問題中取來一張測試圖片：
經過前面的計算過程得到的結果如下圖：
為了找出邊角準確的位置，我們需要選出極點 (一個區域內亮度最高的地方)。這裡需要使用一個非最大值轉化。和我們在坎尼邊界探測中所做的一樣，我們要考察一箇中心畫素周圍的臨近畫素 (從一個畫素半徑開始，半徑可以擴大)，只有當中心畫素的亮度高於它所有臨近畫素時才保留他，否則就將這個畫素變為黑色。這樣一來最後留下的就應該是一片區域中亮度最高的畫素，也就是最可能是邊角的地方。
通過這個過程，我們現在可以從圖片中看到任意不是黑色的畫素都是一個邊角所在的位置：
目前我是使用 CPU 來進行點的提取，這可能會是邊角探測的一個瓶頸，不過在 GPU 上使用柱狀圖金字塔[8]可能會加速這個過程。
哈里斯邊角探測只是在場景中尋找邊角的方法之一。"Machine learning for high-speed corner detection,"[9] 中 Edward Rosten 的 FAST 邊角探測方法是另一個性能更好的方法，甚至可能超越基於 GPU 的哈里斯探測。
霍夫 (Hough) 變換線段探測
筆直的線段是另一種我們會在一個場景需要探測的常見的特徵。尋找筆直的線段可以幫助應用進行文件掃描和條形碼讀取。然而，傳統的線段探測方法並不適合在 GPU 上實現，特別是在移動裝置的 GPU 上。
許多線段探測過程都基於霍夫變換，這是一項將真實世界笛卡爾直角座標空間中的點轉化到另一個座標空間中去的技術。轉化之後在另一個座標空間中進行計算，計算的結果又轉化回正常空間代表線段的位置或者其他特徵資訊。不幸的是，許多已經提出的計算方法都不適合在 GPU 上執行，因為它們在特性上就不太可能充分地並行執行，並且都需要大量的數學計算，比如在每個畫素上進行三角函式計算。
2011年，Dubská 等人 [10] [11] 提出了一種更簡單並更有效的座標空間轉換方法和分析方法，這種方法更合適在 GPU 上執行。他們的方法依賴與一個叫做平行座標空間的概念，聽上去很抽象但是我會展示出它其實很容易理解。
我們首先選擇一條線段和線段上的三個點：
要將這條線段轉化到平行座標空間去，我們需要畫出三個平行的垂直軸。在中間的軸上，我們選取三個點在 X 軸上的值，也就是 1，2，3 處畫一個點。在左邊的軸上，我們選取三個點在 Y 軸上的值，在 3，5，7 處畫一個點。在右邊的軸上我們做同樣的事情，但是取 Y 軸的負值。
接下來我們將代表 Y 軸值的點和它對應的 X 軸值連線起來。連線後的效果像下圖：
你會注意到在右邊的三條線會相交於一點。這個點的座標值代表了在真實空間中線段的斜率和截距。如果我們用一個向下斜的線段，那麼相交會發生在圖的左半邊。
如果我們取交點到中間軸的距離作為 u (在這個例子中是 2)，取豎直方向到 0 的距離作為 v (這裡是 1/3)，將軸之間的距離作為 d (這個例子中我使用的距離是 6)，我們可以用這樣的公式計算斜率和截距
斜率 = −1 + d/u
截距 = d × v/u
斜率是 2，截距是 1，和上面我們所畫的線段一致。
這種簡單有序的線段繪畫非常適合 GPU 進行，所以這種方法是一種利用 GPU 進行線段探測理想的方式。
探測線段的第一步是尋找可能代表一個線段的點。我們尋找的是位於邊界位置的點，並且我們希望將需要分析的點的數量控制在最少，所以之前談論的坎尼邊界探測是一個非常好的起點。
進行邊界探測之後，邊界點被用來在平行座標空間進行畫線。每一個邊界點會畫兩條線，一條在中間軸和左邊軸之間，另一條在中間軸和右邊軸之間。我們使用一種混合新增的方式使線段的交點變得更亮。在一片區域內最亮的點代表了線段。
舉例來說，我們可以從這張測試圖片開始：
下面是我們在平行座標空間中得到的 (我已經將負值對稱過來使圖片高度減半)
圖中的亮點就是我們探測到線段的地方。進行一個非最大值轉化來找到區域最值並將其他地方變為黑色。然後，點被轉化回線段的斜率和截距，得到下面的結果：
我必須指出在 GPUImage 中這個非最大值轉換過程是一個薄弱的環節。它可能會導致錯誤的探測出線段，或者在有噪點的地方將一條線段探測為多條線段。
正如之前所提到的，線段探測有許多有趣的應用。其中一種就是條形碼識別。有關平行座標空間轉換有趣的一點是，在真實空間中平行的線段轉換到平行座標空間中後是垂直對齊的一排點。不論平行線段是怎樣的都一樣。這就意味著你可以通過一排有特定順序間距的點來探測出條形碼無論條形碼是怎樣擺放的。這對於有視力障礙的手機使用者進行條形碼掃描是有巨大幫助的，畢竟他們無法看到盒子也很難將條形碼對齊。
對我而言，這種線段探測過程中的幾何學優雅是令我感到十分著迷的，我希望將它介紹給更多開發者。
小結
這些就是在過去幾年中發展出來的機器視覺方法中的幾個，它們僅僅是適合在 GPU 上工作的方法中的一部分。我個人認為在這個領域還有著令人激動的開創性工作要去做，這將會誕生可以提高許多人生活質量的應用。希望這篇文章至少為你提供了一個機器視覺領域簡要的總體介紹，並且展示了這個領域並不像許多開發者想象的那樣無法進入。
參考文獻
- I. Sobel. An Isotropic 3x3 Gradient Operator, Machine Vision for Three-Dimensional Scenes, Academic Press, 1990.
- J.M.S. Prewitt. Object Enhancement and Extraction, Picture processing and Psychopictorics, Academic Press, 1970.
- J. Canny. A Computational Approach To Edge Detection, IEEE Trans. Pattern Analysis and Machine Intelligence, 8(6):679–698, 1986.
- A. Ensor, S. Hall. GPU-based Image Analysis on Mobile Devices. Proceedings of Image and Vision Computing New Zealand 2011.
- C. Harris and M. Stephens. A Combined Corner and Edge Detector. Proc. Alvey Vision Conf., Univ. Manchester, pp. 147-151, 1988.
- J. Shi and C. Tomasi. Good features to track. Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 593-600, June 1994.
- A. Noble. Descriptions of Image Surfaces. PhD thesis, Department of Engineering Science, Oxford University 1989, p45.
- G. Ziegler, A. Tevs, C. Theobalt, H.-P. Seidel. GPU Point List Generation through HistogramPyramids. Research Report, Max-Planck-Institut fur Informatik, 2006.
- E. Rosten and T. Drummond. Machine learning for high-speed corner detection. European Conference on Computer Vision 2006.
- M. Dubská, J. Havel, and A. Herout. Real-Time Detection of Lines using Parallel Coordinates and OpenGL. Proceedings of SCCG 2011, Bratislava, SK, p. 7.
- M. Dubská, J. Havel, and A. Herout. PClines — Line detection using parallel coordinates. 2011 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), p. 1489- 1494.