CIL Phase contribution

Changes

Add Paganin phase retrieval methods based on the description in https://onlinelibrary.wiley.com/doi/10.1046/j.1365-2818.2002.01010.x

The phase retrieval can be applied using

processor = PaganinProcessor(delta, beta, energy)
processor.set_input(data)
thickness = processor.get_output()

Runs the Paganin method, returning the sample thickness, $ T = -\frac{1}{\mu}\ln(F^{-1}\frac{F(M^2 I_{norm}(x,y,z=\Delta))}{1+\frac{\Delta\lambda\delta}{4\pi\beta}(k_x^2+k_y^2)/M})$

Or with the filter_type argument ='generalised_paganin_method' uses the filter described in https://iopscience.iop.org/article/10.1088/2040-8986/abbab9 $ T = -\frac{1}{\mu}\ln(F^{-1}\frac{F(M^2 I_{norm}(x,y,z=\Delta))}{1-\frac{2\Delta\lambda\delta}{4\pi\beta W^2}(cos(Wk_x)+cos(Wk_y)-2)/M})$

Both methods require normalised transmission data, $I_{norm}$, physical parameters for delta, beta (the real and complex part of the material refractive index) and energy. The propagation distance, $\Delta$, magnification, $M$, and pixel size, $W$, are accessed from the data geometry. If these values are not found in the geometry, the processor prints a warning and uses default values. The geometry values can be over-ridden using

processor = PaganinProcessor(delta, beta, energy)
processor.set_input(data)
processor.update_parameters({'propagation_distance': 0.1})

We can return only the filtered image using,

processor = PaganinProcessor(delta, beta, energy)
processor.set_input(data)
filtered_image = processor.get_output(full_retrieval=False)

which returns

$ I_{filtered} = F^{-1}\frac{F( I(x,y,z=\Delta))}{1+\frac{\Delta\lambda\delta}{4\pi\beta}(k_x^2+k_y^2)/M}$

Testing you performed

I've made a pull request to add an example to CIL-Demos/misc using these methods https://github.com/TomographicImaging/CIL-Demos/pull/151 Below is an output from the notebook showing a comparison of the cross-section of a test image with different phase retrieval methods. The scaled phase retrieval in CIL matches the phase retrieval method in Tomopy when performed on transmission data then converted to absorption, and the filter in CIL matches the Tomopy method performed on absorption data where the Tomopy regularisation parameter argument is tomopy_alpha = 1/(4*np.pi**2*(delta/beta))

Unit tests of the functionality are in test_DataProcessor.py TestPaganinPhaseRetriver() and TestFilter()

Related issues/links

Checklist

[x] I have performed a self-review of my code
[x] I have added docstrings in line with the guidance in the developer guide
[x] I have updated the relevant documentation
[x] I have implemented unit tests that cover any new or modified functionality
[x] CHANGELOG.md has been updated with any functionality change
[x] Request review from all relevant developers
[x] Change pull request label to 'Waiting for review'

Contribution Notes

Please read and adhere to the developer guide and local patterns and conventions.

[x] The content of this Pull Request (the Contribution) is intentionally submitted for inclusion in CIL (the Work) under the terms and conditions of the Apache-2.0 License
[x] I confirm that the contribution does not violate any intellectual property rights of third parties

Mar 01 '24 16:03 hrobarts

We don't normally require dask do we? Can we avoid it?

On Mon, 4 Mar 2024 at 15:58, Casper da Costa-Luis @.***> wrote:

@.**** commented on this pull request.

In Wrappers/Python/cil/processors/PaganinPhaseProcessor.py https://github.com/TomographicImaging/CIL/pull/1737#discussion_r1511298739 :

@@ -0,0 +1,387 @@ +# -- coding: utf-8 --

⬇️ Suggested change

-# -- coding: utf-8 --

In Wrappers/Python/cil/processors/PaganinPhaseProcessor.py https://github.com/TomographicImaging/CIL/pull/1737#discussion_r1511298821 :

+import dask +from dask import delayed

if dask is a new dependency you'll need to add it to:

recipe/meta.yaml

scripts/create_local_env_for_cil_development.sh

scripts/requirements-test.yml

— Reply to this email directly, view it on GitHub https://github.com/TomographicImaging/CIL/pull/1737#pullrequestreview-1914536482, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACMDSCBUL3QNPDBDJ625ERTYWSD2BAVCNFSM6AAAAABECCP67KVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMYTSMJUGUZTMNBYGI . You are receiving this because you are subscribed to this thread.Message ID: @.***>

