Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

SPDHG algorithm from CIL #238

Merged
merged 47 commits into from
Oct 23, 2019
Merged

SPDHG algorithm from CIL #238

merged 47 commits into from
Oct 23, 2019

Conversation

paskino
Copy link
Contributor

@paskino paskino commented Nov 15, 2018

Adds an example on using the SPDHG algorithm with OS from CIL

Requires SyneRBI/Hackathon-SIRF-SuperBuild#1 which I am going to port from the hackathon repo to this.

paskino and others added 30 commits July 27, 2018 10:26
added: copy, conjugate
closes Hackathon-SIRF/#7
notice that it's not taking into consideration
SyneRBI/Hackathon-SIRF#2 (comment)
for MR images dot operator calculates the inner product with conjugate
by itself.
changed AcquisitionModel to inherit from object.
fixes an assertion error
@mehrhardt
Copy link
Contributor

The PR looks good to me.

A general question here is where in STIR/SIRF should all this code go? I.e. the algorithm spdhg, definition of KL divergence etc. Some of the files are in what was intended as a temporary folder "pCIL".

Also there are a few "hacks" to quickly code things up in examples/Python/PET/interactive/spdhgutils.py like the SubsetOperator. These are not spdhg specific and might be useful for other subset based algorithms. One of the purposes of the class OperatorInd is to "crop" the output of the STIR forward operator to a specific subset in order to minimize memory and cpu footprint of spdhg. At the moment STIR would output an array of the same size as the usual data but only with non-zeros for the actual subset. spdhg then takes this output and runs a few other operations on variables of the same size (square root, ...).

@paskino
Copy link
Contributor Author

paskino commented Nov 15, 2018

The actual algorithm is in CIL and the pCIL package should be deleted.

Also #237 should be merged first.

import pSTIR as pet
from ccpi.optimisation.spdhg import spdhg
from ccpi.optimisation.spdhg import KullbackLeibler
from ccpi.optimisation.spdhg import KullbackLeiblerConvexConjugate
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The KullbackLeiblerConvexConjugate is not needed here. It is been used implicitly through KullbackLeibler.

@@ -1202,6 +1212,10 @@ def backward(self, ad):
(self.handle, ad.handle)
check_status(image.handle)
return image
def direct(self, image):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why does CoilImageData have a method direct? This is unintuitive to me.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what about adjoint? Should it be there?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well, I don't know. To me this is a data object and not an operator. So it should neither have direct, forward nor adjoint. But it could be that I did not see any specific CCP-PETMR point of view on this. @KrisThielemans ?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

agree with Matthias - this class should have neither direct nor adjoint

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm quite confused about this. The line where this comment is pinned to (here is part of AcquisitionModel. CoilImageData doesn't have forward/direct/backward/adjoint in this PR, as it should be.

Copy link
Contributor

@mehrhardt mehrhardt left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good overall. A made a few small comments separately from the review.



def PowerMethodNonsquare(op, numiters, x0=None):
# Initialise random
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

better remove all the comments

Copy link
Member

@KrisThielemans KrisThielemans left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this needs more comments on what's going on I'm afraid.

@@ -130,6 +130,10 @@ def main():
print('norm of image*10: %f' % image.norm())
diff = image.clone() - image
print('norm of image.clone() - image: %f' % diff.norm())

# test clone vs copy
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is a test. should it be in an example? I'd rather put it in a test

g_noreg = ZeroFun()


g_reg = FGP_TV_SIRF(lambdaReg=.3,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

some doc please what this does

# Author: Matthias J Ehrhardt, Edoardo Pasca
#
## CCP PETMR Synergistic Image Reconstruction Framework (SIRF)
## Copyright 2015 - 2018 Rutherford Appleton Laboratory STFC
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this will need some change!

STIR.py:
  adds:
    * get_background_term
    * saves background term in self.bt member of AcquisitionModel
    * get_linear_acquisition_model
    * direct
    * adjoint
    * is_linear
    * is_affine
addresses #237 (comment)
added get_constant_term
added get_additive_term

           PET acquisition model that relates an image x to the
           acquisition data y as
           (F)    y = [1/n] (G x + [a]) + [b]
           where:
           G is the geometric (ray tracing) projector from the image voxels
           to the scanner's pairs of detectors (bins);
           a and b are otional additive and background terms representing
           the effects of noise and scattering;

           Returns [a] + [b]
added copy method to SIRF.DatContainer class
@evgueni-ovtchinnikov evgueni-ovtchinnikov merged commit abf11c5 into master Oct 23, 2019
@paskino paskino deleted the spdhg_from_cil branch September 3, 2020 13:39
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants