[Insight-users] Using formulas to imitate result of ResampleImageFilter with itkRayCastInterpolateImageFunction

[Klaus Drechsler]

Dear all,

I have a problem getting the same results when doing the projection 
manaully with formulas compared to a ResampleImageFilter that gets a 
RayCastInterpolateImageFunction and an CenteredEuler3D transform as 
Parameters.

First let me define some abbreviations:

3DV : Original 3D volume
2DP : 2D Projection of the 3DV
3DFP : 3D Feature Points (extracted from 3DV)
2DFP :  2D Projection of the 3D Feature Points

Using a ResampleImageFilter, a CenteredEuler3DTransform (Center = center 
of 3DV) and a RayCastInterpolateImageFunction I get the 2DP.
The CenteredEuler3Dtransform is set as parameter to the 
ResampleImageFilter and the RayCastInterpolateImage Function.
I know  the distance d from focal point to the center of the 3DV and the 
distance e from the center of the 3DV to the projection plane. Thus, the 
focal length is f = d+e.
I also know also the rotation angles and the center of the volume which 
are set as parameters to the CenteredEuler3DTransform.
i know the coordinate of the focal point.

Important: For the sake of simplicity I set the CenteredEuler3DTransform 
to identity (=no rotation).

Now I extract some 3DFP from the 3DV and using the same 
ResampleImageFilter with the same CenteredEuler3DTransform and 
RayCastInterpolateImageFunction I get a 2DFP. 2DP and 2DFP match 
perfectly -> nice!

But i dont like the following:

1. I dont need interpolation
2. I dont want to integrate intensities, because I only have points that 
I want to map to 2D.
3. Too slow
4. I want to remember which 3D point mapped to which 2D point. So i need 
more 'control'.

So, my problem:

Given a list of 3D points p1,p2,...,pn from the 3DV I want to calculate 
the corresponding 2D points q1,q2,...,q3 manually with a formula, so 
that the result matches exactly the 2DFP that I got from the 
ResampleImageFilter.

My try so far:

I assume that the RayCastInterpolateImageFilter produces a 'perspective 
projection' (not with a formulas converting each point, but with rays).
I assume that the projection plane is at z=f;

Using 3D projection (perspective projection) formulas I calculate

p=(x,y,z) -> q=(f/z*x,f/z*y)

But the result is far away from being compareable to the result of the 
ResampleImageFilter.
I played the whole week with different variations of the formula 
(try&error principle) and also with formulas found in the internet, but 
without success so far :-(

Now I am not very far away from giving up, but before that I would like 
to ask this list for help.
I think it could have some to do with coordinate systems and/or using 
physical (mm) values instead of indices, but I am not sure.

What do I oversee? Any ideas and hints are welcome!

If I am missing something in my explanation than ask please.

Thank you very much in advance for your help.

Klaus