5.12 边缘检测
5.12 邊緣檢測
圖像中不連續的灰度值會產生邊緣,圖像的邊緣檢測是基于邊界的圖像分割方法的基礎,例如分水嶺算法,通常是分割原圖的梯度圖像,梯度實際上也是反應的圖像邊緣信息。圖像邊緣一般常用圖像一階導數和二階導數來檢測。
5.12.1 梯度算子
梯度算子對應于圖像一階導數。圖像一階導數計算一般是通過差分運算來近似的。VTK中可以使用vtkImageGradient計算圖像梯度。注意圖像梯度是一個向量,具有方向和大小。因此vtkImageGradient的計算結果是一個梯度場,也就是每個像素值都是一個梯度向量。顯示梯度圖像時需要計算每個像素點的梯度大小,即模值。下面代碼演示了怎么計算圖像梯度。
?
1:???? vtkSmartPointer<vtkJPEGReader>reader =
?? 2:?????????vtkSmartPointer<vtkJPEGReader>::New();
?? 3:?????reader->SetFileName( "lena2.jpg" );
?? 4:?????reader->Update();
?? 5:??
?? 6:?????vtkSmartPointer<vtkImageGradient> gradientFilter =
?? 7:?????????vtkSmartPointer<vtkImageGradient>::New();
?? 8:?????gradientFilter->SetInputConnection(reader->GetOutputPort());
?? 9:?????gradientFilter->SetDimensionality(2);
? 10:??
? 11:?????vtkSmartPointer<vtkImageMagnitude> magnitudeFilter =
? 12:?????????vtkSmartPointer<vtkImageMagnitude>::New();
? 13:?????magnitudeFilter->SetInputConnection(gradientFilter->GetOutputPort());
? 14:?????magnitudeFilter->Update();
? 15:??
? 16:?????double range[2];
? 17:?????magnitudeFilter->GetOutput()->GetScalarRange(range);
? 18:??
? 19:?????vtkSmartPointer<vtkImageShiftScale> ShiftScale =
? 20:?????????vtkSmartPointer<vtkImageShiftScale>::New();
? 21:?????ShiftScale->SetOutputScalarTypeToUnsignedChar();
? 22:?????ShiftScale->SetScale( 255 / range[1] );
? 23:?????ShiftScale->SetInputConnection(magnitudeFilter->GetOutputPort());
? ????24:????? ShiftScale->Update();
vtkImageGradient的使用比較簡單,只需要設置輸入圖像即可。計算梯度時,采用的是中間差分法,即像素在每個方向的差分,都是利用的前后兩個像素值之差。這樣在圖像在邊界處的差分計算需要特殊處理。其內部定義了HandleBoundaries變量,通過函數SetHandleBoundaries()定賦值。當HandleBoundaries為真時算子會特殊處理計算邊界像素的梯度;當為假時不計算邊界像素的梯度值,因此輸出圖像大小要小于輸入圖像。另外函數SetDimensionality()用于設置要計算的圖像維數,默認為二維,此時梯度向量也為二維。
前面也提到過,梯度是一個向量,不能直接顯示。因此上面代碼中定義了vtkImageMagnitude對象來計算梯度向量的2范數,即向量的模。利用vtkImageShiftScale將圖像的數據范圍調整到0-255然后顯示。另外還可以通過vtkImageExtractComponents來提取每個方向的梯度分量進行顯示。注意,彩色圖像不能直接用來計算梯度,需要先轉換為灰度圖像。本例的執行結果如下圖所示。
?
圖5.26 梯度算子
索貝爾算子(Sobel)也是一種常用的梯度算子(圖5.27)。Sobel算子計算稍微復雜,它采用3x3的模板。計算時模板在圖像上移動,并在每個位置上計算對應中心像素的梯度值。VTK中vtkSobel2D計算圖像的sobel算子,使用代碼如下:
?
圖5.27 Sobel算子
?
1:????? vtkSmartPointer<vtkJPEGReader>reader =
?? 2:?????????vtkSmartPointer<vtkJPEGReader>::New();
?? 3:?????reader->SetFileName( "lena2.jpg" );
?? 4:?????reader->Update();
?? 5:??
?? 6:?????vtkSmartPointer<vtkImageSobel2D> sobelFilter =
?? 7:?????????vtkSmartPointer<vtkImageSobel2D>::New();
?? 8:?????sobelFilter->SetInputConnection(reader->GetOutputPort());
?? 9:??
? 10:?????vtkSmartPointer<vtkImageExtractComponents> extractXFilter =
? 11:?????????vtkSmartPointer<vtkImageExtractComponents>::New();
? 12:?????extractXFilter->SetComponents(0);
? 13:?????extractXFilter->SetInputConnection(sobelFilter->GetOutputPort());
? 14:?????extractXFilter->Update();
? 15:??
? 16:?????double xRange[2];
? 17:?????extractXFilter->GetOutput()->GetScalarRange(xRange);
? 18:??
? 19:?????vtkSmartPointer<vtkImageMathematics> xImageAbs =
? 20:?????????vtkSmartPointer<vtkImageMathematics>::New();
? 21:?????xImageAbs->SetOperationToAbsoluteValue();
? 22:?????xImageAbs->SetInputConnection(extractXFilter->GetOutputPort());
? 23:?????xImageAbs->Update();
? 24:??
? 25:?????vtkSmartPointer<vtkImageShiftScale> xShiftScale =
? 26:?????????vtkSmartPointer<vtkImageShiftScale>::New();
? 27:?????xShiftScale->SetOutputScalarTypeToUnsignedChar();
? 28:?????xShiftScale->SetScale( 255 / xRange[1] );
? 29:?????xShiftScale->SetInputConnection(xImageAbs->GetOutputPort());
? 30:?????xShiftScale->Update();
? 31:??
? 32:?????vtkSmartPointer<vtkImageExtractComponents> extractYFilter =
? 33:?????????vtkSmartPointer<vtkImageExtractComponents>::New();
? 34:?????extractYFilter->SetComponents(1);
? 35:? ????extractYFilter->SetInputConnection(sobelFilter->GetOutputPort());
? 36:?????extractYFilter->Update();
? 37:??
? 38:?????double yRange[2];
? 39:?????extractYFilter->GetOutput()->GetScalarRange(yRange);
? 40:??
? 41:?????vtkSmartPointer<vtkImageMathematics> yImageAbs =
? 42:?????????vtkSmartPointer<vtkImageMathematics>::New();
? 43:?????yImageAbs->SetOperationToAbsoluteValue();
? 44:?????yImageAbs->SetInputConnection(extractYFilter->GetOutputPort());
? 45:?????yImageAbs->Update();
? 46:??
? 47:?????vtkSmartPointer<vtkImageShiftScale> yShiftScale =
? 48:?????????vtkSmartPointer<vtkImageShiftScale>::New();
? 49:?????yShiftScale->SetOutputScalarTypeToUnsignedChar();
? 50:?????yShiftScale->SetScale( 255 / yRange[1] );
? 51:?????yShiftScale->SetInputConnection(yImageAbs->GetOutputPort());
? 52:?????yShiftScale->Update();
?
該例中計算利用Sobel算子計算圖像的梯度圖像,然后提取X方向的梯度分量和Y方向的梯度分量。由于計算Sobel算子的值可能存在負值,因此利用vtkImageMathematics對各個分量圖像計算絕對值,再由vtkImageShiftScale將圖像的數值范圍調節到0-255之間再顯示。執行結果如下。
?
圖5.28 vtkSobel2D計算圖像的sobel算子
5.12.2 Canny算子
Canny算子是John Canny于20世紀80年代提出的一種多級邊緣檢測算法。John Canny研究了最優邊緣的特性,即檢測到的邊緣要盡可能跟實際的邊緣接近并盡可能的多,同時,要盡量降低噪聲對邊緣檢測的干擾。其計算步驟如下
1)對源圖像進行高斯平滑以消除圖像中噪聲
2)采用差分法近似計算圖像每一個像素的梯度,并計算梯度的模值和方向
3)對梯度進行"非極大抑制":圖像邊緣點梯度值通常在梯度方向是極大值,因此檢測邊緣需要將非極大值賦值0來抑制非邊緣點。檢測方法就是在一個局部窗口內,如果中心像素點的梯度不比梯度方向上相鄰兩個像素值大,那么該中心像素點梯度值賦0。
4)雙閾值法檢測邊緣和連接邊緣。取兩個梯度閾值high和low,將梯度圖像中小于high的像素賦0得到邊緣圖像I1,該圖像能夠接近圖像邊緣但是可能會存在間斷點;將梯度圖像中小于low的像素賦0得到邊緣圖像I2,該圖中受噪聲影響比較大,但是邊緣信息更多。在連接邊緣時,以I1為基礎,對非零點進行邊緣跟蹤,如果追蹤過程中出現中斷,則從I2對應像素點及其鄰域來尋找可以連接的邊緣,直至結束。
以上是Canny算子的計算步驟。在VTK中沒有實現一個專門的類來做Canny邊緣檢測。但是我們可以根據以上步驟來實現。
?
?? 1: ?????vtkSmartPointer<vtkJPEGReader> reader =
?? 2:?????????vtkSmartPointer<vtkJPEGReader>::New();
?? 3:?????reader->SetFileName( "lena2.jpg" );
?? 4:?????reader->Update();
?? 5:??
?? 6:?????vtkSmartPointer<vtkImageCast> ic =
?? 7:?????????vtkSmartPointer<vtkImageCast>::New();
?? 8:?????ic->SetOutputScalarTypeToFloat();
?? 9:?????ic->SetInputConnection(reader->GetOutputPort());
? 10:??
? 11:?????vtkSmartPointer<vtkImageGaussianSmooth> gs =
? 12:?????????vtkSmartPointer<vtkImageGaussianSmooth>::New();
? 13:?????gs->SetInputConnection(ic->GetOutputPort());
? 14:?????gs->SetDimensionality(2);
? 15:?????gs->SetRadiusFactors(1, 1, 0);
? 16:??
? 17:?????vtkSmartPointer<vtkImageGradient> imgGradient =
? 18:?????????vtkSmartPointer<vtkImageGradient>::New();
? 19:?????imgGradient->SetInputConnection(gs->GetOutputPort());
? 20:?????imgGradient->SetDimensionality(2);
? 21:??
? 22:?????vtkSmartPointer<vtkImageMagnitude> imgMagnitude =
? 23:?????????vtkSmartPointer<vtkImageMagnitude>::New();
? 24:?????imgMagnitude->SetInputConnection(imgGradient->GetOutputPort());
? 25:??
? 26:?????vtkSmartPointer<vtkImageNonMaximumSuppression> nonMax =
? 27:?????????vtkSmartPointer<vtkImageNonMaximumSuppression>::New();
? 28:?????nonMax->SetMagnitudeInput(imgMagnitude->GetOutput());
? 29:?????nonMax->SetVectorInput(imgGradient->GetOutput());
? 30:?????nonMax->SetDimensionality(2);
? 31:??
? 32:?????vtkSmartPointer<vtkImageConstantPad> pad =
? 33:?????????vtkSmartPointer<vtkImageConstantPad>::New();
? 34:?????pad->SetInputConnection(imgGradient->GetOutputPort());
? 35:?????pad->SetOutputNumberOfScalarComponents(3);
? 36:?????pad->SetConstant(0);
? 37:??
? 38:?????vtkSmartPointer<vtkImageToStructuredPoints> i2sp1 =
? 39:?????????vtkSmartPointer<vtkImageToStructuredPoints>::New();
? 40:?????i2sp1->SetInputConnection(nonMax->GetOutputPort());
? 41:?????i2sp1->SetVectorInput(pad->GetOutput());
? 42:??
? 43:?????vtkSmartPointer<vtkLinkEdgels> imgLink =
? 44:?????????vtkSmartPointer<vtkLinkEdgels>::New();
? 45:?????imgLink->SetInput(i2sp1->GetOutput());
? 46:?????imgLink->SetGradientThreshold(2);
? 47:??
? 48:?????vtkSmartPointer<vtkThreshold> thresholdEdgels =
? 49:?????????vtkSmartPointer<vtkThreshold>::New();
? 50:?????thresholdEdgels->SetInputConnection(imgLink->GetOutputPort());
? 51:?????thresholdEdgels->ThresholdByUpper(10);
? 52:?????thresholdEdgels->AllScalarsOff();
? 53:??
? 54:?????vtkSmartPointer<vtkGeometryFilter> gf =
? 55:?????????vtkSmartPointer<vtkGeometryFilter>::New();
? 56:?????gf->SetInputConnection(thresholdEdgels->GetOutputPort());
? 57:??
? 58:?????vtkSmartPointer<vtkImageToStructuredPoints> i2sp =
? 59:?????????vtkSmartPointer<vtkImageToStructuredPoints>::New();
? 60:?????i2sp->SetInputConnection(imgMagnitude->GetOutputPort());
? 61:?????i2sp->SetVectorInput(pad->GetOutput());
? 62:??
? 63:?????vtkSmartPointer<vtkSubPixelPositionEdgels> spe =
? 64:?????????vtkSmartPointer<vtkSubPixelPositionEdgels>::New();
? 65:?????spe->SetInputConnection(gf->GetOutputPort());
? 66:?????spe->SetGradMaps(i2sp->GetStructuredPointsOutput());
? 67:??
? 68:?????vtkSmartPointer<vtkStripper> strip =
? 69:?????????vtkSmartPointer<vtkStripper>::New();
? 70:?????strip->SetInputConnection(spe->GetOutputPort());
? 71:?????strip->Update();
?
該程序比較復雜,處理邊緣時將其作為幾何數據來進行處理。因此涉及了部分幾何數據操作的filter,這里如果不明白可以先放一下,再幾何數據處理部分會做詳細介紹。程序首先讀入圖像,計算圖像的梯度和模值。接下來按照Canny算子的步驟進行處理。我們依次來介紹用到的相應的filter。
vtkImageNonMaximumSuppression將圖像中的非局部峰值設置為0,輸入和輸出類型都是vtkImageData:其中輸入有兩個,模值圖像(magnitude)和向量圖像,一個典型的應用就是輸入梯度模值圖像和梯度向量圖像對梯度做非極大值抑制。
vtkImageConstantPad增加圖像的大小,其輸入和輸出都為vtkImageData。其中函數SetOutputNumberOfScalarComponents(3)用于設置輸出圖像的像素數據組分個數,函數SetConstant(0)用于設置輸出圖像中擴大的區域像素值。而SetOutputWholeExtent()則用于設置輸出圖像的范圍。這里的作用是將梯度圖像像素的組分修改為3,方便下面vtkImageToStructuredPoints使用。
vtkImageToStructuredPoints將vtkImageData格式轉換為規則點集。該類的輸入類型是vtkImageData,另外還有一個可選的RGB三組分向量圖像輸入;輸出類型是vtkStructuredPoints,當輸入向量圖像時,向量圖像像素數據會轉為輸出圖像的對應點的屬性。
vtkLinkEdgels類根據點的相鄰關系連接成連續的折線Polyline。其內部閾值變量GradientThreshold,可以用來排除輸入點中梯度值小于該閾值的點。當使用vtkLinkEdgels進行Canny算子的雙閾值邊緣檢測時,GradientThreshold可以用作較小的閾值。設置該閾值的函數是SetGradientThreshold(2)。
vtkThreshold用于獲取輸入任意類型數據的滿足閾值條件的單元數據。該類的輸入為VTK的任意數據類型,輸出數據類型是不規則網格。閾值設置有:大于閾值,小于閾值和介于兩個閾值之間。內部提供了兩種屬性模式AttributeMode設置,即閾值比較時是采用的點屬性還是單元屬性,默認下是點屬性。而當屬性為多元數據時,還需要設置閾值比較時使用哪個組分的數據。其中提供了三種模式選擇,所有組分都滿足閾值條件,任意一個滿足閾值條件和用戶指定的組分滿足閾值條件。當使用點屬性數據時,如果設置了AllScalars,那么單元滿足閾值條件的前提會是其所有點的屬性都滿足閾值條件。這里將閾值設置為10,即Canny中雙閾值的較大閾值。
vtkGeometryFilter將數據轉換為幾何數據,輸出類型為vtkPolyData。該類從vtkThreshold的輸出中提取圖像邊緣的幾何數據。
vtkSubPixelPositionEdgels接收一系列連續曲線及其對應的梯度系信息作為輸入,利用梯度信息來調整曲線位置。這里對前面提取的圖像邊緣再根據其梯度進行調整。
vtkStripper用來將輸入的多邊形、三角形或者線段生成三角形帶或者折線段。輸入的多邊形數據必須是三角形,否則不會進行帶化處理。因此處理多邊形數據時,可以先用vtkTriangleFilter進行三角化后再使用本類。如果輸入中存在孤立點的話,也不會進行任何處理。默認情況下,該filter處理后會丟棄掉屬性數據。
5.12.3 拉普拉斯算子
拉普拉斯算子是一個二階邊緣算子,即梯度的散度。拉普拉斯算子的實現也是通過模板實現。常用的拉普拉斯模板定義如下:
?
圖5.29拉普拉斯算子
拉普拉斯算子計算圖像的二階導數,對于圖像噪聲比較敏感。拉普拉斯算子的結果為標量,表示邊緣的寬度。但是它常產生雙像素寬邊緣,而且不能提供方向信息,因此較少直接用于邊緣檢測。在VTK中由vtkImageLaplacian實現。
?
1:???? vtkSmartPointer<vtkJPEGReader>reader =
??2:?????????vtkSmartPointer<vtkJPEGReader>::New();
??3:????? reader->SetFileName("lena2.jpg" );
??4:????? reader->Update();
??5:? ?
??6:?????vtkSmartPointer<vtkImageLaplacian> lapFilter =
??7:?????????vtkSmartPointer<vtkImageLaplacian>::New();
??8:?????lapFilter->SetInputConnection(reader->GetOutputPort());
??9:???? ?lapFilter->SetDimensionality(2);
?10:? ?
?11:????? double range[2];
?12:?????lapFilter->GetOutput()->GetScalarRange(range);
?13:? ?
?14:?????vtkSmartPointer<vtkImageShiftScale> ShiftScale =
?15:?????????vtkSmartPointer<vtkImageShiftScale>::New();
?16:?????ShiftScale->SetOutputScalarTypeToUnsignedChar();
?17:????? ShiftScale->SetScale(255 / (range[1]-range[0]) );
?18:?????ShiftScale->SetShift(-range[0]);
?19:?????ShiftScale->SetInputConnection(lapFilter->GetOutputPort());
?20:????? ShiftScale->Update();
?
vtkImageLaplacian輸入和輸出數據都是vtkImageData,與梯度算子不同,該filter的輸出圖像像素為標量。函數SetDimensionality用于設置輸入圖像的維數,默認為2維。計算完畢后,利用vtkImageShiftScale將圖像的數據范圍變換至0-255之間。計算結果如下圖所示:
?
圖5.30 vtkImageLaplacian計算結果
由于拉普拉斯算子對噪聲比較敏感,為了減少噪聲應用,可以先對圖像做高斯濾波來平滑圖像,然后再計算拉普拉斯算子,稱為LoG算子(laplacian of gaussian)。關于高斯平滑我們將在圖像平滑一節中介紹。
==========歡迎轉載,轉載時請保留該聲明信息==========
版權歸@東靈工作室所有,更多信息請訪問東靈工作室
教程系列導航:http://blog.csdn.net/www_doling_net/article/details/8763686
總結
- 上一篇: 5.13 图像平滑
- 下一篇: 统计测序数据reads数和碱基数的几种方