ITK/VTK：绘制二维及三维DICOM图像分布直方图（附源码）

        使用ITK及VTK读取DICOM图像序列，并绘制图像统计直方图。

        注意事项及实现思路：

        1. 必须使用ITK读取DCM序列，使用VTK读取数据会产生失真，详见使用VTK和Matlab读取DICOM图像数据失真问题；

        2. 使用ITK读取数据后，将统计结果存入一维数组中，使用VTK的vtkBarChartActor类绘制直方图；

        3. 使用vtkBarChartActor前需要将数组转换成vtkIntArray、vtkDataObject类型。

绘制结果如图所示：

代码：

#include "vtkDICOMImageReader.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkMarchingCubes.h"
#include "vtkStripper.h"
#include "vtkActor.h"
#include "vtkPolyDataMapper.h"
#include "vtkProperty.h"
#include "vtkCamera.h"
#include "vtkBoxWidget.h"
#include "vtkSmartPointer.h" 
#include "vtkTriangleFilter.h"
#include "vtkMassProperties.h"
#include "vtkSmoothPolyDataFilter.h"
#include "vtkPolyDataNormals.h"
#include "vtkContourFilter.h"
#include "vtkRecursiveDividingCubes.h"
#include "vtkSTLWriter.h"
#include "vtkStringArray.h"


#include <vtkDataObject.h> //
#include <vtkFieldData.h> 
#include <vtkImageAccumulate.h>
#include <vtkImageData.h>
#include <vtkIntArray.h>
#include <vtkBarChartActor.h>
#include <vtkProperty2D.h>
#include <vtkTextProperty.h>
#include <vtkLegendBoxActor.h>
#include <vtkImageActor.h>

#include "itkGDCMImageIO.h"
#include "itkGDCMSeriesFileNames.h"
#include "itkImageSeriesReader.h"
#include "itkImageSeriesWriter.h"
#include <itkImageToVTKImageFilter.h>
#include "itkMedianImageFilter.h"
#include "itkNumericSeriesFileNames.h"

using namespace std;
#define ONE_FRAME 1;
//typedef itk::Image<unsigned char, 2 > ImageType;
#define THREE_DIMENSION 0;//是否绘制三维


//static void CreateImage(ImageType::Pointer image);
using PixelType = int;
#if THREE_DIMENSION
constexpr unsigned int Dimension = 3;
#endif
#if !THREE_DIMENSION
constexpr unsigned int Dimension = 2;
#endif

//初始化待读取序列的格式类型
using ImageType = itk::Image<PixelType, Dimension>;
using ImageType2D = itk::Image<PixelType, 2>;
using ReaderType = itk::ImageSeriesReader<ImageType>;
using ReaderType2D = itk::ImageFileReader<ImageType2D>;


