[Rtk-users] Understanding projection matrices management in rtkCudaBackProjectionImageFilter

[Simon Rit]

Your welcome (I cc the mailing list for info).
Simon

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

> I replaced the original Update by the "UpdateLargestPossibleRegion" to the
> divideFilter in SART (I bypassed displacedDetectorFilter for now) and it
> work as expected, I have now a "clean" matProjIdx.
>
> Thank you for your help !
>
> 2015-11-10 15:03 GMT+01:00 Simon Rit <simon.rit at creatis.insa-lyon.fr>:
>
>> 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/b7eba683/attachment-0010.html>