itk/clitkVectorBSplineDecompositionImageFilter.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 _clitkVectorBSplineDecompositionImageFilter_txx
  19 #define _clitkVectorBSplineDecompositionImageFilter_txx
  20 #include "clitkVectorBSplineDecompositionImageFilter.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 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
  34 ::VectorBSplineDecompositionImageFilter()
  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 VectorBSplineDecompositionImageFilter<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 VectorBSplineDecompositionImageFilter<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 VectorBSplineDecompositionImageFilter<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
 129 template <class TInputImage, class TOutputImage>
 130 void
 131 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
 132 ::SetPoles()
 133 {
 134   /* See Unser, 1997. Part II, Table I for Pole values */
 135   // See also, Handbook of Medical Imaging, Processing and Analysis, Ed. Isaac N. Bankman,
 136   //  2000, pg. 416.
 137   switch (m_SplineOrder)
 138     {
 139     case 3:
 140       m_NumberOfPoles = 1;
 141       m_SplinePoles[0] = vcl_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] = vcl_sqrt(8.0) - 3.0;
 152       break;
 153     case 4:
 154       m_NumberOfPoles = 2;
 155       m_SplinePoles[0] = vcl_sqrt(664.0 - vcl_sqrt(438976.0)) + vcl_sqrt(304.0) - 19.0;
 156       m_SplinePoles[1] = vcl_sqrt(664.0 + vcl_sqrt(438976.0)) - vcl_sqrt(304.0) - 19.0;
 157       break;
 158     case 5:
 159       m_NumberOfPoles = 2;
 160       m_SplinePoles[0] = vcl_sqrt(135.0 / 2.0 - vcl_sqrt(17745.0 / 4.0)) + vcl_sqrt(105.0 / 4.0)
 161         - 13.0 / 2.0;
 162       m_SplinePoles[1] = vcl_sqrt(135.0 / 2.0 + vcl_sqrt(17745.0 / 4.0)) - vcl_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 VectorBSplineDecompositionImageFilter<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     {
 193     horizon = (long)vcl_ceil(log(m_Tolerance) / vcl_log(fabs(z)));
 194     }
 195   if (horizon < m_DataLength[m_IteratorDirection])
 196     {
 197     /* accelerated loop */
 198     sum = m_Scratch[0];   // verify this
 199     for (unsigned int n = 1; n < horizon; n++)
 200       {
 201       sum += zn * m_Scratch[n];
 202       zn *= z;
 203       }
 204     m_Scratch[0] = sum;
 205     }
 206   else {
 207   /* full loop */
 208   iz = 1.0 / z;
 209   z2n = vcl_pow(z, (double)(m_DataLength[m_IteratorDirection] - 1L));
 210   sum = m_Scratch[0] + z2n * m_Scratch[m_DataLength[m_IteratorDirection] - 1L];
 211   z2n *= z2n * iz;
 212   for (unsigned int n = 1; n <= (m_DataLength[m_IteratorDirection] - 2); n++)
 213     {
 214     sum += (zn + z2n) * m_Scratch[n];
 215     zn *= z;
 216     z2n *= iz;
 217     }
 218   m_Scratch[0] = sum / (1.0 - zn * zn);
 219   }
 220 }
 221
 222
 223 template <class TInputImage, class TOutputImage>
 224 void
 225 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
 226 ::SetInitialAntiCausalCoefficient(double z)
 227 {
 228   // this initialization corresponds to mirror boundaries
 229   /* See Unser, 1999, Box 2 for explaination */
 230   //  Also see erratum at http://bigwww.epfl.ch/publications/unser9902.html
 231   m_Scratch[m_DataLength[m_IteratorDirection] - 1] =
 232     (z / (z * z - 1.0)) *
 233     (z * m_Scratch[m_DataLength[m_IteratorDirection] - 2] + m_Scratch[m_DataLength[m_IteratorDirection] - 1]);
 234 }
 235
 236
 237 template <class TInputImage, class TOutputImage>
 238 void
 239 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
 240 ::DataToCoefficientsND()
 241 {
 242   OutputImagePointer output = this->GetOutput();
 243
 244   itk::Size<ImageDimension> size = output->GetBufferedRegion().GetSize();
 245
 246   unsigned int count = output->GetBufferedRegion().GetNumberOfPixels() / size[0] * ImageDimension;
 247
 248   itk::ProgressReporter progress(this, 0, count, 10);
 249
 250   // Initialize coeffient array
 251   this->CopyImageToImage();   // Coefficients are initialized to the input data
 252
 253   for (unsigned int n=0; n < ImageDimension; n++)
 254     {
 255     m_IteratorDirection = n;
 256     // Loop through each dimension
 257
 258     // Initialize iterators
 259     OutputLinearIterator CIterator( output, output->GetBufferedRegion() );
 260     CIterator.SetDirection( m_IteratorDirection );
 261     // For each data vector
 262     while ( !CIterator.IsAtEnd() )
 263       {
 264       // Copy coefficients to scratch
 265       this->CopyCoefficientsToScratch( CIterator );
 266
 267
 268       // Perform 1D BSpline calculations
 269       this->DataToCoefficients1D();
 270
 271       // Copy scratch back to coefficients.
 272       // Brings us back to the end of the line we were working on.
 273       CIterator.GoToBeginOfLine();
 274       this->CopyScratchToCoefficients( CIterator ); // m_Scratch = m_Image;
 275       CIterator.NextLine();
 276       progress.CompletedPixel();
 277       }
 278     }
 279 }
 280
 281
 282 /**
 283  * Copy the input image into the output image
 284  */
 285 template <class TInputImage, class TOutputImage>
 286 void
 287 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
 288 ::CopyImageToImage()
 289 {
 290
 291   typedef itk::ImageRegionConstIteratorWithIndex< TInputImage > InputIterator;
 292   typedef itk::ImageRegionIterator< TOutputImage > OutputIterator;
 293   typedef typename TOutputImage::PixelType OutputPixelType;
 294
 295   InputIterator inIt( this->GetInput(), this->GetInput()->GetBufferedRegion() );
 296   OutputIterator outIt( this->GetOutput(), this->GetOutput()->GetBufferedRegion() );
 297
 298   inIt = inIt.Begin();
 299   outIt = outIt.Begin();
 300   OutputPixelType v;
 301   while ( !outIt.IsAtEnd() )
 302     {
 303       for (unsigned int i=0; i< VectorDimension;i++)
 304         {
 305           v[i]= static_cast<typename OutputPixelType::ComponentType>( inIt.Get()[i] );
 306         }
 307       outIt.Set( v );
 308       ++inIt;
 309     ++outIt;
 310     }
 311
 312 }
 313
 314
 315 /**
 316  * Copy the scratch to one line of the output image
 317  */
 318 template <class TInputImage, class TOutputImage>
 319 void
 320 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
 321 ::CopyScratchToCoefficients(OutputLinearIterator & Iter)
 322 {
 323   typedef typename TOutputImage::PixelType OutputPixelType;
 324   unsigned long j = 0;
 325   OutputPixelType v;
 326   while ( !Iter.IsAtEndOfLine() )
 327     {
 328       for(unsigned int i=0; i<VectorDimension; i++) v[i]=static_cast<typename OutputPixelType::ComponentType>( m_Scratch[j][i]);
 329     Iter.Set( v );
 330     ++Iter;
 331     ++j;
 332     }
 333
 334 }
 335
 336
 337 /**
 338  * Copy one line of the output image to the scratch
 339  */
 340 template <class TInputImage, class TOutputImage>
 341 void
 342 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
 343 ::CopyCoefficientsToScratch(OutputLinearIterator & Iter)
 344 {
 345   unsigned long j = 0;
 346   itk::Vector<double, VectorDimension> v;
 347   while ( !Iter.IsAtEndOfLine() )
 348     {
 349       for(unsigned int i=0; i<VectorDimension; i++)v[i]=static_cast<double>( Iter.Get()[i] );
 350     m_Scratch[j] = v ;
 351     ++Iter;
 352     ++j;
 353     }
 354 }
 355
 356
 357 /**
 358  * GenerateInputRequestedRegion method.
 359  */
 360 template <class TInputImage, class TOutputImage>
 361 void
 362 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
 363 ::GenerateInputRequestedRegion()
 364 {
 365   // this filter requires the all of the input image to be in
 366   // the buffer
 367   InputImagePointer  inputPtr = const_cast< TInputImage * > ( this->GetInput() );
 368   if( inputPtr )
 369     {
 370     inputPtr->SetRequestedRegionToLargestPossibleRegion();
 371     }
 372 }
 373
 374
 375 /**
 376  * EnlargeOutputRequestedRegion method.
 377  */
 378 template <class TInputImage, class TOutputImage>
 379 void
 380 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
 381 ::EnlargeOutputRequestedRegion( itk::DataObject *output )
 382 {
 383
 384   // this filter requires the all of the output image to be in
 385   // the buffer
 386   TOutputImage *imgData;
 387   imgData = dynamic_cast<TOutputImage*>( output );
 388   if( imgData )
 389     {
 390     imgData->SetRequestedRegionToLargestPossibleRegion();
 391     }
 392
 393 }
 394
 395 /**
 396  * Generate data
 397  */
 398 template <class TInputImage, class TOutputImage>
 399 void
 400 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
 401 ::GenerateData()
 402 {
 403   DD("VectorBSplineDecompositionImageFilter GenerateData()");
 404   // Allocate scratch memory
 405   InputImageConstPointer inputPtr = this->GetInput();
 406   m_DataLength = inputPtr->GetBufferedRegion().GetSize();
 407
 408   unsigned long maxLength = 0;
 409   for ( unsigned int n = 0; n < ImageDimension; n++ )
 410     {
 411     if ( m_DataLength[n] > maxLength )
 412       {
 413       maxLength = m_DataLength[n];
 414       }
 415     }
 416   m_Scratch.resize( maxLength );
 417
 418   // Allocate memory for output image
 419   OutputImagePointer outputPtr = this->GetOutput();
 420   outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() );
 421   outputPtr->Allocate();
 422
 423   // Calculate actual output
 424   this->DataToCoefficientsND();
 425
 426   // Clean up
 427   m_Scratch.clear();
 428
 429 }
 430
 431
 432 } // namespace clitk
 433
 434 #endif