+/*=========================================================================
+ 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 <itkMath.h>
-/* =========================================================================
-
-@file itkBSplineInterpolateImageFunctionWithLUT.txx
-@author David Sarrut <david.sarrut@creatis.insa-lyon.fr>
-
-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 <class TImageType, class TCoordRep, class TCoefficientType>
BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
BSplineInterpolateImageFunctionWithLUT():Superclass() {
+
// Set default values
- //for(int i=0; i<TImageType::ImageDimension; i++) mOutputSpacing[i] = -1;
SetLUTSamplingFactor(20);
SetSplineOrder(3);
mWeightsAreUpToDate = false;
- // Following need to be pointer beacause values are updated into
- // "const" function
- // mIntrinsecError = new double;
- // mNumberOfError = new long;
- // mIntrinsecErrorMax = new double;
- // mInputIsCoef = false;
}
//====================================================================
//====================================================================
- // //====================================================================
- // template <class TImageType, class TCoordRep, class TCoefficientType>
- // void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
- // SetOutputSpacing(const SpacingType & s) {
- // for(int i=0; i<TImageType::ImageDimension; i++)
- // mOutputSpacing[i] = s[i];
- // // mWeightsAreUpToDate = false;
- // }
- //====================================================================
-
//====================================================================
template <class TImageType, class TCoordRep, class TCoefficientType>
void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
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<TImageType, ITK_TYPENAME itk::NumericTraits<typename TImageType::PixelType>::RealType, TCoefficientType>::SetInputImage(inputData);
- // }
+
+ // Update the weightproperties
if (!inputData) return;
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);
- for(unsigned int l=0; l<d; l++) mSupportIndex[0][l] = 0;
- for(unsigned int k=0; k<mSupportSize; k++) {
- // Get memory offset
- mSupportOffset[k] = Index2Offset<TImageType::ImageDimension>(mSupportIndex[k], mInputMemoryOffset);
- // next coefficient index
- if (k != mSupportSize-1) {
- for(unsigned int l=0; l<d; l++) mSupportIndex[k+1][l] = mSupportIndex[k][l];
- int l=0;
- bool stop = false;
- while (!stop) {
- mSupportIndex[k+1][l]++;
- if (static_cast<unsigned int>(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<d; l++) mSupportIndex[0][l] = 0;
+ for(unsigned int k=0; k<mSupportSize; k++) {
+ // Get memory offset
+ mSupportOffset[k] = Index2Offset<TImageType::ImageDimension>(mSupportIndex[k], mInputMemoryOffset);
+ // next coefficient index
+ if (k != mSupportSize-1) {
+ for(unsigned int l=0; l<d; l++) mSupportIndex[k+1][l] = mSupportIndex[k][l];
+ int l=0;
+ bool stop = false;
+ while (!stop) {
+ mSupportIndex[k+1][l]++;
+ if (static_cast<unsigned int>(mSupportIndex[k+1][l]) == mSupport[l]) {
+ mSupportIndex[k+1][l] = 0;
+ l++;
}
+ else stop = true;
}
-
- // // Check
- // for(unsigned int l=0; l<d; l++) {
- // if (mOutputSpacing[l] == -1) {
- // std::cerr << "Please use SetOutputSpacing before using BLUT interpolator" << std::endl;
- // exit(0);
- // }
- // }
-
- // Compute BSpline weights if not up to date
- //if (!mWeightsAreUpToDate)
- UpdatePrecomputedWeights();
}
- // Initialisation
- // *mIntrinsecErrorMax = *mIntrinsecError = 0.0;
- // *mNumberOfError = 0;
+ }
+
+ // Compute BSpline weights if not up to date
+ if (!mWeightsAreUpToDate)
+ UpdatePrecomputedWeights();
}
//====================================================================
template <class TImageType, class TCoordRep, class TCoefficientType>
void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
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");
}
//====================================================================
for(int l=0; l<TImageType::ImageDimension; l++)
{
- // bool mChange = false;
-
// Compute t1 = distance to floor
- TCoefficientType t1 = x[l]- vcl_floor(x[l]);
+ TCoefficientType t1 = x[l]- std::floor(x[l]);
// Compute index in precomputed weights table
TCoefficientType t2 = mSamplingFactors[l]*t1;
- index[l] = (IndexValueType)lrint(t2);
+ index[l] = (IndexValueType)Math::Round<TCoefficientType>(t2);
// For even order : test if too close to 0.5 (but lower). In this
// case : take the next coefficient
}
}
- // 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]) {
BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
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;
long indx;
for (unsigned int l=0; l<d; l++) {
if (mSplineOrders[l] & 1) { // Use this index calculation for odd splineOrder (like cubic)
- indx = (long)vcl_floor(x[l]) - mSplineOrders[l] / 2 ; //this->m_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<typename ContinuousIndexType::ValueType>(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;
}
}
// Special case for boundary (to be changed ....)
std::vector<int> correctedSupportOffset;
if (boundaryCase) {
- // return -1000;
- //std::vector<TCoefficientType> coef(mSupportSize);
- // DD(EvaluateIndex);
- //std::vector<int> CorrectedSupportOffset;//(mSupportSize);
+
std::vector<IndexType> correctedSupportIndex;//(mSupportSize);
correctedSupportIndex.resize(mSupportSize);
correctedSupportOffset.resize(mSupportSize);
for(unsigned int i=0; i<mSupportSize; i++) {
- // DD(mSupportIndex[i]);
for (unsigned int l=0; l<d; l++) {
long a = mSupportIndex[i][l] + evaluateIndex[l];
long b = this->m_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;
*/
}
}
- // DD(correctedSupportIndex[i]);
correctedSupportOffset[i] = itk::Index2Offset<TImageType::ImageDimension>(correctedSupportIndex[i], mInputMemoryOffset);
}
- // for (unsigned int l=0; l<d; l++) {
- // EvaluateIndex[l] = EvaluateIndex[l] + correctedSupportIndex[0][l];
- // }
- // DD(EvaluateIndex);
psupport = &correctedSupportOffset[0];
}
else {
const TCoefficientType * pcoef = &(this->m_Coefficients->GetPixel(evaluateIndex));
// Main loop over BSpline support
- TCoefficientType interpolated = 0.0;
+ TCoefficientType interpolated = itk::NumericTraits<TCoefficientType>::Zero;
for (unsigned int p=0; p<mSupportSize; p++) {
interpolated += pcoef[*psupport] * (*pweights);
++psupport;
git://git.creatis.insa-lyon.fr / clitk.git / blobdiff