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 "clitkDD.h"
  22 #include "itkImageRegionConstIteratorWithIndex.h"
  23 #include "itkImageRegionIterator.h"
  24 #include "itkProgressReporter.h"
  25 #include "itkVector.h"
  26
  27 namespace clitk
  28 {
  29
  30 /**
  31  * Constructor
  32  */
  33 template <class TInputImage, class TOutputImage>
  34 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
  35 ::VectorBSplineDecompositionImageFilter()
  36 {
  37   m_SplineOrder = 0;
  38   int SplineOrder = 3;
  39   m_Tolerance = 1e-10;   // Need some guidance on this one...what is reasonable?
  40   m_IteratorDirection = 0;
  41   this->SetSplineOrder(SplineOrder);
  42 }
  43
  44
  45 /**
  46  * Standard "PrintSelf" method
  47  */
  48 template <class TInputImage, class TOutputImage>
  49 void
  50 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
  51 ::PrintSelf(
  52   std::ostream& os,
  53   itk::Indent indent) const
  54 {
  55   Superclass::PrintSelf( os, indent );
  56   os << indent << "Spline Order: " << m_SplineOrder << std::endl;
  57
  58 }
  59
  60
  61 template <class TInputImage, class TOutputImage>
  62 bool
  63 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
  64 ::DataToCoefficients1D()
  65 {
  66
  67   // See Unser, 1993, Part II, Equation 2.5,
  68   //   or Unser, 1999, Box 2. for an explaination.
  69
  70   double c0 = 1.0;
  71
  72   if (m_DataLength[m_IteratorDirection] == 1) { //Required by mirror boundaries
  73     return false;
  74   }
  75
  76   // Compute overall gain
  77   for (int k = 0; k < m_NumberOfPoles; k++) {
  78     // Note for cubic splines lambda = 6
  79     c0 = c0 * (1.0 - m_SplinePoles[k]) * (1.0 - 1.0 / m_SplinePoles[k]);
  80   }
  81
  82   // apply the gain
  83   for (unsigned int n = 0; n < m_DataLength[m_IteratorDirection]; n++) {
  84     m_Scratch[n] *= c0;
  85   }
  86
  87   // loop over all poles
  88   for (int k = 0; k < m_NumberOfPoles; k++) {
  89     // causal initialization
  90     this->SetInitialCausalCoefficient(m_SplinePoles[k]);
  91     // causal recursion
  92     for (unsigned int n = 1; n < m_DataLength[m_IteratorDirection]; n++) {
  93       m_Scratch[n] += m_SplinePoles[k] * m_Scratch[n - 1];
  94     }
  95
  96     // anticausal initialization
  97     this->SetInitialAntiCausalCoefficient(m_SplinePoles[k]);
  98     // anticausal recursion
  99     for ( int n = m_DataLength[m_IteratorDirection] - 2; 0 <= n; n--) {
 100       m_Scratch[n] = m_SplinePoles[k] * (m_Scratch[n + 1] - m_Scratch[n]);
 101     }
 102   }
 103   return true;
 104
 105 }
 106
 107
 108 template <class TInputImage, class TOutputImage>
 109 void
 110 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
 111 ::SetSplineOrder(unsigned int SplineOrder)
 112 {
 113   if (SplineOrder == m_SplineOrder) {
 114     return;
 115   }
 116   m_SplineOrder = SplineOrder;
 117   this->SetPoles();
 118   this->Modified();
 119
 120 }
 121
 122
 123 template <class TInputImage, class TOutputImage>
 124 void
 125 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
 126 ::SetPoles()
 127 {
 128   /* See Unser, 1997. Part II, Table I for Pole values */
 129   // See also, Handbook of Medical Imaging, Processing and Analysis, Ed. Isaac N. Bankman,
 130   //  2000, pg. 416.
 131   switch (m_SplineOrder) {
 132   case 3:
 133     m_NumberOfPoles = 1;
 134     m_SplinePoles[0] = std::sqrt(3.0) - 2.0;
 135     break;
 136   case 0:
 137     m_NumberOfPoles = 0;
 138     break;
 139   case 1:
 140     m_NumberOfPoles = 0;
 141     break;
 142   case 2:
 143     m_NumberOfPoles = 1;
 144     m_SplinePoles[0] = std::sqrt(8.0) - 3.0;
 145     break;
 146   case 4:
 147     m_NumberOfPoles = 2;
 148     m_SplinePoles[0] = std::sqrt(664.0 - std::sqrt(438976.0)) + std::sqrt(304.0) - 19.0;
 149     m_SplinePoles[1] = std::sqrt(664.0 + std::sqrt(438976.0)) - std::sqrt(304.0) - 19.0;
 150     break;
 151   case 5:
 152     m_NumberOfPoles = 2;
 153     m_SplinePoles[0] = std::sqrt(135.0 / 2.0 - std::sqrt(17745.0 / 4.0)) + std::sqrt(105.0 / 4.0)
 154                        - 13.0 / 2.0;
 155     m_SplinePoles[1] = std::sqrt(135.0 / 2.0 + std::sqrt(17745.0 / 4.0)) - std::sqrt(105.0 / 4.0)
 156                        - 13.0 / 2.0;
 157     break;
 158   default:
 159     // SplineOrder not implemented yet.
 160     itk::ExceptionObject err(__FILE__, __LINE__);
 161     err.SetLocation( ITK_LOCATION);
 162     err.SetDescription( "SplineOrder must be between 0 and 5. Requested spline order has not been implemented yet." );
 163     throw err;
 164     break;
 165   }
 166 }
 167
 168
 169 template <class TInputImage, class TOutputImage>
 170 void
 171 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
 172 ::SetInitialCausalCoefficient(double z)
 173 {
 174   /* begining InitialCausalCoefficient */
 175   /* See Unser, 1999, Box 2 for explaination */
 176   //JV
 177   itk::Vector<double, VectorDimension> sum;
 178   double  zn, z2n, iz; //sum
 179   unsigned long  horizon;
 180
 181   /* this initialization corresponds to mirror boundaries */
 182   horizon = m_DataLength[m_IteratorDirection];
 183   zn = z;
 184   if (m_Tolerance > 0.0) {
 185     horizon = (long)std::ceil(log(m_Tolerance) / std::log(fabs(z)));
 186   }
 187   if (horizon < m_DataLength[m_IteratorDirection]) {
 188     /* accelerated loop */
 189     sum = m_Scratch[0];   // verify this
 190     for (unsigned int n = 1; n < horizon; n++) {
 191       sum += zn * m_Scratch[n];
 192       zn *= z;
 193     }
 194     m_Scratch[0] = sum;
 195   } else {
 196     /* full loop */
 197     iz = 1.0 / z;
 198     z2n = std::pow(z, (double)(m_DataLength[m_IteratorDirection] - 1L));
 199     sum = m_Scratch[0] + z2n * m_Scratch[m_DataLength[m_IteratorDirection] - 1L];
 200     z2n *= z2n * iz;
 201     for (unsigned int n = 1; n <= (m_DataLength[m_IteratorDirection] - 2); n++) {
 202       sum += (zn + z2n) * m_Scratch[n];
 203       zn *= z;
 204       z2n *= iz;
 205     }
 206     m_Scratch[0] = sum / (1.0 - zn * zn);
 207   }
 208 }
 209
 210
 211 template <class TInputImage, class TOutputImage>
 212 void
 213 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
 214 ::SetInitialAntiCausalCoefficient(double z)
 215 {
 216   // this initialization corresponds to mirror boundaries
 217   /* See Unser, 1999, Box 2 for explaination */
 218   //  Also see erratum at http://bigwww.epfl.ch/publications/unser9902.html
 219   m_Scratch[m_DataLength[m_IteratorDirection] - 1] =
 220     (z / (z * z - 1.0)) *
 221     (z * m_Scratch[m_DataLength[m_IteratorDirection] - 2] + m_Scratch[m_DataLength[m_IteratorDirection] - 1]);
 222 }
 223
 224
 225 template <class TInputImage, class TOutputImage>
 226 void
 227 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
 228 ::DataToCoefficientsND()
 229 {
 230   OutputImagePointer output = this->GetOutput();
 231
 232   itk::Size<ImageDimension> size = output->GetBufferedRegion().GetSize();
 233
 234   unsigned int count = output->GetBufferedRegion().GetNumberOfPixels() / size[0] * ImageDimension;
 235
 236   itk::ProgressReporter progress(this, 0, count, 10);
 237
 238   // Initialize coeffient array
 239   this->CopyImageToImage();   // Coefficients are initialized to the input data
 240
 241   for (unsigned int n=0; n < ImageDimension; n++) {
 242     m_IteratorDirection = n;
 243     // Loop through each dimension
 244
 245     // Initialize iterators
 246     OutputLinearIterator CIterator( output, output->GetBufferedRegion() );
 247     CIterator.SetDirection( m_IteratorDirection );
 248     // For each data vector
 249     while ( !CIterator.IsAtEnd() ) {
 250       // Copy coefficients to scratch
 251       this->CopyCoefficientsToScratch( CIterator );
 252
 253
 254       // Perform 1D BSpline calculations
 255       this->DataToCoefficients1D();
 256
 257       // Copy scratch back to coefficients.
 258       // Brings us back to the end of the line we were working on.
 259       CIterator.GoToBeginOfLine();
 260       this->CopyScratchToCoefficients( CIterator ); // m_Scratch = m_Image;
 261       CIterator.NextLine();
 262       progress.CompletedPixel();
 263     }
 264   }
 265 }
 266
 267
 268 /**
 269  * Copy the input image into the output image
 270  */
 271 template <class TInputImage, class TOutputImage>
 272 void
 273 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
 274 ::CopyImageToImage()
 275 {
 276
 277   typedef itk::ImageRegionConstIteratorWithIndex< TInputImage > InputIterator;
 278   typedef itk::ImageRegionIterator< TOutputImage > OutputIterator;
 279   typedef typename TOutputImage::PixelType OutputPixelType;
 280
 281   InputIterator inIt( this->GetInput(), this->GetInput()->GetBufferedRegion() );
 282   OutputIterator outIt( this->GetOutput(), this->GetOutput()->GetBufferedRegion() );
 283
 284   inIt.GoToBegin();
 285   outIt.GoToBegin();
 286   OutputPixelType v;
 287   while ( !outIt.IsAtEnd() ) {
 288     for (unsigned int i=0; i< VectorDimension; i++) {
 289       v[i]= static_cast<typename OutputPixelType::ComponentType>( inIt.Get()[i] );
 290     }
 291     outIt.Set( v );
 292     ++inIt;
 293     ++outIt;
 294   }
 295
 296 }
 297
 298
 299 /**
 300  * Copy the scratch to one line of the output image
 301  */
 302 template <class TInputImage, class TOutputImage>
 303 void
 304 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
 305 ::CopyScratchToCoefficients(OutputLinearIterator & Iter)
 306 {
 307   typedef typename TOutputImage::PixelType OutputPixelType;
 308   unsigned long j = 0;
 309   OutputPixelType v;
 310   while ( !Iter.IsAtEndOfLine() ) {
 311     for(unsigned int i=0; i<VectorDimension; i++) v[i]=static_cast<typename OutputPixelType::ComponentType>( m_Scratch[j][i]);
 312     Iter.Set( v );
 313     ++Iter;
 314     ++j;
 315   }
 316
 317 }
 318
 319
 320 /**
 321  * Copy one line of the output image to the scratch
 322  */
 323 template <class TInputImage, class TOutputImage>
 324 void
 325 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
 326 ::CopyCoefficientsToScratch(OutputLinearIterator & Iter)
 327 {
 328   unsigned long j = 0;
 329   itk::Vector<double, VectorDimension> v;
 330   while ( !Iter.IsAtEndOfLine() ) {
 331     for(unsigned int i=0; i<VectorDimension; i++)v[i]=static_cast<double>( Iter.Get()[i] );
 332     m_Scratch[j] = v ;
 333     ++Iter;
 334     ++j;
 335   }
 336 }
 337
 338
 339 /**
 340  * GenerateInputRequestedRegion method.
 341  */
 342 template <class TInputImage, class TOutputImage>
 343 void
 344 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
 345 ::GenerateInputRequestedRegion()
 346 {
 347   // this filter requires the all of the input image to be in
 348   // the buffer
 349   InputImagePointer  inputPtr = const_cast< TInputImage * > ( this->GetInput() );
 350   if( inputPtr ) {
 351     inputPtr->SetRequestedRegionToLargestPossibleRegion();
 352   }
 353 }
 354
 355
 356 /**
 357  * EnlargeOutputRequestedRegion method.
 358  */
 359 template <class TInputImage, class TOutputImage>
 360 void
 361 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
 362 ::EnlargeOutputRequestedRegion( itk::DataObject *output )
 363 {
 364
 365   // this filter requires the all of the output image to be in
 366   // the buffer
 367   TOutputImage *imgData;
 368   imgData = dynamic_cast<TOutputImage*>( output );
 369   if( imgData ) {
 370     imgData->SetRequestedRegionToLargestPossibleRegion();
 371   }
 372
 373 }
 374
 375 /**
 376  * Generate data
 377  */
 378 template <class TInputImage, class TOutputImage>
 379 void
 380 VectorBSplineDecompositionImageFilter<TInputImage, TOutputImage>
 381 ::GenerateData()
 382 {
 383   DD("VectorBSplineDecompositionImageFilter GenerateData()");
 384   // Allocate scratch memory
 385   InputImageConstPointer inputPtr = this->GetInput();
 386   m_DataLength = inputPtr->GetBufferedRegion().GetSize();
 387
 388   unsigned long maxLength = 0;
 389   for ( unsigned int n = 0; n < ImageDimension; n++ ) {
 390     if ( m_DataLength[n] > maxLength ) {
 391       maxLength = m_DataLength[n];
 392     }
 393   }
 394   m_Scratch.resize( maxLength );
 395
 396   // Allocate memory for output image
 397   OutputImagePointer outputPtr = this->GetOutput();
 398   outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() );
 399   outputPtr->Allocate();
 400
 401   // Calculate actual output
 402   this->DataToCoefficientsND();
 403
 404   // Clean up
 405   m_Scratch.clear();
 406
 407 }
 408
 409
 410 } // namespace clitk
 411
 412 #endif