5.8 直方图统计

5.8.1 灰度图像直方图

直方图统计是图像处理中的一个非常重要的操作。VTK中实现直方图统计功能的filter是vtkImageAccumulate。其将每个组分的数值范围划分为离散的间隔，然后统计每个灰度间隔上的像素数目。vtkImageAccumulate输入和输出都是vtkImageData类型，因此直方图也可以看做是一幅图像；对于输入图像的像素数据类型可以是任意的，但是最大支持3个组分像素类型，而输出图像的像素数据类型为int型。一个灰度图像的直方图为一个一维图像。首先来看一下怎么计算灰度图像直方图。

1: vtkSmartPointer<vtkJPEGReader> reader =

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

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

4: reader->Update();

6: int bins = 16;

8: vtkSmartPointer<vtkImageAccumulate> histogram =

9: vtkSmartPointer<vtkImageAccumulate>::New();

10: histogram->SetInput(reader->GetOutput());

11: histogram->SetComponentExtent(0, bins-1, 0, 0, 0, 0);

12: histogram->SetComponentOrigin(0, 0, 0);

13: histogram->SetComponentSpacing(256/bins, 0, 0);

14: histogram->Update();

15:

16: vtkSmartPointer<vtkIntArray> frequencies =

17: vtkSmartPointer<vtkIntArray>::New();

18:

19: frequencies->SetNumberOfComponents(1);

20: frequencies->SetNumberOfTuples(bins);

21: int* output = static_cast<int*>(histogram->GetOutput()->GetScalarPointer());

22:

23: for(int j = 0; j < bins; ++j)

24: {

25: frequencies->SetTuple1(j, *output++);

26: }

27:

28: vtkSmartPointer<vtkDataObject> dataObject =

29: vtkSmartPointer<vtkDataObject>::New();

30: dataObject->GetFieldData()->AddArray( frequencies );

31:

32: vtkSmartPointer<vtkBarChartActor> barChart =

33: vtkSmartPointer<vtkBarChartActor>::New();

34:

35: barChart->SetInput(dataObject);

36: barChart->SetTitle("Histogram");

37: barChart->GetPositionCoordinate()->SetValue(0.05,0.05,0.0);

38: barChart->GetPosition2Coordinate()->SetValue(0.95,0.95,0.0);

39: barChart->GetProperty()->SetColor(1,1,1);

40:

41: barChart->GetLegendActor()->SetNumberOfEntries(

42: dataObject->GetFieldData()->GetArray(0)->GetNumberOfTuples());

43: barChart->LegendVisibilityOff();

44: barChart->LabelVisibilityOff();

45:

46: double red[3] = {1, 0, 0 };

47: for(int i = 0; i < bins; ++i)

48: {

49: barChart->SetBarColor(i, red );

50: }

51:

52: vtkSmartPointer<vtkRenderer> renderer =

53: vtkSmartPointer<vtkRenderer>::New();

54: renderer->AddActor(barChart);

55:

56: vtkSmartPointer<vtkRenderWindow> renderWindow =

57: vtkSmartPointer<vtkRenderWindow>::New();

58: renderWindow->AddRenderer(renderer);

59: renderWindow->SetSize(640, 480);

60:

61: vtkSmartPointer<vtkRenderWindowInteractor> interactor =

62: vtkSmartPointer<vtkRenderWindowInteractor>::New();

63: interactor->SetRenderWindow(renderWindow);

64: interactor->Initialize();

65: interactor->Start();

下面来分析一下代码。首先是读入一副灰度图像，一般的灰度图像的灰度范围为0-255。定义了一个变量bins = 16，表示要图像灰度范围上的间隔数目，也可以理解为直方图一维数组的维数。然后定义vtkImageAccumulate对象，并设置输入数据为我们读入的图像数据，接着调用了三个函数：

