itk/itkBSplineDecompositionImageFilterWithOBD.txx

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