5.12 边缘检测

图像中不连续的灰度值会产生边缘，图像的边缘检测是基于边界的图像分割方法的基础，例如分水岭算法，通常是分割原图的梯度图像，梯度实际上也是反应的图像边缘信息。图像边缘一般常用图像一阶导数和二阶导数来检测。

5.12.1 梯度算子

梯度算子对应于图像一阶导数。图像一阶导数计算一般是通过差分运算来近似的。VTK中可以使用vtkImageGradient计算图像梯度。注意图像梯度是一个向量，具有方向和大小。因此vtkImageGradient的计算结果是一个梯度场，也就是每个像素值都是一个梯度向量。显示梯度图像时需要计算每个像素点的梯度大小，即模值。下面代码演示了怎么计算图像梯度。

1: vtkSmartPointer<vtkJPEGReader>reader =

2: vtkSmartPointer<vtkJPEGReader>::New();

3: reader->SetFileName( "lena2.jpg" );

4: reader->Update();

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算子计算稍微复杂，它采用3×3的模板。计算时模板在图像上移动，并在每个位置上计算对应中心像素的梯度值。VTK中vtkSobel2D计算图像的sobel算子，使用代码如下：

图5.27 Sobel算子

1: vtkSmartPointer<vtkJPEGReader>reader =

2: vtkSmartPointer<vtkJPEGReader>::New();

3: reader->SetFileName( "lena2.jpg" );

4: reader->Update();

6: vtkSmartPointer<vtkImageSobel2D> sobelFilter =

7: vtkSmartPointer<vtkImageSobel2D>::New();

8: sobelFilter->SetInputConnection(reader->GetOutputPort());

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();

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();

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之间。计算结果如下图所示：