SetComponentExtent(0,bins-1, 0, 0, 0, 0)，该函数设置要计算每个组分的直方图的最小和最大值。vtkImageAccumulate最大支持像素值为三个组分（如RGB图像）的直方图，支持共有六个参数。分别表示每个组分的直方图最小和最大值。该例中由于计算的是灰度图像直方图，只有一个组分，因此第二个和第三个组分都设置为0；而第一组分直方图维数为bins = 16，那么其最小和最大范围为0和bins-1。

SetComponentOrigin(0,0, 0)，该函数设置的是统计每个组分直方图时的起始灰度值，这里设置为0，表示灰度从0开始统计直方图。同样，vtkImageAccumulate最大支持像素值为三个组分，这里也要设置三个参数。如果图像的灰度范围为[1000, 2000]，那么计算直方图时，其起始灰度应该设置为1000。

SetComponentSpacing(16,0, 0)，设置直方图每个间隔代表的灰度范围，例如当一个图像灰度范围为[1000, 2000]，统计直方图的间隔数bins为100时，那么对应的space应该设置为SetComponentSpacing(100, 0, 0)。

参数设置完毕后执行Update()即可计算直方图。前面已经提到过，vtkImageAccumulate的输出结果也是一个vtkImageData类型，这样就可以方便的访问图像的每个数据。图像像素访问前面已经介绍过，这里就不在重述。需要注意的是输出直方图图像的数据类型为int。

这里再顺便简单的介绍一下直方图的显示。虽然vtkImageAccumulate的输出类型为vtkImageData但是并不能直接按照图像的方式进行显示。VTK中定义了vtkBarChartActor用来显示条形图，因此可以利用其来显示直方图。但是该类接收的数据类型为vtkDataObject类型，因此需要先将直方图数据进行转换。首先将直方图数组存储到vtkIntArray数组frequencies中，通过函数vtkDataObject函数GetFieldData()->AddArray(frequencies)将其添加到vtkDataObject对象中。vtkBarChartActor对象接收vtkDataObject对象作为输入，另外还需要设置图表的名字，颜色等，需要注意两个函数：

barChart->GetPositionCoordinate()->SetValue(0.05,0.05,0.0);

barChart->GetPosition2Coordinate()->SetValue(0.95,0.95,0.0);

这里设置的是窗口中显示图表的所在矩形的左下角点和右上角点坐标，VTK的坐标系原点位于左下角点，设置时需要格外注意。设置完毕后，即可定义相应的vtkRenderer，vtkRenderWindow和vtkRenderWindowInteractor对象显示图像直方图。本例的显示效果如下：

图5.20 VTK直方图显示

5.8.2 彩色图像直方图

彩色图像由于内部有三个通道，不能直接计算直方图，需要提取RGB三个通道数据，分别计算直方图。每个通道计算直方图的方法与灰度图像直方图计算方法一致。下面给出相关代码：

1: vtkSmartPointer<vtkBMPReader>reader =

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

3: reader->SetFileName ( "lena.bmp" );

4: reader->Update();

6: int numComponents =reader->GetOutput()->GetNumberOfScalarComponents();

8: vtkSmartPointer<vtkXYPlotActor> plot =

9: vtkSmartPointer<vtkXYPlotActor>::New();

10: plot->ExchangeAxesOff();

11: plot->SetLabelFormat( "%g" );

12: plot->SetXTitle( "Value" );

13: plot->SetYTitle( "Frequency" );

14: plot->SetXValuesToValue();

15:

16: double colors[3][3] = {

17: { 1, 0, 0 },

18: { 0, 1, 0 },

19: { 0, 0, 1 }

20: };

21:

22: const char* labels[3] = { "Red", "Green","Blue" };

23:

24: int xmax = 0;

25: int ymax = 0;

26:

27: for( int i = 0; i < numComponents; ++i )

28: {

29: vtkSmartPointer<vtkImageExtractComponents> extract =

30: vtkSmartPointer<vtkImageExtractComponents>::New();

31: extract->SetInputConnection( reader->GetOutputPort() );

32: extract->SetComponents( i );

33: extract->Update();

34:

35: double range[2];

36: extract->GetOutput()->GetScalarRange( range );

37: int extent =static_cast<int>(range[1])-static_cast<int>(range[0])-1;

38:

39: vtkSmartPointer<vtkImageAccumulate> histogram =

40: vtkSmartPointer<vtkImageAccumulate>::New();

41: histogram->SetInputConnection( extract->GetOutputPort() );

42: histogram->SetComponentExtent( 0,extent, 0,0, 0,0);

43: histogram->SetComponentOrigin( range[0],0,0 );

44: histogram->SetComponentSpacing( 1,0,0 );

45: histogram->SetIgnoreZero( 1 );

46: histogram->Update();

47:

48: if( range[1] > xmax )

49: {

50: xmax = range[1];

51: }

52: if( histogram->GetOutput()->GetScalarRange()[1] > ymax )

53: {

54: ymax =histogram->GetOutput()->GetScalarRange()[1];

55: }

56:

57: plot->AddInput( histogram->GetOutput() );

58: plot->SetPlotColor(i,colors[i]);

59: plot->SetPlotLabel(i,labels[i]);

60: plot->LegendOn();

61: }

62:

63: plot->SetXRange( 0, xmax);

64: plot->SetYRange( 0, ymax);

65:

66: vtkSmartPointer<vtkRenderer> renderer =

67: vtkSmartPointer<vtkRenderer>::New();

68: renderer->AddActor(plot);

69:

70: vtkSmartPointer<vtkRenderWindow> renderWindow =

71: vtkSmartPointer<vtkRenderWindow>::New();

72: renderWindow->AddRenderer( renderer );

73: renderWindow->SetSize(640, 480);

74:

75: vtkSmartPointer<vtkRenderWindowInteractor> interactor =

76: vtkSmartPointer<vtkRenderWindowInteractor>::New();

77: interactor->SetRenderWindow( renderWindow );

78:

79: interactor->Initialize();

80: interactor->Start();

上面代码说明了怎样计算彩色图像直方图。计算直方图的主要代码段是27-61行。由于彩色图像不能直接计算直方图，因此需要先通过vtkImageExtractComponents来提取每个通道图像，然后再利用vtkImageAccumulate统计直方图。在本例中计算直方图的间隔取(1, 0, 0)，即每个灰度计算统计一个频率，而且灰度起点为图像的最小灰度值，这样间隔的个数即为：最大灰度值减去最小灰度值，再减1，如第37行代码。同时，设置了SetIgnoreZero()为1，即在统计直方图时，像素值为0的像素不进行统计。

在灰度图像直方图实例中，我们使用的是vtkBarChartActor柱状图来显示直方图，在本例中则使用vtkXYPlotActor曲线来表示直方图，这里做一简单介绍。vtkXYPlotActor类可以用来显示二维曲线，它可以接收多个输入数据，如本例中我们输入了三条曲线，分别是图像红色分量直方图区域，绿色分量直方图曲线和蓝色分量直方图曲线。SetXRange()和SetYRange()用来设置X轴和Y轴的数据范围，另外还可以设置X轴和Y轴的名字，曲线的标题等属性，详细可以查阅vtkXYPlotActor类的文档。vtkXYPlotActor类是一个vtkActor2D的子类，因此定义相应的vtkRenderer，vtkRenderWindow和vtkRenderWindowInteractor对象建立可视化管道来显示图像直方图曲线。本例的显示效果如下：其中，红色曲线代码红色分量的直方图，绿色代表绿色分量的直方图曲线，蓝色代码蓝色分量的直方图曲线。

图5.21 彩色图像直方图

5.9 图像重采样

图像重采样是指对采样后形成的由离散数据组成的数字图像按所需的像元位置或像元问距重新采样,以构成几何变换后的新图像。重采样过程本质上是图像恢复过程,它用输入的离散数字图像重建代表原始图像二维连续函数,再按新的像元间距和像元位置进行采样。其数学过程是根据重建的连续函数(曲面),用周围若干像元点的值估计或内插出新采样点的值。图像重采样在图像处理中应用非常广泛，如SIFT特征提取。

