[Rtk-users] Fwd: High pixel/voxel values in SART

[Chao Wu]

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>:

> 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, 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 listRtk-users at public.kitware.comhttps://public.kitware.com/mailman/listinfo/rtk-users
>
>
>
> _______________________________________________
> Rtk-users mailing list
> Rtk-users at public.kitware.com
> https://public.kitware.com/mailman/listinfo/rtk-users
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20180221/54196d85/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: sart_th0_7.png
Type: image/png
Size: 351349 bytes
Desc: not available
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20180221/54196d85/attachment.png>