itk/clitkVectorBSplineDecompositionImageFilter.txx

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