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     {
  73     return false;
  74     }
  75
  76   // Compute overall gain
  77   for (int k = 0; k < m_NumberOfPoles; k++)
  78     {
  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     {
  86     m_Scratch[n] *= c0;
  87     }
  88
  89   // loop over all poles
  90   for (int k = 0; k < m_NumberOfPoles; k++)
  91     {
  92     // causal initialization
  93     this->SetInitialCausalCoefficient(m_SplinePoles[k]);
  94     // causal recursion
  95     for (unsigned int n = 1; n < m_DataLength[m_IteratorDirection]; n++)
  96       {
  97       m_Scratch[n] += m_SplinePoles[k] * m_Scratch[n - 1];
  98       }
  99
 100     // anticausal initialization
 101     this->SetInitialAntiCausalCoefficient(m_SplinePoles[k]);
 102     // anticausal recursion
 103     for ( int n = m_DataLength[m_IteratorDirection] - 2; 0 <= n; n--)
 104       {
 105       m_Scratch[n] = m_SplinePoles[k] * (m_Scratch[n + 1] - m_Scratch[n]);
 106       }
 107     }
 108   return true;
 109
 110 }
 111
 112
 113 template <class TInputImage, class TOutputImage>
 114 void
 115 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 116 ::SetSplineOrder(unsigned int SplineOrder)
 117 {
 118   if (SplineOrder == m_SplineOrder)
 119     {
 120     return;
 121     }
 122   m_SplineOrder = SplineOrder;
 123   this->SetPoles();
 124   this->Modified();
 125
 126 }
 127
 128 //JV
 129 template <class TInputImage, class TOutputImage>
 130 void
 131 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 132 ::SetSplineOrders(SizeType SplineOrders)
 133 {
 134   m_SplineOrders=SplineOrders;
 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
 148     case 3:
 149       m_NumberOfPoles = 1;
 150       m_SplinePoles[0] = vcl_sqrt(3.0) - 2.0;
 151       break;
 152     case 0:
 153       m_NumberOfPoles = 0;
 154       break;
 155     case 1:
 156       m_NumberOfPoles = 0;
 157       break;
 158     case 2:
 159       m_NumberOfPoles = 1;
 160       m_SplinePoles[0] = vcl_sqrt(8.0) - 3.0;
 161       break;
 162     case 4:
 163       m_NumberOfPoles = 2;
 164       m_SplinePoles[0] = vcl_sqrt(664.0 - vcl_sqrt(438976.0)) + vcl_sqrt(304.0) - 19.0;
 165       m_SplinePoles[1] = vcl_sqrt(664.0 + vcl_sqrt(438976.0)) - vcl_sqrt(304.0) - 19.0;
 166       break;
 167     case 5:
 168       m_NumberOfPoles = 2;
 169       m_SplinePoles[0] = vcl_sqrt(135.0 / 2.0 - vcl_sqrt(17745.0 / 4.0)) + vcl_sqrt(105.0 / 4.0)
 170         - 13.0 / 2.0;
 171       m_SplinePoles[1] = vcl_sqrt(135.0 / 2.0 + vcl_sqrt(17745.0 / 4.0)) - vcl_sqrt(105.0 / 4.0)
 172         - 13.0 / 2.0;
 173       break;
 174     default:
 175       // SplineOrder not implemented yet.
 176       ExceptionObject err(__FILE__, __LINE__);
 177       err.SetLocation( ITK_LOCATION);
 178       err.SetDescription( "SplineOrder must be between 0 and 5. Requested spline order has not been implemented yet." );
 179       throw err;
 180       break;
 181     }
 182 }
 183
 184
 185 template <class TInputImage, class TOutputImage>
 186 void
 187 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 188 ::SetInitialCausalCoefficient(double z)
 189 {
 190   /* begining InitialCausalCoefficient */
 191   /* See Unser, 1999, Box 2 for explaination */
 192
 193   double  sum, zn, z2n, iz;
 194   unsigned long  horizon;
 195
 196   /* this initialization corresponds to mirror boundaries */
 197   horizon = m_DataLength[m_IteratorDirection];
 198   zn = z;
 199   if (m_Tolerance > 0.0)
 200     {
 201     horizon = (long)vcl_ceil(log(m_Tolerance) / vcl_log(fabs(z)));
 202     }
 203   if (horizon < m_DataLength[m_IteratorDirection])
 204     {
 205     /* accelerated loop */
 206     sum = m_Scratch[0];   // verify this
 207     for (unsigned int n = 1; n < horizon; n++)
 208       {
 209       sum += zn * m_Scratch[n];
 210       zn *= z;
 211       }
 212     m_Scratch[0] = sum;
 213     }
 214   else {
 215   /* full loop */
 216   iz = 1.0 / z;
 217   z2n = vcl_pow(z, (double)(m_DataLength[m_IteratorDirection] - 1L));
 218   sum = m_Scratch[0] + z2n * m_Scratch[m_DataLength[m_IteratorDirection] - 1L];
 219   z2n *= z2n * iz;
 220   for (unsigned int n = 1; n <= (m_DataLength[m_IteratorDirection] - 2); n++)
 221     {
 222     sum += (zn + z2n) * m_Scratch[n];
 223     zn *= z;
 224     z2n *= iz;
 225     }
 226   m_Scratch[0] = sum / (1.0 - zn * zn);
 227   }
 228 }
 229
 230
 231 template <class TInputImage, class TOutputImage>
 232 void
 233 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 234 ::SetInitialAntiCausalCoefficient(double z)
 235 {
 236   // this initialization corresponds to mirror boundaries
 237   /* See Unser, 1999, Box 2 for explaination */
 238   //  Also see erratum at http://bigwww.epfl.ch/publications/unser9902.html
 239   m_Scratch[m_DataLength[m_IteratorDirection] - 1] =
 240     (z / (z * z - 1.0)) *
 241     (z * m_Scratch[m_DataLength[m_IteratorDirection] - 2] + m_Scratch[m_DataLength[m_IteratorDirection] - 1]);
 242 }
 243
 244
 245 template <class TInputImage, class TOutputImage>
 246 void
 247 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 248 ::DataToCoefficientsND()
 249 {
 250   OutputImagePointer output = this->GetOutput();
 251
 252   Size<ImageDimension> size = output->GetBufferedRegion().GetSize();
 253
 254   unsigned int count = output->GetBufferedRegion().GetNumberOfPixels() / size[0] * ImageDimension;
 255
 256   ProgressReporter progress(this, 0, count, 10);
 257
 258   // Initialize coeffient array
 259   this->CopyImageToImage();   // Coefficients are initialized to the input data
 260
 261   for (unsigned int n=0; n < ImageDimension; n++)
 262     {
 263     m_IteratorDirection = n;
 264     // Loop through each dimension
 265
 266     //JV Set the correct order by dimension!
 267     SetSplineOrder(m_SplineOrders[n]);
 268
 269     // Initialize iterators
 270     OutputLinearIterator CIterator( output, output->GetBufferedRegion() );
 271     CIterator.SetDirection( m_IteratorDirection );
 272     // For each data vector
 273     while ( !CIterator.IsAtEnd() )
 274       {
 275       // Copy coefficients to scratch
 276       this->CopyCoefficientsToScratch( CIterator );
 277
 278       // Perform 1D BSpline calculations
 279       this->DataToCoefficients1D();
 280
 281       // Copy scratch back to coefficients.
 282       // Brings us back to the end of the line we were working on.
 283       CIterator.GoToBeginOfLine();
 284       this->CopyScratchToCoefficients( CIterator ); // m_Scratch = m_Image;
 285       CIterator.NextLine();
 286       progress.CompletedPixel();
 287       }
 288     }
 289 }
 290
 291
 292 /**
 293  * Copy the input image into the output image
 294  */
 295 template <class TInputImage, class TOutputImage>
 296 void
 297 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 298 ::CopyImageToImage()
 299 {
 300
 301   typedef ImageRegionConstIteratorWithIndex< TInputImage > InputIterator;
 302   typedef ImageRegionIterator< TOutputImage > OutputIterator;
 303   typedef typename TOutputImage::PixelType OutputPixelType;
 304
 305   InputIterator inIt( this->GetInput(), this->GetInput()->GetBufferedRegion() );
 306   OutputIterator outIt( this->GetOutput(), this->GetOutput()->GetBufferedRegion() );
 307
 308   inIt = inIt.Begin();
 309   outIt = outIt.Begin();
 310
 311   while ( !outIt.IsAtEnd() )
 312     {
 313     outIt.Set( static_cast<OutputPixelType>( inIt.Get() ) );
 314     ++inIt;
 315     ++outIt;
 316     }
 317
 318 }
 319
 320
 321 /**
 322  * Copy the scratch to one line of the output image
 323  */
 324 template <class TInputImage, class TOutputImage>
 325 void
 326 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 327 ::CopyScratchToCoefficients(OutputLinearIterator & Iter)
 328 {
 329   typedef typename TOutputImage::PixelType OutputPixelType;
 330   unsigned long j = 0;
 331   while ( !Iter.IsAtEndOfLine() )
 332     {
 333     Iter.Set( static_cast<OutputPixelType>( m_Scratch[j] ) );
 334     ++Iter;
 335     ++j;
 336     }
 337
 338 }
 339
 340
 341 /**
 342  * Copy one line of the output image to the scratch
 343  */
 344 template <class TInputImage, class TOutputImage>
 345 void
 346 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 347 ::CopyCoefficientsToScratch(OutputLinearIterator & Iter)
 348 {
 349   unsigned long j = 0;
 350   while ( !Iter.IsAtEndOfLine() )
 351     {
 352     m_Scratch[j] = static_cast<double>( Iter.Get() ) ;
 353     ++Iter;
 354     ++j;
 355     }
 356 }
 357
 358
 359 /**
 360  * GenerateInputRequestedRegion method.
 361  */
 362 template <class TInputImage, class TOutputImage>
 363 void
 364 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 365 ::GenerateInputRequestedRegion()
 366 {
 367   // this filter requires the all of the input image to be in
 368   // the buffer
 369   InputImagePointer  inputPtr = const_cast< TInputImage * > ( this->GetInput() );
 370   if( inputPtr )
 371     {
 372     inputPtr->SetRequestedRegionToLargestPossibleRegion();
 373     }
 374 }
 375
 376
 377 /**
 378  * EnlargeOutputRequestedRegion method.
 379  */
 380 template <class TInputImage, class TOutputImage>
 381 void
 382 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 383 ::EnlargeOutputRequestedRegion(
 384   DataObject *output )
 385 {
 386
 387   // this filter requires the all of the output image to be in
 388   // the buffer
 389   TOutputImage *imgData;
 390   imgData = dynamic_cast<TOutputImage*>( output );
 391   if( imgData )
 392     {
 393     imgData->SetRequestedRegionToLargestPossibleRegion();
 394     }
 395
 396 }
 397
 398 /**
 399  * Generate data
 400  */
 401 template <class TInputImage, class TOutputImage>
 402 void
 403 BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 404 ::GenerateData()
 405 {
 406
 407   // Allocate scratch memory
 408   InputImageConstPointer inputPtr = this->GetInput();
 409   m_DataLength = inputPtr->GetBufferedRegion().GetSize();
 410
 411   unsigned long maxLength = 0;
 412   for ( unsigned int n = 0; n < ImageDimension; n++ )
 413     {
 414     if ( m_DataLength[n] > maxLength )
 415       {
 416       maxLength = m_DataLength[n];
 417       }
 418     }
 419   m_Scratch.resize( maxLength );
 420
 421   // Allocate memory for output image
 422   OutputImagePointer outputPtr = this->GetOutput();
 423   outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() );
 424   outputPtr->Allocate();
 425
 426   // Calculate actual output
 427   this->DataToCoefficientsND();
 428
 429   // Clean up
 430   m_Scratch.clear();
 431
 432 }
 433
 434
 435 } // namespace itk
 436
 437 #endif