//int main() 
int Histgram()
{
	using ImageIOType = itk::GDCMImageIO;
	using NamesGeneratorType = itk::GDCMSeriesFileNames;

	

	//设置IO，并获取文件名
	ImageIOType::Pointer        gdcmIO = ImageIOType::New();
	

	ReaderType::Pointer reader = ReaderType::New();
	reader->SetImageIO(gdcmIO);

#if THREE_DIMENSION
	/********************************读取dcm序列************************************************/
	string directoryName = "D:\\input";//dcm序列读取地址
	NamesGeneratorType::Pointer namesGenerator = NamesGeneratorType::New();
	namesGenerator->SetInputDirectory(directoryName);

	const ReaderType::FileNamesContainer& filenames =
		namesGenerator->GetInputFileNames();

	std::size_t numberOfFileNames = filenames.size();
	std::cout << numberOfFileNames << std::endl;
	for (unsigned int fni = 0; fni < numberOfFileNames; ++fni)
	{
		std::cout << "filename # " << fni << " = ";
		std::cout << filenames[fni] << std::endl;
	}
	reader->SetFileNames(filenames);
#endif
#
#if !THREE_DIMENSION
/*********************************读取dcm图像***************************************************/
	reader->SetFileName("D:\\file.dcm");//自定义地址
#endif
	
	try
	{
		reader->Update();
	}
	catch (const itk::ExceptionObject& e)
	{
		std::cerr << "exception in file reader " << std::endl;
		std::cerr << e << std::endl;
		return EXIT_FAILURE;
	}

	using MedianFilterType = itk::MedianImageFilter<ImageType, ImageType>;

	MedianFilterType::Pointer filter = MedianFilterType::New();

	ImageType::SizeType indexRadius;

	indexRadius[0] = 2; // radius along x
	indexRadius[1] = 2; // radius along y
#if THREE_DIMENSION
	indexRadius[2] = 1; // radius along z
#endif
	

	filter->SetRadius(indexRadius);

	filter->SetInput(reader->GetOutput());
	filter->Update();

	int x = reader->GetOutput()->GetBufferedRegion().GetSize()[0];
	int y = reader->GetOutput()->GetBufferedRegion().GetSize()[1];
#if THREE_DIMENSION
	int z = reader->GetOutput()->GetBufferedRegion().GetSize()[2];
#endif
	

	int max = 0;
	int min = 0;
	for (int i = 1; i < x; i++) {
		for (int j = 1; j < y; j++) {
#if THREE_DIMENSION
			for (int k = 1; k < z; k++) {
#endif

				ImageType::IndexType id;
				id[0] = i;
				id[1] = j;
#if THREE_DIMENSION
				id[2] = k;
#endif
				
				ImageType::PixelType p = filter->GetOutput()->GetPixel(id);
				if (p > max) {
					max = p;
				}
				else if (p <= min) {
					min = p;
				}
#if THREE_DIMENSION
			}
#endif
		}
	}
	cout << "max is " << max << " min is " << min << endl;
	
	int* hist = new int[max + 1]{0};
	
	for (int i = 1; i < x; i++) {
		for (int j = 1; j < y; j++) {
#if THREE_DIMENSION
			for (int k = 1; k < z; k++) {
#endif

				ImageType::IndexType id;
				id[0] = i;
				id[1] = j;
#if THREE_DIMENSION
				id[2] = k;
#endif

				ImageType::PixelType p = filter->GetOutput()->GetPixel(id);
				
				if (p <= 0) {
					continue;//不统计亮度值小于0的点
				}

				hist[p]++;
#if THREE_DIMENSION
			}
#endif
			

		}
	}

	
#if !THREE_DIMENSION
	typedef itk::ImageToVTKImageFilter<ImageType> ConnectorType;
	ConnectorType::Pointer originalConnector = ConnectorType::New();
	originalConnector->SetInput(filter->GetOutput());
	//originalConnector->UpdateLargestPossibleRegion();
	originalConnector->Update();
#endif
	
	int bins = max - 0;
	//int bins = (max - min) / 10;
	int comps = 1;
	
	vtkSmartPointer<vtkIntArray> frequencies =
		vtkSmartPointer<vtkIntArray>::New();
	frequencies->SetNumberOfComponents(1);
	frequencies->SetArray(hist, max /4, 1);
	
	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.1, 0.1, 0.0);
	barChart->GetPosition2Coordinate()->SetValue(0.9, 0.9, 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->SetBarLabel(1, "Count");
	//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]); //单通道 红色
		}
	}


	double barView[4] = { 0.0, 0.0, 1.0, 1.0 };

	vtkSmartPointer<vtkRenderer> barRender =
		vtkSmartPointer<vtkRenderer>::New();
	barRender->SetViewport(barView);
	barRender->AddActor(barChart);
	barRender->SetBackground(1.0, 1.0, 1.0);

	vtkSmartPointer<vtkRenderWindow> renderWindow =
		vtkSmartPointer<vtkRenderWindow>::New();
	renderWindow->AddRenderer(barRender);

	renderWindow->SetSize(640*3, 320*3);
	renderWindow->Render();
	renderWindow->SetWindowName("Gray-Image Histogram");

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

	interactor->Initialize();
	interactor->Start();

	return 0;

}