Skip to content

Small update of standardize functions. #46

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.

Sign up for GitHub

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

Open
maelclemence wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Small update of standardize functions. #46

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
41e6f68
Small modifications of phase ramp correciton and standarize set
maelclemence Jul 7, 2025
b472163
adding a different window for FRC than for phase ramp
maelclemence Jul 8, 2025
b6c4dba
modification remove amplitude exponent, there was dimension mismatch
maelclemence Jul 8, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
29 changes: 16 additions & 13 deletions src/cdtools/tools/analysis/analysis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1067,7 +1067,7 @@ def remove_phase_ramp(im, window, probe=None):
Is, Js = np.mgrid[:window.shape[0],:window.shape[1]]
def zero_freq_component(freq):
phase_ramp = np.exp(2j * np.pi * (freq[0] * Is + freq[1] * Js))
return -np.abs(np.sum(phase_ramp * window))**2
return 1/np.abs(np.sum(phase_ramp * window))**2 ##Somehow scipy did not like to minize a negative value, so to circumvent that, we can minimize 1/f instead.

x0 = np.array([0,0])
result = opt.minimize(zero_freq_component, x0)
Expand Down Expand Up @@ -1114,7 +1114,8 @@ def rms_error(x):
if weights is not None:
pix_translations = cdtools.tools.interactions.translations_to_pixel(t.as_tensor(basis), t.as_tensor(translations)).numpy()
pix_translations -= np.min(pix_translations,axis=0)
weights = weights * np.exp(growth_rate[0] * pix_translations[:,0] + growth_rate[1] * pix_translations[:,1])
scale = np.exp(growth_rate[0] * pix_translations[:,0] + growth_rate[1] * pix_translations[:,1])
weights = weights * scale[:,np.newaxis,np.newaxis] # taking into account when dm_rank is not zero.
to_return = to_return + (weights,)

if len(to_return) == 1:
Expand Down Expand Up @@ -1153,7 +1154,8 @@ def standardize_reconstruction_set(
correct_phase_offset=True,
correct_phase_ramp=True,
correct_amplitude_exponent=False,
window=np.s_[:,:],
window = np.s_[:,:],
window_frc = np.s_[:,:],
Copy link
Collaborator

@allevitan allevitan Aug 4, 2025

Choose a reason for hiding this comment

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

This will change the default behavior, i.e. if there was a script which set "window" before, they'd be assuming that this also sets what is now window_frc, but the new behavior requires both to be set.

I suggest we set "window_frc=None", and then in the set it equal to window if it is not set explicitly.

Copy link
Collaborator

@allevitan allevitan Aug 4, 2025

Choose a reason for hiding this comment

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

Also, maybe "frc_window" will be clearer than "window_frc", just because it is read out as a more natural phrase, e.g. "let's change the FRC window" vs "let's change the window FRC"?

nbins=50,
frc_limit='side',
):
Expand Down Expand Up @@ -1213,17 +1215,17 @@ def standardize_reconstruction_set(
obj_1, probe_1, weights_1 = remove_amplitude_exponent(
obj_1, window, probe=probe_1,
weights=half_1['weights'],
basis=half_1['basis'],
basis=half_1['obj_basis'],
Copy link
Collaborator

@allevitan allevitan Aug 4, 2025

Choose a reason for hiding this comment

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

Uh oh, so this does stop this code from failing, but it introduces a silent bug. With the latest versions of fancy_ptycho, obj_basis is not necessarily equal to probe_basis. remove_amplitude_exponent assumes they are equal.

So, the behavior before this PR is to simply fail, which is not good, but we should fix remove_amplitude_exponent before merging. Also, I suspect there is an equivalent silent bug in remove_phase_ramp at present, which we also need to fix 😬

translations=half_1['translations'])
obj_2, probe_2, weights_2 = remove_amplitude_exponent(
obj_2, window, probe=probe_2,
weights=half_2['weights'],
basis=half_2['basis'],
basis=half_2['obj_basis'],
translations=half_2['translations'])
obj, probe, weights = remove_amplitude_exponent(
obj, window, probe=probe,
weights=full['weights'],
basis=full['basis'],
basis=full['obj_basis'],
translations=full['translations'])


Expand All @@ -1233,22 +1235,23 @@ def standardize_reconstruction_set(
obj = np.exp(-1j* np.angle(np.sum(obj[window]))) * obj


# Todo update the translations to account for the determined shift
# Todo: update the translations to account for the determined shift
# Using a different window for the FRC than for the other function
shift_1 = ip.find_shift(
t.as_tensor(ip.hann_window(np.abs(obj[window]))),
t.as_tensor(ip.hann_window(np.abs(obj_1[window]))))
t.as_tensor(ip.hann_window(np.abs(obj[window_frc]))),
t.as_tensor(ip.hann_window(np.abs(obj_1[window_frc]))))
obj_1 = ip.sinc_subpixel_shift(
t.as_tensor(obj_1), shift_1).numpy()

shift_2 = ip.find_shift(
t.as_tensor(ip.hann_window(np.abs(obj[window]))),
t.as_tensor(ip.hann_window(np.abs(obj_2[window]))))
t.as_tensor(ip.hann_window(np.abs(obj[window_frc]))),
t.as_tensor(ip.hann_window(np.abs(obj_2[window_frc]))))
obj_2 = ip.sinc_subpixel_shift(
t.as_tensor(obj_2), shift_2).numpy()

freqs, frc, threshold = calc_frc(
ip.hann_window(obj_1[window]),
Copy link
Collaborator

@allevitan allevitan Aug 4, 2025

Choose a reason for hiding this comment

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

This has nothing to do with your code, but at present this apodization is a bit nonstandard - while we're mucking around with this code, I may push a commit to this section to update the apodization strategy to match e.g. ptychoshelves

ip.hann_window(obj_2[window]),
ip.hann_window(obj_1[window_frc]),
ip.hann_window(obj_2[window_frc]),
full['obj_basis'], nbins=nbins, limit=frc_limit)


Expand Down