Skip to content
GitLab
Commit 64af4abb authored by christoph's avatar christoph
Browse files

small bug in roi-finders, bug in E0 in xrs_read

parent 5b2d831f
Hide whitespace changes
Inline Side-by-side
XRStools/roifinder_and_gui.py
View file @ 64af4abb
......@@ -9,7 +9,7 @@ from matplotlib.path import Path
from . import xrs_utilities, xrs_rois, xrs_scans, roiSelectionWidget, math_functions
from matplotlib.widgets import Cursor, Button
from scipy.ndimage import measurements
from scipy import signal
from scipy import signal, stats
import copy
import matplotlib
......@@ -850,7 +850,8 @@ class roi_finder:
data from the sample and which one only has background.
Args:
* hydra_obj (hydra_object): Object from the xrs_read.Hydra class that hold scans to be used for the refinement of the ROIs.
* hydra_obj (hydra_object): Object from the xrs_read.Hydra class that hold scans to be used
for the refinement of the ROIs.
* scan_numbers (int or list): Scan numbers of scans to be used in the refinement.
* n_components (int): Number of components in the decomposition.
* method (string): Keyword describing which decomposition to be used ('pca', 'ica', 'nnma').
......@@ -959,6 +960,8 @@ class roi_finder:
# compact the new red_rois
self.roi_obj.red_rois = xrs_rois.convert_matrix_to_redmatrix(self.roi_obj.roi_matrix, labelformat= 'ROI%02d')
def find_pw_rois(self,roi_obj,pw_data,save_dataset=False):
"""
**find_pw_rois**
......@@ -1383,10 +1386,14 @@ class roi_finder:
# prepare the figure
fig, (ax1, ax2) = plt.subplots( 2, 1 )
plt.subplots_adjust(bottom=0.2)
title_txt = 'Active ROI is No. %d of %d'%(1, len(cw_data)) + '.'
plt.suptitle(title_txt)
# plot the very first column spectrum
ax1.plot(cw_data[roi_keys[0]][:,0], lw=2)
ax1.set_xlabel( 'points along scan' )
ax1.set_ylabel( 'intensity [arb. units]' )
ax1.legend(['Column No. %02d / %02d' %(1, cw_data[roi_keys[0]].shape[1])],
frameon=False, loc=1, fontsize='x-small')
# plot the very first mask image
ax2.plot([0, 0], [-0.5, data_dict_norm['ROI00'].shape[1]+0.5, ], 'k-')
ax2.imshow( np.sum(data_dict_norm['ROI00'], axis=0), interpolation=interpolation, cmap=cmap )
......@@ -1400,14 +1407,15 @@ class roi_finder:
# plot the spectra
ax1.cla()
data = cw_data[roi_keys[roi_ind]]
title_txt = 'Click above to keep/below black line to discard column for ROI %02d'%(roi_ind+1) + '.'
title_txt = 'Active ROI is No. %d of %d.'%(roi_ind+1, len(cw_data))
plt.suptitle(title_txt)
ax1.plot(data[:,column_ind], lw=2)
#offset = np.mean(ax1.get_ylim())
#ax1.plot(np.zeros_like(data[:,column_ind]) + offset ,'k-')
ax1.set_xlabel( 'points along scan' )
ax1.set_ylabel( 'intensity [arb. units]' )
ax1.legend(['Column No. %02d / %02d' %(column_ind+1,data.shape[1])])
ax1.legend(['Column No. %02d / %02d' %(column_ind+1,data.shape[1])],
frameon=False, loc=1, fontsize='x-small')
# plot the mask
ax2.cla()
ax2.plot([column_ind-0.5, column_ind+0.5, column_ind+0.5, column_ind-0.5, column_ind-0.5],
......@@ -1437,7 +1445,7 @@ class roi_finder:
self.num_col_min = 0
def next_roi(self, event):
if self.roi_ind < self.num_ROI_max:
if self.roi_ind < self.num_ROI_max-1:
red_rois_copy[self.roi_key][1] = self.red_mask
self.roi_ind += 1
self.roi_key = roi_keys[self.roi_ind]
......@@ -1464,7 +1472,7 @@ class roi_finder:
pass
def next_column( self, event ):
if self.column_ind < self.num_col_max:
if self.column_ind < self.num_col_max-1:
self.column_ind += 1
self.update_shape()
if self.roi_ind < self.num_ROI_max and self.roi_ind >= self.num_ROI_min:
......@@ -1490,7 +1498,7 @@ class roi_finder:
def take( self, event ):
if self.column_ind > self.num_col_min and self.column_ind < self.num_col_max-1:
if self.column_ind > self.num_col_min and self.column_ind < self.num_col_max:
self.red_mask[:,self.column_ind] = self.roi_ind+1
self.red_rois[self.roi_key][1] = self.red_mask
if self.roi_ind < self.num_ROI_max and self.roi_ind >= self.num_ROI_min:
......@@ -1639,13 +1647,13 @@ class roi_finder:
Define ROIs by matplotlib's zoom/pan functions.
Args:
* input_image (np.array)
* input_image (np.array): 2D numpy array with the image to be displayed.
* logscaling (boolean)
* logscaling (boolean): Display image using logarithmic scale.
* colormap (str)
* colormap (str): matplotlib color map string.
* interpolation (str)
* interpolation (str): matplotlib interpolation for the display.
"""
# unset interactive mode
......@@ -1687,7 +1695,7 @@ class roi_finder:
labelbottom=False )
ax.set_xlim( 0.0, input_image.shape[1] )
ax.set_ylim( input_image.shape[0], 0.0 )
title_txt = 'Active Nr. is %d. of %d'%(1, 0)
title_txt = 'Active ROI Nr. is %d. of %d'%(1, 0)
ax.set_title(title_txt)
# activate the zoom function already
......@@ -1699,7 +1707,7 @@ class roi_finder:
fig.canvas.flush_events()
ax.cla()
title_txt = 'Active Nr. is %d. of %d'%( int(roi_key[-2:])+1, len(contour_lines))
title_txt = 'Active ROI Nr. is %d. of %d'%( int(roi_key[-2:])+1, len(contour_lines))
ax.imshow( input_image, interpolation=interpolation, cmap=cmap)
......@@ -1835,11 +1843,11 @@ class roi_finder:
bprev_roi = Button(axprev_roi, 'prev ROI')
bprev_roi.on_clicked(callback.prev_roi)
# take and next button
axtakeN = plt.axes([0.04, 0.05, 0.15, 0.075])
btakeN = Button(axtakeN, 'take+next')
axtakeN = plt.axes([0.04, 0.05, 0.17, 0.075])
btakeN = Button(axtakeN, 'take+next\n (mouse wheel)')
btakeN.on_clicked(callback.take_and_next)
# take and exit
axtakeE = plt.axes([0.24, 0.05, 0.15, 0.075])
axtakeE = plt.axes([0.24, 0.05, 0.17, 0.075])
btakeE = Button(axtakeE, 'take+exit')
btakeE.on_clicked(callback.take_and_exit)
# discard button
......@@ -1869,7 +1877,7 @@ class roi_finder:
self.roi_obj.red_rois = xrs_rois.convert_matrix_to_redmatrix(self.roi_obj.roi_matrix, labelformat= 'ROI%02d')
#return red_rois_copy
def refine_rois_RW(self, hydra_obj, scan_numbers, interpolation='nearest', cmap='Blues'):
def refine_rois_RW( self, hydra_obj, scan_numbers, interpolation='nearest', cmap='Blues' ):
""" **refine_rois_CW**
Allows for manual refinement of ROIs by plotting the spectra row-wise.
......@@ -1958,10 +1966,14 @@ class roi_finder:
# prepare the figure
fig, (ax1, ax2) = plt.subplots( 2, 1 )
plt.subplots_adjust(bottom=0.2)
title_txt = 'Active ROI is No. %d of %d'%(1, len(rw_data)) + '.'
plt.suptitle(title_txt)
# plot the very first column spectrum
ax1.plot(rw_data[roi_keys[0]][:,0], lw=2)
ax1.set_xlabel( 'points along scan' )
ax1.set_ylabel( 'intensity [arb. units]' )
ax1.legend(['Row No. %02d / %02d' %(1, rw_data[roi_keys[0]].shape[1])],
frameon=False, loc=1, fontsize='x-small')
# plot the very first mask image
ax2.plot([-0.5, data_dict_norm['ROI00'].shape[0]+0.5, data_dict_norm['ROI00'].shape[0]+0.5, -0.5],
[-0.5, -0.5, 0.5, 0.5 ], 'k-')
......@@ -1976,14 +1988,15 @@ class roi_finder:
# plot the spectra
ax1.cla()
data = rw_data[roi_keys[roi_ind]]
title_txt = 'Click above to keep/below black line to discard column for ROI %02d'%(roi_ind+1) + '.'
title_txt = 'Active ROI is No. %d of %d.'%(roi_ind+1, len(rw_data))
plt.suptitle(title_txt)
ax1.plot(data[:,row_ind], lw=2)
#offset = np.mean(ax1.get_ylim())
#ax1.plot(np.zeros_like(data[:,column_ind]) + offset ,'k-')
ax1.set_xlabel( 'points along scan' )
ax1.set_ylabel( 'intensity [arb. units]' )
ax1.legend(['Row No. %02d / %02d' %(row_ind+1,data.shape[0])])
ax1.legend(['Row No. %02d / %02d' %(row_ind+1,data.shape[1])],
frameon=False, loc=1, fontsize='x-small')
# plot the mask
ax2.cla()
ax2.plot([-0.5, data_dict_norm[roi_keys[roi_ind]].shape[0]+0.5, data_dict_norm[roi_keys[roi_ind]].shape[0]+0.5, -0.5, -0.5 ],
......@@ -2013,7 +2026,7 @@ class roi_finder:
self.num_row_min = 0
def next_roi(self, event):
if self.roi_ind < self.num_ROI_max:
if self.roi_ind < self.num_ROI_max-1:
red_rois_copy[self.roi_key][1] = self.red_mask
self.roi_ind += 1
self.roi_key = roi_keys[self.roi_ind]
......@@ -2040,7 +2053,7 @@ class roi_finder:
pass
def next_column( self, event ):
if self.row_ind < self.num_row_max:
if self.row_ind < self.num_row_max-1:
self.row_ind += 1
self.update_shape()
if self.roi_ind < self.num_ROI_max and self.roi_ind >= self.num_ROI_min:
......@@ -2066,7 +2079,7 @@ class roi_finder:
def take( self, event ):
if self.row_ind > self.num_row_min and self.row_ind < self.num_row_max-1:
if self.row_ind > self.num_row_min and self.row_ind < self.num_row_max:
self.red_mask[self.row_ind,:] = self.roi_ind+1
self.red_rois[self.roi_key][1] = self.red_mask
if self.roi_ind < self.num_ROI_max and self.roi_ind >= self.num_ROI_min:
......@@ -2088,7 +2101,7 @@ class roi_finder:
def discard( self, event ):
if self.row_ind >= self.num_row_min and self.column_ind <= self.num_row_max:
if self.row_ind >= self.num_row_min and self.row_ind <= self.num_row_max:
self.red_mask[self.row_ind,:] = 0
self.red_rois[self.roi_key][1] = self.red_mask
if self.roi_ind <= self.num_ROI_max and self.roi_ind >= self.num_ROI_min:
......@@ -2143,8 +2156,213 @@ class roi_finder:
self.roi_obj.red_rois = xrs_rois.convert_matrix_to_redmatrix(self.roi_obj.roi_matrix, labelformat= 'ROI%02d')
#return red_rois_copy
def refine_rois_MF_new( self, hydra_obj, scan_numbers, n_components=2,
method='nnma', cov_thresh=-1, cmap='Blues' ):
""" **refine_rois_MF**
Use decomposition of pixelwise data for each ROI to find which of the pixels holds
data from the sample and which one only has background.
Args:
* hydra_obj (hydra_object): Object from the xrs_read.Hydra class that hold scans to be used
for the refinement of the ROIs.
* scan_numbers (int or list): Scan numbers of scans to be used in the refinement.
* n_components (int): Number of components in the decomposition.
* method (string): Keyword describing which decomposition to be used ('pca', 'ica', 'nnma').
"""
# check if scikit learn is available
try:
from sklearn.decomposition import FastICA, PCA, ProjectedGradientNMF
except ImportError:
from sklearn.decomposition import FastICA, PCA
from sklearn.decomposition import NMF as ProjectedGradientNMF
except:
raise ImportError('Please install the scikit-learn package to use this feature.')
return
# try importing a cursor
try:
from matplotlib.widgets import Cursor
except:
pass
# color map for plotting
color_m = plt.cm.get_cmap( name=cmap )
# check if available method is used
if not method in [ 'pca', 'ica', 'nnma' ]:
print('Please use one of the following methods: ' + str(avail_methods) + '!')
return
# make scan_numbers itarable
if isinstance(scan_numbers,list):
scannums = scan_numbers
elif isinstance(scan_numbers,int):
scannums = [scan_numbers]
# make sure there are rois defined
if not self.roi_obj.red_rois:
print("Please define ROIs first.")
return
# get EDF-files and pw_data
scans = []
for ii in scannums:
scan = xrs_scans.Scan()
scan.load( hydra_obj.path, hydra_obj.SPECfname, hydra_obj.EDFprefix, hydra_obj.EDFname, \
hydra_obj.EDFpostfix, ii, \
direct=True, roi_obj=self.roi_obj, scaling=None, scan_type='generic', \
en_column=None, moni_column=hydra_obj.moni_column, method='pixel', comp_factor=None,\
rot_angles=None, clean_edf_stack=False, cenom_dict=None, storeInsets=False )
scans.append( scan )
data_dict_norm = {}
# sum normalized raw data, if more than one scan
if len( scans ) > 1:
# undo individual normalization
for scan in scans:
for key2 in scan.raw_signals:
scan.raw_signals[key2] *= scan.monitor[:,None,None]
# sum all monitor signals
monitor_sum = np.zeros_like( scans[0].monitor )
for scan in scans:
monitor_sum += scan.monitor
# sum all signals
for key2 in scans[0].raw_signals:
data_dict_norm[key2] = np.zeros_like( scans[0].raw_signals[key2] )
for scan in scans:
data_dict_norm[key2] += scan.raw_signals[key2]
# normalize summed signals
for key2 in data_dict_norm:
data_dict_norm[key2] /= monitor_sum[:,None,None]
elif len( scans ) == 1:
for key in scans[0].raw_signals:
data_dict_norm[key] = scans[0].raw_signals[key]
# fetch all red_roi keys
roi_keys = [key for key in sorted(self.roi_obj.red_rois)]
# deep copy all ROIs, set all to zero
red_rois_copy = {}
for key in sorted( self.roi_obj.red_rois ):
red_rois_copy[key] = copy.deepcopy( self.roi_obj.red_rois[key] )
red_rois_copy[key][1] = np.zeros_like( self.roi_obj.red_rois[key][1] )
# get pixel-wise data in matrix form (n_samples, n_features)
pw_data = {}
for key in data_dict_norm:
s = data_dict_norm[key].shape
pw_data[key] = data_dict_norm[key].reshape( s[0], s[1]*s[2], order='C' )
# plot the components
plt.ioff()
fig = plt.figure()
ax0 = plt.subplot2grid( (16, 16), (0, 0), colspan=16, rowspan=16 )
new_rois = {}
for counter, key in enumerate(sorted(pw_data)): # go through each matrix (one per ROI)
# decompose data, choose method
if method == 'nnma': # non negative matrix factorisation
nnm = ProjectedGradientNMF( n_components=n_components )
N = nnm.fit_transform( pw_data[key] )
elif method == 'pca': # principal component analysis
pca = PCA( n_components=n_components )
N = pca.fit_transform( pw_data[key] )
elif method == 'ica': # independent component analysis
ica = FastICA( n_components=n_components )
N = ica.fit_transform( pw_data[key] )
else:
print('No method: \'' + method + '\' available. Will stop here.' )
return
# let user decide which component belongs to the data:
user_choise = 0
# prepare figure
plt.cla()
title_txt = 'Click component that resembles the sample spectrum for ROI %02d'%(counter+1) + '.'
plt.title(title_txt)
# try:
# cursor1 = Cursor(ax0, useblit=True, color='black', linewidth=1)
# except:
# pass
legendstr = []
for ii in range(n_components):
ax0.plot(N[:,ii], label='Component No. %01d'%ii, color=color_m( (ii+1)/(n_components+1)-0.01 ) )
ax0.legend( frameon=False )
ax0.set_xlabel('points along scan')
ax0.set_ylabel('intensity [arb. units]')
plt.draw()
# which curve was chosen
user_input = np.array( plt.ginput(1,timeout=-1)[0] )
nearest_points = [(np.abs(N[:,ii]-user_input[1])).argmin() for ii in range(n_components)]
user_choice = (np.abs(nearest_points-user_input[0])).argmin()
# find covariance for all pixels with user choice
covariance = np.array([])
for ii in range(len(pw_data[key][0,:])):
covariance = np.append( covariance, np.cov(pw_data[key][:,ii], N[:,user_choice])[0,0] )
# try with qui-squared test
# covariance = np.array([])
# for ii in range(len(pw_data[key][0,:])):
# chi_s = stats.chisquare( pw_data[key][:,ii], N[:,user_choice] )
# print( chi_s )
# covariance = np.append( covariance, chi_s[0] )
# plot covariance, let user choose the the cutoff in y direction
plt.cla()
title_txt = 'Click to define a y-threshold for ROI %02d'%(counter+1) + '.'
plt.title( title_txt )
ax0.plot( covariance, '-o', color=color_m(0.9) )
ax0.set_xlabel( 'pixels in ROI' )
ax0.set_ylabel( 'covariance [arb. units]' )
plt.draw()
try:
cursor1 = Cursor(ax0, useblit=True, color='black', linewidth=1, vertOn=False )
except:
pass
if cov_thresh < 0:
user_cutoff = np.array( plt.ginput(1,timeout=-1)[0] )
elif cov_thresh>0 and isinstance(cov_thresh,int):
if len( covariance ) < cov_thresh:
print('ROI has fewer pixels than cov_thresh, will break here.')
return
else:
user_cutoff = np.array([0.0, np.sort(covariance)[-cov_thresh]])
else:
print('Please provide cov_thresh as positive integer!')
# find the ROI indices above the cutoff, reassign ROI indices
inds = covariance >= user_cutoff[1]
s = data_dict_norm[key].shape
red_rois_copy[key][1][np.reshape(inds, (s[1],s[2]))] = counter+1
# reassign ROI object
self.roi_obj.roi_matrix = xrs_rois.convert_redmatrix_to_matrix( red_rois_copy, np.zeros(self.roi_obj.input_image.shape))
self.roi_obj.indices = xrs_rois.convert_matrix_rois_to_inds(self.roi_obj.roi_matrix)
self.roi_obj.kind = 'refined'
self.roi_obj.x_indices = xrs_rois.convert_inds_to_xinds(self.roi_obj.indices)
self.roi_obj.y_indices = xrs_rois.convert_inds_to_yinds(self.roi_obj.indices)
self.roi_obj.masks = xrs_rois.convert_roi_matrix_to_masks(self.roi_obj.roi_matrix)
self.roi_obj.number_of_rois = int(np.amax(self.roi_obj.roi_matrix))
# compact the new red_rois
# self.roi_obj.red_rois = xrs_rois.convert_matrix_to_redmatrix(self.roi_obj.roi_matrix, labelformat= 'ROI%02d')
self.roi_obj.red_rois = red_rois_copy
# =======
......
XRStools/xrs_calctools.py
View file @ 64af4abb
......@@ -567,6 +567,24 @@ class xyzAtom:
norm = np.trapz(self.spectrum[inds,1], self.spectrum[inds,0])
self.spectrum[:,1] /= norm
def getDist( self, atom ):
return np.linalg.norm(self.coordinates - atom.coordinates)
def getDistPBCarb( self, atom, lattice, lattice_inv ):
return getDistancePBC_arb( self, atom, lattice, lattice_inv )
def getAnglePBCarb( self, atom2, atom3, lattice, lattice_inv, degrees=True):
""" **get_angle**
Return angle between the three given atoms (as seen from atom2).
"""
vec1 = getDistVectorPBC_arb(atom2, self, lattice, lattice_inv)
vec2 = getDistVectorPBC_arb(atom3, self, lattice, lattice_inv)
dotp = np.dot(vec1/np.linalg.norm(vec1), vec2/np.linalg.norm(vec2))
if degrees:
return np.degrees( np.arccos( np.clip( dotp, -1.0, 1.0 ) ) )
else:
return np.arccos( np.clip( dotp, -1.0, 1.0 ) )
class xyzMolecule:
""" **xyzMolecule**
......
XRStools/xrs_read.py
View file @ 64af4abb
......@@ -528,8 +528,8 @@ class Hydra:
# master eloss scale in eV is the one of central pixel in first ROI
self.signals = np.zeros((len(self.energy),len(self.cenom_dict)))
self.errors = np.zeros((len(self.energy),len(self.cenom_dict)))
master_eloss = ( self.energy - np.median([np.median(self.cenom_dict[key]) for key in self.cenom_dict]) )*1.0e3
self.E0 = np.mean(self.cenom_dict[first_key][self.cenom_dict[first_key] > 0.0])
master_eloss = ( self.energy - np.median(self.cenom_dict[first_key][self.cenom_dict[first_key] > 0.0]) )*1.0e3
self.E0 = np.median(self.cenom_dict[first_key][self.cenom_dict[first_key] > 0.0])
for key,ii in zip(sorted(self.cenom_dict), range(len(self.cenom_dict))):
print ('Pixel-by-pixel compensation for ' + key +'.')
signal = np.zeros(len(master_eloss))
......@@ -539,17 +539,18 @@ class Hydra:
x = ( self.energy - self.cenom_dict[key][dim1, dim2] )*1.0e3
# signals
y = self.raw_signals[key][:, dim1, dim2]
#print "Y AMAX", signal.max()
#rbfi = Rbf( x, y, function='linear' )
rbfi = interp1d(x, y,bounds_error=False, fill_value=0.0)
#print "rbf AMAX", Rbf( x, y, function='linear' )
#signal += rbfi( master_eloss )
signal += rbfi( master_eloss )
#print "SIGNAL AMAX", signal.max()
# errors
y = self.raw_errors[key][:, dim1, dim2]
rbfi = Rbf( x, y, function='linear' )
error += rbfi( master_eloss )**2
if np.any(y)>0.0:
#print "Y AMAX", signal.max()
#rbfi = Rbf( x, y, function='linear' )
rbfi = interp1d(x, y,bounds_error=False, fill_value=0.0)
#print "rbf AMAX", Rbf( x, y, function='linear' )
#signal += rbfi( master_eloss )
signal += rbfi( master_eloss )
#print "SIGNAL AMAX", signal.max()
# errors
y = self.raw_errors[key][:, dim1, dim2]
rbfi = Rbf( x, y, function='linear' )
error += rbfi( master_eloss )**2
self.signals[:,ii] = signal
self.errors[:,ii] = np.sqrt(error)
self.eloss = master_eloss
......
XRStools/xrs_rois.py
View file @ 64af4abb
......@@ -318,7 +318,7 @@ class roi_object:
* shiftVal : int
Value by which the ROIs should be shifted.
* direction : string
Description of which direction to shit by.
Description of which direction to shift by (can be 'horiz' or 'vert').
* whichroi : sequence
Sequence (iterable) for which ROIs should be shifted.
"""
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment