5.3.6 直方图统计

1、灰度图像直方图统计

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

#include"vtkSmartPointer.h"
#include"vtkJPEGReader.h"
#include"vtkImageAccumulate.h"
#include"vtkImageData.h"
#include"vtkIntArray.h"
#include"vtkDataObject.h"
#include"vtkFieldData.h"
#include"vtkBarChartActor.h"
#include"vtkProperty2D.h"
#include"vtkTextProperty.h"
#include"vtkLegendBoxActor.h"
#include"vtkImageActor.h"
#include"vtkRenderer.h"
#include"vtkRenderWindow.h"
#include"vtkRenderWindowInteractor.h"

int main()
{
	vtkSmartPointer<vtkJPEGReader>reader = vtkSmartPointer<vtkJPEGReader>::New();//读入一幅灰度图像
	reader->SetFileName("data\\lena-gray.jpg");
	reader->Update();

	int bins = 16; //直方图一维数组的维数
	int comps = 1;
	vtkSmartPointer<vtkImageAccumulate>histogram = vtkSmartPointer<vtkImageAccumulate>::New();
	histogram->SetInputData(reader->GetOutput());
	histogram->SetComponentExtent(0, bins - 1, 0, 0, 0, 0);  //6个参数是RGB3个组分直方图的最大最小值
	histogram->SetComponentOrigin(0, 0, 0);          //函数设置的是统计每个组分直方图时的起始灰度值
	histogram->SetComponentSpacing(256.0 / bins, 0, 0);  //设置直方图每个间隔代表的灰度范围
	histogram->Update();
	int* output = static_cast<int*>(histogram->GetOutput()->GetScalarPointer());
	vtkSmartPointer<vtkIntArray> frequencies = vtkSmartPointer<vtkIntArray>::New();
	frequencies->SetNumberOfComponents(1);
	for (int j = 0; j < bins; ++j)
	{
		for (int i = 0; i < comps; ++i)
		{
		frequencies->InsertNextTuple1(*output++);
		}
	}
	vtkSmartPointer<vtkDataObject> dataObject = vtkSmartPointer<vtkDataObject>::New();
	dataObject->GetFieldData()->AddArray(frequencies);
	/*******************************************************/
	vtkSmartPointer<vtkBarChartActor> barChart = vtkSmartPointer<vtkBarChartActor>::New();
	barChart->SetInput(dataObject);
	barChart->SetTitle("Histogram");
	barChart->GetPositionCoordinate()->SetValue(0.05, 0.05, 0.0);
	barChart->GetPosition2Coordinate()->SetValue(0.95, 0.95, 0.0);
	barChart->GetProperty()->SetColor(0, 0, 0);
	barChart->GetTitleTextProperty()->SetColor(0, 0, 0);
	barChart->GetLabelTextProperty()->SetColor(0, 0, 0);
	barChart->GetLegendActor()->SetNumberOfEntries(dataObject->GetFieldData()->GetArray(0)->GetNumberOfTuples());
	//barChart->LegendVisibilityOff();
	//barChart->LabelVisibilityOff();

	double colors[3][3] = { { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } };
	int count = 0;
	for (int i = 0; i < bins; ++i)
	{
		for (int j = 0; j < comps; ++j)
		{
			barChart->SetBarColor(count++, colors[j]);//单通道红色,colors数组第一个元素的值
		}
	}
	vtkSmartPointer<vtkImageActor> imgActor = vtkSmartPointer<vtkImageActor>::New();
	imgActor->SetInputData(reader->GetOutput());
	/*******************************************************/
	double imgView[4] = { 0.0, 0.0, 0.5, 1.0 };
	double barView[4] = { 0.5, 0.0, 1.0, 1.0 };
	vtkSmartPointer<vtkRenderer> barRender = vtkSmartPointer<vtkRenderer>::New();
	barRender->SetViewport(barView);
	barRender->AddActor(barChart);
	barRender->SetBackground(1, 1, 1);

	vtkSmartPointer<vtkRenderer> imgRender = vtkSmartPointer<vtkRenderer>::New();
	imgRender->SetViewport(imgView);
	imgRender->AddActor(imgActor);
	imgRender->SetBackground(1, 1, 1);

	vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
	renderWindow->AddRenderer(barRender);
	renderWindow->AddRenderer(imgRender);
	renderWindow->SetSize(640, 320);
	renderWindow->Render();
	renderWindow->SetWindowName("Gray-Image Histogram");

	vtkSmartPointer<vtkRenderWindowInteractor> rwi = vtkSmartPointer<vtkRenderWindowInteractor>::New();
	rwi->SetRenderWindow(renderWindow);
	rwi->Initialize();
	rwi->Start();

	return EXIT_SUCCESS;
}

运行结果：

代码分析：

首先是读入一副灰度图像，一般的灰度图像的灰度范围为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对象显示图像直方图。

2、彩色图像直方图统计

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

