X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=itk%2FclitkVectorBSplineInterpolateImageFunction.txx;h=2599b80a0f68bc3dcbf8f8ec3d4aa0a47af913f1;hb=573d80d0f7a17607d2ee883c21c940c0ba020282;hp=d80178bfdbcfa0494527e9eb5765e7ecd6eeb934;hpb=931a42358442f4ee4f314613c991c838d4b4e3b7;p=clitk.git diff --git a/itk/clitkVectorBSplineInterpolateImageFunction.txx b/itk/clitkVectorBSplineInterpolateImageFunction.txx index d80178b..2599b80 100644 --- a/itk/clitkVectorBSplineInterpolateImageFunction.txx +++ b/itk/clitkVectorBSplineInterpolateImageFunction.txx @@ -1,9 +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 _clitkVectorBSplineInterpolateImageFunction_txx #define _clitkVectorBSplineInterpolateImageFunction_txx - -// First, make sure that we include the configuration file. -// This line may be removed once the ThreadSafeTransform gets -// integrated into ITK. #include "itkConfigure.h" // Second, redirect to the optimized version if necessary @@ -48,75 +61,70 @@ template void VectorBSplineInterpolateImageFunction ::PrintSelf( - std::ostream& os, + std::ostream& os, itk::Indent indent) const { Superclass::PrintSelf( os, indent ); os << indent << "Spline Order: " << m_SplineOrder << std::endl; - os << indent << "UseImageDirection = " + os << indent << "UseImageDirection = " << (this->m_UseImageDirection ? "On" : "Off") << std::endl; } template -void +void VectorBSplineInterpolateImageFunction ::SetInputImage(const TImageType * inputData) { - if ( inputData ) - { - - DD("calling decomposition filter"); - m_CoefficientFilter->SetInput(inputData); - - // the Coefficient Filter requires that the spline order and the input data be set. - // TODO: We need to ensure that this is only run once and only after both input and - // spline order have been set. Should we force an update after the - // splineOrder has been set also? - - m_CoefficientFilter->Update(); - m_Coefficients = m_CoefficientFilter->GetOutput(); - - - // Call the Superclass implementation after, in case the filter - // pulls in more of the input image - Superclass::SetInputImage(inputData); - - m_DataLength = inputData->GetBufferedRegion().GetSize(); - - } - else - { - m_Coefficients = NULL; - } + if ( inputData ) { + + DD("calling decomposition filter"); + m_CoefficientFilter->SetInput(inputData); + + // the Coefficient Filter requires that the spline order and the input data be set. + // TODO: We need to ensure that this is only run once and only after both input and + // spline order have been set. Should we force an update after the + // splineOrder has been set also? + + m_CoefficientFilter->Update(); + m_Coefficients = m_CoefficientFilter->GetOutput(); + + + // Call the Superclass implementation after, in case the filter + // pulls in more of the input image + Superclass::SetInputImage(inputData); + + m_DataLength = inputData->GetBufferedRegion().GetSize(); + + } else { + m_Coefficients = NULL; + } } template -void +void VectorBSplineInterpolateImageFunction ::SetSplineOrder(unsigned int SplineOrder) { - if (SplineOrder == m_SplineOrder) - { + if (SplineOrder == m_SplineOrder) { return; - } + } m_SplineOrder = SplineOrder; m_CoefficientFilter->SetSplineOrder( SplineOrder ); //this->SetPoles();/ m_MaxNumberInterpolationPoints = 1; - for (unsigned int n=0; n < ImageDimension; n++) - { + for (unsigned int n=0; n < ImageDimension; n++) { m_MaxNumberInterpolationPoints *= ( m_SplineOrder + 1); - } + } this->GeneratePointsToIndex( ); } template -typename +typename VectorBSplineInterpolateImageFunction ::OutputType VectorBSplineInterpolateImageFunction @@ -124,8 +132,8 @@ VectorBSplineInterpolateImageFunction { vnl_matrix EvaluateIndex(ImageDimension, ( m_SplineOrder + 1 )); - // compute the interpolation indexes - this->DetermineRegionOfSupport(EvaluateIndex, x, m_SplineOrder); + // compute the interpolation indexes + this->DetermineRegionOfSupport(EvaluateIndex, x, m_SplineOrder); // Determine weights vnl_matrix weights(ImageDimension, ( m_SplineOrder + 1 )); @@ -133,59 +141,57 @@ VectorBSplineInterpolateImageFunction // Modify EvaluateIndex at the boundaries using mirror boundary conditions this->ApplyMirrorBoundaryConditions(EvaluateIndex, m_SplineOrder); - + // perform interpolation - //JV - itk::Vector interpolated; - for (unsigned int i=0; i< VectorDimension; i++) interpolated[i]=itk::NumericTraits::Zero; + //JV + itk::Vector interpolated; + for (unsigned int i=0; i< VectorDimension; i++) interpolated[i]=itk::NumericTraits::Zero; IndexType coefficientIndex; // Step through eachpoint in the N-dimensional interpolation cube. - for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++) - { + for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++) { // translate each step into the N-dimensional index. // IndexType pointIndex = PointToIndex( p ); double w = 1.0; - for (unsigned int n = 0; n < ImageDimension; n++ ) - { + for (unsigned int n = 0; n < ImageDimension; n++ ) { - w *= weights[n][ m_PointsToIndex[p][n] ]; - coefficientIndex[n] = EvaluateIndex[n][m_PointsToIndex[p][n]]; // Build up ND index for coefficients. - } + w *= weights[n][ m_PointsToIndex[p][n] ]; + coefficientIndex[n] = EvaluateIndex[n][m_PointsToIndex[p][n]]; // Build up ND index for coefficients. + } // Convert our step p to the appropriate point in ND space in the // m_Coefficients cube. //JV shouldn't be necessary for (unsigned int i=0; iGetPixel(coefficientIndex)[i]; - } - -/* double interpolated = 0.0; - IndexType coefficientIndex; - // Step through eachpoint in the N-dimensional interpolation cube. - for (unsigned int sp = 0; sp <= m_SplineOrder; sp++) - { - for (unsigned int sp1=0; sp1 <= m_SplineOrder; sp1++) - { + } - double w = 1.0; - for (unsigned int n1 = 0; n1 < ImageDimension; n1++ ) + /* double interpolated = 0.0; + IndexType coefficientIndex; + // Step through eachpoint in the N-dimensional interpolation cube. + for (unsigned int sp = 0; sp <= m_SplineOrder; sp++) + { + for (unsigned int sp1=0; sp1 <= m_SplineOrder; sp1++) { - w *= weights[n1][ sp1 ]; - coefficientIndex[n1] = EvaluateIndex[n1][sp]; // Build up ND index for coefficients. - } - interpolated += w * m_Coefficients->GetPixel(coefficientIndex); - } - } -*/ + double w = 1.0; + for (unsigned int n1 = 0; n1 < ImageDimension; n1++ ) + { + w *= weights[n1][ sp1 ]; + coefficientIndex[n1] = EvaluateIndex[n1][sp]; // Build up ND index for coefficients. + } + + interpolated += w * m_Coefficients->GetPixel(coefficientIndex); + } + } + */ return(interpolated); - + } template -typename +typename VectorBSplineInterpolateImageFunction :: CovariantVectorType VectorBSplineInterpolateImageFunction @@ -193,7 +199,7 @@ VectorBSplineInterpolateImageFunction { vnl_matrix EvaluateIndex(ImageDimension, ( m_SplineOrder + 1 )); - // compute the interpolation indexes + // compute the interpolation indexes // TODO: Do we need to revisit region of support for the derivatives? this->DetermineRegionOfSupport(EvaluateIndex, x, m_SplineOrder); @@ -206,7 +212,7 @@ VectorBSplineInterpolateImageFunction // Modify EvaluateIndex at the boundaries using mirror boundary conditions this->ApplyMirrorBoundaryConditions(EvaluateIndex, m_SplineOrder); - + const InputImageType * inputImage = this->GetInputImage(); const typename InputImageType::SpacingType & spacing = inputImage->GetSpacing(); @@ -214,39 +220,32 @@ VectorBSplineInterpolateImageFunction CovariantVectorType derivativeValue; double tempValue; IndexType coefficientIndex; - for (unsigned int n = 0; n < ImageDimension; n++) - { + for (unsigned int n = 0; n < ImageDimension; n++) { derivativeValue[n] = 0.0; - for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++) - { - tempValue = 1.0 ; - for (unsigned int n1 = 0; n1 < ImageDimension; n1++) - { + for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++) { + tempValue = 1.0 ; + for (unsigned int n1 = 0; n1 < ImageDimension; n1++) { //coefficientIndex[n1] = EvaluateIndex[n1][sp]; coefficientIndex[n1] = EvaluateIndex[n1][m_PointsToIndex[p][n1]]; - if (n1 == n) - { + if (n1 == n) { //w *= weights[n][ m_PointsToIndex[p][n] ]; tempValue *= weightsDerivative[n1][ m_PointsToIndex[p][n1] ]; - } - else - { - tempValue *= weights[n1][ m_PointsToIndex[p][n1] ]; - } + } else { + tempValue *= weights[n1][ m_PointsToIndex[p][n1] ]; } - derivativeValue[n] += m_Coefficients->GetPixel(coefficientIndex) * tempValue ; } - derivativeValue[n] /= spacing[n]; // take spacing into account + derivativeValue[n] += m_Coefficients->GetPixel(coefficientIndex) * tempValue ; } + derivativeValue[n] /= spacing[n]; // take spacing into account + } #ifdef ITK_USE_ORIENTED_IMAGE_DIRECTION - if( this->m_UseImageDirection ) - { + if( this->m_UseImageDirection ) { CovariantVectorType orientedDerivative; inputImage->TransformLocalVectorToPhysicalVector( derivativeValue, orientedDerivative ); return orientedDerivative; - } + } #endif return(derivativeValue); @@ -254,107 +253,100 @@ VectorBSplineInterpolateImageFunction template -void +void VectorBSplineInterpolateImageFunction -::SetInterpolationWeights( const ContinuousIndexType & x, const vnl_matrix & EvaluateIndex, +::SetInterpolationWeights( const ContinuousIndexType & x, const vnl_matrix & EvaluateIndex, vnl_matrix & weights, unsigned int splineOrder ) const { // For speed improvements we could make each case a separate function and use // function pointers to reference the correct weight order. // Left as is for now for readability. double w, w2, w4, t, t0, t1; - - switch (splineOrder) - { - case 3: - for (unsigned int n = 0; n < ImageDimension; n++) - { - w = x[n] - (double) EvaluateIndex[n][1]; - weights[n][3] = (1.0 / 6.0) * w * w * w; - weights[n][0] = (1.0 / 6.0) + 0.5 * w * (w - 1.0) - weights[n][3]; - weights[n][2] = w + weights[n][0] - 2.0 * weights[n][3]; - weights[n][1] = 1.0 - weights[n][0] - weights[n][2] - weights[n][3]; - } - break; - case 0: - for (unsigned int n = 0; n < ImageDimension; n++) - { - weights[n][0] = 1; // implements nearest neighbor - } - break; - case 1: - for (unsigned int n = 0; n < ImageDimension; n++) - { - w = x[n] - (double) EvaluateIndex[n][0]; - weights[n][1] = w; - weights[n][0] = 1.0 - w; - } - break; - case 2: - for (unsigned int n = 0; n < ImageDimension; n++) - { - /* x */ - w = x[n] - (double)EvaluateIndex[n][1]; - weights[n][1] = 0.75 - w * w; - weights[n][2] = 0.5 * (w - weights[n][1] + 1.0); - weights[n][0] = 1.0 - weights[n][1] - weights[n][2]; - } - break; - case 4: - for (unsigned int n = 0; n < ImageDimension; n++) - { - /* x */ - w = x[n] - (double)EvaluateIndex[n][2]; - w2 = w * w; - t = (1.0 / 6.0) * w2; - weights[n][0] = 0.5 - w; - weights[n][0] *= weights[n][0]; - weights[n][0] *= (1.0 / 24.0) * weights[n][0]; - t0 = w * (t - 11.0 / 24.0); - t1 = 19.0 / 96.0 + w2 * (0.25 - t); - weights[n][1] = t1 + t0; - weights[n][3] = t1 - t0; - weights[n][4] = weights[n][0] + t0 + 0.5 * w; - weights[n][2] = 1.0 - weights[n][0] - weights[n][1] - weights[n][3] - weights[n][4]; - } - break; - case 5: - for (unsigned int n = 0; n < ImageDimension; n++) - { - /* x */ - w = x[n] - (double)EvaluateIndex[n][2]; - w2 = w * w; - weights[n][5] = (1.0 / 120.0) * w * w2 * w2; - w2 -= w; - w4 = w2 * w2; - w -= 0.5; - t = w2 * (w2 - 3.0); - weights[n][0] = (1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - weights[n][5]; - t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0); - t1 = (-1.0 / 12.0) * w * (t + 4.0); - weights[n][2] = t0 + t1; - weights[n][3] = t0 - t1; - t0 = (1.0 / 16.0) * (9.0 / 5.0 - t); - t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0); - weights[n][1] = t0 + t1; - weights[n][4] = t0 - t1; - } - 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 (splineOrder) { + case 3: + for (unsigned int n = 0; n < ImageDimension; n++) { + w = x[n] - (double) EvaluateIndex[n][1]; + weights[n][3] = (1.0 / 6.0) * w * w * w; + weights[n][0] = (1.0 / 6.0) + 0.5 * w * (w - 1.0) - weights[n][3]; + weights[n][2] = w + weights[n][0] - 2.0 * weights[n][3]; + weights[n][1] = 1.0 - weights[n][0] - weights[n][2] - weights[n][3]; + } + break; + case 0: + for (unsigned int n = 0; n < ImageDimension; n++) { + weights[n][0] = 1; // implements nearest neighbor + } + break; + case 1: + for (unsigned int n = 0; n < ImageDimension; n++) { + w = x[n] - (double) EvaluateIndex[n][0]; + weights[n][1] = w; + weights[n][0] = 1.0 - w; } - + break; + case 2: + for (unsigned int n = 0; n < ImageDimension; n++) { + /* x */ + w = x[n] - (double)EvaluateIndex[n][1]; + weights[n][1] = 0.75 - w * w; + weights[n][2] = 0.5 * (w - weights[n][1] + 1.0); + weights[n][0] = 1.0 - weights[n][1] - weights[n][2]; + } + break; + case 4: + for (unsigned int n = 0; n < ImageDimension; n++) { + /* x */ + w = x[n] - (double)EvaluateIndex[n][2]; + w2 = w * w; + t = (1.0 / 6.0) * w2; + weights[n][0] = 0.5 - w; + weights[n][0] *= weights[n][0]; + weights[n][0] *= (1.0 / 24.0) * weights[n][0]; + t0 = w * (t - 11.0 / 24.0); + t1 = 19.0 / 96.0 + w2 * (0.25 - t); + weights[n][1] = t1 + t0; + weights[n][3] = t1 - t0; + weights[n][4] = weights[n][0] + t0 + 0.5 * w; + weights[n][2] = 1.0 - weights[n][0] - weights[n][1] - weights[n][3] - weights[n][4]; + } + break; + case 5: + for (unsigned int n = 0; n < ImageDimension; n++) { + /* x */ + w = x[n] - (double)EvaluateIndex[n][2]; + w2 = w * w; + weights[n][5] = (1.0 / 120.0) * w * w2 * w2; + w2 -= w; + w4 = w2 * w2; + w -= 0.5; + t = w2 * (w2 - 3.0); + weights[n][0] = (1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - weights[n][5]; + t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0); + t1 = (-1.0 / 12.0) * w * (t + 4.0); + weights[n][2] = t0 + t1; + weights[n][3] = t0 - t1; + t0 = (1.0 / 16.0) * (9.0 / 5.0 - t); + t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0); + weights[n][1] = t0 + t1; + weights[n][4] = t0 - t1; + } + 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; + } + } template -void +void VectorBSplineInterpolateImageFunction -::SetDerivativeWeights( const ContinuousIndexType & x, const vnl_matrix & EvaluateIndex, +::SetDerivativeWeights( const ContinuousIndexType & x, const vnl_matrix & EvaluateIndex, vnl_matrix & weights, unsigned int splineOrder ) const { // For speed improvements we could make each case a separate function and use @@ -364,103 +356,96 @@ VectorBSplineInterpolateImageFunction // Left as is for now for readability. double w, w1, w2, w3, w4, w5, t, t0, t1, t2; int derivativeSplineOrder = (int) splineOrder -1; - - switch (derivativeSplineOrder) - { - - // Calculates B(splineOrder) ( (x + 1/2) - xi) - B(splineOrder -1) ( (x - 1/2) - xi) - case -1: - // Why would we want to do this? - for (unsigned int n = 0; n < ImageDimension; n++) - { - weights[n][0] = 0.0; - } - break; - case 0: - for (unsigned int n = 0; n < ImageDimension; n++) - { - weights[n][0] = -1.0; - weights[n][1] = 1.0; - } - break; - case 1: - for (unsigned int n = 0; n < ImageDimension; n++) - { - w = x[n] + 0.5 - (double)EvaluateIndex[n][1]; - // w2 = w; - w1 = 1.0 - w; - weights[n][0] = 0.0 - w1; - weights[n][1] = w1 - w; - weights[n][2] = w; - } - break; - case 2: - - for (unsigned int n = 0; n < ImageDimension; n++) - { - w = x[n] + .5 - (double)EvaluateIndex[n][2]; - w2 = 0.75 - w * w; - w3 = 0.5 * (w - w2 + 1.0); - w1 = 1.0 - w2 - w3; - - weights[n][0] = 0.0 - w1; - weights[n][1] = w1 - w2; - weights[n][2] = w2 - w3; - weights[n][3] = w3; - } - break; - case 3: - - for (unsigned int n = 0; n < ImageDimension; n++) - { - w = x[n] + 0.5 - (double)EvaluateIndex[n][2]; - w4 = (1.0 / 6.0) * w * w * w; - w1 = (1.0 / 6.0) + 0.5 * w * (w - 1.0) - w4; - w3 = w + w1 - 2.0 * w4; - w2 = 1.0 - w1 - w3 - w4; - - weights[n][0] = 0.0 - w1; - weights[n][1] = w1 - w2; - weights[n][2] = w2 - w3; - weights[n][3] = w3 - w4; - weights[n][4] = w4; - } - break; - case 4: - for (unsigned int n = 0; n < ImageDimension; n++) - { - w = x[n] + .5 - (double)EvaluateIndex[n][3]; - t2 = w * w; - t = (1.0 / 6.0) * t2; - w1 = 0.5 - w; - w1 *= w1; - w1 *= (1.0 / 24.0) * w1; - t0 = w * (t - 11.0 / 24.0); - t1 = 19.0 / 96.0 + t2 * (0.25 - t); - w2 = t1 + t0; - w4 = t1 - t0; - w5 = w1 + t0 + 0.5 * w; - w3 = 1.0 - w1 - w2 - w4 - w5; - - weights[n][0] = 0.0 - w1; - weights[n][1] = w1 - w2; - weights[n][2] = w2 - w3; - weights[n][3] = w3 - w4; - weights[n][4] = w4 - w5; - weights[n][5] = w5; - } - break; - - default: - // SplineOrder not implemented yet. - itk::ExceptionObject err(__FILE__, __LINE__); - err.SetLocation( ITK_LOCATION ); - err.SetDescription( "SplineOrder (for derivatives) must be between 1 and 5. Requested spline order has not been implemented yet." ); - throw err; - break; + switch (derivativeSplineOrder) { + + // Calculates B(splineOrder) ( (x + 1/2) - xi) - B(splineOrder -1) ( (x - 1/2) - xi) + case -1: + // Why would we want to do this? + for (unsigned int n = 0; n < ImageDimension; n++) { + weights[n][0] = 0.0; + } + break; + case 0: + for (unsigned int n = 0; n < ImageDimension; n++) { + weights[n][0] = -1.0; + weights[n][1] = 1.0; + } + break; + case 1: + for (unsigned int n = 0; n < ImageDimension; n++) { + w = x[n] + 0.5 - (double)EvaluateIndex[n][1]; + // w2 = w; + w1 = 1.0 - w; + + weights[n][0] = 0.0 - w1; + weights[n][1] = w1 - w; + weights[n][2] = w; + } + break; + case 2: + + for (unsigned int n = 0; n < ImageDimension; n++) { + w = x[n] + .5 - (double)EvaluateIndex[n][2]; + w2 = 0.75 - w * w; + w3 = 0.5 * (w - w2 + 1.0); + w1 = 1.0 - w2 - w3; + + weights[n][0] = 0.0 - w1; + weights[n][1] = w1 - w2; + weights[n][2] = w2 - w3; + weights[n][3] = w3; } - + break; + case 3: + + for (unsigned int n = 0; n < ImageDimension; n++) { + w = x[n] + 0.5 - (double)EvaluateIndex[n][2]; + w4 = (1.0 / 6.0) * w * w * w; + w1 = (1.0 / 6.0) + 0.5 * w * (w - 1.0) - w4; + w3 = w + w1 - 2.0 * w4; + w2 = 1.0 - w1 - w3 - w4; + + weights[n][0] = 0.0 - w1; + weights[n][1] = w1 - w2; + weights[n][2] = w2 - w3; + weights[n][3] = w3 - w4; + weights[n][4] = w4; + } + break; + case 4: + for (unsigned int n = 0; n < ImageDimension; n++) { + w = x[n] + .5 - (double)EvaluateIndex[n][3]; + t2 = w * w; + t = (1.0 / 6.0) * t2; + w1 = 0.5 - w; + w1 *= w1; + w1 *= (1.0 / 24.0) * w1; + t0 = w * (t - 11.0 / 24.0); + t1 = 19.0 / 96.0 + t2 * (0.25 - t); + w2 = t1 + t0; + w4 = t1 - t0; + w5 = w1 + t0 + 0.5 * w; + w3 = 1.0 - w1 - w2 - w4 - w5; + + weights[n][0] = 0.0 - w1; + weights[n][1] = w1 - w2; + weights[n][2] = w2 - w3; + weights[n][3] = w3 - w4; + weights[n][4] = w4 - w5; + weights[n][5] = w5; + } + break; + + default: + // SplineOrder not implemented yet. + itk::ExceptionObject err(__FILE__, __LINE__); + err.SetLocation( ITK_LOCATION ); + err.SetDescription( "SplineOrder (for derivatives) must be between 1 and 5. Requested spline order has not been implemented yet." ); + throw err; + break; + } + } @@ -473,87 +458,71 @@ VectorBSplineInterpolateImageFunction // m_PointsToIndex is used to convert a sequential location to an N-dimension // index vector. This is precomputed to save time during the interpolation routine. m_PointsToIndex.resize(m_MaxNumberInterpolationPoints); - for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++) - { + for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++) { int pp = p; unsigned long indexFactor[ImageDimension]; indexFactor[0] = 1; - for (int j=1; j< static_cast(ImageDimension); j++) - { + for (int j=1; j< static_cast(ImageDimension); j++) { indexFactor[j] = indexFactor[j-1] * ( m_SplineOrder + 1 ); - } - for (int j = (static_cast(ImageDimension) - 1); j >= 0; j--) - { + } + for (int j = (static_cast(ImageDimension) - 1); j >= 0; j--) { m_PointsToIndex[p][j] = pp / indexFactor[j]; pp = pp % indexFactor[j]; - } } + } } template void VectorBSplineInterpolateImageFunction -::DetermineRegionOfSupport( vnl_matrix & evaluateIndex, - const ContinuousIndexType & x, +::DetermineRegionOfSupport( vnl_matrix & evaluateIndex, + const ContinuousIndexType & x, unsigned int splineOrder ) const -{ +{ long indx; -// compute the interpolation indexes - for (unsigned int n = 0; n< ImageDimension; n++) - { - if (splineOrder & 1) // Use this index calculation for odd splineOrder - { +// compute the interpolation indexes + for (unsigned int n = 0; n< ImageDimension; n++) { + if (splineOrder & 1) { // Use this index calculation for odd splineOrder indx = (long)vcl_floor((float)x[n]) - splineOrder / 2; - for (unsigned int k = 0; k <= splineOrder; k++) - { + for (unsigned int k = 0; k <= splineOrder; k++) { evaluateIndex[n][k] = indx++; - } } - else // Use this index calculation for even splineOrder - { + } else { // Use this index calculation for even splineOrder indx = (long)vcl_floor((float)(x[n] + 0.5)) - splineOrder / 2; - for (unsigned int k = 0; k <= splineOrder; k++) - { + for (unsigned int k = 0; k <= splineOrder; k++) { evaluateIndex[n][k] = indx++; - } } } + } } template void VectorBSplineInterpolateImageFunction -::ApplyMirrorBoundaryConditions(vnl_matrix & evaluateIndex, +::ApplyMirrorBoundaryConditions(vnl_matrix & evaluateIndex, unsigned int splineOrder) const { - for (unsigned int n = 0; n < ImageDimension; n++) - { + for (unsigned int n = 0; n < ImageDimension; n++) { long dataLength2 = 2 * m_DataLength[n] - 2; - // apply the mirror boundary conditions + // apply the mirror boundary conditions // TODO: We could implement other boundary options beside mirror - if (m_DataLength[n] == 1) - { - for (unsigned int k = 0; k <= splineOrder; k++) - { + if (m_DataLength[n] == 1) { + for (unsigned int k = 0; k <= splineOrder; k++) { evaluateIndex[n][k] = 0; - } } - else - { - for (unsigned int k = 0; k <= splineOrder; k++) - { + } else { + for (unsigned int k = 0; k <= splineOrder; k++) { // btw - Think about this couldn't this be replaced with a more elagent modulus method? evaluateIndex[n][k] = (evaluateIndex[n][k] < 0L) ? (-evaluateIndex[n][k] - dataLength2 * ((-evaluateIndex[n][k]) / dataLength2)) - : (evaluateIndex[n][k] - dataLength2 * (evaluateIndex[n][k] / dataLength2)); - if ((long) m_DataLength[n] <= evaluateIndex[n][k]) - { + : (evaluateIndex[n][k] - dataLength2 * (evaluateIndex[n][k] / dataLength2)); + if ((long) m_DataLength[n] <= evaluateIndex[n][k]) { evaluateIndex[n][k] = dataLength2 - evaluateIndex[n][k]; - } } } } + } } } // namespace itk