itk/clitkVectorBSplineDecompositionImageFilterWithOBD.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 _clitkVectorBSplineDecompositionImageFilterWithOBD_txx
  19 #define _clitkVectorBSplineDecompositionImageFilterWithOBD_txx
  20 #include "clitkVectorBSplineDecompositionImageFilterWithOBD.h"
  21 #include "itkImageRegionConstIteratorWithIndex.h"
  22 #include "itkImageRegionIterator.h"
  23 #include "itkProgressReporter.h"
  24 #include "itkVector.h"
  25
  26 namespace clitk
  27 {
  28
  29 /**
  30  * Constructor
  31  */
  32 template <class TInputImage, class TOutputImage>
  33 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
  34 ::VectorBSplineDecompositionImageFilterWithOBD()
  35 {
  36   m_SplineOrder = 0;
  37   int SplineOrder = 3;
  38   m_Tolerance = 1e-10;   // Need some guidance on this one...what is reasonable?
  39   m_IteratorDirection = 0;
  40   this->SetSplineOrder(SplineOrder);
  41 }
  42
  43
  44 /**
  45  * Standard "PrintSelf" method
  46  */
  47 template <class TInputImage, class TOutputImage>
  48 void
  49 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
  50 ::PrintSelf(
  51   std::ostream& os,
  52   itk::Indent indent) const
  53 {
  54   Superclass::PrintSelf( os, indent );
  55   os << indent << "Spline Order: " << m_SplineOrder << std::endl;
  56
  57 }
  58
  59
  60 template <class TInputImage, class TOutputImage>
  61 bool
  62 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
  63 ::DataToCoefficients1D()
  64 {
  65
  66   // See Unser, 1993, Part II, Equation 2.5,
  67   //   or Unser, 1999, Box 2. for an explaination.
  68
  69   double c0 = 1.0;
  70
  71   if (m_DataLength[m_IteratorDirection] == 1) { //Required by mirror boundaries
  72     return false;
  73   }
  74
  75   // Compute overall gain
  76   for (int k = 0; k < m_NumberOfPoles; k++) {
  77     // Note for cubic splines lambda = 6
  78     c0 = c0 * (1.0 - m_SplinePoles[k]) * (1.0 - 1.0 / m_SplinePoles[k]);
  79   }
  80
  81   // apply the gain
  82   for (unsigned int n = 0; n < m_DataLength[m_IteratorDirection]; n++) {
  83     m_Scratch[n] *= c0;
  84   }
  85
  86   // loop over all poles
  87   for (int k = 0; k < m_NumberOfPoles; k++) {
  88     // causal initialization
  89     this->SetInitialCausalCoefficient(m_SplinePoles[k]);
  90     // causal recursion
  91     for (unsigned int n = 1; n < m_DataLength[m_IteratorDirection]; n++) {
  92       m_Scratch[n] += m_SplinePoles[k] * m_Scratch[n - 1];
  93     }
  94
  95     // anticausal initialization
  96     this->SetInitialAntiCausalCoefficient(m_SplinePoles[k]);
  97     // anticausal recursion
  98     for ( int n = m_DataLength[m_IteratorDirection] - 2; 0 <= n; n--) {
  99       m_Scratch[n] = m_SplinePoles[k] * (m_Scratch[n + 1] - m_Scratch[n]);
 100     }
 101   }
 102   return true;
 103
 104 }
 105
 106
 107 template <class TInputImage, class TOutputImage>
 108 void
 109 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 110 ::SetSplineOrder(unsigned int SplineOrder)
 111 {
 112   if (SplineOrder == m_SplineOrder) {
 113     return;
 114   }
 115   m_SplineOrder = SplineOrder;
 116   this->SetPoles();
 117   this->Modified();
 118
 119 }
 120
 121 //JV
 122 template <class TInputImage, class TOutputImage>
 123 void
 124 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 125 ::SetSplineOrders(SizeType SplineOrders)
 126 {
 127   m_SplineOrders=SplineOrders;
 128 }
 129
 130 template <class TInputImage, class TOutputImage>
 131 void
 132 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 133 ::SetPoles()
 134 {
 135   /* See Unser, 1997. Part II, Table I for Pole values */
 136   // See also, Handbook of Medical Imaging, Processing and Analysis, Ed. Isaac N. Bankman,
 137   //  2000, pg. 416.
 138   switch (m_SplineOrder) {
 139   case 3:
 140     m_NumberOfPoles = 1;
 141     m_SplinePoles[0] = std::sqrt(3.0) - 2.0;
 142     break;
 143   case 0:
 144     m_NumberOfPoles = 0;
 145     break;
 146   case 1:
 147     m_NumberOfPoles = 0;
 148     break;
 149   case 2:
 150     m_NumberOfPoles = 1;
 151     m_SplinePoles[0] = std::sqrt(8.0) - 3.0;
 152     break;
 153   case 4:
 154     m_NumberOfPoles = 2;
 155     m_SplinePoles[0] = std::sqrt(664.0 - std::sqrt(438976.0)) + std::sqrt(304.0) - 19.0;
 156     m_SplinePoles[1] = std::sqrt(664.0 + std::sqrt(438976.0)) - std::sqrt(304.0) - 19.0;
 157     break;
 158   case 5:
 159     m_NumberOfPoles = 2;
 160     m_SplinePoles[0] = std::sqrt(135.0 / 2.0 - std::sqrt(17745.0 / 4.0)) + std::sqrt(105.0 / 4.0)
 161                        - 13.0 / 2.0;
 162     m_SplinePoles[1] = std::sqrt(135.0 / 2.0 + std::sqrt(17745.0 / 4.0)) - std::sqrt(105.0 / 4.0)
 163                        - 13.0 / 2.0;
 164     break;
 165   default:
 166     // SplineOrder not implemented yet.
 167     itk::ExceptionObject err(__FILE__, __LINE__);
 168     err.SetLocation( ITK_LOCATION);
 169     err.SetDescription( "SplineOrder must be between 0 and 5. Requested spline order has not been implemented yet." );
 170     throw err;
 171     break;
 172   }
 173 }
 174
 175
 176 template <class TInputImage, class TOutputImage>
 177 void
 178 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 179 ::SetInitialCausalCoefficient(double z)
 180 {
 181   /* begining InitialCausalCoefficient */
 182   /* See Unser, 1999, Box 2 for explaination */
 183   //JV
 184   itk::Vector<double, VectorDimension> sum;
 185   double  zn, z2n, iz; //sum
 186   unsigned long  horizon;
 187
 188   /* this initialization corresponds to mirror boundaries */
 189   horizon = m_DataLength[m_IteratorDirection];
 190   zn = z;
 191   if (m_Tolerance > 0.0) {
 192     horizon = (long)std::ceil(log(m_Tolerance) / std::log(fabs(z)));
 193   }
 194   if (horizon < m_DataLength[m_IteratorDirection]) {
 195     /* accelerated loop */
 196     sum = m_Scratch[0];   // verify this
 197     for (unsigned int n = 1; n < horizon; n++) {
 198       sum += zn * m_Scratch[n];
 199       zn *= z;
 200     }
 201     m_Scratch[0] = sum;
 202   } else {
 203     /* full loop */
 204     iz = 1.0 / z;
 205     z2n = std::pow(z, (double)(m_DataLength[m_IteratorDirection] - 1L));
 206     sum = m_Scratch[0] + z2n * m_Scratch[m_DataLength[m_IteratorDirection] - 1L];
 207     z2n *= z2n * iz;
 208     for (unsigned int n = 1; n <= (m_DataLength[m_IteratorDirection] - 2); n++) {
 209       sum += (zn + z2n) * m_Scratch[n];
 210       zn *= z;
 211       z2n *= iz;
 212     }
 213     m_Scratch[0] = sum / (1.0 - zn * zn);
 214   }
 215 }
 216
 217
 218 template <class TInputImage, class TOutputImage>
 219 void
 220 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 221 ::SetInitialAntiCausalCoefficient(double z)
 222 {
 223   // this initialization corresponds to mirror boundaries
 224   /* See Unser, 1999, Box 2 for explaination */
 225   //  Also see erratum at http://bigwww.epfl.ch/publications/unser9902.html
 226   m_Scratch[m_DataLength[m_IteratorDirection] - 1] =
 227     (z / (z * z - 1.0)) *
 228     (z * m_Scratch[m_DataLength[m_IteratorDirection] - 2] + m_Scratch[m_DataLength[m_IteratorDirection] - 1]);
 229 }
 230
 231
 232 template <class TInputImage, class TOutputImage>
 233 void
 234 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 235 ::DataToCoefficientsND()
 236 {
 237   OutputImagePointer output = this->GetOutput();
 238
 239   itk::Size<ImageDimension> size = output->GetBufferedRegion().GetSize();
 240
 241   unsigned int count = output->GetBufferedRegion().GetNumberOfPixels() / size[0] * ImageDimension;
 242
 243   itk::ProgressReporter progress(this, 0, count, 10);
 244
 245   // Initialize coeffient array
 246   this->CopyImageToImage();   // Coefficients are initialized to the input data
 247
 248   for (unsigned int n=0; n < ImageDimension; n++) {
 249     m_IteratorDirection = n;
 250     // Loop through each dimension
 251
 252     //JV Set the correct order by dimension!
 253     SetSplineOrder(m_SplineOrders[n]);
 254
 255     // Initialize iterators
 256     OutputLinearIterator CIterator( output, output->GetBufferedRegion() );
 257     CIterator.SetDirection( m_IteratorDirection );
 258     // For each data vector
 259     while ( !CIterator.IsAtEnd() ) {
 260       // Copy coefficients to scratch
 261       this->CopyCoefficientsToScratch( CIterator );
 262
 263
 264       // Perform 1D BSpline calculations
 265       this->DataToCoefficients1D();
 266
 267       // Copy scratch back to coefficients.
 268       // Brings us back to the end of the line we were working on.
 269       CIterator.GoToBeginOfLine();
 270       this->CopyScratchToCoefficients( CIterator ); // m_Scratch = m_Image;
 271       CIterator.NextLine();
 272       progress.CompletedPixel();
 273     }
 274   }
 275 }
 276
 277
 278 /**
 279  * Copy the input image into the output image
 280  */
 281 template <class TInputImage, class TOutputImage>
 282 void
 283 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 284 ::CopyImageToImage()
 285 {
 286
 287   typedef itk::ImageRegionConstIteratorWithIndex< TInputImage > InputIterator;
 288   typedef itk::ImageRegionIterator< TOutputImage > OutputIterator;
 289   typedef typename TOutputImage::PixelType OutputPixelType;
 290
 291   InputIterator inIt( this->GetInput(), this->GetInput()->GetBufferedRegion() );
 292   OutputIterator outIt( this->GetOutput(), this->GetOutput()->GetBufferedRegion() );
 293
 294   inIt.GoToBegin();
 295   outIt.GoToBegin();
 296   OutputPixelType v;
 297   while ( !outIt.IsAtEnd() ) {
 298     for (unsigned int i=0; i< VectorDimension; i++) {
 299       v[i]= static_cast<typename OutputPixelType::ComponentType>( inIt.Get()[i] );
 300     }
 301     outIt.Set( v );
 302     ++inIt;
 303     ++outIt;
 304   }
 305
 306 }
 307
 308
 309 /**
 310  * Copy the scratch to one line of the output image
 311  */
 312 template <class TInputImage, class TOutputImage>
 313 void
 314 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 315 ::CopyScratchToCoefficients(OutputLinearIterator & Iter)
 316 {
 317   typedef typename TOutputImage::PixelType OutputPixelType;
 318   unsigned long j = 0;
 319   OutputPixelType v;
 320   while ( !Iter.IsAtEndOfLine() ) {
 321     for(unsigned int i=0; i<VectorDimension; i++) v[i]=static_cast<typename OutputPixelType::ComponentType>( m_Scratch[j][i]);
 322     Iter.Set( v );
 323     ++Iter;
 324     ++j;
 325   }
 326
 327 }
 328
 329
 330 /**
 331  * Copy one line of the output image to the scratch
 332  */
 333 template <class TInputImage, class TOutputImage>
 334 void
 335 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 336 ::CopyCoefficientsToScratch(OutputLinearIterator & Iter)
 337 {
 338   unsigned long j = 0;
 339   itk::Vector<double, VectorDimension> v;
 340   while ( !Iter.IsAtEndOfLine() ) {
 341     for(unsigned int i=0; i<VectorDimension; i++)v[i]=static_cast<double>( Iter.Get()[i] );
 342     m_Scratch[j] = v ;
 343     ++Iter;
 344     ++j;
 345   }
 346 }
 347
 348
 349 /**
 350  * GenerateInputRequestedRegion method.
 351  */
 352 template <class TInputImage, class TOutputImage>
 353 void
 354 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 355 ::GenerateInputRequestedRegion()
 356 {
 357   // this filter requires the all of the input image to be in
 358   // the buffer
 359   InputImagePointer  inputPtr = const_cast< TInputImage * > ( this->GetInput() );
 360   if( inputPtr ) {
 361     inputPtr->SetRequestedRegionToLargestPossibleRegion();
 362   }
 363 }
 364
 365
 366 /**
 367  * EnlargeOutputRequestedRegion method.
 368  */
 369 template <class TInputImage, class TOutputImage>
 370 void
 371 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 372 ::EnlargeOutputRequestedRegion( itk::DataObject *output )
 373 {
 374
 375   // this filter requires the all of the output image to be in
 376   // the buffer
 377   TOutputImage *imgData;
 378   imgData = dynamic_cast<TOutputImage*>( output );
 379   if( imgData ) {
 380     imgData->SetRequestedRegionToLargestPossibleRegion();
 381   }
 382
 383 }
 384
 385 /**
 386  * Generate data
 387  */
 388 template <class TInputImage, class TOutputImage>
 389 void
 390 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 391 ::GenerateData()
 392 {
 393   DD("VectorBSplineDecompositionImageFilterWithOBD GenerateData()");
 394   // Allocate scratch memory
 395   InputImageConstPointer inputPtr = this->GetInput();
 396   m_DataLength = inputPtr->GetBufferedRegion().GetSize();
 397
 398   unsigned long maxLength = 0;
 399   for ( unsigned int n = 0; n < ImageDimension; n++ ) {
 400     if ( m_DataLength[n] > maxLength ) {
 401       maxLength = m_DataLength[n];
 402     }
 403   }
 404   m_Scratch.resize( maxLength );
 405
 406   // Allocate memory for output image
 407   OutputImagePointer outputPtr = this->GetOutput();
 408   outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() );
 409   outputPtr->Allocate();
 410
 411   // Calculate actual output
 412   this->DataToCoefficientsND();
 413
 414   // Clean up
 415   m_Scratch.clear();
 416
 417 }
 418
 419
 420 } // namespace clitk
 421
 422 #endif