X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=itk%2FclitkVectorBSplineDecompositionImageFilter.txx;h=6161d22a8ff9e83390cc5bb4cd1a2b17c170062a;hb=573d80d0f7a17607d2ee883c21c940c0ba020282;hp=a6d99e35093e00df752708b4b65d2d3cf48aab49;hpb=0083c3fb2c66812489631c7551709d121de51625;p=clitk.git diff --git a/itk/clitkVectorBSplineDecompositionImageFilter.txx b/itk/clitkVectorBSplineDecompositionImageFilter.txx index a6d99e3..6161d22 100644 --- a/itk/clitkVectorBSplineDecompositionImageFilter.txx +++ b/itk/clitkVectorBSplineDecompositionImageFilter.txx @@ -1,3 +1,20 @@ +/*========================================================================= + 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 _clitkVectorBSplineDecompositionImageFilter_txx #define _clitkVectorBSplineDecompositionImageFilter_txx #include "clitkVectorBSplineDecompositionImageFilter.h" @@ -31,7 +48,7 @@ template void VectorBSplineDecompositionImageFilter ::PrintSelf( - std::ostream& os, + std::ostream& os, itk::Indent indent) const { Superclass::PrintSelf( os, indent ); @@ -44,50 +61,44 @@ template bool VectorBSplineDecompositionImageFilter ::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; } @@ -98,10 +109,9 @@ void VectorBSplineDecompositionImageFilter ::SetSplineOrder(unsigned int SplineOrder) { - if (SplineOrder == m_SplineOrder) - { + if (SplineOrder == m_SplineOrder) { return; - } + } m_SplineOrder = SplineOrder; this->SetPoles(); this->Modified(); @@ -115,44 +125,43 @@ VectorBSplineDecompositionImageFilter ::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. - itk::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. + itk::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; + } } @@ -171,34 +180,29 @@ VectorBSplineDecompositionImageFilter /* 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); } } @@ -208,11 +212,11 @@ void VectorBSplineDecompositionImageFilter ::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]); } @@ -233,8 +237,7 @@ VectorBSplineDecompositionImageFilter // 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 @@ -242,23 +245,22 @@ VectorBSplineDecompositionImageFilter 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(); - } } + } } @@ -281,17 +283,15 @@ VectorBSplineDecompositionImageFilter inIt = inIt.Begin(); outIt = outIt.Begin(); OutputPixelType v; - while ( !outIt.IsAtEnd() ) - { - for (unsigned int i=0; i< VectorDimension;i++) - { - v[i]= static_cast( inIt.Get()[i] ); - } - outIt.Set( v ); - ++inIt; - ++outIt; + while ( !outIt.IsAtEnd() ) { + for (unsigned int i=0; i< VectorDimension; i++) { + v[i]= static_cast( inIt.Get()[i] ); } - + outIt.Set( v ); + ++inIt; + ++outIt; + } + } @@ -306,13 +306,12 @@ VectorBSplineDecompositionImageFilter typedef typename TOutputImage::PixelType OutputPixelType; unsigned long j = 0; OutputPixelType v; - while ( !Iter.IsAtEndOfLine() ) - { - for(unsigned int i=0; i( m_Scratch[j][i]); + while ( !Iter.IsAtEndOfLine() ) { + for(unsigned int i=0; i( m_Scratch[j][i]); Iter.Set( v ); ++Iter; ++j; - } + } } @@ -327,13 +326,12 @@ VectorBSplineDecompositionImageFilter { unsigned long j = 0; itk::Vector v; - while ( !Iter.IsAtEndOfLine() ) - { - for(unsigned int i=0; i( Iter.Get()[i] ); + while ( !Iter.IsAtEndOfLine() ) { + for(unsigned int i=0; i( Iter.Get()[i] ); m_Scratch[j] = v ; ++Iter; ++j; - } + } } @@ -348,10 +346,9 @@ VectorBSplineDecompositionImageFilter // 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(); - } + } } @@ -368,10 +365,9 @@ VectorBSplineDecompositionImageFilter // the buffer TOutputImage *imgData; imgData = dynamic_cast( output ); - if( imgData ) - { + if( imgData ) { imgData->SetRequestedRegionToLargestPossibleRegion(); - } + } } @@ -389,13 +385,11 @@ VectorBSplineDecompositionImageFilter 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