<div dir="ltr">Hi Cyril,<div><br></div><div>Thanks for your suggestion. </div><div><br></div><div>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).</div><div>To leave some margin I set a threshold of 7 mm. See the attached picture for the results of one SART iteration. </div><div>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.</div><div>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.</div><div>So this works.</div><div><br></div><div>I don't think that incresing the volume and cropping it in the end will
<span style="color:rgb(34,34,34);font-family:arial,sans-serif;font-size:small;font-style:normal;font-variant-ligatures:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;background-color:rgb(255,255,255);text-decoration-style:initial;text-decoration-color:initial;float:none;display:inline">simply</span>
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.</div><div><br></div><div>I believe that OS-SART and SIRT can mitigate this problem too since they are less sensitive to noise, although they are slower.</div><div><br></div><div>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...</div><div><br></div><div>Best regards,</div><div>Chao</div></div><div class="gmail_extra"><br><div class="gmail_quote">2018-02-21 13:38 GMT+01:00 Cyril Mory <span dir="ltr"><<a href="mailto:cyril.mory@creatis.insa-lyon.fr" target="_blank">cyril.mory@creatis.insa-lyon.fr</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div text="#000000" bgcolor="#FFFFFF">
<p>Hi Chao,</p>
<p>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:</p>
<p>- 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</p>
<p>- 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</p>
<p>- 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,
<a class="m_-4283454247097413405moz-txt-link-freetext" href="https://public.kitware.com/pipermail/rtk-users/2017-July/010470.html" target="_blank">https://public.kitware.com/<wbr>pipermail/rtk-users/2017-July/<wbr>010470.html</a>,
but the email doesn't display correctly, so I'm copy-pasting it
below between <<<<<<
>>>>>>>.<br>
</p>
<p>- 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).</p>
<p>Please keep us posted with the results of your experiments,</p>
<p>Cyril<br>
</p>
<p><<<<<<<br>
</p>
<p>Hi Lotte,</p>
<div dir="auto"><br>
</div>
<div dir="auto">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 :</div>
<div dir="auto">- 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 </div>
<div dir="auto">- 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</div>
<div dir="auto">- 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)</div>
<div dir="auto"><br>
</div>
<div dir="auto">Hope it helps,</div>
<div dir="auto">Cyril <br>
>>>>>><br>
<br>
<br>
<br>
</div><div><div class="h5">
<br>
<div class="m_-4283454247097413405moz-cite-prefix">On 21/02/2018 12:57, Chao Wu wrote:<br>
</div>
</div></div><blockquote type="cite"><div><div class="h5">
<div dir="ltr">
<div class="gmail_quote">
<div dir="ltr">L.S.,
<div><br>
</div>
<div>I was working on FDK in the past and interative
reconstruction methods are still new to me.</div>
<div>I understand the concept of iteratvie methods but are
not aware of technical details in implementation.</div>
<div><br>
</div>
<div>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.</div>
<div>When I checked intermediate images in the pipleline I
found that those are introduced in itk::DivideOrZeroOutImageFilte<wbr>r.</div>
<div>You can see from the attached picture: the left half
shows the output of rtk::RayBoxIntersectionImag<wbr>eFilter
and the right half the output of itk::DivideOrZeroOutImageFilte<wbr>r,
both during processing of the first projection in the
first iteration.</div>
<div>Apparently, although it contains the whole object, my
volume is relatively small compared to the size of the
detector images.</div>
<div>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
<span style="color:rgb(34,34,34);font-family:arial,sans-serif;font-size:small;font-style:normal;font-variant-ligatures:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;background-color:rgb(255,255,255);text-decoration-style:initial;text-decoration-color:initial;float:none;display:inline">largely</span>
after division.</div>
<div>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.</div>
<div><br>
</div>
<div>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.</div>
<div><br>
</div>
<div>Regards,</div>
<div>Chao</div>
<div><br>
</div>
</div>
</div>
<br>
</div>
<br>
<fieldset class="m_-4283454247097413405mimeAttachmentHeader"></fieldset>
<br>
</div></div><pre>______________________________<wbr>_________________
Rtk-users mailing list
<a class="m_-4283454247097413405moz-txt-link-abbreviated" href="mailto:Rtk-users@public.kitware.com" target="_blank">Rtk-users@public.kitware.com</a>
<a class="m_-4283454247097413405moz-txt-link-freetext" href="https://public.kitware.com/mailman/listinfo/rtk-users" target="_blank">https://public.kitware.com/<wbr>mailman/listinfo/rtk-users</a>
</pre>
</blockquote>
<br>
</div>
<br>______________________________<wbr>_________________<br>
Rtk-users mailing list<br>
<a href="mailto:Rtk-users@public.kitware.com">Rtk-users@public.kitware.com</a><br>
<a href="https://public.kitware.com/mailman/listinfo/rtk-users" rel="noreferrer" target="_blank">https://public.kitware.com/<wbr>mailman/listinfo/rtk-users</a><br>
<br></blockquote></div><br></div>