From 920374b14480e6fde0c23f99825318dfd5b2bb8a Mon Sep 17 00:00:00 2001
From: linarphy <linarphy@linarphy.net>
Date: Thu, 20 Jul 2023 17:10:16 +0200
Subject: [PATCH] Update polyfit to interpolation

---
 ETA.py | 122 +++++++++++++++++++++++++++++++++++++++++----------------
 1 file changed, 89 insertions(+), 33 deletions(-)

diff --git a/ETA.py b/ETA.py
index 231e250..30c1704 100755
--- a/ETA.py
+++ b/ETA.py
@@ -8,6 +8,7 @@ import shelve
 from scipy.signal import convolve as sp_convolve
 from scipy.signal import find_peaks
 from scipy.ndimage import rotate
+import scipy.interpolate as interpolate
 
 cache , filename , output , calibration = None , None , None , None
 verbose , no_cache = False , False
@@ -466,7 +467,7 @@ else:
         if verbose:
             print( 'cache saved' )
 
-size_x = border[ 'x' ][ 'max' ]                 - border[ 'x' ][ 'min' ]
+size_x = border[ 'x' ][ 'max' ] - border[ 'x' ][ 'min' ]
 size_y = calibrations[ 'down' ][ 'y' ][ 'min' ] - calibrations[ 'top' ][ 'y' ][ 'max' ]
 
 # Calibration
