# -*- coding: utf-8 -*-
#GSASII image calculations: Image calibration, masking & integration routines.
'''
Classes and routines defined in :mod:`GSASIIimage` follow.
'''
from __future__ import division, print_function
import math
import time
import copy
import numpy as np
import numpy.linalg as nl
import numpy.ma as ma
from scipy.optimize import leastsq
import scipy.interpolate as scint
from . import GSASIIpath
GSASIIpath.SetBinaryPath()
from . import GSASIIlattice as G2lat
from . import GSASIIpwd as G2pwd
from . import GSASIIspc as G2spc
#import GSASIImath as G2mth
from . import GSASIIfiles as G2fil
from . import ImageCalibrants as calFile
# trig functions in degrees
sind = lambda x: math.sin(x*math.pi/180.)
asind = lambda x: 180.*math.asin(x)/math.pi
sinhd = lambda x: math.sinh(x*math.pi/180.)
tand = lambda x: math.tan(x*math.pi/180.)
atand = lambda x: 180.*math.atan(x)/math.pi
atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
cosd = lambda x: math.cos(x*math.pi/180.)
acosd = lambda x: 180.*math.acos(x)/math.pi
coshd = lambda x: math.cosh(x*math.pi/180.)
rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
#numpy versions
npsind = lambda x: np.sin(x*np.pi/180.)
npasind = lambda x: 180.*np.arcsin(x)/np.pi
npcosd = lambda x: np.cos(x*np.pi/180.)
npacosd = lambda x: 180.*np.arccos(x)/np.pi
nptand = lambda x: np.tan(x*np.pi/180.)
npatand = lambda x: 180.*np.arctan(x)/np.pi
npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
nxs = np.newaxis
debug = False
[docs]
def peneCorr(tth,dep,dist):
''' Compute empirical position correction due to detector absorption
:param float/array tth: angle between detector normal & scattered ray vector
:param float dep: coefficient
:param float dist: sample to detector surface in mm
:returns: float/array distance: correction for penetration
'''
return dep*(1.-npcosd(tth))*dist**2/1000. #best one
[docs]
def GetTthP(x,y,parmDict):
''' Compute angle between detector normal & sample scattering ray vector
:param float/array x: detector x-position in mm
:param float/array y: detector y-position in mm
:param dict parmDict: dictionary of detector orientation parameters in fitting routine
names are: det-X, det-Y, tilt, dist & phi
:returns: float/array angle: = 2-theta if tilt is zero
'''
def costth(xyz,d0):
''' compute cos of angle between vectors; xyz not normalized, d0 normalized'''
u = xyz/nl.norm(xyz,axis=-1)[:,:,nxs]
return np.dot(u,d0)
dx = x-parmDict['det-X']
dy = y-parmDict['det-Y']
tilt = parmDict['tilt']
dist = parmDict['dist']/npcosd(tilt) #sample-beam intersection point on detector plane
T = makeMat(tilt,0) #detector tilt matrix
R = makeMat(parmDict['phi'],2) #rotation of tilt axis matrix
MN = np.inner(R,np.inner(R,T)) #should be detector transformation matrix; why not np.inner(R,T)
d001 = np.array([0.,0.,1.]) #vector along z (beam direction); normal to untilted detector plane
r001 = np.inner(d001,MN) #should rotate vector same as detector
dxyz0 = np.inner(np.dstack([dx,dy,np.zeros_like(dx)]),MN) #transform detector pixel x,y by tilt/rotate
dxyz0 += np.array([0.,0.,dist]) #shift away from sample
ctth0 = costth(dxyz0,r001) #cos of angle between detector normal & sample-pixel vector
return npacosd(ctth0)[0]
[docs]
def GetTthPmulti(x,y,detX,detY,tilt,dist,phi):
'''Compute angle between detector normal & sample scattering ray vector
:param float/array x: detector x-position in mm
:param float/array y: detector y-position in mm
:param float/array detX: X Position of beam center
:param float/array detY: Y Position of beam center
:param float/array tilt: detector tilt (typically a single value)
:param float/array dist: detector distance
:param float/array phi: detector phi setting (typically a single value)
:returns: float/array angle: = 2-theta if tilt is zero
'''
def costth(xyz,d0):
''' compute cos of angle between vectors; xyz not normalized, d0 normalized'''
u = xyz/nl.norm(xyz,axis=-1)[:,:,nxs]
return np.dot(u,d0)
dx = detX
dy = detY
dist = dist/npcosd(tilt) #sample-beam intersection point on detector plane
T = makeMat(tilt,0) #detector tilt matrix
R = makeMat(phi,2) #rotation of tilt axis matrix
MN = np.inner(R,np.inner(R,T)) #should be detector transformation matrix; why not np.inner(R,T)
d001 = np.array([0.,0.,1.]) #vector along z (beam direction); normal to untilted detector plane
r001 = np.inner(d001,MN) #should rotate vector same as detector
dxyz0 = np.inner(np.dstack([dx,dy,np.zeros_like(dx)]),MN) #transform detector pixel x,y by tilt/rotate
dxyz0[:,:,2] += dist #shift away from sample
ctth0 = costth(dxyz0,r001) #cos of angle between detector normal & sample-pixel vector
return npacosd(ctth0)[0]
[docs]
def SamAbs(data,tax,tay,muT):
'Compute sample absorption correction for images'
if 'Cylind' in data['SampleShape']:
muR = muT*(1.+npsind(tay)**2/2.)/(npcosd(tax)) #adjust for additional thickness off sample normal
tabs = G2pwd.Absorb(data['SampleShape'],muR,tay)
elif 'Fixed' in data['SampleShape']: #assumes flat plate sample normal to beam
tabs = G2pwd.Absorb('Fixed',muT,tay)
else:
tabs = np.ones_like(tax)
return tabs
[docs]
def makeMat(Angle,Axis):
'''Make rotation matrix from Angle and Axis
:param float Angle: in degrees
:param int Axis: 0 for rotation about x, 1 for about y, etc.
'''
cs = npcosd(Angle)
ss = npsind(Angle)
M = np.array(([1.,0.,0.],[0.,cs,-ss],[0.,ss,cs]),dtype=np.float32)
return np.roll(np.roll(M,Axis,axis=0),Axis,axis=1)
def FitEllipse(xy):
def ellipse_center(p):
''' gives ellipse center coordinates
'''
b,c,d,f,a = p[1]/2., p[2], p[3]/2., p[4]/2., p[0]
num = b*b-a*c
x0=(c*d-b*f)/num
y0=(a*f-b*d)/num
return np.array([x0,y0])
def ellipse_angle_of_rotation( p ):
''' gives rotation of ellipse major axis from x-axis
range will be -90 to 90 deg
'''
b,c,a = p[1]/2., p[2], p[0]
return 0.5*npatand(2*b/(a-c))
def ellipse_axis_length( p ):
''' gives ellipse radii in [minor,major] order
'''
b,c,d,f,g,a = p[1]/2., p[2], p[3]/2., p[4]/2, p[5], p[0]
up = 2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g)
down1=(b*b-a*c)*( (c-a)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
down2=(b*b-a*c)*( (a-c)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
res1=np.sqrt(up/down1)
res2=np.sqrt(up/down2)
return np.array([ res2,res1,0.0])
xy = np.array(xy)
x = np.asarray(xy.T[0])[:,np.newaxis]
y = np.asarray(xy.T[1])[:,np.newaxis]
D = np.hstack((x*x, x*y, y*y, x, y, np.ones_like(x)))
S = np.dot(D.T,D)
C = np.zeros([6,6])
C[0,2] = C[2,0] = 2; C[1,1] = -1
E, V = nl.eig(np.dot(nl.inv(S), C))
n = np.argmax(np.abs(E))
a = V[:,n]
cent = ellipse_center(a)
phi = ellipse_angle_of_rotation(a)
radii = ellipse_axis_length(a)
phi += 90.
if radii[0] > radii[1]:
radii = [radii[1],radii[0]]
phi -= 90.
return cent,phi,radii
[docs]
def FitDetector(rings,varyList,parmDict,Print=True,covar=False):
'''Fit detector calibration parameters
:param np.array rings: vector of ring positions
:param list varyList: calibration parameters to be refined
:param dict parmDict: all calibration parameters
:param bool Print: set to True (default) to print the results
:param bool covar: set to True to return the covariance matrix (default is False)
:returns: [chisq,vals,sigList] unless covar is True, then
[chisq,vals,sigList,coVarMatrix] is returned
'''
def CalibPrint(ValSig,chisq,Npts):
print ('Image Parameters: red. chi**2: %12.3g, Np: %d'%(chisq,Npts))
ptlbls = 'names :'
ptstr = 'values:'
sigstr = 'esds :'
for name,value,sig in ValSig:
ptlbls += "%s" % (name.rjust(12))
if name == 'phi':
ptstr += Fmt[name] % (value%360.)
else:
ptstr += Fmt[name] % (value)
if sig:
sigstr += Fmt[name] % (sig)
else:
sigstr += 12*' '
print (ptlbls)
print (ptstr)
print (sigstr)
def ellipseCalcD(B,xyd,varyList,parmDict):
x,y,dsp = xyd
varyDict = dict(zip(varyList,B))
parms = {}
for parm in parmDict:
if parm in varyList:
parms[parm] = varyDict[parm]
else:
parms[parm] = parmDict[parm]
phi = parms['phi']-90. #get rotation of major axis from tilt axis
dtth = GetTthP(x,y,parmDict) #sample-detector ray wrt detector plane != 2-theta if tilted
tth = 2.0*npasind(parms['wave']/(2.*dsp))
phi0 = npatan2d(y-parms['det-Y'],x-parms['det-X'])
dxy = peneCorr(dtth,parms['dep'],parms['dist'])
stth = npsind(tth)
cosb = npcosd(parms['tilt'])
tanb = nptand(parms['tilt'])
tbm = nptand((tth-parms['tilt'])/2.)
tbp = nptand((tth+parms['tilt'])/2.)
d = parms['dist']+dxy
fplus = d*tanb*stth/(cosb+stth)
fminus = d*tanb*stth/(cosb-stth)
vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
R0 = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2. #+minor axis
R1 = (vplus+vminus)/2. #major axis
zdis = (fplus-fminus)/2.
Robs = np.sqrt((x-parms['det-X'])**2+(y-parms['det-Y'])**2)
rsqplus = R0**2+R1**2
rsqminus = R0**2-R1**2
R = rsqminus*npcosd(2.*phi0-2.*phi)+rsqplus
Q = np.sqrt(2.)*R0*R1*np.sqrt(R-2.*zdis**2*npsind(phi0-phi)**2)
P = 2.*R0**2*zdis*npcosd(phi0-phi)
Rcalc = (P+Q)/R
M = (Robs-Rcalc)*25. #why 25? does make "chi**2" more reasonable
return M
names = ['dist','det-X','det-Y','tilt','phi','dep','wave']
fmt = ['%12.3f','%12.3f','%12.3f','%12.3f','%12.3f','%12.4f','%12.6f']
Fmt = dict(zip(names,fmt))
p0 = [parmDict[key] for key in varyList]
result = leastsq(ellipseCalcD,p0,args=(rings.T,varyList,parmDict),full_output=True,ftol=1.e-8)
chisq = np.sum(result[2]['fvec']**2)/(rings.shape[0]-len(p0)) #reduced chi^2 = M/(Nobs-Nvar)
parmDict.update(zip(varyList,result[0]))
vals = list(result[0])
if not len(vals):
sig = []
ValSig = []
sigList = []
else:
sig = list(np.sqrt(chisq*np.diag(result[1])))
sigList = np.zeros(7)
for i,name in enumerate(varyList):
sigList[i] = sig[varyList.index(name)]
ValSig = zip(varyList,vals,sig)
if Print:
if len(sig):
CalibPrint(ValSig,chisq,rings.shape[0])
else:
print(' Nothing refined')
if covar:
return [chisq,vals,sigList,result[1]]
else:
return [chisq,vals,sigList]
[docs]
def FitMultiDist(rings,varyList,parmDict,Print=True,covar=False):
'''Fit detector calibration parameters with multi-distance data
:param np.array rings: vector of ring positions (x,y,dist,d-space)
:param list varyList: calibration parameters to be refined
:param dict parmDict: calibration parameters
:param bool Print: set to True (default) to print the results
:param bool covar: set to True to return the covariance matrix (default is False)
:returns: [chisq,vals,sigDict] unless covar is True, then
[chisq,vals,sigDict,coVarMatrix] is returned
'''
def CalibPrint(parmDict,sigDict,chisq,Npts):
ptlbls = 'names :'
ptstr = 'values:'
sigstr = 'esds :'
for d in sorted(set([i[5:] for i in parmDict.keys() if 'det-X' in i]),key=lambda x:int(x)):
fmt = '%12.3f'
for key in 'det-X','det-Y','delta':
name = key+d
if name not in parmDict: continue
ptlbls += "%12s" % name
ptstr += fmt % (parmDict[name])
if name in sigDict:
sigstr += fmt % (sigDict[name])
else:
sigstr += 12*' '
if len(ptlbls) > 68:
print()
print (ptlbls)
print (ptstr)
print (sigstr)
ptlbls = 'names :'
ptstr = 'values:'
sigstr = 'esds :'
if len(ptlbls) > 8:
print()
print (ptlbls)
print (ptstr)
print (sigstr)
print ('\nImage Parameters: chi**2: %12.3g, Np: %d'%(chisq,Npts))
ptlbls = 'names :'
ptstr = 'values:'
sigstr = 'esds :'
names = ['wavelength', 'dep', 'phi', 'tilt']
if 'deltaDist' in parmDict:
names += ['deltaDist']
for name in names:
if name == 'wavelength':
fmt = '%12.6f'
elif name == 'dep':
fmt = '%12.4f'
else:
fmt = '%12.3f'
ptlbls += "%s" % (name.rjust(12))
if name == 'phi':
ptstr += fmt % (parmDict[name]%360.)
else:
ptstr += fmt % (parmDict[name])
if name in sigDict:
sigstr += fmt % (sigDict[name])
else:
sigstr += 12*' '
print (ptlbls)
print (ptstr)
print (sigstr)
print()
def ellipseCalcMultiD(B,xyd,varyList,parmDict):
x,y,dist,dsp = xyd
varyDict = dict(zip(varyList,B))
parms = {}
for parm in parmDict:
if parm in varyList:
parms[parm] = varyDict[parm]
else:
parms[parm] = parmDict[parm]
# create arrays with detector center values
detX = np.array([parms['det-X'+str(int(d))] for d in dist])
detY = np.array([parms['det-Y'+str(int(d))] for d in dist])
if 'deltaDist' in parms:
deltaDist = parms['deltaDist']
else:
deltaDist = np.array([parms['delta'+str(int(d))] for d in dist])
phi = parms['phi']-90. #get rotation of major axis from tilt axis
tth = 2.0*npasind(parms['wavelength']/(2.*dsp))
dtth = GetTthPmulti(x,y,detX,detY,parms['tilt'],dist,parms['phi'])
phi0 = npatan2d(y-detY,x-detX)
dxy = peneCorr(dtth,parms['dep'],dist-deltaDist)
stth = npsind(tth)
cosb = npcosd(parms['tilt'])
tanb = nptand(parms['tilt'])
tbm = nptand((tth-parms['tilt'])/2.)
tbp = nptand((tth+parms['tilt'])/2.)
d = (dist-deltaDist)+dxy
fplus = d*tanb*stth/(cosb+stth)
fminus = d*tanb*stth/(cosb-stth)
vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
R0 = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2. #+minor axis
R1 = (vplus+vminus)/2. #major axis
zdis = (fplus-fminus)/2.
Robs = np.sqrt((x-detX)**2+(y-detY)**2)
rsqplus = R0**2+R1**2
rsqminus = R0**2-R1**2
R = rsqminus*npcosd(2.*phi0-2.*phi)+rsqplus
Q = np.sqrt(2.)*R0*R1*np.sqrt(R-2.*zdis**2*npsind(phi0-phi)**2)
P = 2.*R0**2*zdis*npcosd(phi0-phi)
Rcalc = (P+Q)/R
return (Robs-Rcalc)*25. #why 25? does make "chi**2" more reasonable
p0 = [parmDict[key] for key in varyList]
result = leastsq(ellipseCalcMultiD,p0,args=(rings.T,varyList,parmDict),full_output=True,ftol=1.e-8)
chisq = np.sum(result[2]['fvec']**2)/(rings.shape[0]-len(p0)) #reduced chi^2 = M/(Nobs-Nvar)
parmDict.update(zip(varyList,result[0]))
vals = list(result[0])
if chisq > 1:
sig = list(np.sqrt(chisq*np.diag(result[1])))
else:
sig = list(np.sqrt(np.diag(result[1])))
sigDict = {name:s for name,s in zip(varyList,sig)}
if Print:
CalibPrint(parmDict,sigDict,chisq,rings.shape[0])
if covar:
return [chisq,vals,sigDict,result[1]]
else:
return [chisq,vals,sigDict]
[docs]
def ImageLocalMax(image,w,Xpix,Ypix):
'Needs a doc string'
w2 = w*2
sizey,sizex = image.shape
xpix = int(Xpix) #get reference corner of pixel chosen
ypix = int(Ypix)
if not w:
ZMax = np.sum(image[ypix-2:ypix+2,xpix-2:xpix+2])
return xpix,ypix,ZMax,0.0001
if (w2 < xpix < sizex-w2) and (w2 < ypix < sizey-w2) and image[ypix,xpix]:
ZMax = image[ypix-w:ypix+w,xpix-w:xpix+w]
Zmax = np.argmax(ZMax)
ZMin = image[ypix-w2:ypix+w2,xpix-w2:xpix+w2]
Zmin = np.argmin(ZMin)
xpix += Zmax%w2-w
ypix += Zmax//w2-w
return xpix,ypix,np.ravel(ZMax)[Zmax],max(0.0001,np.ravel(ZMin)[Zmin]) #avoid neg/zero minimum
else:
return 0,0,0,0
[docs]
def makeRing(dsp,ellipse,pix,reject,scalex,scaley,image,mul=1):
'Needs a doc string'
def ellipseC():
'compute estimate of ellipse circumference'
if radii[0] <= 0: #hyperbola
if debug:
theta = npacosd(1./np.sqrt(1.+(radii[0]/radii[1])**2))
print ('hyperbola:',theta)
return 720.
apb = radii[1]+radii[0]
amb = radii[1]-radii[0]
return np.pi*apb*(1+3*(amb/apb)**2/(10+np.sqrt(4-3*(amb/apb)**2)))
cent,phi,radii = ellipse
if radii[1] < 0.:
return None
cphi = cosd(phi-90.) #convert to major axis rotation
sphi = sind(phi-90.)
ring = []
C = int(ellipseC())*mul #ring circumference in mm
azm = []
for i in range(0,C,1): #step around ring in 1mm increments
a = 360.*i/C
if radii[0] <= 0: #parabola or hyperbola
x = -radii[0]/cosd(a-phi) #major axis
y = radii[1]*tand(a-phi)
else:
x = radii[1]*cosd(a-phi+90.) #major axis
y = radii[0]*sind(a-phi+90.)
X = (cphi*x-sphi*y+cent[0])*scalex #convert mm to pixels
Y = (sphi*x+cphi*y+cent[1])*scaley
X,Y,I,J = ImageLocalMax(image,pix,X,Y)
if I and J and float(I)/J > reject:
X += .5 #set to center of pixel
Y += .5
X /= scalex #convert back to mm
Y /= scaley
if [X,Y,dsp] not in ring: #no duplicates!
ring.append([X,Y,dsp])
azm.append(a)
if len(ring) < 10:
ring = []
azm = []
return ring,azm
[docs]
def GetEllipse2(tth,dxy,dist,cent,tilt,phi):
'''uses Dandelin spheres to find ellipse or hyperbola parameters from detector geometry
on output
radii[0] (b-minor axis) set < 0. for hyperbola
'''
radii = [0,0,0]
stth = sind(tth)
cosb = cosd(tilt)
tanb = tand(tilt)
tbm = tand((tth-tilt)/2.)
tbp = tand((tth+tilt)/2.)
sinb = sind(tilt)
d = dist+dxy
if tth+abs(tilt) < 90.: #ellipse
fplus = d*tanb*stth/(cosb+stth)
fminus = d*tanb*stth/(cosb-stth)
vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
radii[0] = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2. #+minor axis
radii[1] = (vplus+vminus)/2. #major axis
radii[2] = tth #save for ellipse; might be useful
zdis = (fplus-fminus)/2.
else: #hyperbola!
f = d*abs(tanb)*stth/(cosb+stth)
v = d*(abs(tanb)+tand(tth-abs(tilt)))
delt = d*stth*(1.+stth*cosb)/(abs(sinb)*cosb*(stth+cosb))
eps = (v-f)/(delt-v)
radii[0] = -eps*(delt-f)/np.sqrt(eps**2-1.) #-minor axis
radii[1] = eps*(delt-f)/(eps**2-1.) #major axis
radii[2] = tth #save 2-theta for hyperbola
if tilt > 0:
zdis = f+radii[1]*eps
else:
zdis = -f-radii[1]*eps
#NB: zdis is || to major axis & phi is rotation of minor axis
#thus shift from beam to ellipse center is [Z*sin(phi),-Z*cos(phi)]
elcent = [cent[0]+zdis*sind(phi),cent[1]-zdis*cosd(phi)]
return elcent,phi,radii
[docs]
def GetEllipse(dsp,data):
'''uses Dandelin spheres to find ellipse or hyperbola parameters from detector geometry
as given in image controls dictionary (data) and a d-spacing (dsp)
'''
cent = data['center']
tilt = data['tilt']
phi = data['rotation']
dep = data.get('DetDepth',0.0)
tth = 2.0*asind(data['wavelength']/(2.*dsp))
dist = data['distance']
dxy = peneCorr(tth,dep,dist)
return GetEllipse2(tth,dxy,dist,cent,tilt,phi)
def LinePlaneCollision(planeNormal, planePoint, rayDirection, rayPoint, epsilon=1e-6):
ndotu = planeNormal.dot(rayDirection)
if ndotu < epsilon:
return None
w = rayPoint - planePoint
si = -planeNormal.dot(w) / ndotu
Psi = w + si * rayDirection + planePoint
return Psi
def GetDetSamAngle(x,y,data):
def costth(xyz,d0):
''' compute cos of angle between vectors; xyz not normalized, d0 normalized'''
u = xyz/nl.norm(xyz,axis=-1)[:,:,nxs]
return np.dot(u,d0)
#zero detector 2-theta: tested with tilted images - perfect integrations
dx = x-data['center'][0]
dy = y-data['center'][1]
tilt = data['tilt']
dist = data['distance']/npcosd(tilt) #sample-beam intersection point on detector plane
T = makeMat(tilt,0) #detector tilt matrix
R = makeMat(data['rotation'],2) #rotation of tilt axis matrix
MN = np.inner(R,np.inner(R,T)) #should be detector transformation matrix; why not np.inner(R,T)
d001 = np.array([0.,0.,1.]) #vector along z (beam direction); normal to untilted detector plane
r001 = np.inner(d001,MN) #should rotate vector same as detector
dxyz0 = np.inner(np.dstack([dx,dy,np.zeros_like(dx)]),MN) #transform detector pixel x,y by tilt/rotate
dxyz0 += np.array([0.,0.,dist]) #shift away from sample
ctth0 = costth(dxyz0,r001) #cos of angle between detector normal & sample-pixel vector
return ctth0
[docs]
def GetDetectorXY(dsp,azm,data):
'''Get detector x,y position from d-spacing (dsp), azimuth (azm,deg)
& image controls dictionary (data) - new version
'''
dist = data['distance']/npcosd(data['tilt']) #sample-beam intersection point on detector plane
cent = data['center']
phi = data['rotation']
T = makeMat(data['tilt'],0) #rotate about X
R = makeMat(phi,2) #rotate about Z
MN = np.inner(R,np.inner(R,T))
iMN= nl.inv(MN)
tth = 2.0*npasind(data['wavelength']/(2.*dsp))
vect = np.array([npsind(tth)*npcosd(azm),npsind(tth)*npsind(azm),npcosd(tth)])
dxyzN = np.inner(np.array([0.,0.,1.0]),MN) #tilt detector normal
dxyzO = np.array([0.,0.,dist]) #translate to distance
dorig = np.zeros(3)
xyz = LinePlaneCollision(dxyzN,dxyzO,vect,dorig)
if xyz is None:
return np.zeros(2)
xyz = np.inner(xyz,makeMat(data['det2theta'],1).T)
xyz -= np.array([0.,0.,dist]) #translate back
xyz = np.inner(xyz,iMN)
return np.squeeze(xyz)[:2]+cent
[docs]
def GetTthAzmDsp2(x,y,data): #expensive
'''Computes a 2theta, etc. from a detector position and calibration constants - checked
OK for ellipses & hyperbola.
Use only for detector 2-theta = 0
:returns: np.array(tth,azm,dsp) where tth is 2theta, azm is the azimutal angle,
and dsp is the d-space - not used in integration
'''
wave = data['wavelength']
cent = data['center']
tilt = data['tilt']
dist = data['distance']/cosd(tilt)
phi = data['rotation']
dep = data.get('DetDepth',0.)
azmthoff = data['azmthOff']
dx = np.array(x-cent[0],dtype=np.float32)
dy = np.array(y-cent[1],dtype=np.float32)
X = np.array(([dx,dy,np.zeros_like(dx)]),dtype=np.float32).T
X = np.dot(X,makeMat(phi,2))
Z = np.dot(X,makeMat(tilt,0)).T[2]
tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z))
dxy = peneCorr(tth,dep,data['distance'])
DX = dist-Z+dxy
DY = np.sqrt(dx**2+dy**2-Z**2)
tth = npatan2d(DY,DX)
dsp = wave/(2.*npsind(tth/2.))
azm = (npatan2d(dy,dx)+azmthoff+720.)%360.
return np.array([tth,azm,dsp])
[docs]
def GetTthAzmDsp(x,y,data): #expensive
'''Computes a 2theta, etc. from a detector position and calibration constants - checked
OK for ellipses & hyperbola.
Use for detector 2-theta != 0.
:returns: np.array(tth,azm,dsp) where tth is 2theta, azm is the azimutal angle,
and dsp is the d-space - only used for non-zero detector 2-thetas
'''
def costth(xyz):
u = xyz/nl.norm(xyz,axis=-1)[:,:,nxs]
return np.dot(u,np.array([0.,0.,1.]))
#zero detector 2-theta: tested with tilted images - perfect integrations
wave = data['wavelength']
dx = x-data['center'][0]
dy = y-data['center'][1]
tilt = data['tilt']
dist = data['distance']/npcosd(tilt) #sample-beam intersection point
T = makeMat(tilt,0)
R = makeMat(data['rotation'],2)
MN = np.inner(R,np.inner(R,T))
dxyz0 = np.inner(np.dstack([dx,dy,np.zeros_like(dx)]),MN) #correct for 45 deg tilt
dxyz0 += np.array([0.,0.,dist])
if data['DetDepth']:
ctth0 = costth(dxyz0)
tth0 = npacosd(ctth0)
dzp = peneCorr(tth0,data['DetDepth'],data['distance'])
dxyz0[:,:,2] += dzp
#non zero detector 2-theta:
if data['det2theta']:
tthMat = makeMat(data['det2theta'],1)
dxyz = np.inner(dxyz0,tthMat.T)
else:
dxyz = dxyz0
ctth = costth(dxyz)
tth = npacosd(ctth)
dsp = wave/(2.*npsind(tth/2.))
azm = (npatan2d(dxyz[:,:,1],dxyz[:,:,0])+data['azmthOff']+720.)%360.
return [tth,azm,dsp]
[docs]
def GetTth(x,y,data):
'Give 2-theta value for detector x,y position; calibration info in data'
if data['det2theta']:
return GetTthAzmDsp(x,y,data)[0]
else:
return GetTthAzmDsp2(x,y,data)[0]
[docs]
def GetTthAzm(x,y,data):
'Give 2-theta, azimuth values for detector x,y position; calibration info in data'
if data['det2theta']:
return GetTthAzmDsp(x,y,data)[0:2]
else:
return GetTthAzmDsp2(x,y,data)[0:2]
[docs]
def GetDsp(x,y,data):
'Give d-spacing value for detector x,y position; calibration info in data'
if data['det2theta']:
return GetTthAzmDsp(x,y,data)[2]
else:
return GetTthAzmDsp2(x,y,data)[2]
[docs]
def GetAzm(x,y,data):
'Give azimuth value for detector x,y position; calibration info in data'
if data['det2theta']:
return GetTthAzmDsp(x,y,data)[1]
else:
return GetTthAzmDsp2(x,y,data)[1]
# these two are used only for integration & finding pixel masks
[docs]
def GetTthAzmG(x,y,data):
'''Give 2-theta, azimuth & geometric corr. values for detector x,y position;
calibration info in data - only used in integration for detector 2-theta != 0.
checked OK for ellipses & hyperbola
This is the slow step in image integration
'''
def costth(xyz,d0):
''' compute cos of angle between vectors; xyz not normalized, d0 normalized'''
u = xyz/nl.norm(xyz,axis=-1)[:,:,nxs]
return np.dot(u,d0)
#zero detector 2-theta: tested with tilted images - perfect integrations
dx = x-data['center'][0]
dy = y-data['center'][1]
tilt = data['tilt']
dist = data['distance']/npcosd(tilt) #sample-beam intersection point on detector plane
T = makeMat(tilt,0) #detector tilt matrix
R = makeMat(data['rotation'],2) #rotation of tilt axis matrix
MN = np.inner(R,np.inner(R,T)) #should be detector transformation matrix; why not np.inner(R,T)
d001 = np.array([0.,0.,1.]) #vector along z (beam direction); normal to untilted detector plane
r001 = np.inner(d001,MN) #should rotate vector same as detector
dxyz0 = np.inner(np.dstack([dx,dy,np.zeros_like(dx)]),MN) #transform detector pixel x,y by tilt/rotate
dxyz0 += np.array([0.,0.,dist]) #shift away from sample
ctth0 = costth(dxyz0,r001) #cos of angle between detector normal & sample-pixel vector
if data['DetDepth']:
tth0 = npacosd(ctth0)
dzp = peneCorr(tth0,data['DetDepth'],data['distance'])
dxyz0[:,:,2] += dzp
#non zero detector 2-theta:
if data.get('det2theta',0):
tthMat = makeMat(data['det2theta'],1)
dxyz = np.inner(dxyz0,tthMat.T)
else:
dxyz = dxyz0
ctth = costth(dxyz,d001)
tth = npacosd(ctth)
azm = (npatan2d(dxyz[:,:,1],dxyz[:,:,0])+data['azmthOff']+720.)%360.
# new G calculation - parallax & distance (in meters)
G = 1./ctth0*np.sum(dxyz0**2,axis=2)*1.e-6 #parallax*distance^2; convert to m^2 - correct!
return tth,azm,G
def meanAzm(a,b):
AZM = lambda a,b: npacosd(0.5*(npsind(2.*b)-npsind(2.*a))/(np.pi*(b-a)/180.))/2.
azm = AZM(a,b)
return azm
[docs]
def ImageCompress(image,scale):
''' Reduces size of image by selecting every n'th point
param: image array: original image
param: scale int: intervsl between selected points
returns: array: reduced size image
'''
if scale == 1:
return image
else:
return image[::scale,::scale]
[docs]
def checkEllipse(Zsum,distSum,xSum,ySum,dist,x,y):
'Needs a doc string'
avg = np.array([distSum/Zsum,xSum/Zsum,ySum/Zsum])
curr = np.array([dist,x,y])
return abs(avg-curr)/avg < .02
def GetLineScan(image,data):
Nx,Ny = data['size']
pixelSize = data['pixelSize']
scalex = 1000./pixelSize[0] #microns --> 1/mm
scaley = 1000./pixelSize[1]
wave = data['wavelength']
numChans = data['outChannels']
LUtth = np.array(data['IOtth'],dtype=float)
azm = data['linescan'][1]-data['azmthOff']
Tx = np.array([tth for tth in np.linspace(LUtth[0],LUtth[1],numChans+1)])
Ty = np.zeros_like(Tx)
dsp = wave/(2.0*npsind(Tx/2.0))
xy = [GetDetectorXY(d,azm,data) for d in dsp]
xy = np.array(xy).T
xy[1] *= scalex
xy[0] *= scaley
xy = np.array(xy,dtype=int)
Xpix = ma.masked_outside(xy[1],0,Ny-1)
Ypix = ma.masked_outside(xy[0],0,Nx-1)
xpix = Xpix[~(Xpix.mask+Ypix.mask)].compressed()
ypix = Ypix[~(Xpix.mask+Ypix.mask)].compressed()
Ty = np.array(image[xpix,ypix],dtype=float)
Tx = ma.array(Tx,mask=Xpix.mask+Ypix.mask).compressed()
Ty /= npcosd(Tx)**2 #do parallax
Ty *= (data['distance']/1000.)**2 #do dist correction
return [Tx,Ty]
[docs]
def EdgeFinder(image,data):
'''this makes list of all x,y where I>edgeMin suitable for an ellipse search?
Not currently used but might be useful in future?
'''
import numpy.ma as ma
Nx,Ny = data['size']
pixelSize = data['pixelSize']
edgemin = data['edgemin']
scalex = pixelSize[0]/1000.
scaley = pixelSize[1]/1000.
tay,tax = np.mgrid[0:Nx,0:Ny]
tax = np.asarray(tax*scalex,dtype=np.float32)
tay = np.asarray(tay*scaley,dtype=np.float32)
tam = ma.getmask(ma.masked_less(image.flatten(),edgemin))
tax = ma.compressed(ma.array(tax.flatten(),mask=tam))
tay = ma.compressed(ma.array(tay.flatten(),mask=tam))
return zip(tax,tay)
def CalcRings(G2frame,ImageZ,data,masks):
''' '''
pixelSize = data['pixelSize']
scalex = 1000./pixelSize[0]
scaley = 1000./pixelSize[1]
data['rings'] = []
data['ellipses'] = []
if not data['calibrant']:
G2fil.G2Print ('warning: no calibration material selected')
return
skip = data['calibskip']
dmin = data['calibdmin']
if data['calibrant'] not in calFile.Calibrants:
G2fil.G2Print('Warning: %s not in local copy of image calibrants file'%data['calibrant'])
return
calibrant = calFile.Calibrants[data['calibrant']]
Bravais,SGs,Cells = calibrant[:3]
HKL = []
for bravais,sg,cell in zip(Bravais,SGs,Cells):
A = G2lat.cell2A(cell)
if sg:
SGData = G2spc.SpcGroup(sg)[1]
hkl = G2pwd.getHKLpeak(dmin,SGData,A,Inst=None,nodup=True)
HKL += list(hkl)
else:
hkl = G2lat.GenHBravais(dmin,bravais,A)
HKL += list(hkl)
if len(calibrant) > 5:
absent = calibrant[5]
else:
absent = ()
HKL = G2lat.sortHKLd(HKL,True,False)
frame = masks['Frames']
tam = ma.make_mask_none(ImageZ.shape)
if frame:
tam = ma.mask_or(tam,ma.make_mask(np.abs(polymask(data,[frame,])-255)))
for iH,H in enumerate(HKL):
if debug: print (H)
dsp = H[3]
ellipse = GetEllipse(dsp,data)
if iH not in absent and iH >= skip:
Ring = makeRing(dsp,ellipse,0,-1.,scalex,scaley,ma.array(ImageZ,mask=tam))[0]
else:
Ring = makeRing(dsp,ellipse,0,-1.,scalex,scaley,ma.array(ImageZ,mask=tam))[0]
if Ring:
if iH not in absent and iH >= skip:
data['rings'].append(np.array(Ring))
data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
[docs]
def ImageRecalibrate(G2frame,ImageZ,data,masks,getRingsOnly=False):
'''Called to repeat the calibration on an image, usually called after
calibration is done initially to improve the fit, but also
can be used after reading approximate calibration parameters,
if they are close enough that the first ring can be found.
:param G2frame: The top-level GSAS-II frame or None, to skip plotting
:param np.Array ImageZ: the image to calibrate
:param dict data: the Controls dict for the image
:param dict masks: a dict with masks
:returns: a list containing vals,varyList,sigList,parmDict,covar or rings
(with an array of x, y, and d-space values) if getRingsOnly is True
or an empty list, in case of an error
'''
if not getRingsOnly:
G2fil.G2Print ('Image recalibration:')
time0 = time.time()
pixelSize = data['pixelSize']
scalex = 1000./pixelSize[0]
scaley = 1000./pixelSize[1]
pixLimit = data['pixLimit']
cutoff = data['cutoff']
data['rings'] = []
data['ellipses'] = []
if data['DetDepth'] > 0.5: #patch - redefine DetDepth
data['DetDepth'] /= data['distance']
if not data['calibrant']:
G2fil.G2Print ('warning: no calibration material selected')
return []
skip = data['calibskip']
dmin = data['calibdmin']
if data['calibrant'] not in calFile.Calibrants:
G2fil.G2Print('Warning: %s not in local copy of image calibrants file'%data['calibrant'])
return []
calibrant = calFile.Calibrants[data['calibrant']]
Bravais,SGs,Cells = calibrant[:3]
HKL = []
for bravais,sg,cell in zip(Bravais,SGs,Cells):
A = G2lat.cell2A(cell)
if sg:
SGData = G2spc.SpcGroup(sg)[1]
hkl = G2pwd.getHKLpeak(dmin,SGData,A,Inst=None,nodup=True)
HKL += list(hkl)
else:
hkl = G2lat.GenHBravais(dmin,bravais,A)
HKL += list(hkl)
if len(calibrant) > 5:
absent = calibrant[5]
else:
absent = ()
HKL = G2lat.sortHKLd(HKL,True,False)
varyList = [item for item in data['varyList'] if data['varyList'][item]]
parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
'setdist':data.get('setdist',data['distance']),
'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
Found = False
frame = masks['Frames']
tam = ma.make_mask_none(ImageZ.shape)
if frame:
tam = ma.mask_or(tam,ma.make_mask(np.abs(polymask(data,[frame,])-255)))
hyperbola = False
for iH,H in enumerate(HKL):
if debug: print (H)
dsp = H[3]
ellipse = GetEllipse(dsp,data)
if iH not in absent and iH >= skip:
Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,ma.array(ImageZ,mask=tam))[0]
else:
Ring = makeRing(dsp,ellipse,pixLimit,1000.0,scalex,scaley,ma.array(ImageZ,mask=tam))[0]
if Ring:
if iH not in absent and iH >= skip:
data['rings'].append(np.array(Ring))
data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
if debug: print('added ellipse:',ellipse)
Found = True
elif not Found: #skipping inner rings, keep looking until ring found
continue
else: #no more rings beyond edge of detector
data['ellipses'].append([])
if Ring is None:
hyperbola = True
continue
if not data['rings']:
G2fil.G2Print ('no rings found; try lower Min ring I/Ib',mode='warn')
return []
if hyperbola:
print('Hyperbola found, outer rings not fitted')
rings = np.concatenate((data['rings']),axis=0)
if getRingsOnly:
return rings,HKL
[chisq,vals,sigList,covar] = FitDetector(rings,varyList,parmDict,True,True)
data['wavelength'] = parmDict['wave']
data['distance'] = parmDict['dist']
data['center'] = [parmDict['det-X'],parmDict['det-Y']]
data['rotation'] = np.mod(parmDict['phi'],360.0)
data['tilt'] = parmDict['tilt']
data['DetDepth'] = parmDict['dep']
data['chisq'] = chisq
N = len(data['ellipses'])
data['ellipses'] = [] #clear away individual ellipse fits
for H in HKL[:N]:
ellipse = GetEllipse(H[3],data)
data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
G2fil.G2Print ('calibration time = %.3f'%(time.time()-time0))
if G2frame:
from . import GSASIIplot as G2plt
G2plt.PlotImage(G2frame,newImage=True)
return [vals,varyList,sigList,parmDict,covar]
[docs]
def ImageCalibrate(G2frame,data):
'''Called to perform an initial image calibration after points have been
selected for the inner ring.
Called only from ``OnImRelease`` (mouse release) in
:func:`GSASIIplot.PlotImage`, thus expected to be used from GUI
only (not scripted)
'''
from . import GSASIIplot as G2plt
G2fil.G2Print ('Image calibration:')
time0 = time.time()
ring = data['ring']
pixelSize = data['pixelSize']
scalex = 1000./pixelSize[0]
scaley = 1000./pixelSize[1]
pixLimit = data['pixLimit']
cutoff = data['cutoff']
varyDict = data['varyList']
if varyDict['dist'] and varyDict['wave']:
G2fil.G2Print ('ERROR - you can not simultaneously calibrate distance and wavelength')
return False
if len(ring) < 5:
G2fil.G2Print ('ERROR - not enough inner ring points for ellipse')
return False
#fit start points on inner ring
data['ellipses'] = []
data['rings'] = []
outE = FitEllipse(ring)
fmt = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f'
fmt2 = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f, chi**2: %.3f, Np: %d'
if outE:
G2fil.G2Print (fmt%('start ellipse: ',outE[0][0],outE[0][1],outE[1],outE[2][0],outE[2][1]))
ellipse = outE
else:
return False
#setup 360 points on that ring for "good" fit
data['ellipses'].append(ellipse[:]+('g',))
Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0]
if Ring:
ellipse = FitEllipse(Ring)
Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0] #do again
ellipse = FitEllipse(Ring)
else:
G2fil.G2Print ('1st ring not sufficiently complete to proceed',mode='warn')
return False
if debug:
G2fil.G2Print (fmt2%('inner ring: ',ellipse[0][0],ellipse[0][1],ellipse[1],
ellipse[2][0],ellipse[2][1],0.,len(Ring))) #cent,phi,radii
data['ellipses'].append(ellipse[:]+('r',))
data['rings'].append(np.array(Ring))
G2plt.PlotImage(G2frame,newImage=True)
#setup for calibration
data['rings'] = []
if not data['calibrant']:
G2fil.G2Print ('Warning: no calibration material selected')
return True
skip = data['calibskip']
dmin = data['calibdmin']
#generate reflection set
calibrant = calFile.Calibrants[data['calibrant']]
Bravais,SGs,Cells = calibrant[:3]
HKL = []
for bravais,sg,cell in zip(Bravais,SGs,Cells):
A = G2lat.cell2A(cell)
if sg:
SGData = G2spc.SpcGroup(sg)[1]
hkl = G2pwd.getHKLpeak(dmin,SGData,A,Inst=None,nodup=True)
HKL += list(hkl)
else:
hkl = G2lat.GenHBravais(dmin,bravais,A)
HKL += list(hkl)
HKL = G2lat.sortHKLd(HKL,True,False)[skip:]
#set up 1st ring
elcent,phi,radii = ellipse #from fit of 1st ring
dsp = HKL[0][3]
G2fil.G2Print ('1st ring: try %.4f'%(dsp))
if varyDict['dist']:
wave = data['wavelength']
tth = 2.0*asind(wave/(2.*dsp))
else: #varyDict['wave']!
dist = data['distance']
tth = npatan2d(radii[0],dist)
data['wavelength'] = wave = 2.0*dsp*sind(tth/2.0)
Ring0 = makeRing(dsp,ellipse,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
ttth = nptand(tth)
ctth = npcosd(tth)
#1st estimate of tilt; assume ellipse - don't know sign though
if varyDict['tilt']:
tilt = npasind(np.sqrt(max(0.,1.-(radii[0]/radii[1])**2))*ctth)
if not tilt:
G2fil.G2Print ('WARNING - selected ring was fitted as a circle')
G2fil.G2Print (' - if detector was tilted we suggest you skip this ring - WARNING')
else:
tilt = data['tilt']
#1st estimate of dist: sample to detector normal to plane
if varyDict['dist']:
data['distance'] = dist = radii[0]**2/(ttth*radii[1])
else:
dist = data['distance']
if varyDict['tilt']:
#ellipse to cone axis (x-ray beam); 2 choices depending on sign of tilt
zdisp = radii[1]*ttth*tand(tilt)
zdism = radii[1]*ttth*tand(-tilt)
#cone axis position; 2 choices. Which is right?
#NB: zdisp is || to major axis & phi is rotation of minor axis
#thus shift from beam to ellipse center is [Z*sin(phi),-Z*cos(phi)]
centp = [elcent[0]+zdisp*sind(phi),elcent[1]-zdisp*cosd(phi)]
centm = [elcent[0]+zdism*sind(phi),elcent[1]-zdism*cosd(phi)]
#check get same ellipse parms either way
#now do next ring; estimate either way & do a FitDetector each way; best fit is correct one
fail = True
i2 = 1
while fail:
dsp = HKL[i2][3]
G2fil.G2Print ('2nd ring: try %.4f'%(dsp))
tth = 2.0*asind(wave/(2.*dsp))
ellipsep = GetEllipse2(tth,0.,dist,centp,tilt,phi)
G2fil.G2Print (fmt%('plus ellipse :',ellipsep[0][0],ellipsep[0][1],ellipsep[1],ellipsep[2][0],ellipsep[2][1]))
Ringp = makeRing(dsp,ellipsep,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
parmDict = {'dist':dist,'det-X':centp[0],'det-Y':centp[1],
'tilt':tilt,'phi':phi,'wave':wave,'dep':0.0}
varyList = [item for item in varyDict if varyDict[item]]
if len(Ringp) > 10:
chip = FitDetector(np.array(Ring0+Ringp),varyList,parmDict,True)[0]
tiltp = parmDict['tilt']
phip = parmDict['phi']
centp = [parmDict['det-X'],parmDict['det-Y']]
fail = False
else:
chip = 1e6
ellipsem = GetEllipse2(tth,0.,dist,centm,-tilt,phi)
G2fil.G2Print (fmt%('minus ellipse:',ellipsem[0][0],ellipsem[0][1],ellipsem[1],ellipsem[2][0],ellipsem[2][1]))
Ringm = makeRing(dsp,ellipsem,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
if len(Ringm) > 10:
parmDict['tilt'] *= -1
chim = FitDetector(np.array(Ring0+Ringm),varyList,parmDict,True)[0]
tiltm = parmDict['tilt']
phim = parmDict['phi']
centm = [parmDict['det-X'],parmDict['det-Y']]
fail = False
else:
chim = 1e6
if fail:
i2 += 1
if chip < chim:
data['tilt'] = tiltp
data['center'] = centp
data['rotation'] = phip
else:
data['tilt'] = tiltm
data['center'] = centm
data['rotation'] = phim
data['ellipses'].append(ellipsep[:]+('b',))
data['rings'].append(np.array(Ringp))
data['ellipses'].append(ellipsem[:]+('r',))
data['rings'].append(np.array(Ringm))
G2plt.PlotImage(G2frame,newImage=True)
if data['DetDepth'] > 0.5: #patch - redefine DetDepth
data['DetDepth'] /= data['distance']
parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
varyList = [item for item in varyDict if varyDict[item]]
data['rings'] = []
data['ellipses'] = []
for i,H in enumerate(HKL):
dsp = H[3]
tth = 2.0*asind(wave/(2.*dsp))
if debug: print ('HKLD:'+str(H[:4])+'2-theta: %.4f'%(tth))
elcent,phi,radii = ellipse = GetEllipse(dsp,data)
data['ellipses'].append(copy.deepcopy(ellipse+('g',)))
if debug: print (fmt%('predicted ellipse:',elcent[0],elcent[1],phi,radii[0],radii[1]))
Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0]
if Ring:
data['rings'].append(np.array(Ring))
rings = np.concatenate((data['rings']),axis=0)
if i:
chisq = FitDetector(rings,varyList,parmDict,False)[0]
data['distance'] = parmDict['dist']
data['center'] = [parmDict['det-X'],parmDict['det-Y']]
data['rotation'] = parmDict['phi']
data['tilt'] = parmDict['tilt']
data['DetDepth'] = parmDict['dep']
data['chisq'] = chisq
elcent,phi,radii = ellipse = GetEllipse(dsp,data)
if debug: print (fmt2%('fitted ellipse: ',elcent[0],elcent[1],phi,radii[0],radii[1],chisq,len(rings)))
data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
# G2plt.PlotImage(G2frame,newImage=True)
else:
if debug: print ('insufficient number of points in this ellipse to fit')
# break
G2plt.PlotImage(G2frame,newImage=True)
fullSize = len(G2frame.ImageZ)/scalex
if 2*radii[1] < .9*fullSize:
G2fil.G2Print ('Are all usable rings (>25% visible) used? Try reducing Min ring I/Ib')
N = len(data['ellipses'])
if N > 2:
FitDetector(rings,varyList,parmDict)[0]
data['wavelength'] = parmDict['wave']
data['distance'] = parmDict['dist']
data['center'] = [parmDict['det-X'],parmDict['det-Y']]
data['rotation'] = parmDict['phi']
data['tilt'] = parmDict['tilt']
data['DetDepth'] = parmDict['dep']
for H in HKL[:N]:
ellipse = GetEllipse(H[3],data)
data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
G2fil.G2Print ('calibration time = %.3f'%(time.time()-time0))
G2plt.PlotImage(G2frame,newImage=True)
return True
[docs]
def Make2ThetaAzimuthMap(data,iLim,jLim): #most expensive part of integration!
'''Makes a set of matrices that provide the 2-theta, azimuth and geometric
correction values for each pixel in an image taking into account the
detector orientation. Can be used for the entire image or a rectangular
section of an image (determined by iLim and jLim).
This is used in two ways. For image integration, the computation is done
over blocks of fixed size (typically 128 or 256 pixels) but for pixel mask
generation, the two-theta matrix for all pixels is computed. Note that
for integration, this routine will be called to generate sections as needed
or may be called by :func:`MakeUseTA`, which creates all sections at
once, so they can be reused multiple times.
:param dict data: GSAS-II image data object (describes the image)
:param list iLim: boundary along x-pixels
:param list jLim: boundary along y-pixels
:returns: TA, array[4,nI,nJ]: 2-theta, azimuth, 2 geometric corrections
'''
pixelSize = data['pixelSize']
scalex = pixelSize[0]/1000.
scaley = pixelSize[1]/1000.
tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5] #bin centers not corners
tax = np.asarray(tax*scalex,dtype=np.float32).flatten()
tay = np.asarray(tay*scaley,dtype=np.float32).flatten()
nI = iLim[1]-iLim[0]
nJ = jLim[1]-jLim[0]
TA = np.empty((4,nI,nJ))
TA[:3] = np.array(GetTthAzmG(np.reshape(tax,(nI,nJ)),np.reshape(tay,(nI,nJ)),data)) #includes geom. corr. as cos(detang)*dist**2/10^6 - most expensive step
TA[1] = np.where(TA[1]<0,TA[1]+360,TA[1])
TA[3] = G2pwd.Polarization(data['PolaVal'][0],TA[0],TA[1]-90.)[0]
return TA #2-theta, azimuth & geom. corr. arrays
[docs]
def polymask(data,Poly,Spots=[]):
''' Applies spot(point), polygon & frame masks via calls to matplotlib routines;
should be called only once each during image processing. A separate call is used for a frame.
Individual masked blocks are then pulled from the output array.
:param dict data: GSAS-II image data object (describes the image)
:param list Poly: list of polygons; if empty, returns None
:param list Spots: list of spots/points; if empty, returns None
:returns: Zimg, array[Nx,Ny] size of full image mask for all polygons/spots or frame considered
'''
import matplotlib.figure as mplfig
from matplotlib.patches import Circle
try:
from matplotlib.backends.backend_agg import FigureCanvasAgg as hcCanvas
except ImportError:
from matplotlib.backends.backend_agg import FigureCanvas as hcCanvas # standard name
if not Poly and not Spots:
return []
outmask = 'black'
inmask = 'white'
Nx,Ny = data['size']
pixelSize = data['pixelSize']
scalex = pixelSize[0]/1000.
scaley = pixelSize[1]/1000.
figure = mplfig.Figure(figsize=(Nx/400.,Ny/400.),dpi=400,facecolor=outmask)
canvas = hcCanvas(figure)
figure.clf()
ax0 = figure.add_subplot()
ax0.axis("off")
figure.subplots_adjust(bottom=0.,top=1.,left=0.,right=1.,wspace=0.,hspace=0.)
for poly in Poly:
px = np.array(poly).T[1]/scalex
py = np.array(poly).T[0]/scaley
ax0.fill(px,py,inmask)
for spot in Spots:
px = np.array(spot).T[1]/scalex
py = np.array(spot).T[0]/scaley
rad = 0.5*np.array(spot).T[2]/scaley
psp = Circle((px,py),radius=rad,fc=inmask,ec='none')
ax0.add_artist(psp)
ax0.set_xbound(0,Nx)
ax0.set_ybound(0,Ny)
img, (width,height) = canvas.print_to_buffer()
Zimg = np.frombuffer(img, np.uint8).reshape((height,width,4))
Z = np.swapaxes(Zimg[:,:,0],0,1)
A = np.flip(Z,axis=1)
return A
[docs]
def MakeMaskMap(data,masks,iLim,jLim):
'''Makes a mask array from masking parameters that are not determined by
image calibration parameters or the image intensities. Thus this uses
mask Frames, Polygons, Points, and Lines settings (but not Thresholds, Rings or
Arcs). Used on a rectangular section of an image (must be 1024x1024 or
smaller) where the size is determined
by iLim and jLim.
:param dict data: GSAS-II image data object (describes the image)
:param list iLim: boundary along x-pixels
:param list jLim: boundary along y-pixels
:returns: array[nI,nJ] TA: X, Y
'''
pixelSize = data['pixelSize']
scalex = pixelSize[0]/1000.
scaley = pixelSize[1]/1000.
frame = np.zeros(data['size'],dtype='uint8')
poly = np.zeros(data['size'],dtype='uint8')
if iLim[0] == jLim[0] == 0:
if masks['Frames']:
frame = np.abs(polymask(data,[masks['Frames'],])-255) #turn inner to outer mask
if masks['Polygons'] or masks['Points']:
poly = polymask(data,masks['Polygons'],masks['Points'])
masks['Pmask'] = frame+poly
tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5] #bin centers not corners
tax = np.asarray(tax*scalex,dtype=np.float32).flatten()
tay = np.asarray(tay*scaley,dtype=np.float32).flatten()
nI = iLim[1]-iLim[0]
nJ = jLim[1]-jLim[0]
#make position masks here
tam = ma.make_mask_none((nI*nJ))
if len(masks['Pmask']):
tam = ma.mask_or(tam,ma.make_mask(masks['Pmask'][iLim[0]:iLim[1],jLim[0]:jLim[1]].flatten()))
if tam.shape:
tam = np.reshape(tam,(nI,nJ))
else:
tam = ma.make_mask_none((nI,nJ))
for xline in masks.get('Xlines',[]): #a y pixel position
if iLim[0] <= xline <= iLim[1]:
tam[xline-iLim[0],:] = True
for yline in masks.get('Ylines',[]): #a x pixel position
if jLim[0] <= yline <= jLim[1]:
tam[:,yline-jLim[0]] = True
return tam #position mask
[docs]
def Fill2ThetaAzimuthMap(masks,TAr,tam,image,ringMask=False):
'''Makes masked intensity correction arrays that depend on image
intensity, 2theta and azimuth. Masking is generated from the
combination of the following:
an array previously generated by :func:`MakeMaskMap` combined with
Thresholds, Rings and Arcs mask input.
These correction arrays are generated for a rectangular section
of an image (must be 1024x1024 or smaller) where the size is
determined the input arrays.
Note that older, less optimized, code has been left commented out below
in case there are future problems or questions.
:param dict masks: GSAS-II mask settings
:param np.array TAr: 2theta/azimuth/correction arrays, reshaped
:param np.array tam: mask array from :func:`MakeMaskMap`
:param np.array image: image array
:returns: a list of 4 masked arrays with values for: azimuth, 2-theta,
intensity/polarization, dist**2/d0**2
'''
tax,tay,tad,pol = TAr #azimuth, 2-theta, dist**2/d0**2, pol
if ma.is_masked(image):
# get prev. masks
mask = ma.getmask(image) # Pixel mask (N.B. mask is True if pixel is masked)
mask |= tam.reshape(image.shape)
image = image.data
else:
mask = tam.reshape(image.shape)
# apply Ring & Arc masks. Note that this could be done in advance
# and be cached (like tam & TAr) but it is not clear this is needed
# often or that a lot of time is saved.
for tth,thick in masks['Rings']:
#tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.)))
mask |= ((tay > max(0.01,tth-thick/2.)) & (tay < (tth+thick/2.)))
for tth,azm,thick in masks['Arcs']:
# tamt = ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.))
# tama = ma.getmask(ma.masked_inside(tax.flatten(),azm[0],azm[1]))
# tam = ma.mask_or(tam.flatten(),tamt*tama)
mask |= ((tay > max(0.01,tth-thick/2.)) & (tay < (tth+thick/2.))
& (tax > azm[0]) & (tax < azm[1]))
# apply threshold masks
Zlim = masks['Thresholds'][1]
#taz = ma.masked_outside(image.flatten(),int(Zlim[0]),Zlim[1])
mask |= (image < Zlim[0]) | (image > Zlim[1])
#tam = ma.mask_or(tam.flatten(),ma.getmask(taz))
#tax = ma.compressed(ma.array(tax.flatten(),mask=tam)) #azimuth
#tay = ma.compressed(ma.array(tay.flatten(),mask=tam)) #2-theta
#taz = ma.compressed(ma.array(taz.flatten(),mask=tam)) #intensity
#tad = ma.compressed(ma.array(tad.flatten(),mask=tam)) #dist**2/d0**2
#tabs = ma.compressed(ma.array(tabs.flatten(),mask=tam)) #ones - later used for absorption corr.
#pol = ma.compressed(ma.array(pol.flatten(),mask=tam)) #polarization
#return tax,tay,taz/pol,tad,tabs
if ringMask: #returns contents of ring mask only
return tax[mask],tay[mask],image[mask]/pol[mask],tad[mask]
else:
mask = ~mask
return tax[mask],tay[mask],image[mask]/pol[mask],tad[mask]
[docs]
def MakeUseTA(data,blkSize=128):
'''Precomputes the set of blocked arrays for 2theta-azimuth mapping from
the controls settings of the current image for image integration.
This computation is done optionally, but provides speed as the results
from this can be cached to avoid recomputation for a series of images
with the same calibration parameters.
:param np.array data: specifies parameters for an image
:param int blkSize: a blocksize that is selected for speed
:returns: a list of TA blocks
'''
Nx,Ny = data['size']
nXBlks = (Nx-1)//blkSize+1
nYBlks = (Ny-1)//blkSize+1
useTA = []
for iBlk in range(nYBlks):
iBeg = iBlk*blkSize
iFin = min(iBeg+blkSize,Ny)
useTAj = []
for jBlk in range(nXBlks):
jBeg = jBlk*blkSize
jFin = min(jBeg+blkSize,Nx)
TA = Make2ThetaAzimuthMap(data,(iBeg,iFin),(jBeg,jFin)) #2-theta & azimuth arrays & create position mask
TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0]),ma.getdata(TA[2]),ma.getdata(TA[3]))) #azimuth, 2-theta, dist, pol
TAr = [i[:,:,0] for i in np.dsplit(TA,4)] #azimuth, 2-theta, dist**2/d0**2, pol
useTAj.append(TAr)
useTA.append(useTAj)
return useTA
[docs]
def MakeUseMask(data,masks,blkSize=128):
'''Precomputes a set of blocked mask arrays for the mask elements
that do not depend on the instrument controls (see :func:`MakeMaskMap`).
This computation is done optionally, but provides speed as the results
from this can be cached to avoid recomputation for a series of images
with the same mask parameters.
:param np.array data: specifies mask parameters for an image
:param int blkSize: a blocksize that is selected for speed
:returns: a list of TA blocks
'''
Masks = copy.deepcopy(masks)
Nx,Ny = data['size']
nXBlks = (Nx-1)//blkSize+1
nYBlks = (Ny-1)//blkSize+1
useMask = []
for iBlk in range(nYBlks):
iBeg = iBlk*blkSize
iFin = min(iBeg+blkSize,Ny)
useMaskj = []
for jBlk in range(nXBlks):
jBeg = jBlk*blkSize
jFin = min(jBeg+blkSize,Nx)
mask = MakeMaskMap(data,Masks,(iBeg,iFin),(jBeg,jFin)) #2-theta & azimuth arrays & create position mask
useMaskj.append(mask)
useMask.append(useMaskj)
return useMask
def MakeGainMap(image,Ix,Iy,data,mask,blkSize=128):
Iy *= npcosd(Ix[:-1])**2 #undo parallax
Iy /= np.array(G2pwd.Polarization(data['PolaVal'][0],Ix[:-1],0.)[0]) #undo polarization
if data['Oblique'][1]:
Iy *= G2pwd.Oblique(data['Oblique'][0],Ix[:-1]) #undo penetration
IyInt = scint.interp1d(Ix[:-1],Iy[0],bounds_error=False)
GainMap = np.zeros_like(image,dtype=float)
Mask = np.zeros_like(image,dtype=bool)
#do interpolation on all points - fills in the empty bins; leaves others the same
Nx,Ny = data['size']
nXBlks = (Nx-1)//blkSize+1
nYBlks = (Ny-1)//blkSize+1
for iBlk in range(nYBlks):
iBeg = iBlk*blkSize
iFin = min(iBeg+blkSize,Ny)
for jBlk in range(nXBlks):
jBeg = jBlk*blkSize
jFin = min(jBeg+blkSize,Nx)
TA = Make2ThetaAzimuthMap(data,(iBeg,iFin),(jBeg,jFin)) #2-theta & azimuth arrays & create position mask
tam = MakeMaskMap(data,mask,(iBeg,iFin),(jBeg,jFin))
Ipix = IyInt(TA[0])
Mask[iBeg:iFin,jBeg:jFin] = tam
GainMap[iBeg:iFin,jBeg:jFin] = np.array(image[iBeg:iFin,jBeg:jFin]/(Ipix*TA[3]))
GainMap /= np.nanmedian(GainMap)
return 1./GainMap,Mask
[docs]
def AzimuthIntegrate(image,data,masks,ringId,blkSize=1024):
''' Integrate by azimuth around the ring masked region in 0.5 deg steps
'''
if GSASIIpath.binaryPath:
import histogram2d as h2d
else:
from . import histogram2d as h2d
LRazm = np.array(data['LRazimuth'],dtype=np.float64)
numAzms = 720
ring = masks['Rings'][ringId]
Dazm = 360./numAzms
Azms = Dazm*np.arange(numAzms)
NST = np.zeros(shape=(numAzms,1),order='F',dtype=np.float32)
H0 = np.zeros(shape=(numAzms,1),order='F',dtype=np.float32)
AMasks = {'Points':[],'Rings':[ring,],'Arcs':[],'Polygons':[],'Frames':[],
'Thresholds':data['range'],'SpotMask':{'esdMul':3.,'spotMask':None}}
muT = data.get('SampleAbs',[0.0,''])[0]
LUtth = [ring[0],ring[0]+ring[1]]
dtth = ring[1]
Nx,Ny = data['size']
nXBlks = (Nx-1)//blkSize+1
nYBlks = (Ny-1)//blkSize+1
for iBlk in range(nYBlks):
iBeg = iBlk*blkSize
iFin = min(iBeg+blkSize,Ny)
for jBlk in range(nXBlks):
jBeg = jBlk*blkSize
jFin = min(jBeg+blkSize,Nx)
TA = Make2ThetaAzimuthMap(data,(iBeg,iFin),(jBeg,jFin)) #2-theta & azimuth arrays & create position mask (none here)
TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0]),ma.getdata(TA[2]),ma.getdata(TA[3]))) #azimuth, 2-theta, dist, pol
TAr = [i[:,:,0] for i in np.dsplit(TA,4)] #azimuth, 2-theta, dist**2/d0**2, pol
tam = MakeMaskMap(data,AMasks,(iBeg,iFin),(jBeg,jFin))
Block = image[iBeg:iFin,jBeg:jFin] # image Pixel mask has been applied here
tax,tay,taz,tad = Fill2ThetaAzimuthMap(AMasks,TAr,tam,Block,ringMask=True) #applies Ring masks only & returns contents
if data.get('SampleAbs',[0.0,''])[1]:
tabs = SamAbs(data,tax,tay,muT)
else:
tabs = np.ones_like(taz)
taz = np.array((taz*tad*tabs),dtype='float32')
if any([tax.shape[0],tay.shape[0],taz.shape[0]]):
NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,
numAzms,1,LRazm,LUtth,Dazm,dtth,NST,H0)
H0 /= NST
H0 = np.nan_to_num(np.array(H0))
return np.stack((Azms,H0.T[0]))
[docs]
def ImageIntegrate(image,data,masks,blkSize=128,returnN=False,useTA=None,useMask=None):
'''Integrate an image; called from ``OnIntegrate()`` and
``OnIntegrateAll()`` inside :func:`GSASIIimgGUI.UpdateImageControls`
as well as :meth:`GSASIIscriptable.G2Image.Integrate`.
:param np.array image: contains the 2-D image
:param np.array data: specifies controls/calibration parameters for an image
:param np.array masks: specifies masks parameters for an image
:param int blkSize: a blocksize that is selected for speed
:param bool returnN: If True, causes an extra matrix (NST) to be
returned. The default is False.
:param np.array useTA: contains a cached set of blocked
2theta/azimuth/correction matrices (see :func:`MakeUseTA`) for the
current image.
The default, None, causes this to be computed as needed.
:param np.array useMask: contains a cached set of blocked masks (see
:func:`MakeUseMask`) for the current image.
The default, None, causes this to be computed as needed.
:returns: list ints, azms, Xvals, cancel (or ints, azms, Xvals,
NST, cancel if returnN is True), where azms is a list of ``M`` azimuth
values that were requested for integration, ints is a list ``M`` arrays
of diffraction intensities (where each array of diffraction data is
length ``N``), Xvals is an array of "x" values, 2theta, Q, log(q)
(determined by data['binType']), also of length ``N``. Variable cancel
will always be False, since a status window is no longer supported.
'''
#for q, log(q) bins need data['binType']
if GSASIIpath.binaryPath:
import histogram2d as h2d
else:
from . import histogram2d as h2d
G2fil.G2Print ('Beginning image integration; image range: %d %d'%(np.min(image),np.max(image)))
CancelPressed = False
LUtth = np.array(data['IOtth'])
LRazm = np.array(data['LRazimuth'],dtype=np.float64)
numAzms = data['outAzimuths']
numChans = (data['outChannels']//4)*4
Dazm = (LRazm[1]-LRazm[0])/numAzms
if '2-theta' in data.get('binType','2-theta'):
lutth = LUtth
elif 'log(q)' in data['binType']:
lutth = np.log(4.*np.pi*npsind(LUtth/2.)/data['wavelength'])
elif 'q' == data['binType'].lower():
lutth = 4.*np.pi*npsind(LUtth/2.)/data['wavelength']
dtth = (lutth[1]-lutth[0])/numChans
muT = data.get('SampleAbs',[0.0,''])[0]
if data['DetDepth'] > 0.5: #patch - redefine DetDepth
data['DetDepth'] /= data['distance']
if 'SASD' in data['type']:
muT = -np.log(muT)/2. #Transmission to 1/2 thickness muT
Masks = copy.deepcopy(masks)
NST = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
H2 = np.linspace(lutth[0],lutth[1],numChans+1)
Nx,Ny = data['size']
nXBlks = (Nx-1)//blkSize+1
nYBlks = (Ny-1)//blkSize+1
tbeg = time.time()
times = [0,0,0,0,0]
for iBlk in range(nYBlks):
iBeg = iBlk*blkSize
iFin = min(iBeg+blkSize,Ny)
for jBlk in range(nXBlks):
jBeg = jBlk*blkSize
jFin = min(jBeg+blkSize,Nx)
# next is most expensive step!
t0 = time.time()
if useTA:
TAr = useTA[iBlk][jBlk]
else:
TA = Make2ThetaAzimuthMap(data,(iBeg,iFin),(jBeg,jFin)) #2-theta & azimuth arrays & create position mask
TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0]),ma.getdata(TA[2]),ma.getdata(TA[3]))) #azimuth, 2-theta, dist, pol
TAr = [i[:,:,0] for i in np.dsplit(TA,4)] #azimuth, 2-theta, cos(detang)*d**2/10^6, pol
times[1] += time.time()-t0 #xy->th,azm
t0 = time.time()
if useMask:
tam = useMask[iBlk][jBlk]
else:
tam = MakeMaskMap(data,Masks,(iBeg,iFin),(jBeg,jFin))
Block = image[iBeg:iFin,jBeg:jFin] # image Pixel mask has been applied here
tax,tay,taz,tad = Fill2ThetaAzimuthMap(Masks,TAr,tam,Block) #applies remaining masks
times[0] += time.time()-t0 # time mask application
t0 = time.time()
tax = np.where(tax > LRazm[1],tax-360.,tax) #put azm inside limits if possible
tax = np.where(tax < LRazm[0],tax+360.,tax) #are these really needed?
if data.get('SampleAbs',[0.0,''])[1]:
tabs = SamAbs(data,tax,tay,muT)
else:
tabs = np.ones_like(taz)
if 'log(q)' in data.get('binType',''):
tay = np.log(4.*np.pi*npsind(tay/2.)/data['wavelength'])
elif 'q' == data.get('binType','').lower():
tay = 4.*np.pi*npsind(tay/2.)/data['wavelength']
times[2] += time.time()-t0 #fill map
t0 = time.time()
taz = np.array((taz*tad*tabs),dtype='float32')
if any([tax.shape[0],tay.shape[0],taz.shape[0]]):
NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,
numAzms,numChans,LRazm,lutth,Dazm,dtth,NST,H0)
times[3] += time.time()-t0 #binning
G2fil.G2Print('End integration loops')
t0 = time.time()
#prepare masked arrays of bins with pixels for interpolation setup
H2msk = [ma.array(H2[:-1],mask=np.logical_not(nst)) for nst in NST]
# H0msk = [ma.array(np.divide(h0,nst),mask=np.logical_not(nst)) for nst,h0 in zip(NST,H0)]
H0msk = [ma.array(h0/nst,mask=np.logical_not(nst)) for nst,h0 in zip(NST,H0)]
#make linear interpolators; outside limits give NaN
H0 = []
problemEntries = []
for i,(h0msk,h2msk) in enumerate(zip(H0msk,H2msk)):
try:
h0int = scint.interp1d(h2msk.compressed(),h0msk.compressed(),bounds_error=False)
#do interpolation on all points - fills in the empty bins; leaves others the same
h0 = h0int(H2[:-1])
H0.append(h0)
except ValueError:
problemEntries.append(i)
H0.append(np.zeros(numChans,order='F',dtype=np.float32))
H0 = np.nan_to_num(np.array(H0))
if 'log(q)' in data.get('binType',''):
H2 = 2.*npasind(np.exp(H2)*data['wavelength']/(4.*np.pi))
elif 'q' == data.get('binType','').lower():
H2 = 2.*npasind(H2*data['wavelength']/(4.*np.pi))
if Dazm:
H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)])
else:
H1 = LRazm
if 'SASD' not in data['type']:
H0 *= np.array(G2pwd.Polarization(data['PolaVal'][0],H2[:-1],0.)[0])
if 'SASD' in data['type']:
H0 /= npcosd(H2[:-1]) #one more for small angle scattering data?
if data['Oblique'][1]:
H0 /= G2pwd.Oblique(data['Oblique'][0],H2[:-1])
times[4] += time.time()-t0 #cleanup
G2fil.G2Print ('Step times: \n apply masks %8.3fs xy->th,azm %8.3fs fill map %8.3fs \
\n binning %8.3fs cleanup %8.3fs'%(times[0],times[1],times[2],times[3],times[4]))
G2fil.G2Print ("Integration complete. Elapsed time:","%8.3fs"%(time.time()-tbeg))
if problemEntries:
msg = ""
for i in problemEntries:
if msg: msg += ', '
msg += "{:.2f}".format((H1[i+1]+H1[i])/2.)
print('\nWarning, integrations have no pixels at these azimuth(s)\n\t',msg,'\n')
if returnN: #As requested by Steven Weigand
return H0,H1,H2,NST,CancelPressed
else:
return H0,H1,H2,CancelPressed
def MakeStrStaRing(ring,Image,Controls):
ellipse = GetEllipse(ring['Dset'],Controls)
pixSize = Controls['pixelSize']
scalex = 1000./pixSize[0]
scaley = 1000./pixSize[1]
Ring = np.array(makeRing(ring['Dset'],ellipse,ring['pixLimit'],ring['cutoff'],scalex,scaley,Image)[0]).T #returns x,y,dsp for each point in ring
if len(Ring):
ring['ImxyObs'] = copy.copy(Ring[:2])
TA = GetTthAzm(Ring[0],Ring[1],Controls) #convert x,y to tth,azm
TA[0] = Controls['wavelength']/(2.*npsind(TA[0]/2.)) #convert 2th to d
ring['ImtaObs'] = TA
ring['ImtaCalc'] = np.zeros_like(ring['ImtaObs'])
Ring[0] = TA[0]
Ring[1] = TA[1]
return Ring,ring
else:
ring['ImxyObs'] = [[],[]]
ring['ImtaObs'] = [[],[]]
ring['ImtaCalc'] = [[],[]]
return [],[] #bad ring; no points found
[docs]
def FitStrSta(Image,StrSta,Controls):
'Needs a doc string'
StaControls = copy.deepcopy(Controls)
phi = StrSta['Sample phi']
wave = Controls['wavelength']
pixelSize = Controls['pixelSize']
scalex = 1000./pixelSize[0]
scaley = 1000./pixelSize[1]
StaType = StrSta['Type']
StaControls['distance'] += StrSta['Sample z']*cosd(phi)
for ring in StrSta['d-zero']: #get observed x,y,d points for the d-zeros
dset = ring['Dset']
Ring,R = MakeStrStaRing(ring,Image,StaControls)
if len(Ring):
ring.update(R)
p0 = ring['Emat']
val,esd,covMat = FitStrain(Ring,p0,dset,wave,phi,StaType)
ring['Emat'] = val
ring['Esig'] = esd
ellipse = FitEllipse(R['ImxyObs'].T)
ringxy,ringazm = makeRing(ring['Dcalc'],ellipse,0,0.,scalex,scaley,Image)
ring['ImxyCalc'] = np.array(ringxy).T[:2]
ringint = np.array([float(Image[int(x*scalex),int(y*scaley)]) for y,x in np.array(ringxy)[:,:2]])
ringint /= np.mean(ringint)
ring['Ivar'] = np.var(ringint)
ring['covMat'] = covMat
G2fil.G2Print ('Variance in normalized ring intensity: %.3f'%(ring['Ivar']))
CalcStrSta(StrSta,Controls)
def IntStrSta(Image,StrSta,Controls):
StaControls = copy.deepcopy(Controls)
pixelSize = Controls['pixelSize']
scalex = 1000./pixelSize[0]
scaley = 1000./pixelSize[1]
phi = StrSta['Sample phi']
StaControls['distance'] += StrSta['Sample z']*cosd(phi)
RingsAI = []
for ring in StrSta['d-zero']: #get observed x,y,d points for the d-zeros
Ring,R = MakeStrStaRing(ring,Image,StaControls)
if len(Ring):
ellipse = FitEllipse(R['ImxyObs'].T)
ringxy,ringazm = makeRing(ring['Dcalc'],ellipse,0,0.,scalex,scaley,Image,5)
XY = np.array(ringxy).T
Th,Azm = GetTthAzm(XY[0],XY[1],Controls)
pola = G2pwd.Polarization(Controls['PolaVal'][0],Th,Azm-90.)[0] #get pola not dpola
ring['ImxyCalc'] = np.array(ringxy).T[:2]
ringint = np.array([float(Image[int(x*scalex),int(y*scaley)]) for y,x in np.array(ringxy)[:,:2]])
ringint /= np.mean(ringint)
if Controls.get('SampleAbs',[0.0,''])[1]:
muT = Controls.get('SampleAbs',[0.0,''])[0]
tabs = SamAbs(Controls,Th,Azm,muT)
ringint *= tabs
ringint /= pola[0] #just 1st column
G2fil.G2Print (' %s %.3f %s %.3f %s %d'%('d-spacing',ring['Dcalc'],'sig(MRD):',np.sqrt(np.var(ringint)),'# points:',len(ringint)))
RingsAI.append(np.array(list(zip(ringazm,ringint))).T)
return RingsAI
def CalcStrSta(StrSta,Controls):
wave = Controls['wavelength']
phi = StrSta['Sample phi']
StaType = StrSta['Type']
for ring in StrSta['d-zero']:
Eij = ring['Emat']
E = [[Eij[0],Eij[1],0],[Eij[1],Eij[2],0],[0,0,0]]
th,azm = ring['ImtaObs']
th0 = np.ones_like(azm)*npasind(wave/(2.*ring['Dset']))
V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th0).T/1.e6,axis=2),axis=1)
if StaType == 'True':
ring['ImtaCalc'] = np.array([np.exp(V)*ring['Dset'],azm])
else:
ring['ImtaCalc'] = np.array([(V+1.)*ring['Dset'],azm])
dmin = np.min(ring['ImtaCalc'][0])
dmax = np.max(ring['ImtaCalc'][0])
if ring.get('fixDset',True):
if abs(Eij[0]) < abs(Eij[2]): #tension
ring['Dcalc'] = dmin+(dmax-dmin)/4.
else: #compression
ring['Dcalc'] = dmin+3.*(dmax-dmin)/4.
else:
ring['Dcalc'] = np.mean(ring['ImtaCalc'][0])
[docs]
def calcFij(omg,phi,azm,th):
''' Uses parameters as defined by Bob He & Kingsley Smith, Adv. in X-Ray Anal. 41, 501 (1997)
:param omg: his omega = sample omega rotation; 0 when incident beam || sample surface,
90 when perp. to sample surface
:param phi: his phi = sample phi rotation; usually = 0, axis rotates with omg.
:param azm: his chi = azimuth around incident beam
:param th: his theta = theta
'''
a = npsind(th)*npcosd(omg)+npsind(azm)*npcosd(th)*npsind(omg)
b = -npcosd(azm)*npcosd(th)
c = npsind(th)*npsind(omg)-npsind(azm)*npcosd(th)*npcosd(omg)
d = a*npsind(phi)+b*npcosd(phi)
e = a*npcosd(phi)-b*npsind(phi)
Fij = np.array([
[d**2,d*e,c*d],
[d*e,e**2,c*e],
[c*d,c*e,c**2]])
return -Fij
[docs]
def FitStrain(rings,p0,dset,wave,phi,StaType):
'Fits external strain tensor from distortion of Bragg rings in images'
def StrainPrint(ValSig,dset):
print ('Strain tensor for Dset: %.6f'%(dset))
ptlbls = 'names :'
ptstr = 'values:'
sigstr = 'esds :'
for name,fmt,value,sig in ValSig:
ptlbls += "%s" % (name.rjust(12))
ptstr += fmt % (value)
if sig:
sigstr += fmt % (sig)
else:
sigstr += 12*' '
print (ptlbls)
print (ptstr)
print (sigstr)
def strainCalc(p,xyd,dset,wave,phi,StaType):
E = np.array([[p[0],p[1],0],[p[1],p[2],0],[0,0,0]])
dspo,azm,dsp = xyd
th = npasind(wave/(2.0*dspo))
V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th).T/1.e6,axis=2),axis=1)
if StaType == 'True':
dspc = dset*np.exp(V)
else:
dspc = dset*(V+1.)
return dspo-dspc
names = ['e11','e12','e22']
fmt = ['%12.2f','%12.2f','%12.2f']
result = leastsq(strainCalc,p0,args=(rings,dset,wave,phi,StaType),full_output=True)
vals = list(result[0])
chisq = np.sum(result[2]['fvec']**2)/(rings.shape[1]-3) #reduced chi^2 = M/(Nobs-Nvar)
covM = result[1]
covMatrix = covM*chisq
sig = list(np.sqrt(chisq*np.diag(result[1])))
ValSig = zip(names,fmt,vals,sig)
StrainPrint(ValSig,dset)
return vals,sig,covMatrix
[docs]
def FitImageSpots(Image,ImMax,ind,pixSize,nxy,spotSize=1.0):
'''Used with "s" key in image plots to search for spot masks
'''
def calcMean(nxy,pixSize,img):
gdx,gdy = np.mgrid[0:nxy,0:nxy]
gdx = ma.array((gdx-nxy//2)*pixSize[0]/1000.,mask=~ma.getmaskarray(ImBox))
gdy = ma.array((gdy-nxy//2)*pixSize[1]/1000.,mask=~ma.getmaskarray(ImBox))
posx = ma.sum(gdx)/ma.count(gdx)
posy = ma.sum(gdy)/ma.count(gdy)
return posx,posy
def calcPeak(values,nxy,pixSize,img):
back,mag,px,py,sig = values
peak = np.zeros([nxy,nxy])+back
nor = 1./(2.*np.pi*sig**2)
gdx,gdy = np.mgrid[0:nxy,0:nxy]
gdx = (gdx-nxy//2)*pixSize[0]/1000.
gdy = (gdy-nxy//2)*pixSize[1]/1000.
arg = (gdx-px)**2+(gdy-py)**2
peak += mag*nor*np.exp(-arg/(2.*sig**2))
return ma.compressed(img-peak)/np.sqrt(ma.compressed(img))
def calc2Peak(values,nxy,pixSize,img):
back,mag,px,py,sigx,sigy,rho = values
peak = np.zeros([nxy,nxy])+back
nor = 1./(2.*np.pi*sigx*sigy*np.sqrt(1.-rho**2))
gdx,gdy = np.mgrid[0:nxy,0:nxy]
gdx = (gdx-nxy//2)*pixSize[0]/1000.
gdy = (gdy-nxy//2)*pixSize[1]/1000.
argnor = -1./(2.*(1.-rho**2))
arg = (gdx-px)**2/sigx**2+(gdy-py)**2/sigy**2-2.*rho*(gdx-px)*(gdy-py)/(sigx*sigy)
peak += mag*nor*np.exp(argnor*arg)
return ma.compressed(img-peak)/np.sqrt(ma.compressed(img))
nxy2 = nxy//2
ImBox = Image[ind[1]-nxy2:ind[1]+nxy2+1,ind[0]-nxy2:ind[0]+nxy2+1]
back = np.min(ImBox)
mag = np.sum(ImBox-back)
vals = [back,mag,0.,0.,0.2,0.2,0.]
ImBox = ma.array(ImBox,dtype=float,mask=ImBox>0.75*ImMax)
px = (ind[0]+.5)*pixSize[0]/1000.
py = (ind[1]+.5)*pixSize[1]/1000.
if ma.any(ma.getmaskarray(ImBox)):
vals = calcMean(nxy,pixSize,ImBox)
posx,posy = [px+vals[0],py+vals[1]]
return [posx,posy,spotSize]
else:
result = leastsq(calc2Peak,vals,args=(nxy,pixSize,ImBox),full_output=True)
if result[1] is None:
return [px,py,spotSize]
vals = result[0]
ratio = vals[4]/vals[5]
if 0.5 < ratio < 2.0 and vals[2] < 2. and vals[3] < 2.:
posx,posy = [px+vals[2],py+vals[3]]
return [posx,posy,min(6.*vals[4],spotSize)]
else:
return [px,py,spotSize]
[docs]
def TestFastPixelMask():
'''Test if the fast (C) version of Auto Pixel Masking is available.
:returns: True if the airxd.mask package can be imported; False otherwise.
'''
try:
if GSASIIpath.binaryPath:
import fmask
else:
from . import fmask
except ModuleNotFoundError:
if GSASIIpath.GetConfigValue('debug'): print('fmask not found')
return False
return True
[docs]
def FastAutoPixelMask(Image, Masks, Controls, numChans, dlg=None):
'''Find "bad" regions on an image and create a pixel mask to remove them.
This works by masking pixels that are m*sigma outside the range of the
median at that radial distance using the using the fmask C module (based on the
AIRXD C++ code https://github.com/AdvancedPhotonSource/AIRXD-ML-PUB, developed
by Howard Yanxon, Wenqian Xu and James Weng.)
This is much faster than AutoPixelMask, which does pretty much the
same computation, but uses pure Python/numpy code.
Called from GSASIIimgGUI.UpdateMasks.OnFindPixelMask (single image)
and GSASIIimgGUI.UpdateMasks.OnAutoFindPixelMask (multiple images)
[see :func:`GSASIIimgGUI.UpdateMasks`]
:param np.array Image: 2D data structure describing a diffaction image
:param dict Masks: contents of Masks data tree
:param dict Controls: diffraction & calibration parameters for image from
IMG data tree entry
:param int numChans: number of channels in eventual 2theta pattern
after integration
:returns: a bool mask array with the same shape as Image
'''
try:
if GSASIIpath.binaryPath:
import fmask
else:
from . import fmask
if GSASIIpath.GetConfigValue('debug'): print('Loaded fmask from',fmask.__file__)
except:
return None
frame = Masks['Frames']
tam = ma.make_mask_none(Image.shape)
if frame:
tam = ma.mask_or(tam,ma.make_mask(np.abs(polymask(Controls,[frame,])-255)))
ttmin = float(Masks['SpotMask'].get('SearchMin',0.0))
ttmax = float(Masks['SpotMask'].get('SearchMax',180.0))
esdMul = float(Masks['SpotMask']['esdMul'])
TA = Make2ThetaAzimuthMap(Controls, (0, Image.shape[0]), (0, Image.shape[1]))[0]
LUtth = np.array(Controls['IOtth'])
TThs = np.linspace(LUtth[0], LUtth[1], numChans, False)
# fp = open('/tmp/maskdump.pickle','wb') # make external test file
# import pickle
# pickle.dump(tam,fp) # frame mask
# pickle.dump(TA,fp) # 2theta values
# pickle.dump(Image,fp) # image values
# pickle.dump(TThs,fp) # 2theta bins
# fp.close()
print(' Fast mask: Spots greater or less than %.1f of median abs deviation are masked'%esdMul)
outMask = np.zeros_like(tam,dtype=bool).ravel()
if dlg: dlg.Update(10,"Fast scan in progress")
try:
masked = fmask.mask(esdMul, tam.ravel(), TA.ravel(), Image.ravel(), TThs, outMask, ttmin, ttmax)
except Exception as msg:
print('Exception in fmask.mask\n\t',msg)
return outMask.reshape(Image.shape)
[docs]
def AutoPixelMask(Image, Masks, Controls, numChans, dlg=None):
'''Find "bad" regions on an image and creata a pixel mask to remove them.
This works by masking pixels that are well outside the range of the
median at that radial distance.
This is ~4x faster than the original version from RBVD.
Developed by Howard Yanxon, Wenqian Xu and James Weng.
Called from GSASIIimgGUI.UpdateMasks.OnFindPixelMask (single image)
and GSASIIimgGUI.UpdateMasks.OnAutoFindPixelMask (multiple images)
[see :func:`GSASIIimgGUI.UpdateMasks`]
:param np.array Image: 2D data structure describing a diffaction image
:param dict Masks: contents of Masks data tree
:param dict Controls: diffraction & calibration parameters for image from
IMG data tree entry
:param int numChans: number of channels in eventual 2theta pattern
after integration
:param wx.Dialog dlg: a widget that can be used to show the status of
the pixel mask scan and can optionally be used to cancel the scan. If
dlg=None then this is ignored (for non-GUI use).
:returns: a mask array with the same shape as Image or None if the
the scan is cancelled from the dlg Dialog.
'''
#if GSASIIpath.GetConfigValue('debug'): print('faster all-Python AutoPixelMask')
try:
from scipy.stats import median_abs_deviation as newMAD # new in 1.5.0
def MAD(args,**kwargs):
kwargs['scale'] = 1.4826
return newMAD(args,**kwargs)
except ImportError:
try:
from scipy.stats import median_absolute_deviation as MAD
if GSASIIpath.GetConfigValue('debug'):
print('Using deprecated scipy.stats.median_absolute_deviation routine')
except:
print('Unable to load scipy.stats.median_abs_deviation')
return None
frame = Masks['Frames']
tam = ma.make_mask_none(Image.shape)
if frame:
tam = ma.mask_or(tam,ma.make_mask(np.abs(polymask(Controls,[frame,])-255)))
LUtth = np.array(Controls['IOtth'])
dtth = (LUtth[1]-LUtth[0])/numChans
esdMul = Masks['SpotMask']['esdMul']
print(' Spots greater or less than %.1f of median abs deviation are masked'%esdMul)
band = np.array(Image)
TA = Make2ThetaAzimuthMap(Controls, (0, Image.shape[0]), (0, Image.shape[1]))[0]
TThs = np.linspace(LUtth[0], LUtth[1], numChans, False)
mask = (TA >= LUtth[1]) | (TA < LUtth[0]) | tam
for it in range(len(TThs)):
masker = (TA >= TThs[it]) & (TA < TThs[it]+dtth) & ~tam
if (TThs[it] < Masks['SpotMask'].get('SearchMin',0.0) or
TThs[it] > Masks['SpotMask'].get('SearchMax',180.0)):
mask |= masker
continue
bin = band[masker]
if np.all(np.isnan(bin)):
continue
if bin.size < 1:
continue
median = np.nanmedian(bin)
mad = MAD(bin, nan_policy='omit')
anom = np.abs(band-median)/mad <= esdMul
mask |= (anom & masker)
#print(TThs[it],sum(sum(masker))- sum(sum(anom & masker)))
if not dlg is None:
GoOn = dlg.Update(it,newmsg='Processed 2-theta rings = %d'%(it))
if not GoOn[0]:
return None
#break
return ~mask
[docs]
def DoPolaCalib(ImageZ,imageData,arcTth):
''' Determine image polarization by successive integrations with & without preset arc mask.
After initial search, does a set of five with offset azimuth to get mean(std) result.
'''
from scipy.optimize import minimize_scalar
print(' Polarization fit for mask at %.1f deg 2-theta'%arcTth)
data = copy.deepcopy(imageData)
data['IOtth'] = [arcTth-2.,arcTth+2.]
data['fullIntegrate'] = True
data['LRazimuth'] = [0.,360.]
data['outChannels'] = 200
data['outAzimuths'] = 1
Arc = [arcTth,[85.,95.],2.0]
print(' Integration 2-theta test range %.1f - %.1f in 200 steps'%(data['IOtth'][0],data['IOtth'][1]))
print(' Mask azimuth range: %.1f - %.1f'%(Arc[1][0],Arc[1][1]))
Masks = {'Points':[],'Rings':[],'Arcs':[],'Polygons':[],'Frames':[],
'Thresholds':imageData['range'],'SpotMask':{'esdMul':3.,'spotMask':None}}
AMasks = {'Points':[],'Rings':[],'Arcs':[Arc,],'Polygons':[],'Frames':[],
'Thresholds':imageData['range'],'SpotMask':{'esdMul':3.,'spotMask':None}}
def func(p):
p = min(1.0,max(p,0.0))
data['PolaVal'][0] = p
H0 = ImageIntegrate(ImageZ,data,Masks)[0]
H0A = ImageIntegrate(ImageZ,data,AMasks)[0]
M = np.sum(H0)-np.sum(H0A)
print(' Polarization %.4f, fxn: %.1f'%(p,M))
return M**2
time0 = time.time()
res = minimize_scalar(func,bracket=[1.,.999],tol=.0001)
print(res)
pola = min(1.0,max(res.x,.0))
Pola = []
for arc in [75,80,85,90,95]:
Arc = [arcTth,[arc,arc+10.],2.0]
AMasks = {'Points':[],'Rings':[],'Arcs':[Arc,],'Polygons':[],'Frames':[],
'Thresholds':imageData['range'],'SpotMask':{'esdMul':3.,'spotMask':None}}
res = minimize_scalar(func,bracket=[pola-.001,pola],tol=.0001)
print(' Mask azimuth range: %.1f - %.1f'%(Arc[1][0],Arc[1][1]))
print(' pola: %.5f'%res.x)
Pola.append(res.x)
Pola = np.array(Pola)
mean = np.mean(Pola)
std = int(10000*np.std(Pola))
print(' Polarization: %.4f(%d)'%(mean,std))
print(' time: %.2fs'%(time.time()-time0))
imageData['PolaVal'][0] = mean