diff --git a/src/cdtools/tools/analysis/analysis.py b/src/cdtools/tools/analysis/analysis.py index 9844edb4..25dcf630 100644 --- a/src/cdtools/tools/analysis/analysis.py +++ b/src/cdtools/tools/analysis/analysis.py @@ -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) @@ -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: @@ -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_[:,:], nbins=50, frc_limit='side', ): @@ -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'], 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']) @@ -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]), - 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)