Skip to content

fixed bugs in PaganinProcessor #2214 - #2225

Merged
hrobarts merged 12 commits into
TomographicImaging:masterfrom
adamdoc:paganin-fix
Nov 6, 2025
Merged

fixed bugs in PaganinProcessor #2214#2225
hrobarts merged 12 commits into
TomographicImaging:masterfrom
adamdoc:paganin-fix

Conversation

@adamdoc

@adamdoc adamdoc commented Oct 21, 2025

Copy link
Copy Markdown
Contributor

Description

Fixed bugs in Paganin Processor that gave incorrect output for cone-beam phase retrieval

Example Usage

sim_image3.tif

from cil.framework import AcquisitionGeometry, AcquisitionData, DataContainer
from cil.processors import PaganinProcessor
import os
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image

# Load simulated paganin image 
corrected =  Image.open('sim_image3.tif')
corrected = np.array(corrected)

# Constants for the retrieval
Z_so = 15
Z_sd = 45
lx = 812
ly = 812
num_proj = 1
px = 8e-3
E = 13000
M = Z_sd/Z_so
x = np.linspace(0,lx*px,lx)

# Set up geometry
ag = AcquisitionGeometry.create_Cone3D(source_position=[0, -Z_so, 0],
                                        detector_position=[0, Z_sd - Z_so, 0],
                                        units='mm').set_panel(num_pixels=[lx, ly], pixel_size=px)\
    .set_angles(angles=[0])
Acq_data = AcquisitionData(corrected, geometry=ag, deep_copy=True)

# Phase retrieval using "paganin_method"
processor = PaganinProcessor(delta = 1e-6, beta = 1e-9, energy = E)
processor.set_input(Acq_data)
Acq_data_ret = processor.get_output()

phase = Image.fromarray(Acq_data_ret.as_array())
phase_np = np.array(phase)
plt.plot(x,phase_np[:,1])

# Phase retrieval using "generalised_paganin_method"
processor = PaganinProcessor(delta = 1e-6, beta = 1e-9, energy = E, filter_type = 'generalised_paganin_method')
processor.set_input(Acq_data)
Acq_data_ret = processor.get_output()

phase = Image.fromarray(Acq_data_ret.as_array())
phase_np = np.array(phase)
plt.plot(x,phase_np[:,1])

Contribution Notes

I made three changes to Paganin_processor.py

  1. Removed mag2 variable as a scaling factor on the retrieved thickness
  2. Defined a new variable as the demagnified pixel size pixel_size_dmag which defines the fourier mesh and filter, replacing use of self.pixel_size
  3. Added self.magnification to the definition of alpha
  • [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

❤️ Thanks for your contribution!

@lauramurgatroyd

lauramurgatroyd commented Oct 22, 2025

Copy link
Copy Markdown
Member

Hi @adamdoc thanks for opening this PR! 😄

@adamdoc

adamdoc commented Oct 24, 2025

Copy link
Copy Markdown
Contributor Author

Hi @lauramurgatroyd
First time doing a PR, I will try again with the tests on my local env and update the PR once it passes

@adamdoc

adamdoc commented Oct 24, 2025

Copy link
Copy Markdown
Contributor Author

Hi @lauramurgatroyd
I've now updated the tests. There are a few lines that assert alpha and filter are calculated correctly, but since I changed the definition of these in PaganinProcessor I've had to update the tests accordingly. Passed on my local branch
numpy.testing.assert_allclose(processor.filter, filter)

@hrobarts hrobarts linked an issue Oct 29, 2025 that may be closed by this pull request
Incorrect output from Paganin Processor with cone beam #2214
Closed
@hrobarts

hrobarts commented Oct 29, 2025

Copy link
Copy Markdown
Contributor

Hi @adamdoc, thank you for looking at this it looks great to me! My only question is about the change that scales alpha by 1/magnification. As I understand it you are now accounting for the magnification in the Fourier-space pixel size. Do we also need to scale alpha as well? The propagation distance we use is the physical sample-to-detector distance and not scaled by magnification, I think in some implementations propagation distance is also scaled so need to scale alpha as well.

@adamdoc

adamdoc commented Nov 3, 2025

Copy link
Copy Markdown
Contributor Author

Hi @adamdoc, thank you for looking at this it looks great to me! My only question is about the change that scales alpha by 1/magnification. As I understand it you are now accounting for the magnification in the Fourier-space pixel size. Do we also need to scale alpha as well? The propagation distance we use is the physical sample-to-detector distance and not scaled by magnification, I think in some implementations propagation distance is also scaled so need to scale alpha as well.

Hi @hrobarts

Yes it's related to Fresnel scaling theorem, which is always something I've struggled to understand why it exists but I know at least it has to be used in the Paganin equation for cone beam setups

Fresnel scaling theorem says you can model a cone beam system using a parallel beam system with a shorter propagation distance (R2/M). One of the ways you can use this is apply theory designed for parallel beam systems to cone beam systems that have been scaled accordingly. The derivation in Paganin's paper is designed for a parallel beam, and states this can also be applied to a cone beam setup by first applying Fresnel scaling theorem (going from Equation 10 to Equation 12 in his 2002 paper). The PaganinProcessor.py function I believe to now implement Equation 12 in that paper, which you will notice has a R2/M term on the filter

So to convert a cone beam image to one that would be captured by a parallel beam through the scaling theorem you scale:
i) the intensity by a factor of M**2, which is just the inverse square law, but we account for this by normalising with the flat field
ii) the transverse coordinates to that on the sample plane, which we do through defining the Fourier grid with pixel_size/M
iii) then finally the propagation distance has to be reduced to R2/M. It's related to that fact a cone beam has a spherical wavefront which is different than a parallel beam with a plane wavefront, and hence phase terms propagates differently (somehow)

I've validated the retrieval with a code we have for simulating phase propagation. Although this does use Fresnel scaling theorem too, so it makes sense that the retrieval un-does this propagation. But happy to chat more !

@hrobarts
hrobarts self-requested a review November 3, 2025 16:59
hrobarts
hrobarts reviewed Nov 4, 2025
View reviewed changes

@hrobarts hrobarts left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Thank you for the explanation @adamdoc . I'm happy to approve the code as it is, I've pushed some changes to the documentation and made a suggestion. Let me know what you think.

Comment thread Wrappers/Python/cil/processors/PaganinProcessor.py
Comment thread Wrappers/Python/cil/processors/PaganinProcessor.py Outdated
hrobarts
hrobarts approved these changes Nov 6, 2025
View reviewed changes

@hrobarts hrobarts left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Thank you very much for this bug fix @adamdoc :) happy to approve and merge now

@hrobarts
hrobarts merged commit 65b66c5 into TomographicImaging:master Nov 6, 2025
11 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Reviewers

@hrobarts hrobarts hrobarts approved these changes

Assignees

No one assigned

Labels

None yet

Projects

None yet

Milestone

No milestone

Development

Successfully merging this pull request may close these issues.

Incorrect output from Paganin Processor with cone beam

3 participants

@adamdoc @lauramurgatroyd @hrobarts