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.m...@creatis.insa-lyon.fr
<mailto:cyril.m...@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@public.kitware.com <mailto:Rtk-users@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@public.kitware.com <mailto:Rtk-users@public.kitware.com>
https://public.kitware.com/mailman/listinfo/rtk-users
<https://public.kitware.com/mailman/listinfo/rtk-users>
