From a87fe1bfb8c6b85456952f0b584d44b1da6cb3fa Mon Sep 17 00:00:00 2001
From: srit <srit>
Date: Wed, 15 Sep 2010 08:16:51 +0000
Subject: [PATCH] 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!
---
...orBSplineInterpolateImageFunctionWithLUT.h | 29 +-
...BSplineInterpolateImageFunctionWithLUT.txx | 95 +--
...BSplineDecompositionImageFilterWithOBD.txx | 7 +
...tkBSplineInterpolateImageFunctionWithLUT.h | 51 +-
...BSplineInterpolateImageFunctionWithLUT.txx | 611 ++++++++----------
5 files changed, 284 insertions(+), 509 deletions(-)
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<int> mCorrectedSupportOffset;
- //std::vector<IndexType> 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<TImageType,TCoordRep,TCoefficientTy
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;
}
//====================================================================
@@ -78,16 +71,6 @@ SetLUTSamplingFactors(const SizeType& s)
//====================================================================
-// //====================================================================
-// 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>::
@@ -140,27 +123,11 @@ template <class TImageType, class TCoordRep, class 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();
@@ -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<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;
-}
//====================================================================
//====================================================================
@@ -227,18 +181,10 @@ template <class TImageType, class TCoordRep, class TCoefficientType>
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");
}
//====================================================================
@@ -259,8 +205,6 @@ GetSampleIndexOfPixelPosition(const ContinuousIndexType & x, IndexType & Evaluat
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]);
@@ -280,13 +224,6 @@ GetSampleIndexOfPixelPosition(const ContinuousIndexType & x, IndexType & Evaluat
}
}
- // 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]) {
@@ -309,6 +246,8 @@ VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientTy
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;
@@ -347,20 +286,13 @@ EvaluateAtContinuousIndex(const ContinuousIndexType & x) const
// 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;
@@ -383,10 +315,6 @@ EvaluateAtContinuousIndex(const ContinuousIndexType & x) const
// 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];
@@ -398,7 +326,7 @@ EvaluateAtContinuousIndex(const ContinuousIndexType & x) const
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;
@@ -425,7 +353,6 @@ EvaluateWeightsAtContinuousIndex(const ContinuousIndexType& x, const TCoeffici
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)
diff --git a/itk/itkBSplineDecompositionImageFilterWithOBD.txx b/itk/itkBSplineDecompositionImageFilterWithOBD.txx
index 3812bb5..e268ed7 100644
--- a/itk/itkBSplineDecompositionImageFilterWithOBD.txx
+++ b/itk/itkBSplineDecompositionImageFilterWithOBD.txx
@@ -34,6 +34,8 @@ BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
::BSplineDecompositionImageFilterWithOBD()
{
m_SplineOrder = 0;
+ for(unsigned int i=0; i<ImageDimension; i++)
+ m_SplineOrders[i]=3;
int SplineOrder = 3;
m_Tolerance = 1e-10; // Need some guidance on this one...what is reasonable?
m_IteratorDirection = 0;
@@ -124,7 +126,12 @@ void
BSplineDecompositionImageFilterWithOBD<TInputImage, TOutputImage>
::SetSplineOrders(SizeType SplineOrders)
{
+ if (SplineOrders == m_SplineOrders)
+ {
+ return;
+ }
m_SplineOrders=SplineOrders;
+ this->Modified();
}
template <class TInputImage, class TOutputImage>
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 <itkBSplineInterpolateImageFunction.h>
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<int> mCorrectedSupportOffset;
- //std::vector<IndexType> 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 <david.sarrut@creatis.insa-lyon.fr>
-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 <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>
+ BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
+ BSplineInterpolateImageFunctionWithLUT():Superclass() {
+
+ // Set default values
+ SetLUTSamplingFactor(20);
+ SetSplineOrder(3);
+ mWeightsAreUpToDate = false;
+ }
-//====================================================================
-template <class TImageType, class TCoordRep, class TCoefficientType>
-void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
-SetLUTSamplingFactor(const int& s)
-{
- for(int i=0; i<TImageType::ImageDimension; i++) mSamplingFactors[i] = s;
- mWeightsAreUpToDate = false;
-}
+ //====================================================================
-//====================================================================
+ //====================================================================
+ template <class TImageType, class TCoordRep, class TCoefficientType>
+ void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
+ SetLUTSamplingFactor(const int& s) {
+ for(int i=0; i<TImageType::ImageDimension; i++) mSamplingFactors[i] = s;
+ mWeightsAreUpToDate = false;
+ }
-//====================================================================
-template <class TImageType, class TCoordRep, class TCoefficientType>
-void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
-SetLUTSamplingFactors(const SizeType& s)
-{
- for(int i=0; i<TImageType::ImageDimension; i++) mSamplingFactors[i] = s[i];
- mWeightsAreUpToDate = 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>::
-SetSplineOrder(const unsigned int& SplineOrder)
-{
- Superclass::SetSplineOrder(SplineOrder);
- // Compute support and half size
- static const int d = TImageType::ImageDimension;
- for(int l=0; l<d; l++) {
- mSplineOrders[l]= SplineOrder;
- mSupport[l] = SplineOrder+1;
- if (mSupport[l] % 2 == 0) { // support is even
- mHalfSupport[l] = mSupport[l]/2-1;
- } else mHalfSupport[l] = mSupport[l]/2; // support is odd (like cubic spline)
+ //====================================================================
+ template <class TImageType, class TCoordRep, class TCoefficientType>
+ void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
+ SetLUTSamplingFactors(const SizeType& s) {
+ for(int i=0; i<TImageType::ImageDimension; i++) mSamplingFactors[i] = s[i];
+ mWeightsAreUpToDate = false;
}
- mSupportSize = 1;
- for(int l=0; l<d; l++) {
- mSupportSize *= mSupport[l];
- }
- mWeightsAreUpToDate = false;
-}
-//====================================================================
-//JV
-//====================================================================
-template <class TImageType, class TCoordRep, class TCoefficientType>
-void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
-SetSplineOrders(const SizeType& SplineOrders)
-{
- mSplineOrders=SplineOrders;
-
- // Compute support and half size
- static const int d = TImageType::ImageDimension;
- for(int l=0; l<d; l++) {
- mSupport[l] = mSplineOrders[l]+1;
- if (mSupport[l] % 2 == 0) { // support is even
- mHalfSupport[l] = mSupport[l]/2-1;
- } else mHalfSupport[l] = mSupport[l]/2; // support is odd (like cubic spline)
+ //====================================================================
+
+ //====================================================================
+ template <class TImageType, class TCoordRep, class TCoefficientType>
+ void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
+ SetSplineOrder(const unsigned int& SplineOrder) {
+ Superclass::SetSplineOrder(SplineOrder);
+ // Compute support and half size
+ static const int d = TImageType::ImageDimension;
+ for(int l=0; l<d; l++) {
+ mSplineOrders[l]= SplineOrder;
+ mSupport[l] = SplineOrder+1;
+ if (mSupport[l] % 2 == 0) { // support is even
+ mHalfSupport[l] = mSupport[l]/2-1;
+ }
+ else mHalfSupport[l] = mSupport[l]/2; // support is odd (like cubic spline)
+ }
+ mSupportSize = 1;
+ for(int l=0; l<d; l++) {
+ mSupportSize *= mSupport[l];
+ }
+ mWeightsAreUpToDate = false;
}
- mSupportSize = 1;
- for(int l=0; l<d; l++) {
- mSupportSize *= mSupport[l];
+ //====================================================================
+
+ //JV
+ //====================================================================
+ template <class TImageType, class TCoordRep, class TCoefficientType>
+ void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
+ SetSplineOrders(const SizeType& SplineOrders) {
+ mSplineOrders=SplineOrders;
+
+ // Compute support and half size
+ static const int d = TImageType::ImageDimension;
+ for(int l=0; l<d; l++) {
+ mSupport[l] = mSplineOrders[l]+1;
+ if (mSupport[l] % 2 == 0) { // support is even
+ mHalfSupport[l] = mSupport[l]/2-1;
+ }
+ else mHalfSupport[l] = mSupport[l]/2; // support is odd (like cubic spline)
+ }
+ mSupportSize = 1;
+ for(int l=0; l<d; l++) {
+ mSupportSize *= mSupport[l];
+ }
+ mWeightsAreUpToDate = false;
}
- mWeightsAreUpToDate = false;
-}
-//====================================================================
-
-//====================================================================
-template <class TImageType, class TCoordRep, class TCoefficientType>
-void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
-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<TImageType, ITK_TYPENAME itk::NumericTraits<typename TImageType::PixelType>::RealType, TCoefficientType>::SetInputImage(inputData);
- // }
- if (!inputData) return;
- UpdateWeightsProperties();
-
-}
-
-//====================================================================
-template <class TImageType, class TCoordRep, class TCoefficientType>
-void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
-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; l<d; l++) {
- mInputMemoryOffset[l] =
- mInputMemoryOffset[l-1]*this->m_Coefficients->GetLargestPossibleRegion().GetSize(l-1);
+ //====================================================================
+
+ //====================================================================
+ template <class TImageType, class TCoordRep, class TCoefficientType>
+ void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
+ 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 <class TImageType, class TCoordRep, class TCoefficientType>
+ void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
+ 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; l<d; l++) {
+ mInputMemoryOffset[l] =
+ mInputMemoryOffset[l-1]*this->m_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<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;
- }
+ 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();
+ if (!mWeightsAreUpToDate)
+ UpdatePrecomputedWeights();
}
- // Initialisation
- // *mIntrinsecErrorMax = *mIntrinsecError = 0.0;
- // *mNumberOfError = 0;
-}
-//====================================================================
+ //====================================================================
-//====================================================================
-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");
-}
-//====================================================================
+ //====================================================================
+ template <class TImageType, class TCoordRep, class TCoefficientType>
+ void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
+ UpdatePrecomputedWeights() {
-//====================================================================
-template <class TImageType, class TCoordRep, class TCoefficientType>
-typename BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::IndexType
-BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
-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<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);
-
- // For even order : test if too close to 0.5 (but lower). In this
- // case : take the next coefficient
- if (!(mSplineOrders[l] & 1)) {
- if (t1<0.5) {
- if (mSamplingFactors[l] & 1) {
- if (index[l] == (int) mSamplingFactors[l]/2+1) EvaluateIndex[l] = EvaluateIndex[l]+1;
- }
-
- else if (index[l] == (int) mSamplingFactors[l]/2) EvaluateIndex[l] = EvaluateIndex[l]+1;
- }
- }
+ //====================================================================
+ template <class TImageType, class TCoordRep, class TCoefficientType>
+ typename BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::IndexType
+ BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
+ 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<TImageType::ImageDimension; l++)
+ {
+ // 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);
+
+ // For even order : test if too close to 0.5 (but lower). In this
+ // case : take the next coefficient
+ if (!(mSplineOrders[l] & 1)) {
+ if (t1<0.5) {
+ if (mSamplingFactors[l] & 1) {
+ if (index[l] == (int) mSamplingFactors[l]/2+1) EvaluateIndex[l] = EvaluateIndex[l]+1;
+ }
+
+ else if (index[l] == (int) mSamplingFactors[l]/2) EvaluateIndex[l] = EvaluateIndex[l]+1;
+ }
+ }
+
+ // When to close to 1, take the next coefficient for odd order, but
+ // only change index for odd
+ if (index[l] == (int)mSamplingFactors[l]) {
+ index[l] = 0;
+ if (mSplineOrders[l] & 1) EvaluateIndex[l] = EvaluateIndex[l]+1;
+ }
+ }
- // When to close to 1, take the next coefficient for odd order, but
- // only change index for odd
- if (index[l] == (int)mSamplingFactors[l]) {
- index[l] = 0;
- if (mSplineOrders[l] & 1) EvaluateIndex[l] = EvaluateIndex[l]+1;
- }
+ // The end
+ return index;
}
- // The end
- return index;
-}
+ //====================================================================
-//====================================================================
+ //====================================================================
+ template <class TImageType, class TCoordRep, class TCoefficientType>
+ typename BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::OutputType
+ BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
+ EvaluateAtContinuousIndex(const ContinuousIndexType & x) const {
-//====================================================================
-template <class TImageType, class TCoordRep, class TCoefficientType>
-typename BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::OutputType
-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;
- // 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; 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]);
- 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; 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]);
+ 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<d; l++) {
- if ((evaluateIndex[l] < 0) ||
- (evaluateIndex[l]+mSupport[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<d; l++) {
+ if ((evaluateIndex[l] < 0) ||
+ (evaluateIndex[l]+mSupport[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<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;
- } 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<int> correctedSupportOffset;
+ if (boundaryCase) {
+
+ std::vector<IndexType> correctedSupportIndex;//(mSupportSize);
+ correctedSupportIndex.resize(mSupportSize);
+ correctedSupportOffset.resize(mSupportSize);
+ for(unsigned int i=0; i<mSupportSize; i++) {
+ for (unsigned int l=0; l<d; l++) {
+ long a = mSupportIndex[i][l] + evaluateIndex[l];
+ long b = this->m_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<TImageType::ImageDimension>(correctedSupportIndex[i], mInputMemoryOffset);
}
- // DD(correctedSupportIndex[i]);
- correctedSupportOffset[i] = itk::Index2Offset<TImageType::ImageDimension>(correctedSupportIndex[i], mInputMemoryOffset);
+ psupport = &correctedSupportOffset[0];
+ }
+ else {
+ psupport = &mSupportOffset[0];
}
- // for (unsigned int l=0; l<d; l++) {
- // EvaluateIndex[l] = EvaluateIndex[l] + correctedSupportIndex[0][l];
- // }
- // DD(EvaluateIndex);
- psupport = &correctedSupportOffset[0];
- } else {
- psupport = &mSupportOffset[0];
- }
- // 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));
+ // 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<mSupportSize; p++) {
- interpolated += pcoef[*psupport] * (*pweights);
- ++psupport;
- ++pweights;
- }
+ // Main loop over BSpline support
+ TCoefficientType interpolated = itk::NumericTraits<TCoefficientType>::Zero;
+ for (unsigned int p=0; p<mSupportSize; p++) {
+ interpolated += pcoef[*psupport] * (*pweights);
+ ++psupport;
+ ++pweights;
+ }
- // Return interpolated value
- return(interpolated);
-}
-//====================================================================
+ // Return interpolated value
+ return(interpolated);
+ }
+ //====================================================================
}
#endif //_ITKINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX
--
2.49.0