X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=itk%2FitkBSplineDecompositionImageFilterWithOBD.txx;h=80414092681c333ecf58ee0004c8990b9af85237;hb=595624eb4297e747630105d45017de69733caef2;hp=1ed152bf0e9d372b69cb9f92f8c9dafad5cad92d;hpb=931a42358442f4ee4f314613c991c838d4b4e3b7;p=clitk.git diff --git a/itk/itkBSplineDecompositionImageFilterWithOBD.txx b/itk/itkBSplineDecompositionImageFilterWithOBD.txx index 1ed152b..8041409 100644 --- a/itk/itkBSplineDecompositionImageFilterWithOBD.txx +++ b/itk/itkBSplineDecompositionImageFilterWithOBD.txx @@ -1,6 +1,22 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://www.centreleonberard.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html +===========================================================================**/ #ifndef _itkBSplineDecompositionImageFilterWithOBD_txx #define _itkBSplineDecompositionImageFilterWithOBD_txx - #include "itkBSplineDecompositionImageFilterWithOBD.h" #include "itkImageRegionConstIteratorWithIndex.h" #include "itkImageRegionIterator.h" @@ -18,6 +34,8 @@ BSplineDecompositionImageFilterWithOBD ::BSplineDecompositionImageFilterWithOBD() { m_SplineOrder = 0; + for(unsigned int i=0; i void BSplineDecompositionImageFilterWithOBD ::PrintSelf( - std::ostream& os, + std::ostream& os, Indent indent) const { Superclass::PrintSelf( os, indent ); @@ -45,50 +63,44 @@ template bool BSplineDecompositionImageFilterWithOBD ::DataToCoefficients1D() -{ +{ - // See Unser, 1993, Part II, Equation 2.5, - // or Unser, 1999, Box 2. for an explaination. + // See Unser, 1993, Part II, Equation 2.5, + // or Unser, 1999, Box 2. for an explaination. - double c0 = 1.0; - - if (m_DataLength[m_IteratorDirection] == 1) //Required by mirror boundaries - { + double c0 = 1.0; + + if (m_DataLength[m_IteratorDirection] == 1) { //Required by mirror boundaries return false; - } + } // Compute overall gain - for (int k = 0; k < m_NumberOfPoles; k++) - { - // Note for cubic splines lambda = 6 + for (int k = 0; k < m_NumberOfPoles; k++) { + // Note for cubic splines lambda = 6 c0 = c0 * (1.0 - m_SplinePoles[k]) * (1.0 - 1.0 / m_SplinePoles[k]); - } + } - // apply the gain - for (unsigned int n = 0; n < m_DataLength[m_IteratorDirection]; n++) - { + // apply the gain + for (unsigned int n = 0; n < m_DataLength[m_IteratorDirection]; n++) { m_Scratch[n] *= c0; - } + } - // loop over all poles - for (int k = 0; k < m_NumberOfPoles; k++) - { - // causal initialization + // loop over all poles + for (int k = 0; k < m_NumberOfPoles; k++) { + // causal initialization this->SetInitialCausalCoefficient(m_SplinePoles[k]); - // causal recursion - for (unsigned int n = 1; n < m_DataLength[m_IteratorDirection]; n++) - { + // causal recursion + for (unsigned int n = 1; n < m_DataLength[m_IteratorDirection]; n++) { m_Scratch[n] += m_SplinePoles[k] * m_Scratch[n - 1]; - } + } - // anticausal initialization + // anticausal initialization this->SetInitialAntiCausalCoefficient(m_SplinePoles[k]); - // anticausal recursion - for ( int n = m_DataLength[m_IteratorDirection] - 2; 0 <= n; n--) - { + // anticausal recursion + for ( int n = m_DataLength[m_IteratorDirection] - 2; 0 <= n; n--) { m_Scratch[n] = m_SplinePoles[k] * (m_Scratch[n + 1] - m_Scratch[n]); - } } + } return true; } @@ -99,10 +111,9 @@ void BSplineDecompositionImageFilterWithOBD ::SetSplineOrder(unsigned int SplineOrder) { - if (SplineOrder == m_SplineOrder) - { + if (SplineOrder == m_SplineOrder) { return; - } + } m_SplineOrder = SplineOrder; this->SetPoles(); this->Modified(); @@ -115,7 +126,12 @@ void BSplineDecompositionImageFilterWithOBD ::SetSplineOrders(SizeType SplineOrders) { + if (SplineOrders == m_SplineOrders) + { + return; + } m_SplineOrders=SplineOrders; + this->Modified(); } template @@ -124,45 +140,44 @@ BSplineDecompositionImageFilterWithOBD ::SetPoles() { /* See Unser, 1997. Part II, Table I for Pole values */ - // See also, Handbook of Medical Imaging, Processing and Analysis, Ed. Isaac N. Bankman, + // See also, Handbook of Medical Imaging, Processing and Analysis, Ed. Isaac N. Bankman, // 2000, pg. 416. - switch (m_SplineOrder) - { - - case 3: - m_NumberOfPoles = 1; - m_SplinePoles[0] = vcl_sqrt(3.0) - 2.0; - break; - case 0: - m_NumberOfPoles = 0; - break; - case 1: - m_NumberOfPoles = 0; - break; - case 2: - m_NumberOfPoles = 1; - m_SplinePoles[0] = vcl_sqrt(8.0) - 3.0; - break; - case 4: - m_NumberOfPoles = 2; - m_SplinePoles[0] = vcl_sqrt(664.0 - vcl_sqrt(438976.0)) + vcl_sqrt(304.0) - 19.0; - m_SplinePoles[1] = vcl_sqrt(664.0 + vcl_sqrt(438976.0)) - vcl_sqrt(304.0) - 19.0; - break; - case 5: - m_NumberOfPoles = 2; - m_SplinePoles[0] = vcl_sqrt(135.0 / 2.0 - vcl_sqrt(17745.0 / 4.0)) + vcl_sqrt(105.0 / 4.0) - - 13.0 / 2.0; - m_SplinePoles[1] = vcl_sqrt(135.0 / 2.0 + vcl_sqrt(17745.0 / 4.0)) - vcl_sqrt(105.0 / 4.0) - - 13.0 / 2.0; - break; - default: - // SplineOrder not implemented yet. - ExceptionObject err(__FILE__, __LINE__); - err.SetLocation( ITK_LOCATION); - err.SetDescription( "SplineOrder must be between 0 and 5. Requested spline order has not been implemented yet." ); - throw err; - break; - } + switch (m_SplineOrder) { + + case 3: + m_NumberOfPoles = 1; + m_SplinePoles[0] = vcl_sqrt(3.0) - 2.0; + break; + case 0: + m_NumberOfPoles = 0; + break; + case 1: + m_NumberOfPoles = 0; + break; + case 2: + m_NumberOfPoles = 1; + m_SplinePoles[0] = vcl_sqrt(8.0) - 3.0; + break; + case 4: + m_NumberOfPoles = 2; + m_SplinePoles[0] = vcl_sqrt(664.0 - vcl_sqrt(438976.0)) + vcl_sqrt(304.0) - 19.0; + m_SplinePoles[1] = vcl_sqrt(664.0 + vcl_sqrt(438976.0)) - vcl_sqrt(304.0) - 19.0; + break; + case 5: + m_NumberOfPoles = 2; + m_SplinePoles[0] = vcl_sqrt(135.0 / 2.0 - vcl_sqrt(17745.0 / 4.0)) + vcl_sqrt(105.0 / 4.0) + - 13.0 / 2.0; + m_SplinePoles[1] = vcl_sqrt(135.0 / 2.0 + vcl_sqrt(17745.0 / 4.0)) - vcl_sqrt(105.0 / 4.0) + - 13.0 / 2.0; + break; + default: + // SplineOrder not implemented yet. + ExceptionObject err(__FILE__, __LINE__); + err.SetLocation( ITK_LOCATION); + err.SetDescription( "SplineOrder must be between 0 and 5. Requested spline order has not been implemented yet." ); + throw err; + break; + } } @@ -180,34 +195,29 @@ BSplineDecompositionImageFilterWithOBD /* this initialization corresponds to mirror boundaries */ horizon = m_DataLength[m_IteratorDirection]; zn = z; - if (m_Tolerance > 0.0) - { + if (m_Tolerance > 0.0) { horizon = (long)vcl_ceil(log(m_Tolerance) / vcl_log(fabs(z))); - } - if (horizon < m_DataLength[m_IteratorDirection]) - { + } + if (horizon < m_DataLength[m_IteratorDirection]) { /* accelerated loop */ sum = m_Scratch[0]; // verify this - for (unsigned int n = 1; n < horizon; n++) - { + for (unsigned int n = 1; n < horizon; n++) { sum += zn * m_Scratch[n]; zn *= z; - } - m_Scratch[0] = sum; } - else { - /* full loop */ - iz = 1.0 / z; - z2n = vcl_pow(z, (double)(m_DataLength[m_IteratorDirection] - 1L)); - sum = m_Scratch[0] + z2n * m_Scratch[m_DataLength[m_IteratorDirection] - 1L]; - z2n *= z2n * iz; - for (unsigned int n = 1; n <= (m_DataLength[m_IteratorDirection] - 2); n++) - { - sum += (zn + z2n) * m_Scratch[n]; - zn *= z; - z2n *= iz; + m_Scratch[0] = sum; + } else { + /* full loop */ + iz = 1.0 / z; + z2n = vcl_pow(z, (double)(m_DataLength[m_IteratorDirection] - 1L)); + sum = m_Scratch[0] + z2n * m_Scratch[m_DataLength[m_IteratorDirection] - 1L]; + z2n *= z2n * iz; + for (unsigned int n = 1; n <= (m_DataLength[m_IteratorDirection] - 2); n++) { + sum += (zn + z2n) * m_Scratch[n]; + zn *= z; + z2n *= iz; } - m_Scratch[0] = sum / (1.0 - zn * zn); + m_Scratch[0] = sum / (1.0 - zn * zn); } } @@ -217,11 +227,11 @@ void BSplineDecompositionImageFilterWithOBD ::SetInitialAntiCausalCoefficient(double z) { - // this initialization corresponds to mirror boundaries + // this initialization corresponds to mirror boundaries /* See Unser, 1999, Box 2 for explaination */ // Also see erratum at http://bigwww.epfl.ch/publications/unser9902.html m_Scratch[m_DataLength[m_IteratorDirection] - 1] = - (z / (z * z - 1.0)) * + (z / (z * z - 1.0)) * (z * m_Scratch[m_DataLength[m_IteratorDirection] - 2] + m_Scratch[m_DataLength[m_IteratorDirection] - 1]); } @@ -242,8 +252,7 @@ BSplineDecompositionImageFilterWithOBD // Initialize coeffient array this->CopyImageToImage(); // Coefficients are initialized to the input data - for (unsigned int n=0; n < ImageDimension; n++) - { + for (unsigned int n=0; n < ImageDimension; n++) { m_IteratorDirection = n; // Loop through each dimension @@ -254,22 +263,21 @@ BSplineDecompositionImageFilterWithOBD OutputLinearIterator CIterator( output, output->GetBufferedRegion() ); CIterator.SetDirection( m_IteratorDirection ); // For each data vector - while ( !CIterator.IsAtEnd() ) - { + while ( !CIterator.IsAtEnd() ) { // Copy coefficients to scratch this->CopyCoefficientsToScratch( CIterator ); // Perform 1D BSpline calculations this->DataToCoefficients1D(); - + // Copy scratch back to coefficients. // Brings us back to the end of the line we were working on. CIterator.GoToBeginOfLine(); this->CopyScratchToCoefficients( CIterator ); // m_Scratch = m_Image; CIterator.NextLine(); progress.CompletedPixel(); - } } + } } @@ -292,13 +300,12 @@ BSplineDecompositionImageFilterWithOBD inIt = inIt.Begin(); outIt = outIt.Begin(); - while ( !outIt.IsAtEnd() ) - { + while ( !outIt.IsAtEnd() ) { outIt.Set( static_cast( inIt.Get() ) ); ++inIt; ++outIt; - } - + } + } @@ -312,12 +319,11 @@ BSplineDecompositionImageFilterWithOBD { typedef typename TOutputImage::PixelType OutputPixelType; unsigned long j = 0; - while ( !Iter.IsAtEndOfLine() ) - { + while ( !Iter.IsAtEndOfLine() ) { Iter.Set( static_cast( m_Scratch[j] ) ); ++Iter; ++j; - } + } } @@ -331,12 +337,11 @@ BSplineDecompositionImageFilterWithOBD ::CopyCoefficientsToScratch(OutputLinearIterator & Iter) { unsigned long j = 0; - while ( !Iter.IsAtEndOfLine() ) - { + while ( !Iter.IsAtEndOfLine() ) { m_Scratch[j] = static_cast( Iter.Get() ) ; ++Iter; ++j; - } + } } @@ -351,10 +356,9 @@ BSplineDecompositionImageFilterWithOBD // this filter requires the all of the input image to be in // the buffer InputImagePointer inputPtr = const_cast< TInputImage * > ( this->GetInput() ); - if( inputPtr ) - { + if( inputPtr ) { inputPtr->SetRequestedRegionToLargestPossibleRegion(); - } + } } @@ -372,10 +376,9 @@ BSplineDecompositionImageFilterWithOBD // the buffer TOutputImage *imgData; imgData = dynamic_cast( output ); - if( imgData ) - { + if( imgData ) { imgData->SetRequestedRegionToLargestPossibleRegion(); - } + } } @@ -393,13 +396,11 @@ BSplineDecompositionImageFilterWithOBD m_DataLength = inputPtr->GetBufferedRegion().GetSize(); unsigned long maxLength = 0; - for ( unsigned int n = 0; n < ImageDimension; n++ ) - { - if ( m_DataLength[n] > maxLength ) - { + for ( unsigned int n = 0; n < ImageDimension; n++ ) { + if ( m_DataLength[n] > maxLength ) { maxLength = m_DataLength[n]; - } } + } m_Scratch.resize( maxLength ); // Allocate memory for output image