diff --git a/ETA.py b/ETA.py index 0d03b84..3cca47c 100644 --- a/ETA.py +++ b/ETA.py @@ -1,6 +1,9 @@ import numpy as np import utils import sys +from scipy.signal import convolve as sp_convolve +from scipy.signal import find_peaks +from scipy.ndimage import rotate if len( sys.argv ) < 2: raise Exception( 'this command must have a filename of an ETA fits as an argument' ) @@ -141,9 +144,97 @@ elif label_x[ 'max' ] > data.shape[1] // 2: else: raise Exception( 'for an unkown reason, label_x[ \'min\' ] > label_x[ \'max\' ]' ) +""" +Rotation +""" + +index = border[ 'x' ][ 'min' ] +gradient = np.gradient( + data[ + border[ 'y' ][ 'min' ] : border[ 'y' ][ 'min' ] + ( + border[ 'y' ][ 'max' ] - border[ 'y' ][ 'min' ] + ) // 2, + index + ] +) +while np.max( gradient ) - np.min( gradient ) > 5500: + index += 1 + gradient = np.gradient( + data[ + border[ 'y' ][ 'min' ] : border[ 'y' ][ 'min' ] + ( + border[ 'y' ][ 'max' ] - border[ 'y' ][ 'min' ] + ) // 2, + index + ] + ) + +positions = np.argmax( + sp_convolve( + np.gradient( + data[ + border[ 'y' ][ 'min' ] : border[ 'y' ][ 'min' ] + ( + border[ 'y' ][ 'max' ] - border[ 'y' ][ 'min' ] + ) // 2 , + border[ 'x' ][ 'min' ] : index + ] , + axis = 0 + ) , + np.ones( ( 100 , 1 ) ), + 'valid' + ) , + axis = 0 +) + +list_ = np.arange( 0 , index - border[ 'x' ][ 'min' ] , 1 ) +polyval = np.polyfit( list_ , positions , 1 ) + +angle = np.arctan( polyval[0] ) + +data = rotate( data , angle * ( 180 / np.pi ) ) # utils.rotate does not keep intensity absolute value ? TODO + +diff_y = int( np.tan( angle ) * ( border[ 'x' ][ 'max' ] - border[ 'x' ][ 'min' ] ) ) + +border[ 'y' ][ 'min' ] -= diff_y +border[ 'y' ][ 'max' ] -= diff_y + +""" +Calibration +""" + +tot_avg = np.mean( + data[ + border[ 'y' ][ 'min' ] : border[ 'y' ][ 'max' ], + border[ 'x' ][ 'min' ] : border[ 'x' ][ 'max' ] + ] +) +def indicator( list_ ): + if np.mean( list_ ) > 0.75 * tot_avg: + return 0 + if np.mean( list_ ) < 0.25 * tot_avg: + return 1 + list_ -= np.min( list_ ) + list_ /= np.max( list_ ) + positions = np.where( list_ > 0.5 )[0] + if len( positions ) < 10: + return 2 + if len( positions ) > 400: + return 3 + distance = np.mean( positions[ 1 : ] - positions[ : -1 ] ) + if distance < 10: + return 4 + return 10 + +indicators = np.array( [ indicator( data[ i , border[ 'x' ][ 'min' ] : border[ 'x' ][ 'max' ] ] ) for i in range( border[ 'y' ][ 'min' ] , border[ 'y' ][ 'max' ] , 1 ) ] ) + +calibration_areas = utils.consecutive( np.where( indicators == 10 )[0] ) + import matplotlib.pyplot as plt +plt.plot( indicators ) +plt.savefig( 'asset/indicator.png' ) +""" plt.imshow( data[ border[ 'y' ][ 'min' ] : border[ 'y' ][ 'max' ], border[ 'x' ][ 'min' ] : border[ 'x' ][ 'max' ] ] ) -plt.savefig( 'asset/test.png' ) +plt.savefig( 'asset/test_rotated.png' ) +"""