linarphy/naroo_reader
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Add intensity step correction

Browse source 
- Finished curve correction
- Add intensity correction for each x
This commit is contained in:
linarphy 2023-05-15 16:43:59 +02:00
parent 2840f86646
commit 15325e0e8f
No known key found for this signature in database
GPG key ID: 3D4AAAC3AD16E79C
1 changed files with 281 additions and 233 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

230
ETA.py
View file

@ -1,6 +1,9 @@
import numpy as np import numpy as np
import matplotlib.pyplot as plt
import utils import utils
import sys import sys
import pathlib
import shelve
from scipy.signal import convolve as sp_convolve from scipy.signal import convolve as sp_convolve
from scipy.signal import find_peaks from scipy.signal import find_peaks
from scipy.ndimage import rotate from scipy.ndimage import rotate
@ -9,56 +12,64 @@ if len( sys.argv ) < 2:
raise Exception( 'this command must have a filename of an ETA fits as an argument' ) raise Exception( 'this command must have a filename of an ETA fits as an argument' )
data = utils.load( sys.argv[1] ) data = utils.load( sys.argv[1] )
""" cache_file = pathlib.Path( 'asset/points_' + sys.argv[1].split( '/' )[-1][:-5] + '.pag' )
find fill points
"""
points = []
points += utils.find_point( data[ : , 0 ] , 0 ) # x_min if cache_file.is_file():
points += utils.find_point( with shelve.open( str( cache_file ) ) as cache:
data = cache[ 'rotated_data' ]
border = cache[ 'border' ]
calibrations = cache[ 'calibrations' ]
else:
"""
find fill points
"""
points = []
points += utils.find_point( data[ : , 0 ] , 0 ) # x_min
points += utils.find_point(
data[ : , data.shape[1] - 1 ], data[ : , data.shape[1] - 1 ],
data.shape[1] - 1 , data.shape[1] - 1 ,
) # x_max ) # x_max
index_min = 0 index_min = 0
while data.shape[0] - 1 > index_min: while data.shape[0] - 1 > index_min:
index_min += 1 index_min += 1
if len( utils.find_point( data[ index_min , : ] , index_min , 'y' ) ) == 3: if len( utils.find_point( data[ index_min , : ] , index_min , 'y' ) ) == 3:
break break
points.append( points.append(
utils.find_point( data[ index_min , : ] , index_min , 'y' )[1] utils.find_point( data[ index_min , : ] , index_min , 'y' )[1]
) # y_min ) # y_min
index_max = data.shape[0] - 1 index_max = data.shape[0] - 1
while index_min < index_max: while index_min < index_max:
index_max -= 1 index_max -= 1
if len( utils.find_point( data[ index_max , : ] , index_max , 'y' ) ) == 3: if len( utils.find_point( data[ index_max , : ] , index_max , 'y' ) ) == 3:
break break
points.append( points.append(
utils.find_point( data[ index_max , : ] , index_max , 'y' )[1] utils.find_point( data[ index_max , : ] , index_max , 'y' )[1]
) # y_max ) # y_max
small_data = utils.compress( data , 5 ) small_data = utils.compress( data , 5 )
points = utils.big_to_small( points , 5 ) points = utils.big_to_small( points , 5 )
# size - 1 # size - 1
points[ 1 ][ 1 ] -= 1 points[ 1 ][ 1 ] -= 1
points[ 3 ][ 0 ] -= 1 points[ 3 ][ 0 ] -= 1
# little shift to be inside the light # little shift to be inside the light
points[ 2 ][ 1 ] += 3 points[ 2 ][ 1 ] += 3
points[ 3 ][ 1 ] += 3 points[ 3 ][ 1 ] += 3
""" """
fill data fill data
""" """
extremum = [] extremum = []
for point in points: for point in points:
if point[0] < points[2][0]: if point[0] < points[2][0]:
point[0] = points[2][0] point[0] = points[2][0]
if point[1] < points[0][1]: if point[1] < points[0][1]:
@ -87,7 +98,7 @@ for point in points:
np.max( taken_points[ : , 0 ] ), np.max( taken_points[ : , 0 ] ),
] ) ] )
border = { border = {
'x': { 'x': {
'min': points[0][1] + extremum[0][1] + 1, 'min': points[0][1] + extremum[0][1] + 1,
'max': points[0][1] + extremum[1][0] , 'max': points[0][1] + extremum[1][0] ,
@ -96,13 +107,13 @@ border = {
'min': points[2][0] + extremum[2][3] + 1, 'min': points[2][0] + extremum[2][3] + 1,
'max': points[2][0] + extremum[3][2] , 'max': points[2][0] + extremum[3][2] ,
}, },
} }
""" """
label deletion label deletion
""" """
mean_data = np.convolve( mean_data = np.convolve(
np.gradient( np.gradient(
np.mean( np.mean(
data[ data[
@ -118,46 +129,46 @@ mean_data = np.convolve(
) ) ) )
), ),
'same' 'same'
) )
mean_data -= np.min( mean_data ) mean_data -= np.min( mean_data )
mean_data /= np.max( mean_data ) mean_data /= np.max( mean_data )
top = utils.consecutive( np.where( mean_data > 0.75 )[0] ) top = utils.consecutive( np.where( mean_data > 0.75 )[0] )
down = utils.consecutive( np.where( mean_data < 0.25 )[0] ) down = utils.consecutive( np.where( mean_data < 0.25 )[0] )
size_top = [ len( list_ ) for list_ in top ] size_top = [ len( list_ ) for list_ in top ]
size_down = [ len( list_ ) for list_ in down ] size_down = [ len( list_ ) for list_ in down ]
label_x = { label_x = {
'min': border[ 'x' ][ 'min' ] + top[ np.argmax( size_top ) ][0] , 'min': border[ 'x' ][ 'min' ] + top[ np.argmax( size_top ) ][0] ,
'max': border[ 'x' ][ 'min' ] + down[ np.argmax( size_down ) ][-1] 'max': border[ 'x' ][ 'min' ] + down[ np.argmax( size_down ) ][-1]
} }
if label_x[ 'min' ] < data.shape[1] // 2: if label_x[ 'min' ] < data.shape[1] // 2:
if label_x[ 'max' ] < data.shape[1] // 2: if label_x[ 'max' ] < data.shape[1] // 2:
border[ 'x' ][ 'min' ] = label_x[ 'max' ] border[ 'x' ][ 'min' ] = label_x[ 'max' ]
else: else:
raise Exception( 'the label seems to be in the middle of the picture' ) raise Exception( 'the label seems to be in the middle of the picture' )
elif label_x[ 'max' ] > data.shape[1] // 2: elif label_x[ 'max' ] > data.shape[1] // 2:
border[ 'x' ][ 'max' ] = label_x[ 'min' ] border[ 'x' ][ 'max' ] = label_x[ 'min' ]
else: else:
raise Exception( 'for an unkown reason, label_x[ \'min\' ] > label_x[ \'max\' ]' ) raise Exception( 'for an unkown reason, label_x[ \'min\' ] > label_x[ \'max\' ]' )
""" """
Rotation Rotation
""" """
index = border[ 'x' ][ 'min' ] index = border[ 'x' ][ 'min' ]
gradient = np.gradient( gradient = np.gradient(
data[ data[
border[ 'y' ][ 'min' ] : border[ 'y' ][ 'min' ] + ( border[ 'y' ][ 'min' ] : border[ 'y' ][ 'min' ] + (
border[ 'y' ][ 'max' ] - border[ 'y' ][ 'min' ] border[ 'y' ][ 'max' ] - border[ 'y' ][ 'min' ]
) // 2, ) // 2,
index index
] ]
) )
while np.max( gradient ) - np.min( gradient ) > 5500: while np.max( gradient ) - np.min( gradient ) > 5500:
index += 1 index += 1
gradient = np.gradient( gradient = np.gradient(
data[ data[
@ -168,7 +179,7 @@ while np.max( gradient ) - np.min( gradient ) > 5500:
] ]
) )
positions = np.argmax( positions = np.argmax(
sp_convolve( sp_convolve(
np.gradient( np.gradient(
data[ data[
@ -183,31 +194,31 @@ positions = np.argmax(
'valid' 'valid'
) , ) ,
axis = 0 axis = 0
) )
list_ = np.arange( 0 , index - border[ 'x' ][ 'min' ] , 1 ) list_ = np.arange( 0 , index - border[ 'x' ][ 'min' ] , 1 )
polyval = np.polyfit( list_ , positions , 1 ) polyval = np.polyfit( list_ , positions , 1 )
angle = np.arctan( polyval[0] ) angle = np.arctan( polyval[0] )
data = rotate( data , angle * ( 180 / np.pi ) ) # utils.rotate does not keep intensity absolute value ? TODO 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' ] ) ) diff_y = int( np.tan( angle ) * ( border[ 'x' ][ 'max' ] - border[ 'x' ][ 'min' ] ) )
border[ 'y' ][ 'min' ] -= diff_y border[ 'y' ][ 'min' ] -= diff_y
border[ 'y' ][ 'max' ] -= diff_y border[ 'y' ][ 'max' ] -= diff_y
""" """
Calibration Calibration
""" """
tot_avg = np.mean( tot_avg = np.mean(
data[ data[
border[ 'y' ][ 'min' ] : border[ 'y' ][ 'max' ], border[ 'y' ][ 'min' ] : border[ 'y' ][ 'max' ],
border[ 'x' ][ 'min' ] : border[ 'x' ][ 'max' ] border[ 'x' ][ 'min' ] : border[ 'x' ][ 'max' ]
] ]
) )
def indicator( list_ ): def indicator( list_ ):
if np.mean( list_ ) > 0.75 * tot_avg: if np.mean( list_ ) > 0.75 * tot_avg:
return 0 return 0
if np.mean( list_ ) < 0.25 * tot_avg: if np.mean( list_ ) < 0.25 * tot_avg:
@ -224,13 +235,13 @@ def indicator( list_ ):
return 4 return 4
return 10 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 ) ] ) indicators = np.array( [ indicator( data[ i , border[ 'x' ][ 'min' ] : border[ 'x' ][ 'max' ] ].copy() ) for i in range( border[ 'y' ][ 'min' ] , border[ 'y' ][ 'max' ] , 1 ) ] )
calibration_areas = utils.consecutive( np.where( indicators == 10 )[0] ) calibration_areas = utils.consecutive( np.where( indicators == 10 )[0] )
calibration_sizes = [ len( calibration_area ) for calibration_area in calibration_areas ] calibration_sizes = [ len( calibration_area ) for calibration_area in calibration_areas ]
y_calibrations = [ calibration_areas[ i ] for i in np.argsort( calibration_sizes ) ][ -2 : ] y_calibrations = [ calibration_areas[ i ] for i in np.argsort( calibration_sizes ) ][ -2 : ]
calibrations = { calibrations = {
'top': { 'top': {
'x': { 'x': {
'min': border['x']['min'], 'min': border['x']['min'],
@ -251,38 +262,75 @@ calibrations = {
'max': border['y']['min'] + y_calibrations[1][-1], 'max': border['y']['min'] + y_calibrations[1][-1],
}, },
}, },
} }
with shelve.open( str( cache_file ) ) as cache:
cache[ 'rotated_data' ] = data
cache[ 'border' ] = border
cache[ 'calibrations'] = calibrations
""" """
stripes curves detection stripes curves detection
""" """
list_ = data[ list_ = data[
calibrations[ 'top' ][ 'y' ][ 'max' ] : calibrations[ 'down' ][ 'y' ][ 'min' ], calibrations[ 'top' ][ 'y' ][ 'max' ] : calibrations[ 'down' ][ 'y' ][ 'min' ],
border[ 'x' ][ 'min' ] : border[ 'x' ][ 'max' ] border[ 'x' ][ 'min' ] : border[ 'x' ][ 'max' ]
].copy() ].copy()
list_ -= np.min( list_ , axis = 0 ) list_ -= np.min( list_ , axis = 0 )
list_ /= np.max( list_ , axis = 0 ) list_ /= np.max( list_ , axis = 0 )
size = list_.shape[1]
size = border[ 'x' ][ 'max' ] - border[ 'x' ][ 'min' ] y_stripe = np.argmax( list_ , axis = 0 )
x_stripe = np.arange( border[ 'x' ][ 'min' ] + 1 * size // 4 , border[ 'x' ][ 'min' ] + 3 * size // 4 , 1 ).astype( int ) good_x = np.where( y_stripe < 2 * np.mean( y_stripe ) )[0]
y_stripe = np.array( [ x_stripe = np.arange( 0 , size , 1 ).astype( int )[ good_x ]
np.where( y_stripe = y_stripe[ good_x ]
list_[ : , x ] > 0.8
)[0][0] for x in x_stripe
] )
stripes = [ # list of polyval result for each stripe stripes = [ # list of polyval result for each stripe
np.polyfit( x_stripe , y_stripe , 2 ) np.polyfit( x_stripe , y_stripe , 3 )
] ]
# First deformation # First deformation
y_diff = np.polyval( stripes[0] , np.arange( 0 , size , 1 ) ).astype( int )
results = np.zeros( ( list_.shape[1] , list_.shape[0] - np.max( y_diff ) ) )
for i in range( list_.shape[1] ):
results[i] = list_[ y_diff[ i ] : list_.shape[0] + y_diff[ i ] - np.max( y_diff ) , i ]
results = results.transpose()
import matplotlib.pyplot as plt y_diff = ( np.polyval( stripes[0] , np.arange( 0 , size , 1 ) ) ).astype( int )
plt.imshow( results ) y_diff[ np.where( y_diff < 0 ) ] = 0
plt.savefig( 'asset/deformation.png' ) results = np.zeros( ( list_.shape[0] + np.max( y_diff ) , list_.shape[1] ) )
for i in range( list_.shape[1] ):
results[ : , i ] = np.concatenate( ( np.zeros( np.max( y_diff ) - y_diff[ i ] ) , list_[ : , i ] , np.zeros( y_diff[i] ) ) )
list_results = np.convolve(
np.gradient(
np.mean( results , axis = 1 ),
) ,
np.ones( 50 ),
'same' ,
)
fall = utils.consecutive( np.where( list_results < - 0.02 )[0] )
fall = np.array( [
np.argmax( list_results )
] + [
consecutive[0] + np.argmin(
list_results[ consecutive[0] : consecutive[-1] ]
) for consecutive in fall
] ).astype( int )
"""
plt.imshow( results , aspect = 'auto' )
plt.hlines( fall , 0 , size )
plt.show()
"""
temp = np.convolve( results[ : , 10000 ] , np.ones( 50 ) , 'same' )
for i in range( len( fall ) - 1 ):
temp[ fall[ i ] : fall[ i + 1 ] ] = np.mean( temp[ fall[ i ] : fall[ i + 1 ] ] )
plt.plot( temp )
plt.plot(
np.convolve(
results[ : , 10000 ],
np.ones( 50 ) ,
'same' ,
),
)
plt.vlines( fall , 0 , 50 , colors = 'red' )
plt.show()