[Insight-developers] Threads and WarpImageFilter

[Vincent A. Magnotta]

--=-XSSh52yu3rhFumce6er3
Content-Type: text/plain
Content-Transfer-Encoding: 7bit

I am trying to load deformation fields from a file (actually three
files) and then apply the deformation fields to an image. I started out
by modifying the Thirion Demons image registration application and
saving the deformation fields as three analyze files (one for each
direction x, y and z). If I use the WarpImageFilter in this program I
can warp the images and the output looks fine (register+warp.jpg). If I
load the deformation fields from disk and warp the image the result
looks like noise gets into the data (warp-only.jpg). If I make set the
maximum number of threads to be 1 before warping my output image does
not have this artifact. Does anyone have ideas why this may occur?

Attached is a little stand along program that contains the warping
procedure that exhibits the problem. 

Thanks,
Vince

-- 
Assistant Professor
Department of Radiology
0437 JCP
E-mail: vincent-magnotta at uiowa . edu
Phone: 319-356-8255
Fax: 319-353-6275
Website: http://www . radiology . uiowa . edu


--=-XSSh52yu3rhFumce6er3
Content-Disposition: attachment; filename=TestRegistration.cxx
Content-Type: text/plain; name=TestRegistration.cxx; charset=UTF-8
Content-Transfer-Encoding: 7bit

#include <iostream>
#include <fstream>
#include <stdio.h>
#include "itkImage.h"
#include "itkExceptionObject.h"
#include "itkMultiThreader.h"
#include "itkCommand.h"
#include "itkVectorIndexSelectionCastImageFilter.h"
#include "itkImageFileWriter.h"
#include "itkImageFileReader.h"
#include "itkWarpImageFilter.h"
#include "itkImage.h"
#include "itkComposeRGBImageFilter.h"
#include <itkAdaptImageFilter.h>
#include <itkRGBPixel.h>
#include <itkRGBToVectorPixelAccessor.h>
#include "itkRGBToVectorImageAdaptor.h"
#include "itkMetaDataObject.h"
#include "itkIOCommon.h"
#include "itkStatisticsImageFilter.h"


int main (int argc, char **argv)
{
	
	std::string DeformFieldName(argv[1]);
	
	itk::ImageFileReader<itk::Image<float, 3> >::Pointer xDisplacementImageReader=itk::ImageFileReader<itk::Image<float, 3> >::New();
	std::string CurrentComponentFilename;

	CurrentComponentFilename=DeformFieldName+"_xdisp.img";
	xDisplacementImageReader->SetFileName(CurrentComponentFilename.c_str());
	xDisplacementImageReader->Update();				
	itk::Image<float, 3>::Pointer xDisplacementImage = xDisplacementImageReader->GetOutput();
	std::cout << "Load X deformation field: " << CurrentComponentFilename << std::endl;

	itk::ImageFileReader<itk::Image<float, 3> >::Pointer yDisplacementImageReader=itk::ImageFileReader<itk::Image<float, 3> >::New();
	CurrentComponentFilename=DeformFieldName+"_ydisp.img";
	yDisplacementImageReader->SetFileName(CurrentComponentFilename.c_str());
	yDisplacementImageReader->Update();				
	itk::Image<float, 3>::Pointer yDisplacementImage = yDisplacementImageReader->GetOutput();
	std::cout << "Load Y deformation field: " << CurrentComponentFilename << std::endl;


	itk::ImageFileReader<itk::Image<float, 3> >::Pointer zDisplacementImageReader=itk::ImageFileReader<itk::Image<float, 3> >::New();
	CurrentComponentFilename=DeformFieldName+"_zdisp.img";
	zDisplacementImageReader->SetFileName(CurrentComponentFilename.c_str());
	zDisplacementImageReader->Update();				
	itk::Image<float, 3>::Pointer zDisplacementImage = zDisplacementImageReader->GetOutput();
	std::cout << "Load Z deformation field: " << CurrentComponentFilename << std::endl;



	typedef itk::ComposeRGBImageFilter< itk::Image<float, 3> >  ComposeFilterType;
    ComposeFilterType::Pointer composeRGBfilter = ComposeFilterType::New();
	typedef ComposeFilterType::OutputImageType  RGBImageType;

	composeRGBfilter->SetInput1(xDisplacementImage);
	composeRGBfilter->SetInput2(yDisplacementImage);
	composeRGBfilter->SetInput3(zDisplacementImage);

	typedef itk::Accessor::RGBToVectorPixelAccessor<float>    AccessorType;
	typedef AccessorType::ExternalType              		  VectorPixelType;
	typedef itk::Image< VectorPixelType,   3 >      		  VectorImageType;
	typedef itk::RGBPixel< float >                  		  RGBPixelType;
	typedef itk::AdaptImageFilter< RGBImageType,
					    		   VectorImageType,
					    		   AccessorType       >  RGBToVectorAdaptorType;
	
	RGBToVectorAdaptorType::Pointer rgbToVectorFilter = RGBToVectorAdaptorType::New();

	rgbToVectorFilter->SetInput(composeRGBfilter->GetOutput());				
	rgbToVectorFilter->Update();

	VectorImageType::Pointer DeformationField = rgbToVectorFilter->GetOutput();


	typedef itk::Image<unsigned char, 3> 					MovingImageType;
	typedef itk::WarpImageFilter< MovingImageType, MovingImageType, VectorImageType > WarperType;
	WarperType::Pointer warper = WarperType::New();

	itk::ImageFileReader< MovingImageType >::Pointer MovingImageReader=itk::ImageFileReader< MovingImageType >::New();
	MovingImageReader->SetFileName( argv[2] );
	MovingImageReader->Update();		
	
	MovingImageType::Pointer MovingImage = MovingImageReader->GetOutput();

	warper->SetInput( MovingImage );
	warper->SetDeformationField( DeformationField );
	warper->SetOutputOrigin( DeformationField->GetOrigin() );
	warper->SetOutputSpacing( DeformationField->GetSpacing() );
	warper->Update();
	
	MovingImageType::Pointer TempImageSptr=warper->GetOutput();
	TempImageSptr->SetMetaDataDictionary(MovingImage->GetMetaDataDictionary());
	typedef itk::ImageFileWriter<MovingImageType> ImageWriterType;
	ImageWriterType::Pointer Imagewriter = ImageWriterType::New();
	Imagewriter->SetFileName( argv[3] );
	Imagewriter->SetInput(TempImageSptr);
	Imagewriter->Update();
	
	return (0);
		
}
			


--=-XSSh52yu3rhFumce6er3--