+++ /dev/null
-/*=========================================================================
- Program: vv http://www.creatis.insa-lyon.fr/rio/vv
-
- Authors belong to:
- - University of LYON http://www.universite-lyon.fr/
- - Léon Bérard cancer center http://www.centreleonberard.fr
- - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
-
- This software is distributed WITHOUT ANY WARRANTY; without even
- the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
- PURPOSE. See the copyright notices for more information.
-
- It is distributed under dual licence
-
- - BSD See included LICENSE.txt file
- - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
-===========================================================================**/
-#ifndef clitkPointTrajectoryGenericFilter_cxx
-#define clitkPointTrajectoryGenericFilter_cxx
-
-/* =================================================
- * @file clitkPointTrajectoryGenericFilter.cxx
- * @author
- * @date
- *
- * @brief
- *
- ===================================================*/
-
-#include "clitkPointTrajectoryGenericFilter.h"
-
-
-namespace clitk
-{
-
-
- //-----------------------------------------------------------
- // Constructor
- //-----------------------------------------------------------
- PointTrajectoryGenericFilter::PointTrajectoryGenericFilter()
- {
- m_Verbose=false;
- // m_InputFileName="";
- }
-
-
- //-----------------------------------------------------------
- // Update
- //-----------------------------------------------------------
- void PointTrajectoryGenericFilter::Update()
- {
- // ImageTypes
- const unsigned int ImageDimension=4;
- const unsigned int SpaceDimension=3;
- typedef itk::Vector<double, SpaceDimension> CoefficientPixelType;
- typedef itk::Vector<float, SpaceDimension> VectorPixelType;
- typedef itk::Image<CoefficientPixelType, ImageDimension> CoefficientImageType;
- typedef itk::Image<VectorPixelType, ImageDimension> VectorFieldType;
-
-
- // -----------------------------------------------
- // Reference Point List 3D
- // -----------------------------------------------
- typedef itk::Point<double, SpaceDimension> SpacePointType;
- typedef clitk::List<SpacePointType> PointListType;
- PointListType referencePointList;
- if (m_Verbose) std::cout<<"Reference point list:"<<std::endl;
- referencePointList=PointListType(m_ArgsInfo.ref_arg, m_Verbose);
-
-
- // -----------------------------------------------
- // Transform: based on points, 4DVF, spatio-Temporal transform
- // -----------------------------------------------
- typedef itk::Transform<double, 4, 4> TransformType;
- TransformType::Pointer transform;
- switch (m_ArgsInfo.transform_arg)
- {
- // ==========================
- // List of points
- // ==========================
- case 0:
- {
- //-----------------------------
- // Input point lists
- //-----------------------------
- typedef itk::Point<double, SpaceDimension> PointType;
- typedef clitk::List<PointType> PointListType;
- typedef clitk::Lists<PointType> PointListsType;
- PointListsType inputPointLists, sortedPointLists;
-
- // Read the lists
- for (unsigned int i=0; i<m_ArgsInfo.points_given; i++)
- inputPointLists.push_back(PointListType(m_ArgsInfo.points_arg[i], m_Verbose) );
-
- // Convert/sort the lists
- sortedPointLists.resize(inputPointLists[0].size());
- for (unsigned int i=0; i<inputPointLists[0].size(); i++)
- {
- sortedPointLists[i].push_back(referencePointList[i]);
- for (unsigned int j=0; j<inputPointLists.size(); j++)
- {
- sortedPointLists[i].push_back(inputPointLists[j][i]);
- }
- }
-
- // Point list Transform
- typedef clitk::PointListTransform<double, ImageDimension,ImageDimension> PointListTransformType;
- PointListTransformType::Pointer pointListTransform=PointListTransformType::New();
- pointListTransform->SetPointLists(sortedPointLists);
-
- // Vector Interpolator
- typedef PointListTransformType::PointListImageType PointListImageType;
- typedef clitk::GenericVectorInterpolator<args_info_clitkPointTrajectory,PointListImageType, double> GenericVectorInterpolatorType;
- GenericVectorInterpolatorType::Pointer genericInterpolator=GenericVectorInterpolatorType::New();
- genericInterpolator->SetArgsInfo(m_ArgsInfo);
- typedef itk::VectorInterpolateImageFunction<PointListImageType, double> InterpolatorType;
- InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
- pointListTransform->SetInterpolator(interpolator);
- transform=pointListTransform;
-
- break;
- }
-
- // ==========================
- // 4D vector field
- // ==========================
- case 1:
- {
- // Deformation field transform
- typedef clitk::DeformationFieldTransform<double, ImageDimension,ImageDimension, SpaceDimension> DeformationFieldTransformType;
- DeformationFieldTransformType::Pointer deformationFieldTransform=DeformationFieldTransformType::New();
-
- // The deformation field
- typedef DeformationFieldTransformType::DeformationFieldType DeformationFieldType;
- typedef itk::ImageFileReader<DeformationFieldType> InputReaderType;
- InputReaderType::Pointer reader = InputReaderType::New();
- reader->SetFileName( m_ArgsInfo.input_arg);
- reader->Update();
- DeformationFieldType::Pointer input= reader->GetOutput();
- deformationFieldTransform->SetDeformationField(input);
-
- // Vector Interpolator
- typedef clitk::GenericVectorInterpolator<args_info_clitkPointTrajectory,DeformationFieldType, double> GenericVectorInterpolatorType;
- GenericVectorInterpolatorType::Pointer genericInterpolator=GenericVectorInterpolatorType::New();
- genericInterpolator->SetArgsInfo(m_ArgsInfo);
- typedef itk::VectorInterpolateImageFunction<DeformationFieldType, double> InterpolatorType;
- InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
- deformationFieldTransform->SetInterpolator(interpolator);
- transform=deformationFieldTransform;
-
- break;
- }
-
- // ==========================
- // Spatio-Temporal transform
- // ==========================
- case 2:
- {
- // S-T transform
- typedef clitk::ShapedBLUTSpatioTemporalDeformableTransform< double, ImageDimension, ImageDimension > TransformType;
- TransformType::Pointer spatioTemporalTransform = TransformType::New();
-
-
- // Spline orders: Default is cubic splines
- CoefficientImageType::RegionType::SizeType splineOrders ;
- splineOrders.Fill(3);
- if (m_ArgsInfo.order_given)
- for(unsigned int i=0; i<ImageDimension;i++)
- splineOrders[i]=m_ArgsInfo.order_arg[i];
- if (m_Verbose) std::cout<<"Setting the spline orders to "<<splineOrders<<"..."<<std::endl;
-
- // Coefficient image
- typedef itk::ImageFileReader<CoefficientImageType> InputReaderType;
- InputReaderType::Pointer reader = InputReaderType::New();
- reader->SetFileName( m_ArgsInfo.input_arg);
- reader->Update();
- CoefficientImageType::Pointer input= reader->GetOutput();
- // itk::Vector<double,3> vector;
- // vector.Fill(0.);
- // vector[2]=100;
- // input->FillBuffer(vector);
-
- // Mask
- typedef itk::ImageMaskSpatialObject< ImageDimension > MaskType;
- MaskType::Pointer spatialObjectMask=NULL;
- if (m_ArgsInfo.mask_given)
- {
- typedef itk::Image< unsigned char, ImageDimension > ImageMaskType;
- typedef itk::ImageFileReader< ImageMaskType > MaskReaderType;
- MaskReaderType::Pointer maskReader = MaskReaderType::New();
- maskReader->SetFileName(m_ArgsInfo.mask_arg);
-
- try
- {
- maskReader->Update();
- }
- catch( itk::ExceptionObject & err )
- {
- std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
- std::cerr << err << std::endl;
- return;
- }
- if (m_Verbose)std::cout <<"Mask was read..." <<std::endl;
-
- // Set the image to the spatialObject
- spatialObjectMask = MaskType::New();
- spatialObjectMask->SetImage( maskReader->GetOutput() );
- }
-
- // Samplingfactors
- CoefficientImageType::SizeType samplingFactors;
- for (unsigned int i=0; i< ImageDimension-1; i++)
- {
- samplingFactors[i]= (int) ( input->GetSpacing()[i]/ m_ArgsInfo.spacing_arg);
- if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<samplingFactors[i]<<"..."<<std::endl;
- }
- samplingFactors[ImageDimension-1]= (int) ( input->GetSpacing()[ImageDimension-1]/ m_ArgsInfo.phaseIncrement_arg);
- if (m_Verbose) std::cout<<"Setting sampling factor "<<ImageDimension-1<<" to "<<samplingFactors[ImageDimension-1]<<"..."<<std::endl;
- if( (m_ArgsInfo.shape_arg==3) |
- (m_ArgsInfo.shape_arg==4) |
- (m_ArgsInfo.shape_arg==6) |
- (m_ArgsInfo.shape_arg==8)
- ) samplingFactors[ImageDimension-1]*=2.5;
-
- // Set
- spatioTemporalTransform->SetTransformShape(m_ArgsInfo.shape_arg);
- spatioTemporalTransform->SetSplineOrders(splineOrders);
- spatioTemporalTransform->SetMask(spatialObjectMask);
- spatioTemporalTransform->SetLUTSamplingFactors(samplingFactors);
- spatioTemporalTransform->SetCoefficientImage(input);
- transform=spatioTemporalTransform;
-
- break;
- }
-
-// // ==========================
-// // Spatio-Temporal transform
-// // ==========================
-// case 3:
-// {
-// // S-T transform
-// typedef clitk::BSplineSpatioTemporalDeformableTransform< double, ImageDimension, ImageDimension > TransformType;
-// TransformType::Pointer spatioTemporalTransform = TransformType::New();
-
-
-// // Spline orders: Default is cubic splines
-// CoefficientImageType::RegionType::SizeType splineOrders ;
-// splineOrders.Fill(3);
-// if (m_ArgsInfo.order_given)
-// for(unsigned int i=0; i<ImageDimension;i++)
-// splineOrders[i]=m_ArgsInfo.order_arg[i];
-// if (m_Verbose) std::cout<<"Setting the spline orders to "<<splineOrders<<"..."<<std::endl;
-
-// // Coefficient image
-// typedef itk::ImageFileReader<CoefficientImageType> InputReaderType;
-// InputReaderType::Pointer reader = InputReaderType::New();
-// reader->SetFileName( m_ArgsInfo.input_arg);
-// reader->Update();
-// CoefficientImageType::Pointer input= reader->GetOutput();
-// // itk::Vector<double,3> vector;
-// // vector.Fill(0.);
-// // vector[2]=100;
-// // input->FillBuffer(vector);
-
-// // Mask
-// typedef itk::ImageMaskSpatialObject< ImageDimension > MaskType;
-// MaskType::Pointer spatialObjectMask=NULL;
-// if (m_ArgsInfo.mask_given)
-// {
-// typedef itk::Image< unsigned char, ImageDimension > ImageMaskType;
-// typedef itk::ImageFileReader< ImageMaskType > MaskReaderType;
-// MaskReaderType::Pointer maskReader = MaskReaderType::New();
-// maskReader->SetFileName(m_ArgsInfo.mask_arg);
-
-// try
-// {
-// maskReader->Update();
-// }
-// catch( itk::ExceptionObject & err )
-// {
-// std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
-// std::cerr << err << std::endl;
-// return;
-// }
-// if (m_Verbose)std::cout <<"Mask was read..." <<std::endl;
-
-// // Set the image to the spatialObject
-// spatialObjectMask = MaskType::New();
-// spatialObjectMask->SetImage( maskReader->GetOutput() );
-// }
-
-// // Samplingfactors
-// CoefficientImageType::SizeType samplingFactors;
-// for (unsigned int i=0; i< ImageDimension-1; i++)
-// {
-// samplingFactors[i]= (int) ( input->GetSpacing()[i]/ m_ArgsInfo.spacing_arg);
-// if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<samplingFactors[i]<<"..."<<std::endl;
-// }
-// samplingFactors[ImageDimension-1]= (int) ( input->GetSpacing()[ImageDimension-1]/ m_ArgsInfo.phaseIncrement_arg);
-// if (m_Verbose) std::cout<<"Setting sampling factor "<<ImageDimension-1<<" to "<<samplingFactors[ImageDimension-1]<<"..."<<std::endl;
-
-// // Set
-// //spatioTemporalTransform->SetTransformShape(m_ArgsInfo.shape_arg);
-// spatioTemporalTransform->SetSplineOrders(splineOrders);
-// spatioTemporalTransform->SetMask(spatialObjectMask);
-// spatioTemporalTransform->SetLUTSamplingFactors(samplingFactors);
-// spatioTemporalTransform->SetCoefficientImage(input);
-// transform=spatioTemporalTransform;
-
-// break;
-// }
- }
-
-
- // -----------------------------------------------
- // Construct Spatio-Temporal Point lists 4D
- // -----------------------------------------------
- typedef itk::Point<double, ImageDimension> SpaceTimePointType;
- typedef clitk::Lists<SpaceTimePointType> PointListsType;
- PointListsType pointLists(referencePointList.size());
- SpaceTimePointType spaceTimePoint;
- double phase;
- for (unsigned int i=0; i<referencePointList.size(); i++)
- {
- for (unsigned int d=0; d<SpaceDimension; d++)
- {
- spaceTimePoint[d]=referencePointList[i][d];
- }
- phase=0;
- while (phase<10.)
- {
- spaceTimePoint[ImageDimension-1]=phase;
- pointLists[i].push_back(spaceTimePoint);
- phase+=m_ArgsInfo.phaseIncrement_arg;
- }
- }
-
-
- // -----------------------------------------------
- // Transform Points
- // -----------------------------------------------
- typedef itk::Vector<double, ImageDimension> VectorType;
- typedef clitk::List<VectorType> VectorListType;
- typedef clitk::Lists<VectorType> VectorListsType;
- VectorListsType displacementLists(pointLists.size());
- VectorType displacement;
- for (unsigned int i=0; i<pointLists.size(); i++)
- {
- if (m_Verbose) std::cout<<"Transforming point "<<pointLists[i][0]<<"..."<<std::endl;
- for (unsigned int j=0; j<pointLists[i].size(); j++)
- {
- spaceTimePoint= transform->TransformPoint(pointLists[i][j]);
- if (m_Verbose) std::cout<<"Transformed point "<<spaceTimePoint<<"..."<<std::endl;
- displacement=spaceTimePoint-pointLists[i][j];
- displacementLists[i].push_back(displacement);
- }
- }
-
-
- // -----------------------------------------------
- // Write displacements
- // -----------------------------------------------
- std::vector<std::string> filenames;
- for (unsigned int i=0;i<displacementLists.size();i++)
- {
- std::ostringstream number_ostr;
- number_ostr << i;
- std::string number_str;
- if (i<10) number_str= "0"+number_ostr.str();
- else number_str = number_ostr.str();
- filenames.push_back(m_ArgsInfo.trajectory_arg+number_str);
- }
- displacementLists.Write(filenames, m_Verbose );
-
- }
-
-
-} //end clitk
-
-#endif //#define clitkPointTrajectoryGenericFilter_cxx
git://git.creatis.insa-lyon.fr / clitk.git / blobdiff