itk/itkBSplineDecompositionImageFilterWithOBD.txx

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