[Rtk-users] SIRT Reconstruction

[Lotte Schyns]

Hi Cyril. I performed a SIRT reconstruction with 100 iterations as you 
suggested. After a whopping 80 hours of calculation time on 64 cores, 
the result was a bit disappointing (in the dropbox: 
https://www.dropbox.com/sh/7a4bkyjg43zjmy7/AAAp4tJrk9HedMEEuIKsGYDSa?dl=0). 
If anything, it looks worse than 3 iterations. Do you think we need even 
more iterations? Another question, you said the conjugate gradient 
method minimizes the same cost function as the SIRT method, but in the 
doxygen documentation I find that the conjugate gradient is similar to 
SART (not SIRT). Could you elaborate a bit on this please? Thanks again!


On 18-07-17 14:28, Cyril Mory wrote:
> The weights are a way to adapt the cost function if you know that some 
> pixels of your projection data are more reliable than others (you then 
> give them a larger weight than the others). This is called Weighted 
> Least Squares (WLS), you can easily find a lot of theory on it. IYou 
> should only do this if you have a reliable way to estimate the 
> variance of the noise on your projection data, and use the inverse 
> variances as weights. There is one weight per pixel, so the weights 
> image should be exactly the same size as the full projection dataset. 
> Note that this input is mandatory for the 
> rtk::ConjugateGradientConeBeamReconstructionFilter, so if you do not 
> want to use it, you still have to generate an image full of ones and 
> pass it as the weights. Clearly that behavior is not optimal, but that 
> is the way it works at the moment.
>
> The support mask enforces a support constraint. It should be a binary 
> volume (ones in the ROI you are interested in reconstructing, zeros 
> outside). If it is provided, over the course of the iterations, the 
> iterates are multiplied by the mask, leaving the voxels inside the ROI 
> unchanged, and setting the ones outside the ROI to zero. 
> Mathematically, it is a special case of preconditioner (a diagonal, 
> binary preconditioner). Not passing a mask yields the same result (but 
> faster) as passing a mask full of ones.
>
> Hope this helps,
> Cyril
>
> PS: If at some point you decide to use the "Weights" input of the 
> conjugate gradient filter, in order to perform WLS reconstructions, do 
> write another email on this list: I have a lot of advice to share on 
> this topic.
>
> On 18/07/2017 14:02, Lotte Schyns wrote:
>> Thanks! I will try both the SIRT with 100 iterations and the conjugate
>> gradient method. What exactly are the weights (input 2) and input
>> support mask in rtkConjugateGradientConeBeamReconstructionFilter?
>>
>> On 18-07-17 11:31, Cyril Mory wrote:
>>> Hi Lotte,
>>>
>>> I had a very quick look at what you are doing. It looks fine to me,
>>> you just need to perform way more iterations in the SIRT case. 100
>>> would be a good start. SIRT is more stable than SART when there are
>>> inconsistencies in the projection data, but converges slowly.
>>>
>>> An alternative to SIRT, which minimizes the same cost function with a
>>> faster algorithm, is the conjugate gradient algorithm. You should
>>> obtain nice results with something like 30-40 iterations (look for
>>> rtk::ConjugateGradientConeBeamReconstructionFilter if you want to give
>>> it a try).
>>>
>>> Best regards,
>>> Cyril
>>>
>>>
>>> On 18/07/2017 11:17, Lotte Schyns wrote:
>>>> Hello,
>>>>
>>>> We are having problems with the SIRT reconstructions. They seem very
>>>> strange and blurry. However, other reconstructions (SART, FDK,
>>>> iterativeFDK) look perfect. I uploaded an example of a SART and a SIRT
>>>> reconstruction (same parameters) to
>>>> https://www.dropbox.com/sh/7a4bkyjg43zjmy7/AAAp4tJrk9HedMEEuIKsGYDSa?dl=0. 
>>>>
>>>>
>>>> I also uploaded a minimalistic version of the code that I used (the
>>>> paths in CMakeLists.txt probably need to be adapted to your 
>>>> system). You
>>>> can alternate between SART and SIRT by changing line 143. The raw data
>>>> is also available on dropbox. I didn't use
>>>> rtkXRadRawToAttenuationImageFilter, because our projections are 
>>>> already
>>>> corrected for the dark field and flood field, so
>>>> (signal-dark)/(flood-dark). I just take the natural logarithm and
>>>> multiply by -1. Do you know what could be the problem with the SIRT
>>>> reconstructions? Are we using wrong parameters? Thanks for your time.
>>>>
>>>> Lotte
>>>> _______________________________________________
>>>> Rtk-users mailing list
>>>> Rtk-users at public.kitware.com
>>>> http://public.kitware.com/mailman/listinfo/rtk-users
>> _______________________________________________
>> Rtk-users mailing list
>> Rtk-users at public.kitware.com
>> http://public.kitware.com/mailman/listinfo/rtk-users
>