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