[Insight-developers] What is LINEAR_INTERPOLATOR_FIXED?
Luis Ibanez
luis.ibanez at kitware.com
Tue Jul 28 12:48:37 EDT 2009
Hi Tom,
Thanks a lot for looking at this problem and for fixing it.
The patch look fine to me,
and the Green submission that you produced is very reassuring. :-)
Please feel free to commit it at your convenience.
Thanks
Luis
-----------------------
Tom Vercauteren wrote:
> Hi,
>
> The attached patch enables optimized linear interpolation when
> ITK_USE_OPTIMIZED_REGISTRATION_METHODS is turned on. It also fixes the
> behavior of the optimized linear interpolation. On my setup (optimized
> registration turned on & centered pixels turned on), an experimental
> build turns out as green as before:
> http://www.cdash.org/CDash/buildSummary.php?buildid=388451
>
> Anyone against me committing this patch?
>
> Cheers,
> Tom
>
> On Tue, Jul 28, 2009 at 00:13, Tom Vercauteren<tom.vercauteren at m4x.org> wrote:
>
>>Ok. Thanks Luis. I'll give it a try.
>>
>>Tom
>>
>>On Mon, Jul 27, 2009 at 22:43, Luis Ibanez<luis.ibanez at kitware.com> wrote:
>>
>>>Hi Tom,
>>>
>>>
>>>The linear interpolator in the Code/Review directory has a bug.
>>>Currently, it doesn't produce the same results as the linear
>>>interpolator in Code/Common.
>>>
>>>
>>>When you enable that CMake flag, you should see a test failing...
>>>
>>>
>>>We should fix its behavior before we rely on this code.
>>>
>>>
>>>
>>> Luis
>>>
>>>
>>>
>>>-----------------------
>>>Tom Vercauteren wrote:
>>>
>>>>Hi Luis,
>>>>
>>>>I just realized that the optimized registration framework does not use
>>>>the optimized linear interpolator. Apparently the optimized linear
>>>>interpolator is not used at all. This is because of the following line
>>>>in itkLinearInterpolateImageFunction.h:
>>>>
>>>>#if defined( ITK_USE_OPTIMIZED_REGISTRATION_METHODS ) && defined(
>>>>LINEAR_INTERPOLATOR_FIXED )
>>>>#include "itkOptLinearInterpolateImageFunction.h"
>>>>#else
>>>>
>>>>Any reason for not having
>>>>
>>>>#ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
>>>>#include "itkOptLinearInterpolateImageFunction.h"
>>>>#else
>>>>
>>>>instead? I saw that this modification was done for ITK 3.6
>>>>
>>>>
>>>> http://www.itk.org/cgi-bin/viewcvs.cgi/Code/Common/itkLinearInterpolateImageFunction.h?root=Insight&r1=1.36&r2=1.33&sortby=date
>>>>
>>>>but the only commit message is
>>>>
>>>> ENH: Temporarily disabling the redirection of this interpolator to the
>>>> OptLinearInterpolateImageFunction class.
>>>>
>>>>Should we proceed to re-enabling this functionality?
>>>>
>>>>Regards,
>>>>Tom
>>>>_______________________________________________
>>>>Powered by www.kitware.com
>>>>
>>>>Visit other Kitware open-source projects at
>>>>http://www.kitware.com/opensource/opensource.html
>>>>
>>>>Please keep messages on-topic and check the ITK FAQ at:
>>>>http://www.itk.org/Wiki/ITK_FAQ
>>>>
>>>>Follow this link to subscribe/unsubscribe:
>>>>http://www.itk.org/mailman/listinfo/insight-developers
>>>>
>>>
>>
>>------------------------------------------------------------------------
>>
>>Index: Code/Common/itkLinearInterpolateImageFunction.h
>>===================================================================
>>RCS file: /cvsroot/Insight/Insight/Code/Common/itkLinearInterpolateImageFunction.h,v
>>retrieving revision 1.36
>>diff -u -r1.36 itkLinearInterpolateImageFunction.h
>>--- Code/Common/itkLinearInterpolateImageFunction.h 20 Mar 2009 10:25:36 -0000 1.36
>>+++ Code/Common/itkLinearInterpolateImageFunction.h 28 Jul 2009 09:06:30 -0000
>>@@ -22,7 +22,7 @@
>> // gets integrated into the main directories.
>> #include "itkConfigure.h"
>>
>>-#if defined( ITK_USE_OPTIMIZED_REGISTRATION_METHODS ) && defined( LINEAR_INTERPOLATOR_FIXED )
>>+#ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
>> #include "itkOptLinearInterpolateImageFunction.h"
>> #else
>>
>>Index: Code/Common/itkLinearInterpolateImageFunction.txx
>>===================================================================
>>RCS file: /cvsroot/Insight/Insight/Code/Common/itkLinearInterpolateImageFunction.txx,v
>>retrieving revision 1.42
>>diff -u -r1.42 itkLinearInterpolateImageFunction.txx
>>--- Code/Common/itkLinearInterpolateImageFunction.txx 7 May 2009 14:03:45 -0000 1.42
>>+++ Code/Common/itkLinearInterpolateImageFunction.txx 28 Jul 2009 09:06:30 -0000
>>@@ -23,7 +23,7 @@
>> #include "itkConfigure.h"
>>
>> // Second, redirect to the optimized version if necessary
>>-#if defined( ITK_USE_OPTIMIZED_REGISTRATION_METHODS ) && defined( LINEAR_INTERPOLATOR_FIXED )
>>+#ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
>> #include "itkOptLinearInterpolateImageFunction.txx"
>> #else
>>
>>@@ -85,27 +85,11 @@
>> */
>> signed long baseIndex[ImageDimension];
>> double distance[ImageDimension];
>>- long tIndex;
>>
>> for( dim = 0; dim < ImageDimension; dim++ )
>> {
>>- // The following "if" block is equivalent to the following line without
>>- // having to call floor.
>>- // baseIndex[dim] = (long) vcl_floor(index[dim] );
>>- if (index[dim] >= 0.0)
>>- {
>>- baseIndex[dim] = (long) index[dim];
>>- }
>>- else
>>- {
>>- tIndex = (long) index[dim];
>>- if (double(tIndex) != index[dim])
>>- {
>>- tIndex--;
>>- }
>>- baseIndex[dim] = tIndex;
>>- }
>>- distance[dim] = index[dim] - double( baseIndex[dim] );
>>+ baseIndex[dim] = Math::Floor( index[dim] );
>>+ distance[dim] = index[dim] - static_cast< double >( baseIndex[dim] );
>> }
>>
>> /**
>>Index: Code/Review/itkOptLinearInterpolateImageFunction.h
>>===================================================================
>>RCS file: /cvsroot/Insight/Insight/Code/Review/itkOptLinearInterpolateImageFunction.h,v
>>retrieving revision 1.6
>>diff -u -r1.6 itkOptLinearInterpolateImageFunction.h
>>--- Code/Review/itkOptLinearInterpolateImageFunction.h 20 Mar 2009 10:25:38 -0000 1.6
>>+++ Code/Review/itkOptLinearInterpolateImageFunction.h 28 Jul 2009 09:06:31 -0000
>>@@ -39,7 +39,7 @@
>> *
>> * \ingroup ImageFunctions ImageInterpolators
>> */
>>-template <class TInputImage, class TCoordRep = float>
>>+template <class TInputImage, class TCoordRep = double>
>> class ITK_EXPORT LinearInterpolateImageFunction :
>> public InterpolateImageFunction<TInputImage,TCoordRep>
>> {
>>@@ -118,21 +118,29 @@
>> {
>> IndexType basei;
>>
>>- double i = index[0];
>>- basei[0] = (long)i;
>>- if( i < 0.0 && double(basei[0]) != i)
>>+ basei[0] = Math::Floor(index[0]);
>>+ if( basei[0] < this->m_StartIndex[0] )
>> {
>>- basei[0]--;
>>+ basei[0] = this->m_StartIndex[0];
>> }
>>
>>- double distance = i - double(basei[0]);
>>+ const double distance = index[0] - static_cast<double>(basei[0]);
>>+
>>+ const RealType val0 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ if(distance <= 0.)
>>+ {
>>+ return( static_cast<OutputType>( val0 ) );
>>+ }
>>
>>- double val0 = this->GetInputImage()->GetPixel( basei );
>>- double val1 = val0;
>> ++basei[0];
>>- val1 = this->GetInputImage()->GetPixel( basei );
>>+ if(basei[0]>this->m_EndIndex[0])
>>+ {
>>+ return( static_cast<OutputType>( val0 ) );
>>+ }
>>+ const RealType val1 = this->GetInputImage()->GetPixel( basei );
>>
>>- return( static_cast<OutputType>( val0 + distance * ( val1 - val0 ) ) );
>>+ return( static_cast<OutputType>( val0 + ( val1 - val0 ) * distance ) );
>> }
>>
>> inline OutputType EvaluateOptimized( const Dispatch<2>&,
>>@@ -140,61 +148,76 @@
>> {
>> IndexType basei;
>>
>>- double i = index[0];
>>- basei[0] = (long)i;
>>- if( i < 0.0 && double(basei[0]) != i)
>>+ basei[0] = Math::Floor(index[0]);
>>+ if( basei[0] < this->m_StartIndex[0] )
>> {
>>- basei[0]--;
>>+ basei[0] = this->m_StartIndex[0];
>> }
>>- double distance0 = i - double(basei[0]);
>>+ const double distance0 = index[0] - static_cast<double>(basei[0]);
>>
>>- i = index[1];
>>- basei[1] = (long)i;
>>- if( i < 0.0 && double(basei[1]) != i)
>>+ basei[1] = Math::Floor(index[1]);
>>+ if( basei[1] < this->m_StartIndex[1] )
>> {
>>- basei[1]--;
>>+ basei[1] = this->m_StartIndex[1];
>> }
>>- double distance1 = i - double(basei[1]);
>>+ const double distance1 = index[1] - static_cast<double>(basei[1]);
>>
>>
>>- double val00 = this->GetInputImage()->GetPixel( basei );
>>- if(distance0+distance1 == 0)
>>+ const RealType val00 = this->GetInputImage()->GetPixel( basei );
>>+ if(distance0 <= 0. && distance1 <= 0.)
>> {
>> return( static_cast<OutputType>( val00 ) );
>> }
>>- else if(distance1 == 0) // if they have the same "y"
>>+ else if(distance1 <= 0.) // if they have the same "y"
>> {
>> ++basei[0]; // then interpolate across "x"
>> if(basei[0]>this->m_EndIndex[0])
>> {
>> return( static_cast<OutputType>( val00 ) );
>> }
>>- double val10 = this->GetInputImage()->GetPixel( basei );
>>- return( static_cast<OutputType>(val00 + distance0 * (val10 - val00)) );
>>+ const RealType val10 = this->GetInputImage()->GetPixel( basei );
>>+ return( static_cast<OutputType>(val00 + (val10 - val00) * distance0) );
>> }
>>- else if(distance0 == 0) // if they have the same "x"
>>+ else if(distance0 <= 0.) // if they have the same "x"
>> {
>> ++basei[1]; // then interpolate across "y"
>> if(basei[1]>this->m_EndIndex[1])
>> {
>> return( static_cast<OutputType>( val00 ) );
>> }
>>- double val01 = this->GetInputImage()->GetPixel( basei );
>>- return( static_cast<OutputType>(val00 + distance1 * (val01 - val00)) );
>>+ const RealType val01 = this->GetInputImage()->GetPixel( basei );
>>+ return( static_cast<OutputType>(val00 + (val01 - val00) * distance1) );
>> }
>>- else
>>+ else // interpolate across "xy"
>> {
>> ++basei[0];
>>- double val10 = this->GetInputImage()->GetPixel( basei );
>>+ if(basei[0]>this->m_EndIndex[0]) // interpolate across "y"
>>+ {
>>+ --basei[0];
>>+ ++basei[1];
>>+ if(basei[1]>this->m_EndIndex[1])
>>+ {
>>+ return( static_cast<OutputType>( val00 ) );
>>+ }
>>+ const RealType val01 = this->GetInputImage()->GetPixel( basei );
>>+ return( static_cast<OutputType>(val00 + (val01 - val00) * distance1) );
>>+ }
>>+ const RealType val10 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ const RealType valx0 = val00 + (val10 - val00) * distance0;
>>+
>> ++basei[1];
>>- double val11 = this->GetInputImage()->GetPixel( basei );
>>+ if(basei[1]>this->m_EndIndex[1]) // interpolate across "x"
>>+ {
>>+ return( static_cast<OutputType>( valx0 ) );
>>+ }
>>+ const RealType val11 = this->GetInputImage()->GetPixel( basei );
>> --basei[0];
>>- double val01 = this->GetInputImage()->GetPixel( basei );
>>-
>>- double val0 = val00 + distance1 * ( val01 - val00 ); // interpolate across "y"
>>- double val1 = val10 + distance1 * ( val11 - val10 ); // interpolate across "y"
>>+ const RealType val01 = this->GetInputImage()->GetPixel( basei );
>>
>>- return( static_cast<OutputType>( val0 + distance0 * (val1-val0) ) ); // interpolate across "X"
>>+ const RealType valx1 = val01 + (val11 - val01) * distance0;
>>+
>>+ return( static_cast<OutputType>( valx0 + (valx1-valx0) * distance1 ) );
>> }
>> }
>>
>>@@ -202,166 +225,267 @@
>> const ContinuousIndexType & index) const
>> {
>> IndexType basei;
>>- double val[8];
>>
>>- unsigned long min0;
>>- unsigned long max0;
>>- unsigned long min1;
>>- unsigned long max1;
>>- unsigned long min2;
>>- unsigned long max2;
>>-
>>- double i = index[0];
>>- basei[0] = (long)i;
>>- if( i < 0.0 && double(basei[0]) != i)
>>+ basei[0] = Math::Floor(index[0]);
>>+ if( basei[0] < this->m_StartIndex[0] )
>> {
>>- basei[0]--;
>>+ basei[0] = this->m_StartIndex[0];
>> }
>>- double distance0 = i - double(basei[0]);
>>+ const double distance0 = index[0] - static_cast<double>(basei[0]);
>>
>>- i = index[1];
>>- basei[1] = (long)i;
>>- if( i < 0.0 && double(basei[1]) != i)
>>+ basei[1] = Math::Floor(index[1]);
>>+ if( basei[1] < this->m_StartIndex[1] )
>> {
>>- basei[1]--;
>>+ basei[1] = this->m_StartIndex[1];
>> }
>>- double distance1 = i - double(basei[1]);
>>+ const double distance1 = index[1] - static_cast<double>(basei[1]);
>>
>>- i = index[2];
>>- basei[2] = (long)i;
>>- if( i < 0.0 && double(basei[2]) != i)
>>+ basei[2] = Math::Floor(index[2]);
>>+ if( basei[2] < this->m_StartIndex[2] )
>> {
>>- basei[2]--;
>>+ basei[2] = this->m_StartIndex[2];
>> }
>>- double distance2 = i - double(basei[2]);
>>+ const double distance2 = index[2] - static_cast<double>(basei[2]);
>>
>>- val[0] = this->GetInputImage()->GetPixel( basei );
>>- if(distance0+distance1+distance2 == 0)
>>- {
>>- return( static_cast<OutputType>( val[0] ) );
>>- }
>>- if(distance0 > 0.0)
>>- {
>>- min0 = basei[0];
>>- max0 = basei[0]+1;
>>- if(max0>this->m_EndIndex[0])
>>- {
>>- max0 = this->m_EndIndex[0];
>>- }
>>- }
>>- if(distance1 > 0.0)
>>- {
>>- min1 = basei[1];
>>- max1 = basei[1]+1;
>>- if(max1>this->m_EndIndex[1])
>>- {
>>- max1 = this->m_EndIndex[1];
>>- }
>>- }
>>- if(distance2 > 0.0)
>>+ if(distance0<=0. && distance1<=0. && distance2<=0.)
>> {
>>- min2 = basei[2];
>>- max2 = basei[2]+1;
>>- if(max2>this->m_EndIndex[2])
>>- {
>>- max2 = this->m_EndIndex[2];
>>- }
>>+ return( static_cast<OutputType>( this->GetInputImage()->GetPixel( basei ) ) );
>> }
>>- if(distance2 == 0)
>>+
>>+ typedef typename IndexType::IndexValueType IndexValueType;
>>+
>>+ const RealType val000 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ if(distance2 <= 0.)
>> {
>>- if(distance1 == 0)
>>+ if(distance1 <= 0.) // interpolate across "x"
>> {
>>- basei[0] = max0;
>>- val[1] = this->GetInputImage()->GetPixel( basei );
>>-
>>- double val0 = val[0] + distance0 * (val[1]-val[0]);
>>+ ++basei[0];
>>+ if(basei[0]>this->m_EndIndex[0])
>>+ {
>>+ return( static_cast<OutputType>( val000 ) );
>>+ }
>>+ const RealType val100 = this->GetInputImage()->GetPixel( basei );
>>
>>- return( static_cast<OutputType>( val0 ) );
>>+ return static_cast<OutputType>( val000 + (val100-val000) * distance0 );
>> }
>>- else if(distance0 == 0)
>>+ else if(distance0 <= 0.) // interpolate across "y"
>> {
>>- basei[1] = max1;
>>- val[2] = this->GetInputImage()->GetPixel( basei );
>>-
>>- double val0 = val[0] + distance1 * (val[2]-val[0]);
>>+ ++basei[1];
>>+ if(basei[1]>this->m_EndIndex[1])
>>+ {
>>+ return( static_cast<OutputType>( val000 ) );
>>+ }
>>+ const RealType val010 = this->GetInputImage()->GetPixel( basei );
>>
>>- return( static_cast<OutputType>( val0 ) );
>>+ return static_cast<OutputType>( val000 + (val010-val000) * distance1 );
>> }
>>- else
>>+ else // interpolate across "xy"
>> {
>>- basei[0] = max0;
>>- val[1] = this->GetInputImage()->GetPixel( basei );
>>- basei[1] = max1;
>>- val[3] = this->GetInputImage()->GetPixel( basei );
>>- basei[0] = min0;
>>- val[2] = this->GetInputImage()->GetPixel( basei );
>>-
>>- double val0 = val[0] + distance0 * (val[1]-val[0]);
>>- double val1 = val[2] + distance0 * (val[3]-val[2]);
>>+ ++basei[0];
>>+ if(basei[0]>this->m_EndIndex[0]) // interpolate across "y"
>>+ {
>>+ --basei[0];
>>+ ++basei[1];
>>+ if(basei[1]>this->m_EndIndex[1])
>>+ {
>>+ return( static_cast<OutputType>( val000 ) );
>>+ }
>>+ const RealType val010 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ return static_cast<OutputType>( val000 + (val010-val000) * distance1 );
>>+ }
>>+ const RealType val100 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ const RealType valx00 = val000 + (val100-val000) * distance0;
>>
>>- return( static_cast<OutputType>( val0 + distance1 * (val1 - val0) ) );
>>+ ++basei[1];
>>+ if(basei[1]>this->m_EndIndex[1]) // interpolate across "x"
>>+ {
>>+ return( static_cast<OutputType>( valx00 ) );
>>+ }
>>+ const RealType val110 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ --basei[0];
>>+ const RealType val010 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ const RealType valx10 = val010 + (val110-val010) * distance0;
>>+
>>+ return static_cast<OutputType>( valx00 + (valx10-valx00) * distance1 );
>> }
>> }
>> else
>> {
>>- basei[2] = max2;
>>- val[4] = this->GetInputImage()->GetPixel( basei );
>>- if(distance1 == 0)
>>+ if(distance1 <= 0.)
>> {
>>- if(distance0 == 0)
>>+ if(distance0 <= 0.) // interpolate across "z"
>>+ {
>>+ ++basei[2];
>>+ if(basei[2]>this->m_EndIndex[2])
>> {
>>- return( static_cast<OutputType>( val[0] + distance2
>>- * (val[4] - val[0]) ) );
>>+ return( static_cast<OutputType>( val000 ) );
>> }
>>- else
>>+ const RealType val001 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ return static_cast<OutputType>( val000 + (val001-val000) * distance2 );
>>+ }
>>+ else // interpolate across "xz"
>> {
>>- basei[0] = max0;
>>- val[5] = this->GetInputImage()->GetPixel( basei );
>>- basei[2] = min2;
>>- val[1] = this->GetInputImage()->GetPixel( basei );
>>-
>>- double val0 = val[0] + distance0 * (val[1]-val[0]);
>>- double val1 = val[4] + distance0 * (val[5]-val[4]);
>>-
>>- return( static_cast<OutputType>( val0 + distance2 * (val1 - val0) ) );
>>+ ++basei[0];
>>+ if(basei[0]>this->m_EndIndex[0]) // interpolate across "z"
>>+ {
>>+ --basei[0];
>>+ ++basei[2];
>>+ if(basei[2]>this->m_EndIndex[2])
>>+ {
>>+ return( static_cast<OutputType>( val000 ) );
>>+ }
>>+ const RealType val001 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ return static_cast<OutputType>( val000 + (val001-val000) * distance2 );
>>+ }
>>+ const RealType val100 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ const RealType valx00 = val000 + (val100-val000) * distance0;
>>+
>>+ ++basei[2];
>>+ if(basei[2]>this->m_EndIndex[2]) // interpolate across "x"
>>+ {
>>+ return( static_cast<OutputType>( valx00 ) );
>>+ }
>>+ const RealType val101 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ --basei[0];
>>+ const RealType val001 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ const RealType valx01 = val001 + (val101-val001) * distance0;
>>+
>>+ return static_cast<OutputType>( valx00 + (valx01-valx00) * distance2 );
>> }
>> }
>>- else if(distance0 == 0)
>>+ else if(distance0 <= 0.) // interpolate across "yz"
>> {
>>- basei[1] = max1;
>>- val[6] = this->GetInputImage()->GetPixel( basei );
>>- basei[2] = min2;
>>- val[2] = this->GetInputImage()->GetPixel( basei );
>>-
>>- double val0 = val[0] + distance1 * (val[2]-val[0]);
>>- double val1 = val[4] + distance1 * (val[6]-val[4]);
>>-
>>- return( static_cast<OutputType>( val0 + distance2 * (val1 - val0) ) );
>>+ ++basei[1];
>>+ if(basei[1]>this->m_EndIndex[1]) // interpolate across "z"
>>+ {
>>+ --basei[1];
>>+ ++basei[2];
>>+ if(basei[2]>this->m_EndIndex[2])
>>+ {
>>+ return( static_cast<OutputType>( val000 ) );
>>+ }
>>+ const RealType val001 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ return static_cast<OutputType>( val000 + (val001-val000) * distance2 );
>>+ }
>>+ const RealType val010 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ const RealType val0x0 = val000 + (val010-val000) * distance1;
>>+
>>+ ++basei[2];
>>+ if(basei[2]>this->m_EndIndex[2]) // interpolate across "y"
>>+ {
>>+ return( static_cast<OutputType>( val0x0 ) );
>>+ }
>>+ const RealType val011 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ --basei[1];
>>+ const RealType val001 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ const RealType val0x1 = val001 + (val011-val001) * distance1;
>>+
>>+ return static_cast<OutputType>( val0x0 + (val0x1-val0x0) * distance2 );
>> }
>>- else
>>+ else // interpolate across "xyz"
>> {
>>- basei[0] = max0;
>>- val[5] = this->GetInputImage()->GetPixel( basei );
>>- basei[1] = max1;
>>- val[7] = this->GetInputImage()->GetPixel( basei );
>>- basei[0] = min0;
>>- val[6] = this->GetInputImage()->GetPixel( basei );
>>- basei[2] = min2;
>>- val[2] = this->GetInputImage()->GetPixel( basei );
>>- basei[0] = max0;
>>- val[3] = this->GetInputImage()->GetPixel( basei );
>>- basei[1] = min1;
>>- val[1] = this->GetInputImage()->GetPixel( basei );
>>-
>>- double val00 = val[0] + distance0 * (val[1]-val[0]);
>>- double val01 = val[2] + distance0 * (val[3]-val[2]);
>>- double val0 = val00 + distance1 * (val01-val00);
>>-
>>- double val10 = val[4] + distance0 * (val[5]-val[4]);
>>- double val11 = val[6] + distance0 * (val[7]-val[6]);
>>- double val1 = val10 + distance1 * (val11-val10);
>>-
>>- return( static_cast<OutputType>( val0 + distance2 * (val1-val0) ) );
>>+ ++basei[0];
>>+ if(basei[0]>this->m_EndIndex[0]) // interpolate across "yz"
>>+ {
>>+ --basei[0];
>>+ ++basei[1];
>>+ if(basei[1]>this->m_EndIndex[1]) // interpolate across "z"
>>+ {
>>+ --basei[1];
>>+ ++basei[2];
>>+ if(basei[2]>this->m_EndIndex[2])
>>+ {
>>+ return( static_cast<OutputType>( val000 ) );
>>+ }
>>+ const RealType val001 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ return static_cast<OutputType>( val000 + (val001-val000) * distance2 );
>>+ }
>>+ const RealType val010 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ const RealType val0x0 = val000 + (val010-val000) * distance1;
>>+
>>+ ++basei[2];
>>+ if(basei[2]>this->m_EndIndex[2]) // interpolate across "y"
>>+ {
>>+ return( static_cast<OutputType>( val0x0 ) );
>>+ }
>>+ const RealType val011 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ --basei[1];
>>+ const RealType val001 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ const RealType val0x1 = val001 + (val011-val001) * distance1;
>>+
>>+ return static_cast<OutputType>( val0x0 + (val0x1-val0x0) * distance2 );
>>+ }
>>+ const RealType val100 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ const RealType valx00 = val000 + (val100-val000) * distance0;
>>+
>>+ ++basei[1];
>>+ if(basei[1]>this->m_EndIndex[1]) // interpolate across "xz"
>>+ {
>>+ --basei[1];
>>+ ++basei[2];
>>+ if(basei[2]>this->m_EndIndex[2]) // interpolate across "x"
>>+ {
>>+ return( static_cast<OutputType>( valx00 ) );
>>+ }
>>+ const RealType val101 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ --basei[0];
>>+ const RealType val001 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ const RealType valx01 = val001 + (val101-val001) * distance0;
>>+
>>+ return static_cast<OutputType>( valx00 + (valx01-valx00) * distance2 );
>>+ }
>>+ const RealType val110 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ --basei[0];
>>+ const RealType val010 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ const RealType valx10 = val010 + (val110-val010) * distance0;
>>+
>>+ const RealType valxx0 = valx00 + (valx10-valx00) * distance1;
>>+
>>+
>>+ ++basei[2];
>>+ if(basei[2]>this->m_EndIndex[2]) // interpolate across "xy"
>>+ {
>>+ return( static_cast<OutputType>( valxx0 ) );
>>+ }
>>+ const RealType val011 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ ++basei[0];
>>+ const RealType val111 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ --basei[1];
>>+ const RealType val101 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ --basei[0];
>>+ const RealType val001 = this->GetInputImage()->GetPixel( basei );
>>+
>>+ const RealType valx01 = val001 + (val101-val001) * distance0;
>>+ const RealType valx11 = val011 + (val111-val011) * distance0;
>>+ const RealType valxx1 = valx01 + (valx11-valx01) * distance1;
>>+
>>+ return( static_cast<OutputType>( valxx0 + (valxx1-valxx0) * distance2 ) );
>> }
>> }
>> }
>>Index: Code/Review/itkOptLinearInterpolateImageFunction.txx
>>===================================================================
>>RCS file: /cvsroot/Insight/Insight/Code/Review/itkOptLinearInterpolateImageFunction.txx,v
>>retrieving revision 1.6
>>diff -u -r1.6 itkOptLinearInterpolateImageFunction.txx
>>--- Code/Review/itkOptLinearInterpolateImageFunction.txx 20 Mar 2009 10:25:37 -0000 1.6
>>+++ Code/Review/itkOptLinearInterpolateImageFunction.txx 28 Jul 2009 09:06:31 -0000
>>@@ -17,7 +17,7 @@
>> #ifndef __itkOptLinearInterpolateImageFunction_txx
>> #define __itkOptLinearInterpolateImageFunction_txx
>>
>>-#include "itkLinearInterpolateImageFunction.h"
>>+#include "itkOptLinearInterpolateImageFunction.h"
>>
>> #include "vnl/vnl_math.h"
>>
>>@@ -77,27 +77,11 @@
>> */
>> signed long baseIndex[ImageDimension];
>> double distance[ImageDimension];
>>- long tIndex;
>>
>> for( dim = 0; dim < ImageDimension; dim++ )
>> {
>>- // The following "if" block is equivalent to the following line without
>>- // having to call floor.
>>- // baseIndex[dim] = (long) vcl_floor(index[dim] );
>>- if (index[dim] >= 0.0)
>>- {
>>- baseIndex[dim] = (long) index[dim];
>>- }
>>- else
>>- {
>>- tIndex = (long) index[dim];
>>- if (double(tIndex) != index[dim])
>>- {
>>- tIndex--;
>>- }
>>- baseIndex[dim] = tIndex;
>>- }
>>- distance[dim] = index[dim] - static_cast< RealType >( baseIndex[dim] );
>>+ baseIndex[dim] = Math::Floor( index[dim] );
>>+ distance[dim] = index[dim] - static_cast< double >( baseIndex[dim] );
>> }
>>
>> /**
>>@@ -124,11 +108,27 @@
>> if ( upper & 1 )
>> {
>> neighIndex[dim] = baseIndex[dim] + 1;
>>+#ifdef ITK_USE_CENTERED_PIXEL_COORDINATES_CONSISTENTLY
>>+ // Take care of the case where the pixel is just
>>+ // in the outer upper boundary of the image grid.
>>+ if( neighIndex[dim] > this->m_EndIndex[dim] )
>>+ {
>>+ neighIndex[dim] = this->m_EndIndex[dim];
>>+ }
>>+#endif
>> overlap *= distance[dim];
>> }
>> else
>> {
>> neighIndex[dim] = baseIndex[dim];
>>+#ifdef ITK_USE_CENTERED_PIXEL_COORDINATES_CONSISTENTLY
>>+ // Take care of the case where the pixel is just
>>+ // in the outer lower boundary of the image grid.
>>+ if( neighIndex[dim] < this->m_StartIndex[dim] )
>>+ {
>>+ neighIndex[dim] = this->m_StartIndex[dim];
>>+ }
>>+#endif
>> overlap *= 1.0 - distance[dim];
>> }
>>
More information about the Insight-developers
mailing list