[Rtk-users] Fwd: High pixel/voxel values in SART
Cyril Mory
cyril.mory at creatis.insa-lyon.fr
Wed Feb 21 11:07:29 EST 2018
Chao,
Thanks for the feedback. It's a very encouraging first result. I'll mark
the trick you used as an issue on github, so as to remember to implement
it as an option and test it on more data.
Indeed, with iterative methods, you have to process full slices. And it
is even worse when you apply spatial regularization, because then
reconstructing slices independently is possible, but less relevant than
reconstructing a full volume.
Cyril
On 21/02/2018 16:52, Chao Wu wrote:
> Hi Cyril,
>
> Thanks for your suggestion.
>
> I have tried increasing the threshold. My reconstruced slices are
> 32x32 mm so any rays travelling through the volume shorter than 13 mm
> won't cross the 32 mm diameter cylinderical object region (except for
> the two ends which is not of interest).
> To leave some margin I set a threshold of 7 mm. See the attached
> picture for the results of one SART iteration.
> The left one is with the default threshold. You can see dark and
> bright dots at the corners and some streaks coming from the topleft
> corner.
> The right one is with 7 mm threshold and the slice is clean except for
> a trace of a circle outside which is easy to remove afterwards.
> So this works.
>
> I don't think that incresing the volume and cropping it in the end
> will simply work unless the enlarged volume's projection is bigger
> than the detector image; becasue the problematic values are not only
> at edge and corner voxels but are also spread in the volume as streaks
> by the forward projector as shown in the left picture.
>
> I believe that OS-SART and SIRT can mitigate this problem too since
> they are less sensitive to noise, although they are slower.
>
> I will move to CG once I have a good SART implementation for the big
> datasets in my group. There are still a lot of challenges to me.
> Unlike in FDK you can reconstruct a small subvolume directly, with
> iterative methods (I believe) I have to always reconstruct full slices
> which results in memory issues especially with CUDA. I need to stream
> the reconstruction pipeline somehow...
>
> Best regards,
> Chao
>
> 2018-02-21 13:38 GMT+01:00 Cyril Mory <cyril.mory at creatis.insa-lyon.fr
> <mailto:cyril.mory at creatis.insa-lyon.fr>>:
>
> Hi Chao,
>
> Indeed, you identified the problem quite well. That division is
> required from the maths of SART, but it brings its set of
> problems. To make a long story short, I don't know of any best
> practice in order to solve this problem. My suggestions:
>
> - increasing the threshold to the size of a few voxels could do
> the trick. We've never tried it, and I'm curious about the result
>
> - increasing the size of your volume, if you can, and cropping it
> in the end, is also a good idea, and could work, but it would
> increase the memory and time requirements, so I'd try it only if
> the rest fails
>
> - the theoretical origin of these artifacts is that in SART,
> projections are back-projected one by one instead of all together,
> so when its turn comes, each projection can have a strong
> influence on the volume. Try the --nprojspersubset argument. I've
> explained its role in details in an earlier email,
> https://public.kitware.com/pipermail/rtk-users/2017-July/010470.html
> <https://public.kitware.com/pipermail/rtk-users/2017-July/010470.html>,
> but the email doesn't display correctly, so I'm copy-pasting it
> below between <<<<<< >>>>>>>.
>
> - use conjugate gradient instead, removing the lambda and
> increasing the number of iterations (at least 30). CG requires
> more iterations, but each iteration is shorter, and it can run
> fully on GPU (switch --cudacg on if your GPU has enough memory,
> off otherwise).
>
> Please keep us posted with the results of your experiments,
>
> Cyril
>
> <<<<<<
>
> Hi Lotte,
>
>
> I'm on vacation, with very limited access to the Internet, so I
> can't look at your SIRT result, but I can answer your question on
> SART, SIRT and CG : all of those (as well as ART, and another
> method called OS-SART) minimize the same cost function, which only
> consists of a least-squares data-attachment term, i.e. || R f - p
> ||^2, with f the sought volume, p the projections and R the
> forward projection, but with different algorithms :
> - SIRT does a simple gradient descent. Since the gradient of the
> cost function is 2 R* ( R f - p ), with R* the transpose of R,
> i.e. the back projection, this means that at each iteration, the
> algorithm needs one forward and one back projection from ALL
> angles, and one "update" of the volume
> - ART, SART and OS-SART all use the same strategy: they split the
> cost function into smaller bits (individual rays for ART,
> individual projections for SART, sets of several projections for
> OS-SART, so ART splits the most, and SART the least), and
> alternately minimize the cost for each bit. We count one iteration
> when each of the smaller bits has triggered an "update" of the
> volume. This means that, per iteration, the smaller you split, the
> more updates of the volume the algorithm performs, so the faster
> (in terms of number of iterations) you get to convergence.
> Obviously it does have a dangerous drawback: if data is
> inconsistent (noise, scatter, truncation, ...), such strategies
> may not converge
> - Conjugate gradient minimizes the same cost function, without
> splitting it (so like SIRT), but using the conjugate gradient
> algorithm, which converges faster than a simple gradient descent,
> for two reasons : first, the step size is calculated analytically
> at each iteration and is optimal, and second, the descent
> direction is a combination of the gradient at the current
> iteration and the descent direction at the previous iteration (a
> "conjugate" direction, thus the algorithm's name)
>
> Hope it helps,
> Cyril
> >>>>>>
>
>
>
>
> On 21/02/2018 12:57, Chao Wu wrote:
>> L.S.,
>>
>> I was working on FDK in the past and interative reconstruction
>> methods are still new to me.
>> I understand the concept of iteratvie methods but are not aware
>> of technical details in implementation.
>>
>> Recently I am trying SART but got streak artefacts in
>> reconstructed slices, as well as dots with very high value (both
>> negative and positive) at corners of slices.
>> When I checked intermediate images in the pipleline I found that
>> those are introduced in itk::DivideOrZeroOutImageFilter.
>> You can see from the attached picture: the left half shows the
>> output of rtk::RayBoxIntersectionImageFilter and the right half
>> the output of itk::DivideOrZeroOutImageFilter, both during
>> processing of the first projection in the first iteration.
>> Apparently, although it contains the whole object, my volume is
>> relatively small compared to the size of the detector images.
>> Then the rays intersecting the volume near corners and edges
>> result in small values in the output of the raybox filter, and
>> subsequently magnify the pixel values largely after division.
>> This may not be a problem if the detector images are noiseless,
>> but in practice this will magnify the noise and they will stay as
>> streaks and dots in slices.
>>
>> To correct for this I have something in mind, such as making the
>> volume bigger and cropping the detector images so that corners
>> and edges of the volume do not project to the cropped detector;
>> or increasing the threshold in the divide filter so that low
>> values from edge/corner rays wll be zero out. Since I am lack of
>> experiences in interative methods, my question is what the best
>> or common practice will be to handle this? Thanks a lot.
>>
>> Regards,
>> Chao
>>
>>
>>
>>
>> _______________________________________________
>> Rtk-users mailing list
>> Rtk-users at public.kitware.com <mailto:Rtk-users at public.kitware.com>
>> https://public.kitware.com/mailman/listinfo/rtk-users
>> <https://public.kitware.com/mailman/listinfo/rtk-users>
>
>
> _______________________________________________
> Rtk-users mailing list
> Rtk-users at public.kitware.com <mailto:Rtk-users at public.kitware.com>
> https://public.kitware.com/mailman/listinfo/rtk-users
> <https://public.kitware.com/mailman/listinfo/rtk-users>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://public.kitware.com/pipermail/rtk-users/attachments/20180221/d849aa6e/attachment.html>
More information about the Rtk-users
mailing list