From ff6f48de881c8f25b6616dc8215c69e8dcf23d36 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Lo=C3=AFc=20Le=20Guyader?= <loic.le.guyader@xfel.eu>
Date: Thu, 16 May 2024 14:06:58 +0200
Subject: [PATCH] Adds function to get first estimate on polyline flat field

---
 src/toolbox_scs/routines/boz.py | 79 ++++++++++++++++++++++++++-------
 1 file changed, 63 insertions(+), 16 deletions(-)

diff --git a/src/toolbox_scs/routines/boz.py b/src/toolbox_scs/routines/boz.py
index 348d993..732f4b4 100644
--- a/src/toolbox_scs/routines/boz.py
+++ b/src/toolbox_scs/routines/boz.py
@@ -929,6 +929,51 @@ def compute_plane_flat_field_correction(rois, params, plot=False):
 
     return flat_field
 
+def initialize_polyline_ff_correction(avg, rois, params, plot=False):
+    """Initialize the polyline flat field correction.
+
+    Inputs
+    ------
+    avg: 2D array, average module image
+    rois: dictionnary of ROIs.
+    plot: boolean, plot initialized polyline versus data projection
+
+    Returns
+    -------
+    fig: handle to figure or None
+    """
+    refn = avg[rois['n']['yl']:rois['n']['yh'],
+             rois['n']['xl']:rois['n']['xh']]
+    refp = avg[rois['p']['yl']:rois['p']['yh'],
+             rois['p']['xl']:rois['p']['xh']]
+    mid = avg[rois['0']['yl']:rois['0']['yh'],
+             rois['0']['xl']:rois['0']['xh']]
+
+    mref = 0.5*(refn + refp)
+
+    inv_signal = mref/mid # normalization
+    projection = inv_signal[5:-5, :].mean(axis=0) #skip first 5 and last 5 px
+    x = np.arange(0, len(projection))
+
+    z = np.polyfit(x, projection, 6)
+
+    if plot:
+        fig, ax = plt.subplots(1, 1, figsize=(4,3))
+        ax.plot(x, projection, label='projection (n+p)/2x0')
+        p = np.poly1d(z)
+        ax.plot(x, p(x), label='poly')
+        ax.legend()
+        ax.set_xlabel('x (px)')
+        ax.set_ylabel('projection')
+    else:
+        fig = None
+
+    # scaling on polynom coefficients for better fitting
+    params.set_flat_field(z/np.logspace(-(z.shape[0]-1), 0, z.shape[0]))
+    params.ff_type = 'polyline'
+
+    return fig
+
 def compute_polyline_flat_field_correction(rois, params, plot=False):
     """Compute the 1D polyline field correction on beam rois.
 
@@ -1154,7 +1199,8 @@ def inspect_ff_fitting_sk(avg, rois, ff, domain=None, vmin=None, vmax=None):
              rois['p']['xl']:rois['p']['xh']]
     ffmid = ff[rois['0']['yl']:rois['0']['yh'],
                rois['0']['xl']:rois['0']['xh']]
-    norm = (ffn+ffp) 
+    np_norm = 0.5*(ffn+ffp)
+    mid_norm = ffmid
 
     fig, axs = plt.subplots(3, 3, sharex=True, sharey=True,
         figsize=(8, 4))
@@ -1167,32 +1213,33 @@ def inspect_ff_fitting_sk(avg, rois, ff, domain=None, vmin=None, vmax=None):
     fig.colorbar(im, ax=axs[1, 0])
 
     im = axs[2, 0].imshow(mid/mref-1, cmap='RdBu_r', vmin=-1, vmax=1)
-    axs[2, 0].set_title('02/(n+p) - 1')
+    axs[2, 0].set_title('2x0/(n+p) - 1')
     fig.colorbar(im, ax=axs[2, 0])
 
-    im = axs[0, 1].imshow(norm)
-    axs[0, 1].set_title('n+p norm')
+    im = axs[0, 1].imshow(np_norm)
+    axs[0, 1].set_title('norm: (n+p)/2')
     fig.colorbar(im, ax=axs[0, 1])
 
-    im = axs[1, 1].imshow(ffmid)
-    axs[1, 1].set_title('0 norm')
+    im = axs[1, 1].imshow(mid_norm)
+    axs[1, 1].set_title('norm: 0')
     fig.colorbar(im, ax=axs[1, 1])
 
-    im = axs[2, 1].imshow(ffmid/norm-1, cmap='RdBu_r', vmin=-1, vmax=1)
-    axs[2, 1].set_title('0/(n+p)norm-1')
+    im = axs[2, 1].imshow(mid_norm/np_norm-1, cmap='RdBu_r', vmin=-1, vmax=1)
+    axs[2, 1].set_title('norm: 2x0/(n+p) - 1')
     fig.colorbar(im, ax=axs[2, 1])
 
 
-    im = axs[0, 2].imshow(mref/norm)
-    axs[0, 2].set_title('(n+p)/norm')
+    im = axs[0, 2].imshow(mref/np_norm)
+    axs[0, 2].set_title('(n+p)/2 /norm')
     fig.colorbar(im, ax=axs[0, 2])
 
-    im = axs[1, 2].imshow(mid)
-    axs[1, 2].set_title('0')
+    im = axs[1, 2].imshow(mid/mid_norm)
+    axs[1, 2].set_title('0 /norm')
     fig.colorbar(im, ax=axs[1, 2])
 
-    im = axs[2, 2].imshow(mid/(mref/norm)-1, cmap='RdBu_r', vmin=-1, vmax=1)
-    axs[2, 2].set_title('0/(n+p) norm')
+    im = axs[2, 2].imshow((mid/mid_norm)/(mref/np_norm)-1,
+                          cmap='RdBu_r', vmin=-1, vmax=1)
+    axs[2, 2].set_title('2x0/(n+p) - 1 /norm')
     fig.colorbar(im, ax=axs[2, 2])
 
     # fig.suptitle(f'{proposalNB}-run{runNB}-dark{darkrunNB} sat={sat_level}')
@@ -1447,7 +1494,7 @@ def ff_refine_fit(params, crit=ff_refine_crit):
         J1 = crit(x, *temp)
         fit_callback.res.append([J0, Jalpha, J1])
         print(f'{fit_callback.counter-1}: {time_delta} '
-                f'({J0}, {Jalpha}, {J1}), {x}')
+                f'(reg. term: {J0}, {Jalpha}, err. term: {J1}), {x}')
 
         return False
 
@@ -2219,7 +2266,7 @@ def process_module(arr, tid, dark, rois, mask=None, sat_level=511,
         v[n] = r_ff[n].sum(axis=(2, 3))
     
     # np_mean roi where we normalize the sum of flat_field
-    np_mean  = 0.5*(r['n'] + r['p'])/(ff['n'] + ff['p'])
+    np_mean  = (r['n'] + r['p'])/(ff['n'] + ff['p'])
     v['np_mean'] = np_mean.sum(axis=(2,3))
 
     res = xr.Dataset()
-- 
GitLab