[Insight-users] Interpolation and centered pixel coordinates

[Dženan Zukić]

Hi Luis,

I found out what the problem is.

The way I call evaluate is OK:

image->TransformPointFromPhysicalToContinouosIndex(p, ind, ok);
if (ok)
  pixel=interpolator->EvaluateAtContinuousIndex(ind);

but my original point p had NaN values in it, which were caused by code
which looked something like this:

double cosFi=vec1*vec2/(vec1.length()*vec2.length());
//cosine of angle between vectors
p+=acos(cosFi)*...;

Due to numerical round-off errors it would rarely occur that cosFi>1
(something like 1.00000002534), and acos would produce NaN, which would then
propagate all the way to poor interpolator which would then complain.
Shortly after starting to use interpolator I also started using different
sort of input files, which made my program do much more iterations of the
above numerically unstable vector code and made manifestations of that
problem more often.

Since the problem was manifesting itself in interpolator, and just before
calling interpolator I was getting its input data from
TransformPhysicalPointToContinuousIndex, and code for this transformation is
very convoluted, I suspected it might be wrong and I had to ask.

I am also quite busy these last few weeks (I barely managed to find time to
fix that bug which I adopted :D), so I soon forgot I asked this question,
and of course didn't write that I found out the problem (which was quite
far, in calculation of normals and curvatures of the mesh).

Thanks for the answer,
Dženan

2010/6/16 Luis Ibanez <luis.ibanez at kitware.com>

>
> Dženan,
>
> We did quite a heavy testing when we fixed
> the problem of centered pixel coordinates.
>
> The Dashboard was green when we released
> ITK 3.18, which give us some confidence on
> the correctness of the code.
>
> --
>
> Please note that it is your responsibility to call
>
>      interpolator->IsInsideBuffer( point )
>
> before you call
>
>      interpolator->Evaluate( point );
>
>
> That is, your code should take the form:
>
>
>    if(    interpolator->IsInsideBuffer( p1 ) )
>       {
>       p2 = interpolator->Evaluate( p1 );
>       }
>
>
> --
>
> If you still find run-time errors, please
> post to the list a minimal code example
> that displays the problem.
>
>
>     Thanks
>
>
>         Luis
>
>
> ----------------------------------------------------------
> 2010/6/11 Dženan Zukić <dzenanz at gmail.com>
>
>> Hi everyone,
>>
>> I am running into some crashes in
>> \InsightToolkit-3.18.0\Code\Common\itkLinearInterpolateImageFunction.txx
>> (line 150), and I wonder could it be due to recent move to consistent use of
>> centered pixels (physical coordinates going from -spacing/2 to
>> (n-1)*spacing+spacing/2, instead of 0 to n*spacing). Did anyone have
>> encounters with similar problems? More specifically, does
>> TransformPhysicalPointToContinuousIndex in
>> \InsightToolkit-3.18.0\Code\Common\itkImageBase.h honor the centered pixel
>> coordinates?
>>
>> Regards,
>> Dženan
>>
>> _____________________________________
>> Powered by www.kitware.com
>>
>> Visit other Kitware open-source projects at
>> http://www.kitware.com/opensource/opensource.html
>>
>> Kitware offers ITK Training Courses, for more information visit:
>> http://www.kitware.com/products/protraining.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-users
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20100616/e1ef4ce7/attachment.htm>