itk/clitkVectorBSplineDecompositionImageFilter.txx

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