<html>
<head>
<meta content="text/html; charset=utf-8" http-equiv="Content-Type">
</head>
<body bgcolor="#FFFFFF" text="#000000">
<p>Hi Vahid,</p>
<p>It's becoming unclear to me, but I don't think you want to
perform a dot product between a matrix and a vector. My advice: do
the math, turn the equations into a pipeline (that's the tricky
part) and try to copy some pieces of existing pipelines in RTK
filters. And I don't recommend the rtkConjugateGradientFilter for
this.<br>
</p>
<p>Then run your pipeline with known data, and at every step where
you have a known reference (e.g. from python code), compare the
intermediate image you get with that reference. <br>
</p>
<p>As for the python wrapping: I personally have no experience with
writing python wrappings for RTK. It is supposed to be easily
done, but only for full RTK filters. You cannot wrap part of a
filter. So you will have to create your own filter, and only then
wrap it in python. If I'm mistaken, please, RTK users, do correct
me :)</p>
<p>Regards,</p>
<p>Cyril<br>
</p>
<div class="moz-cite-prefix">On 11/03/2016 09:10 PM, vahid ettehadi
wrote:<br>
</div>
<blockquote
cite="mid:1775591917.651990.1478203831679@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_1478202307292_6193">OK. good. Thank
Cyril to clarify it to me. To implement the Newton's methods
instead of multiplication of R with R^T, we need the
multiplication of Jacobian/Sensitivity (R here) matrix with
vector variable (here x) and the gradient vector (here
multiplication of R with -r). I think it is possible to
implement the Newton optimization methods with the current
modules as I think we could extract the R and f dot product
and R and r dot product. Am I right?</div>
<div id="yui_3_16_0_ym19_1_1478202307292_6193" dir="ltr">About
the displayed algorithms, I think the gradient vector is dot
product of matrix A (Jacobian) with the minus residual vector
(-r). </div>
<div id="yui_3_16_0_ym19_1_1478202307292_6193" dir="ltr"><br>
</div>
<div id="yui_3_16_0_ym19_1_1478202307292_6193" dir="ltr">One
more question:</div>
<div id="yui_3_16_0_ym19_1_1478202307292_6193" dir="ltr">I
already developed my codes in python. Do you think it is easy
to wrap your modules in python? I don't have experience in
these subject.</div>
<div id="yui_3_16_0_ym19_1_1478202307292_6193" dir="ltr"><br>
</div>
<div id="yui_3_16_0_ym19_1_1478202307292_6193" dir="ltr">Regards,</div>
<div id="yui_3_16_0_ym19_1_1478202307292_6193" dir="ltr">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 Thursday,
November 3, 2016 2:23 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="yiv0060408286">
<div>
<div>Hello Vahid,</div>
<div>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 clear="none">
</div>
<div> 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):</div>
<div><img moz-do-not-send="true"
class="yiv0060408286mwe-math-fallback-image-inline"
src="cid:2ELmSQtDCBHEiAzhdYFX"
style="vertical-align:-23.598ex;margin-bottom:-0.24ex;width:45.967ex;height:48.843ex;"></div>
<div>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 clear="none">
</div>
<div>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.</div>
<div>Best,</div>
<div>Cyril<br clear="none">
</div>
<br clear="none">
<div class="yiv0060408286yqt7446797867"
id="yiv0060408286yqt19943">
<div class="yiv0060408286moz-cite-prefix">On
11/03/2016 03:38 AM, vahid ettehadi 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="yiv0060408286yui_3_16_0_ym19_1_1478133746044_18490"><span
id="yiv0060408286yui_3_16_0_ym19_1_1478133746044_18491">Hello Simon and
Cyril,</span></div>
<div
id="yiv0060408286yui_3_16_0_ym19_1_1478133746044_18492"><span
id="yiv0060408286yui_3_16_0_ym19_1_1478133746044_18493">Thanks for the
reply.</span></div>
<div dir="ltr"
id="yiv0060408286yui_3_16_0_ym19_1_1478133746044_18494"><span
id="yiv0060408286yui_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="yiv0060408286yui_3_16_0_ym19_1_1478133746044_18494">Thanks
again for the reply.</div>
<div dir="ltr"
id="yiv0060408286yui_3_16_0_ym19_1_1478133746044_18494"><br
clear="none">
</div>
<div dir="ltr"
id="yiv0060408286yui_3_16_0_ym19_1_1478133746044_18494">@Cyril:</div>
<div dir="ltr"
id="yiv0060408286yui_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="yiv0060408286yui_3_16_0_ym19_1_1478133746044_18494"><br
clear="none">
</div>
<div dir="ltr"
id="yiv0060408286yui_3_16_0_ym19_1_1478133746044_18494">Regards,</div>
<div dir="ltr"
id="yiv0060408286yui_3_16_0_ym19_1_1478133746044_18494">Vahid</div>
<div class="yiv0060408286qtdSeparateBR"><br
clear="none">
<br clear="none">
</div>
<div class="yiv0060408286yahoo_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
moz-do-not-send="true"
rel="nofollow" shape="rect"
class="yiv0060408286moz-txt-link-rfc2396E"
ymailto="mailto:cyril.mory@creatis.insa-lyon.fr" target="_blank"
href="mailto:cyril.mory@creatis.insa-lyon.fr"><cyril.mory@creatis.insa-lyon.fr></a>
wrote:<br clear="none">
</font></div>
<br clear="none">
<br clear="none">
<div
class="yiv0060408286y_msg_container">
<div id="yiv0060408286">
<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="yiv0060408286moz-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="yiv0060408286yqt2869888754"
id="yiv0060408286yqt14261">
<div
class="yiv0060408286moz-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="yiv0060408286yui_3_16_0_ym19_1_1478018323282_15982">Hello
RTK users and developers,</div>
<div
id="yiv0060408286yui_3_16_0_ym19_1_1478018323282_15982"><br
clear="none">
</div>
<div
id="yiv0060408286yui_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="yiv0060408286yui_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="yiv0060408286yui_3_16_0_ym19_1_1478018323282_17709"><br
clear="none">
</div>
<div dir="ltr"
id="yiv0060408286yui_3_16_0_ym19_1_1478018323282_17709">Best
Regards,</div>
<div dir="ltr"
id="yiv0060408286yui_3_16_0_ym19_1_1478018323282_17709">Vahid</div>
<div dir="ltr"
id="yiv0060408286yui_3_16_0_ym19_1_1478018323282_17710"><br
id="yiv0060408286yui_3_16_0_ym19_1_1478018323282_17711" clear="none">
</div>
</div>
<br clear="none">
<fieldset
class="yiv0060408286mimeAttachmentHeader"></fieldset>
<br clear="none">
<pre>_______________________________________________
Rtk-users mailing list
<a moz-do-not-send="true" rel="nofollow" shape="rect" class="yiv0060408286moz-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="yiv0060408286moz-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 clear="none">
<div
class="yiv0060408286yqt2869888754"
id="yiv0060408286yqt28907">_______________________________________________<br
clear="none">
Rtk-users mailing list<br
clear="none">
<a moz-do-not-send="true"
rel="nofollow" shape="rect"
ymailto="mailto:Rtk-users@public.kitware.com"
target="_blank"
href="mailto:Rtk-users@public.kitware.com">Rtk-users@public.kitware.com</a><br
clear="none">
<a moz-do-not-send="true"
rel="nofollow" shape="rect"
target="_blank"
href="http://public.kitware.com/mailman/listinfo/rtk-users">http://public.kitware.com/mailman/listinfo/rtk-users</a><br
clear="none">
</div>
<br clear="none">
<br clear="none">
</div>
</div>
</div>
</div>
</div>
</blockquote>
</div>
<br clear="none">
</div>
</div>
<br>
<div class="yqt7446797867" id="yqt43941">_______________________________________________<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>