<html>
<head>
<meta content="text/html; charset=windows-1252"
http-equiv="Content-Type">
</head>
<body bgcolor="#FFFFFF" text="#000000">
<p>Hello Vahid,</p>
<p>Thank you for this insight on Newton's methods. Yes, the output
of the backprojection filter, in SART, is the gradient of the cost
function. You may have noticed that the pipeline performs a
division of the forward projection by something coming from
"RayBoxIntersectionFilter". This is to normalize the forward
projection, so that R^T R f ~= blurry f. If you don't do it,
you'll have R^T R f ~= alpha * blurry f, with alpha that can reach
200 or so, and your algorithm will quickly diverge.<br>
</p>
<p> You could try to extract the gradient from the conjugate
gradient filter as well, but it is significantly more tricky:
first, the CG filter was implemented from the following algorithm,
taken from wikipedia (which embeds the normalization I was
mentioning earlier):</p>
<p><img src="cid:part1.00987626.E3135833@creatis.insa-lyon.fr"
class="mwe-math-fallback-image-inline" aria-hidden="true"
style="vertical-align: -23.598ex; margin-bottom: -0.24ex;
width:45.967ex; height:48.843ex;"
alt="{\begin{aligned}&\mathbf {r} _{0}:=\mathbf {b} -\mathbf
{Ax} _{0}\\&\mathbf {p} _{0}:=\mathbf {r}
_{0}\\&k:=0\\&{\hbox{repeat}}\\&\qquad \alpha
_{k}:={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r}
_{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {Ap}
_{k}}}\\&\qquad \mathbf {x} _{k+1}:=\mathbf {x} _{k}+\alpha
_{k}\mathbf {p} _{k}\\&\qquad \mathbf {r} _{k+1}:=\mathbf
{r} _{k}-\alpha _{k}\mathbf {Ap} _{k}\\&\qquad {\hbox{if
}}r_{k+1}{\hbox{ is sufficiently small then exit
loop}}\\&\qquad \beta _{k}:={\frac {\mathbf {r}
_{k+1}^{\mathsf {T}}\mathbf {r} _{k+1}}{\mathbf {r}
_{k}^{\mathsf {T}}\mathbf {r} _{k}}}\\&\qquad \mathbf {p}
_{k+1}:=\mathbf {r} _{k+1}+\beta _{k}\mathbf {p}
_{k}\\&\qquad k:=k+1\\&{\hbox{end
repeat}}\\&{\hbox{The result is }}\mathbf {x}
_{k+1}\end{aligned}}"></p>
<p>In this algorithm, it is not clear to me what exactly is the
gradient of the cost function. I would say it is something like "-
r_k", but I'm not sure. Second, as you see, the CG filter needs an
operator A, which may differ from one problem to another, so this
operator is implemented in a separate filter, which in your case
would be rtkReconstructionConjugateGradientOperator, with the
laplacian regularization parameter gamma set to 0. <br>
</p>
<p>Note that we never actually store the system matrix R. Instead,
the interpolation coefficient it contains are re-computed on the
fly everytime we forward project. And the same holds for
backprojection, i.e the matrix R^T.</p>
<p>Best,</p>
<p>Cyril<br>
</p>
<br>
<div class="moz-cite-prefix">On 11/03/2016 03:38 AM, vahid ettehadi
wrote:<br>
</div>
<blockquote cite="mid:24189134.77160.1478140718379@mail.yahoo.com"
type="cite">
<div style="color:#000; background-color:#fff;
font-family:HelveticaNeue, Helvetica Neue, Helvetica, Arial,
Lucida Grande, sans-serif;font-size:14px">
<div id="yui_3_16_0_ym19_1_1478133746044_18490"><span
id="yui_3_16_0_ym19_1_1478133746044_18491">Hello Simon and
Cyril,</span></div>
<div id="yui_3_16_0_ym19_1_1478133746044_18492"><span
id="yui_3_16_0_ym19_1_1478133746044_18493">Thanks for the
reply.</span></div>
<div dir="ltr" id="yui_3_16_0_ym19_1_1478133746044_18494"><span
id="yui_3_16_0_ym19_1_1478133746044_18495">You are right
Simon. I did not notice it too in the literature. The main
problem as you said is the storage. Actually I developed
the conjugate gradient (CG), quasi-Newton and Newton
optimization methods for optical tomography and I intended
to apply them to the CT reconstruction as well. I
implemented the Newton's methods (Gauss-Newton and
Levenberg-Marquardt) in a Jacobian-Free-Newton-Krylov
approaches to avoid the matrix multiplication of Jacobians
(sensitivity). It means we only need to store the Jacobian
matrix for the these methods (the matrix R that Cyril was
mentioned), that is still a big matrix for practical
problems in CT reconstruction. For the quasi-Newton I
adapted an L-BFGS algorithm that only need the 3 or 8 last
iterations of the gradient vector to calculate the Hessian
matrix. In my case, the L-BFGS and Newton's methods was much
faster than </span>the CG as you know because of using the
second order derivative (hessian matrix). I saw in your last
paper you implement the conjugate gradient method, so I
thought it might be easy to extract the gradient vector from
CG modules and solve the cost function within the
quasi-Newton/Newton methods. I will look at the codes to see
what I can do.</div>
<div dir="ltr" id="yui_3_16_0_ym19_1_1478133746044_18494">Thanks
again for the reply.</div>
<div dir="ltr" id="yui_3_16_0_ym19_1_1478133746044_18494"><br>
</div>
<div dir="ltr" id="yui_3_16_0_ym19_1_1478133746044_18494">@Cyril:</div>
<div dir="ltr" id="yui_3_16_0_ym19_1_1478133746044_18494">Please
correct me if I am wrong. you mean the output of
backProjectionFilter is the gradient of defined cost function?</div>
<div dir="ltr" id="yui_3_16_0_ym19_1_1478133746044_18494"><br>
</div>
<div dir="ltr" id="yui_3_16_0_ym19_1_1478133746044_18494">Regards,</div>
<div dir="ltr" id="yui_3_16_0_ym19_1_1478133746044_18494">Vahid</div>
<div class="qtdSeparateBR"><br>
<br>
</div>
<div class="yahoo_quoted" style="display: block;">
<div style="font-family: HelveticaNeue, Helvetica Neue,
Helvetica, Arial, Lucida Grande, sans-serif; font-size:
14px;">
<div style="font-family: HelveticaNeue, Helvetica Neue,
Helvetica, Arial, Lucida Grande, sans-serif; font-size:
16px;">
<div dir="ltr"><font face="Arial" size="2"> On Wednesday,
November 2, 2016 2:53 AM, Cyril Mory
<a class="moz-txt-link-rfc2396E" href="mailto:cyril.mory@creatis.insa-lyon.fr"><cyril.mory@creatis.insa-lyon.fr></a> wrote:<br>
</font></div>
<br>
<br>
<div class="y_msg_container">
<div id="yiv2652688545">
<div>
<div>Hi Vahid,</div>
<div>Welcome to RTK :)</div>
Indeed, there are several iterative methods already
implemented in RTK, but none of the filters allows
you to easily extract the gradient of the least
squares function there are minimizing. <br
clear="none">
If you need to minimize the classical
non-regularized tomographic cost function, ie || R f
- p ||², with R the forward projection operator, f
the volume you are looking for, and p the measured
projections, my best advice would be to copy some
part of the pipeline of
rtkSARTConeBeamReconstructionFilter to get the job
done, ie the following part (copy-paste this into
webgraphviz.com)<br clear="none">
<br clear="none">
digraph SARTConeBeamReconstructionFilter {<br
clear="none">
<br clear="none">
Input0 [ label="Input 0 (Volume)"];<br clear="none">
Input0 [shape=Mdiamond];<br clear="none">
Input1 [label="Input 1 (Projections)"];<br
clear="none">
Input1 [shape=Mdiamond];<br clear="none">
<br clear="none">
node [shape=box];<br clear="none">
ForwardProject [
label="rtk::ForwardProjectionImageFilter" URL="\ref
rtk::ForwardProjectionImageFilter"];<br clear="none">
Extract [ label="itk::ExtractImageFilter" URL="\ref
itk::ExtractImageFilter"];<br clear="none">
MultiplyByZero [ label="itk::MultiplyImageFilter (by
zero)" URL="\ref itk::MultiplyImageFilter"];<br
clear="none">
AfterExtract [label="", fixedsize="false", width=0,
height=0, shape=none];<br clear="none">
Subtract [ label="itk::SubtractImageFilter"
URL="\ref itk::SubtractImageFilter"];<br
clear="none">
MultiplyByLambda [ label="itk::MultiplyImageFilter
(by lambda)" URL="\ref itk::MultiplyImageFilter"];<br
clear="none">
Divide [ label="itk::DivideOrZeroOutImageFilter"
URL="\ref itk::DivideOrZeroOutImageFilter"];<br
clear="none">
GatingWeight [ label="itk::MultiplyImageFilter (by
gating weight)" URL="\ref itk::MultiplyImageFilter",
style=dashed];<br clear="none">
Displaced [
label="rtk::DisplacedDetectorImageFilter" URL="\ref
rtk::DisplacedDetectorImageFilter"];<br clear="none">
ConstantProjectionStack [
label="rtk::ConstantImageSource" URL="\ref
rtk::ConstantImageSource"];<br clear="none">
ExtractConstantProjection [
label="itk::ExtractImageFilter" URL="\ref
itk::ExtractImageFilter"];<br clear="none">
RayBox [ label="rtk::RayBoxIntersectionImageFilter"
URL="\ref rtk::RayBoxIntersectionImageFilter"];<br
clear="none">
ConstantVolume [ label="rtk::ConstantImageSource"
URL="\ref rtk::ConstantImageSource"];<br
clear="none">
BackProjection [
label="rtk::BackProjectionImageFilter" URL="\ref
rtk::BackProjectionImageFilter"];<br clear="none">
OutofInput0 [label="", fixedsize="false", width=0,
height=0, shape=none];<br clear="none">
OutofBP [label="", fixedsize="false", width=0,
height=0, shape=none];<br clear="none">
BeforeBP [label="", fixedsize="false", width=0,
height=0, shape=none];<br clear="none">
BeforeAdd [label="", fixedsize="false", width=0,
height=0, shape=none];<br clear="none">
Input0 -> OutofInput0 [arrowhead=none];<br
clear="none">
OutofInput0 -> ForwardProject;<br clear="none">
ConstantVolume -> BeforeBP [arrowhead=none];<br
clear="none">
BeforeBP -> BackProjection;<br clear="none">
Extract -> AfterExtract[arrowhead=none];<br
clear="none">
AfterExtract -> MultiplyByZero;<br clear="none">
AfterExtract -> Subtract;<br clear="none">
MultiplyByZero -> ForwardProject;<br clear="none">
Input1 -> Extract;<br clear="none">
ForwardProject -> Subtract;<br clear="none">
Subtract -> MultiplyByLambda;<br clear="none">
MultiplyByLambda -> Divide;<br clear="none">
Divide -> GatingWeight;<br clear="none">
GatingWeight -> Displaced;<br clear="none">
ConstantProjectionStack ->
ExtractConstantProjection;<br clear="none">
ExtractConstantProjection -> RayBox;<br
clear="none">
RayBox -> Divide;<br clear="none">
Displaced -> BackProjection;<br clear="none">
BackProjection -> OutofBP [arrowhead=none];<br
clear="none">
}<br clear="none">
<br clear="none">
As you can see, it is a very large part of the SART
reconstruction filter, so yoiu might be better off
just copying the whole
SARTConeBeamReconstructionFilter and modifying it. <br
clear="none">
<br clear="none">
Of course, you could also look into ITK's cost
function class, and see if one of the classes
inherited from it suits your needs, implement your
cost function this way, and use ITK's off-the-shelf
solvers to minimize it. See the inheritance diagram
in <a moz-do-not-send="true" rel="nofollow"
shape="rect"
class="yiv2652688545moz-txt-link-freetext"
target="_blank"
href="https://itk.org/Doxygen/html/classitk_1_1CostFunctionTemplate.html">https://itk.org/Doxygen/html/classitk_1_1CostFunctionTemplate.html</a>
if you want to try this approach.<br clear="none">
<br clear="none">
Best regards,<br clear="none">
Cyril<br clear="none">
<br clear="none">
<div class="yiv2652688545yqt2869888754"
id="yiv2652688545yqt14261">
<div class="yiv2652688545moz-cite-prefix">On
11/01/2016 05:50 PM, vahid ettehadi via
Rtk-users wrote:<br clear="none">
</div>
<blockquote type="cite">
<div
style="color:#000;background-color:#fff;font-family:HelveticaNeue,
Helvetica Neue, Helvetica, Arial, Lucida
Grande, sans-serif;font-size:14px;">
<div
id="yiv2652688545yui_3_16_0_ym19_1_1478018323282_15982">Hello
RTK users and developers,</div>
<div
id="yiv2652688545yui_3_16_0_ym19_1_1478018323282_15982"><br
clear="none">
</div>
<div
id="yiv2652688545yui_3_16_0_ym19_1_1478018323282_15982">I
already implemented the RTK and
reconstructed some images with the FDK
algorithm implemented in RTK. It works well.
Thanks to RTK developers.<br clear="none">
</div>
<div dir="ltr"
id="yiv2652688545yui_3_16_0_ym19_1_1478018323282_17709">Now,
I am trying to develop a model-based image
reconstruction for our cone-beam micro-CT. I
see already that some iterative algorithms
like ART and its modifications and
conjugate-gradient (CG) method are
implemented in the RTK. I want to develop a
model-based reconstruction through the
Newton/quasi-Newton optimizations methods. I
was wondering is it possible to extract the
gradient of least square function from
implemented algorithms like CG module? Any
recommendation will be appreciated. </div>
<div dir="ltr"
id="yiv2652688545yui_3_16_0_ym19_1_1478018323282_17709"><br
clear="none">
</div>
<div dir="ltr"
id="yiv2652688545yui_3_16_0_ym19_1_1478018323282_17709">Best
Regards,</div>
<div dir="ltr"
id="yiv2652688545yui_3_16_0_ym19_1_1478018323282_17709">Vahid</div>
<div dir="ltr"
id="yiv2652688545yui_3_16_0_ym19_1_1478018323282_17710"><br
id="yiv2652688545yui_3_16_0_ym19_1_1478018323282_17711" clear="none">
</div>
</div>
<br clear="none">
<fieldset
class="yiv2652688545mimeAttachmentHeader"></fieldset>
<br clear="none">
<pre>_______________________________________________
Rtk-users mailing list
<a moz-do-not-send="true" rel="nofollow" shape="rect" class="yiv2652688545moz-txt-link-abbreviated" ymailto="mailto:Rtk-users@public.kitware.com" target="_blank" href="mailto:Rtk-users@public.kitware.com">Rtk-users@public.kitware.com</a>
<a moz-do-not-send="true" rel="nofollow" shape="rect" class="yiv2652688545moz-txt-link-freetext" target="_blank" href="http://public.kitware.com/mailman/listinfo/rtk-users">http://public.kitware.com/mailman/listinfo/rtk-users</a>
</pre>
</blockquote>
</div>
<br clear="none">
</div>
</div>
<br>
<div class="yqt2869888754" id="yqt28907">_______________________________________________<br
clear="none">
Rtk-users mailing list<br clear="none">
<a moz-do-not-send="true" shape="rect"
ymailto="mailto:Rtk-users@public.kitware.com"
href="mailto:Rtk-users@public.kitware.com">Rtk-users@public.kitware.com</a><br
clear="none">
<a moz-do-not-send="true" shape="rect"
href="http://public.kitware.com/mailman/listinfo/rtk-users"
target="_blank">http://public.kitware.com/mailman/listinfo/rtk-users</a><br
clear="none">
</div>
<br>
<br>
</div>
</div>
</div>
</div>
</div>
</blockquote>
<br>
</body>
</html>