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