X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=itk%2FitkBSplineInterpolateImageFunctionWithLUT.txx;h=58a173a10cea02189c190a61ec9458b226e9d44c;hb=378ee630dce37a3e15baf3a8027542c2f8cf43de;hp=275bc5396c1c6dd2b02c4a61687ef0236384651f;hpb=931a42358442f4ee4f314613c991c838d4b4e3b7;p=clitk.git diff --git a/itk/itkBSplineInterpolateImageFunctionWithLUT.txx b/itk/itkBSplineInterpolateImageFunctionWithLUT.txx index 275bc53..58a173a 100644 --- a/itk/itkBSplineInterpolateImageFunctionWithLUT.txx +++ b/itk/itkBSplineInterpolateImageFunctionWithLUT.txx @@ -1,40 +1,35 @@ +/*========================================================================= + 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 _ITKINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX #define _ITKINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX +#include -/* ========================================================================= - -@file itkBSplineInterpolateImageFunctionWithLUT.txx -@author David Sarrut - -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 - // void BSplineInterpolateImageFunctionWithLUT:: - // SetOutputSpacing(const SpacingType & s) { - // for(int i=0; i void BSplineInterpolateImageFunctionWithLUT:: @@ -121,28 +106,13 @@ namespace itk void BSplineInterpolateImageFunctionWithLUT:: SetInputImage(const TImageType * inputData) { - //============================== - // if (!mInputIsCoef) - // { + // JV Call superclass (decomposition filter is executeed each time!) + // JV Should call itkBSplineDecompositionFilterWithOBD 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(); - } //==================================================================== @@ -158,47 +128,32 @@ namespace itk 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); - for(unsigned int l=0; l(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; - } + // Compute mSupportOffset according to input size + mSupportOffset.resize(mSupportSize); + mSupportIndex.resize(mSupportSize); + for(unsigned int l=0; l(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; } - - // // 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"); } //==================================================================== @@ -238,14 +186,12 @@ namespace itk for(int l=0; l(t2); // For even order : test if too close to 0.5 (but lower). In this // case : take the next coefficient @@ -259,13 +205,6 @@ namespace itk } } - // Statistics (to be removed) - /* - *mIntrinsecError += fabs(index[l]-t2); - (*mNumberOfError)++; - if (fabs(index[l]-t2)> *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]) { @@ -287,6 +226,9 @@ namespace itk 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; @@ -295,13 +237,13 @@ namespace itk long indx; for (unsigned int l=0; lm_SplineOrder / 2; + indx = (long)std::floor(x[l]) - mSplineOrders[l] / 2 ; //this->m_SplineOrder / 2; evaluateIndex[l] = indx; } else { // Use this index calculation for even splineOrder - if (mSplineOrders[l] == 0) evaluateIndex[l] = (long)rint(x[l]); + if (mSplineOrders[l] == 0) evaluateIndex[l] = Math::Round(x[l]); else { - indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2; + indx = (long)std::floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2; evaluateIndex[l] = indx; } } @@ -326,20 +268,14 @@ namespace itk // 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; @@ -361,13 +297,8 @@ namespace itk */ } } - // DD(correctedSupportIndex[i]); correctedSupportOffset[i] = itk::Index2Offset(correctedSupportIndex[i], mInputMemoryOffset); } - // for (unsigned int l=0; lm_Coefficients->GetPixel(evaluateIndex)); // Main loop over BSpline support - TCoefficientType interpolated = 0.0; + TCoefficientType interpolated = itk::NumericTraits::Zero; for (unsigned int p=0; p