itk/clitkVectorBSplineDecompositionImageFilterWithOBD.txx

   1 #ifndef _clitkVectorBSplineDecompositionImageFilterWithOBD_txx
   2 #define _clitkVectorBSplineDecompositionImageFilterWithOBD_txx
   3
   4 #include "clitkVectorBSplineDecompositionImageFilterWithOBD.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 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
  18 ::VectorBSplineDecompositionImageFilterWithOBD()
  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 VectorBSplineDecompositionImageFilterWithOBD<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 VectorBSplineDecompositionImageFilterWithOBD<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 VectorBSplineDecompositionImageFilterWithOBD<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 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 116 ::SetSplineOrders(SizeType SplineOrders)
 117 {
 118   m_SplineOrders=SplineOrders;
 119 }
 120
 121 template <class TInputImage, class TOutputImage>
 122 void
 123 VectorBSplineDecompositionImageFilterWithOBD<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     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       itk::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 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 171 ::SetInitialCausalCoefficient(double z)
 172 {
 173   /* begining InitialCausalCoefficient */
 174   /* See Unser, 1999, Box 2 for explaination */
 175   //JV
 176   itk::Vector<double, VectorDimension> sum;
 177   double  zn, z2n, iz; //sum
 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 VectorBSplineDecompositionImageFilterWithOBD<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 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 232 ::DataToCoefficientsND()
 233 {
 234   OutputImagePointer output = this->GetOutput();
 235
 236   itk::Size<ImageDimension> size = output->GetBufferedRegion().GetSize();
 237
 238   unsigned int count = output->GetBufferedRegion().GetNumberOfPixels() / size[0] * ImageDimension;
 239
 240   itk::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
 263       // Perform 1D BSpline calculations
 264       this->DataToCoefficients1D();
 265
 266       // Copy scratch back to coefficients.
 267       // Brings us back to the end of the line we were working on.
 268       CIterator.GoToBeginOfLine();
 269       this->CopyScratchToCoefficients( CIterator ); // m_Scratch = m_Image;
 270       CIterator.NextLine();
 271       progress.CompletedPixel();
 272       }
 273     }
 274 }
 275
 276
 277 /**
 278  * Copy the input image into the output image
 279  */
 280 template <class TInputImage, class TOutputImage>
 281 void
 282 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 283 ::CopyImageToImage()
 284 {
 285
 286   typedef itk::ImageRegionConstIteratorWithIndex< TInputImage > InputIterator;
 287   typedef itk::ImageRegionIterator< TOutputImage > OutputIterator;
 288   typedef typename TOutputImage::PixelType OutputPixelType;
 289
 290   InputIterator inIt( this->GetInput(), this->GetInput()->GetBufferedRegion() );
 291   OutputIterator outIt( this->GetOutput(), this->GetOutput()->GetBufferedRegion() );
 292
 293   inIt = inIt.Begin();
 294   outIt = outIt.Begin();
 295   OutputPixelType v;
 296   while ( !outIt.IsAtEnd() )
 297     {
 298       for (unsigned int i=0; i< VectorDimension;i++)
 299         {
 300           v[i]= static_cast<typename OutputPixelType::ComponentType>( inIt.Get()[i] );
 301         }
 302       outIt.Set( v );
 303       ++inIt;
 304     ++outIt;
 305     }
 306
 307 }
 308
 309
 310 /**
 311  * Copy the scratch to one line of the output image
 312  */
 313 template <class TInputImage, class TOutputImage>
 314 void
 315 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 316 ::CopyScratchToCoefficients(OutputLinearIterator & Iter)
 317 {
 318   typedef typename TOutputImage::PixelType OutputPixelType;
 319   unsigned long j = 0;
 320   OutputPixelType v;
 321   while ( !Iter.IsAtEndOfLine() )
 322     {
 323       for(unsigned int i=0; i<VectorDimension; i++) v[i]=static_cast<typename OutputPixelType::ComponentType>( m_Scratch[j][i]);
 324     Iter.Set( v );
 325     ++Iter;
 326     ++j;
 327     }
 328
 329 }
 330
 331
 332 /**
 333  * Copy one line of the output image to the scratch
 334  */
 335 template <class TInputImage, class TOutputImage>
 336 void
 337 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 338 ::CopyCoefficientsToScratch(OutputLinearIterator & Iter)
 339 {
 340   unsigned long j = 0;
 341   itk::Vector<double, VectorDimension> v;
 342   while ( !Iter.IsAtEndOfLine() )
 343     {
 344       for(unsigned int i=0; i<VectorDimension; i++)v[i]=static_cast<double>( Iter.Get()[i] );
 345     m_Scratch[j] = v ;
 346     ++Iter;
 347     ++j;
 348     }
 349 }
 350
 351
 352 /**
 353  * GenerateInputRequestedRegion method.
 354  */
 355 template <class TInputImage, class TOutputImage>
 356 void
 357 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 358 ::GenerateInputRequestedRegion()
 359 {
 360   // this filter requires the all of the input image to be in
 361   // the buffer
 362   InputImagePointer  inputPtr = const_cast< TInputImage * > ( this->GetInput() );
 363   if( inputPtr )
 364     {
 365     inputPtr->SetRequestedRegionToLargestPossibleRegion();
 366     }
 367 }
 368
 369
 370 /**
 371  * EnlargeOutputRequestedRegion method.
 372  */
 373 template <class TInputImage, class TOutputImage>
 374 void
 375 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 376 ::EnlargeOutputRequestedRegion( itk::DataObject *output )
 377 {
 378
 379   // this filter requires the all of the output image to be in
 380   // the buffer
 381   TOutputImage *imgData;
 382   imgData = dynamic_cast<TOutputImage*>( output );
 383   if( imgData )
 384     {
 385     imgData->SetRequestedRegionToLargestPossibleRegion();
 386     }
 387
 388 }
 389
 390 /**
 391  * Generate data
 392  */
 393 template <class TInputImage, class TOutputImage>
 394 void
 395 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 396 ::GenerateData()
 397 {
 398   DD("VectorBSplineDecompositionImageFilterWithOBD GenerateData()");
 399   // Allocate scratch memory
 400   InputImageConstPointer inputPtr = this->GetInput();
 401   m_DataLength = inputPtr->GetBufferedRegion().GetSize();
 402
 403   unsigned long maxLength = 0;
 404   for ( unsigned int n = 0; n < ImageDimension; n++ )
 405     {
 406     if ( m_DataLength[n] > maxLength )
 407       {
 408       maxLength = m_DataLength[n];
 409       }
 410     }
 411   m_Scratch.resize( maxLength );
 412
 413   // Allocate memory for output image
 414   OutputImagePointer outputPtr = this->GetOutput();
 415   outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() );
 416   outputPtr->Allocate();
 417
 418   // Calculate actual output
 419   this->DataToCoefficientsND();
 420
 421   // Clean up
 422   m_Scratch.clear();
 423
 424 }
 425
 426
 427 } // namespace clitk
 428
 429 #endif