图像重采样后图像的维数会发生改变。当重采样图像小于原图像维数时，称为降采样；当重采样图像维数大于原图像时，称为升采样。VTK中可以方便的对图像进行重采样。vtkImageShrink3D类实现图像降采样。降采样需要设置每个方向的采样率，降采样率越大，图像越模糊。升采样的原理与降采样原理一致，只是增加采样点数来增加图像的维数。VTK中vtkImageMagnify来实现图像的升采样。下面看一个实例。

1: vtkSmartPointer<vtkBMPReader>reader =

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

3: reader->SetFileName (" lena.bmp" );

4: reader->Update();

6: vtkSmartPointer<vtkImageShrink3D> shrinkFilter =

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

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

9: shrinkFilter->SetShrinkFactors(20,20,1);

10: shrinkFilter->Update();

11:

12: vtkSmartPointer<vtkImageMagnify> magnifyFilter =

13: vtkSmartPointer<vtkImageMagnify>::New();

14: magnifyFilter->SetInputConnection(reader->GetOutputPort());

15: magnifyFilter->SetMagnificationFactors(10,10,1);

16: magnifyFilter->Update();

17:

18: int originalDims[3];

19: reader->GetOutput()->GetDimensions(originalDims);

20:

21: double originalSpace[3];

22: reader->GetOutput()->GetSpacing(originalSpace);

23:

24: int shrinkDims[3];

25: shrinkFilter->GetOutput()->GetDimensions(shrinkDims);

26:

27: double shrinkSpace[3];

28: shrinkFilter->GetOutput()->GetSpacing(shrinkSpace);

29:

30: int magnifyDims[3];

31: magnifyFilter->GetOutput()->GetDimensions(magnifyDims);

32:

33: double magnifySpace[3];

34: magnifyFilter->GetOutput()->GetSpacing(magnifySpace);

35:

36: std::cout<<"原图图像维数："<<originalDims[0]<< " "<<originalDims[1]<<""<<originalDims[2]<<std::endl;

37: std::cout<<"原图图像像素间隔："<<originalSpace[0]<< " "<<originalSpace[1]<<""<<originalSpace[2]<<std::endl;

38: std::cout<<"降采样图像维数："<<shrinkDims[0]<< " "<<shrinkDims[1]<<""<<shrinkDims[2]<<std::endl;

39: std::cout<<"降采样图像像素间隔："<<shrinkSpace[0]<< " "<<shrinkSpace[1]<<""<<shrinkSpace[2]<<std::endl;

40: std::cout<<"升采样图像维数："<<magnifyDims[0]<< " "<<magnifyDims[1]<<""<<magnifyDims[2]<<std::endl;

41: std::cout<<"升采样图像像素间隔："<<magnifySpace[0]<< " "<<magnifySpace[1]<<""<<magnifySpace[2]<<std::endl;

vtkImageShrink3D使用时通过SetShrinkFactors()设置X、Y和Z方向的采样率，参数为三个int类型数据。vtkImageMagnify的使用也是类似，通过SetMagnificationFactors设置相应的放大采样率。上面例子中分别对原图像维数缩小原来的1/20和放大为原来的10倍，执行结果如下：

图5.22 图像重采样

从图中可以看成，重采样图像已经变得十分模糊了。而升采样图像则变化不大。程序中输出了原图和重采样图像的维数和像素间隔，从输出来看，原图像的维数为512×512，缩小20倍后，图像维数变为25×25。需要注意的是，原图的像素间隔为(1,1,1)，而重采样后像素间隔变为(20, 20, 1)，这是因为采样是在世界坐标系下进行的，当采样图像的维数减小时，采样的像素间隔也相应的变大，这样保持图像的区域不会发生改变。同样的道理，当升采样后图像的维数变为原来的10倍，而像素间隔变为原来的十分之一。

另外，图像的重采样类有：vtkImageResample。