ITK(Insight Tool Kit) 医疗影像数据后处理软件模块使用和算法02

本文链接：https://blog.csdn.net/weixin_43904073/article/details/120565669

上节说到如何在windows（10）的VS（2017）平台下导入ITK的library，这次咱们开始说一下咱们一开始最关心的内容，如何用ITK来做3D图像（e.g. CT）的2D投影（X-ray）。

我们先来看看结果吧，先已一个脊柱节段的CT影像为例，右边为轴向的投影。下面就具体说说怎么使用ITK。

下面主要是改编ITK标准example的代码，如果根据我上一篇的内容成功编译了HelloWorld项目，下面代码应该可以成功跑通。

首先我们以需求出发，假设我们手头上是DICOM格式的CT数据，我们希望：

“读取”（ITK模块读取并通过容器存放DICOM序列（比如以列表或者其它iterative的数据格式）） -> “处理” （通过输入的参数，ray source 位置角度，output 图像的尺寸格式，投影到平面的阈值 etc）-> “导出” （导出成固定格式的2D图像）

首先有几个概念在manual中已经提出来了，要事先了解一下:

ITK的IO framework：sources（Associated with the data processing pipeline，reader）；初始化filter（initiate the pipeline, and mappers,）；终止filter（terminate the pipeline）；mapper （output data from pipeline，writer）；

ITK创建图像对象的时候，会创建一个智能指针（Pointer /SmartPointer），官方说明作用就是更好的管理内存，因为图像都比较大，根据使用情况（reference counting，受到referred之后，进行register（），如果count==0，则delete（）），及时释放内容，以免出现内存溢出的情况。

图像的读取

图片类型定义, 如果定义的格式不够，系统会隐式转化：

using ImageType = itk::Image<PixelType, Dimension>;

图片读取类定义（reader）：

using ReaderType = itk::ImageFileReader<ImageType>;

初始化reader类，每次读取后，智能指针会创建：

ReaderType::Pointer reader = ReaderType::New();

为了保证唯一性，至少每个reader pointer都会有一个名字：

reader->SetFileName(filename);

并进行更新，保证在pipeline中该image的reader class被处理：

reader->Update();

获取reader pipeline最新的一张图片（image reader），特别是在update（）之后；也可以在update之前使用，如果pipeline是空，会等到reader处理后进行output：

ImageType::Pointer image = reader->GetOutput();

举一反三，如果是DICOM series如下，就应该是每个图片生成对应的image reader类，幸好ITK已经把它们放在#include "itkGDCMImageIO.h" 不用自己写了，#include "itkGDCMSeriesFileNames.h" 用于生成批量的image name，因为前面已经说了，filename是image reader的标识物。（放在一个类似于[]的数据结构中）。图像读取后，对于图像中区域的计算和像素的处理，也应该比较简单了吧~

	if (input_name)
	{
		//typedef itk::ImageFileReader< InputImageType >  ReaderType;
		//ReaderType::Pointer reader = ReaderType::New();
		//reader->SetFileName(input_name);

		
		using ImageIOType = itk::GDCMImageIO;
		using NamesGeneratorType = itk::GDCMSeriesFileNames;
		ImageIOType::Pointer gdcmIO = ImageIOType::New();
		NamesGeneratorType::Pointer namesGenerator = NamesGeneratorType::New();
		namesGenerator->SetInputDirectory(input_name);

		using ReaderType = itk::ImageSeriesReader<InputImageType>;
		const ReaderType::FileNamesContainer& filenames = namesGenerator->GetInputFileNames();

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

		try
		{
			reader->Update();
		}
		catch (itk::ExceptionObject & err)
		{
			std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
			std::cerr << err << std::endl;
			return EXIT_FAILURE;
		}
		image = reader->GetOutput();

	}

//因为内容比较多
//处理和导出的部分再加一个章节来写吧~~

