[Insight-developers] bug 6558

[Tom Vercauteren]

Hi Simon,

I have written a new patch from the patches that you proposed. It now
includes the flag to turn consistent use of centered pixels on or off.

Note that many tests will start failing or crashing. For example
linear interpolation is not yet capable of handling the interpolation
of points that are outside the "strict range". Thus, asking for an
interpolation at continuous index [-0.4  5.2  8.8] will likely lead to
a segfault.

Anyhow, this patch might be a good starting point.

Tom

P.S. The patch can now also be found on the bug tracker:
http://www.itk.org/Bug/view.php?id=6558

On Wed, Nov 19, 2008 at 11:36 AM, Tom Vercauteren
<tom.vercauteren at m4x.org> wrote:
> Hi Simon,
>
> Thanks for this patch which I believe is important. For the record, I
> have just proposed a patch  to vxl to get a fast implementation of
> vnl_math_rnd:
> http://sourceforge.net/mailarchive/forum.php?thread_name=28392e8b0811190226i20ab46a5m2a72a5d9518ce557%40mail.gmail.com&forum_name=vxl-maintainers
>
> I hope both patches get accepted.
>
> Best regards,
> Tom Vercauteren
>
> On Tue, Nov 18, 2008 at 9:36 PM, Simon Warfield
> <simon.warfield at childrens.harvard.edu> wrote:
>>
>> Attached are two patches that modify the transformation of physical points
>> to indexes to use rounding.
>> This ensures the index selected corresponds to the extent of the voxel.
>>
>>>>>> ------------------------
>>>>>> Tom Vercauteren wrote:
>>>>>>
>>>>>>>
>>>>>>> Hi Luis (and others),
>>>>>>>
>>>>>>> I saw the recent refactoring done on the Index <--> Point transforms.
>>>>>>>
>>>>>>> It looks great and easier to maintain (once we get rid of the ifdef
>>>>>>> blocks). I'm just wondering whether this couldn't be a great
>>>>>>> opportunity to fix this bug:
>>>>>>> http://www.itk.org/Bug/view.php?id=6558
>>>>>>>
>>>>>>> What bothers me most is the casts that are similar to:
>>>>>>>  index[i] = static_cast<typename IndexType::IndexValueType>(
>>>>>>> continuous_index[i] );
>>>>>>>
>>>>>>> Using something similar to
>>>>>>>  index[i] = static_cast<typename IndexType::IndexValueType>(
>>>>>>> vnl_math_rnd( continuous_index[i] ) );
>>>>>>> would make more sense to me.
>>>>>>>
>>>>>>> Then could do a performance optimization on vnl_math_rnd. I did such a
>>>>>>> proposition to the vnl team but for some reason, it remained a
>>>>>>> discussion only:
>>>>>>>
>>>>>>>
>>>>>>> http://sourceforge.net/mailarchive/forum.php?thread_name=28392e8b0808050138r31574c90wddd19c0072a7da3c%40mail.gmail.com&forum_name=vxl-maintainers
>>>>>>>
>>>>>>> Best,
>>>>>>> Tom
>>>>>>>
>>>>>>> P.S.: For the record, optimizing vnl_math_rnd & co would also allow
>>>>>>> for fixing these bugs
>>>>>>> http://www.itk.org/Bug/view.php?id=2078
>>>>>>> http://www.itk.org/Bug/view.php?id=5692
>>>>>>> which have long been waiting for a fix
>>>>>>>
>>>>>>>
>>>
>>>
>>>
>>
>>
>> --
>> Simon
>>
>>
>> --- itkImageTransformHelper.h.orig      2008-11-17 12:10:37.000000000 -0500
>> +++ itkImageTransformHelper.h   2008-11-17 13:54:09.000000000 -0500
>> @@ -20,6 +20,7 @@
>>  #include "itkConceptChecking.h"
>>  #include "itkPoint.h"
>>  #include "itkMatrix.h"
>> +#include "vnl/vnl_math.h"
>>
>>  namespace itk
>>  {
>> @@ -173,7 +174,8 @@
>>     const DoublePoint &, DoublePoint &rindex, IndexType &index,
>>     const UniqueTypeBoolTrue& )
>>     {
>> -     index[R] = static_cast<typename IndexType::IndexValueType>(rindex[R]);
>> +     index[R] = static_cast<typename IndexType::IndexValueType>(
>> +                    vnl_math_rnd(rindex[R]) );
>>     }
>>
>>   //
>> @@ -308,7 +310,8 @@
>>     const FloatPoint &, FloatPoint &rindex, IndexType &index,
>>     const UniqueTypeBoolTrue& )
>>     {
>> -     index[R] = static_cast<typename IndexType::IndexValueType>(rindex[R]);
>> +     index[R] = static_cast<typename IndexType::IndexValueType>(
>> +                    vnl_math_rnd(rindex[R]) );
>>     }
>>
>>  };
>>
>> --- itkImageBase.h.orig 2008-11-17 12:15:07.000000000 -0500
>> +++ itkImageBase.h      2008-11-17 12:17:25.000000000 -0500
>> @@ -344,8 +344,8 @@
>>   virtual void SetSpacing (const float spacing[VImageDimension]);
>>
>>
>> -  /** Get the index (discrete) from a physical point.
>> -   * Floating point index results are truncated to integers.
>> +  /** Get the index (discrete) of a voxel from a physical point.
>> +   * Floating point index results are rounded to integers.
>>    * Returns true if the resulting index is within the image, false
>> otherwise
>>    * \sa Transform */
>>  #ifdef ITK_USE_TEMPLATE_META_PROGRAMMING_LOOP_UNROLLING
>> @@ -374,7 +374,7 @@
>>         {
>>         sum += this->m_PhysicalPointToIndex[i][j] * (point[j] -
>> this->m_Origin[j]);
>>         }
>> -      index[i] = static_cast< IndexValueType>( sum );
>> +      index[i] = static_cast< IndexValueType>( vnl_math_round(sum) );
>>       }
>>
>>     // Now, check to see if the index is within allowed bounds
>>
>> _______________________________________________
>> Insight-developers mailing list
>> Insight-developers at itk.org
>> http://www.itk.org/mailman/listinfo/insight-developers
>>
>>
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: itk-bug6558.patch
Type: text/x-diff
Size: 8646 bytes
Desc: not available
URL: <http://www.itk.org/mailman/private/insight-developers/attachments/20081127/7eaa79d7/attachment.patch>