PDHG algorithm for PET recon with SIRF

[epapoutsellis]

Demo that compares the PET reconstruction with PDHG from CIL using data fidelity:

1. KullbackLeibler from CIL
   1. pet.PoissonLogLikelihoodWithLinearModelForMeanAndProjData

Hacks to make PDHG from CIL to work with SIRF

1. We need to compute the norm of the AcquisitionModel. Currently it uses the PowerMethod from the LinearOperator class from CIL, however we suggest to add the PowerMethod in the AcquisitionModel class.

2. For the PoissonLogLikelihoodWithLinearModelForMeanAndProjData, we need

   1. convex_conjugate method to compute the convex conjugate of PoissonLogLikelihoodWithLinearModelForMeanAndProjData. In CIL, we are using scipy which is the convex conjugate of scipy.special.kl_div

   2. compute the proximal of the convex conjugate

      1. the __call__ method for PoissonLogLikelihoodWithLinearModelForMeanAndProjData. Currently, we have this method. In the STIR KullbackLeibler I believe that get_value computes

\sum \overline{y_b} - y_{b}\log\overline{y_b}     ( 1 )

and not

KL = sum \overline{y_b} - y_{b} log \overline{y_b}  - y_{b} + y_{b} log y_{b}  ( 2 )

3) From PoissonLogLikelihoodWithLinearModelForMeanAndProjData we would like

get_background_term get_acquisition_data

so we can avoid writing these 2 lines

[KrisThielemans]

A few comments:

• adding PowerMethod in the AcquisitionModel class seems weird (a generic method inside an acquisition model?) I guess that's not what you meant.

• convex_conjugate and its proximal I don't know really what it'd involve. Maybe someone can elaborate on what is actually needed? Probably best in separate issue?

• adding __call__ seems trivial to do in the add_cil branch.

• point 3. please create a separate issue for this such that we can comment on where the appropriate place would be.

[paskino]

@KrisThielemans we have the PowerMethod as static method of the LinearOperator class, for this reason we thought to add it to the AcquisitionModel. However, if it seems not right, we can put it in SIRF.py as function, maybe?

[KrisThielemans]

done in CIL, with some remaining issues listed separately already. #804