[Rtk-users] Model-based image reconstruction based on the RTK modules

[vahid ettehadi]

Hello Simon and Cyril,Thanks for the reply.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 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.Thanks again for the reply.
@Cyril:Please correct me if I am wrong. you mean the output of backProjectionFilter is the gradient of defined cost function?
Regards,Vahid 

    On Wednesday, November 2, 2016 2:53 AM, Cyril Mory <cyril.mory at creatis.insa-lyon.fr> wrote:
 

  Hi Vahid, Welcome to RTK :) 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. 
 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)
 
 digraph SARTConeBeamReconstructionFilter {
 
 Input0 [ label="Input 0 (Volume)"];
 Input0 [shape=Mdiamond];
 Input1 [label="Input 1 (Projections)"];
 Input1 [shape=Mdiamond];
 
 node [shape=box];
 ForwardProject [ label="rtk::ForwardProjectionImageFilter" URL="\ref rtk::ForwardProjectionImageFilter"];
 Extract [ label="itk::ExtractImageFilter" URL="\ref itk::ExtractImageFilter"];
 MultiplyByZero [ label="itk::MultiplyImageFilter (by zero)" URL="\ref itk::MultiplyImageFilter"];
 AfterExtract [label="", fixedsize="false", width=0, height=0, shape=none];
 Subtract [ label="itk::SubtractImageFilter" URL="\ref itk::SubtractImageFilter"];
 MultiplyByLambda [ label="itk::MultiplyImageFilter (by lambda)" URL="\ref itk::MultiplyImageFilter"];
 Divide [ label="itk::DivideOrZeroOutImageFilter" URL="\ref itk::DivideOrZeroOutImageFilter"];
 GatingWeight [ label="itk::MultiplyImageFilter (by gating weight)" URL="\ref itk::MultiplyImageFilter", style=dashed];
 Displaced [ label="rtk::DisplacedDetectorImageFilter" URL="\ref rtk::DisplacedDetectorImageFilter"];
 ConstantProjectionStack [ label="rtk::ConstantImageSource" URL="\ref rtk::ConstantImageSource"];
 ExtractConstantProjection [ label="itk::ExtractImageFilter" URL="\ref itk::ExtractImageFilter"];
 RayBox [ label="rtk::RayBoxIntersectionImageFilter" URL="\ref rtk::RayBoxIntersectionImageFilter"];
 ConstantVolume [ label="rtk::ConstantImageSource" URL="\ref rtk::ConstantImageSource"];
 BackProjection [ label="rtk::BackProjectionImageFilter" URL="\ref rtk::BackProjectionImageFilter"];
 OutofInput0 [label="", fixedsize="false", width=0, height=0, shape=none];
 OutofBP [label="", fixedsize="false", width=0, height=0, shape=none];
 BeforeBP [label="", fixedsize="false", width=0, height=0, shape=none];
 BeforeAdd [label="", fixedsize="false", width=0, height=0, shape=none];
 Input0 -> OutofInput0 [arrowhead=none];
 OutofInput0 -> ForwardProject;
 ConstantVolume -> BeforeBP [arrowhead=none];
 BeforeBP -> BackProjection;
 Extract -> AfterExtract[arrowhead=none];
 AfterExtract -> MultiplyByZero;
 AfterExtract -> Subtract;
 MultiplyByZero -> ForwardProject;
 Input1 -> Extract;
 ForwardProject -> Subtract;
 Subtract -> MultiplyByLambda;
 MultiplyByLambda -> Divide;
 Divide -> GatingWeight;
 GatingWeight -> Displaced;
 ConstantProjectionStack -> ExtractConstantProjection;
 ExtractConstantProjection -> RayBox;
 RayBox -> Divide;
 Displaced -> BackProjection;
 BackProjection -> OutofBP [arrowhead=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. 
 
 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 https://itk.org/Doxygen/html/classitk_1_1CostFunctionTemplate.html if you want to try this approach.
 
 Best regards,
 Cyril
 
 On 11/01/2016 05:50 PM, vahid ettehadi via Rtk-users wrote:
  
  Hello RTK users and developers, 
  I already implemented the RTK and reconstructed some images with the FDK algorithm implemented in RTK. It works well. Thanks to RTK developers.
  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.  
  Best Regards, Vahid 
   
  
 _______________________________________________
Rtk-users mailing list
Rtk-users at public.kitware.com
http://public.kitware.com/mailman/listinfo/rtk-users
 
 
 
_______________________________________________
Rtk-users mailing list
Rtk-users at public.kitware.com
http://public.kitware.com/mailman/listinfo/rtk-users


   
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/rtk-users/attachments/20161103/a7ddd25a/attachment-0010.html>