CCPi Regularisers do not work with SIRF ImageData

[paskino]

From @mehrhardt

Here is "another" issue. x_CIL and x_SIRF in the end should be the same, not equal to zero and a "denoised" version of "image". They are not equal: x_SIRF = image and x_CIL is equal to zero. Any idea what is going wrong here?

Below it follows a full example, but before a little explanation:

To access the CCPi regularisers within the CCPi-Framework for optimisation and image reconstruction, we wrap the regularisers in so-called plugins. Although, this structure seems a bit involved, it allows to keep the regularisers and the CCPi-Framework independent. In particular, the regularisers work on numpy arrays without any notion of ImageData or DataContainer.

The role of the plugin is, therefore, to get a DataContainer or sibling and pass the relevant data to the regulariser and return the appropriate DataContainer.

However, the CCPi-FrameworkPlugins for the regularisers are CCPi-Framework specific so they cannot be simply used with SIRF DataContainers. @paskino proposed a trick to do this. One would create a class FGP_TV_SIRF inheriting from the CCPi plugin FGP_TV, call the prox method of the parent class of the plugin, i.e. FGP_TV which extracts the data from the DataContainer with the as_array() method'; get the result and put it back into a proper SIRF DataContainer:

class FGP_TV_SIRF(FGP_TV):
    def prox(self, x, sigma):
       out = super(FGP_TV, self).prox(x, sigma)
       y = x.copy()
       y.fill(out.as_array())
       return y

Full example which fails.

import numpy
import matplotlib.pyplot as plt
import os
import shutil
import pSTIR as pet
os.chdir(pet.petmr_data_path('pet'))
shutil.rmtree('working_folder/thorax_single_slice',True)
shutil.copytree('thorax_single_slice','working_folder/thorax_single_slice')
os.chdir('working_folder/thorax_single_slice')

image = pet.ImageData('emission.hv');
image_array=image.as_array()*.05
image.fill(image_array);
from plugins.regularisers import FGP_TV

class FGP_TV_SIRF(FGP_TV):
    def prox(self, x, sigma):
       print("calling FGP")
       out = super(FGP_TV, self).prox(x, sigma)
       y = x.copy()
       y.fill(out.as_array())
       return y

g_reg = FGP_TV_SIRF(lambdaReg=1,
                iterationsTV=5000,
                tolerance=1e-16,
                methodTV=0,
                nonnegativity=1,
                printing=0,
                device='cpu')

g = FGP_TV(lambdaReg=1,
                iterationsTV=5000,
                tolerance=1e-16,
                methodTV=0,
                nonnegativity=1,
                printing=0,
                device='cpu')

x_SIRF = g_reg.prox(image, 2)
x_CIL = g.prox(image, 2)

(x_SIRF - image).norm()
numpy.linalg.norm(x_CIL.as_array() - image.as_array())
numpy.linalg.norm(image.as_array())
numpy.linalg.norm(x_CIL.as_array())

[paskino]

x_CIL is zero

This is fixed by https://github.com/vais-ral/CCPi-Regularisation-Toolkit/pull/62

And additionally by https://github.com/vais-ral/CCPi-Framework/pull/152 for Python 2.7

[paskino]

x_SIRF = image

is fixed by correcting the use of super

class FGP_TV_SIRF(FGP_TV):
    def prox(self, x, sigma):
       out = super(FGP_TV_SIRF, self).prox(x, sigma)
       y = x.copy()
       y.fill(out.as_array())
       return y

In the previous implementation the class wasn't doing what expected and described in the issue. I'm still puzzled that it didn't get into a recursion loop yielding a RuntimeError: maximum recursion depth exceeded

[mehrhardt]

This is interesting. Before with out = super(FGP_TV, self).prox(x, sigma) we called the parent of FGP_TV which is Function and we were surprised about its behavior. This actually reveals dangerous coding currently being done in https://github.com/vais-ral/CCPi-Framework/blob/master/Wrappers/Python/ccpi/optimisation/funcs.py where the prox (and other functions) has a default behavior that makes no sense for a generic Function. Perhaps instead of doing the "identity", it should throw an error as one should not be allowed to call this if not overwritten.