X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=itk%2FitkBSplineInterpolateImageFunctionWithLUT.txx;h=b9354a5e13ce2d13762de9a7f5eefef4f7d410e8;hb=dca456c52728e92190a48a2823c4bb9d7738c8d5;hp=3c9dff72d4e96d550d91075f3e01e755b313ae2d;hpb=0b7c9b1e1215634b02cbd38d4e4ba101d6111ba8;p=clitk.git diff --git a/itk/itkBSplineInterpolateImageFunctionWithLUT.txx b/itk/itkBSplineInterpolateImageFunctionWithLUT.txx index 3c9dff7..b9354a5 100644 --- a/itk/itkBSplineInterpolateImageFunctionWithLUT.txx +++ b/itk/itkBSplineInterpolateImageFunctionWithLUT.txx @@ -1,9 +1,9 @@ /*========================================================================= 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 @@ -14,43 +14,22 @@ - 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). -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. - -========================================================================= */ +#include + 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:: @@ -137,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(); - } //==================================================================== @@ -174,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"); } //==================================================================== @@ -254,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 @@ -275,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]) { @@ -303,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; @@ -315,7 +241,7 @@ namespace itk 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; evaluateIndex[l] = indx; @@ -342,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; @@ -377,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