#include <vtkSmartPointer.h>
#include <vtkBMPReader.h>
#include <vtkImageData.h>
#include <vtkXYPlotActor.h>  //vtkXYPlotActor类可以用来显示二维曲线，它可以接收多个输入数据
#include <vtkAxisActor2D.h>
#include <vtkImageAccumulate.h>
#include <vtkImageExtractComponents.h>
#include <vtkProperty2D.h>
#include <vtkTextProperty.h>
#include <vtkActor.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>

//测试图像：../data/lena.bmp
int main()
{
	
	vtkSmartPointer<vtkBMPReader> reader = vtkSmartPointer<vtkBMPReader>::New();
	reader->SetFileName("data//lena.bmp");
	reader->Update();

	int numComponents = reader->GetOutput()->GetNumberOfScalarComponents();
//**********************************************************************/
	vtkSmartPointer<vtkXYPlotActor> plot = vtkSmartPointer<vtkXYPlotActor>::New();//使用vtkBarChartActor柱状图来显示直方图，本例使用vtkXYPlotActor曲线来表示直方图。
	plot->ExchangeAxesOff();
	plot->SetLabelFormat("%g");
	plot->SetXTitle("Intensity");
	plot->SetYTitle("Frequency");
	plot->SetXValuesToValue();
	plot->GetProperty()->SetColor(0.0, 0.0, 0.0);
	plot->GetAxisLabelTextProperty()->SetColor(0.0, 0.0, 0.0);
	plot->GetAxisTitleTextProperty()->SetColor(0.0, 0.0, 0.0);

	double colors[3][3] = {{ 1, 0, 0 },{ 0, 1, 0 },{ 0, 0, 1 }};

	const char* labels[3] = { "Red", "Green", "Blue" };
	int xmax = 0;
	int ymax = 0;

	for (int i = 0; i < numComponents; ++i)
	{
		vtkSmartPointer<vtkImageExtractComponents> extract = vtkSmartPointer<vtkImageExtractComponents>::New();
		extract->SetInputConnection(reader->GetOutputPort()); //彩色图像不能直接计算直方图，因此需要先通过vtkImageExtractComponents来提取每个通道图像
		extract->SetComponents(i);
		extract->Update();

		double range[2];
		extract->GetOutput()->GetScalarRange(range);
		int extent = static_cast<int>(range[1]) - static_cast<int>(range[0]) - 1;

		vtkSmartPointer<vtkImageAccumulate> histogram =	vtkSmartPointer<vtkImageAccumulate>::New();
		histogram->SetInputConnection(extract->GetOutputPort());
		histogram->SetComponentExtent(0, extent, 0, 0, 0, 0);
		histogram->SetComponentOrigin(range[0], 0, 0);  //灰度起点为图像的最小灰度值
		histogram->SetComponentSpacing(1, 0, 0); //直方图的间隔取(1, 0, 0),每个灰度计算统计一个频率
		histogram->SetIgnoreZero(1);  //统计直方图时，像素值为0的像素不进行统计
		histogram->Update();
**********************************************************************/
		if (range[1] > xmax)
		{
			xmax = range[1];
		}
		if (histogram->GetOutput()->GetScalarRange()[1] > ymax)
		{
			ymax = histogram->GetOutput()->GetScalarRange()[1];
		}

		plot->AddDataSetInput(histogram->GetOutput());
		plot->SetPlotColor(i, colors[i]);
		plot->SetPlotLabel(i, labels[i]);
		plot->LegendOn(); //red\green\blue注释
	}

	plot->SetXRange(0, xmax);
	plot->SetYRange(0, ymax);

	vtkSmartPointer<vtkRenderer> renderer =	vtkSmartPointer<vtkRenderer>::New();
	renderer->AddActor(plot);
	renderer->SetBackground(1.0, 1.0, 1.0);

	vtkSmartPointer<vtkRenderWindow> renderWindow =	vtkSmartPointer<vtkRenderWindow>::New();
	renderWindow->AddRenderer(renderer);
	renderWindow->SetSize(640, 480);
	renderWindow->Render();
	renderWindow->SetWindowName("ImageAccumulateExample2");

	vtkSmartPointer<vtkRenderWindowInteractor> interactor =	vtkSmartPointer<vtkRenderWindowInteractor>::New();
	interactor->SetRenderWindow(renderWindow);
	interactor->Initialize();
	interactor->Start();

	return EXIT_SUCCESS;
} 

运行结果：

代码分析：

由于彩色图像不能直接计算直方图，因此需要先通过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对象建立可视化管道来显示图像直方图曲线。

注：此文知识学习笔记，仅记录完整程序和实现结果，具体原理参见：

https://blog.csdn.net/www_doling_net/article/details/8541534

https://blog.csdn.net/shenziheng1/article/category/6114053/4

参考资料：

1.《The Visualization Toolkit – AnObject-Oriented Approach To 3D Graphics (4th Edition)》
2. 张晓东, 罗火灵. VTK图形图像开发进阶[M]. 机械工业出版社, 2015.

所用软件：vtk7.0+visual studio 2013