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://oncora1.lyon.fnclcc.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   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 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
  50 ::PrintSelf(
  51   std::ostream& os,
  52   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 BSplineDecompositionImageFilterWithOBD<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 BSplineDecompositionImageFilterWithOBD<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 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 125 ::SetSplineOrders(SizeType SplineOrders)
 126 {
 127   m_SplineOrders=SplineOrders;
 128 }
 129
 130 template <class TInputImage, class TOutputImage>
 131 void
 132 BSplineDecompositionImageFilterWithOBD<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
 140   case 3:
 141     m_NumberOfPoles = 1;
 142     m_SplinePoles[0] = vcl_sqrt(3.0) - 2.0;
 143     break;
 144   case 0:
 145     m_NumberOfPoles = 0;
 146     break;
 147   case 1:
 148     m_NumberOfPoles = 0;
 149     break;
 150   case 2:
 151     m_NumberOfPoles = 1;
 152     m_SplinePoles[0] = vcl_sqrt(8.0) - 3.0;
 153     break;
 154   case 4:
 155     m_NumberOfPoles = 2;
 156     m_SplinePoles[0] = vcl_sqrt(664.0 - vcl_sqrt(438976.0)) + vcl_sqrt(304.0) - 19.0;
 157     m_SplinePoles[1] = vcl_sqrt(664.0 + vcl_sqrt(438976.0)) - vcl_sqrt(304.0) - 19.0;
 158     break;
 159   case 5:
 160     m_NumberOfPoles = 2;
 161     m_SplinePoles[0] = vcl_sqrt(135.0 / 2.0 - vcl_sqrt(17745.0 / 4.0)) + vcl_sqrt(105.0 / 4.0)
 162                        - 13.0 / 2.0;
 163     m_SplinePoles[1] = vcl_sqrt(135.0 / 2.0 + vcl_sqrt(17745.0 / 4.0)) - vcl_sqrt(105.0 / 4.0)
 164                        - 13.0 / 2.0;
 165     break;
 166   default:
 167     // SplineOrder not implemented yet.
 168     ExceptionObject err(__FILE__, __LINE__);
 169     err.SetLocation( ITK_LOCATION);
 170     err.SetDescription( "SplineOrder must be between 0 and 5. Requested spline order has not been implemented yet." );
 171     throw err;
 172     break;
 173   }
 174 }
 175
 176
 177 template <class TInputImage, class TOutputImage>
 178 void
 179 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 180 ::SetInitialCausalCoefficient(double z)
 181 {
 182   /* begining InitialCausalCoefficient */
 183   /* See Unser, 1999, Box 2 for explaination */
 184
 185   double  sum, zn, z2n, iz;
 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)vcl_ceil(log(m_Tolerance) / vcl_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 = vcl_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 BSplineDecompositionImageFilterWithOBD<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 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 235 ::DataToCoefficientsND()
 236 {
 237   OutputImagePointer output = this->GetOutput();
 238
 239   Size<ImageDimension> size = output->GetBufferedRegion().GetSize();
 240
 241   unsigned int count = output->GetBufferedRegion().GetNumberOfPixels() / size[0] * ImageDimension;
 242
 243   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       // Perform 1D BSpline calculations
 264       this->DataToCoefficients1D();
 265
 266       // Copy scratch back to coefficients.
 267       // Brings us back to the end of the line we were working on.
 268       CIterator.GoToBeginOfLine();
 269       this->CopyScratchToCoefficients( CIterator ); // m_Scratch = m_Image;
 270       CIterator.NextLine();
 271       progress.CompletedPixel();
 272     }
 273   }
 274 }
 275
 276
 277 /**
 278  * Copy the input image into the output image
 279  */
 280 template <class TInputImage, class TOutputImage>
 281 void
 282 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 283 ::CopyImageToImage()
 284 {
 285
 286   typedef ImageRegionConstIteratorWithIndex< TInputImage > InputIterator;
 287   typedef ImageRegionIterator< TOutputImage > OutputIterator;
 288   typedef typename TOutputImage::PixelType OutputPixelType;
 289
 290   InputIterator inIt( this->GetInput(), this->GetInput()->GetBufferedRegion() );
 291   OutputIterator outIt( this->GetOutput(), this->GetOutput()->GetBufferedRegion() );
 292
 293   inIt = inIt.Begin();
 294   outIt = outIt.Begin();
 295
 296   while ( !outIt.IsAtEnd() ) {
 297     outIt.Set( static_cast<OutputPixelType>( inIt.Get() ) );
 298     ++inIt;
 299     ++outIt;
 300   }
 301
 302 }
 303
 304
 305 /**
 306  * Copy the scratch to one line of the output image
 307  */
 308 template <class TInputImage, class TOutputImage>
 309 void
 310 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 311 ::CopyScratchToCoefficients(OutputLinearIterator & Iter)
 312 {
 313   typedef typename TOutputImage::PixelType OutputPixelType;
 314   unsigned long j = 0;
 315   while ( !Iter.IsAtEndOfLine() ) {
 316     Iter.Set( static_cast<OutputPixelType>( m_Scratch[j] ) );
 317     ++Iter;
 318     ++j;
 319   }
 320
 321 }
 322
 323
 324 /**
 325  * Copy one line of the output image to the scratch
 326  */
 327 template <class TInputImage, class TOutputImage>
 328 void
 329 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 330 ::CopyCoefficientsToScratch(OutputLinearIterator & Iter)
 331 {
 332   unsigned long j = 0;
 333   while ( !Iter.IsAtEndOfLine() ) {
 334     m_Scratch[j] = static_cast<double>( Iter.Get() ) ;
 335     ++Iter;
 336     ++j;
 337   }
 338 }
 339
 340
 341 /**
 342  * GenerateInputRequestedRegion method.
 343  */
 344 template <class TInputImage, class TOutputImage>
 345 void
 346 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 347 ::GenerateInputRequestedRegion()
 348 {
 349   // this filter requires the all of the input image to be in
 350   // the buffer
 351   InputImagePointer  inputPtr = const_cast< TInputImage * > ( this->GetInput() );
 352   if( inputPtr ) {
 353     inputPtr->SetRequestedRegionToLargestPossibleRegion();
 354   }
 355 }
 356
 357
 358 /**
 359  * EnlargeOutputRequestedRegion method.
 360  */
 361 template <class TInputImage, class TOutputImage>
 362 void
 363 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 364 ::EnlargeOutputRequestedRegion(
 365   DataObject *output )
 366 {
 367
 368   // this filter requires the all of the output image to be in
 369   // the buffer
 370   TOutputImage *imgData;
 371   imgData = dynamic_cast<TOutputImage*>( output );
 372   if( imgData ) {
 373     imgData->SetRequestedRegionToLargestPossibleRegion();
 374   }
 375
 376 }
 377
 378 /**
 379  * Generate data
 380  */
 381 template <class TInputImage, class TOutputImage>
 382 void
 383 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 384 ::GenerateData()
 385 {
 386
 387   // Allocate scratch memory
 388   InputImageConstPointer inputPtr = this->GetInput();
 389   m_DataLength = inputPtr->GetBufferedRegion().GetSize();
 390
 391   unsigned long maxLength = 0;
 392   for ( unsigned int n = 0; n < ImageDimension; n++ ) {
 393     if ( m_DataLength[n] > maxLength ) {
 394       maxLength = m_DataLength[n];
 395     }
 396   }
 397   m_Scratch.resize( maxLength );
 398
 399   // Allocate memory for output image
 400   OutputImagePointer outputPtr = this->GetOutput();
 401   outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() );
 402   outputPtr->Allocate();
 403
 404   // Calculate actual output
 405   this->DataToCoefficientsND();
 406
 407   // Clean up
 408   m_Scratch.clear();
 409
 410 }
 411
 412
 413 } // namespace itk
 414
 415 #endif