Hello Vahid,
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.
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):
{\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}}
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.
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.
Best,
Cyril
On 11/03/2016 03:38 AM, vahid ettehadi wrote:
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
<[email protected]> 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
[email protected] <mailto:[email protected]>
http://public.kitware.com/mailman/listinfo/rtk-users
_______________________________________________
Rtk-users mailing list
[email protected] <mailto:[email protected]>
http://public.kitware.com/mailman/listinfo/rtk-users
_______________________________________________
Rtk-users mailing list
[email protected]
http://public.kitware.com/mailman/listinfo/rtk-users
