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
- BSD See included LICENSE.txt file
- CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
-======================================================================-====*/
+===========================================================================**/
#ifndef _CLITKVECTORBSPLINEINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX
#define _CLITKVECTORBSPLINEINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX
/* =========================================================================
VectorBSplineInterpolateImageFunctionWithLUT():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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
void VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
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<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);
if (static_cast<unsigned int>(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<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();
+ if (!mWeightsAreUpToDate)
+ UpdatePrecomputedWeights();
}
- // Initialisation
- // *mIntrinsecErrorMax = *mIntrinsecError = 0.0;
- // *mNumberOfError = 0;
-}
//====================================================================
//====================================================================
void VectorBSplineInterpolateImageFunctionWithLUT<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");
}
//====================================================================
IndexType index;
for(int l=0; l<TImageType::ImageDimension; l++) {
- // bool mChange = false;
-
// Compute t1 = distance to floor
TCoefficientType t1 = x[l]- vcl_floor(x[l]);
// Compute index in precomputed weights table
TCoefficientType t2 = mSamplingFactors[l]*t1;
- index[l] = (IndexValueType)lrint(t2);
+ index[l] = (IndexValueType)itk::Math::Floor<IndexValueType,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]) {
VectorBSplineInterpolateImageFunctionWithLUT<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;
indx = (long)vcl_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] = itk::Math::Round<long, typename ContinuousIndexType::ValueType>(x[l]);
else {
indx = (long)vcl_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 {
psupport = &mSupportOffset[0];
const CoefficientImagePixelType * pcoef = &(this->m_Coefficients->GetPixel(evaluateIndex));
// Main loop over BSpline support
- CoefficientImagePixelType interpolated = 0.0;
+ CoefficientImagePixelType interpolated = itk::NumericTraits<CoefficientImagePixelType>::Zero;;
for (unsigned int p=0; p<mSupportSize; p++) {
interpolated += pcoef[*psupport] * (*pweights);
++psupport;
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; 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;
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] = itk::Math::Round<long, typename ContinuousIndexType::ValueType>(x[l]);
else {
indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2;
evaluateIndex[l] = indx;
git://git.creatis.insa-lyon.fr / clitk.git / blobdiff