itk/clitkVectorBSplineDecompositionImageFilterWithOBD.txx

   1 #ifndef _clitkVectorBSplineDecompositionImageFilterWithOBD_txx
   2 #define _clitkVectorBSplineDecompositionImageFilterWithOBD_txx
   3 #include "clitkVectorBSplineDecompositionImageFilterWithOBD.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 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
  17 ::VectorBSplineDecompositionImageFilterWithOBD()
  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 VectorBSplineDecompositionImageFilterWithOBD<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 VectorBSplineDecompositionImageFilterWithOBD<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 VectorBSplineDecompositionImageFilterWithOBD<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 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 115 ::SetSplineOrders(SizeType SplineOrders)
 116 {
 117   m_SplineOrders=SplineOrders;
 118 }
 119
 120 template <class TInputImage, class TOutputImage>
 121 void
 122 VectorBSplineDecompositionImageFilterWithOBD<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     case 3:
 131       m_NumberOfPoles = 1;
 132       m_SplinePoles[0] = vcl_sqrt(3.0) - 2.0;
 133       break;
 134     case 0:
 135       m_NumberOfPoles = 0;
 136       break;
 137     case 1:
 138       m_NumberOfPoles = 0;
 139       break;
 140     case 2:
 141       m_NumberOfPoles = 1;
 142       m_SplinePoles[0] = vcl_sqrt(8.0) - 3.0;
 143       break;
 144     case 4:
 145       m_NumberOfPoles = 2;
 146       m_SplinePoles[0] = vcl_sqrt(664.0 - vcl_sqrt(438976.0)) + vcl_sqrt(304.0) - 19.0;
 147       m_SplinePoles[1] = vcl_sqrt(664.0 + vcl_sqrt(438976.0)) - vcl_sqrt(304.0) - 19.0;
 148       break;
 149     case 5:
 150       m_NumberOfPoles = 2;
 151       m_SplinePoles[0] = vcl_sqrt(135.0 / 2.0 - vcl_sqrt(17745.0 / 4.0)) + vcl_sqrt(105.0 / 4.0)
 152         - 13.0 / 2.0;
 153       m_SplinePoles[1] = vcl_sqrt(135.0 / 2.0 + vcl_sqrt(17745.0 / 4.0)) - vcl_sqrt(105.0 / 4.0)
 154         - 13.0 / 2.0;
 155       break;
 156     default:
 157       // SplineOrder not implemented yet.
 158       itk::ExceptionObject err(__FILE__, __LINE__);
 159       err.SetLocation( ITK_LOCATION);
 160       err.SetDescription( "SplineOrder must be between 0 and 5. Requested spline order has not been implemented yet." );
 161       throw err;
 162       break;
 163     }
 164 }
 165
 166
 167 template <class TInputImage, class TOutputImage>
 168 void
 169 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 170 ::SetInitialCausalCoefficient(double z)
 171 {
 172   /* begining InitialCausalCoefficient */
 173   /* See Unser, 1999, Box 2 for explaination */
 174   //JV
 175   itk::Vector<double, VectorDimension> sum;
 176   double  zn, z2n, iz; //sum
 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 VectorBSplineDecompositionImageFilterWithOBD<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 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 231 ::DataToCoefficientsND()
 232 {
 233   OutputImagePointer output = this->GetOutput();
 234
 235   itk::Size<ImageDimension> size = output->GetBufferedRegion().GetSize();
 236
 237   unsigned int count = output->GetBufferedRegion().GetNumberOfPixels() / size[0] * ImageDimension;
 238
 239   itk::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
 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 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 282 ::CopyImageToImage()
 283 {
 284
 285   typedef itk::ImageRegionConstIteratorWithIndex< TInputImage > InputIterator;
 286   typedef itk::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   OutputPixelType v;
 295   while ( !outIt.IsAtEnd() )
 296     {
 297       for (unsigned int i=0; i< VectorDimension;i++)
 298         {
 299           v[i]= static_cast<typename OutputPixelType::ComponentType>( inIt.Get()[i] );
 300         }
 301       outIt.Set( v );
 302       ++inIt;
 303     ++outIt;
 304     }
 305
 306 }
 307
 308
 309 /**
 310  * Copy the scratch to one line of the output image
 311  */
 312 template <class TInputImage, class TOutputImage>
 313 void
 314 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 315 ::CopyScratchToCoefficients(OutputLinearIterator & Iter)
 316 {
 317   typedef typename TOutputImage::PixelType OutputPixelType;
 318   unsigned long j = 0;
 319   OutputPixelType v;
 320   while ( !Iter.IsAtEndOfLine() )
 321     {
 322       for(unsigned int i=0; i<VectorDimension; i++) v[i]=static_cast<typename OutputPixelType::ComponentType>( m_Scratch[j][i]);
 323     Iter.Set( v );
 324     ++Iter;
 325     ++j;
 326     }
 327
 328 }
 329
 330
 331 /**
 332  * Copy one line of the output image to the scratch
 333  */
 334 template <class TInputImage, class TOutputImage>
 335 void
 336 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 337 ::CopyCoefficientsToScratch(OutputLinearIterator & Iter)
 338 {
 339   unsigned long j = 0;
 340   itk::Vector<double, VectorDimension> v;
 341   while ( !Iter.IsAtEndOfLine() )
 342     {
 343       for(unsigned int i=0; i<VectorDimension; i++)v[i]=static_cast<double>( Iter.Get()[i] );
 344     m_Scratch[j] = v ;
 345     ++Iter;
 346     ++j;
 347     }
 348 }
 349
 350
 351 /**
 352  * GenerateInputRequestedRegion method.
 353  */
 354 template <class TInputImage, class TOutputImage>
 355 void
 356 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 357 ::GenerateInputRequestedRegion()
 358 {
 359   // this filter requires the all of the input image to be in
 360   // the buffer
 361   InputImagePointer  inputPtr = const_cast< TInputImage * > ( this->GetInput() );
 362   if( inputPtr )
 363     {
 364     inputPtr->SetRequestedRegionToLargestPossibleRegion();
 365     }
 366 }
 367
 368
 369 /**
 370  * EnlargeOutputRequestedRegion method.
 371  */
 372 template <class TInputImage, class TOutputImage>
 373 void
 374 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 375 ::EnlargeOutputRequestedRegion( itk::DataObject *output )
 376 {
 377
 378   // this filter requires the all of the output image to be in
 379   // the buffer
 380   TOutputImage *imgData;
 381   imgData = dynamic_cast<TOutputImage*>( output );
 382   if( imgData )
 383     {
 384     imgData->SetRequestedRegionToLargestPossibleRegion();
 385     }
 386
 387 }
 388
 389 /**
 390  * Generate data
 391  */
 392 template <class TInputImage, class TOutputImage>
 393 void
 394 VectorBSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
 395 ::GenerateData()
 396 {
 397   DD("VectorBSplineDecompositionImageFilterWithOBD GenerateData()");
 398   // Allocate scratch memory
 399   InputImageConstPointer inputPtr = this->GetInput();
 400   m_DataLength = inputPtr->GetBufferedRegion().GetSize();
 401
 402   unsigned long maxLength = 0;
 403   for ( unsigned int n = 0; n < ImageDimension; n++ )
 404     {
 405     if ( m_DataLength[n] > maxLength )
 406       {
 407       maxLength = m_DataLength[n];
 408       }
 409     }
 410   m_Scratch.resize( maxLength );
 411
 412   // Allocate memory for output image
 413   OutputImagePointer outputPtr = this->GetOutput();
 414   outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() );
 415   outputPtr->Allocate();
 416
 417   // Calculate actual output
 418   this->DataToCoefficientsND();
 419
 420   // Clean up
 421   m_Scratch.clear();
 422
 423 }
 424
 425
 426 } // namespace clitk
 427
 428 #endif