今天想做一個影象分析，內容是把影象分成一個個的小區域，然後兩幅影象中，對應區域的均方差。最先計算出來的兩影象方差始終為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;
}