Add auto calibration based on peakless signal

wavelength_calibration.py

@@ -13,20 +13,109 @@ reference  = np.load( 'asset/reference.npy' )
 mean_calib_top  = np.mean( calib_top  , axis = 0 )
 mean_calib_down = np.mean( calib_down , axis = 0 )
 
-peaks_calib_top = utils.retrieve_peaks( mean_calib_top )
+peaks_calib_top = np.array(
-peaks_reference = utils.retrieve_peaks( reference[1] , window_size = 1 , max_window_size = 1 )
+    utils.retrieve_peaks( mean_calib_top )
-
+).astype( int )
-polyval = np.polyfit( peaks_calib_top , peaks_reference[ 2 : ] , 1 )
+peaks_reference = np.array(
-
+    utils.retrieve_peaks(
-peaks_values = [ reference[ 1 , i ] for i in np.array( peaks_reference ).astype( int ) ]
+        reference[1]       ,
-sorting = np.argsort( peaks_values )[ :: -1 ]
+        window_size = 1    ,
-wavelength = [ reference[ 0 , i ] for i in np.array( [ peaks_reference[ j ] for j in sorting ] ).astype( int ) ]
+        max_window_size = 1,
-
+    )
-polyval_wavelength = np.polyfit( peaks_reference , np.sort( wavelength ) , 1 )
+).astype( int )
-
-wavelength = np.polyval( polyval_wavelength , np.polyval( polyval , np.arange( 0 , len( mean_calib_top ) , 1 ) ) )
 
 """
+Signal without peaks
+"""
+
+peakless_calib_top = mean_calib_top.copy()
+for peak in peaks_calib_top:
+    first = peak
+    peak , old = first - 2 , first - 1
+    while peakless_calib_top[ peak ] <= peakless_calib_top[ old ]:
+        old   = peak
+        peak -= 1
+    peakless_calib_top[ peak : first ] = peakless_calib_top[ peak ] * np.ones( first - peak )
+
+    peak , old = first + 2 , first + 1
+    while peakless_calib_top[ peak ] <= peakless_calib_top[ old ]:
+        old   = peak
+        peak += 1
+    peakless_calib_top[ first : peak ] = peakless_calib_top[ peak ] * np.ones( peak - first )
+
+peakless_calib_top = sig.medfilt( peakless_calib_top , 111 )
+
+"""
+It's possible to have an idea of the part the spectrum begins with the
+small large peak at the start of it, the peak at 3810 A should be inside
+it. The same goes for the area at the end.
+ONLY TRUE FOR THE CURRENT CALIBRATION (Hg)
+"""
+argmin = np.argmin(
+    peakless_calib_top[ : len( peakless_calib_top ) // 2 ]
+)
+argmax = np.argmax( peakless_calib_top[ : argmin ] )
+peaks_inside = np.where(
+    np.logical_and(
+        argmax < peaks_calib_top,
+        argmin > peaks_calib_top,
+    )
+)[0]
+if len( peaks_inside ) < 1:
+    raise Error( 'unknown plage, cannot autocalibrate' )
+first_peak_cal = peaks_inside[0]
+first_peak_ref = 3 # hard-coded
+
+"""
+The next peak after the minimum at the end of the spectrum should be
+5079 A.
+ONLY TRUE FOR THE CURRENT CALIBRATION (Hg)
+"""
+argmin_1 = np.argmin(
+    peakless_calib_top[ len( peakless_calib_top ) // 2 : ]
+) + len( peakless_calib_top ) // 2
+
+peaks_inside = np.where(
+    argmin_1 < peaks_calib_top
+)[0]
+if len( peaks_inside ) < 1:
+    raise Error( 'unknown plage, cannot autocalibrate' )
+last_peak_cal = peaks_inside[0]
+last_peak_ref = 20 # hard-coded
+
+print( first_peak_cal , last_peak_cal )
+polyval = np.polyfit(
+    peaks_calib_top[ first_peak_cal : last_peak_cal ],
+    peaks_reference[ first_peak_ref : last_peak_ref ],
+    1                                                ,
+) # We suppose there is as much calib peak than ref peak for now
+
+peaks_values = [ reference[ 1 , i ] for i in peaks_reference ]
+sorting = np.argsort( peaks_values )[ :: -1 ]
+wavelength = np.array( [
+    reference[ 0 , i ] for i in np.array( [
+        peaks_reference[ j ] for j in sorting
+    ] ).astype( int )
+] )
+
+polyval_wavelength = np.polyfit(
+    peaks_reference,
+    np.sort( wavelength ),
+    1,
+)
+
+wavelength = np.polyval(
+    polyval_wavelength,
+    np.polyval(
+        polyval,
+        np.arange(
+            0                    ,
+            len( mean_calib_top ),
+            1                    ,
+        )
+    )
+)
+
 manager = qt.FigureManager(
     qt.FigureCanvas(
         fig.Figure(
@@ -35,73 +124,22 @@ manager = qt.FigureManager(
     ),
     0,
 )
-
+manager.canvas.figure.gca().plot(
-manager.canvas.figure.add_subplot( 2 , 1 , 1 , xlim = ( 3800 , 4000 ) , xmargin = 0 ).plot(
     wavelength    ,
     mean_calib_top,
 )
-manager.canvas.figure.axes[0].set_title( 'raw data' )
-
-manager.canvas.figure.add_subplot( 2 , 1 , 2 , xlim = ( 3800 , 4000 ) , xmargin = 0 ).plot(
-    reference[0],
-    reference[1],
-)
-manager.canvas.figure.axes[1].set_title( 'reference' )
-
-manager.show()
-manager.start_main_loop()
-"""
-for peak in np.array( peaks_calib_top ).astype( int ):
-    first = peak
-    peak , old = first - 2 , first - 1
-    while mean_calib_top[ peak ] <= mean_calib_top[ old ]:
-        old   = peak
-        peak -= 1
-    mean_calib_top[ peak : first ] = mean_calib_top[ peak ] * np.ones( first - peak )
-
-    peak , old = first + 2 , first + 1
-    while mean_calib_top[ peak ] <= mean_calib_top[ old ]:
-        old   = peak
-        peak += 1
-    mean_calib_top[ first : peak ] = mean_calib_top[ peak ] * np.ones( peak - first )
-
-manager = qt.FigureManager(
-    qt.FigureCanvas(
-        fig.Figure(
-            figsize = ( 10 , 5 )
-        ),
-    ),
-    0,
-)
 manager.canvas.figure.gca().plot(
-    mean_calib_top
+    wavelength  ,
-)
+    peakless_calib_top,
-manager.canvas.figure.gca().plot(
-    sig.medfilt( mean_calib_top , 71 )
-)
-manager.show()
-manager.start_main_loop()
-
-manager = qt.FigureManager(
-    qt.FigureCanvas(
-        fig.Figure(),
-    ),
-    1,
-)
-signal_calib = sig.convolve(
-    sig.medfilt( mean_calib_top , 71 ),
-    np.ones( 50 ),
-    'same'
-)
-manager.canvas.figure.gca().plot(
-    signal
 )
 manager.canvas.figure.gca().vlines(
     [
-        np.argmin( signal[ 25 : -25 ] )
+        wavelength[ argmin ]  ,
+        wavelength[ argmax ]  ,
+        wavelength[ argmin_1 ],
     ]                           ,
-    np.min( signal[ 25 : -25 ] ),
+    np.min( peakless_calib_top ),
-    np.max( signal[ 25 : -25 ] ),
+    np.max( peakless_calib_top ),
     color = 'red'               ,
 )
 manager.show()