[Rtk-users] Understanding projection matrices management in rtkCudaBackProjectionImageFilter

[Simon Rit]

A small correction: it's actually the amount of data pre-processed (read
from disk + all the processing before backprojection) but that should
change nothing to what is backprojected.
I remembered something that can help you: you may want to run
UpdateLargestPossibleRegion() to the projections (input 1 of the
BackProjectImageFilter) before backprojection to avoid this and make things
clearer instead of using a large volume as I was suggesting before.
Simon

On Tue, Nov 10, 2015 at 2:34 PM, Notargiacomo Thibault <gnthibault at gmail.com
> wrote:

> Ok, using geometric information to reduce the amount of data to be
> backprojected is pretty clever ;)
> This explains the "random" index I was asking about in my previous
> question, I now know what it is used for and continue my debugging process.
>
> Kind regards
>
> Thibault Notargiacomo
>
> 2015-11-10 14:21 GMT+01:00 Simon Rit <simon.rit at creatis.insa-lyon.fr>:
>
>> Hi Thibault,
>> You started digging in the code so that gets complicated. First thing
>> first, yes, the link you provide is still valid.
>> The offset in X is not random but computed in
>> BackProjectionImageFilter::GenerateInputRequestedRegion
>> <https://github.com/SimonRit/RTK/blob/master/code/rtkBackProjectionImageFilter.txx>.
>> The idea is that we request only the part of the projections we need to do
>> the backprojection. This optimizes the streaming (we only read from disk
>> what we need). If you use a huge volume (the test I'd suggest to validate
>> that I'm correct and there is nothing wrong with your code), then you'll
>> need the full projection and you will always get 0 in X as in Y. If you use
>> a small volume and do a computation projection per projection, you will get
>> a varying index depending on your input volume and the geometry for that
>> projection.
>> I think you understood the rest pretty well. Does that make sense or am I
>> missing something else?
>> Simon
>>
>> On Tue, Nov 10, 2015 at 1:39 PM, Notargiacomo Thibault <
>> gnthibault at gmail.com> wrote:
>>
>>> Hi all,
>>>
>>> I am currently trying to use the "AddProjectionInRadians" method from
>>> the ThreeDCircularProjectionGeometry class in RTK to import my geometry.
>>>
>>> First things first, I would like to know if the information on this pdf
>>> are still valid: http://www.openrtk.org/Doxygen/geometry.pdf
>>>
>>> Because I used them to deduce the source offset and Projection offset
>>> from my own matrix decomposition.
>>>
>>> Notice that I don't want to use geometric informations coming from the
>>> projection images themselves, so their geometric origin is assumed to be
>>> (0,0,0) in my case.
>>>
>>> In order to check if my geometry import method is correct, I print out
>>> one of the intermediate form of the projection matrix used in the
>>> backprojection process, in order to check if it is the same as in my own
>>> implementation.
>>>
>>> To do so, I print out one of the intermediate matrix that is generated
>>> in rtkCudaBackProjectionImageFilter while launching a SART reconstruction
>>> for instance.
>>>
>>> I was surpised by the fact that, the matrix allowing the transformation
>>> from 2D geometric world to 2D image index world, namely matrixIdxProj in
>>> rtkCudaBackProjectionImageFilte.cxx
>>> <http://www.openrtk.org/Doxygen/rtkCudaBackProjectionImageFilter_8cxx_source.html>
>>> is not always the same along the reconstruction process ?
>>>
>>> Here is the line that generates the matrix I don't understand:
>>>
>>> for(unsigned int i=0; i<2; i++)
>>>  //SR: 0.5 for 2D texture
>>>   matrixIdxProj[i][2] =
>>> -1*(this->GetInput(1)->GetBufferedRegion().GetIndex()[i])+0.5;
>>>
>>> What I understand here, is that, we are generating a 2D translation
>>> matrix embedded in a 3*3 matrix that allows to take into account two things
>>> related to the 2D projection:
>>>
>>> -If the internal behaviour of the imageContainer, allows buffering of
>>> only a part of the image, beginning at index (i,j), this information will
>>> be embeded in the projection matrix
>>>
>>> -the ITK floating point addressing definition seems to differ from the
>>> one from cuda (and openCV):
>>> When I take a look at this page
>>> <http://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#texture-fetching>
>>> I understand that the floating point indexing (0.5,0.5) in Cuda points to
>>> the center of the pixel (0,0) hence a linear contribution of 100%, while it
>>> seems that we must add a translation of (0.5,0.5) to the ITK geometric
>>> framework in order to get the same behaviour.
>>>
>>> In practice matrixIdxProj is based on a 3*3 identity matrix, where the
>>> inner 2D translation (ie 2 first term of the last column) in my case are:
>>> -in X, a random number between 10 and 150 plus 0.5
>>> -in Y: always 0.5
>>>
>>> I understand the 0.5 part, but why is there always a 0 offset in Y but a
>>> random offset in X remains a mystery to me.
>>>
>>> At first, I believed that each projection was splitted into multiple
>>> part, each having a different X and Y offset.
>>> But reconstructing using only one projection generates only one
>>> backprojection operation, with a non zero offset in X, 123 in my case, for
>>> an image of size 780 in X.
>>>
>>> Do you think there is something wrong whith the values of matrixIdxProjI
>>> am experiencing ?
>>>
>>> Thank you in advance.
>>>
>>> Kind regards
>>>
>>> Thibault Notargiacomo
>>>
>>> _______________________________________________
>>> Rtk-users mailing list
>>> Rtk-users at public.kitware.com
>>> http://public.kitware.com/mailman/listinfo/rtk-users
>>>
>>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/rtk-users/attachments/20151110/d8913dc3/attachment-0010.html>