X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=itk%2FclitkVectorBSplineInterpolateImageFunctionWithLUT.txx;h=a237a35871b8866b2238c61eebea1e5386de4752;hb=5f801bf0b07486889123e941d6913d4369dfc86f;hp=4c6f3eb99efd1fb4acb65258159c9c2503ac7588;hpb=1e034c70105f0926939acaaa27ddb46e904ae8bf;p=clitk.git

diff --git a/itk/clitkVectorBSplineInterpolateImageFunctionWithLUT.txx b/itk/clitkVectorBSplineInterpolateImageFunctionWithLUT.txx
index 4c6f3eb..a237a35 100644
--- a/itk/clitkVectorBSplineInterpolateImageFunctionWithLUT.txx
+++ b/itk/clitkVectorBSplineInterpolateImageFunctionWithLUT.txx
@@ -3,7 +3,7 @@
 
   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,7 +14,7 @@
 
   - 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
 /* =========================================================================
@@ -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,14 +205,12 @@ 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]);
 
     // 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
@@ -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]) {
@@ -308,7 +245,8 @@ typename VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoef
 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;
 
@@ -320,7 +258,7 @@ EvaluateAtContinuousIndex(const ContinuousIndexType & x) const
       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;
@@ -347,20 +285,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 +314,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 +325,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,14 +352,13 @@ 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)
       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;