/*=========================================================================
 *
 *  Copyright Insight Software Consortium
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
 // Software Guide : BeginLatex
 //
 // This example illustrates the use of the
 // \code{RayCastInterpolateImageFunction} class to generate digitally
 // reconstructed radiographs (DRRs) from a 3D image volume such as CT
 // or MR.
 //
 // Software Guide : EndLatex
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCenteredEuler3DTransform.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkRescaleIntensityImageFilter.h"
// Software Guide : BeginLatex
//
// The \code{RayCastInterpolateImageFunction} class definition for
// this example is contained in the following header file.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkRayCastInterpolateImageFunction.h"
// Software Guide : EndCodeSnippet
//#define WRITE_CUBE_IMAGE_TO_FILE

#include "itkGDCMImageIO.h"
#include "itkGDCMSeriesFileNames.h"  
#include "itkImageSeriesWriter.h"
#include "itkImageSeriesReader.h" 
#include "itkNumericSeriesFileNames.h"

void usage()
{
	std::cerr << "\n";
	std::cerr << "Usage: DRR <options> [input]\n";
	std::cerr << "       calculates the Digitally Reconstructed Radiograph from  \n";
	std::cerr << "       a volume. \n\n";
	std::cerr << "   where <options> is one or more of the following:\n\n";
	std::cerr << "       <-h>                    Display (this) usage information\n";
	std::cerr << "       <-v>                    Verbose output [default: no]\n";
	std::cerr << "       <-res float float>      Pixel spacing of the output image [default: 1x1mm]  \n";
	std::cerr << "       <-size int int>         Dimension of the output image [default: 501x501]  \n";
	std::cerr << "       <-sid float>            Distance of ray source (focal point) [default: 400mm]\n";
	std::cerr << "       <-t float float float>  Translation parameter of the camera \n";
	std::cerr << "       <-rx float>             Rotation around x,y,z axis in degrees \n";
	std::cerr << "       <-ry float>\n";
	std::cerr << "       <-rz float>\n";
	std::cerr << "       <-normal float float>   The 2D projection normal position [default: 0x0mm]\n";
	std::cerr << "       <-cor float float float> The centre of rotation relative to centre of volume\n";
	std::cerr << "       <-threshold float>      Threshold [default: 0]\n";
	std::cerr << "       <-o file>               Output image filename\n\n";
	std::cerr << "                               by  thomas@hartkens.de\n";
	std::cerr << "                               and john.hipwell@kcl.ac.uk (CISG London)\n\n";
	exit(1);
}
int main(int argc, char *argv[])
{
	//reading mhd file
	//char *input_name = "输入3D数据文件夹目录";
	char *output_name = "输出2D数据\\12.png";
	bool ok;
	bool verbose = true;
	float rx = 0.;
	float ry = 0.;
	float rz = 0.;
	float tx = 0.;
	float ty = 0.;
	float tz = 0.;
	float cx = 0.;
	float cy = 0.;
	float cz = 0.;
	float sid = 500.;
	float sx = 0.6;
	float sy = 0.6;
	int dx = 512;
	int dy = 512;
	float o2Dx = 0;
	float o2Dy = 0;
	double threshold = -100.;
	// Parse command line parameters
	while (argc > 1)
	{
		ok = false;
		if ((ok == false) && (strcmp(argv[1], "-h") == 0))
		{
			argc--; argv++;
			ok = true;
			usage();
		}
		if ((ok == false) && (strcmp(argv[1], "-v") == 0))
		{
			argc--; argv++;
			ok = true;
			verbose = true;
		}
		if ((ok == false) && (strcmp(argv[1], "-rx") == 0))
		{
			argc--; argv++;
			ok = true;
			rx = atof(argv[1]);
			argc--; argv++;
		}
		if ((ok == false) && (strcmp(argv[1], "-ry") == 0))
		{
			argc--; argv++;
			ok = true;
			ry = atof(argv[1]);
			argc--; argv++;
		}
		if ((ok == false) && (strcmp(argv[1], "-rz") == 0))
		{
			argc--; argv++;
			ok = true;
			rz = atof(argv[1]);
			argc--; argv++;
		}
		if ((ok == false) && (strcmp(argv[1], "-threshold") == 0))
		{
			argc--; argv++;
			ok = true;
			threshold = atof(argv[1]);
			argc--; argv++;
		}
		if ((ok == false) && (strcmp(argv[1], "-t") == 0))
		{
			argc--; argv++;
			ok = true;
			tx = atof(argv[1]);
			argc--; argv++;
			ty = atof(argv[1]);
			argc--; argv++;
			tz = atof(argv[1]);
			argc--; argv++;
		}
		if ((ok == false) && (strcmp(argv[1], "-cor") == 0))
		{
			argc--; argv++;
			ok = true;
			cx = atof(argv[1]);
			argc--; argv++;
			cy = atof(argv[1]);
			argc--; argv++;
			cz = atof(argv[1]);
			argc--; argv++;
		}
		if ((ok == false) && (strcmp(argv[1], "-res") == 0))
		{
			argc--; argv++;
			ok = true;
			sx = atof(argv[1]);
			argc--; argv++;
			sy = atof(argv[1]);
			argc--; argv++;
		}
		if ((ok == false) && (strcmp(argv[1], "-size") == 0))
		{
			argc--; argv++;
			ok = true;
			dx = atoi(argv[1]);
			argc--; argv++;
			dy = atoi(argv[1]);
			argc--; argv++;
		}
		if ((ok == false) && (strcmp(argv[1], "-sid") == 0))
		{
			argc--; argv++;
			ok = true;
			sid = atof(argv[1]);
			argc--; argv++;
		}
		if ((ok == false) && (strcmp(argv[1], "-normal") == 0))
		{
			argc--; argv++;
			ok = true;
			o2Dx = atof(argv[1]);
			argc--; argv++;
			o2Dy = atof(argv[1]);
			argc--; argv++;
		}
		if ((ok == false) && (strcmp(argv[1], "-o") == 0))
		{
			argc--; argv++;
			ok = true;
			output_name = argv[1];
			argc--; argv++;
		}
		if (ok == false)
		{
			if (input_name == NULL)
			{
				input_name = argv[1];
				argc--;
				argv++;
			}
			else
			{
				std::cerr << "ERROR: Can not parse argument " << argv[1] << std::endl;
				usage();
			}
		}
	}
	if (verbose)
	{
		if (input_name)  std::cout << "Input image: " << input_name << std::endl;
		if (output_name) std::cout << "Output image: " << output_name << std::endl;
	}
	// Software Guide : BeginLatex
	//
	// Although we generate a 2D projection of the 3D volume for the
	// purposes of the interpolator both images must be three dimensional.
	//
	// Software Guide : EndLatex
	// Software Guide : BeginCodeSnippet
	const     unsigned int   Dimension = 3;
	typedef   short         InputPixelType;
	typedef   unsigned char OutputPixelType;
	typedef itk::Image< InputPixelType, Dimension >   InputImageType;
	typedef itk::Image< OutputPixelType, Dimension >   OutputImageType;
	InputImageType::Pointer image;
	// Software Guide : EndCodeSnippet
	// Software Guide : BeginLatex
	//
	// For the purposes of this example we assume the input volume has
	// been loaded into an \code{itk::Image image}.
	//
	// Software Guide : EndLatex
	if (input_name)
	{
		//typedef itk::ImageFileReader< InputImageType >  ReaderType;
		//ReaderType::Pointer reader = ReaderType::New();
		//reader->SetFileName(input_name);

		
		using ImageIOType = itk::GDCMImageIO;
		using NamesGeneratorType = itk::GDCMSeriesFileNames;
		ImageIOType::Pointer gdcmIO = ImageIOType::New();
		NamesGeneratorType::Pointer namesGenerator = NamesGeneratorType::New();
		namesGenerator->SetInputDirectory(input_name);

		using ReaderType = itk::ImageSeriesReader<InputImageType>;
		const ReaderType::FileNamesContainer& filenames = namesGenerator->GetInputFileNames();

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

		try
		{
			reader->Update();
		}
		catch (itk::ExceptionObject & err)
		{
			std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
			std::cerr << err << std::endl;
			return EXIT_FAILURE;
		}
		image = reader->GetOutput();

	}
	else
	{   // No input image specified so create a cube
		image = InputImageType::New();
		InputImageType::SpacingType spacing;
		spacing[0] = 3.;
		spacing[1] = 3.;
		spacing[2] = 3.;
		image->SetSpacing(spacing);
		InputImageType::PointType origin;
		origin[0] = 0.;
		origin[1] = 0.;
		origin[2] = 0.;
		image->SetOrigin(origin);
		InputImageType::IndexType start;
		start[0] = 0;  // first index on X
		start[1] = 0;  // first index on Y
		start[2] = 0;  // first index on Z
		InputImageType::SizeType  size;
		size[0] = 61;  // size along X
		size[1] = 61;  // size along Y
		size[2] = 61;  // size along Z
		InputImageType::RegionType region;
		region.SetSize(size);
		region.SetIndex(start);
		image->SetRegions(region);
		image->Allocate(true); // initialize to zero.
		image->Update();
		typedef itk::ImageRegionIteratorWithIndex< InputImageType > IteratorType;
		IteratorType iterate(image, image->GetLargestPossibleRegion());
		while (!iterate.IsAtEnd())
		{
			InputImageType::IndexType idx = iterate.GetIndex();
			if ((idx[0] >= 6) && (idx[0] <= 54)
				&& (idx[1] >= 6) && (idx[1] <= 54)
				&& (idx[2] >= 6) && (idx[2] <= 54)
				&& ((((idx[0] <= 11) || (idx[0] >= 49))
					&& ((idx[1] <= 11) || (idx[1] >= 49)))
					|| (((idx[0] <= 11) || (idx[0] >= 49))
						&& ((idx[2] <= 11) || (idx[2] >= 49)))
					|| (((idx[1] <= 11) || (idx[1] >= 49))
						&& ((idx[2] <= 11) || (idx[2] >= 49)))))
			{
				iterate.Set(10);
			}
			else if ((idx[0] >= 18) && (idx[0] <= 42)
				&& (idx[1] >= 18) && (idx[1] <= 42)
				&& (idx[2] >= 18) && (idx[2] <= 42)
				&& ((((idx[0] <= 23) || (idx[0] >= 37))
					&& ((idx[1] <= 23) || (idx[1] >= 37)))
					|| (((idx[0] <= 23) || (idx[0] >= 37))
						&& ((idx[2] <= 23) || (idx[2] >= 37)))
					|| (((idx[1] <= 23) || (idx[1] >= 37))
						&& ((idx[2] <= 23) || (idx[2] >= 37)))))
			{
				iterate.Set(60);
			}
			else if ((idx[0] == 30) && (idx[1] == 30) && (idx[2] == 30))
			{
				iterate.Set(100);
			}
			++iterate;
		}
#ifdef WRITE_CUBE_IMAGE_TO_FILE
		char *filename = "cube.gipl";
		typedef itk::ImageFileWriter< OutputImageType >  WriterType;
		WriterType::Pointer writer = WriterType::New();
		writer->SetFileName(filename);
		writer->SetInput(image);
		try
		{
			std::cout << "Writing image: " << filename << std::endl;
			writer->Update();
		}
		catch (itk::ExceptionObject & err)
		{
			std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
			std::cerr << err << std::endl;
			return EXIT_FAILURE;
		}
#endif
	}
	// Print out the details of the input volume
	if (verbose)
	{
		unsigned int i;
		const InputImageType::SpacingType spacing = image->GetSpacing();
		std::cout << std::endl << "Input ";
		InputImageType::RegionType region = image->GetBufferedRegion();
		region.Print(std::cout);
		std::cout << "  Resolution: [";
		for (i = 0; i < Dimension; i++)
		{
			std::cout << spacing[i];
			if (i < Dimension - 1) std::cout << ", ";
		}
		std::cout << "]" << std::endl;
		const InputImageType::PointType origin = image->GetOrigin();
		std::cout << "  Origin: [";
		for (i = 0; i < Dimension; i++)
		{
			std::cout << origin[i];
			if (i < Dimension - 1) std::cout << ", ";
		}
		std::cout << "]" << std::endl << std::endl;
	}
	// Software Guide : BeginLatex
	//
	// Creation of a \code{ResampleImageFilter} enables coordinates for
	// each of the pixels in the DRR image to be generated. These
	// coordinates are used by the \code{RayCastInterpolateImageFunction}
	// to determine the equation of each corresponding ray which is cast
	// through the input volume.
	//
	// Software Guide : EndLatex
	// Software Guide : BeginCodeSnippet
	typedef itk::ResampleImageFilter<InputImageType, InputImageType > FilterType;
	FilterType::Pointer filter = FilterType::New();
	filter->SetInput(image);
	filter->SetDefaultPixelValue(0);
	// Software Guide : EndCodeSnippet
	// Software Guide : BeginLatex
	//
	// An Euler transformation is defined to position the input volume.
	// The \code{ResampleImageFilter} uses this transform to position the
	// output DRR image for the desired view.
	//
	// Software Guide : EndLatex
	// Software Guide : BeginCodeSnippet
	typedef itk::CenteredEuler3DTransform< double >  TransformType;
	TransformType::Pointer transform = TransformType::New();
	transform->SetComputeZYX(true);
	TransformType::OutputVectorType translation;
	translation[0] = tx;
	translation[1] = ty;
	translation[2] = tz;
	// constant for converting degrees into radians
	const double dtr = (std::atan(1.0) * 4.0) / 180.0;
	transform->SetTranslation(translation);
	transform->SetRotation(dtr*rx, dtr*ry, dtr*rz);
	InputImageType::PointType   imOrigin = image->GetOrigin();
	InputImageType::SpacingType imRes = image->GetSpacing();
	typedef InputImageType::RegionType     InputImageRegionType;
	typedef InputImageRegionType::SizeType InputImageSizeType;
	InputImageRegionType imRegion = image->GetBufferedRegion();
	InputImageSizeType   imSize = imRegion.GetSize();
	imOrigin[0] += imRes[0] * static_cast<double>(imSize[0]) / 2.0;
	imOrigin[1] += imRes[1] * static_cast<double>(imSize[1]) / 2.0;
	imOrigin[2] += imRes[2] * static_cast<double>(imSize[2]) / 2.0;
	TransformType::InputPointType center;
	center[0] = cx + imOrigin[0];
	center[1] = cy + imOrigin[1];
	center[2] = cz + imOrigin[2];
	transform->SetCenter(center);
	if (verbose)
	{
		std::cout << "Image size: "
			<< imSize[0] << ", " << imSize[1] << ", " << imSize[2]
			<< std::endl << "   resolution: "
			<< imRes[0] << ", " << imRes[1] << ", " << imRes[2]
			<< std::endl << "   origin: "
			<< imOrigin[0] << ", " << imOrigin[1] << ", " << imOrigin[2]
			<< std::endl << "   center: "
			<< center[0] << ", " << center[1] << ", " << center[2]
			<< std::endl << "Transform: " << transform << std::endl;
	}
	// Software Guide : EndCodeSnippet
	// Software Guide : BeginLatex
	//
	// The \code{RayCastInterpolateImageFunction} is instantiated and passed the transform
	// object. The \code{RayCastInterpolateImageFunction} uses this
	// transform to reposition the x-ray source such that the DRR image
	// and x-ray source move as one around the input volume. This coupling
	// mimics the rigid geometry of the x-ray gantry.
	//
	// Software Guide : EndLatex
	// Software Guide : BeginCodeSnippet
	typedef itk::RayCastInterpolateImageFunction<InputImageType, double>
		InterpolatorType;
	InterpolatorType::Pointer interpolator = InterpolatorType::New();
	interpolator->SetTransform(transform);
	// Software Guide : EndCodeSnippet
	// Software Guide : BeginLatex
	//
	// We can then specify a threshold above which the volume's
	// intensities will be integrated.
	//
	// Software Guide : EndLatex
	// Software Guide : BeginCodeSnippet
	interpolator->SetThreshold(threshold);
	// Software Guide : EndCodeSnippet
	// Software Guide : BeginLatex
	//
	// The ray-cast interpolator needs to know the initial position of the
	// ray source or focal point. In this example we place the input
	// volume at the origin and halfway between the ray source and the
	// screen. The distance between the ray source and the screen
	// is the "source to image distance" \code{sid} and is specified by
	// the user.
	//
	// Software Guide : EndLatex
	// Software Guide : BeginCodeSnippet
	InterpolatorType::InputPointType focalpoint;
	focalpoint[0] = imOrigin[0];
	focalpoint[1] = imOrigin[1];
	focalpoint[2] = imOrigin[2] - sid / 2.;
	interpolator->SetFocalPoint(focalpoint);
	// Software Guide : EndCodeSnippet
	if (verbose)
	{
		std::cout << "Focal Point: "
			<< focalpoint[0] << ", "
			<< focalpoint[1] << ", "
			<< focalpoint[2] << std::endl;
	}
	// Software Guide : BeginLatex
	//
	// Having initialised the interpolator we pass the object to the
	// resample filter.
	//
	// Software Guide : EndLatex
	// Software Guide : BeginCodeSnippet
	interpolator->Print(std::cout);
	filter->SetInterpolator(interpolator);
	filter->SetTransform(transform);
	// Software Guide : EndCodeSnippet
	// Software Guide : BeginLatex
	//
	// The size and resolution of the output DRR image is specified via the
	// resample filter.
	//
	// Software Guide : EndLatex
	// Software Guide : BeginCodeSnippet
	  // setup the scene
	InputImageType::SizeType   size;
	size[0] = dx;  // number of pixels along X of the 2D DRR image
	size[1] = dy;  // number of pixels along Y of the 2D DRR image
	size[2] = 1;   // only one slice
	filter->SetSize(size);
	InputImageType::SpacingType spacing;
	spacing[0] = sx;  // pixel spacing along X of the 2D DRR image [mm]
	spacing[1] = sy;  // pixel spacing along Y of the 2D DRR image [mm]
	spacing[2] = 1.0; // slice thickness of the 2D DRR image [mm]
	filter->SetOutputSpacing(spacing);
	// Software Guide : EndCodeSnippet
	if (verbose)
	{
		std::cout << "Output image size: "
			<< size[0] << ", "
			<< size[1] << ", "
			<< size[2] << std::endl;
		std::cout << "Output image spacing: "
			<< spacing[0] << ", "
			<< spacing[1] << ", "
			<< spacing[2] << std::endl;
	}
	// Software Guide : BeginLatex
	//
	// In addition the position of the DRR is specified. The default
	// position of the input volume, prior to its transformation is
	// half-way between the ray source and screen and unless specified
	// otherwise the normal from the "screen" to the ray source passes
	// directly through the centre of the DRR.
	//
	// Software Guide : EndLatex
	// Software Guide : BeginCodeSnippet
	double origin[Dimension];
	origin[0] = imOrigin[0] + o2Dx - sx * ((double)dx - 1.) / 2.;
	origin[1] = imOrigin[1] + o2Dy - sy * ((double)dy - 1.) / 2.;
	origin[2] = imOrigin[2] + sid / 2.;
	filter->SetOutputOrigin(origin);
	// Software Guide : EndCodeSnippet
	if (verbose)
	{
		std::cout << "Output image origin: "
			<< origin[0] << ", "
			<< origin[1] << ", "
			<< origin[2] << std::endl;
	}
	// create writer
	if (output_name)
	{
		// Software Guide : BeginLatex
		//
		// The output of the resample filter can then be passed to a writer to
		// save the DRR image to a file.
		//
		// Software Guide : EndLatex
		// Software Guide : BeginCodeSnippet
		typedef itk::RescaleIntensityImageFilter<
			InputImageType, OutputImageType > RescaleFilterType;
		RescaleFilterType::Pointer rescaler = RescaleFilterType::New();
		rescaler->SetOutputMinimum(0);
		rescaler->SetOutputMaximum(255);
		rescaler->SetInput(filter->GetOutput());
		typedef itk::ImageFileWriter< OutputImageType >  WriterType;
		WriterType::Pointer writer = WriterType::New();
		writer->SetFileName(output_name);
		writer->SetInput(rescaler->GetOutput());
		try
		{
			std::cout << "Writing image: " << output_name << std::endl;
			writer->Update();
		}
		catch (itk::ExceptionObject & err)
		{
			std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
			std::cerr << err << std::endl;
		}
		// Software Guide : EndCodeSnippet
	}
	else
	{
		filter->Update();
	}
	return EXIT_SUCCESS;
}