/*=========================================================================
Program: vv http://www.creatis.insa-lyon.fr/rio/vv
- Authors belong to:
+ Authors belong to:
- University of LYON http://www.universite-lyon.fr/
- - Léon Bérard cancer center http://oncora1.lyon.fnclcc.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
- 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"
::BSplineDecompositionImageFilterWithOBD()
{
m_SplineOrder = 0;
+ for(unsigned int i=0; i<ImageDimension; i++)
+ m_SplineOrders[i]=3;
int SplineOrder = 3;
m_Tolerance = 1e-10; // Need some guidance on this one...what is reasonable?
m_IteratorDirection = 0;
void
BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
::PrintSelf(
- std::ostream& os,
+ std::ostream& os,
Indent indent) const
{
Superclass::PrintSelf( os, indent );
bool
BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
::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;
}
BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
::SetSplineOrder(unsigned int SplineOrder)
{
- if (SplineOrder == m_SplineOrder)
- {
+ if (SplineOrder == m_SplineOrder) {
return;
- }
+ }
m_SplineOrder = SplineOrder;
this->SetPoles();
this->Modified();
BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
::SetSplineOrders(SizeType SplineOrders)
{
+ if (SplineOrders == m_SplineOrders)
+ {
+ return;
+ }
m_SplineOrders=SplineOrders;
+ this->Modified();
}
template <class TInputImage, class TOutputImage>
::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] = std::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] = std::sqrt(8.0) - 3.0;
+ break;
+ case 4:
+ m_NumberOfPoles = 2;
+ m_SplinePoles[0] = std::sqrt(664.0 - std::sqrt(438976.0)) + std::sqrt(304.0) - 19.0;
+ m_SplinePoles[1] = std::sqrt(664.0 + std::sqrt(438976.0)) - std::sqrt(304.0) - 19.0;
+ break;
+ case 5:
+ m_NumberOfPoles = 2;
+ m_SplinePoles[0] = std::sqrt(135.0 / 2.0 - std::sqrt(17745.0 / 4.0)) + std::sqrt(105.0 / 4.0)
+ - 13.0 / 2.0;
+ m_SplinePoles[1] = std::sqrt(135.0 / 2.0 + std::sqrt(17745.0 / 4.0)) - std::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;
+ }
}
/* this initialization corresponds to mirror boundaries */
horizon = m_DataLength[m_IteratorDirection];
zn = z;
- if (m_Tolerance > 0.0)
- {
- horizon = (long)vcl_ceil(log(m_Tolerance) / vcl_log(fabs(z)));
- }
- if (horizon < m_DataLength[m_IteratorDirection])
- {
+ if (m_Tolerance > 0.0) {
+ horizon = (long)std::ceil(log(m_Tolerance) / std::log(fabs(z)));
+ }
+ 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 = std::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);
}
}
BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
::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]);
}
// 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
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();
- }
}
+ }
}
inIt = inIt.Begin();
outIt = outIt.Begin();
- while ( !outIt.IsAtEnd() )
- {
+ while ( !outIt.IsAtEnd() ) {
outIt.Set( static_cast<OutputPixelType>( inIt.Get() ) );
++inIt;
++outIt;
- }
-
+ }
+
}
{
typedef typename TOutputImage::PixelType OutputPixelType;
unsigned long j = 0;
- while ( !Iter.IsAtEndOfLine() )
- {
+ while ( !Iter.IsAtEndOfLine() ) {
Iter.Set( static_cast<OutputPixelType>( m_Scratch[j] ) );
++Iter;
++j;
- }
+ }
}
::CopyCoefficientsToScratch(OutputLinearIterator & Iter)
{
unsigned long j = 0;
- while ( !Iter.IsAtEndOfLine() )
- {
+ while ( !Iter.IsAtEndOfLine() ) {
m_Scratch[j] = static_cast<double>( Iter.Get() ) ;
++Iter;
++j;
- }
+ }
}
// 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();
- }
+ }
}
// the buffer
TOutputImage *imgData;
imgData = dynamic_cast<TOutputImage*>( output );
- if( imgData )
- {
+ if( imgData ) {
imgData->SetRequestedRegionToLargestPossibleRegion();
- }
+ }
}
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