Mar 04 '24 15:03 jakobsj

Hi Jakob, yes we can avoid it. I thought it already was a dependency but I can remove it if not

Mar 04 '24 15:03 hrobarts

Thanks for the discussion @jakobsj , @paskino , @gfardell , @casperdcl We discussed:

It makes more sense to keep the Paganin-like filter functionality in the same class as the full Paganin method
We considered having a PhaseRetriever type class which could contain other phase retrieval methods in the future. But for now it makes sense to have a single Paganin class which could in future be wrapped into a bigger phase module.

Here is a possible structure:

class PaganinProcessor(Processor):

    def __init__(self, delta=1, beta = 1e-2, energy = 40000, propagation_distance=None, magnification=None, pixel_size=None, geometry_unit_multiplier = 1,  filter_type='paganin_method'):

    def filter(self):
          data = self.get_input()
          self.create_filter(filter_shape[0], filter_shape[1])

          for i in len(data.geometry.angles):
              projection = data.take(indices = i, axis = out.get_dimension_axis('angle'))
              iffI = ifft2(fft2(projection)*self.filter)
              out.fill(np.squeeze(iffI), angle = i)

         return out

    def create_filter(self):
        ...
    def process(self):
        ...

Suggested usage for the full phase retrieval would be:

processor = PaganinProcessor(delta = 1, beta = 1e-2, energy = 30000, propagation_distance = 0.5)
processor.set_input(data)
thickness = processor.get_output()

And for calling the filter method:

processor = PaganinProcessor()# using default parameters
processor.set_input(data) 
filtered_image = processor.filter()

This uses separate parameters for delta and beta but we could allow the default to be over-ridden with a single parameter when we call filter.

def filter(self, delta_beta):
    '''
    Parameters
    ----------
    delta_beta: float (optional)
        Filter strength. From the ratio of the real and complex part of the material refractive index delta/beta, where refractive index n = (1 - delta) + i beta. Energy-dependent refractive index information for x-ray wavelengths can be found https://henke.lbl.gov/optical_constants/getdb2.html , default is 1e2
    '''

May 07 '24 13:05 hrobarts

After chatting with Gemma, latest suggestion is to add a new argument to get_output so the user can specify if they want the full retrieval or just the filter

def get_output(self, out=None, full_retrieval=True):
            self.full_retrieval = full_retrieval
            return super().get_output(out)

The usage would look like this

processor = PaganinProcessor(delta = 1, beta = 1e-2, energy = 30000)
processor.set_input(data)
thickness = processor.get_output()

And to filter

processor = PaganinProcessor() # using default parameters
processor.set_input(data)
filtered_image = processor.get_output(full_retrieval=False)

We also discussed that we should require the propagation distance, pixel size and magnification to be correctly set. In the geometry, so these should no longer be arguments to the processor. We could allow the values to be over-ridden with set methods

processor = PaganinProcessor() # using default parameters
processor.set_input(data)
processor.update_propagation_distance(0.1)
filtered_image = processor.get_output(full_retrieval=False)

May 09 '24 11:05 hrobarts

Updates from discussion with @gfardell

Could remove the geometry_unit_multiplier argument and require the user to specify the units in the geometry, as we already have a parameter for this in AcquisitionGeometry - create a new PR to make changes to framework.
We discussed that check_input() should raise an error rather than a warning if the distance values in geometry are not valid. This would mean users should call the update_parameters() function before set/check_input() e.g.

processor = PaganinProcessor()
processor.update_parameters({'propagation_distance':new_distance_value})
processor.set_input(data)
thickness = processor.get_output()

But it's confusing to have different behaviour depending on the order the functions are called, and this doesn't follow the behaviour of other processors An alternative could be to have an override_parameters argument to get_output() then check the validity of the distance parameters when process is called. This could look like:

processor = PaganinProcessor()
processor.set_input(data)
processor.get_output(override_parameters = {'pixel_size':1e-6, 'magnification':1})

I think this would mean always using the value from geometry unless an over-ride is passed whereas it would not always be clear if that's the case using update_parameters(). This could also be used to re-configure the processor if you want to check the result for different values of a parameter e.g.

processor = PaganinProcessor(delta=1)
processor.set_input(data)
thickness1 = processor.get_output()
thickness2 = processor.get_output(override_parameters = {'delta':10})
thickness3 = processor.get_output(override_parameters = {'delta':100})

Should update __call__() so arguments to get_output() can also be called in-line

May 23 '24 13:05 hrobarts