[ITK] [ITK-users] Generating a warping field with ThinPlateSplineKernelTransform

Andras Lasso lasso at queensu.ca
Thu Aug 6 13:33:09 EDT 2015


Have you tested if the registration works well with those landmarks? I would suggest to first experiment with the registration in an interactive environment where you can see the registration result in real-time as you define/move landmarks (for example, http://www.slicer.org/slicerWiki/index.php/Documentation/Nightly/Modules/LandmarkRegistration). After you've figured out how many landmarks, at what distance, what configuration are needed for a good-quality alignment you can dive into coding your own implementation.

Andras

-----Original Message-----
From: Community [mailto:community-bounces at itk.org] On Behalf Of Olivier Comtois
Sent: Thursday, August 06, 2015 11:54 AM
To: insight-users at itk.org
Subject: Re: [ITK] [ITK-users] Generating a warping field with ThinPlateSplineKernelTransform

Thank you for your reply Siavash,

Here is my code :

It differs from the example in that I take my landmark points from 2 files (one for the source landmarks and another for the target). The source landmarks come from the segmentation of a curved spinal cord while the straight ones come from a template (which is straightened).

std::vector<double> ReadLabelValueFromFile(std::string filename) {
	std::ifstream file(filename.c_str());
	string line;
	char* token;
	std::vector<double> result;
	std::vector<double>::iterator it;

	it = result.begin();
	if (file.is_open()) {
		while (file.good()) {
			getline(file, line);
			char * charLine = new char[line.length() + 1];
			strcpy(charLine, line.c_str());
			token = std::strtok(charLine, ",");
			while (token != NULL) {
				result.push_back(std::atof(token));
				token = strtok(NULL, ",");
			}
			delete[] charLine;
		}
		file.close();
	} else {
		throw std::runtime_error(std::string("Cannot open file : ") + filename);
	}
	return result;
}

void GetRealValuePointSetFromFile(itk::ThinPlateSplineKernelTransform<
double, 3>::PointSetType::Pointer &source,
		itk::ThinPlateSplineKernelTransform< double,  
3>::PointSetType::Pointer &destination,
		string curvedFilename, string straightFilename) {
	//define variables needed for this function
	string line;
	string token;
	std::vector<double> straightValues;
	std::vector<double> curvedValues;

	//fetching coordinates from files
	try {
		straightValues = ReadLabelValueFromFile(straightFilename);
		curvedValues = ReadLabelValueFromFile(curvedFilename);
		//check if not empty
		if (curvedValues.size() == 0) {
			throw std::runtime_error(
					"No value found in the curved labels text file.");
		}
		if (straightValues.size() == 0) {
			throw std::runtime_error(
					"No value found in straight the labels text file.");
		}
		if ((straightValues.size() % 3) != 0) {
			throw std::runtime_error(
					"The straight file size must be divisible by 3 (since there is 3 dimensions to the image)");
		}
		if ((curvedValues.size() % 3) != 0) {
			throw std::runtime_error(
					"The curved file size must be divisible by 3 (since there is 3 dimensions to the image)");
		}
		if (straightValues.size() != curvedValues.size()) {
			throw std::runtime_error(
					"Curved and Straight files do not have the same size.");
		}
	} catch (const std::exception &e) {
		throw;
	}

	//Produce the Straight points
	typedef double Realtype;
	typedef itk::ThinPlateSplineKernelTransform< double, ImageDimension> TransformType;
	typedef TransformType::PointSetType PointSetType;

	std::vector<double>::iterator itS = straightValues.begin();
	std::vector<double>::iterator itC = curvedValues.begin();
	typedef PointSetType::PointIdentifier PointIdType;
	PointIdType id = itk::NumericTraits< PointIdType >::ZeroValue();

	while (itS != straightValues.end()) {
		typedef itk::Point<double, ImageDimension> PointType;
		PointType pS;
		PointType pC;
		for (int i = 0; i < 3; i++) {
			pS[i] = *itS;
			pC[i] = *itC;
			itS++;
			itC++;
		}
		//Place points in the point sets
		source->SetPoint(id, pC);
		source->SetPointData(id, id);
		destination->SetPoint(id, pS);
		destination->SetPointData(id, id);
		id++;
	}
	printf("RETURNING POINT VALUES");
	return;

}

int LandmarkBasedWithTextDisplacementFieldTransformInitializer(int
argc, char *argv[]) {


	////////////////////////////////////////////////////////////////////////////////////////////
	/*CREATING USEFUL TYPES*/
	////////////////////////////////////////////////////////////////////////////////////////////
	printf("\nCREATING USEFUL TYPES\n");
	typedef double RealType;
	typedef unsigned int LabelType;
	typedef itk::Image<LabelType, ImageDimension> LabelImageType;
	typedef itk::Image<RealType, ImageDimension> RealImageType;
	typedef itk::Vector<RealType, ImageDimension> VectorType;
	typedef itk::Vector< double, ImageDimension > FieldVectorType;
	typedef itk::ResampleImageFilter<LabelImageType, LabelImageType> ResamplerType;
	typedef itk::LinearInterpolateImageFunction<LabelImageType, RealType  > InterpolatorType;
	typedef itk::ImageFileWriter< LabelImageType >  DeformedImageWriterType;
	typedef itk::Image< FieldVectorType,  ImageDimension > DisplacementFieldType;
	typedef itk::ImageFileWriter< DisplacementFieldType > FieldWriterType;
	typedef itk::ThinPlateSplineKernelTransform< RealType,  
ImageDimension> TransformType;
	typedef itk::ImageFileReader<LabelImageType> ImageReaderType;
	typedef TransformType::PointSetType PointSetType;
	typedef PointSetType::PointIdentifier PointIdType;
	////////////////////////////////////////////////////////////////////////////////////////////

	////////////////////////////////////////////////////////////////////////////////////////////
	/*LOADING FIXED IMAGE*/
	////////////////////////////////////////////////////////////////////////////////////////////
	printf("\nLOADING DESTINATION IMAGE\n");

	ImageReaderType::Pointer destinationReader = ImageReaderType::New();
	destinationReader->SetFileName(argv[2]);
	destinationReader->Update();
	LabelImageType::Pointer destinationImage = destinationReader->GetOutput();
	////////////////////////////////////////////////////////////////////////////////////////////

	////////////////////////////////////////////////////////////////////////////////////////////
	/*LOADING MOVING IMAGE*/
	////////////////////////////////////////////////////////////////////////////////////////////
	printf("\nLOADING SOURCE IMAGE\n");

	ImageReaderType::Pointer sourceReader = ImageReaderType::New();
	sourceReader->SetFileName(argv[1]);
	sourceReader->Update();
	LabelImageType::Pointer sourceImage = sourceReader->GetOutput();
	////////////////////////////////////////////////////////////////////////////////////////////

	////////////////////////////////////////////////////////////////////////////////////////////
	/*INIT IMAGE CONTAINERS*/
	////////////////////////////////////////////////////////////////////////////////////////////
	typedef itk::PointSet<RealType, ImageDimension> RealPointSetType;
	PointSetType::Pointer destinationPts = PointSetType::New();
	PointSetType::Pointer sourcePts = PointSetType::New();
	RegistrationType::Pointer registration = RegistrationType::New();

	destinationPts->Initialize();
	sourcePts->Initialize();
	////////////////////////////////////////////////////////////////////////////////////////////

	////////////////////////////////////////////////////////////////////////////////////////////
	/*IMPORTING LANDMARKS FROM TEXT FILES*/
	////////////////////////////////////////////////////////////////////////////////////////////
	printf("\nIMPORTING LANDMARKS FROM TEXT FILES\n");
	if (string(argv[3]).compare("") != 0) {
		if (string(argv[4]).compare("") != 0) {
			try {
				string sourceFilename = argv[3];
				string destinationFilename = argv[4];
				GetRealValuePointSetFromFile(sourcePts, destinationPts, sourceFilename, destinationFilename);
			} catch (exception &e) {
				throw;
			}
		} else
			throw std::runtime_error(
					"Both arguments the straight and curved file must be specified");
	} else
		throw std::runtime_error(
				"A file must be specified in order to take input real labels.");
	////////////////////////////////////////////////////////////////////////////////////////////

	////////////////////////////////////////////////////////////////////////////////////////////
	/*DEFINING KERNEL TRANSFORM*/
	////////////////////////////////////////////////////////////////////////////////////////////
	printf("\nDEFINING KERNEL TRANSFORM\n");
	TransformType::Pointer tps = TransformType::New();

	tps->SetSourceLandmarks(sourcePts);
	tps->SetTargetLandmarks(destinationPts);
	tps->ComputeWMatrix();
	////////////////////////////////////////////////////////////////////////////////////////////

	////////////////////////////////////////////////////////////////////////////////////////////
	/*DEFINING RESAMPLER TO EXECUTE TRANSFORM*/
	////////////////////////////////////////////////////////////////////////////////////////////
	printf("\nDEFINING RESAMPLER\n");
	LabelImageType::ConstPointer inputImage = sourceReader->GetOutput();

	LabelImageType::SpacingType spacing = destinationImage->GetSpacing();
	LabelImageType::PointType   origin  = destinationImage->GetOrigin();
	LabelImageType::DirectionType direction  = destinationImage->GetDirection();
	LabelImageType::RegionType region =  
destinationImage->GetLargestPossibleRegion();
	LabelImageType::SizeType   size =  region.GetSize();
	cout << spacing <<std::endl;
	cout << origin <<std::endl;
	cout << direction <<std::endl;
	cout << region <<std::endl;
	cout << size <<std::endl;



	if (argc > 6)
	{
		ResamplerType::Pointer resampler = ResamplerType::New();
		InterpolatorType::Pointer interpolator = InterpolatorType::New();
		resampler->SetInterpolator( interpolator );
		resampler->SetOutputSpacing( spacing );
		resampler->SetOutputDirection( direction );
		resampler->SetOutputOrigin(  origin  );
		resampler->SetSize( size );
		resampler->SetTransform( tps );

		resampler->SetOutputStartIndex(  region.GetIndex() );
		resampler->SetInput( sourceReader->GetOutput() );
		////////////////////////////////////////////////////////////////////////////////////////////

		DeformedImageWriterType::Pointer deformedImageWriter = DeformedImageWriterType::New();
		deformedImageWriter->SetInput( resampler->GetOutput() );
		deformedImageWriter->SetFileName( argv[6] );

		try
		{
		   deformedImageWriter->Update();
		}
		catch( itk::ExceptionObject & excp )
		{
		   std::cerr << "Exception thrown " << std::endl;
		   std::cerr << excp << std::endl;
		   return EXIT_FAILURE;
		}
	}

	////////////////////////////////////////////////////////////////////////////////////////////
	/*DEFINING FIELD*/
	////////////////////////////////////////////////////////////////////////////////////////////
	printf("\nDEFINING FIELD\n");
     DisplacementFieldType::Pointer field = DisplacementFieldType::New();
     field->SetRegions( region );
     field->SetOrigin( origin );
     field->SetDirection( direction );
     field->SetSpacing( spacing );
     field->Allocate();

	cout << field->GetLargestPossibleRegion() <<std::endl;
	cout << field->GetOrigin() <<std::endl;
	cout << field->GetDirection() <<std::endl;
	cout << field->GetSpacing() <<std::endl;
	cout << field->GetLargestPossibleRegion().GetSize() <<std::endl;
	////////////////////////////////////////////////////////////////////////////////////////////

	////////////////////////////////////////////////////////////////////////////////////////////
	/*APPLYING TRANSFORM TO POINTS*/
	////////////////////////////////////////////////////////////////////////////////////////////
     printf("\nAPPLYING TRANSFORM TO POINTS\n");
     typedef itk::ImageRegionIterator< DisplacementFieldType > FieldIterator;
     FieldIterator fi( field, region );
     fi.GoToBegin();
     TransformType::InputPointType  point1;
     TransformType::OutputPointType point2;
     DisplacementFieldType::IndexType index;

     FieldVectorType displacement;
     int fieldCount
     while( ! fi.IsAtEnd() )
     {
	  index = fi.GetIndex();
	  field->TransformIndexToPhysicalPoint( index, point1 );
	  point2 = tps->TransformPoint( point1 );

	  for ( int i = 0;i < ImageDimension;i++)
	  {
		displacement[i] = point2[i] - point1[i];
	  }
	  fi.Set( displacement );
	  ++fi;
     }
	////////////////////////////////////////////////////////////////////////////////////////////

	////////////////////////////////////////////////////////////////////////////////////////////
	/*GENERATING FIELD FILE*/
	////////////////////////////////////////////////////////////////////////////////////////////
     //Write computed deformation field
     printf("\nGENERATING FIELD FILE\n");
     FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
     fieldWriter->SetFileName( argv[5] );
     fieldWriter->SetInput( field );
     try
       {
       fieldWriter->Update();
       }
     catch( itk::ExceptionObject & excp )
       {
       std::cerr << "Exception thrown " << std::endl;
       std::cerr << excp << std::endl;
       return EXIT_FAILURE;
       }
     return EXIT_SUCCESS;
	////////////////////////////////////////////////////////////////////////////////////////////

}


insight-users-request at itk.org a écrit :

> Send Insight-users mailing list submissions to
> 	insight-users at itk.org
>
> To subscribe or unsubscribe via the World Wide Web, visit
> 	http://public.kitware.com/mailman/listinfo/insight-users
> or, via email, send a message with subject or body 'help' to
> 	insight-users-request at itk.org
>
> You can reach the person managing the list at
> 	insight-users-owner at itk.org
>
> When replying, please edit your Subject line so it is more specific 
> than "Re: Contents of Insight-users digest..."
>
>
> Today's Topics:
>
>    1. Re: Generating a warping field with
>       ThinPlateSplineKernelTransform (Siavash Khallaghi)
>    2. Re: 3D/2D registration with
>       itkRayCastInterpolateImageFunction (Siavash Khallaghi)
>    3. Re: [ITK] Generating a warping field with
>       ThinPlateSplineKernelTransform (Andras Lasso)
>
>
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Tue, 4 Aug 2015 14:13:16 -0700 (MST)
> From: Siavash Khallaghi <siavashk at ece.ubc.ca>
> To: insight-users at itk.org
> Subject: Re: [ITK-users] Generating a warping field with
> 	ThinPlateSplineKernelTransform
> Message-ID: <1438722796207-35972.post at n7.nabble.com>
> Content-Type: text/plain; charset=us-ascii
>
> Having your exact code (on the small chance that it differs from the 
> example), your data, and the argument list would help in debugging 
> your problem. But in the mean time, you can try the following to see 
> where the problem lies:
>
> 1. If you write out the defomation field and view it in the context of 
> your volume in a 3rd party software (e.g., Paraview), what would you 
> see? Would it make sense in terms of smoothness and direction?
>
> 2. Is the deformed volume (output of itk::ResmapleImageFilter) 
> different compared to the input volume? In other words, is the 
> resampler doing anything?
>
> 3. If you invert the deformation field, would your deformed volume be 
> the output that you expect? You can do this by swapping source and 
> target landmarks.
>
> Siavash
>
>
>
> --
> View this message in context:  
> http://itk-users.7.n7.nabble.com/ITK-users-Generating-a-warping-field-
> with-ThinPlateSplineKernelTransform-tp35971p35972.html
> Sent from the ITK - Users mailing list archive at Nabble.com.
>
>
> ------------------------------
>
> Message: 2
> Date: Tue, 4 Aug 2015 14:53:46 -0700 (MST)
> From: Siavash Khallaghi <siavashk at ece.ubc.ca>
> To: insight-users at itk.org
> Subject: Re: [ITK-users] 3D/2D registration with
> 	itkRayCastInterpolateImageFunction
> Message-ID: <1438725226436-35973.post at n7.nabble.com>
> Content-Type: text/plain; charset=us-ascii
>
> Is the metric decreasing? (i.e. similarity value) It could be that 
> everything is correct, it is just that you have an ill-posed problem. 
> If the metric is decreasing, it could be that you are converging to a local minimum.
>
> Siavash
>
>
>
> --
> View this message in context:  
> http://itk-users.7.n7.nabble.com/ITK-users-3D-2D-registration-with-itk
> RayCastInterpolateImageFunction-tp35970p35973.html
> Sent from the ITK - Users mailing list archive at Nabble.com.
>
>
> ------------------------------
>
> Message: 3
> Date: Tue, 4 Aug 2015 22:38:37 +0000
> From: Andras Lasso <lasso at queensu.ca>
> To: Siavash Khallaghi <siavashk at ece.ubc.ca>, "insight-users at itk.org"
> 	<insight-users at itk.org>
> Subject: Re: [ITK-users] [ITK] Generating a warping field with
> 	ThinPlateSplineKernelTransform
> Message-ID:
> 	
> <CY1PR0701MB1613B15AF8B413FC7C21C855D8760 at CY1PR0701MB1613.namprd07.pro
> d.outlook.com>
>
> Content-Type: text/plain; charset="us-ascii"
>
> You can also load non-linear transforms into 3D Slicer for interactive 
> visualization, application to your images, comparing 
> fixed/moving/registered images, etc. You can also invert, compose, 
> decompose transforms, apply it to marked points, segmented surfaces, 
> etc.
>
> See some examples here:
> http://www.slicer.org/slicerWiki/index.php/Documentation/Nightly/Modul
> es/Transforms
>
> Paraview is very good but not for analyzing transforms: it can only 
> load exported displacement fields (but for that you would need to be 
> sure that your region of interest is correct and you lose information 
> when you sample the field), it only supports axis-aligned images, and 
> dense vector field visualization has been broken in recent versions (I 
> haven't checked the latest version).
>
> Andras
>
> -----Original Message-----
> From: Community [mailto:community-bounces at itk.org] On Behalf Of 
> Siavash Khallaghi
> Sent: Tuesday, August 4, 2015 5:13 PM
> To: insight-users at itk.org
> Subject: Re: [ITK] [ITK-users] Generating a warping field with 
> ThinPlateSplineKernelTransform
>
> Having your exact code (on the small chance that it differs from the 
> example), your data, and the argument list would help in debugging 
> your problem. But in the mean time, you can try the following to see 
> where the problem lies:
>
> 1. If you write out the defomation field and view it in the context of 
> your volume in a 3rd party software (e.g., Paraview), what would you 
> see? Would it make sense in terms of smoothness and direction?
>
> 2. Is the deformed volume (output of itk::ResmapleImageFilter) 
> different compared to the input volume? In other words, is the 
> resampler doing anything?
>
> 3. If you invert the deformation field, would your deformed volume be 
> the output that you expect? You can do this by swapping source and 
> target landmarks.
>
> Siavash
>
>
>
> --
> View this message in context:  
> http://itk-users.7.n7.nabble.com/ITK-users-Generating-a-warping-field-
> with-ThinPlateSplineKernelTransform-tp35971p35972.html
> Sent from the ITK - Users mailing list archive at Nabble.com.
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at 
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.php
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://public.kitware.com/mailman/listinfo/insight-users
> _______________________________________________
> Community mailing list
> Community at itk.org
> http://public.kitware.com/mailman/listinfo/community
>
>
> ------------------------------
>
> Subject: Digest Footer
>
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://public.kitware.com/mailman/listinfo/insight-users
>
>
> ------------------------------
>
> End of Insight-users Digest, Vol 136, Issue 2
> *********************************************



_____________________________________
Powered by www.kitware.com

Visit other Kitware open-source projects at http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/listinfo/insight-users
_______________________________________________
Community mailing list
Community at itk.org
http://public.kitware.com/mailman/listinfo/community
_____________________________________
Powered by www.kitware.com

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/listinfo/insight-users


More information about the Community mailing list