[Rtk-users] SIRT Reconstruction

[Cyril Mory]

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