<html><head></head><body><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 size="2" face="Arial"> On Thursday, November 3, 2016 2:23 AM, Cyril Mory <cyril.mory@creatis.insa-lyon.fr> 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 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 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 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 clear="none" id="yiv0060408286yui_3_16_0_ym19_1_1478018323282_17711">
                          </div>
                        </div>
                        <br clear="none">
                        <fieldset class="yiv0060408286mimeAttachmentHeader"></fieldset>
                        <br clear="none">
                        <pre>_______________________________________________
Rtk-users mailing list
<a 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 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 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 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 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 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></body></html>