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