template \
itk::Image<PixelType, Dim>::Pointer clitk::ImageToImageGenericFilterBase::GetInput<itk::Image<PixelType, Dim> >(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<itk::Image<itk::Vector<float, Compo>, Dim> >(itk::Image<itk::Vector<float, Compo>,Dim>::Pointer output); \
+ void clitk::ImageToImageGenericFilterBase::SetNextOutput<itk::Image<itk::Vector<PixelType, Compo>, Dim> >(itk::Image<itk::Vector<PixelType, Compo>,Dim>::Pointer output); \
template \
- itk::Image<itk::Vector<float,Compo>, Dim>::Pointer clitk::ImageToImageGenericFilterBase::GetInput<itk::Image<itk::Vector<float, Compo>, Dim> >(unsigned int n);
+ itk::Image<itk::Vector<PixelType,Compo>, Dim>::Pointer clitk::ImageToImageGenericFilterBase::GetInput<itk::Image<itk::Vector<PixelType, Compo>, Dim> >(unsigned int n);
DEF_SetNextOutput_And_GetInput(char, 2);
DEF_SetNextOutput_And_GetInput(unsigned char, 2);
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);
#include <itkImageSeriesReader.h>
#include <itkImageToVTKImageFilter.h>
#include <itkAnalyzeImageIO.h>
-#include <itkVectorCastImageFilter.h>
+#include <itkFlexibleVectorCastImageFilter.h>
#include <vtkTransform.h>
void vvImageReader::UpdateWithDim(std::string InputPixelType)
{
if (mType == VECTORFIELD || mType == VECTORFIELDWITHTIME)
- UpdateWithDimAndInputVectorPixelType<itk::Vector<float,VImageDimension>,VImageDimension>();
+ {
+ if (VImageDimension == 4)
+ UpdateWithDimAndInputVectorPixelType<itk::Vector<float,3>,VImageDimension>();
+ else
+ UpdateWithDimAndInputVectorPixelType<itk::Vector<float,VImageDimension>,VImageDimension>();
+ }
else if (InputPixelType == "short")
UpdateWithDimAndInputPixelType<short,VImageDimension>();
else if (InputPixelType == "unsigned_short")
}
analyzeImageIO = dynamic_cast<itk::AnalyzeImageIO*>( reader->GetImageIO() );
}
-
+
typedef itk::Image< itk::Vector<float , 3>, VImageDimension > VectorImageType;
- typedef itk::VectorCastImageFilter<InputImageType, VectorImageType> CasterType;
+ typedef itk::FlexibleVectorCastImageFilter<InputImageType, VectorImageType> CasterType;
typename VectorImageType::Pointer casted_input;
typename CasterType::Pointer caster = CasterType::New();
caster->SetInput(input);
+ caster->Update();
casted_input = caster->GetOutput();
-
+
mImage = vvImageFromITK<VImageDimension, itk::Vector<float, 3> >(casted_input, mType == IMAGEWITHTIME || mType == VECTORFIELDWITHTIME);
// For unknown analyze orientations, we set identity
===========================================================================**/
#ifndef __clitkSetBackgroundImageFilter_h
#define __clitkSetBackgroundImageFilter_h
-#include "itkBinaryFunctorImageFilter.h"
+#include "itkFlexibleBinaryFunctorImageFilter.h"
#include "itkNumericTraits.h"
template <class TInputImage, class TMaskImage, class TOutputImage=TInputImage>
class ITK_EXPORT SetBackgroundImageFilter :
public
- itk::BinaryFunctorImageFilter<TInputImage,TMaskImage,TOutputImage,
+ itk::FlexibleBinaryFunctorImageFilter<TInputImage,TMaskImage,TOutputImage,
Functor::SetBackground<
typename TInputImage::PixelType,
typename TMaskImage::PixelType,
public:
/** Standard class typedefs. */
typedef SetBackgroundImageFilter Self;
- typedef itk::BinaryFunctorImageFilter<TInputImage,TMaskImage,TOutputImage,
+ typedef itk::FlexibleBinaryFunctorImageFilter<TInputImage,TMaskImage,TOutputImage,
Functor::SetBackground<
typename TInputImage::PixelType,
typename TMaskImage::PixelType,
/** Runtime information support. */
itkTypeMacro(SetBackgroundImageFilter,
- BinaryFunctorImageFilter);
+ FlexibleBinaryFunctorImageFilter);
/** Method to explicitly set the outside value of the mask. Defaults to 0 */
void SetOutsideValue( const typename TOutputImage::PixelType & outsideValue )
--- /dev/null
+/*=========================================================================
+
+ Program: Insight Segmentation & Registration Toolkit
+ Module: $RCSfile: itkFlexibleBinaryFunctorImageFilter.h,v $
+ Language: C++
+ Date: $Date: 2012-02-17 11:01:02 $
+ Version: $Revision: 1.37 $
+
+ 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_h
+#define __itkFlexibleBinaryFunctorImageFilter_h
+
+#include "itkInPlaceImageFilter.h"
+#include "itkImageRegionIteratorWithIndex.h"
+
+namespace itk
+{
+
+/** \class FlexibleBinaryFunctorImageFilter
+ * \brief Implements pixel-wise generic operation of two images.
+ *
+ * This class is parameterized over the types of the two input images
+ * and the type of the output image. It is also parameterized by the
+ * operation to be applied. A Functor style is used. The main difference
+ * wrt the standard BinaryFunctorImageFilter is that the second input's
+ * dimensions/spacing/origin can differ from those of the first. The
+ * functor is then only called if the pixel from input1 is inside
+ * input2's largest region.
+ *
+ * \sa UnaryFunctorImageFilter BinaryFunctorImageFilter TernaryFunctorImageFilter
+ *
+ * \ingroup IntensityImageFilters Multithreaded
+ */
+template <class TInputImage1, class TInputImage2,
+ class TOutputImage, class TFunction >
+class ITK_EXPORT FlexibleBinaryFunctorImageFilter :
+ public InPlaceImageFilter<TInputImage1,TOutputImage>
+{
+public:
+ /** Standard class typedefs. */
+ typedef FlexibleBinaryFunctorImageFilter Self;
+ typedef InPlaceImageFilter<TInputImage1,TOutputImage> Superclass;
+ typedef SmartPointer<Self> Pointer;
+ typedef SmartPointer<const Self> 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<itkGetStaticConstMacro(InputImage1Dimension),
+ itkGetStaticConstMacro(InputImage2Dimension)>));
+ itkConceptMacro(SameDimensionCheck2,
+ (Concept::SameDimension<itkGetStaticConstMacro(InputImage1Dimension),
+ itkGetStaticConstMacro(OutputImageDimension)>));
+ /** 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
--- /dev/null
+/*=========================================================================
+
+ 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 <class TInputImage1, class TInputImage2,
+ class TOutputImage, class TFunction >
+FlexibleBinaryFunctorImageFilter<TInputImage1,TInputImage2,TOutputImage,TFunction>
+::FlexibleBinaryFunctorImageFilter()
+{
+ this->SetNumberOfRequiredInputs( 1 );
+ this->SetInPlace(false);
+}
+
+
+/**
+ * Connect one of the operands for pixel-wise addition
+ */
+template <class TInputImage1, class TInputImage2,
+ class TOutputImage, class TFunction >
+void
+FlexibleBinaryFunctorImageFilter<TInputImage1,TInputImage2,TOutputImage,TFunction>
+::SetInput1( const TInputImage1 * image1 )
+{
+ // Process object is not const-correct so the const casting is required.
+ this->SetNthInput(0, const_cast<TInputImage1 *>( image1 ));
+}
+
+
+/**
+ * Connect one of the operands for pixel-wise addition
+ */
+template <class TInputImage1, class TInputImage2,
+ class TOutputImage, class TFunction >
+void
+FlexibleBinaryFunctorImageFilter<TInputImage1,TInputImage2,TOutputImage,TFunction>
+::SetInput2( const TInputImage2 * image2 )
+{
+ // Process object is not const-correct so the const casting is required.
+ //this->SetNthInput(1, const_cast<TInputImage2 *>( image2 ));
+
+ m_Input2 = image2;
+}
+
+template <class TInputImage1, class TInputImage2,
+ class TOutputImage, class TFunction >
+void
+FlexibleBinaryFunctorImageFilter<TInputImage1,TInputImage2,TOutputImage,TFunction>
+::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<TInputImage2 *>( m_Input2.GetPointer() );
+ image2->Update();
+}
+
+template <class TInputImage1, class TInputImage2,
+ class TOutputImage, class TFunction >
+void
+FlexibleBinaryFunctorImageFilter<TInputImage1,TInputImage2,TOutputImage,TFunction>
+::GenerateOutputInformation()
+{
+ Superclass::GenerateOutputInformation() ;
+}
+
+/**
+ * ThreadedGenerateData Performs the pixel-wise addition
+ */
+template <class TInputImage1, class TInputImage2, class TOutputImage, class TFunction >
+void
+FlexibleBinaryFunctorImageFilter<TInputImage1, TInputImage2, TOutputImage, TFunction>
+::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<const TInputImage1*>(ProcessObject::GetInput(0));
+ Input2ImagePointer inputPtr2 = m_Input2;
+/* = dynamic_cast<const TInputImage2*>(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<TInputImage1> inputIt1(inputPtr1, outputRegionForThread);
+ ImageRegionConstIterator<TInputImage2> inputIt2(inputPtr2, outputRegionForThread);
+
+ ImageRegionIterator<TOutputImage> 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
--- /dev/null
+/*=========================================================================
+
+ 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<OutputValueType>( A[k] );
+ }
+
+ for (; k < TOutput::Dimension; k++)
+ value[k] = 0;
+
+ return value;
+ }
+};
+}
+
+template <class TInputImage, class TOutputImage = itk::Vector<float, 3> >
+class ITK_EXPORT FlexibleVectorCastImageFilter :
+ public
+UnaryFunctorImageFilter<TInputImage,TOutputImage,
+ Functor::FlexibleVectorCast< typename TInputImage::PixelType,
+ typename TOutputImage::PixelType> >
+{
+public:
+ /** Standard class typedefs. */
+ typedef FlexibleVectorCastImageFilter Self;
+ typedef UnaryFunctorImageFilter<
+ TInputImage,TOutputImage,
+ Functor::FlexibleVectorCast< typename TInputImage::PixelType,
+ typename TOutputImage::PixelType> > Superclass;
+ typedef SmartPointer<Self> Pointer;
+ typedef SmartPointer<const Self> 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
# 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.*
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
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?)
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..."
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
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
# 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
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..."
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.*
+}
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
#define CLITKIMAGEARITHMGENERICFILTER_CXX
#include "clitkImageArithmGenericFilter.h"
-#include "clitkImageArithm_ggo.h"
namespace clitk {
// Specialisation
+// template<>
+// class ImageArithmGenericFilter<args_info_clitkImageArithm>;
+
+////////////////////////////////////////////////////////////////////
+ // specializations for itk::Vector<float, 3u>, 3u
+ template<>
+ template<>
+ void ImageArithmGenericFilter<args_info_clitkImageArithm>::UpdateWithInputImageType< itk::Image< itk::Vector<float, 3u>, 3u > >()
+ {
+ typedef itk::Image< itk::Vector<float, 3u>, 3u > ImageType;
+
+ // Read input1
+ ImageType::Pointer input1 = this->GetInput<ImageType>(0);
+
+ // Set input image iterator
+ typedef itk::ImageRegionIterator<ImageType> IteratorType;
+ IteratorType it(input1, input1->GetLargestPossibleRegion());
+
+ // typedef input2
+ ImageType::Pointer input2 = NULL;
+ IteratorType it2;
+
+ /*
+ // Special case for normalisation
+ if (mTypeOfOperation == 12) {
+ typedef itk::MinimumMaximumImageCalculator<ImageType> 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<ImageType>(1);
+ // Set input image iterator
+ it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
+ // Check dimension
+ if (!clitk::HaveSameSize<ImageType, ImageType>(input1, input2)) {
+ itkExceptionMacro(<< "The images (input and input2) must have the same size");
+ }
+ if(!clitk::HaveSameSpacing<ImageType, ImageType>(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<ImageType>(input1);
+ }
+ // ---------------- Create new output Image ---------------------
+ else {
+ /*if (mOutputIsFloat) {
+ // Create output image
+ typedef itk::Image<float,ImageType::ImageDimension> 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<OutputImageType> IteratorOutputType;
+ IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
+ if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
+ else ComputeImage(it, ito);
+ this->template SetNextOutput<OutputImageType>(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<OutputImageType> IteratorOutputType;
+ IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
+ if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
+ else ComputeImage(it, ito);
+ this->SetNextOutput<OutputImageType>(output);
+ }
+ }
+ }
+
+ template<>
+ template<>
+ void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
+ itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > >
+ (itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > ito)
+ {
+ typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 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<double, PixelType>(mScalar));
+ else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
+ ++it;
+ ++ito;
+ }
+ break;
+ case 4: // Min
+ while (!it.IsAtEnd()) {
+ if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
+ else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
+ ++it;
+ ++ito;
+ }
+ break;
+ case 5: // Absolute value
+ while (!it.IsAtEnd()) {
+ if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
+ // <= zero to avoid warning for unsigned types
+ else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
+ ++it;
+ ++ito;
+ }
+ break;
+ case 6: // Squared value
+ while (!it.IsAtEnd()) {
+ ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
+ ++it;
+ ++ito;
+ }
+ break;
+ case 7: // Log
+ while (!it.IsAtEnd()) {
+ if (it.Get() > 0)
+ ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
+ else ito.Set(mDefaultPixelValue);
+ ++it;
+ ++ito;
+ }
+ break;
+ case 8: // exp
+ while (!it.IsAtEnd()) {
+ ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
+ ++it;
+ ++ito;
+ }
+ break;
+ case 9: // sqrt
+ while (!it.IsAtEnd()) {
+ if (it.Get() > 0)
+ ito.Set(PixelTypeDownCast<double, PixelType>(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<double, PixelType>((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<args_info_clitkImageArithm>::ComputeImage<
+ itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >
+ >
+ (itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it1,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it2,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > ito)
+ {
+ typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter1;
+ typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter2;
+ typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 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, PixelType>(((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<double, 3u>, 3u
+ template<>
+ template<>
+ void ImageArithmGenericFilter<args_info_clitkImageArithm>::UpdateWithInputImageType< itk::Image< itk::Vector<double, 3u>, 3u > >()
+ {
+ typedef itk::Image< itk::Vector<double, 3u>, 3u > ImageType;
+
+ // Read input1
+ ImageType::Pointer input1 = this->GetInput<ImageType>(0);
+
+ // Set input image iterator
+ typedef itk::ImageRegionIterator<ImageType> IteratorType;
+ IteratorType it(input1, input1->GetLargestPossibleRegion());
+
+ // typedef input2
+ ImageType::Pointer input2 = NULL;
+ IteratorType it2;
+
+ /*
+ // Special case for normalisation
+ if (mTypeOfOperation == 12) {
+ typedef itk::MinimumMaximumImageCalculator<ImageType> 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<ImageType>(1);
+ // Set input image iterator
+ it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
+ // Check dimension
+ if (!clitk::HaveSameSize<ImageType, ImageType>(input1, input2)) {
+ itkExceptionMacro(<< "The images (input and input2) must have the same size");
+ }
+ if(!clitk::HaveSameSpacing<ImageType, ImageType>(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<ImageType>(input1);
+ }
+ // ---------------- Create new output Image ---------------------
+ else {
+ /*if (mOutputIsFloat) {
+ // Create output image
+ typedef itk::Image<float,ImageType::ImageDimension> 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<OutputImageType> IteratorOutputType;
+ IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
+ if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
+ else ComputeImage(it, ito);
+ this->template SetNextOutput<OutputImageType>(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<OutputImageType> IteratorOutputType;
+ IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
+ if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
+ else ComputeImage(it, ito);
+ this->SetNextOutput<OutputImageType>(output);
+ }
+ }
+ }
+
template<>
- class ImageArithmGenericFilter<args_info_clitkImageArithm>;
+ template<>
+ void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
+ itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > >
+ (itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > ito)
+ {
+ typedef itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 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<double, PixelType>(mScalar));
+ else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
+ ++it;
+ ++ito;
+ }
+ break;
+ case 4: // Min
+ while (!it.IsAtEnd()) {
+ if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
+ else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
+ ++it;
+ ++ito;
+ }
+ break;
+ case 5: // Absolute value
+ while (!it.IsAtEnd()) {
+ if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
+ // <= zero to avoid warning for unsigned types
+ else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
+ ++it;
+ ++ito;
+ }
+ break;
+ case 6: // Squared value
+ while (!it.IsAtEnd()) {
+ ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
+ ++it;
+ ++ito;
+ }
+ break;
+ case 7: // Log
+ while (!it.IsAtEnd()) {
+ if (it.Get() > 0)
+ ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
+ else ito.Set(mDefaultPixelValue);
+ ++it;
+ ++ito;
+ }
+ break;
+ case 8: // exp
+ while (!it.IsAtEnd()) {
+ ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
+ ++it;
+ ++ito;
+ }
+ break;
+ case 9: // sqrt
+ while (!it.IsAtEnd()) {
+ if (it.Get() > 0)
+ ito.Set(PixelTypeDownCast<double, PixelType>(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<double, PixelType>((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<args_info_clitkImageArithm>::ComputeImage<
+ itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >
+ >
+ (itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it1,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it2,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > ito)
+ {
+ typedef itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > Iter1;
+ typedef itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > Iter2;
+ typedef itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 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, PixelType>(((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
// clitk include
#include "clitkCommon.h"
#include "clitkImageToImageGenericFilter.h"
+#include "clitkImageArithm_ggo.h"
// itk include
#include "itkImage.h"
//--------------------------------------------------------------------
}; // end class ImageArithmGenericFilter
+
+ // specializations for itk::Vector<float, 3u>, 3u
+ template<> template<>
+ void ImageArithmGenericFilter<args_info_clitkImageArithm>::UpdateWithInputImageType< itk::Image< itk::Vector<float, 3u>, 3u > >();
+
+ template<> template<>
+ void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
+ itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >
+ >
+ (itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > ito);
+
+ template<> template<>
+ void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
+ itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >
+ >
+ (itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it1,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it2,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > ito);
+
+ // specializations for itk::Vector<double, 3u>, 3u
+ template<> template<>
+ void ImageArithmGenericFilter<args_info_clitkImageArithm>::UpdateWithInputImageType< itk::Image< itk::Vector<double, 3u>, 3u > >();
+
+ template<> template<>
+ void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
+ itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >
+ >
+ (itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > ito);
+
+ template<> template<>
+ void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
+ itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >
+ >
+ (itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it1,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it2,
+ itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > ito);
} // end namespace
//--------------------------------------------------------------------
+
#ifndef ITK_MANUAL_INSTANTIATION
#include "clitkImageArithmGenericFilter.txx"
#endif
void ImageArithmGenericFilter<args_info_type>::InitializeImageType()
{
ADD_DEFAULT_IMAGE_TYPES(Dim);
+ ADD_VEC_IMAGE_TYPE(3u,3u,float);
+ ADD_VEC_IMAGE_TYPE(3u,3u,double);
}
//--------------------------------------------------------------------
}
//--------------------------------------------------------------------
+
+
} // end namespace
#endif //#define CLITKIMAGEARITHMGENERICFILTER_TXX
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!!!"<<std::endl ;
return;
}
+
+ if (Components > 3) {
+ std::cout<<"Error, Only 1, 2, or 3-component images are supported!!!"<<std::endl ;
+ return;
+ }
+
+ if (Components > 1 && PixelType != "float") {
+ std::cout<<"Error, Only float multi-component images are supported!!!"<<std::endl ;
+ return;
+ }
+
+ switch (Components) {
+ case 1:
+ // Call UpdateWithDim
+ if(Dimension==2) UpdateWithDim<2>(PixelType);
+ else if(Dimension==3) UpdateWithDim<3>(PixelType);
+ else if (Dimension==4)UpdateWithDim<4>(PixelType);
+ break;
+ case 2:
+ {
+ typedef itk::Vector<float, 2> 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<float, 3> TPixelType;
+ if(Dimension==2) UpdateWithDimAndPixelType<2, TPixelType>();
+ else if(Dimension==3) UpdateWithDimAndPixelType<3, TPixelType>();
+ else if (Dimension==4)UpdateWithDimAndPixelType<4, TPixelType>();
+ break;
+ }
+ }
}
}
}
-
//-------------------------------------------------------------------
// Update with the number of dimensions and the pixeltype
//-------------------------------------------------------------------
// ImageTypes
typedef itk::Image<PixelType, Dimension> InputImageType;
- typedef itk::Image<float, Dimension> MaskImageType;
+ typedef itk::Image<unsigned char, Dimension> MaskImageType;
// Read the input
typedef itk::ImageFileReader<InputImageType> InputReaderType;
<item>
<layout class="QHBoxLayout" name="horizontalLayout_4">
<item>
- <widget class="QLineEdit" name="Zval"/>
+ <widget class="QLineEdit" name="Xval"/>
</item>
<item>
<widget class="QLineEdit" name="Yval"/>
</item>
<item>
- <widget class="QLineEdit" name="Xval"/>
+ <widget class="QLineEdit" name="Zval"/>
</item>
</layout>
</item>
imageorigin=mInput->GetImage()->GetOrigin();
std::vector<int> imageSize = mInput->GetImage()->GetSize();
std::vector<double> 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);
euler->SetMatrix(rotMat);
euler->SetOffset(transVec);
-
// Modify GUI according to the new parameters
std::vector<QSlider *> transSliders, rotSliders;
std::vector<QDoubleSpinBox *> transSBs, rotSBs;