From: srit Date: Wed, 15 Sep 2010 08:16:51 +0000 (+0000) Subject: There were some f***** differences between clitk2 and clitk3 in the BLUT. Tried to... X-Git-Tag: v1.2.0~393 X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=a87fe1bfb8c6b85456952f0b584d44b1da6cb3fa;p=clitk.git There were some f***** differences between clitk2 and clitk3 in the BLUT. Tried to take clitk2 version which is more recent according to cvs. They don't look big but they have an effect. Please check! --- diff --git a/itk/clitkVectorBSplineInterpolateImageFunctionWithLUT.h b/itk/clitkVectorBSplineInterpolateImageFunctionWithLUT.h index eb96f69..e28e768 100644 --- a/itk/clitkVectorBSplineInterpolateImageFunctionWithLUT.h +++ b/itk/clitkVectorBSplineInterpolateImageFunctionWithLUT.h @@ -36,7 +36,6 @@ ========================================================================= */ #include "itkBSplineWeightsCalculator.h" -//#include "clitkTimer.h" #include "clitkVectorBSplineInterpolateImageFunction.h" namespace clitk { @@ -69,7 +68,7 @@ namespace clitk { itkNewMacro(Self); /** Setting LUT sampling (one parameters by dimension or a single - one for all dim); Default is 10 (for each dim) **/ + one for all dim); Default is 20 (for each dim) **/ void SetLUTSamplingFactor(const int& s); void SetLUTSamplingFactors(const SizeType& s); @@ -94,15 +93,6 @@ namespace clitk { various order, dimension, sampling ... **/ static void ComputeBlendingWeights(int dim, int order, int sampling, TCoefficientType * weights); - /** Timer giving computation time for coefficients computation **/ - // const clitk::Timer & GetCoefTimer() const { return mCoefficientTimer; } - - /** Get estimated error **/ - - double GetIntrinsicError() const { return *mIntrinsecError; } - long GetNumberOfError() const { return *mNumberOfError; } - double GetIntrinsicErrorMax() const { return *mIntrinsecErrorMax; } - protected: VectorBSplineInterpolateImageFunctionWithLUT(); ~VectorBSplineInterpolateImageFunctionWithLUT(){;} @@ -117,14 +107,6 @@ namespace clitk { /** Sampling factors for LUT weights **/ SizeType mSamplingFactors; bool mWeightsAreUpToDate; - //SpacingType mOutputSpacing; - - double * mIntrinsecError; - double * mIntrinsecErrorMax; - long * mNumberOfError; - - //JV add iscoeff, and splineorders - // bool mInputIsCoef; SizeType mSplineOrders; // Filter to compute weights @@ -135,15 +117,6 @@ namespace clitk { void UpdateWeightsProperties(); IndexType GetSampleIndexOfPixelPosition(const ContinuousIndexType & x, IndexType & EvaluateIndex) const; - // Timing options - // clitk::Timer mCoefficientTimer; - // clitk::Timer mLUTTimer; - bool mTimerEnabled; - - //JV threadsafety: everything on the stack - //std::vector mCorrectedSupportOffset; - //std::vector mCorrectedSupportIndex; - TCoefficientType * coef; }; // end class clitkVectorBSplineInterpolateImageFunctionWithLUT } // end namespace diff --git a/itk/clitkVectorBSplineInterpolateImageFunctionWithLUT.txx b/itk/clitkVectorBSplineInterpolateImageFunctionWithLUT.txx index 4c6f3eb..337b4c0 100644 --- a/itk/clitkVectorBSplineInterpolateImageFunctionWithLUT.txx +++ b/itk/clitkVectorBSplineInterpolateImageFunctionWithLUT.txx @@ -42,16 +42,9 @@ VectorBSplineInterpolateImageFunctionWithLUT -// void VectorBSplineInterpolateImageFunctionWithLUT:: -// SetOutputSpacing(const SpacingType & s) { -// for(int i=0; i void VectorBSplineInterpolateImageFunctionWithLUT:: @@ -140,27 +123,11 @@ template void VectorBSplineInterpolateImageFunctionWithLUT:: SetInputImage(const TImageType * inputData) { - - //============================== - // if (!mInputIsCoef) - // { + // JV Call superclass (decomposition filter is executeed each time!) + // JV Should call clitkVectorBSplineDecompositionFilterWithOBD to allow different order by dimension Superclass::SetInputImage(inputData); - // } - - //============================== - // //JV Don't call superclass (decomposition filter is executeed each time!) - // else - // { - // this->m_Coefficients = inputData; - // if ( this->m_Coefficients.IsNotNull()) - // { - // this->m_DataLength = this->m_Coefficients->GetBufferedRegion().GetSize(); - // } - - // //Call super-superclass in case more input arrives - // itk::ImageFunction::RealType, TCoefficientType>::SetInputImage(inputData); - // } - + + // Update the weightproperties if (!inputData) return; UpdateWeightsProperties(); @@ -180,8 +147,6 @@ UpdateWeightsProperties() mInputMemoryOffset[l-1]*this->m_Coefficients->GetLargestPossibleRegion().GetSize(l-1); } - //JV Put here? - if (!mWeightsAreUpToDate) { // Compute mSupportOffset according to input size mSupportOffset.resize(mSupportSize); mSupportIndex.resize(mSupportSize); @@ -199,27 +164,16 @@ UpdateWeightsProperties() if (static_cast(mSupportIndex[k+1][l]) == mSupport[l]) { mSupportIndex[k+1][l] = 0; l++; - } else stop = true; + } + else stop = true; } } } - // // Check - // for(unsigned int l=0; l void VectorBSplineInterpolateImageFunctionWithLUT:: UpdatePrecomputedWeights() { - // mLUTTimer.Reset(); - // mLUTTimer.Start(); mWeightsCalculator.SetSplineOrders(mSplineOrders); mWeightsCalculator.SetSamplingFactors(mSamplingFactors); mWeightsCalculator.ComputeTensorProducts(); mWeightsAreUpToDate = true; - //DS - // coef = new TCoefficientType[mSupportSize]; - // mCorrectedSupportIndex.resize(mSupportSize); - // mCorrectedSupportOffset.resize(mSupportSize); - // mLUTTimer.Stop(); - // mLUTTimer.Print("LUT \t"); } //==================================================================== @@ -259,8 +205,6 @@ GetSampleIndexOfPixelPosition(const ContinuousIndexType & x, IndexType & Evaluat IndexType index; for(int l=0; l *mIntrinsecErrorMax) *mIntrinsecErrorMax = fabs(index[l]-t2); - */ - // When to close to 1, take the next coefficient for odd order, but // only change index for odd if (index[l] == (int)mSamplingFactors[l]) { @@ -309,6 +246,8 @@ VectorBSplineInterpolateImageFunctionWithLUT correctedSupportOffset; if (boundaryCase) { - // return -1000; - //std::vector coef(mSupportSize); - // DD(EvaluateIndex); - //std::vector CorrectedSupportOffset;//(mSupportSize); std::vector correctedSupportIndex;//(mSupportSize); correctedSupportIndex.resize(mSupportSize); correctedSupportOffset.resize(mSupportSize); for(unsigned int i=0; im_Coefficients->GetLargestPossibleRegion().GetSize(l); - // DD(a); - // DD(b); long d2 = 2 * b - 2; if (a < 0) { correctedSupportIndex[i][l] = -a - d2*(-a/d2) - evaluateIndex[l];//mSupportIndex[i][l]-a; @@ -383,10 +315,6 @@ EvaluateAtContinuousIndex(const ContinuousIndexType & x) const // DD(correctedSupportIndex[i]); correctedSupportOffset[i] = itk::Index2Offset(correctedSupportIndex[i], mInputMemoryOffset); } - // for (unsigned int l=0; lm_Coefficients->GetPixel(evaluateIndex)); // Main loop over BSpline support - CoefficientImagePixelType interpolated = 0.0; + CoefficientImagePixelType interpolated = itk::NumericTraits::Zero;; for (unsigned int p=0; p ::BSplineDecompositionImageFilterWithOBD() { m_SplineOrder = 0; + for(unsigned int i=0; i ::SetSplineOrders(SizeType SplineOrders) { + if (SplineOrders == m_SplineOrders) + { + return; + } m_SplineOrders=SplineOrders; + this->Modified(); } template diff --git a/itk/itkBSplineInterpolateImageFunctionWithLUT.h b/itk/itkBSplineInterpolateImageFunctionWithLUT.h index 1656bef..2bbc130 100644 --- a/itk/itkBSplineInterpolateImageFunctionWithLUT.h +++ b/itk/itkBSplineInterpolateImageFunctionWithLUT.h @@ -1,22 +1,6 @@ -/*========================================================================= - 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://oncora1.lyon.fnclcc.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 ITKBSPLINEINTERPOLATEIMAGEFUNCTIONWITHLUT_H #define ITKBSPLINEINTERPOLATEIMAGEFUNCTIONWITHLUT_H + /* ========================================================================= @file itkBSplineInterpolateImageFunctionWithLUT.h @@ -36,7 +20,6 @@ ========================================================================= */ #include "itkBSplineWeightsCalculator.h" -//#include "clitkTimer.h" #include namespace itk { @@ -65,7 +48,7 @@ namespace itk { itkNewMacro(Self); /** Setting LUT sampling (one parameters by dimension or a single - one for all dim); Default is 10 (for each dim) **/ + one for all dim); Default is 20 (for each dim) **/ void SetLUTSamplingFactor(const int& s); void SetLUTSamplingFactors(const SizeType& s); @@ -78,9 +61,7 @@ namespace itk { /** Set the input image. This must be set by the user. */ virtual void SetInputImage(const TImageType * inputData); - //void SetOutputSpacing(const SpacingType & s); - //void SetInputImageIsCoefficient(bool inputIsCoef) { mInputIsCoef = inputIsCoef; } - + /** Evaluate the function at a ContinuousIndex position. Overwritten for taking LUT into account */ virtual OutputType EvaluateAtContinuousIndex(const ContinuousIndexType & index ) const; @@ -89,15 +70,6 @@ namespace itk { various order, dimension, sampling ... **/ static void ComputeBlendingWeights(int dim, int order, int sampling, TCoefficientType * weights); - /** Timer giving computation time for coefficients computation **/ - // const clitk::Timer & GetCoefTimer() const { return mCoefficientTimer; } - - /** Get estimated error **/ - - double GetIntrinsicError() const { return *mIntrinsecError; } - long GetNumberOfError() const { return *mNumberOfError; } - double GetIntrinsicErrorMax() const { return *mIntrinsecErrorMax; } - protected: BSplineInterpolateImageFunctionWithLUT(); ~BSplineInterpolateImageFunctionWithLUT(){;} @@ -112,14 +84,6 @@ namespace itk { /** Sampling factors for LUT weights **/ SizeType mSamplingFactors; bool mWeightsAreUpToDate; - //SpacingType mOutputSpacing; - - double * mIntrinsecError; - double * mIntrinsecErrorMax; - long * mNumberOfError; - - //JV add iscoeff, and splineorders - // bool mInputIsCoef; SizeType mSplineOrders; // Filter to compute weights @@ -130,15 +94,6 @@ namespace itk { void UpdateWeightsProperties(); IndexType GetSampleIndexOfPixelPosition(const ContinuousIndexType & x, IndexType & EvaluateIndex) const; - // Timing options - // clitk::Timer mCoefficientTimer; - // clitk::Timer mLUTTimer; - bool mTimerEnabled; - - //JV threadsafety: everything on the stack - //std::vector mCorrectedSupportOffset; - //std::vector mCorrectedSupportIndex; - TCoefficientType * coef; }; // end class itkBSplineInterpolateImageFunctionWithLUT } // end namespace diff --git a/itk/itkBSplineInterpolateImageFunctionWithLUT.txx b/itk/itkBSplineInterpolateImageFunctionWithLUT.txx index b059bc3..c776afa 100644 --- a/itk/itkBSplineInterpolateImageFunctionWithLUT.txx +++ b/itk/itkBSplineInterpolateImageFunctionWithLUT.txx @@ -1,186 +1,132 @@ -/*========================================================================= - 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://oncora1.lyon.fnclcc.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 _ITKINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX #define _ITKINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX -/* ========================================================================= +/* ========================================================================= + @file itkBSplineInterpolateImageFunctionWithLUT.txx @author David Sarrut -Copyright (c) -* CREATIS (Centre de Recherche et d'Applications en Traitement de l'Image). +Copyright (c) +* CREATIS (Centre de Recherche et d'Applications en Traitement de l'Image). All rights reserved. See Doc/License.txt or http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details. * Léon Bérard cancer center, 28 rue Laënnec, 69373 Lyon cedex 08, France * http://www.creatis.insa-lyon.fr/rio - + This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. - + ========================================================================= */ namespace itk { -//==================================================================== -template -BSplineInterpolateImageFunctionWithLUT:: -BSplineInterpolateImageFunctionWithLUT():Superclass() -{ - // Set default values - //for(int i=0; i + BSplineInterpolateImageFunctionWithLUT:: + BSplineInterpolateImageFunctionWithLUT():Superclass() { + + // Set default values + SetLUTSamplingFactor(20); + SetSplineOrder(3); + mWeightsAreUpToDate = false; + } -//==================================================================== -template -void BSplineInterpolateImageFunctionWithLUT:: -SetLUTSamplingFactor(const int& s) -{ - for(int i=0; i + void BSplineInterpolateImageFunctionWithLUT:: + SetLUTSamplingFactor(const int& s) { + for(int i=0; i -void BSplineInterpolateImageFunctionWithLUT:: -SetLUTSamplingFactors(const SizeType& s) -{ - for(int i=0; i -// void BSplineInterpolateImageFunctionWithLUT:: -// SetOutputSpacing(const SpacingType & s) { -// for(int i=0; i -void BSplineInterpolateImageFunctionWithLUT:: -SetSplineOrder(const unsigned int& SplineOrder) -{ - Superclass::SetSplineOrder(SplineOrder); - // Compute support and half size - static const int d = TImageType::ImageDimension; - for(int l=0; l + void BSplineInterpolateImageFunctionWithLUT:: + SetLUTSamplingFactors(const SizeType& s) { + for(int i=0; i -void BSplineInterpolateImageFunctionWithLUT:: -SetSplineOrders(const SizeType& SplineOrders) -{ - mSplineOrders=SplineOrders; - - // Compute support and half size - static const int d = TImageType::ImageDimension; - for(int l=0; l + void BSplineInterpolateImageFunctionWithLUT:: + SetSplineOrder(const unsigned int& SplineOrder) { + Superclass::SetSplineOrder(SplineOrder); + // Compute support and half size + static const int d = TImageType::ImageDimension; + for(int l=0; l + void BSplineInterpolateImageFunctionWithLUT:: + SetSplineOrders(const SizeType& SplineOrders) { + mSplineOrders=SplineOrders; + + // Compute support and half size + static const int d = TImageType::ImageDimension; + for(int l=0; l -void BSplineInterpolateImageFunctionWithLUT:: -SetInputImage(const TImageType * inputData) -{ - - //============================== - // if (!mInputIsCoef) - // { - Superclass::SetInputImage(inputData); - // } - - //============================== - // //JV Don't call superclass (decomposition filter is executeed each time!) - // else - // { - // this->m_Coefficients = inputData; - // if ( this->m_Coefficients.IsNotNull()) - // { - // this->m_DataLength = this->m_Coefficients->GetBufferedRegion().GetSize(); - // } - - // //Call super-superclass in case more input arrives - // itk::ImageFunction::RealType, TCoefficientType>::SetInputImage(inputData); - // } - if (!inputData) return; - UpdateWeightsProperties(); - -} - -//==================================================================== -template -void BSplineInterpolateImageFunctionWithLUT:: -UpdateWeightsProperties() -{ - - // Compute Memory offset inside coefficients images (for looping over coefficients) - static const unsigned int d = TImageType::ImageDimension; - mInputMemoryOffset[0] = 1; - for(unsigned int l=1; lm_Coefficients->GetLargestPossibleRegion().GetSize(l-1); + //==================================================================== + + //==================================================================== + template + void BSplineInterpolateImageFunctionWithLUT:: + SetInputImage(const TImageType * inputData) { + + // JV Call superclass (decomposition filter is executeed each time!) + // JV Should call itkBSplineDecompositionFilterWithOBD to allow different order by dimension + Superclass::SetInputImage(inputData); + + // Update the weightproperties + if (!inputData) return; + UpdateWeightsProperties(); } - - //JV Put here? - if (!mWeightsAreUpToDate) { + + //==================================================================== + template + void BSplineInterpolateImageFunctionWithLUT:: + UpdateWeightsProperties() { + + // Compute Memory offset inside coefficients images (for looping over coefficients) + static const unsigned int d = TImageType::ImageDimension; + mInputMemoryOffset[0] = 1; + for(unsigned int l=1; lm_Coefficients->GetLargestPossibleRegion().GetSize(l-1); + } + // Compute mSupportOffset according to input size mSupportOffset.resize(mSupportSize); mSupportIndex.resize(mSupportSize); @@ -190,223 +136,190 @@ UpdateWeightsProperties() mSupportOffset[k] = Index2Offset(mSupportIndex[k], mInputMemoryOffset); // next coefficient index if (k != mSupportSize-1) { - for(unsigned int l=0; l(mSupportIndex[k+1][l]) == mSupport[l]) { - mSupportIndex[k+1][l] = 0; - l++; - } else stop = true; - } + for(unsigned int l=0; l(mSupportIndex[k+1][l]) == mSupport[l]) { + mSupportIndex[k+1][l] = 0; + l++; + } + else stop = true; + } } } - - // // Check - // for(unsigned int l=0; l -void BSplineInterpolateImageFunctionWithLUT:: -UpdatePrecomputedWeights() -{ - // mLUTTimer.Reset(); - // mLUTTimer.Start(); - mWeightsCalculator.SetSplineOrders(mSplineOrders); - mWeightsCalculator.SetSamplingFactors(mSamplingFactors); - mWeightsCalculator.ComputeTensorProducts(); - mWeightsAreUpToDate = true; - //DS - // coef = new TCoefficientType[mSupportSize]; - // mCorrectedSupportIndex.resize(mSupportSize); - // mCorrectedSupportOffset.resize(mSupportSize); - // mLUTTimer.Stop(); - // mLUTTimer.Print("LUT \t"); -} -//==================================================================== + //==================================================================== + template + void BSplineInterpolateImageFunctionWithLUT:: + UpdatePrecomputedWeights() { -//==================================================================== -template -typename BSplineInterpolateImageFunctionWithLUT::IndexType -BSplineInterpolateImageFunctionWithLUT:: -GetSampleIndexOfPixelPosition(const ContinuousIndexType & x, IndexType & EvaluateIndex) const -{ + mWeightsCalculator.SetSplineOrders(mSplineOrders); + mWeightsCalculator.SetSamplingFactors(mSamplingFactors); + mWeightsCalculator.ComputeTensorProducts(); + mWeightsAreUpToDate = true; + } + //==================================================================== - /* - WARNING : sometimes, a floating number x could not really be - represented in memory. In this case, the difference between x and - floor(x) can be almost 1 (instead of 0). So we take into account - such special case, otherwise it could lead to annoying results. - */ - // static const TCoefficientType tiny = 1.0e-7; - IndexType index; - - for(int l=0; l + typename BSplineInterpolateImageFunctionWithLUT::IndexType + BSplineInterpolateImageFunctionWithLUT:: + GetSampleIndexOfPixelPosition(const ContinuousIndexType & x, IndexType & EvaluateIndex) const { - // Statistics (to be removed) /* - *mIntrinsecError += fabs(index[l]-t2); - (*mNumberOfError)++; - if (fabs(index[l]-t2)> *mIntrinsecErrorMax) *mIntrinsecErrorMax = fabs(index[l]-t2); + WARNING : sometimes, a floating number x could not really be + represented in memory. In this case, the difference between x and + floor(x) can be almost 1 (instead of 0). So we take into account + such special case, otherwise it could lead to annoying results. */ + // static const TCoefficientType tiny = 1.0e-7; + IndexType index; + + for(int l=0; l + typename BSplineInterpolateImageFunctionWithLUT::OutputType + BSplineInterpolateImageFunctionWithLUT:: + EvaluateAtContinuousIndex(const ContinuousIndexType & x) const { -//==================================================================== -template -typename BSplineInterpolateImageFunctionWithLUT::OutputType -BSplineInterpolateImageFunctionWithLUT:: -EvaluateAtContinuousIndex(const ContinuousIndexType & x) const -{ + // JV Compute BSpline weights if not up to date! Problem const: pass image as last + // if (!mWeightsAreUpToDate) UpdatePrecomputedWeights(); + + // For shorter coding + static const unsigned int d = TImageType::ImageDimension; - // For shorter coding - static const unsigned int d = TImageType::ImageDimension; - - // Compute the index of the first interpolation coefficient in the coefficient image - IndexType evaluateIndex; - long indx; - for (unsigned int l=0; lm_SplineOrder / 2; - evaluateIndex[l] = indx; - } else { // Use this index calculation for even splineOrder - if (mSplineOrders[l] == 0) evaluateIndex[l] = (long)rint(x[l]); - else { - indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2; - evaluateIndex[l] = indx; + // Compute the index of the first interpolation coefficient in the coefficient image + IndexType evaluateIndex; + long indx; + for (unsigned int l=0; lm_SplineOrder / 2; + evaluateIndex[l] = indx; + } + else { // Use this index calculation for even splineOrder + if (mSplineOrders[l] == 0) evaluateIndex[l] = (long)rint(x[l]); + else { + indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2; + evaluateIndex[l] = indx; + } } } - } - - // Compute index of precomputed weights and get pointer to first weights - const IndexType weightIndex = GetSampleIndexOfPixelPosition(x, evaluateIndex); - const TCoefficientType * pweights = mWeightsCalculator.GetFirstTensorProduct(weightIndex); - - // Check boundaries - bool boundaryCase = false; - for (unsigned int l=0; l= this->m_Coefficients->GetLargestPossibleRegion().GetSize(l)) { - boundaryCase = true; + + // Compute index of precomputed weights and get pointer to first weights + const IndexType weightIndex = GetSampleIndexOfPixelPosition(x, evaluateIndex); + const TCoefficientType * pweights = mWeightsCalculator.GetFirstTensorProduct(weightIndex); + + // Check boundaries + bool boundaryCase = false; + for (unsigned int l=0; l= this->m_Coefficients->GetLargestPossibleRegion().GetSize(l)) { + boundaryCase = true; + } } - } - // Pointer to support offset - const int * psupport; - - // Special case for boundary (to be changed ....) - std::vector correctedSupportOffset; - if (boundaryCase) { - // return -1000; - //std::vector coef(mSupportSize); - // DD(EvaluateIndex); - //std::vector CorrectedSupportOffset;//(mSupportSize); - std::vector correctedSupportIndex;//(mSupportSize); - correctedSupportIndex.resize(mSupportSize); - correctedSupportOffset.resize(mSupportSize); - for(unsigned int i=0; im_Coefficients->GetLargestPossibleRegion().GetSize(l); - // DD(a); - // DD(b); - long d2 = 2 * b - 2; - if (a < 0) { - correctedSupportIndex[i][l] = -a - d2*(-a/d2) - evaluateIndex[l];//mSupportIndex[i][l]-a; - } else { - if (a>=b) { - correctedSupportIndex[i][l] = d2 - a - evaluateIndex[l]; - } else { - correctedSupportIndex[i][l] = mSupportIndex[i][l]; //a - d2*(a/d2) - EvaluateIndex[l]; - } - /* - if (a>=b) { - correctedSupportIndex[i][l] = d2 - a - EvaluateIndex[l];//mSupportIndex[i][l] - (a-(b-1)); - } - else { - correctedSupportIndex[i][l] = mSupportIndex[i][l]; - } - */ - } + // Pointer to support offset + const int * psupport; + + // Special case for boundary (to be changed ....) + std::vector correctedSupportOffset; + if (boundaryCase) { + + std::vector correctedSupportIndex;//(mSupportSize); + correctedSupportIndex.resize(mSupportSize); + correctedSupportOffset.resize(mSupportSize); + for(unsigned int i=0; im_Coefficients->GetLargestPossibleRegion().GetSize(l); + long d2 = 2 * b - 2; + if (a < 0) { + correctedSupportIndex[i][l] = -a - d2*(-a/d2) - evaluateIndex[l];//mSupportIndex[i][l]-a; + } + else { + if (a>=b) { + correctedSupportIndex[i][l] = d2 - a - evaluateIndex[l]; + } + else { + correctedSupportIndex[i][l] = mSupportIndex[i][l]; //a - d2*(a/d2) - EvaluateIndex[l]; + } + /* + if (a>=b) { + correctedSupportIndex[i][l] = d2 - a - EvaluateIndex[l];//mSupportIndex[i][l] - (a-(b-1)); + } + else { + correctedSupportIndex[i][l] = mSupportIndex[i][l]; + } + */ + } + } + correctedSupportOffset[i] = itk::Index2Offset(correctedSupportIndex[i], mInputMemoryOffset); } - // DD(correctedSupportIndex[i]); - correctedSupportOffset[i] = itk::Index2Offset(correctedSupportIndex[i], mInputMemoryOffset); + psupport = &correctedSupportOffset[0]; + } + else { + psupport = &mSupportOffset[0]; } - // for (unsigned int l=0; lm_Coefficients->GetPixel(evaluateIndex)); + // Get pointer to first coefficient. Even if in some boundary cases, + // EvaluateIndex is out of the coefficient image, + const TCoefficientType * pcoef = &(this->m_Coefficients->GetPixel(evaluateIndex)); - // Main loop over BSpline support - TCoefficientType interpolated = 0.0; - for (unsigned int p=0; p::Zero; + for (unsigned int p=0; p