ITK求两图像局部均方差

今天想做一个图像分析，内容是把图像分成一个个的小区域，然后两幅图像中，对应区域的均方差。最先计算出来的两图像方差始终为0，最后发现是ITK的Update方法是在调用方法结束后才会写入硬盘的。也就是说，对每一个管道，最好都单独定义的自己的滤波器。比如读入两幅图像，最好每一幅图像都单独定义一个reader。我最先只定义了一个reader，所以虽然每次都给reader写入了新的文件名，即调用了两次setfilename方法，并调用了两次update方法，最后得到的两图图像是一样的。最后修改的代码如下。

#include "itkImage.h"
#include "itkListSample.h"
#include "itkVector.h"
#include "itkCastImageFilter.h"
#include "itkMeanSquaresImageToImageMetric.h"
#include "itkTranslationTransform.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkImageFileWriter.h"
#include "itkImageFileReader.h"

typedef itk::Vector<float,256> MeasurementVectorType;
typedef itk::Image<unsigned char,2> InputImageType;
typedef itk::Image<float,2> FloatImageType;
typedef itk::CastImageFilter<InputImageType,FloatImageType> CasterFilterType;

int GetMeasureVector(FloatImageType::Pointer image1,FloatImageType::Pointer image2,MeasurementVectorType mv);

int main()
{
	typedef itk::ImageFileReader<InputImageType> ImageReaderType;
	ImageReaderType::Pointer reader1 = ImageReaderType::New();

	reader1->SetFileName("input1.bmp");
	reader1->Update();
	InputImageType::Pointer image1 = reader1->GetOutput();

	ImageReaderType::Pointer reader2 = ImageReaderType::New();
	reader2->SetFileName("input2.bmp");
	reader2->Update();
	InputImageType::Pointer image2 = reader2->GetOutput();

	MeasurementVectorType mv ;

	CasterFilterType::Pointer caster1 = CasterFilterType::New();
	caster1->SetInput(image1);
	caster1->Update();
	FloatImageType::Pointer floatImage1 = caster1->GetOutput();
	CasterFilterType::Pointer caster2 = CasterFilterType::New();
	caster2->SetInput(image2);
	caster2->Update();
	FloatImageType::Pointer floatImage2 = caster2->GetOutput();

	GetMeasureVector(floatImage1,floatImage2,mv);
	
	return EXIT_SUCCESS;
}

int GetMeasureVector(FloatImageType::Pointer image1,FloatImageType::Pointer image2,MeasurementVectorType mv)
{

	typedef itk::MeanSquaresImageToImageMetric<FloatImageType,FloatImageType> MetricType;
	typedef itk::TranslationTransform<double,2> TransformType;
	typedef itk::LinearInterpolateImageFunction<FloatImageType,double> InterpolatorType;
	MetricType::Pointer metric = MetricType::New();
	

	FloatImageType::SizeType processSize ;
	const int HorizontalSize = 32;
	const int VerticalSize = 32;
	processSize[0] = HorizontalSize;
	processSize[1] = VerticalSize;//定义区域大小
	FloatImageType::IndexType processIndex;
	processIndex[0] = 0;
	processIndex[1] = 0;
	FloatImageType::SizeType imageSize = image1->GetLargestPossibleRegion().GetSize();
    
	int HorizontalTimes = imageSize[0]/HorizontalSize;
	int VerticalTimes = imageSize[1]/VerticalSize;
	
	FloatImageType::RegionType processRegion(processIndex,processSize);
	MetricType::TransformParametersType transParam(2) ;//这里不指定维数2，将抛出异常
    transParam.Fill(0.0);

	metric->SetFixedImage(image1);
	metric->SetMovingImage(image2);
	metric->SetTransform(TransformType::New());
	InterpolatorType::Pointer interpolator = InterpolatorType::New();
	interpolator->SetInputImage(image1);
	metric->SetInterpolator(interpolator);
	
	metric->SetFixedImageRegion(processRegion);
	metric->Initialize();
	for ( int i = 0; i < HorizontalTimes; i++ )
	{
		processIndex[0] = i*HorizontalSize;
		for ( int j = 0; j < VerticalTimes; j++ )
		{
			processIndex[1] = j*VerticalSize;
		
			processRegion.SetIndex(processIndex);
			metric->SetFixedImageRegion(processRegion);
			metric->Initialize();
			mv[i*VerticalTimes+j] = metric->GetValue(transParam);
			std::cout<<"i"<<i<<","<<"j"<<j<<std::endl;
			std::cout<<mv[i*VerticalTimes+j]<<std::endl;

			//测试，通过迭代器来将不同的区域设置为不同的灰度值
			//int grayvalue = 0;
			//typedef itk::ImageRegionIterator<FloatImageType> IteratorType;
			//image1->SetRequestedRegion(processRegion);
			//IteratorType it(image1,image1->GetRequestedRegion());
			//it.GoToBegin();
			//while(!it.IsAtEnd())
			//{
			//	std::cout<<it.Get()<<std::endl;
			//	it.Set(((i+1)%2)*((j+1)%2)*255);//255是白色，0是黑色
			//	++it;//别忘了这句。
			//}
		}
	}
	
	typedef itk::ImageFileWriter<InputImageType> ImageWriterType;
	ImageWriterType::Pointer writer1 = ImageWriterType::New();
	ImageWriterType::Pointer writer2 = ImageWriterType::New();
	itk::CastImageFilter<FloatImageType,InputImageType>::Pointer caster1 = itk::CastImageFilter<FloatImageType,InputImageType>::New();
	itk::CastImageFilter<FloatImageType,InputImageType>::Pointer caster2 = itk::CastImageFilter<FloatImageType,InputImageType>::New();
	
	writer1->SetFileName("result1.bmp");//这两段代码用来将读入的图像写入到硬盘，如果前面只定义了一个reader，那么这里得到的图像都是image2
   	caster1->SetInput(image1);
	writer1->SetInput(caster1->GetOutput());
	writer1->Update();
	writer2->SetFileName("result2.bmp");//如果指定了一个writer来写出两图，那么将出错
	caster2->SetInput(image2);
	writer2->SetInput(caster2->GetOutput());
	writer2->Update();
	
	system("pause");
	return EXIT_SUCCESS;
}