import numpy as np
from scipy import stats
#from scipy import interpolate
from scipy.interpolate import Rbf
#import cv2
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from ..exp import ebsdconst
#from ..cli import config_xcds
[docs]def pcxyz_to_brkr(pc_xyz, top_clip=0.0):
""" Convert detector pc coordinates into Bruker convention
"""
image_width =ebsdconst.BRKR_WIDTH_MICRONS
image_height=(1.0-top_clip)*ebsdconst.BRKR_HEIGHT_MICRONS
pcx= pc_xyz[:,0]/image_width
pcy=1.0 - pc_xyz[:,1]/image_height
dd = np.abs(pc_xyz[:,2]/image_height)
pc_brkr=np.array([pcx,pcy,dd])
return pc_brkr
[docs]def brkr_to_pcxyz(pcbrkr,top_clip=0.0):
""" convert Bruker pc coordinates PCX,PCY,DD into microns in detector system
"""
image_width =ebsdconst.BRKR_WIDTH_MICRONS
image_height=(1.0-top_clip)*ebsdconst.BRKR_HEIGHT_MICRONS
x_microns= +image_width * pcbrkr[:,0]
y_microns= +image_height * (1.0-pcbrkr[:,1])
z_microns= -image_height * pcbrkr[:,2]
pc_xyz=np.array([x_microns,y_microns,z_microns]).T
return pc_xyz
[docs]def brkr_to_gnom(pcx,pcy,dd,aspect):
"""
convert Bruker PCX,PCY,DD,Aspect to gnomonic
"""
y_gn_max= + ( pcy) /dd
y_gn_min= - (1.0-pcy) /dd
x_gn_max= +((1.0-pcx)*aspect)/dd
x_gn_min= -(( pcx)*aspect)/dd
pc_gnom=np.array([x_gn_min,x_gn_max,y_gn_min,y_gn_max],dtype=np.float32)
return pc_gnom
[docs]def pcxyz_to_gnom(pc_xyz, top_clip=0.0):
""" convert detector pc coordinates to gnomonic projection values
"""
image_width =ebsdconst.BRKR_WIDTH_MICRONS
image_height=(1.0-top_clip)*ebsdconst.BRKR_HEIGHT_MICRONS
# coordinates of image borders relative to PC
x_min = - ( pc_xyz[:,0])
y_min = - (image_height - pc_xyz[:,1])
x_max = (image_width - pc_xyz[:,0])
y_max = ( pc_xyz[:,1])
# gnomonic projection: scale by z
# z is minus in pc_xyz
img_xy =np.array([x_min,x_max,y_min,y_max])
img_z = np.array(- pc_xyz[:,2])
img_gnom = img_xy / img_z
return img_gnom.T[0]
[docs]def calibratePC(ScanPointList,bcf_filename,mapinfo,XTilt=-20.0):
"""
This function calibrates the Projection Center from the current PC fit values
in ScanPointList assuming that the step size and the sample tilt is known.
This makes it possible to extrapolate all experimentally determined PC values
to a common reference point (0,0) assuming a regular x-y grid of steps tilted by
around the detector X-axis.
The extrapolated reference PC values are then averaged to result
in the PC estimation.
Notes:
* the SEM beam X-scan is exactly parallel to the X-Tilt detector axis
* Bruker: XTilt=(SampleTilt-90)-DetectorTilt (deg)
* i.e. SampleTilt=70 usually leads to NEGATIVE XTilts!!!
* the method produces PC values in a tilted rectangular region
* no trapezoidal/projective distortion is considered
"""
print('SCALING TRANSFORMATION PC Calibration...')
print('Assuming pure XTilt, no trapez:', XTilt)
f = open(bcf_filename+'_DSCALIBPC_LOG.TXT', "w")
f.write('DynamicS Pattern Center calibration\n')
StepSize=mapinfo['hstep']
print("scanPointList: ", ScanPointList[0,:])
ImageWidthPX =ScanPointList[0,17]
ImageHeightPX=ScanPointList[0,18]
f.write(str(ImageWidthPX) + ' # PatternWIDTH\n')
f.write(str(ImageHeightPX) + ' # PatternHEIGHT\n')
alpha=-XTilt*np.pi/180.0
f.write(str(XTilt) + ' # total tilt angle between sample surface normal and direction to PC\n')
#relative step size
hStepRel=StepSize/ebsdconst.BRKR_WIDTH_MICRONS
# how many pixels moves the pattern center per scan point
hStepPX= hStepRel*ImageWidthPX
vStepPX= hStepPX*np.cos(alpha) # just tilted, no dependence of hStep on vStep
zStepPX= hStepPX*np.sin(alpha)
f.write('assumed PC beam steps in PX:\n')
f.write(str(-hStepPX)+ ' # horizontal PC step size in pixels\n')
f.write(str(vStepPX)+ ' # vertical PC step size in pixels\n')
f.write(str(zStepPX)+ ' # z PC step size in pixels\n')
#meanBX=np.round(np.mean(ScanPointList[:,15]))
#meanBY=np.round(np.mean(ScanPointList[:,16]))
meanBX=0.0
meanBY=0.0
f.write('Reference beam indices:\n')
f.write(str(meanBX)+ '\n')
f.write(str(meanBY)+ '\n')
dBX=ScanPointList[:,12]-meanBX
dBY=ScanPointList[:,13]-meanBY
# make absolute values in pixels
# get away from 2 different units for PCX,PCY,DD given in Bruker software
pPCX= ScanPointList[:, 9] * ImageWidthPX
pPCY= (1.0-ScanPointList[:,10]) * ImageHeightPX
pDD = ScanPointList[:,11] * ImageHeightPX
# back-interpolate to mean reference PC
pPCX0= np.mean(pPCX + dBX*hStepPX)
pPCY0= np.mean(pPCY - dBY*vStepPX)
pDD0 = np.mean(pDD - dBY*zStepPX)
f.write('reference PC (in pixels, PCX from left, PCY from bottom):\n')
f.write(str(pPCX0)+ '\n')
f.write(str(pPCY0)+ '\n')
f.write(str(pDD0)+ '\n')
f.write('reference PC in PCX,PCY,DD , PCY from top/Height, DD =/Height):\n')
f.write(str(pPCX0/ImageWidthPX)+ ' # measured from LEFT of pattern, in units of PatternWIDTH\n')
f.write(str(1.0-pPCY0/ImageHeightPX)+ ' # measured from TOP of pattern, in units of PatternHEIGHT \n')
f.write(str(pDD0 /ImageHeightPX)+ ' # to sample direction, in units of PatternHEIGHT \n')
f.close()
# explicit maps with the PC values
map_width =mapinfo['nwidth']
map_height=mapinfo['nheight']
# note iz is useless here, dd changes with iy
ix,iy,_ = make_map_indices(map_width,map_height)
pcx_px = pPCX0 + ix * (-hStepPX)
pcy_px = pPCY0 + iy * (+vStepPX)
dd_px = pDD0 + iy * (+zStepPX)
pcx_brkr = pcx_px/ImageWidthPX
pcy_brkr = (1.0-pcy_px/ImageHeightPX)
dd_brkr = dd_px/ImageHeightPX
np.savetxt(bcf_filename+'_PCX.map',pcx_brkr.reshape((map_height,map_width)), fmt='%10.6f')
np.savetxt(bcf_filename+'_PCY.map',pcy_brkr.reshape((map_height,map_width)), fmt='%10.6f')
np.savetxt(bcf_filename+'_DD.map' , dd_brkr.reshape((map_height,map_width)), fmt='%10.6f')
return
[docs]def calibratePC_BRKR(pcdata,mapinfo,XTilt=-20.0):
"""
This function calibrates the Projection Center from the current PC fit values
in ScanPointList assuming that the step size and the sample tilt is known.
This makes it possible to extrapolate all experimentally determined PC values
to a common reference point (0,0) assuming a regular x-y grid of steps tilted by
around the detector X-axis.
The extrapolated reference PC values are then averaged to result
in the PC estimation.
Notes:
* the SEM beam X-scan is exactly parallel to the X-Tilt detector axis
* Bruker: XTilt=(SampleTilt-90)-DetectorTilt (deg)
* i.e. SampleTilt=70 usually leads to NEGATIVE XTilts!!!
* the method produces PC values in a tilted rectangular region
* no trapezoidal/projective distortion is considered
"""
print('SCALING TRANSFORMATION PC Calibration...')
print('Assuming pure XTilt, no trapez:', XTilt)
f = open('XCDSCALIBPC_LOG.TXT', "w")
f.write('SEM SCAN BASED Pattern Center calibration\n')
step_size = 0.094 #mapinfo['hstep']
ImageWidthPX = 160
ImageHeightPX = 120
f.write(str(step_size) + ' # step size (microns)\n')
f.write(str(ImageWidthPX) + ' # PatternWIDTH\n')
f.write(str(ImageHeightPX) + ' # PatternHEIGHT\n')
alpha=-XTilt*np.pi/180.0
f.write(str(XTilt) + ' # total tilt angle between sample surface normal and direction to PC\n')
#relative step size
hStepRel=StepSize/ebsdconst.BRKR_WIDTH_MICRONS
# how many pixels moves the pattern center per scan point
hStepPX= hStepRel*ImageWidthPX
vStepPX= hStepPX*np.cos(alpha) # just tilted, no dependence of hStep on vStep
zStepPX= hStepPX*np.sin(alpha)
f.write('assumed PC beam steps in PX:\n')
f.write(str(-hStepPX)+ ' # horizontal PC step size in pixels\n')
f.write(str(vStepPX)+ ' # vertical PC step size in pixels\n')
f.write(str(zStepPX)+ ' # z PC step size in pixels\n')
#meanBX=np.round(np.mean(ScanPointList[:,15]))
#meanBY=np.round(np.mean(ScanPointList[:,16]))
meanBX=0.0
meanBY=0.0
f.write('Reference beam indices:\n')
f.write(str(meanBX)+ '\n')
f.write(str(meanBY)+ '\n')
dBX=ScanPointList[:,12]-meanBX
dBY=ScanPointList[:,13]-meanBY
# make absolute values in pixels
# get away from 2 different units for PCX,PCY,DD given in Bruker software
pPCX= ScanPointList[:, 9] * ImageWidthPX
pPCY= (1.0-ScanPointList[:,10]) * ImageHeightPX
pDD = ScanPointList[:,11] * ImageHeightPX
# back-interpolate to mean reference PC
pPCX0= np.mean(pPCX + dBX*hStepPX)
pPCY0= np.mean(pPCY - dBY*vStepPX)
pDD0 = np.mean(pDD - dBY*zStepPX)
f.write('reference PC (in pixels, PCX from left, PCY from bottom):\n')
f.write(str(pPCX0)+ '\n')
f.write(str(pPCY0)+ '\n')
f.write(str(pDD0)+ '\n')
f.write('reference PC in PCX,PCY,DD , PCY from top/Height, DD =/Height):\n')
f.write(str(pPCX0/ImageWidthPX)+ ' # measured from LEFT of pattern, in units of PatternWIDTH\n')
f.write(str(1.0-pPCY0/ImageHeightPX)+ ' # measured from TOP of pattern, in units of PatternHEIGHT \n')
f.write(str(pDD0 /ImageHeightPX)+ ' # to sample direction, in units of PatternHEIGHT \n')
f.close()
# explicit maps with the PC values
map_width =mapinfo['nwidth']
map_height=mapinfo['nheight']
# note iz is useless here, dd changes with iy
ix,iy,_ = make_map_indices(map_width,map_height)
pcx_px = pPCX0 + ix * (-hStepPX)
pcy_px = pPCY0 + iy * (+vStepPX)
dd_px = pDD0 + iy * (+zStepPX)
pcx_brkr = pcx_px/ImageWidthPX
pcy_brkr = (1.0-pcy_px/ImageHeightPX)
dd_brkr = dd_px/ImageHeightPX
np.savetxt(bcf_filename+'_PCX.map',pcx_brkr.reshape((map_height,map_width)), fmt='%10.6f')
np.savetxt(bcf_filename+'_PCY.map',pcy_brkr.reshape((map_height,map_width)), fmt='%10.6f')
np.savetxt(bcf_filename+'_DD.map' , dd_brkr.reshape((map_height,map_width)), fmt='%10.6f')
return
[docs]def plotPC(secondary,fit,pc3=None,units='Bruker',plotdir=''):
"""
plots the pattern center data and shows parameters
"""
plt.rc('xtick', labelsize=14)
plt.rc('ytick', labelsize=14)
#plt.title(title,fontsize=26)
plt.ylabel('PCY (pattern height from top)', color='k', fontsize=22)
plt.gca().invert_yaxis() # consistent with y measured from top of pattern
plt.xlabel('PCX (pattern width from left)', color='k',fontsize=22)
plt.scatter(fit.T[0],fit.T[1],s=22,c='red',alpha=1.0, lw=0)
plt.scatter(secondary.T[0],secondary.T[1],s=22,c='blue',alpha=1.0,lw=0)
plt.grid(True)
plt.axes().set_aspect('equal', 'datalim')
plt.savefig(plotdir+'/PCXY_FIT.png',dpi=300,bbox_inches = 'tight')
#plt.show()
plt.close()
plt.rc('xtick', labelsize=14)
plt.rc('ytick', labelsize=14)
#plt.title(title,fontsize=26)
plt.ylabel('DD (pattern height)', color='k', fontsize=22)
plt.xlabel('PCX (pattern width from left)', color='k',fontsize=22)
plt.scatter(fit.T[0],fit.T[2],s=22,c='red',alpha=1.0,lw=0)
plt.scatter(secondary.T[0],secondary.T[2],s=22,c='blue',alpha=1.0,lw=0)
plt.grid(True)
plt.axes().set_aspect('equal', 'datalim')
plt.savefig(plotdir+'/PCXZ_FIT.png',dpi=300,bbox_inches = 'tight')
#plt.show()
plt.close()
# linear fit to YZ-curve to estimate tilt angle
slope, intercept, r_value, p_value, std_err_slope = stats.linregress(fit.T[2],fit.T[1])
print(slope, intercept, r_value, p_value, std_err_slope)
XTilt=-np.abs(np.arctan(1.0/slope)*180.0/np.pi)
print('approximated total XTilt=(SampleTilt-90)-DetectorTilt (deg):',XTilt)
fittedY=intercept+slope*fit.T[2]
plt.rc('xtick', labelsize=14)
plt.rc('ytick', labelsize=14)
#plt.title(title,fontsize=26)
plt.xlabel('DD (pattern height)', color='k', fontsize=22)
plt.ylabel('PCY (pattern height from top)', color='k',fontsize=22)
plt.gca().invert_yaxis() # consistent with y measured from top of pattern
plt.scatter(fit.T[2],fit.T[1],s=22,c='red',alpha=1.0,lw=0)
plt.scatter(secondary.T[2],secondary.T[1],s=22,c='blue',alpha=1.0,lw=0)
plt.plot(fit.T[2],fittedY,'r-' )
plt.grid(True)
plt.axes().set_aspect('equal', 'datalim')
plt.savefig(plotdir+'/PCYZ_FIT.png',dpi=300,bbox_inches = 'tight')
#plt.show()
plt.close()
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
plt.gca().invert_yaxis() # consistent with y measured from top of pattern
ax.scatter(secondary.T[0],secondary.T[2],secondary.T[1],s=16, c='b', marker='o',lw=0,alpha=1.0)
ax.scatter(fit.T[0],fit.T[2],fit.T[1], c='red', s=16, marker='o',lw=0,alpha=1.0)
plt.grid(True)
ax.set_xlabel('PCX')
ax.set_ylabel('DD')
ax.set_zlabel('PCY')
plt.savefig(plotdir+'/PCXYZ_FIT.png',dpi=300,bbox_inches = 'tight')
#plt.show()
plt.close()
return XTilt
[docs]def normrange(Y):
"""normalize the range of the Y values
"""
Ymax=np.amax(Y,axis=0)
Ymin=np.amin(Y,axis=0)
dY =Ymax-Ymin
Yn=(Y-Ymin)/dY - 0.5
return Yn
[docs]def make_map_indices(w,h):
"""
makes 3 1D arrays of x and y, z=1 indices of a map w*h
for projective transformations
"""
ix = np.tile(np.arange(0, w), h)
iy = np.reshape(np.tile(np.arange(0, h), [w,1]).T, -1)
iz = np.ones_like(ix)
return ix,iy,iz
[docs]def project_points(pts_src_3d,h):
""" make 2D plane-projected points from homography matrix
"""
# to get same result with numpy
fit = np.dot(pts_src_3d,h.T)
projected = fit / fit[:,2,np.newaxis] # dehomogenize
projected[:,2]=0.0
#print(fit)
#using OpenCV, assumes 2d points
#pts_src=np.array([pts_src]) # 2d
#projected = cv2.perspectiveTransform(pts_src, h)[0]
return projected
[docs]def make_projective_PCdata(A,MapWidth,MapHeight,outfile=None):
"""
apply projective transformation matrix to all map indices
TODO: 1. save projection center values for direct import into DynamicS
2. save PC values in HDF5
"""
# list of all beam index coordinates
# add z=1 for homogeneous coordinates in 2D
ix,iy,iz=make_map_indices(MapWidth,MapHeight)
map_coo=np.array([ix,iy,iz]).T
# apply projective transformation
#pad = lambda x: np.hstack([x, np.ones((x.shape[0], 1))])
#unpad = lambda x: x[:,:-1]
#transform = lambda x: unpad(np.dot(pad(x), A))
pc_coo=np.dot(map_coo,A)
#pc_coo=transform(map_coo)
if outfile:
np.savetxt(outfile+'_A_proj.dat',A)
np.savetxt(outfile+'_PCX.map',pc_coo.T[0].reshape((MapHeight,MapWidth)))
np.savetxt(outfile+'_PCY.map',pc_coo.T[1].reshape((MapHeight,MapWidth)))
np.savetxt(outfile+'_DD.map' ,pc_coo.T[2].reshape((MapHeight,MapWidth)))
return
[docs]def make_interpolated_PCdata(xc_beam_indices,xc_pc_coords,nwidth,nheight,outfile=''):
"""
interpolates between projection center values, input assumes column vectors
Radial Basis Functions for interpolation of scattered data:
http://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html
"""
print('Interpolating Projection Center Map Data...')
ix=xc_beam_indices[0]
iy=xc_beam_indices[1]
f_pcx = Rbf(ix, iy, xc_pc_coords[0,:])
f_pcy = Rbf(ix, iy, xc_pc_coords[1,:])
f_dd = Rbf(ix, iy, xc_pc_coords[2,:])
xind = np.arange(0,nwidth)
yind = np.arange(0,nheight)
ix_map, iy_map = np.meshgrid(xind, yind)
pcx= f_pcx(ix_map, iy_map)
pcy= f_pcy(ix_map, iy_map)
dd = f_dd(ix_map, iy_map)
if outfile>'':
np.savetxt(outfile+'_PCX_interp.map',pcx, fmt='%10.6f')
np.savetxt(outfile+'_PCY_interp.map',pcy, fmt='%10.6f')
np.savetxt(outfile+'_DD_interp.map' ,dd , fmt='%10.6f')
return
[docs]def PCstats(DataMatrix0,mean_dim,title):
""" Statistics of pattern center positions
PCX, PCY, and DD should be constant along rows or colums respectively
thus the stdDev along a column gives an error estimate
"""
if (mean_dim==1):
DataMatrix=DataMatrix0.T
else:
DataMatrix=DataMatrix0
# mask all zero values (skipped map points)
#DataMatrix=np.ma.masked_values(DataMatrix0,0.0)
print(DataMatrix)
print(title)
print('Mean in best fit:');
# make data sets for row/column means and std
PCMean=np.zeros(DataMatrix.shape[0])
PCStd=np.zeros(DataMatrix.shape[0])
for i in range(DataMatrix.shape[0]):
PCMean[i]=np.mean(DataMatrix[i,1:-1]) # exclude first and last point
PCStd[i]=np.std(DataMatrix[i,1:-1]) # exclude first and last point
if not (np.isnan(PCMean[i]) or np.isnan(PCStd[i])):
print(i+1, PCMean[i], PCStd[i])
print('Mean error bar (all valid points):',np.mean(np.ma.masked_invalid(PCStd)[1:-1]))
PCMean=np.ma.masked_invalid(PCMean)
BeamPos=np.ma.masked_where(np.ma.getmask(PCMean), np.arange(DataMatrix.shape[0]))
#print(BeamPos,PCMeanMasked)
#print(len(BeamPos),len(PCMeanMasked))
mask = np.isfinite(PCMean) & np.isfinite(BeamPos)
slope, intercept, r_value, p_value, std_err_slope = stats.linregress(BeamPos[mask][1:-1],PCMean[mask][1:-1])
print(slope, intercept, r_value, p_value, std_err_slope)
fittedY=intercept+slope*BeamPos
residuals_y=PCMean-fittedY
print('StdDev of Y residuals:',np.std(residuals_y))
#print(fittedY)
plt.rc('xtick', labelsize=18)
plt.rc('ytick', labelsize=18)
plt.title(title,fontsize=26)
plt.ylabel('position (microns)', color='k', fontsize=22)
plt.xlabel('beam scan index', color='k',fontsize=22)
plt.plot(BeamPos[mask],fittedY[mask],'r-' )
plt.errorbar(BeamPos,PCMean,yerr=PCStd,fmt='o')
plt.savefig(title+'.png',dpi=300,bbox_inches = 'tight')
plt.show()
plt.close()
print('')