itk/clitkComposeVFFilter.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 __clitkComposeVFFilter_txx
  19 #define __clitkComposeVFFilter_txx
  20 #include "clitkComposeVFFilter.h"
  21
  22
  23 namespace clitk
  24 {
  25
  26   //=========================================================================================================================
  27   //constructor
  28   template <class InputImageType, class OutputImageType>
  29   ComposeVFFilter<InputImageType, OutputImageType>::ComposeVFFilter()
  30   {
  31     for (unsigned int i=0; i<ImageDimension; i++) m_EdgePaddingValue[i]=0.0; //no other reasonable value?
  32     m_Verbose=false;
  33   }
  34
  35
  36   //=========================================================================================================================
  37   //update the output for the outputRegionForThread
  38 #if ITK_VERSION_MAJOR >= 4
  39   template<class InputImageType, class OutputImageType>
  40   void ComposeVFFilter<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId )
  41 #else
  42   template<class InputImageType, class OutputImageType>
  43   void ComposeVFFilter<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId )
  44 #endif
  45   {
  46
  47     //Get pointer to the output
  48     typename OutputImageType::Pointer outputPtr = this->GetOutput();
  49
  50     //Iterator over output
  51     typedef itk::ImageRegionConstIteratorWithIndex<InputImageType> Input1ImageIteratorType;
  52     typedef itk::ImageRegionIterator<OutputImageType> OutputImageIteratorType;
  53
  54     //define the output and input1 iterator over the outputRegionForThread
  55     OutputImageIteratorType outputIt(outputPtr, outputRegionForThread);
  56     Input1ImageIteratorType input1It(m_Input1, outputRegionForThread);
  57
  58     //Initialize
  59     typename InputImageType::IndexType index;
  60     itk::ContinuousIndex<double,ImageDimension> contIndex;
  61     typename InputImageType::PointType point1;
  62     typename InputImageType::PointType point2;
  63     typedef typename OutputImageType::PixelType DisplacementType;
  64     DisplacementType displacement;
  65     DisplacementType totalDisplacement;
  66     input1It.GoToBegin();
  67
  68     //define some temp variables outside the loop
  69     signed long baseIndex[ImageDimension];
  70     //    long tIndex;
  71     double distance[ImageDimension];
  72     unsigned int dim;
  73     double totalOverlap;
  74     double overlap;
  75     unsigned int upper;
  76     unsigned int counter;
  77     typename OutputImageType::IndexType neighIndex;
  78
  79
  80     //Find the number of neighbors
  81     unsigned int neighbors =  1 << ImageDimension;
  82
  83     //==================================================================================================
  84     //1. loop over the region of the deformationfield and compose the displacements
  85     while( !input1It.IsAtEnd() )
  86       {
  87         // get the input image index
  88         index = input1It.GetIndex();
  89         m_Input1->TransformIndexToPhysicalPoint( index, point1 );
  90
  91         // get the displacement of input 1
  92         displacement = input1It.Get();
  93
  94         // compute the required image point in input2
  95         for(unsigned int j = 0; j < ImageDimension; j++ ) point2[j] = point1[j] + (double)displacement[j];
  96
  97         //JV TODO add something for the borders
  98         //true if inside image
  99         if(m_Input2->TransformPhysicalPointToContinuousIndex(point2, contIndex ) )
 100           {
 101
 102             for(dim = 0; dim < ImageDimension; dim++)
 103               {
 104                 // The following  block is equivalent to the following line without
 105                 // having to call floor. (Only for positive inputs, we already now that is in the image)
 106                 // baseIndex[dim] = (long) vcl_floor(contIndex[dim] );
 107
 108                 baseIndex[dim] = (long) contIndex[dim];
 109                 distance[dim] = contIndex[dim] - double( baseIndex[dim] );
 110
 111                 //Reset the accumulator
 112                 totalDisplacement[dim]=0.0;
 113               }
 114
 115             //Add contribution for each neighbor
 116             totalOverlap = itk::NumericTraits<double>::Zero;
 117             for( counter = 0; counter < neighbors ; counter++ )
 118               {
 119                 overlap = 1.0;     // fraction overlap
 120                 upper = counter;  // each bit indicates upper/lower neighbour
 121
 122
 123                 // get neighbor index and overlap fraction
 124                 bool neighbIndexSupZero = 1;
 125                 for( dim = 0; dim < ImageDimension; dim++ )
 126                   {
 127                     if ( upper & 1 )
 128                       {
 129                         neighIndex[dim] = baseIndex[dim] + 1;
 130                         overlap *= distance[dim];
 131                       }
 132                     else
 133                       {
 134                         neighIndex[dim] = baseIndex[dim];
 135                         overlap *= 1.0 - distance[dim];
 136                       }
 137                     if (neighIndex[dim] < 0)
 138                       neighbIndexSupZero = 0;
 139                     upper >>= 1;
 140                   }
 141
 142                 //JV shouldn't we verify that the index is not over the upper border instead of zero?
 143                 // Set neighbor value only if overlap is not zero and index is still in image
 144                 if ( overlap>0.0 && neighbIndexSupZero )
 145                   {
 146                     //what to store? the weighted displacement vector of Input2?
 147                     totalDisplacement+=m_Input2->GetPixel(neighIndex)*overlap;
 148                     //JV don't we loose more by adding and verifyning each time, it will never be 1.0 before the end no?
 149                     totalOverlap += overlap;
 150                   }
 151
 152                 if( totalOverlap == 1.0 )
 153                   {
 154                     // finished
 155                     break;
 156                   }
 157               }
 158             //add the displacement of input1 and the interpolated displacement of input2 to the output
 159             outputIt.Set(static_cast<PixelType>(input1It.Get())+totalDisplacement);
 160           }
 161         //displacement from input1 goes outside input2
 162         else outputIt.Set(m_EdgePaddingValue);
 163
 164         ++input1It;
 165         ++outputIt;
 166       }
 167   }
 168
 169 }
 170
 171 #endif