itk/clitkInvertVFFilter.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 __clitkInvertVFFilter_txx
  19 #define __clitkInvertVFFilter_txx
  20
  21 namespace
  22 {
  23
  24 //=========================================================================================================================
  25 //helper class 1 to allow a threaded execution: add contributions of input to output and update weights
  26 //=========================================================================================================================
  27 template<class InputImageType, class OutputImageType> class ITK_EXPORT HelperClass1 : public itk::ImageToImageFilter<InputImageType, OutputImageType>
  28 {
  29
  30 public:
  31   /** Standard class typedefs. */
  32   typedef HelperClass1  Self;
  33   typedef itk::ImageToImageFilter<InputImageType,OutputImageType> Superclass;
  34   typedef itk::SmartPointer<Self>         Pointer;
  35   typedef itk::SmartPointer<const Self>   ConstPointer;
  36
  37   /** Method for creation through the object factory. */
  38   itkNewMacro(Self);
  39
  40   /** Run-time type information (and related methods) */
  41   itkTypeMacro( HelperClass1, ImageToImageFilter );
  42
  43   /** Constants for the image dimensions */
  44   itkStaticConstMacro(ImageDimension, unsigned int,InputImageType::ImageDimension);
  45
  46
  47   //Typedefs
  48   typedef typename OutputImageType::PixelType        PixelType;
  49   typedef itk::Image<double, ImageDimension > WeightsImageType;
  50   typedef itk::Image<itk::SimpleFastMutexLock, ImageDimension > MutexImageType;
  51
  52   //===================================================================================
  53   //Set methods
  54   void SetWeights(const typename WeightsImageType::Pointer input) {
  55     m_Weights = input;
  56     this->Modified();
  57   }
  58   void SetMutexImage(const typename MutexImageType::Pointer input) {
  59     m_MutexImage=input;
  60     this->Modified();
  61     m_ThreadSafe=true;
  62   }
  63
  64   //Get methods
  65   typename  WeightsImageType::Pointer GetWeights() {
  66     return m_Weights;
  67   }
  68
  69   /** Typedef to describe the output image region type. */
  70   typedef typename OutputImageType::RegionType OutputImageRegionType;
  71
  72 protected:
  73   HelperClass1();
  74   ~HelperClass1() {};
  75
  76   //the actual processing
  77   void BeforeThreadedGenerateData() ITK_OVERRIDE;
  78   void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId ) ITK_OVERRIDE;
  79
  80   //member data
  81   typename  WeightsImageType::Pointer m_Weights;
  82   typename MutexImageType::Pointer m_MutexImage;
  83   bool m_ThreadSafe;
  84
  85 };
  86
  87
  88
  89 //=========================================================================================================================
  90 //Member functions of the helper class 1
  91 //=========================================================================================================================
  92
  93
  94 //=========================================================================================================================
  95 //Empty constructor
  96 template<class InputImageType, class OutputImageType >
  97 HelperClass1<InputImageType, OutputImageType>::HelperClass1()
  98 {
  99   m_ThreadSafe=false;
 100 }
 101
 102 //=========================================================================================================================
 103 //Before threaded data
 104 template<class InputImageType, class OutputImageType >
 105 void HelperClass1<InputImageType, OutputImageType>::BeforeThreadedGenerateData()
 106 {
 107   //std::cout << "HelperClass1::BeforeThreadedGenerateData - IN" << std::endl;
 108   //Since we will add, put to zero!
 109   this->GetOutput()->FillBuffer(itk::NumericTraits<double>::Zero);
 110   this->GetWeights()->FillBuffer(itk::NumericTraits<double>::Zero);
 111 }
 112
 113 //=========================================================================================================================
 114 //update the output for the outputRegionForThread
 115 template<class InputImageType, class OutputImageType>
 116 void HelperClass1<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId )
 117 {
 118 //   std::cout << "HelperClass1::ThreadedGenerateData - IN " << threadId << std::endl;
 119   //Get pointer to the input
 120   typename InputImageType::ConstPointer inputPtr = this->GetInput();
 121
 122   //Get pointer to the output
 123   typename OutputImageType::Pointer outputPtr = this->GetOutput();
 124   //typename OutputImageType::SizeType size=outputPtr->GetLargestPossibleRegion().GetSize();
 125
 126   //Iterator over input
 127   typedef itk::ImageRegionConstIteratorWithIndex<InputImageType> InputImageIteratorType;
 128
 129   //define them over the outputRegionForThread
 130   InputImageIteratorType inputIt(inputPtr, outputRegionForThread);
 131
 132   //Initialize
 133   typename InputImageType::IndexType index;
 134   itk::ContinuousIndex<double,ImageDimension> contIndex, inContIndex;
 135   typename InputImageType::PointType ipoint;
 136   typename OutputImageType::PointType opoint;
 137   typedef typename OutputImageType::PixelType DisplacementType;
 138   DisplacementType displacement;
 139   inputIt.GoToBegin();
 140
 141   typename OutputImageType::SizeType size = outputPtr->GetLargestPossibleRegion().GetSize();
 142
 143   //define some temp variables
 144   signed long baseIndex[ImageDimension];
 145   double distance[ImageDimension];
 146   unsigned int dim, counter, upper;
 147   double totalOverlap,overlap;
 148   typename OutputImageType::IndexType neighIndex;
 149
 150   //Find the number of neighbors
 151   unsigned int neighbors =  1 << ImageDimension;
 152
 153   //==================================================================================================
 154   //Loop over the output region and add the intensities from the input to the output and the weight to the weights
 155   //==================================================================================================
 156   while( !inputIt.IsAtEnd() ) {
 157     // get the input image index
 158     index = inputIt.GetIndex();
 159     inputPtr->TransformIndexToPhysicalPoint( index,ipoint );
 160
 161     // get the required displacement
 162     displacement = inputIt.Get();
 163
 164     // compute the required output image point
 165     for(unsigned int j = 0; j < ImageDimension; j++ ) opoint[j] = ipoint[j] + (double)displacement[j];
 166
 167     // Update the output and the weights
 168     if(outputPtr->TransformPhysicalPointToContinuousIndex(opoint, contIndex ) ) {
 169       for(dim = 0; dim < ImageDimension; dim++) {
 170         // The following  block is equivalent to the following line without
 171         // having to call floor. (Only for positive inputs, we already now that is in the image)
 172         // baseIndex[dim] = (long) std::floor(contIndex[dim] );
 173
 174         baseIndex[dim] = (long) contIndex[dim];
 175         distance[dim] = contIndex[dim] - double( baseIndex[dim] );
 176       }
 177
 178       //Add contribution for each neighbor
 179       totalOverlap = itk::NumericTraits<double>::Zero;
 180       for( counter = 0; counter < neighbors ; counter++ ) {
 181         overlap = 1.0;          // fraction overlap
 182         upper = counter;  // each bit indicates upper/lower neighbour
 183
 184         // get neighbor index and overlap fraction
 185         for( dim = 0; dim < ImageDimension; dim++ ) {
 186           if ( upper & 1 ) {
 187             neighIndex[dim] = baseIndex[dim] + 1;
 188             overlap *= distance[dim];
 189           } else {
 190             neighIndex[dim] = baseIndex[dim];
 191             overlap *= 1.0 - distance[dim];
 192           }
 193           upper >>= 1;
 194
 195           if (neighIndex[dim] >= size[dim])
 196             neighIndex[dim] = size[dim] - 1;
 197         }
 198
 199
 200         //Set neighbor value only if overlap is not zero
 201         if( (overlap>0.0)) // &&
 202           //                    (static_cast<unsigned int>(neighIndex[0])<size[0]) &&
 203           //                    (static_cast<unsigned int>(neighIndex[1])<size[1]) &&
 204           //                    (static_cast<unsigned int>(neighIndex[2])<size[2]) &&
 205           //                    (neighIndex[0]>=0) &&
 206           //                    (neighIndex[1]>=0) &&
 207           //                    (neighIndex[2]>=0) )
 208         {
 209           //what to store? the original displacement vector?
 210           if (! m_ThreadSafe) {
 211             //Set the pixel and weight at neighIndex
 212             outputPtr->SetPixel(neighIndex, outputPtr->GetPixel(neighIndex) - (displacement*overlap));
 213             m_Weights->SetPixel(neighIndex, m_Weights->GetPixel(neighIndex) + overlap);
 214           }
 215
 216           else {
 217             //Entering critilal section: shared memory
 218             m_MutexImage->GetPixel(neighIndex).Lock();
 219
 220             //Set the pixel and weight at neighIndex
 221             outputPtr->SetPixel(neighIndex, outputPtr->GetPixel(neighIndex) - (displacement*overlap));
 222             m_Weights->SetPixel(neighIndex, m_Weights->GetPixel(neighIndex) + overlap);
 223
 224             //Unlock
 225             m_MutexImage->GetPixel(neighIndex).Unlock();
 226
 227           }
 228           //Add to total overlap
 229           totalOverlap += overlap;
 230         }
 231
 232         if( totalOverlap == 1.0 ) {
 233           // finished
 234           break;
 235         }
 236       }
 237     }
 238
 239     ++inputIt;
 240   }
 241
 242 //   std::cout << "HelperClass1::ThreadedGenerateData - OUT " << threadId << std::endl;
 243 }
 244
 245
 246
 247 //=========================================================================================================================
 248 //helper class 2 to allow a threaded execution of normalisation by the weights
 249 //=========================================================================================================================
 250 template<class InputImageType, class OutputImageType> class HelperClass2 : public itk::ImageToImageFilter<InputImageType, OutputImageType>
 251 {
 252
 253 public:
 254   /** Standard class typedefs. */
 255   typedef HelperClass2  Self;
 256   typedef itk::ImageToImageFilter<InputImageType,OutputImageType> Superclass;
 257   typedef itk::SmartPointer<Self>         Pointer;
 258   typedef itk::SmartPointer<const Self>   ConstPointer;
 259
 260   /** Method for creation through the object factory. */
 261   itkNewMacro(Self);
 262
 263   /** Run-time type information (and related methods) */
 264   itkTypeMacro( HelperClass2, ImageToImageFilter );
 265
 266   /** Constants for the image dimensions */
 267   itkStaticConstMacro(ImageDimension, unsigned int,InputImageType::ImageDimension);
 268
 269   //Typedefs
 270   typedef typename OutputImageType::PixelType        PixelType;
 271   typedef itk::Image<double,ImageDimension> WeightsImageType;
 272
 273   //Set methods
 274   void SetWeights(const typename WeightsImageType::Pointer input) {
 275     m_Weights = input;
 276     this->Modified();
 277   }
 278   void SetEdgePaddingValue(PixelType value) {
 279     m_EdgePaddingValue = value;
 280     this->Modified();
 281   }
 282
 283   /** Typedef to describe the output image region type. */
 284   typedef typename OutputImageType::RegionType OutputImageRegionType;
 285
 286 protected:
 287   HelperClass2();
 288   ~HelperClass2() {};
 289
 290
 291   //the actual processing
 292   void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId ) ITK_OVERRIDE;
 293
 294   //member data
 295   typename     WeightsImageType::Pointer m_Weights;
 296   PixelType m_EdgePaddingValue;
 297
 298 } ;
 299
 300
 301
 302 //=========================================================================================================================
 303 //Member functions of the helper class 2
 304 //=========================================================================================================================
 305
 306
 307 //=========================================================================================================================
 308 //Empty constructor
 309 template<class InputImageType, class OutputImageType > HelperClass2<InputImageType, OutputImageType>::HelperClass2()
 310 {
 311   PixelType zero;
 312   for(unsigned int i=0;i <PixelType::Dimension; i++) zero[i] = 0.0;
 313   m_EdgePaddingValue=zero;
 314   //m_EdgePaddingValue=itk::NumericTraits<PixelType>::Zero;
 315 }
 316
 317
 318 //=========================================================================================================================
 319 //update the output for the outputRegionForThread
 320 template<class InputImageType, class OutputImageType > void HelperClass2<InputImageType, OutputImageType>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId )
 321 {
 322 //   std::cout << "HelperClass2::ThreadedGenerateData - IN " << threadId << std::endl;
 323
 324   //Get pointer to the input
 325   typename InputImageType::ConstPointer inputPtr = this->GetInput();
 326
 327   //Get pointer to the output
 328   typename OutputImageType::Pointer outputPtr = this->GetOutput();
 329
 330   //Iterators over input, weigths  and output
 331   typedef itk::ImageRegionConstIterator<InputImageType> InputImageIteratorType;
 332   typedef itk::ImageRegionIterator<OutputImageType> OutputImageIteratorType;
 333   typedef itk::ImageRegionIterator<WeightsImageType> WeightsImageIteratorType;
 334
 335   //define them over the outputRegionForThread
 336   OutputImageIteratorType outputIt(outputPtr, outputRegionForThread);
 337   InputImageIteratorType inputIt(inputPtr, outputRegionForThread);
 338   WeightsImageIteratorType weightsIt(m_Weights, outputRegionForThread);
 339
 340
 341   //==================================================================================================
 342   //loop over the output and normalize the input, remove holes
 343   PixelType neighValue;
 344   double  zero = itk::NumericTraits<double>::Zero;
 345   while (!outputIt.IsAtEnd()) {
 346     //the weight is not zero
 347     if (weightsIt.Get() != zero) {
 348       //divide by the weight
 349       outputIt.Set(static_cast<PixelType>(inputIt.Get()/weightsIt.Get()));
 350     }
 351
 352     //copy the value of the  neighbour that was just processed
 353     else {
 354       if(!outputIt.IsAtBegin()) {
 355         //go back
 356         --outputIt;
 357
 358         //Neighbour cannot have zero weight because it should be filled already
 359         neighValue=outputIt.Get();
 360         ++outputIt;
 361         outputIt.Set(neighValue);
 362         //DD("hole filled");
 363       } else {
 364         //DD("is at begin, setting edgepadding value");
 365         outputIt.Set(m_EdgePaddingValue);
 366       }
 367     }
 368     ++weightsIt;
 369     ++outputIt;
 370     ++inputIt;
 371
 372   }//end while
 373
 374 //   std::cout << "HelperClass2::ThreadedGenerateData - OUT " << threadId << std::endl;
 375
 376 }//end member
 377
 378
 379 }//end nameless namespace
 380
 381
 382
 383 namespace clitk
 384 {
 385
 386 //=========================================================================================================================
 387 // The rest is the InvertVFFilter
 388 //=========================================================================================================================
 389
 390 //=========================================================================================================================
 391 //constructor
 392 template <class InputImageType, class OutputImageType>
 393 InvertVFFilter<InputImageType, OutputImageType>::InvertVFFilter()
 394 {
 395
 396   //m_EdgePaddingValue=itk::NumericTraits<PixelType>::Zero; //no other reasonable value?
 397   PixelType zero;
 398   for(unsigned int i=0;i <PixelType::Dimension; i++) zero[i] = 0.0;
 399   m_EdgePaddingValue=zero; //no other reasonable value?
 400
 401   m_ThreadSafe=false;
 402   m_Verbose=false;
 403 }
 404
 405 //=========================================================================================================================
 406 //Update
 407 template <class InputImageType, class OutputImageType> void InvertVFFilter<InputImageType, OutputImageType>::GenerateData()
 408 {
 409   // std::cout << "InvertVFFilter::GenerateData - IN" << std::endl;
 410
 411   //Get the properties of the input
 412   typename InputImageType::ConstPointer inputPtr=this->GetInput();
 413   typename WeightsImageType::RegionType region = inputPtr->GetLargestPossibleRegion();
 414
 415   //Allocate the weights
 416   typename WeightsImageType::Pointer weights=WeightsImageType::New();
 417   weights->SetOrigin(inputPtr->GetOrigin());
 418   weights->SetRegions(region);
 419   weights->Allocate();
 420   weights->SetSpacing(inputPtr->GetSpacing());
 421
 422   //===========================================================================
 423   //Inversion is divided in in two loops, for each we will call a threaded helper class
 424   //1. add contribution of input to output and update weights
 425   //2. normalize the output by the weight and remove holes
 426   //===========================================================================
 427
 428
 429   //===========================================================================
 430   //1. add contribution of input to output and update weights
 431
 432   //Define an internal image type
 433
 434   typedef itk::Image<itk::Vector<double,ImageDimension>, ImageDimension > InternalImageType;
 435
 436   //Call threaded helper class 1
 437   typedef HelperClass1<InputImageType, InternalImageType > HelperClass1Type;
 438   typename HelperClass1Type::Pointer helper1=HelperClass1Type::New();
 439
 440   //Set input
 441   if(m_NumberOfThreadsIsGiven) {
 442 #if ITK_VERSION_MAJOR <= 4
 443     helper1->SetNumberOfThreads(m_NumberOfThreads);
 444 #else
 445     helper1->SetNumberOfWorkUnits(m_NumberOfWorkUnits);
 446 #endif
 447   }
 448   helper1->SetInput(inputPtr);
 449   helper1->SetWeights(weights);
 450
 451   //Threadsafe?
 452   if(m_ThreadSafe) {
 453     //Allocate the mutex image
 454     typename MutexImageType::Pointer mutex=InvertVFFilter::MutexImageType::New();
 455     mutex->SetRegions(region);
 456     mutex->Allocate();
 457     mutex->SetSpacing(inputPtr->GetSpacing());
 458     helper1->SetMutexImage(mutex);
 459     if (m_Verbose) std::cout <<"Inverting using a thread-safe algorithm" <<std::endl;
 460   } else  if(m_Verbose)std::cout <<"Inverting using a thread-unsafe algorithm" <<std::endl;
 461
 462   //Execute helper class
 463   helper1->Update();
 464
 465   //Get the output
 466   typename InternalImageType::Pointer temp= helper1->GetOutput();
 467   weights=helper1->GetWeights();
 468
 469
 470   //===========================================================================
 471   //2. Normalize the output by the weights and remove holes
 472   //Call threaded helper class
 473   typedef HelperClass2<InternalImageType, OutputImageType> HelperClass2Type;
 474   typename HelperClass2Type::Pointer helper2=HelperClass2Type::New();
 475
 476   //Set temporary output as input
 477   if(m_NumberOfThreadsIsGiven) {
 478 #if ITK_VERSION_MAJOR <= 4
 479     helper2->SetNumberOfThreads(m_NumberOfThreads);
 480 #else
 481     helper2->SetNumberOfWorkUnits(m_NumberOfWorkUnits);
 482 #endif
 483   }
 484   helper2->SetInput(temp);
 485   helper2->SetWeights(weights);
 486   helper2->SetEdgePaddingValue(m_EdgePaddingValue);
 487
 488   //Execute helper class
 489   if (m_Verbose) std::cout << "Normalizing the output VF..."<<std::endl;
 490   helper2->Update();
 491
 492   //Set the output
 493   this->SetNthOutput(0, helper2->GetOutput());
 494
 495   //std::cout << "InvertVFFilter::GenerateData - OUT" << std::endl;
 496 }
 497
 498
 499
 500 }
 501
 502 #endif