From: David Sarrut Date: Thu, 23 Feb 2012 15:32:44 +0000 (+0100) Subject: Merge branch 'master' of git.creatis.insa-lyon.fr:clitk X-Git-Tag: v1.3.0~52^2~58 X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=5312281fb385245fff3d64693e1549ffc27206e4;hp=268292ed7f629f715a0aa7ae9e51de6dfacbc591;p=clitk.git Merge branch 'master' of git.creatis.insa-lyon.fr:clitk add qvtk Conflicts: cmake/dependencies.cmake --- diff --git a/common/clitkImageToImageGenericFilterBase.cxx b/common/clitkImageToImageGenericFilterBase.cxx index 5f9cfac..8a713d2 100644 --- a/common/clitkImageToImageGenericFilterBase.cxx +++ b/common/clitkImageToImageGenericFilterBase.cxx @@ -284,11 +284,11 @@ void clitk::ImageToImageGenericFilterBase::SetNextOutput::Pointer clitk::ImageToImageGenericFilterBase::GetInput >(unsigned int n); -#define DEF_SetNextOutput_And_GetInput_WithCompo(Compo, Dim) \ +#define DEF_SetNextOutput_And_GetInput_WithCompo(PixelType, Compo, Dim) \ template \ - void clitk::ImageToImageGenericFilterBase::SetNextOutput, Dim> >(itk::Image,Dim>::Pointer output); \ + void clitk::ImageToImageGenericFilterBase::SetNextOutput, Dim> >(itk::Image,Dim>::Pointer output); \ template \ - itk::Image, Dim>::Pointer clitk::ImageToImageGenericFilterBase::GetInput, Dim> >(unsigned int n); + itk::Image, Dim>::Pointer clitk::ImageToImageGenericFilterBase::GetInput, Dim> >(unsigned int n); DEF_SetNextOutput_And_GetInput(char, 2); DEF_SetNextOutput_And_GetInput(unsigned char, 2); @@ -306,15 +306,25 @@ DEF_SetNextOutput_And_GetInput(int, 3); DEF_SetNextOutput_And_GetInput(float, 3); DEF_SetNextOutput_And_GetInput(double, 3); -DEF_SetNextOutput_And_GetInput_WithCompo(2, 2); -DEF_SetNextOutput_And_GetInput_WithCompo(2, 3); -DEF_SetNextOutput_And_GetInput_WithCompo(2, 4); -DEF_SetNextOutput_And_GetInput_WithCompo(3, 2); -DEF_SetNextOutput_And_GetInput_WithCompo(3, 3); -DEF_SetNextOutput_And_GetInput_WithCompo(3, 4); -DEF_SetNextOutput_And_GetInput_WithCompo(4, 2); -DEF_SetNextOutput_And_GetInput_WithCompo(4, 3); -DEF_SetNextOutput_And_GetInput_WithCompo(4, 4); +DEF_SetNextOutput_And_GetInput_WithCompo(float, 2, 2); +DEF_SetNextOutput_And_GetInput_WithCompo(float, 2, 3); +DEF_SetNextOutput_And_GetInput_WithCompo(float, 2, 4); +DEF_SetNextOutput_And_GetInput_WithCompo(float, 3, 2); +DEF_SetNextOutput_And_GetInput_WithCompo(float, 3, 3); +DEF_SetNextOutput_And_GetInput_WithCompo(float, 3, 4); +DEF_SetNextOutput_And_GetInput_WithCompo(float, 4, 2); +DEF_SetNextOutput_And_GetInput_WithCompo(float, 4, 3); +DEF_SetNextOutput_And_GetInput_WithCompo(float, 4, 4); + +DEF_SetNextOutput_And_GetInput_WithCompo(double, 2, 2); +DEF_SetNextOutput_And_GetInput_WithCompo(double, 2, 3); +DEF_SetNextOutput_And_GetInput_WithCompo(double, 2, 4); +DEF_SetNextOutput_And_GetInput_WithCompo(double, 3, 2); +DEF_SetNextOutput_And_GetInput_WithCompo(double, 3, 3); +DEF_SetNextOutput_And_GetInput_WithCompo(double, 3, 4); +DEF_SetNextOutput_And_GetInput_WithCompo(double, 4, 2); +DEF_SetNextOutput_And_GetInput_WithCompo(double, 4, 3); +DEF_SetNextOutput_And_GetInput_WithCompo(double, 4, 4); DEF_SetNextOutput_And_GetInput(char, 4); DEF_SetNextOutput_And_GetInput(unsigned char, 4); diff --git a/common/vvImageReader.txx b/common/vvImageReader.txx index e15e6d7..6a4333f 100644 --- a/common/vvImageReader.txx +++ b/common/vvImageReader.txx @@ -24,7 +24,7 @@ #include #include #include -#include +#include #include @@ -36,7 +36,12 @@ template void vvImageReader::UpdateWithDim(std::string InputPixelType) { if (mType == VECTORFIELD || mType == VECTORFIELDWITHTIME) - UpdateWithDimAndInputVectorPixelType,VImageDimension>(); + { + if (VImageDimension == 4) + UpdateWithDimAndInputVectorPixelType,VImageDimension>(); + else + UpdateWithDimAndInputVectorPixelType,VImageDimension>(); + } else if (InputPixelType == "short") UpdateWithDimAndInputPixelType(); else if (InputPixelType == "unsigned_short") @@ -235,14 +240,15 @@ void vvImageReader::UpdateWithDimAndInputVectorPixelType() } analyzeImageIO = dynamic_cast( reader->GetImageIO() ); } - + typedef itk::Image< itk::Vector, VImageDimension > VectorImageType; - typedef itk::VectorCastImageFilter CasterType; + typedef itk::FlexibleVectorCastImageFilter CasterType; typename VectorImageType::Pointer casted_input; typename CasterType::Pointer caster = CasterType::New(); caster->SetInput(input); + caster->Update(); casted_input = caster->GetOutput(); - + mImage = vvImageFromITK >(casted_input, mType == IMAGEWITHTIME || mType == VECTORFIELDWITHTIME); // For unknown analyze orientations, we set identity diff --git a/itk/clitkSetBackgroundImageFilter.h b/itk/clitkSetBackgroundImageFilter.h index 97b0bed..e5ac118 100644 --- a/itk/clitkSetBackgroundImageFilter.h +++ b/itk/clitkSetBackgroundImageFilter.h @@ -17,7 +17,7 @@ ===========================================================================**/ #ifndef __clitkSetBackgroundImageFilter_h #define __clitkSetBackgroundImageFilter_h -#include "itkBinaryFunctorImageFilter.h" +#include "itkFlexibleBinaryFunctorImageFilter.h" #include "itkNumericTraits.h" @@ -125,7 +125,7 @@ private: template class ITK_EXPORT SetBackgroundImageFilter : public - itk::BinaryFunctorImageFilter +class ITK_EXPORT FlexibleBinaryFunctorImageFilter : + public InPlaceImageFilter +{ +public: + /** Standard class typedefs. */ + typedef FlexibleBinaryFunctorImageFilter Self; + typedef InPlaceImageFilter Superclass; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(FlexibleBinaryFunctorImageFilter, InPlaceImageFilter); + + /** Some convenient typedefs. */ + typedef TFunction FunctorType; + typedef TInputImage1 Input1ImageType; + typedef typename Input1ImageType::ConstPointer Input1ImagePointer; + typedef typename Input1ImageType::RegionType Input1ImageRegionType; + typedef typename Input1ImageType::PixelType Input1ImagePixelType; + + typedef TInputImage2 Input2ImageType; + typedef typename Input2ImageType::ConstPointer Input2ImagePointer; + typedef typename Input2ImageType::RegionType Input2ImageRegionType; + typedef typename Input2ImageType::PixelType Input2ImagePixelType; + + typedef TOutputImage OutputImageType; + typedef typename OutputImageType::Pointer OutputImagePointer; + typedef typename OutputImageType::RegionType OutputImageRegionType; + typedef typename OutputImageType::PixelType OutputImagePixelType; + + /** Connect one of the operands for pixel-wise addition */ + void SetInput1( const TInputImage1 * image1); + + /** Connect one of the operands for pixel-wise addition */ + void SetInput2( const TInputImage2 * image2); + + /** Get the functor object. The functor is returned by reference. + * (Functors do not have to derive from itk::LightObject, so they do + * not necessarily have a reference count. So we cannot return a + * SmartPointer.) */ + FunctorType& GetFunctor() { return m_Functor; } + + /** Get the functor object. The functor is returned by reference. + * (Functors do not have to derive from itk::LightObject, so they do + * not necessarily have a reference count. So we cannot return a + * SmartPointer.) */ + const FunctorType& GetFunctor() const + { + return m_Functor; + } + + /** Set the functor object. This replaces the current Functor with a + * copy of the specified Functor. This allows the user to specify a + * functor that has ivars set differently than the default functor. + * This method requires an operator!=() be defined on the functor + * (or the compiler's default implementation of operator!=() being + * appropriate). */ + void SetFunctor(const FunctorType& functor) + { + if (m_Functor != functor) + { + m_Functor = functor; + this->Modified(); + } + } + + virtual void GenerateOutputInformation(); + virtual void GenerateInputRequestedRegion(); + + /** ImageDimension constants */ + itkStaticConstMacro( + InputImage1Dimension, unsigned int, TInputImage1::ImageDimension); + itkStaticConstMacro( + InputImage2Dimension, unsigned int, TInputImage2::ImageDimension); + itkStaticConstMacro( + OutputImageDimension, unsigned int, TOutputImage::ImageDimension); + +#ifdef ITK_USE_CONCEPT_CHECKING + /** Begin concept checking */ + itkConceptMacro(SameDimensionCheck1, + (Concept::SameDimension)); + itkConceptMacro(SameDimensionCheck2, + (Concept::SameDimension)); + /** End concept checking */ +#endif + +protected: + FlexibleBinaryFunctorImageFilter(); + virtual ~FlexibleBinaryFunctorImageFilter() {} + + /** FlexibleBinaryFunctorImageFilter can be implemented as a multithreaded filter. + * Therefore, this implementation provides a ThreadedGenerateData() routine + * which is called for each processing thread. The output image data is + * allocated automatically by the superclass prior to calling + * ThreadedGenerateData(). ThreadedGenerateData can only write to the + * portion of the output image specified by the parameter + * "outputRegionForThread" + * + * \sa ImageToImageFilter::ThreadedGenerateData(), + * ImageToImageFilter::GenerateData() */ +#if ITK_VERSION_MAJOR >= 4 + void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, + itk::ThreadIdType threadId ); +#else + void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, + int threadId ); +#endif + +private: + FlexibleBinaryFunctorImageFilter(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + FunctorType m_Functor; + Input2ImagePointer m_Input2; +}; + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkFlexibleBinaryFunctorImageFilter.txx" +#endif + +#endif diff --git a/itk/itkFlexibleBinaryFunctorImageFilter.txx b/itk/itkFlexibleBinaryFunctorImageFilter.txx new file mode 100644 index 0000000..f8c7a75 --- /dev/null +++ b/itk/itkFlexibleBinaryFunctorImageFilter.txx @@ -0,0 +1,163 @@ +/*========================================================================= + + Program: Insight Segmentation & Registration Toolkit + Module: $RCSfile: itkFlexibleBinaryFunctorImageFilter.txx,v $ + Language: C++ + Date: $Date: 2008-10-07 17:31:02 $ + Version: $Revision: 1.40 $ + + Copyright (c) Insight Software Consortium. All rights reserved. + See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#ifndef __itkFlexibleBinaryFunctorImageFilter_txx +#define __itkFlexibleBinaryFunctorImageFilter_txx + +#include "itkFlexibleBinaryFunctorImageFilter.h" +#include "itkImageRegionIterator.h" +#include "itkImageRegionConstIterator.h" +#include "itkProgressReporter.h" + +namespace itk +{ + +/** + * Constructor + */ +template +FlexibleBinaryFunctorImageFilter +::FlexibleBinaryFunctorImageFilter() +{ + this->SetNumberOfRequiredInputs( 1 ); + this->SetInPlace(false); +} + + +/** + * Connect one of the operands for pixel-wise addition + */ +template +void +FlexibleBinaryFunctorImageFilter +::SetInput1( const TInputImage1 * image1 ) +{ + // Process object is not const-correct so the const casting is required. + this->SetNthInput(0, const_cast( image1 )); +} + + +/** + * Connect one of the operands for pixel-wise addition + */ +template +void +FlexibleBinaryFunctorImageFilter +::SetInput2( const TInputImage2 * image2 ) +{ + // Process object is not const-correct so the const casting is required. + //this->SetNthInput(1, const_cast( image2 )); + + m_Input2 = image2; +} + +template +void +FlexibleBinaryFunctorImageFilter +::GenerateInputRequestedRegion() +{ + Superclass::GenerateInputRequestedRegion(); + + // Process object is not const-correct so the const casting is required. + // This "manual" update step is necessary because m_Input2 is not in the pipeline, since + // it's dimensions can be different from input 1. + TInputImage2* image2 = const_cast( m_Input2.GetPointer() ); + image2->Update(); +} + +template +void +FlexibleBinaryFunctorImageFilter +::GenerateOutputInformation() +{ + Superclass::GenerateOutputInformation() ; +} + +/** + * ThreadedGenerateData Performs the pixel-wise addition + */ +template +void +FlexibleBinaryFunctorImageFilter +::ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, +#if ITK_VERSION_MAJOR >= 4 + itk::ThreadIdType threadId ) +#else + int threadId) +#endif +{ + const unsigned int dim = Input1ImageType::ImageDimension; + + // We use dynamic_cast since inputs are stored as DataObjects. The + // ImageToImageFilter::GetInput(int) always returns a pointer to a + // TInputImage1 so it cannot be used for the second input. + Input1ImagePointer inputPtr1 + = dynamic_cast(ProcessObject::GetInput(0)); + Input2ImagePointer inputPtr2 = m_Input2; +/* = dynamic_cast(ProcessObject::GetInput(1));*/ + OutputImagePointer outputPtr = this->GetOutput(0); + + typename Input1ImageType::RegionType region2 = inputPtr2->GetLargestPossibleRegion(); + + typename Input1ImageType::IndexType index1; + typename Input2ImageType::IndexType index2; + typename Input1ImageType::PointType point1; + typename Input2ImageType::PointType point2; + + ImageRegionConstIterator inputIt1(inputPtr1, outputRegionForThread); + ImageRegionConstIterator inputIt2(inputPtr2, outputRegionForThread); + + ImageRegionIterator outputIt(outputPtr, outputRegionForThread); + + ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels()); + + inputIt1.GoToBegin(); + index1 = inputIt1.GetIndex(); + inputPtr1->TransformIndexToPhysicalPoint(index1, point1); + for (unsigned int i = 0; i < dim; i++) + point2[i] = point1[i]; + inputPtr2->TransformPhysicalPointToIndex(point2, index2); + inputIt2.SetIndex(index2); + outputIt.GoToBegin(); + + while( !inputIt1.IsAtEnd() ) { + + if (region2.IsInside(index2)) + outputIt.Set( m_Functor( inputIt1.Get(), inputIt2.Get() ) ); + else + outputIt.Set(inputIt1.Get()); + + ++inputIt1; + index1 = inputIt1.GetIndex(); + inputPtr1->TransformIndexToPhysicalPoint(index1, point1); + for (unsigned int i = 0; i < dim; i++) + point2[i] = point1[i]; + inputPtr2->TransformPhysicalPointToIndex(point2, index2); + inputIt2.SetIndex(index2); + ++outputIt; + + progress.CompletedPixel(); // potential exception thrown here + } +} + +} // end namespace itk + +#endif diff --git a/itk/itkFlexibleVectorCastImageFilter.h b/itk/itkFlexibleVectorCastImageFilter.h new file mode 100644 index 0000000..1e73a92 --- /dev/null +++ b/itk/itkFlexibleVectorCastImageFilter.h @@ -0,0 +1,113 @@ +/*========================================================================= + + Program: Insight Segmentation & Registration Toolkit + Module: $RCSfile: itkFlexibleVectorCastImageFilter.h,v $ + Language: C++ + Date: $Date: 2008-10-17 16:30:53 $ + Version: $Revision: 1.16 $ + + Copyright (c) Insight Software Consortium. All rights reserved. + See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#ifndef __itkFlexibleVectorCastImageFilter_h +#define __itkFlexibleVectorCastImageFilter_h + +#include "itkUnaryFunctorImageFilter.h" +#include "itkNumericTraitsFixedArrayPixel.h" + +namespace itk +{ + +/** \class FlexibleVectorCastImageFilter + * + * \brief Casts input vector pixels to output vector pixel type. + * + * This filter is templated over the input image type and + * output image type. + * + * The filter expect both images to have the same number of dimensions, + * and that both the input and output have itk::Vector pixel types + * of the same VectorDimension. + * + * \sa Vector + * + * \ingroup IntensityImageFilters Multithreaded + */ +namespace Functor { + +template< class TInput, class TOutput> +class FlexibleVectorCast +{ +public: + FlexibleVectorCast() {} + ~FlexibleVectorCast() {} + bool operator!=( const FlexibleVectorCast & ) const + { + return false; + } + bool operator==( const FlexibleVectorCast & other ) const + { + return !(*this != other); + } + inline TOutput operator()( const TInput & A ) const + { + typedef typename TOutput::ValueType OutputValueType; + + TOutput value; + unsigned int k; + for( k = 0; k < TOutput::Dimension && k < TInput::Dimension; k++ ) + { + value[k] = static_cast( A[k] ); + } + + for (; k < TOutput::Dimension; k++) + value[k] = 0; + + return value; + } +}; +} + +template > +class ITK_EXPORT FlexibleVectorCastImageFilter : + public +UnaryFunctorImageFilter > +{ +public: + /** Standard class typedefs. */ + typedef FlexibleVectorCastImageFilter Self; + typedef UnaryFunctorImageFilter< + TInputImage,TOutputImage, + Functor::FlexibleVectorCast< typename TInputImage::PixelType, + typename TOutputImage::PixelType> > Superclass; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(FlexibleVectorCastImageFilter, + UnaryFunctorImageFilter); + +protected: + FlexibleVectorCastImageFilter() {} + virtual ~FlexibleVectorCastImageFilter() {} + +private: + FlexibleVectorCastImageFilter(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + +}; + +} // end namespace itk + + +#endif diff --git a/scripts/create_midP-2.0.sh b/scripts/create_midP-2.0.sh index 3b943d3..7af8b09 100755 --- a/scripts/create_midP-2.0.sh +++ b/scripts/create_midP-2.0.sh @@ -89,19 +89,22 @@ registration() # combine in and out results out_result=$output_dir/result_${ref_phase_nb}_$phase_nb.mhd - clitkCombineImage -i $result_in -j $result_out -m $motion_mask -o $out_result + # clitkCombineImage -i $result_in -j $result_out -m $motion_mask -o $out_result + combine_image $result_in $result_out $out_result $motion_mask abort_on_error registration $? clean_up_registration # combine in and out vf vf_result=$vf_dir/vf_${ref_phase_nb}_$phase_nb.mhd - clitkCombineImage -i $vf_in -j $vf_out -m $motion_mask -o $vf_result + #clitkCombineImage -i $vf_in -j $vf_out -m $motion_mask -o $vf_result + combine_image $vf_in $vf_out $vf_result $motion_mask abort_on_error registration $? clean_up_registration clitkZeroVF -i $vf_in -o vf_zero.mhd abort_on_error registration $? clean_up_registration patient_mask=$mask_dir/patient_mask_$phase_nb.mhd - clitkCombineImage -i $vf_result -j vf_zero.mhd -m $patient_mask -o $vf_result + #clitkCombineImage -i $vf_result -j vf_zero.mhd -m $patient_mask -o $vf_result + combine_image $vf_result vf_zero.mhd $vf_result $patient_mask abort_on_error registration $? clean_up_registration rm vf_zero.* @@ -123,7 +126,8 @@ registration() reference_in=$mask_dir/${banded}inside_${ref_phase_nb}.mhd reference_out=$mask_dir/${banded}outside_$ref_phase_nb.mhd out_result=$output_dir/result_${ref_phase_nb}_${ref_phase_nb}.mhd - clitkCombineImage -i $reference_in -j $reference_out -m $motion_mask -o $out_result + #clitkCombineImage -i $reference_in -j $reference_out -m $motion_mask -o $out_result + combine_image $reference_in $reference_out $out_result $motion_mask abort_on_error registration $? clean_up_registration # create 4D vf @@ -164,7 +168,8 @@ midp() vf_midp=$midp_dir/vf_midp_$phase_nb.mhd midp=$midp_dir/midp_$phase_nb.mhd # average the vf's from reference phase to phase - clitkAverageTemporalDimension -i $vf_dir/vf_${ref_phase_nb}_4D.mhd -o $vf_midp + #clitkAverageTemporalDimension -i $vf_dir/vf_${ref_phase_nb}_4D.mhd -o $vf_midp + average_temporal_dimension $vf_dir/vf_${ref_phase_nb}_4D.mhd $vf_midp abort_on_error midp $? clean_up_midp # invert the vf (why?) @@ -207,7 +212,8 @@ midp() create_mhd_4D_pattern.sh $midp_dir/vf_midp_ echo "Calculating midp_avg.mhd..." - clitkAverageTemporalDimension -i $midp_dir/midp_4D.mhd -o $midp_dir/midp_avg.mhd + #clitkAverageTemporalDimension -i $midp_dir/midp_4D.mhd -o $midp_dir/midp_avg.mhd + average_temporal_dimension $midp_dir/midp_4D.mhd $midp_dir/midp_avg.mhd abort_on_error midp $? clean_up_midp echo "Calculating midp_med.mhd..." @@ -245,10 +251,12 @@ midp_in_out() coeff_midp_out=$midp_dir/coeff_outside_midp_$phase_nb.mhd # average the vf's from reference phase to phase # clitkAverageTemporalDimension -i $vf_dir/vf_inside_${ref_phase_nb}_4D.mhd -o $vf_midp_in - clitkAverageTemporalDimension -i $vf_dir/coeff_inside_${ref_phase_nb}_4D_0.mhd -o $coeff_midp_in + #clitkAverageTemporalDimension -i $vf_dir/coeff_inside_${ref_phase_nb}_4D_0.mhd -o $coeff_midp_in + average_temporal_dimension $vf_dir/coeff_inside_${ref_phase_nb}_4D_0.mhd $coeff_midp_in abort_on_error midp $? clean_up_midp # clitkAverageTemporalDimension -i $vf_dir/vf_outside_${ref_phase_nb}_4D.mhd -o $vf_midp_out - clitkAverageTemporalDimension -i $vf_dir/coeff_outside_${ref_phase_nb}_4D_0.mhd -o $coeff_midp_out + #clitkAverageTemporalDimension -i $vf_dir/coeff_outside_${ref_phase_nb}_4D_0.mhd -o $coeff_midp_out + average_temporal_dimension $vf_dir/coeff_outside_${ref_phase_nb}_4D_0.mhd $coeff_midp_out abort_on_error midp $? clean_up_midp # invert the vf @@ -263,9 +271,11 @@ midp_in_out() ref_vf_midp_in=$vf_midp_in ref_vf_midp_out=$vf_midp_out vf_midp=$midp_dir/vf_midp_$phase_nb.mhd - clitkCombineImage -i $vf_midp_in -j $vf_midp_out -o $vf_midp -m $mask_dir/mm_$phase_nb.mhd + #clitkCombineImage -i $vf_midp_in -j $vf_midp_out -o $vf_midp -m $mask_dir/mm_$phase_nb.mhd + combine_image $vf_midp_in $vf_midp_out $vf_midp $mask_dir/mm_$phase_nb.mhd clitkZeroVF -i $vf_midp -o vf_zero.mhd - clitkCombineImage -i $vf_midp -j vf_zero.mhd -o $vf_midp -m $mask_dir/patient_mask_$phase_nb.mhd + #clitkCombineImage -i $vf_midp -j vf_zero.mhd -o $vf_midp -m $mask_dir/patient_mask_$phase_nb.mhd + combine_image $vf_midp vf_zero.mhd $vf_midp $mask_dir/patient_mask_$phase_nb.mhd # create the midp by warping the reference phase with the reference vf midp=$midp_dir/midp_$phase_nb.mhd @@ -293,9 +303,11 @@ midp_in_out() # combine in and out VF's vf_midp=$midp_dir/vf_midp_$phase_nb.mhd - clitkCombineImage -i $vf_midp_in -j $vf_midp_out -o $vf_midp -m $mask_dir/mm_$phase_nb.mhd + #clitkCombineImage -i $vf_midp_in -j $vf_midp_out -o $vf_midp -m $mask_dir/mm_$phase_nb.mhd + combine_image $vf_midp_in $vf_midp_out $vf_midp $mask_dir/mm_$phase_nb.mhd clitkZeroVF -i $vf_midp -o vf_zero.mhd - clitkCombineImage -i $vf_midp -j vf_zero.mhd -o $vf_midp -m $mask_dir/patient_mask_$phase_nb.mhd + #clitkCombineImage -i $vf_midp -j vf_zero.mhd -o $vf_midp -m $mask_dir/patient_mask_$phase_nb.mhd + combine_image $vf_midp vf_zero.mhd $vf_midp $mask_dir/patient_mask_$phase_nb.mhd midp=$midp_dir/midp_$phase_nb.mhd clitkWarpImage -i $phase_file -o $midp --vf=$vf_midp -s 1 @@ -315,7 +327,8 @@ midp_in_out() create_mhd_4D_pattern.sh $midp_dir/vf_outside_midp_ echo "Calculating midp_avg.mhd..." - clitkAverageTemporalDimension -i $midp_dir/midp_4D.mhd -o $midp_dir/midp_avg.mhd + #clitkAverageTemporalDimension -i $midp_dir/midp_4D.mhd -o $midp_dir/midp_avg.mhd + average_temporal_dimension $midp_dir/midp_4D.mhd $midp_dir/midp_avg.mhd abort_on_error midp $? clean_up_midp echo "Calculating midp_med.mhd..." diff --git a/scripts/midp_common.sh b/scripts/midp_common.sh index aa500ae..16535a1 100755 --- a/scripts/midp_common.sh +++ b/scripts/midp_common.sh @@ -173,3 +173,37 @@ extract_4d_phases_ref() fi done } + +# +# replacement for clitkCombineImage +combine_image() +{ +# eg: -i $result_in -j $result_out -o $out_result -m $motion_mask + + clitkSetBackground -i $1 -o temp1.mhd -m $4 + clitkSetBackground -i $2 -o temp2.mhd -m $4 --fg + + clitkImageArithm -i temp1.mhd -j temp2.mhd -o $3 + rm temp?.* +} + +# +# replacement for clitkAverageTemporalDimension +average_temporal_dimension() +{ + # eg: -i $midp_dir/midp_4D.mhd -o $midp_dir/midp_avg.mhd + + local tot=tot.mhd + + local dir=`dirname $1` + local first=`grep raw $1 | sed 's/raw/mhd/g' | head -n 1` + clitkImageArithm -i $dir/$first -o $tot -t 1 -s 0 + + local nbphases=`grep raw $1 | sed 's/raw/mhd/g' | wc -l` + for i in $(grep raw $1 | sed 's/raw/mhd/g'); do + clitkImageArithm -i $dir/$i -j $tot -o $tot + done + + clitkImageArithm -i $tot -o $2 -t 11 -s $nbphases + rm tot.* +} diff --git a/tools/clitkImageArithm.ggo b/tools/clitkImageArithm.ggo index ca7e916..fd01863 100644 --- a/tools/clitkImageArithm.ggo +++ b/tools/clitkImageArithm.ggo @@ -13,7 +13,7 @@ option "input2" j "Input second image filename" string no option "output" o "Output image filename" string yes option "scalar" s "Scalar value" double no -option "operation" t "Type of operation : \n With another image : 0=add, 1=multiply, 2=divide,\n 3=max, 4=min, 5=absdiff, 6=squareddiff, 7=difference, 8=relativ diff\n For 'scalar' : 0=add, 1=multiply, 2=inverse,\n 3=max, 4=min 5=absval 6=squareval\n 7=log 8=exp 9=sqrt 10=EPID 11=divide 12=normalize (divide by max)" int default="0" no +option "operation" t "Type of operation : \n With another image : 0=add*, 1=multiply, 2=divide,\n 3=max, 4=min, 5=absdiff, 6=squareddiff, 7=difference*, 8=relativ diff\n; For 'scalar' : 0=add*, 1=multiply*, 2=inverse,\n 3=max, 4=min 5=absval 6=squareval\n 7=log 8=exp 9=sqrt 10=EPID 11=divide* 12=normalize (divide by max); \n* operations supported with vector fields as inputs." int default="0" no option "pixelValue" - "Default value for NaN/Inf" double default="0.0" no option "setFloatOutput" f "Set output to float pixel type" flag off diff --git a/tools/clitkImageArithmGenericFilter.cxx b/tools/clitkImageArithmGenericFilter.cxx index 30b570a..f0d55c2 100644 --- a/tools/clitkImageArithmGenericFilter.cxx +++ b/tools/clitkImageArithmGenericFilter.cxx @@ -19,12 +19,676 @@ #define CLITKIMAGEARITHMGENERICFILTER_CXX #include "clitkImageArithmGenericFilter.h" -#include "clitkImageArithm_ggo.h" namespace clitk { // Specialisation +// template<> +// class ImageArithmGenericFilter; + +//////////////////////////////////////////////////////////////////// + // specializations for itk::Vector, 3u + template<> + template<> + void ImageArithmGenericFilter::UpdateWithInputImageType< itk::Image< itk::Vector, 3u > >() + { + typedef itk::Image< itk::Vector, 3u > ImageType; + + // Read input1 + ImageType::Pointer input1 = this->GetInput(0); + + // Set input image iterator + typedef itk::ImageRegionIterator IteratorType; + IteratorType it(input1, input1->GetLargestPossibleRegion()); + + // typedef input2 + ImageType::Pointer input2 = NULL; + IteratorType it2; + + /* + // Special case for normalisation + if (mTypeOfOperation == 12) { + typedef itk::MinimumMaximumImageCalculator MinMaxFilterType; + typename MinMaxFilterType::Pointer ff = MinMaxFilterType::New(); + ff->SetImage(input1); + ff->ComputeMaximum(); + mScalar = ff->GetMaximum(); + mTypeOfOperation = 11; // divide + } + */ + + if (mIsOperationUseASecondImage) { + // Read input2 + input2 = this->GetInput(1); + // Set input image iterator + it2 = IteratorType(input2, input2->GetLargestPossibleRegion()); + // Check dimension + if (!clitk::HaveSameSize(input1, input2)) { + itkExceptionMacro(<< "The images (input and input2) must have the same size"); + } + if(!clitk::HaveSameSpacing(input1, input2)) { + itkWarningMacro(<< "The images (input and input2) do not have the same spacing. " + << "Using first input's information."); + } + } + + /* + // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite + if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) { + // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is " + // << typeid(PixelType).name() + // << std::endl; + mOverwriteInputImage = false; + } + */ + + // ---------------- Overwrite input Image --------------------- + if (mOverwriteInputImage) { + // Set output iterator (to input1) + IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion()); + if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito); + else ComputeImage(it, ito); + this->SetNextOutput(input1); + } + // ---------------- Create new output Image --------------------- + else { + /*if (mOutputIsFloat) { + // Create output image + typedef itk::Image OutputImageType; + typename OutputImageType::Pointer output = OutputImageType::New(); + output->SetRegions(input1->GetLargestPossibleRegion()); + output->SetOrigin(input1->GetOrigin()); + output->SetSpacing(input1->GetSpacing()); + output->Allocate(); + // Set output iterator + typedef itk::ImageRegionIterator IteratorOutputType; + IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion()); + if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito); + else ComputeImage(it, ito); + this->template SetNextOutput(output); + } else*/ { + // Create output image + typedef ImageType OutputImageType; + OutputImageType::Pointer output = OutputImageType::New(); + output->SetRegions(input1->GetLargestPossibleRegion()); + output->SetOrigin(input1->GetOrigin()); + output->SetSpacing(input1->GetSpacing()); + output->Allocate(); + // Set output iterator + typedef itk::ImageRegionIterator IteratorOutputType; + IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion()); + if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito); + else ComputeImage(it, ito); + this->SetNextOutput(output); + } + } + } + + template<> + template<> + void ImageArithmGenericFilter::ComputeImage< + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > > + (itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > ito) + { + typedef itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > Iter2; + + ito.GoToBegin(); + it.GoToBegin(); + + typedef Iter2::PixelType PixelType; + + PixelType scalar_vector; + scalar_vector.Fill(mScalar); + + // Perform operation + switch (mTypeOfOperation) { + case 0: // Addition + while (!it.IsAtEnd()) { + ito.Set(it.Get() + scalar_vector); + ++it; + ++ito; + } + break; + case 1: // Multiply + while (!it.IsAtEnd()) { + ito.Set(it.Get() * mScalar); + ++it; + ++ito; + } + break; + /* + case 2: // Inverse + while (!it.IsAtEnd()) { + if (it.Get() != 0) + ito.Set(mScalar / it.Get())); + else ito.Set(mDefaultPixelValue); + ++it; + ++ito; + } + break; + case 3: // Max + while (!it.IsAtEnd()) { + if (it.Get() < mScalar) ito.Set(PixelTypeDownCast(mScalar)); + else ito.Set(PixelTypeDownCast(it.Get())); + ++it; + ++ito; + } + break; + case 4: // Min + while (!it.IsAtEnd()) { + if (it.Get() > mScalar) ito.Set(PixelTypeDownCast(mScalar)); + else ito.Set(PixelTypeDownCast(it.Get())); + ++it; + ++ito; + } + break; + case 5: // Absolute value + while (!it.IsAtEnd()) { + if (it.Get() <= 0) ito.Set(PixelTypeDownCast(-it.Get())); + // <= zero to avoid warning for unsigned types + else ito.Set(PixelTypeDownCast(it.Get())); + ++it; + ++ito; + } + break; + case 6: // Squared value + while (!it.IsAtEnd()) { + ito.Set(PixelTypeDownCast((double)it.Get()*(double)it.Get())); + ++it; + ++ito; + } + break; + case 7: // Log + while (!it.IsAtEnd()) { + if (it.Get() > 0) + ito.Set(PixelTypeDownCast(log((double)it.Get()))); + else ito.Set(mDefaultPixelValue); + ++it; + ++ito; + } + break; + case 8: // exp + while (!it.IsAtEnd()) { + ito.Set(PixelTypeDownCast(exp((double)it.Get()))); + ++it; + ++ito; + } + break; + case 9: // sqrt + while (!it.IsAtEnd()) { + if (it.Get() > 0) + ito.Set(PixelTypeDownCast(sqrt((double)it.Get()))); + else { + if (it.Get() ==0) ito.Set(0); + else ito.Set(mDefaultPixelValue); + } + ++it; + ++ito; + } + break; + case 10: // exp + while (!it.IsAtEnd()) { + ito.Set(PixelTypeDownCast((0x10000 - (double)it.Get())/mScalar)); + ++it; + ++ito; + } + break; + */ + case 11: // divide + while (!it.IsAtEnd()) { + ito.Set(it.Get() / mScalar); + ++it; + ++ito; + } + break; + default: // error ? + std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl; + exit(-1); + } + + } + + template<> + template<> + void ImageArithmGenericFilter::ComputeImage< + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > + > + (itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it1, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it2, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > ito) + { + typedef itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > Iter1; + typedef itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > Iter2; + typedef itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > Iter3; + + it1.GoToBegin(); + it2.GoToBegin(); + ito.GoToBegin(); + typedef Iter3::PixelType PixelType; + + switch (mTypeOfOperation) { + case 0: // Addition + while (!ito.IsAtEnd()) { + ito.Set(it1.Get() + it2.Get()); + ++it1; + ++it2; + ++ito; + } + break; + /* + case 1: // Multiply + while (!ito.IsAtEnd()) { + ito.Set(it1.Get() * it2.Get()) ); + ++it1; + ++it2; + ++ito; + } + break; + case 2: // Divide + while (!ito.IsAtEnd()) { + if (it1.Get() != 0) + ito.Set(it1.Get() / it2.Get())); + else ito.Set(mDefaultPixelValue); + ++it1; + ++it2; + ++ito; + } + break; + case 3: // Max + while (!ito.IsAtEnd()) { + if (it1.Get() < it2.Get()) ito.Set(it2.Get()); + ++it1; + ++it2; + ++ito; + } + break; + case 4: // Min + while (!ito.IsAtEnd()) { + if (it1.Get() > it2.Get()) ito.Set(it2.Get()); + ++it1; + ++it2; + ++ito; + } + break; + */ + case 5: // Absolute difference + while (!ito.IsAtEnd()) { + ito.Set(it2.Get()-it1.Get()); + ++it1; + ++it2; + ++ito; + } + break; + /* + case 6: // Squared differences + while (!ito.IsAtEnd()) { + ito.Set(pow(it1.Get()-it2.Get(),2))); + ++it1; + ++it2; + ++ito; + } + break; + */ + case 7: // Difference + while (!ito.IsAtEnd()) { + ito.Set(it1.Get()-it2.Get()); + ++it1; + ++it2; + ++ito; + } + break; + /* + case 8: // Relative Difference + while (!ito.IsAtEnd()) { + if (it1.Get() != 0) ito.Set(PixelTypeDownCast(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get()); + else ito.Set(0.0); + ++it1; + ++it2; + ++ito; + } + break; + */ + default: // error ? + std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl; + exit(-1); + } + } + +//////////////////////////////////////////////////////////////////// + // specializations for itk::Vector, 3u + template<> + template<> + void ImageArithmGenericFilter::UpdateWithInputImageType< itk::Image< itk::Vector, 3u > >() + { + typedef itk::Image< itk::Vector, 3u > ImageType; + + // Read input1 + ImageType::Pointer input1 = this->GetInput(0); + + // Set input image iterator + typedef itk::ImageRegionIterator IteratorType; + IteratorType it(input1, input1->GetLargestPossibleRegion()); + + // typedef input2 + ImageType::Pointer input2 = NULL; + IteratorType it2; + + /* + // Special case for normalisation + if (mTypeOfOperation == 12) { + typedef itk::MinimumMaximumImageCalculator MinMaxFilterType; + typename MinMaxFilterType::Pointer ff = MinMaxFilterType::New(); + ff->SetImage(input1); + ff->ComputeMaximum(); + mScalar = ff->GetMaximum(); + mTypeOfOperation = 11; // divide + } + */ + + if (mIsOperationUseASecondImage) { + // Read input2 + input2 = this->GetInput(1); + // Set input image iterator + it2 = IteratorType(input2, input2->GetLargestPossibleRegion()); + // Check dimension + if (!clitk::HaveSameSize(input1, input2)) { + itkExceptionMacro(<< "The images (input and input2) must have the same size"); + } + if(!clitk::HaveSameSpacing(input1, input2)) { + itkWarningMacro(<< "The images (input and input2) do not have the same spacing. " + << "Using first input's information."); + } + } + + /* + // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite + if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) { + // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is " + // << typeid(PixelType).name() + // << std::endl; + mOverwriteInputImage = false; + } + */ + + // ---------------- Overwrite input Image --------------------- + if (mOverwriteInputImage) { + // Set output iterator (to input1) + IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion()); + if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito); + else ComputeImage(it, ito); + this->SetNextOutput(input1); + } + // ---------------- Create new output Image --------------------- + else { + /*if (mOutputIsFloat) { + // Create output image + typedef itk::Image OutputImageType; + typename OutputImageType::Pointer output = OutputImageType::New(); + output->SetRegions(input1->GetLargestPossibleRegion()); + output->SetOrigin(input1->GetOrigin()); + output->SetSpacing(input1->GetSpacing()); + output->Allocate(); + // Set output iterator + typedef itk::ImageRegionIterator IteratorOutputType; + IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion()); + if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito); + else ComputeImage(it, ito); + this->template SetNextOutput(output); + } else*/ { + // Create output image + typedef ImageType OutputImageType; + OutputImageType::Pointer output = OutputImageType::New(); + output->SetRegions(input1->GetLargestPossibleRegion()); + output->SetOrigin(input1->GetOrigin()); + output->SetSpacing(input1->GetSpacing()); + output->Allocate(); + // Set output iterator + typedef itk::ImageRegionIterator IteratorOutputType; + IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion()); + if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito); + else ComputeImage(it, ito); + this->SetNextOutput(output); + } + } + } + template<> - class ImageArithmGenericFilter; + template<> + void ImageArithmGenericFilter::ComputeImage< + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > > + (itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > ito) + { + typedef itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > Iter2; + + ito.GoToBegin(); + it.GoToBegin(); + + typedef Iter2::PixelType PixelType; + + PixelType scalar_vector; + scalar_vector.Fill(mScalar); + + // Perform operation + switch (mTypeOfOperation) { + case 0: // Addition + while (!it.IsAtEnd()) { + ito.Set(it.Get() + scalar_vector); + ++it; + ++ito; + } + break; + case 1: // Multiply + while (!it.IsAtEnd()) { + ito.Set(it.Get() * mScalar); + ++it; + ++ito; + } + break; + /* + case 2: // Inverse + while (!it.IsAtEnd()) { + if (it.Get() != 0) + ito.Set(mScalar / it.Get())); + else ito.Set(mDefaultPixelValue); + ++it; + ++ito; + } + break; + case 3: // Max + while (!it.IsAtEnd()) { + if (it.Get() < mScalar) ito.Set(PixelTypeDownCast(mScalar)); + else ito.Set(PixelTypeDownCast(it.Get())); + ++it; + ++ito; + } + break; + case 4: // Min + while (!it.IsAtEnd()) { + if (it.Get() > mScalar) ito.Set(PixelTypeDownCast(mScalar)); + else ito.Set(PixelTypeDownCast(it.Get())); + ++it; + ++ito; + } + break; + case 5: // Absolute value + while (!it.IsAtEnd()) { + if (it.Get() <= 0) ito.Set(PixelTypeDownCast(-it.Get())); + // <= zero to avoid warning for unsigned types + else ito.Set(PixelTypeDownCast(it.Get())); + ++it; + ++ito; + } + break; + case 6: // Squared value + while (!it.IsAtEnd()) { + ito.Set(PixelTypeDownCast((double)it.Get()*(double)it.Get())); + ++it; + ++ito; + } + break; + case 7: // Log + while (!it.IsAtEnd()) { + if (it.Get() > 0) + ito.Set(PixelTypeDownCast(log((double)it.Get()))); + else ito.Set(mDefaultPixelValue); + ++it; + ++ito; + } + break; + case 8: // exp + while (!it.IsAtEnd()) { + ito.Set(PixelTypeDownCast(exp((double)it.Get()))); + ++it; + ++ito; + } + break; + case 9: // sqrt + while (!it.IsAtEnd()) { + if (it.Get() > 0) + ito.Set(PixelTypeDownCast(sqrt((double)it.Get()))); + else { + if (it.Get() ==0) ito.Set(0); + else ito.Set(mDefaultPixelValue); + } + ++it; + ++ito; + } + break; + case 10: // exp + while (!it.IsAtEnd()) { + ito.Set(PixelTypeDownCast((0x10000 - (double)it.Get())/mScalar)); + ++it; + ++ito; + } + break; + */ + case 11: // divide + while (!it.IsAtEnd()) { + ito.Set(it.Get() / mScalar); + ++it; + ++ito; + } + break; + default: // error ? + std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl; + exit(-1); + } + + } + + template<> + template<> + void ImageArithmGenericFilter::ComputeImage< + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > + > + (itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it1, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it2, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > ito) + { + typedef itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > Iter1; + typedef itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > Iter2; + typedef itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > Iter3; + + it1.GoToBegin(); + it2.GoToBegin(); + ito.GoToBegin(); + typedef Iter3::PixelType PixelType; + + switch (mTypeOfOperation) { + case 0: // Addition + while (!ito.IsAtEnd()) { + ito.Set(it1.Get() + it2.Get()); + ++it1; + ++it2; + ++ito; + } + break; + /* + case 1: // Multiply + while (!ito.IsAtEnd()) { + ito.Set(it1.Get() * it2.Get()) ); + ++it1; + ++it2; + ++ito; + } + break; + case 2: // Divide + while (!ito.IsAtEnd()) { + if (it1.Get() != 0) + ito.Set(it1.Get() / it2.Get())); + else ito.Set(mDefaultPixelValue); + ++it1; + ++it2; + ++ito; + } + break; + case 3: // Max + while (!ito.IsAtEnd()) { + if (it1.Get() < it2.Get()) ito.Set(it2.Get()); + ++it1; + ++it2; + ++ito; + } + break; + case 4: // Min + while (!ito.IsAtEnd()) { + if (it1.Get() > it2.Get()) ito.Set(it2.Get()); + ++it1; + ++it2; + ++ito; + } + break; + */ + case 5: // Absolute difference + while (!ito.IsAtEnd()) { + ito.Set(it2.Get()-it1.Get()); + ++it1; + ++it2; + ++ito; + } + break; + /* + case 6: // Squared differences + while (!ito.IsAtEnd()) { + ito.Set(pow(it1.Get()-it2.Get(),2))); + ++it1; + ++it2; + ++ito; + } + break; + */ + case 7: // Difference + while (!ito.IsAtEnd()) { + ito.Set(it1.Get()-it2.Get()); + ++it1; + ++it2; + ++ito; + } + break; + /* + case 8: // Relative Difference + while (!ito.IsAtEnd()) { + if (it1.Get() != 0) ito.Set(PixelTypeDownCast(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get()); + else ito.Set(0.0); + ++it1; + ++it2; + ++ito; + } + break; + */ + default: // error ? + std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl; + exit(-1); + } + } + } #endif //define CLITKIMAGEARITHMGENERICFILTER_CXX diff --git a/tools/clitkImageArithmGenericFilter.h b/tools/clitkImageArithmGenericFilter.h index 46ae08c..39bfeb9 100644 --- a/tools/clitkImageArithmGenericFilter.h +++ b/tools/clitkImageArithmGenericFilter.h @@ -29,6 +29,7 @@ // clitk include #include "clitkCommon.h" #include "clitkImageToImageGenericFilter.h" +#include "clitkImageArithm_ggo.h" // itk include #include "itkImage.h" @@ -95,9 +96,54 @@ namespace clitk { //-------------------------------------------------------------------- }; // end class ImageArithmGenericFilter + + // specializations for itk::Vector, 3u + template<> template<> + void ImageArithmGenericFilter::UpdateWithInputImageType< itk::Image< itk::Vector, 3u > >(); + + template<> template<> + void ImageArithmGenericFilter::ComputeImage< + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > + > + (itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > ito); + + template<> template<> + void ImageArithmGenericFilter::ComputeImage< + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > + > + (itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it1, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it2, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > ito); + + // specializations for itk::Vector, 3u + template<> template<> + void ImageArithmGenericFilter::UpdateWithInputImageType< itk::Image< itk::Vector, 3u > >(); + + template<> template<> + void ImageArithmGenericFilter::ComputeImage< + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > + > + (itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > ito); + + template<> template<> + void ImageArithmGenericFilter::ComputeImage< + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > + > + (itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it1, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it2, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > ito); } // end namespace //-------------------------------------------------------------------- + #ifndef ITK_MANUAL_INSTANTIATION #include "clitkImageArithmGenericFilter.txx" #endif diff --git a/tools/clitkImageArithmGenericFilter.txx b/tools/clitkImageArithmGenericFilter.txx index 24ebbbd..1afb77f 100644 --- a/tools/clitkImageArithmGenericFilter.txx +++ b/tools/clitkImageArithmGenericFilter.txx @@ -45,6 +45,8 @@ template void ImageArithmGenericFilter::InitializeImageType() { ADD_DEFAULT_IMAGE_TYPES(Dim); + ADD_VEC_IMAGE_TYPE(3u,3u,float); + ADD_VEC_IMAGE_TYPE(3u,3u,double); } //-------------------------------------------------------------------- @@ -398,6 +400,8 @@ void clitk::ImageArithmGenericFilter::ComputeImage(Iter1 it, Ite } //-------------------------------------------------------------------- + + } // end namespace #endif //#define CLITKIMAGEARITHMGENERICFILTER_TXX diff --git a/tools/clitkSetBackgroundGenericFilter.cxx b/tools/clitkSetBackgroundGenericFilter.cxx index 038fcca..2cb3c00 100644 --- a/tools/clitkSetBackgroundGenericFilter.cxx +++ b/tools/clitkSetBackgroundGenericFilter.cxx @@ -50,19 +50,50 @@ SetBackgroundGenericFilter::SetBackgroundGenericFilter() void SetBackgroundGenericFilter::Update() { // Read the Dimension and PixelType - int Dimension; + int Dimension, Components; std::string PixelType; - ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType); + //ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType); + ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType, Components); - - // Call UpdateWithDim - if(Dimension==2) UpdateWithDim<2>(PixelType); - else if(Dimension==3) UpdateWithDim<3>(PixelType); - else if (Dimension==4)UpdateWithDim<4>(PixelType); - else { + if (Dimension > 4) { std::cout<<"Error, Only for 2 , 3 or 4 Dimensions!!!"< 3) { + std::cout<<"Error, Only 1, 2, or 3-component images are supported!!!"< 1 && PixelType != "float") { + std::cout<<"Error, Only float multi-component images are supported!!!"<(PixelType); + else if(Dimension==3) UpdateWithDim<3>(PixelType); + else if (Dimension==4)UpdateWithDim<4>(PixelType); + break; + case 2: + { + typedef itk::Vector TPixelType; + if(Dimension==2) UpdateWithDimAndPixelType<2, TPixelType>(); + else if(Dimension==3) UpdateWithDimAndPixelType<3, TPixelType>(); + else if (Dimension==4)UpdateWithDimAndPixelType<4, TPixelType>(); + break; + } + case 3: + { + typedef itk::Vector TPixelType; + if(Dimension==2) UpdateWithDimAndPixelType<2, TPixelType>(); + else if(Dimension==3) UpdateWithDimAndPixelType<3, TPixelType>(); + else if (Dimension==4)UpdateWithDimAndPixelType<4, TPixelType>(); + break; + } + } } diff --git a/tools/clitkSetBackgroundGenericFilter.txx b/tools/clitkSetBackgroundGenericFilter.txx index 8a825ad..63ccfb0 100644 --- a/tools/clitkSetBackgroundGenericFilter.txx +++ b/tools/clitkSetBackgroundGenericFilter.txx @@ -54,7 +54,6 @@ SetBackgroundGenericFilter::UpdateWithDim(std::string PixelType) } } - //------------------------------------------------------------------- // Update with the number of dimensions and the pixeltype //------------------------------------------------------------------- @@ -65,7 +64,7 @@ SetBackgroundGenericFilter::UpdateWithDimAndPixelType() // ImageTypes typedef itk::Image InputImageType; - typedef itk::Image MaskImageType; + typedef itk::Image MaskImageType; // Read the input typedef itk::ImageFileReader InputReaderType; diff --git a/vv/qt_ui/vvToolRigidReg.ui b/vv/qt_ui/vvToolRigidReg.ui index 10e77d6..02eb98e 100644 --- a/vv/qt_ui/vvToolRigidReg.ui +++ b/vv/qt_ui/vvToolRigidReg.ui @@ -369,13 +369,13 @@ p, li { white-space: pre-wrap; } - + - + diff --git a/vv/vvToolRigidReg.cxx b/vv/vvToolRigidReg.cxx index 57d7bcb..28591ba 100644 --- a/vv/vvToolRigidReg.cxx +++ b/vv/vvToolRigidReg.cxx @@ -108,9 +108,9 @@ void vvToolRigidReg::InputIsSelected(vvSlicerManager *input) imageorigin=mInput->GetImage()->GetOrigin(); std::vector imageSize = mInput->GetImage()->GetSize(); std::vector imageSpacing = mInput->GetImage()->GetSpacing(); - xcord=xcord.setNum(imageorigin[0]+imageSize[0]*imageSpacing[0]/2, 'g', 3); - ycord=ycord.setNum(imageorigin[1]+imageSize[1]*imageSpacing[1]/2, 'g', 3); - zcord=zcord.setNum(imageorigin[2]+imageSize[2]*imageSpacing[2]/2, 'g', 3); + xcord=xcord.setNum(imageorigin[0]+(imageSize[0]-1)*imageSpacing[0]*0.5, 'g', 3); + ycord=ycord.setNum(imageorigin[1]+(imageSize[1]-1)*imageSpacing[1]*0.5, 'g', 3); + zcord=zcord.setNum(imageorigin[2]+(imageSize[2]-1)*imageSpacing[2]*0.5, 'g', 3); Xval->setText(xcord); Yval->setText(ycord); Zval->setText(zcord); @@ -370,7 +370,6 @@ void vvToolRigidReg::SetTransform(vtkMatrix4x4 *matrix) euler->SetMatrix(rotMat); euler->SetOffset(transVec); - // Modify GUI according to the new parameters std::vector transSliders, rotSliders; std::vector transSBs, rotSBs;