registration/clitkExtractImageFilter.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://www.centreleonberard.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ===========================================================================**/
  18 #ifndef _clitkExtractImageFilter_txx
  19 #define _clitkExtractImageFilter_txx
  20 #include "clitkExtractImageFilter.h"
  21 #include "itkImageRegionIterator.h"
  22 #include "itkImageRegionConstIterator.h"
  23 #include "itkObjectFactory.h"
  24 #include "itkExtractImageFilterRegionCopier.h"
  25 #include "itkProgressReporter.h"
  26
  27
  28 namespace clitk
  29 {
  30
  31 /**
  32  *
  33  */
  34 template <class TInputImage, class TOutputImage>
  35 ExtractImageFilter<TInputImage,TOutputImage>
  36 ::ExtractImageFilter()
  37 {
  38 }
  39
  40
  41 /**
  42  *
  43  */
  44 template <class TInputImage, class TOutputImage>
  45 void
  46 ExtractImageFilter<TInputImage,TOutputImage>
  47 ::PrintSelf(std::ostream& os, itk::Indent indent) const
  48 {
  49   Superclass::PrintSelf(os,indent);
  50
  51   os << indent << "ExtractionRegion: " << m_ExtractionRegion << std::endl;
  52   os << indent << "OutputImageRegion: " << m_OutputImageRegion << std::endl;
  53 }
  54
  55
  56 template<class TInputImage, class TOutputImage>
  57 void
  58 ExtractImageFilter<TInputImage,TOutputImage>
  59 ::CallCopyOutputRegionToInputRegion(InputImageRegionType &destRegion,
  60                                     const OutputImageRegionType &srcRegion)
  61 {
  62   ExtractImageFilterRegionCopierType extractImageRegionCopier;
  63   extractImageRegionCopier(destRegion, srcRegion, m_ExtractionRegion);
  64 }
  65
  66
  67 template <class TInputImage, class TOutputImage>
  68 void
  69 ExtractImageFilter<TInputImage,TOutputImage>
  70 ::SetExtractionRegion(InputImageRegionType extractRegion)
  71 {
  72   m_ExtractionRegion = extractRegion;
  73
  74   unsigned int nonzeroSizeCount = 0;
  75   InputImageSizeType inputSize = extractRegion.GetSize();
  76   OutputImageSizeType outputSize;
  77   OutputImageIndexType outputIndex;
  78
  79   /**
  80    * check to see if the number of non-zero entries in the extraction region
  81    * matches the number of dimensions in the output image.
  82    **/
  83   for (unsigned int i = 0; i < InputImageDimension; ++i)
  84     {
  85     if (inputSize[i])
  86       {
  87       outputSize[nonzeroSizeCount] = inputSize[i];
  88       outputIndex[nonzeroSizeCount] = extractRegion.GetIndex()[i];
  89       nonzeroSizeCount++;
  90       }
  91     }
  92
  93   if (nonzeroSizeCount != OutputImageDimension)
  94     {
  95     itkExceptionMacro("Extraction Region not consistent with output image");
  96     }
  97
  98   m_OutputImageRegion.SetSize(outputSize);
  99   m_OutputImageRegion.SetIndex(outputIndex);
 100   this->Modified();
 101 }
 102
 103
 104
 105 /**
 106  * ExtractImageFilter can produce an image which is a different resolution
 107  * than its input image.  As such, ExtractImageFilter needs to provide an
 108  * implementation for GenerateOutputInformation() in order to inform
 109  * the pipeline execution model.  The original documentation of this
 110  * method is below.
 111  *
 112  * \sa ProcessObject::GenerateOutputInformaton()
 113  */
 114 template <class TInputImage, class TOutputImage>
 115 void
 116 ExtractImageFilter<TInputImage,TOutputImage>
 117 ::GenerateOutputInformation()
 118 {
 119   // do not call the superclass' implementation of this method since
 120   // this filter allows the input and the output to be of different dimensions
 121
 122   // get pointers to the input and output
 123   typename Superclass::OutputImagePointer      outputPtr = this->GetOutput();
 124   typename Superclass::InputImageConstPointer  inputPtr  = this->GetInput();
 125
 126   if ( !outputPtr || !inputPtr)
 127     {
 128     return;
 129     }
 130
 131   // Set the output image size to the same value as the extraction region.
 132   outputPtr->SetLargestPossibleRegion( m_OutputImageRegion );
 133
 134   // Set the output spacing and origin
 135   const itk::ImageBase<InputImageDimension> *phyData;
 136
 137   phyData
 138     = dynamic_cast<const itk::ImageBase<InputImageDimension>*>(this->GetInput());
 139
 140   if (phyData)
 141     {
 142     // Copy what we can from the image from spacing and origin of the input
 143     // This logic needs to be augmented with logic that select which
 144     // dimensions to copy
 145
 146     unsigned int i;
 147     const typename InputImageType::SpacingType&
 148       inputSpacing = inputPtr->GetSpacing();
 149     const typename InputImageType::DirectionType&
 150       inputDirection = inputPtr->GetDirection();
 151     const typename InputImageType::PointType&
 152       inputOrigin = inputPtr->GetOrigin();
 153
 154     typename OutputImageType::SpacingType outputSpacing;
 155     typename OutputImageType::DirectionType outputDirection;
 156     typename OutputImageType::PointType outputOrigin;
 157
 158     if ( static_cast<unsigned int>(OutputImageDimension) >
 159          static_cast<unsigned int>(InputImageDimension )    )
 160       {
 161       // copy the input to the output and fill the rest of the
 162       // output with zeros.
 163       for (i=0; i < InputImageDimension; ++i)
 164         {
 165         outputSpacing[i] = inputSpacing[i];
 166         outputOrigin[i] = inputOrigin[i];
 167         for (unsigned int dim = 0; dim < InputImageDimension; ++dim)
 168           {
 169           outputDirection[i][dim] = inputDirection[i][dim];
 170           }
 171         }
 172       for (; i < OutputImageDimension; ++i)
 173         {
 174         outputSpacing[i] = 1.0;
 175         outputOrigin[i] = 0.0;
 176         for (unsigned int dim = 0; dim < InputImageDimension; ++dim)
 177           {
 178           outputDirection[i][dim] = 0.0;
 179           }
 180         outputDirection[i][i] = 1.0;
 181         }
 182       }
 183     else
 184       {
 185       // copy the non-collapsed part of the input spacing and origing to the output
 186       int nonZeroCount = 0;
 187       for (i=0; i < InputImageDimension; ++i)
 188         {
 189         if (m_ExtractionRegion.GetSize()[i])
 190           {
 191           outputSpacing[nonZeroCount] = inputSpacing[i];
 192           outputOrigin[nonZeroCount] = inputOrigin[i];
 193           int nonZeroCount2 = 0;
 194           for (unsigned int dim = 0; dim < InputImageDimension; ++dim)
 195             {
 196             if (m_ExtractionRegion.GetSize()[dim])
 197               {
 198               outputDirection[nonZeroCount][nonZeroCount2] =
 199                 inputDirection[i][dim];
 200               ++nonZeroCount2;
 201               }
 202             }
 203           nonZeroCount++;
 204           }
 205         }
 206       }
 207
 208     // set the spacing and origin
 209     outputPtr->SetSpacing( outputSpacing );
 210     outputPtr->SetDirection( outputDirection );
 211     outputPtr->SetOrigin( outputOrigin );
 212     outputPtr->SetNumberOfComponentsPerPixel(inputPtr->GetNumberOfComponentsPerPixel() );
 213     }
 214   else
 215     {
 216     // pointer could not be cast back down
 217     itkExceptionMacro(<< "itk::ExtractImageFilter::GenerateOutputInformation "
 218                       << "cannot cast input to "
 219                       << typeid(itk::ImageBase<InputImageDimension>*).name() );
 220     }
 221 }
 222
 223 /**
 224    * ExtractImageFilter can be implemented as a multithreaded filter.
 225    * Therefore, this implementation provides a ThreadedGenerateData()
 226    * routine which is called for each processing thread. The output
 227    * image data is allocated automatically by the superclass prior to
 228    * calling ThreadedGenerateData().  ThreadedGenerateData can only
 229    * write to the portion of the output image specified by the
 230    * parameter "outputRegionForThread"
 231    *
 232    * \sa ImageToImageFilter::ThreadedGenerateData(),
 233    *     ImageToImageFilter::GenerateData()
 234    */
 235 template <class TInputImage, class TOutputImage>
 236 void
 237 ExtractImageFilter<TInputImage,TOutputImage>
 238 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
 239                        int threadId)
 240 {
 241   itkDebugMacro(<<"Actually executing");
 242
 243   // Get the input and output pointers
 244   typename Superclass::InputImageConstPointer  inputPtr = this->GetInput();
 245   typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
 246
 247   // support progress methods/callbacks
 248  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
 249
 250   // Define the portion of the input to walk for this thread
 251   InputImageRegionType inputRegionForThread;
 252   this->CallCopyOutputRegionToInputRegion(inputRegionForThread, outputRegionForThread);
 253
 254   // Define the iterators.
 255   typedef itk::ImageRegionIterator<TOutputImage> OutputIterator;
 256   typedef itk::ImageRegionConstIterator<TInputImage> InputIterator;
 257
 258   OutputIterator outIt(outputPtr, outputRegionForThread);
 259   InputIterator inIt(inputPtr, inputRegionForThread);
 260
 261   // walk the output region, and sample the input image
 262   while( !outIt.IsAtEnd() )
 263     {
 264     // copy the input pixel to the output
 265     outIt.Set( static_cast<OutputImagePixelType>(inIt.Get()));
 266     ++outIt;
 267     ++inIt;
 268     progress.CompletedPixel();
 269     }
 270 }
 271
 272 } // end namespace clitk
 273
 274 #endif