[ITK-users] Generating a warping field with ThinPlateSplineKernelTransform
Olivier Comtois
olivier.comtois at polymtl.ca
Thu Aug 6 11:53:32 EDT 2015
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-itkRayCastInterpolateImageFunction-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.prod.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/Modules/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
> *********************************************
More information about the Insight-users
mailing list