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://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 _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     {
  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 VectorBSplineDecompositionImageFilterWithOBD<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 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 132 ::SetSplineOrders(SizeType SplineOrders)
 133 {
 134   m_SplineOrders=SplineOrders;
 135 }
 136
 137 template <class TInputImage, class TOutputImage>
 138 void
 139 VectorBSplineDecompositionImageFilterWithOBD<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] = vcl_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] = vcl_sqrt(8.0) - 3.0;
 160       break;
 161     case 4:
 162       m_NumberOfPoles = 2;
 163       m_SplinePoles[0] = vcl_sqrt(664.0 - vcl_sqrt(438976.0)) + vcl_sqrt(304.0) - 19.0;
 164       m_SplinePoles[1] = vcl_sqrt(664.0 + vcl_sqrt(438976.0)) - vcl_sqrt(304.0) - 19.0;
 165       break;
 166     case 5:
 167       m_NumberOfPoles = 2;
 168       m_SplinePoles[0] = vcl_sqrt(135.0 / 2.0 - vcl_sqrt(17745.0 / 4.0)) + vcl_sqrt(105.0 / 4.0)
 169         - 13.0 / 2.0;
 170       m_SplinePoles[1] = vcl_sqrt(135.0 / 2.0 + vcl_sqrt(17745.0 / 4.0)) - vcl_sqrt(105.0 / 4.0)
 171         - 13.0 / 2.0;
 172       break;
 173     default:
 174       // SplineOrder not implemented yet.
 175       itk::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 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 187 ::SetInitialCausalCoefficient(double z)
 188 {
 189   /* begining InitialCausalCoefficient */
 190   /* See Unser, 1999, Box 2 for explaination */
 191   //JV
 192   itk::Vector<double, VectorDimension> sum;
 193   double  zn, z2n, iz; //sum
 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 VectorBSplineDecompositionImageFilterWithOBD<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 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 248 ::DataToCoefficientsND()
 249 {
 250   OutputImagePointer output = this->GetOutput();
 251
 252   itk::Size<ImageDimension> size = output->GetBufferedRegion().GetSize();
 253
 254   unsigned int count = output->GetBufferedRegion().GetNumberOfPixels() / size[0] * ImageDimension;
 255
 256   itk::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
 279       // Perform 1D BSpline calculations
 280       this->DataToCoefficients1D();
 281
 282       // Copy scratch back to coefficients.
 283       // Brings us back to the end of the line we were working on.
 284       CIterator.GoToBeginOfLine();
 285       this->CopyScratchToCoefficients( CIterator ); // m_Scratch = m_Image;
 286       CIterator.NextLine();
 287       progress.CompletedPixel();
 288       }
 289     }
 290 }
 291
 292
 293 /**
 294  * Copy the input image into the output image
 295  */
 296 template <class TInputImage, class TOutputImage>
 297 void
 298 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 299 ::CopyImageToImage()
 300 {
 301
 302   typedef itk::ImageRegionConstIteratorWithIndex< TInputImage > InputIterator;
 303   typedef itk::ImageRegionIterator< TOutputImage > OutputIterator;
 304   typedef typename TOutputImage::PixelType OutputPixelType;
 305
 306   InputIterator inIt( this->GetInput(), this->GetInput()->GetBufferedRegion() );
 307   OutputIterator outIt( this->GetOutput(), this->GetOutput()->GetBufferedRegion() );
 308
 309   inIt = inIt.Begin();
 310   outIt = outIt.Begin();
 311   OutputPixelType v;
 312   while ( !outIt.IsAtEnd() )
 313     {
 314       for (unsigned int i=0; i< VectorDimension;i++)
 315         {
 316           v[i]= static_cast<typename OutputPixelType::ComponentType>( inIt.Get()[i] );
 317         }
 318       outIt.Set( v );
 319       ++inIt;
 320     ++outIt;
 321     }
 322
 323 }
 324
 325
 326 /**
 327  * Copy the scratch to one line of the output image
 328  */
 329 template <class TInputImage, class TOutputImage>
 330 void
 331 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 332 ::CopyScratchToCoefficients(OutputLinearIterator & Iter)
 333 {
 334   typedef typename TOutputImage::PixelType OutputPixelType;
 335   unsigned long j = 0;
 336   OutputPixelType v;
 337   while ( !Iter.IsAtEndOfLine() )
 338     {
 339       for(unsigned int i=0; i<VectorDimension; i++) v[i]=static_cast<typename OutputPixelType::ComponentType>( m_Scratch[j][i]);
 340     Iter.Set( v );
 341     ++Iter;
 342     ++j;
 343     }
 344
 345 }
 346
 347
 348 /**
 349  * Copy one line of the output image to the scratch
 350  */
 351 template <class TInputImage, class TOutputImage>
 352 void
 353 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 354 ::CopyCoefficientsToScratch(OutputLinearIterator & Iter)
 355 {
 356   unsigned long j = 0;
 357   itk::Vector<double, VectorDimension> v;
 358   while ( !Iter.IsAtEndOfLine() )
 359     {
 360       for(unsigned int i=0; i<VectorDimension; i++)v[i]=static_cast<double>( Iter.Get()[i] );
 361     m_Scratch[j] = v ;
 362     ++Iter;
 363     ++j;
 364     }
 365 }
 366
 367
 368 /**
 369  * GenerateInputRequestedRegion method.
 370  */
 371 template <class TInputImage, class TOutputImage>
 372 void
 373 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 374 ::GenerateInputRequestedRegion()
 375 {
 376   // this filter requires the all of the input image to be in
 377   // the buffer
 378   InputImagePointer  inputPtr = const_cast< TInputImage * > ( this->GetInput() );
 379   if( inputPtr )
 380     {
 381     inputPtr->SetRequestedRegionToLargestPossibleRegion();
 382     }
 383 }
 384
 385
 386 /**
 387  * EnlargeOutputRequestedRegion method.
 388  */
 389 template <class TInputImage, class TOutputImage>
 390 void
 391 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 392 ::EnlargeOutputRequestedRegion( itk::DataObject *output )
 393 {
 394
 395   // this filter requires the all of the output image to be in
 396   // the buffer
 397   TOutputImage *imgData;
 398   imgData = dynamic_cast<TOutputImage*>( output );
 399   if( imgData )
 400     {
 401     imgData->SetRequestedRegionToLargestPossibleRegion();
 402     }
 403
 404 }
 405
 406 /**
 407  * Generate data
 408  */
 409 template <class TInputImage, class TOutputImage>
 410 void
 411 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 412 ::GenerateData()
 413 {
 414   DD("VectorBSplineDecompositionImageFilterWithOBD GenerateData()");
 415   // Allocate scratch memory
 416   InputImageConstPointer inputPtr = this->GetInput();
 417   m_DataLength = inputPtr->GetBufferedRegion().GetSize();
 418
 419   unsigned long maxLength = 0;
 420   for ( unsigned int n = 0; n < ImageDimension; n++ )
 421     {
 422     if ( m_DataLength[n] > maxLength )
 423       {
 424       maxLength = m_DataLength[n];
 425       }
 426     }
 427   m_Scratch.resize( maxLength );
 428
 429   // Allocate memory for output image
 430   OutputImagePointer outputPtr = this->GetOutput();
 431   outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() );
 432   outputPtr->Allocate();
 433
 434   // Calculate actual output
 435   this->DataToCoefficientsND();
 436
 437   // Clean up
 438   m_Scratch.clear();
 439
 440 }
 441
 442
 443 } // namespace clitk
 444
 445 #endif