@@ -478,12 +479,16 @@ if calibration != None:
     peaks_calib = np.loadtxt( calibration ) # sorted list
     peaks_calib = np.sort( peaks_calib )
     mean_up   = np.mean( data[
-        calibrations[ 'top' ][ 'y' ][ 'min' ] : calibrations[ 'top' ][ 'y' ][ 'max' ],
-        calibrations[ 'top' ][ 'x' ][ 'min' ] : calibrations[ 'top' ][ 'x' ][ 'max' ]
+        calibrations[ 'top' ][ 'y' ][ 'min' ]:
+        calibrations[ 'top' ][ 'y' ][ 'max' ],
+        calibrations[ 'top' ][ 'x' ][ 'min' ]:
+        calibrations[ 'top' ][ 'x' ][ 'max' ]
     ] , axis = 0 )
     mean_down = np.mean( data[
-        calibrations[ 'down' ][ 'y' ][ 'min' ] : calibrations[ 'down' ][ 'y' ][ 'max' ],
-        calibrations[ 'down' ][ 'x' ][ 'min' ] : calibrations[ 'down' ][ 'x' ][ 'max' ]
+        calibrations[ 'down' ][ 'y' ][ 'min' ]:
+        calibrations[ 'down' ][ 'y' ][ 'max' ],
+        calibrations[ 'down' ][ 'x' ][ 'min' ]:
+        calibrations[ 'down' ][ 'x' ][ 'max' ]
     ] , axis = 0 )
     peaks_up   = np.array(
         utils.retrieve_peaks(
@@ -559,7 +564,9 @@ if calibration != None:
         cons_diff , i = utils.same_value( diff ) , 0
         for list_same in cons_diff:
             new_data[ i : i + len( list_same ) , x ] = data[
-                calibrations[ 'top' ][ 'y' ][ 'max' ] + i : calibrations[ 'top' ][ 'y' ][ 'max' ] + i + len( list_same ),
+                calibrations[ 'top' ][ 'y' ][ 'max' ] + i :
+                calibrations[ 'top' ][ 'y' ][ 'max' ] + i
+                + len( list_same ),
                 calibrations[ 'top' ][ 'x' ][ 'min' ] + list_same[0]
             ]
             i += len( list_same )
@@ -572,7 +579,7 @@ if calibration != None:
         print( 'y to x pixel calibrated' )
 
     polyval_wavelength = np.polyfit( # x = 0 begins at the start of
-        peaks_up   ,                 # calibrations[ 'top' ][ 'x' ][ 'min' ]
+        peaks_up   ,        # calibrations[ 'top' ][ 'x' ][ 'min' ]
         peaks_calib,
         1          ,
     )
@@ -600,38 +607,88 @@ if verbose:
 """
 remove artefact
 """
-list_ = np.convolve(
-    np.gradient(
-        np.mean(
-            data[
-                calibrations[ 'top' ][ 'y' ][ 'max' ] : calibrations[ 'down' ][ 'y' ][ 'min' ],
-                border[ 'x' ][ 'min' ] : border[ 'x' ][ 'max' ]
-            ]       ,
-            axis = 1,
-        ),
-    )            ,
-    np.ones( 50 ),
-    'same'       ,
+first_stripe = np.argmax(
+    np.convolve(
+        np.gradient(
+            np.mean(
+                data[
+                    calibrations[ 'top' ][ 'y' ][ 'max' ] :
+                    calibrations[ 'down' ][ 'y' ][ 'min' ],
+                    border[ 'x' ][ 'min' ] : border[ 'x' ][ 'max' ]
+                ]       ,
+                axis = 1,
+            ),
+        )            ,
+        np.ones( 50 ),
+        'same'       ,
+    ),
 )
-first_stripe = np.argmax( list_ )
-range_       = range(
+range_       = np.arange(
     calibrations[ 'top'  ][ 'y' ][ 'max' ] + first_stripe,
-    calibrations[ 'down' ][ 'y' ][ 'max' ]
+    calibrations[ 'down' ][ 'y' ][ 'max' ]               ,
+    100                                                  ,
 )
 
 for line in range_:
-    abciss = range(
+    convolve_point = 500
+    abciss = np.arange(
         border[ 'x' ][ 'min' ],
         border[ 'x' ][ 'max' ],
     )
-    list_  = data[
-        line                                           ,
-        border[ 'x' ][ 'min' ] : border[ 'x' ][ 'max' ]
+    list_  = np.convolve(
+        data[
+            line                                           ,
+            border[ 'x' ][ 'min' ] : border[ 'x' ][ 'max' ]
+        ]      ,
+        np.ones(
+            convolve_point,
+        )      ,
+        'valid', # Need to append value at extremities
+    ) / convolve_point # keep correct value
+    list_ = np.concatenate(
+        (
+            data[
+                line                                        ,
+                border[ 'x' ][ 'min' ]                      :
+                border[ 'x' ][ 'min' ] + convolve_point //
+                2 - 1
+            ]    ,
+            list_,
+            data[
+                line                                                ,
+                border[ 'x' ][ 'max' ]
+                - convolve_point // 2 + convolve_point % 2          :
+                border[ 'x' ][ 'max' ]
+            ]    ,
+        ),
+    )
+    number_point  = 10
+    points_taken  = abciss[
+        convolve_point // 2 + convolve_point % 2  :
+        - convolve_point // 2 + convolve_point % 2:
+        ( abciss.shape[0] - convolve_point + 1 ) // ( number_point - 1 )
     ]
-    polyfit_ = np.polyfit(
-        abciss,
-        list_ ,
-        8     ,
+    values_taken  = list_[
+        convolve_point // 2 + convolve_point % 2  :
+        - convolve_point // 2 + convolve_point % 2:
+        ( abciss.shape[0] - convolve_point + 1 ) // ( number_point - 1 )
+    ]
+    sigma         = np.std( list_ )
+    interpolation = interpolate.PchipInterpolator( # Seems the best
+                    # (Have to check to be sure) -> TODO find reference
+                    # value to determine which interpolation is better
+        points_taken,
+        values_taken,
+    )
+    wrong_points = np.where(
+        np.abs(
+            interpolation(
+                abciss,
+            ) - list_ ,
+        ) > sigma
+    )[0]
+    list_[ wrong_points ] = interpolation( # no border point
+        abciss[ wrong_points ],
     )
 
     import matplotlib.pyplot as plt
@@ -641,9 +698,8 @@ for line in range_:
     )
     plt.plot(
         abciss,
-        np.polyval(
-            polyfit_,
-            abciss  ,
+        interpolation(
+            abciss,
         )     ,
     )
     plt.show()