source: trunk/GSASIIimage.py @ 4009

Last change on this file since 4009 was 4009, checked in by toby, 3 years ago

add convariance to recalibrate all, so esd's are computed properly seq. ref. table

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 54.6 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASII image calculations: ellipse fitting & image integration       
3########### SVN repository information ###################
4# $Date: 2019-05-31 20:25:20 +0000 (Fri, 31 May 2019) $
5# $Author: toby $
6# $Revision: 4009 $
7# $URL: trunk/GSASIIimage.py $
8# $Id: GSASIIimage.py 4009 2019-05-31 20:25:20Z toby $
9########### SVN repository information ###################
10'''
11*GSASIIimage: Image calc module*
12================================
13
14Ellipse fitting & image integration
15
16'''
17from __future__ import division, print_function
18import math
19import time
20import numpy as np
21import numpy.linalg as nl
22import numpy.ma as ma
23from scipy.optimize import leastsq
24import scipy.interpolate as scint
25import copy
26import GSASIIpath
27GSASIIpath.SetVersionNumber("$Revision: 4009 $")
28try:
29    import GSASIIplot as G2plt
30except ImportError: # expected in scriptable w/o matplotlib and/or wx
31    pass
32import GSASIIlattice as G2lat
33import GSASIIpwd as G2pwd
34import GSASIIspc as G2spc
35import GSASIImath as G2mth
36
37# trig functions in degrees
38sind = lambda x: math.sin(x*math.pi/180.)
39asind = lambda x: 180.*math.asin(x)/math.pi
40tand = lambda x: math.tan(x*math.pi/180.)
41atand = lambda x: 180.*math.atan(x)/math.pi
42atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
43cosd = lambda x: math.cos(x*math.pi/180.)
44acosd = lambda x: 180.*math.acos(x)/math.pi
45rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
46#numpy versions
47npsind = lambda x: np.sin(x*np.pi/180.)
48npasind = lambda x: 180.*np.arcsin(x)/np.pi
49npcosd = lambda x: np.cos(x*np.pi/180.)
50npacosd = lambda x: 180.*np.arccos(x)/np.pi
51nptand = lambda x: np.tan(x*np.pi/180.)
52npatand = lambda x: 180.*np.arctan(x)/np.pi
53npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
54nxs = np.newaxis
55debug = False
56   
57def pointInPolygon(pXY,xy):
58    'Needs a doc string'
59    #pXY - assumed closed 1st & last points are duplicates
60    Inside = False
61    N = len(pXY)
62    p1x,p1y = pXY[0]
63    for i in range(N+1):
64        p2x,p2y = pXY[i%N]
65        if (max(p1y,p2y) >= xy[1] > min(p1y,p2y)) and (xy[0] <= max(p1x,p2x)):
66            if p1y != p2y:
67                xinters = (xy[1]-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
68            if p1x == p2x or xy[0] <= xinters:
69                Inside = not Inside
70        p1x,p1y = p2x,p2y
71    return Inside
72   
73def peneCorr(tth,dep,dist,tilt=0.,azm=0.):
74    'Needs a doc string'
75#    return dep*(1.-npcosd(abs(tilt*npsind(azm))-tth*npcosd(azm)))  #something wrong here
76    return dep*(1.-npcosd(tth))*dist**2/1000.         #best one
77#    return dep*npsind(tth)             #not as good as 1-cos2Q
78       
79def makeMat(Angle,Axis):
80    '''Make rotation matrix from Angle and Axis
81
82    :param float Angle: in degrees
83    :param int Axis: 0 for rotation about x, 1 for about y, etc.
84    '''
85    cs = npcosd(Angle)
86    ss = npsind(Angle)
87    M = np.array(([1.,0.,0.],[0.,cs,-ss],[0.,ss,cs]),dtype=np.float32)
88    return np.roll(np.roll(M,Axis,axis=0),Axis,axis=1)
89                   
90def FitEllipse(xy):
91   
92    def ellipse_center(p):
93        ''' gives ellipse center coordinates
94        '''
95        b,c,d,f,a = p[1]/2., p[2], p[3]/2., p[4]/2., p[0]
96        num = b*b-a*c
97        x0=(c*d-b*f)/num
98        y0=(a*f-b*d)/num
99        return np.array([x0,y0])
100   
101    def ellipse_angle_of_rotation( p ):
102        ''' gives rotation of ellipse major axis from x-axis
103        range will be -90 to 90 deg
104        '''
105        b,c,a = p[1]/2., p[2], p[0]
106        return 0.5*npatand(2*b/(a-c))
107   
108    def ellipse_axis_length( p ):
109        ''' gives ellipse radii in [minor,major] order
110        '''
111        b,c,d,f,g,a = p[1]/2., p[2], p[3]/2., p[4]/2, p[5], p[0]
112        up = 2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g)
113        down1=(b*b-a*c)*( (c-a)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
114        down2=(b*b-a*c)*( (a-c)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
115        res1=np.sqrt(up/down1)
116        res2=np.sqrt(up/down2)
117        return np.array([ res2,res1])
118   
119    xy = np.array(xy)
120    x = np.asarray(xy.T[0])[:,np.newaxis]
121    y = np.asarray(xy.T[1])[:,np.newaxis]
122    D =  np.hstack((x*x, x*y, y*y, x, y, np.ones_like(x)))
123    S = np.dot(D.T,D)
124    C = np.zeros([6,6])
125    C[0,2] = C[2,0] = 2; C[1,1] = -1
126    E, V =  nl.eig(np.dot(nl.inv(S), C))
127    n = np.argmax(np.abs(E))
128    a = V[:,n]
129    cent = ellipse_center(a)
130    phi = ellipse_angle_of_rotation(a)
131    radii = ellipse_axis_length(a)
132    phi += 90.
133    if radii[0] > radii[1]:
134        radii = [radii[1],radii[0]]
135        phi -= 90.
136    return cent,phi,radii
137
138def FitDetector(rings,varyList,parmDict,Print=True,covar=False):
139    '''Fit detector calibration paremeters
140
141    :param np.array rings: vector of ring positions
142    :param list varyList: calibration parameters to be refined
143    :param dict parmDict: all calibration parameters
144    :param bool Print: set to True (default) to print the results
145    :param bool covar: set to True to return the covariance matrix (default is False)
146    :returns: [chisq,vals,sigList] unless covar is True, then
147        [chisq,vals,sigList,coVarMatrix] is returned
148    '''
149       
150    def CalibPrint(ValSig,chisq,Npts):
151        print ('Image Parameters: chi**2: %12.3g, Np: %d'%(chisq,Npts))
152        ptlbls = 'names :'
153        ptstr =  'values:'
154        sigstr = 'esds  :'
155        for name,value,sig in ValSig:
156            ptlbls += "%s" % (name.rjust(12))
157            if name == 'phi':
158                ptstr += Fmt[name] % (value%360.)
159            else:
160                ptstr += Fmt[name] % (value)
161            if sig:
162                sigstr += Fmt[name] % (sig)
163            else:
164                sigstr += 12*' '
165        print (ptlbls)
166        print (ptstr)
167        print (sigstr)       
168       
169    def ellipseCalcD(B,xyd,varyList,parmDict):
170       
171        x,y,dsp = xyd
172        varyDict = dict(zip(varyList,B))
173        parms = {}
174        for parm in parmDict:
175            if parm in varyList:
176                parms[parm] = varyDict[parm]
177            else:
178                parms[parm] = parmDict[parm]
179        phi = parms['phi']-90.               #get rotation of major axis from tilt axis
180        tth = 2.0*npasind(parms['wave']/(2.*dsp))
181        phi0 = npatan2d(y-parms['det-Y'],x-parms['det-X'])
182        dxy = peneCorr(tth,parms['dep'],parms['dist'],parms['tilt'],phi0)
183        stth = npsind(tth)
184        cosb = npcosd(parms['tilt'])
185        tanb = nptand(parms['tilt'])       
186        tbm = nptand((tth-parms['tilt'])/2.)
187        tbp = nptand((tth+parms['tilt'])/2.)
188        d = parms['dist']+dxy
189        fplus = d*tanb*stth/(cosb+stth)
190        fminus = d*tanb*stth/(cosb-stth)
191        vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
192        vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
193        R0 = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2.      #+minor axis
194        R1 = (vplus+vminus)/2.                                    #major axis
195        zdis = (fplus-fminus)/2.
196        Robs = np.sqrt((x-parms['det-X'])**2+(y-parms['det-Y'])**2)
197        rsqplus = R0**2+R1**2
198        rsqminus = R0**2-R1**2
199        R = rsqminus*npcosd(2.*phi0-2.*phi)+rsqplus
200        Q = np.sqrt(2.)*R0*R1*np.sqrt(R-2.*zdis**2*npsind(phi0-phi)**2)
201        P = 2.*R0**2*zdis*npcosd(phi0-phi)
202        Rcalc = (P+Q)/R
203        M = (Robs-Rcalc)*10.        #why 10? does make "chi**2" more reasonable
204        return M
205       
206    names = ['dist','det-X','det-Y','tilt','phi','dep','wave']
207    fmt = ['%12.3f','%12.3f','%12.3f','%12.3f','%12.3f','%12.4f','%12.6f']
208    Fmt = dict(zip(names,fmt))
209    p0 = [parmDict[key] for key in varyList]
210    result = leastsq(ellipseCalcD,p0,args=(rings.T,varyList,parmDict),full_output=True,ftol=1.e-8)
211    chisq = np.sum(result[2]['fvec']**2)/(rings.shape[0]-len(p0))   #reduced chi^2 = M/(Nobs-Nvar)
212    parmDict.update(zip(varyList,result[0]))
213    vals = list(result[0])
214    sig = list(np.sqrt(chisq*np.diag(result[1])))
215    sigList = np.zeros(7)
216    for i,name in enumerate(varyList):
217        sigList[i] = sig[varyList.index(name)]
218    ValSig = zip(varyList,vals,sig)
219    if Print:
220        CalibPrint(ValSig,chisq,rings.shape[0])
221    if covar:
222        return [chisq,vals,sigList,result[1]]
223    else:
224        return [chisq,vals,sigList]
225
226def ImageLocalMax(image,w,Xpix,Ypix):
227    'Needs a doc string'
228    w2 = w*2
229    sizey,sizex = image.shape
230    xpix = int(Xpix)            #get reference corner of pixel chosen
231    ypix = int(Ypix)
232    if not w:
233        ZMax = np.sum(image[ypix-2:ypix+2,xpix-2:xpix+2])
234        return xpix,ypix,ZMax,0.0001
235    if (w2 < xpix < sizex-w2) and (w2 < ypix < sizey-w2) and image[ypix,xpix]:
236        ZMax = image[ypix-w:ypix+w,xpix-w:xpix+w]
237        Zmax = np.argmax(ZMax)
238        ZMin = image[ypix-w2:ypix+w2,xpix-w2:xpix+w2]
239        Zmin = np.argmin(ZMin)
240        xpix += Zmax%w2-w
241        ypix += Zmax//w2-w
242        return xpix,ypix,np.ravel(ZMax)[Zmax],max(0.0001,np.ravel(ZMin)[Zmin])   #avoid neg/zero minimum
243    else:
244        return 0,0,0,0     
245   
246def makeRing(dsp,ellipse,pix,reject,scalex,scaley,image,mul=1):
247    'Needs a doc string'
248    def ellipseC():
249        'compute estimate of ellipse circumference'
250        if radii[0] < 0:        #hyperbola
251#            theta = npacosd(1./np.sqrt(1.+(radii[0]/radii[1])**2))
252#            print (theta)
253            return 0
254        apb = radii[1]+radii[0]
255        amb = radii[1]-radii[0]
256        return np.pi*apb*(1+3*(amb/apb)**2/(10+np.sqrt(4-3*(amb/apb)**2)))
257       
258    cent,phi,radii = ellipse
259    cphi = cosd(phi-90.)        #convert to major axis rotation
260    sphi = sind(phi-90.)
261    ring = []
262    C = int(ellipseC())*mul         #ring circumference in mm
263    azm = []
264    for i in range(0,C,1):      #step around ring in 1mm increments
265        a = 360.*i/C
266        x = radii[1]*cosd(a-phi+90.)        #major axis
267        y = radii[0]*sind(a-phi+90.)
268        X = (cphi*x-sphi*y+cent[0])*scalex      #convert mm to pixels
269        Y = (sphi*x+cphi*y+cent[1])*scaley
270        X,Y,I,J = ImageLocalMax(image,pix,X,Y)
271        if I and J and float(I)/J > reject:
272            X += .5                             #set to center of pixel
273            Y += .5
274            X /= scalex                         #convert back to mm
275            Y /= scaley
276            if [X,Y,dsp] not in ring:           #no duplicates!
277                ring.append([X,Y,dsp])
278                azm.append(a)
279    if len(ring) < 10:
280        ring = []
281        azm = []
282    return ring,azm
283   
284def GetEllipse2(tth,dxy,dist,cent,tilt,phi):
285    '''uses Dandelin spheres to find ellipse or hyperbola parameters from detector geometry
286    on output
287    radii[0] (b-minor axis) set < 0. for hyperbola
288   
289    '''
290    radii = [0,0]
291    stth = sind(tth)
292    cosb = cosd(tilt)
293    tanb = tand(tilt)
294    tbm = tand((tth-tilt)/2.)
295    tbp = tand((tth+tilt)/2.)
296    sinb = sind(tilt)
297    d = dist+dxy
298    if tth+abs(tilt) < 90.:      #ellipse
299        fplus = d*tanb*stth/(cosb+stth)
300        fminus = d*tanb*stth/(cosb-stth)
301        vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
302        vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
303        radii[0] = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2.      #+minor axis
304        radii[1] = (vplus+vminus)/2.                                    #major axis
305        zdis = (fplus-fminus)/2.
306    else:   #hyperbola!
307        f = d*abs(tanb)*stth/(cosb+stth)
308        v = d*(abs(tanb)+tand(tth-abs(tilt)))
309        delt = d*stth*(1.+stth*cosb)/(abs(sinb)*cosb*(stth+cosb))
310        eps = (v-f)/(delt-v)
311        radii[0] = -eps*(delt-f)/np.sqrt(eps**2-1.)                     #-minor axis
312        radii[1] = eps*(delt-f)/(eps**2-1.)                             #major axis
313        if tilt > 0:
314            zdis = f+radii[1]*eps
315        else:
316            zdis = -f
317#NB: zdis is || to major axis & phi is rotation of minor axis
318#thus shift from beam to ellipse center is [Z*sin(phi),-Z*cos(phi)]
319    elcent = [cent[0]+zdis*sind(phi),cent[1]-zdis*cosd(phi)]
320    return elcent,phi,radii
321   
322def GetEllipse(dsp,data):
323    '''uses Dandelin spheres to find ellipse or hyperbola parameters from detector geometry
324    as given in image controls dictionary (data) and a d-spacing (dsp)
325    '''
326    cent = data['center']
327    tilt = data['tilt']
328    phi = data['rotation']
329    dep = data.get('DetDepth',0.0)
330    tth = 2.0*asind(data['wavelength']/(2.*dsp))
331    dist = data['distance']
332    dxy = peneCorr(tth,dep,dist,tilt)
333    return GetEllipse2(tth,dxy,dist,cent,tilt,phi)
334       
335def GetDetectorXY(dsp,azm,data):
336    '''Get detector x,y position from d-spacing (dsp), azimuth (azm,deg)
337    & image controls dictionary (data)
338    it seems to be only used in plotting
339    '''   
340    elcent,phi,radii = GetEllipse(dsp,data)
341    phi = data['rotation']-90.          #to give rotation of major axis
342    tilt = data['tilt']
343    dist = data['distance']
344    cent = data['center']
345    tth = 2.0*asind(data['wavelength']/(2.*dsp))
346    stth = sind(tth)
347    cosb = cosd(tilt)
348    if radii[0] > 0.:
349        sinb = sind(tilt)
350        tanb = tand(tilt)
351        fplus = dist*tanb*stth/(cosb+stth)
352        fminus = dist*tanb*stth/(cosb-stth)
353        zdis = (fplus-fminus)/2.
354        rsqplus = radii[0]**2+radii[1]**2
355        rsqminus = radii[0]**2-radii[1]**2
356        R = rsqminus*cosd(2.*azm-2.*phi)+rsqplus
357        Q = np.sqrt(2.)*radii[0]*radii[1]*np.sqrt(R-2.*zdis**2*sind(azm-phi)**2)
358        P = 2.*radii[0]**2*zdis*cosd(azm-phi)
359        radius = (P+Q)/R
360        xy = np.array([radius*cosd(azm),radius*sind(azm)])
361        xy += cent
362    else:   #hyperbola - both branches (one is way off screen!)
363        sinb = abs(sind(tilt))
364        tanb = abs(tand(tilt))
365        f = dist*tanb*stth/(cosb+stth)
366        v = dist*(tanb+tand(tth-abs(tilt)))
367        delt = dist*stth*(1+stth*cosb)/(sinb*cosb*(stth+cosb))
368        ecc = (v-f)/(delt-v)
369        R = radii[1]*(ecc**2-1)/(1-ecc*cosd(azm))
370        if tilt > 0.:
371            offset = 2.*radii[1]*ecc+f      #select other branch
372            xy = [-R*cosd(azm)-offset,-R*sind(azm)]
373        else:
374            offset = -f
375            xy = [-R*cosd(azm)-offset,R*sind(azm)]
376        xy = -np.array([xy[0]*cosd(phi)+xy[1]*sind(phi),xy[0]*sind(phi)-xy[1]*cosd(phi)])
377        xy += cent
378    return xy
379   
380def GetDetXYfromThAzm(Th,Azm,data):
381    '''Computes a detector position from a 2theta angle and an azimultal
382    angle (both in degrees) - apparently not used!
383    '''
384    dsp = data['wavelength']/(2.0*npsind(Th))   
385    return GetDetectorXY(dsp,Azm,data)
386                   
387def GetTthAzmDsp(x,y,data): #expensive
388    '''Computes a 2theta, etc. from a detector position and calibration constants - checked
389    OK for ellipses & hyperbola.
390
391    :returns: np.array(tth,azm,G,dsp) where tth is 2theta, azm is the azimutal angle,
392       G is ? and dsp is the d-space
393    '''
394    wave = data['wavelength']
395    cent = data['center']
396    tilt = data['tilt']
397    dist = data['distance']/cosd(tilt)
398    x0 = dist*tand(tilt)
399    phi = data['rotation']
400    dep = data.get('DetDepth',0.)
401    azmthoff = data['azmthOff']
402    dx = np.array(x-cent[0],dtype=np.float32)
403    dy = np.array(y-cent[1],dtype=np.float32)
404    D = ((dx-x0)**2+dy**2+dist**2)      #sample to pixel distance
405    X = np.array(([dx,dy,np.zeros_like(dx)]),dtype=np.float32).T
406    X = np.dot(X,makeMat(phi,2))
407    Z = np.dot(X,makeMat(tilt,0)).T[2]
408    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z))
409    dxy = peneCorr(tth,dep,dist,tilt,npatan2d(dy,dx))
410    DX = dist-Z+dxy
411    DY = np.sqrt(dx**2+dy**2-Z**2)
412    tth = npatan2d(DY,DX) 
413    dsp = wave/(2.*npsind(tth/2.))
414    azm = (npatan2d(dy,dx)+azmthoff+720.)%360.
415    G = D/dist**2       #for geometric correction = 1/cos(2theta)^2 if tilt=0.
416    return np.array([tth,azm,G,dsp])
417   
418def GetTth(x,y,data):
419    'Give 2-theta value for detector x,y position; calibration info in data'
420    return GetTthAzmDsp(x,y,data)[0]
421   
422def GetTthAzm(x,y,data):
423    'Give 2-theta, azimuth values for detector x,y position; calibration info in data'
424    return GetTthAzmDsp(x,y,data)[0:2]
425   
426def GetTthAzmG(x,y,data):
427    '''Give 2-theta, azimuth & geometric corr. values for detector x,y position;
428     calibration info in data - only used in integration
429    '''
430    'Needs a doc string - checked OK for ellipses & hyperbola'
431    tilt = data['tilt']
432    dist = data['distance']/npcosd(tilt)
433    x0 = data['distance']*nptand(tilt)
434    MN = -np.inner(makeMat(data['rotation'],2),makeMat(tilt,0))
435    distsq = data['distance']**2
436    dx = x-data['center'][0]
437    dy = y-data['center'][1]
438    G = ((dx-x0)**2+dy**2+distsq)/distsq       #for geometric correction = 1/cos(2theta)^2 if tilt=0.
439    Z = np.dot(np.dstack([dx.T,dy.T,np.zeros_like(dx.T)]),MN).T[2]
440    xyZ = dx**2+dy**2-Z**2   
441    tth = npatand(np.sqrt(xyZ)/(dist-Z))
442    dxy = peneCorr(tth,data['DetDepth'],dist,tilt,npatan2d(dy,dx))
443    tth = npatan2d(np.sqrt(xyZ),dist-Z+dxy) 
444    azm = (npatan2d(dy,dx)+data['azmthOff']+720.)%360.
445    return tth,azm,G
446
447def GetDsp(x,y,data):
448    'Give d-spacing value for detector x,y position; calibration info in data'
449    return GetTthAzmDsp(x,y,data)[3]
450       
451def GetAzm(x,y,data):
452    'Give azimuth value for detector x,y position; calibration info in data'
453    return GetTthAzmDsp(x,y,data)[1]
454   
455def meanAzm(a,b):
456    AZM = lambda a,b: npacosd(0.5*(npsind(2.*b)-npsind(2.*a))/(np.pi*(b-a)/180.))/2.
457    azm = AZM(a,b)
458#    quad = int((a+b)/180.)
459#    if quad == 1:
460#        azm = 180.-azm
461#    elif quad == 2:
462#        azm += 180.
463#    elif quad == 3:
464#        azm = 360-azm
465    return azm     
466       
467def ImageCompress(image,scale):
468    ''' Reduces size of image by selecting every n'th point
469    param: image array: original image
470    param: scale int: intervsl between selected points
471    returns: array: reduced size image
472    '''
473    if scale == 1:
474        return image
475    else:
476        return image[::scale,::scale]
477       
478def checkEllipse(Zsum,distSum,xSum,ySum,dist,x,y):
479    'Needs a doc string'
480    avg = np.array([distSum/Zsum,xSum/Zsum,ySum/Zsum])
481    curr = np.array([dist,x,y])
482    return abs(avg-curr)/avg < .02
483
484def GetLineScan(image,data):
485    Nx,Ny = data['size']
486    pixelSize = data['pixelSize']
487    scalex = 1000./pixelSize[0]         #microns --> 1/mm
488    scaley = 1000./pixelSize[1]
489    wave = data['wavelength']
490    numChans = data['outChannels']
491    LUtth = np.array(data['IOtth'],dtype=np.float)
492    azm = data['linescan'][1]-data['azmthOff']
493    Tx = np.array([tth for tth in np.linspace(LUtth[0],LUtth[1],numChans+1)])
494    Ty = np.zeros_like(Tx)
495    dsp = wave/(2.0*npsind(Tx/2.0))
496    xy = np.array([GetDetectorXY(d,azm,data) for d in dsp]).T
497    xy[1] *= scalex
498    xy[0] *= scaley
499    xy = np.array(xy,dtype=int)
500    Xpix = ma.masked_outside(xy[1],0,Ny-1)
501    Ypix = ma.masked_outside(xy[0],0,Nx-1)
502    xpix = Xpix[~(Xpix.mask+Ypix.mask)].compressed()
503    ypix = Ypix[~(Xpix.mask+Ypix.mask)].compressed()
504    Ty = image[xpix,ypix]
505    Tx = ma.array(Tx,mask=Xpix.mask+Ypix.mask).compressed()
506    return [Tx,Ty]
507
508def EdgeFinder(image,data):
509    '''this makes list of all x,y where I>edgeMin suitable for an ellipse search?
510    Not currently used but might be useful in future?
511    '''
512    import numpy.ma as ma
513    Nx,Ny = data['size']
514    pixelSize = data['pixelSize']
515    edgemin = data['edgemin']
516    scalex = pixelSize[0]/1000.
517    scaley = pixelSize[1]/1000.   
518    tay,tax = np.mgrid[0:Nx,0:Ny]
519    tax = np.asfarray(tax*scalex,dtype=np.float32)
520    tay = np.asfarray(tay*scaley,dtype=np.float32)
521    tam = ma.getmask(ma.masked_less(image.flatten(),edgemin))
522    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))
523    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))
524    return zip(tax,tay)
525   
526def MakeFrameMask(data,frame):
527    import polymask as pm
528    pixelSize = data['pixelSize']
529    scalex = pixelSize[0]/1000.
530    scaley = pixelSize[1]/1000.
531    blkSize = 512
532    Nx,Ny = data['size']
533    nXBlks = (Nx-1)//blkSize+1
534    nYBlks = (Ny-1)//blkSize+1
535    tam = ma.make_mask_none(data['size'])
536    for iBlk in range(nXBlks):
537        iBeg = iBlk*blkSize
538        iFin = min(iBeg+blkSize,Nx)
539        for jBlk in range(nYBlks):
540            jBeg = jBlk*blkSize
541            jFin = min(jBeg+blkSize,Ny)               
542            nI = iFin-iBeg
543            nJ = jFin-jBeg
544            tax,tay = np.mgrid[iBeg+0.5:iFin+.5,jBeg+.5:jFin+.5]         #bin centers not corners
545            tax = np.asfarray(tax*scalex,dtype=np.float32)
546            tay = np.asfarray(tay*scaley,dtype=np.float32)
547            tamp = ma.make_mask_none((1024*1024))
548            tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
549                tay.flatten(),len(frame),frame,tamp)[:nI*nJ])-True  #switch to exclude around frame
550            if tamp.shape:
551                tamp = np.reshape(tamp[:nI*nJ],(nI,nJ))
552                tam[iBeg:iFin,jBeg:jFin] = ma.mask_or(tamp[0:nI,0:nJ],tam[iBeg:iFin,jBeg:jFin])
553            else:
554                tam[iBeg:iFin,jBeg:jFin] = True
555    return tam.T
556   
557def ImageRecalibrate(G2frame,ImageZ,data,masks):
558    '''Called to repeat the calibration on an image, usually called after
559    calibration is done initially to improve the fit.
560
561    :param G2frame: The top-level GSAS-II frame or None, to skip plotting
562    :param np.Array ImageZ: the image to calibrate
563    :param dict data: the Controls dict for the image
564    :param dict masks: a dict with masks
565    :returns: a list containing vals,varyList,sigList,parmDict,covar
566    '''
567    import ImageCalibrants as calFile
568    print ('Image recalibration:')
569    time0 = time.time()
570    pixelSize = data['pixelSize']
571    scalex = 1000./pixelSize[0]
572    scaley = 1000./pixelSize[1]
573    pixLimit = data['pixLimit']
574    cutoff = data['cutoff']
575    data['rings'] = []
576    data['ellipses'] = []
577    if data['DetDepth'] > 0.5:          #patch - redefine DetDepth
578        data['DetDepth'] /= data['distance']
579    if not data['calibrant']:
580        print ('no calibration material selected')
581        return []   
582    skip = data['calibskip']
583    dmin = data['calibdmin']
584    if data['calibrant'] not in calFile.Calibrants:
585        print(' %s not in local copy of image calibrants file'%data['calibrant'])
586        return []
587    Bravais,SGs,Cells = calFile.Calibrants[data['calibrant']][:3]
588    HKL = []
589    for bravais,sg,cell in zip(Bravais,SGs,Cells):
590        A = G2lat.cell2A(cell)
591        if sg:
592            SGData = G2spc.SpcGroup(sg)[1]
593            hkl = G2pwd.getHKLpeak(dmin,SGData,A)
594            HKL += list(hkl)
595        else:
596            hkl = G2lat.GenHBravais(dmin,bravais,A)
597            HKL += list(hkl)
598    HKL = G2lat.sortHKLd(HKL,True,False)
599    varyList = [item for item in data['varyList'] if data['varyList'][item]]
600    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
601        'setdist':data.get('setdist',data['distance']),
602        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
603    Found = False
604    wave = data['wavelength']
605    frame = masks['Frames']
606    tam = ma.make_mask_none(ImageZ.shape)
607    if frame:
608        tam = ma.mask_or(tam,MakeFrameMask(data,frame))
609    for iH,H in enumerate(HKL):
610        if debug:   print (H) 
611        dsp = H[3]
612        tth = 2.0*asind(wave/(2.*dsp))
613        if tth+abs(data['tilt']) > 90.:
614            print ('next line is a hyperbola - search stopped')
615            break
616        ellipse = GetEllipse(dsp,data)
617        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,ma.array(ImageZ,mask=tam))[0]
618        if Ring:
619            if iH >= skip:
620                data['rings'].append(np.array(Ring))
621            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
622            Found = True
623        elif not Found:         #skipping inner rings, keep looking until ring found
624            continue
625        else:                   #no more rings beyond edge of detector
626            data['ellipses'].append([])
627            continue
628    if not data['rings']:
629        print ('no rings found; try lower Min ring I/Ib')
630        return []   
631       
632    rings = np.concatenate((data['rings']),axis=0)
633    [chisq,vals,sigList,covar] = FitDetector(rings,varyList,parmDict,True,True)
634    data['wavelength'] = parmDict['wave']
635    data['distance'] = parmDict['dist']
636    data['center'] = [parmDict['det-X'],parmDict['det-Y']]
637    data['rotation'] = np.mod(parmDict['phi'],360.0)
638    data['tilt'] = parmDict['tilt']
639    data['DetDepth'] = parmDict['dep']
640    data['chisq'] = chisq
641    N = len(data['ellipses'])
642    data['ellipses'] = []           #clear away individual ellipse fits
643    for H in HKL[:N]:
644        ellipse = GetEllipse(H[3],data)
645        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))   
646    print ('calibration time = %.3f'%(time.time()-time0))
647    if G2frame:
648        G2plt.PlotImage(G2frame,newImage=True)       
649    return [vals,varyList,sigList,parmDict,covar]
650
651def ImageCalibrate(G2frame,data):
652    '''Called to perform an initial image calibration after points have been
653    selected for the inner ring.
654    '''
655    import ImageCalibrants as calFile
656    print ('Image calibration:')
657    time0 = time.time()
658    ring = data['ring']
659    pixelSize = data['pixelSize']
660    scalex = 1000./pixelSize[0]
661    scaley = 1000./pixelSize[1]
662    pixLimit = data['pixLimit']
663    cutoff = data['cutoff']
664    varyDict = data['varyList']
665    if varyDict['dist'] and varyDict['wave']:
666        print ('ERROR - you can not simultaneously calibrate distance and wavelength')
667        return False
668    if len(ring) < 5:
669        print ('ERROR - not enough inner ring points for ellipse')
670        return False
671       
672    #fit start points on inner ring
673    data['ellipses'] = []
674    data['rings'] = []
675    outE = FitEllipse(ring)
676    fmt  = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f'
677    fmt2 = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f, chi**2: %.3f, Np: %d'
678    if outE:
679        print (fmt%('start ellipse: ',outE[0][0],outE[0][1],outE[1],outE[2][0],outE[2][1]))
680        ellipse = outE
681    else:
682        return False
683       
684    #setup 360 points on that ring for "good" fit
685    data['ellipses'].append(ellipse[:]+('g',))
686    Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0]
687    if Ring:
688        ellipse = FitEllipse(Ring)
689        Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0]    #do again
690        ellipse = FitEllipse(Ring)
691    else:
692        print ('1st ring not sufficiently complete to proceed')
693        return False
694    if debug:
695        print (fmt2%('inner ring:    ',ellipse[0][0],ellipse[0][1],ellipse[1],
696            ellipse[2][0],ellipse[2][1],0.,len(Ring)))     #cent,phi,radii
697    data['ellipses'].append(ellipse[:]+('r',))
698    data['rings'].append(np.array(Ring))
699    G2plt.PlotImage(G2frame,newImage=True)
700   
701#setup for calibration
702    data['rings'] = []
703    if not data['calibrant']:
704        print ('no calibration material selected')
705        return True
706   
707    skip = data['calibskip']
708    dmin = data['calibdmin']
709#generate reflection set
710    Bravais,SGs,Cells = calFile.Calibrants[data['calibrant']][:3]
711    HKL = []
712    for bravais,sg,cell in zip(Bravais,SGs,Cells):
713        A = G2lat.cell2A(cell)
714        if sg:
715            SGData = G2spc.SpcGroup(sg)[1]
716            hkl = G2pwd.getHKLpeak(dmin,SGData,A)
717            print(hkl)
718            HKL += list(hkl)
719        else:
720            hkl = G2lat.GenHBravais(dmin,bravais,A)
721            HKL += list(hkl)
722    HKL = G2lat.sortHKLd(HKL,True,False)[skip:]
723#set up 1st ring
724    elcent,phi,radii = ellipse              #from fit of 1st ring
725    dsp = HKL[0][3]
726    print ('1st ring: try %.4f'%(dsp))
727    if varyDict['dist']:
728        wave = data['wavelength']
729        tth = 2.0*asind(wave/(2.*dsp))
730    else:   #varyDict['wave']!
731        dist = data['distance']
732        tth = npatan2d(radii[0],dist)
733        data['wavelength'] = wave =  2.0*dsp*sind(tth/2.0)
734    Ring0 = makeRing(dsp,ellipse,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
735    ttth = nptand(tth)
736    ctth = npcosd(tth)
737#1st estimate of tilt; assume ellipse - don't know sign though
738    if varyDict['tilt']:
739        tilt = npasind(np.sqrt(max(0.,1.-(radii[0]/radii[1])**2))*ctth)
740        if not tilt:
741            print ('WARNING - selected ring was fitted as a circle')
742            print (' - if detector was tilted we suggest you skip this ring - WARNING')
743    else:
744        tilt = data['tilt']
745#1st estimate of dist: sample to detector normal to plane
746    if varyDict['dist']:
747        data['distance'] = dist = radii[0]**2/(ttth*radii[1])
748    else:
749        dist = data['distance']
750    if varyDict['tilt']:
751#ellipse to cone axis (x-ray beam); 2 choices depending on sign of tilt
752        zdisp = radii[1]*ttth*tand(tilt)
753        zdism = radii[1]*ttth*tand(-tilt)
754#cone axis position; 2 choices. Which is right?     
755#NB: zdisp is || to major axis & phi is rotation of minor axis
756#thus shift from beam to ellipse center is [Z*sin(phi),-Z*cos(phi)]
757        centp = [elcent[0]+zdisp*sind(phi),elcent[1]-zdisp*cosd(phi)]
758        centm = [elcent[0]+zdism*sind(phi),elcent[1]-zdism*cosd(phi)]
759#check get same ellipse parms either way
760#now do next ring; estimate either way & do a FitDetector each way; best fit is correct one
761        fail = True
762        i2 = 1
763        while fail:
764            dsp = HKL[i2][3]
765            print ('2nd ring: try %.4f'%(dsp))
766            tth = 2.0*asind(wave/(2.*dsp))
767            ellipsep = GetEllipse2(tth,0.,dist,centp,tilt,phi)
768            print (fmt%('plus ellipse :',ellipsep[0][0],ellipsep[0][1],ellipsep[1],ellipsep[2][0],ellipsep[2][1]))
769            Ringp = makeRing(dsp,ellipsep,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
770            parmDict = {'dist':dist,'det-X':centp[0],'det-Y':centp[1],
771                'tilt':tilt,'phi':phi,'wave':wave,'dep':0.0}       
772            varyList = [item for item in varyDict if varyDict[item]]
773            if len(Ringp) > 10:
774                chip = FitDetector(np.array(Ring0+Ringp),varyList,parmDict,True)[0]
775                tiltp = parmDict['tilt']
776                phip = parmDict['phi']
777                centp = [parmDict['det-X'],parmDict['det-Y']]
778                fail = False
779            else:
780                chip = 1e6
781            ellipsem = GetEllipse2(tth,0.,dist,centm,-tilt,phi)
782            print (fmt%('minus ellipse:',ellipsem[0][0],ellipsem[0][1],ellipsem[1],ellipsem[2][0],ellipsem[2][1]))
783            Ringm = makeRing(dsp,ellipsem,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
784            if len(Ringm) > 10:
785                parmDict['tilt'] *= -1
786                chim = FitDetector(np.array(Ring0+Ringm),varyList,parmDict,True)[0]
787                tiltm = parmDict['tilt']
788                phim = parmDict['phi']
789                centm = [parmDict['det-X'],parmDict['det-Y']]
790                fail = False
791            else:
792                chim = 1e6
793            if fail:
794                i2 += 1
795        if chip < chim:
796            data['tilt'] = tiltp
797            data['center'] = centp
798            data['rotation'] = phip
799        else:
800            data['tilt'] = tiltm
801            data['center'] = centm
802            data['rotation'] = phim
803        data['ellipses'].append(ellipsep[:]+('b',))
804        data['rings'].append(np.array(Ringp))
805        data['ellipses'].append(ellipsem[:]+('r',))
806        data['rings'].append(np.array(Ringm))
807        G2plt.PlotImage(G2frame,newImage=True)
808    if data['DetDepth'] > 0.5:          #patch - redefine DetDepth
809        data['DetDepth'] /= data['distance']
810    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
811        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
812    varyList = [item for item in varyDict if varyDict[item]]
813    data['rings'] = []
814    data['ellipses'] = []
815    for i,H in enumerate(HKL):
816        dsp = H[3]
817        tth = 2.0*asind(wave/(2.*dsp))
818        if tth+abs(data['tilt']) > 90.:
819            print ('next line is a hyperbola - search stopped')
820            break
821        if debug:   print ('HKLD:'+str(H[:4])+'2-theta: %.4f'%(tth))
822        elcent,phi,radii = ellipse = GetEllipse(dsp,data)
823        data['ellipses'].append(copy.deepcopy(ellipse+('g',)))
824        if debug:   print (fmt%('predicted ellipse:',elcent[0],elcent[1],phi,radii[0],radii[1]))
825        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0]
826        if Ring:
827            data['rings'].append(np.array(Ring))
828            rings = np.concatenate((data['rings']),axis=0)
829            if i:
830                chisq = FitDetector(rings,varyList,parmDict,False)[0]
831                data['distance'] = parmDict['dist']
832                data['center'] = [parmDict['det-X'],parmDict['det-Y']]
833                data['rotation'] = parmDict['phi']
834                data['tilt'] = parmDict['tilt']
835                data['DetDepth'] = parmDict['dep']
836                data['chisq'] = chisq
837                elcent,phi,radii = ellipse = GetEllipse(dsp,data)
838                if debug:   print (fmt2%('fitted ellipse:   ',elcent[0],elcent[1],phi,radii[0],radii[1],chisq,len(rings)))
839            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
840#            G2plt.PlotImage(G2frame,newImage=True)
841        else:
842            if debug:   print ('insufficient number of points in this ellipse to fit')
843#            break
844    G2plt.PlotImage(G2frame,newImage=True)
845    fullSize = len(G2frame.ImageZ)/scalex
846    if 2*radii[1] < .9*fullSize:
847        print ('Are all usable rings (>25% visible) used? Try reducing Min ring I/Ib')
848    N = len(data['ellipses'])
849    if N > 2:
850        FitDetector(rings,varyList,parmDict)[0]
851        data['wavelength'] = parmDict['wave']
852        data['distance'] = parmDict['dist']
853        data['center'] = [parmDict['det-X'],parmDict['det-Y']]
854        data['rotation'] = parmDict['phi']
855        data['tilt'] = parmDict['tilt']
856        data['DetDepth'] = parmDict['dep']
857    for H in HKL[:N]:
858        ellipse = GetEllipse(H[3],data)
859        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
860    print ('calibration time = %.3f'%(time.time()-time0))
861    G2plt.PlotImage(G2frame,newImage=True)       
862    return True
863   
864def Make2ThetaAzimuthMap(data,iLim,jLim): #most expensive part of integration!
865    'Needs a doc string'
866    #transforms 2D image from x,y space to 2-theta,azimuth space based on detector orientation
867    pixelSize = data['pixelSize']
868    scalex = pixelSize[0]/1000.
869    scaley = pixelSize[1]/1000.
870   
871    tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5]         #bin centers not corners
872    tax = np.asfarray(tax*scalex,dtype=np.float32).flatten()
873    tay = np.asfarray(tay*scaley,dtype=np.float32).flatten()
874    nI = iLim[1]-iLim[0]
875    nJ = jLim[1]-jLim[0]
876    TA = np.array(GetTthAzmG(np.reshape(tax,(nI,nJ)),np.reshape(tay,(nI,nJ)),data))     #includes geom. corr. as dist**2/d0**2 - most expensive step
877    TA[1] = np.where(TA[1]<0,TA[1]+360,TA[1])
878    return TA           #2-theta, azimuth & geom. corr. arrays
879
880def MakeMaskMap(data,masks,iLim,jLim,tamp):
881    import polymask as pm
882    pixelSize = data['pixelSize']
883    scalex = pixelSize[0]/1000.
884    scaley = pixelSize[1]/1000.
885   
886    tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5]         #bin centers not corners
887    tax = np.asfarray(tax*scalex,dtype=np.float32).flatten()
888    tay = np.asfarray(tay*scaley,dtype=np.float32).flatten()
889    nI = iLim[1]-iLim[0]
890    nJ = jLim[1]-jLim[0]
891    #make position masks here
892    frame = masks['Frames']
893    tam = ma.make_mask_none((nI*nJ))
894    if frame:
895        tam = ma.mask_or(tam,ma.make_mask(pm.polymask(nI*nJ,tax,
896            tay,len(frame),frame,tamp)[:nI*nJ])^True)
897    polygons = masks['Polygons']
898    for polygon in polygons:
899        if polygon:
900            tam = ma.mask_or(tam,ma.make_mask(pm.polymask(nI*nJ,tax,
901                tay,len(polygon),polygon,tamp)[:nI*nJ]))
902    if True:
903        for X,Y,rsq in masks['Points'].T:
904            tam = ma.mask_or(tam,ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,rsq)))
905    if tam.shape: tam = np.reshape(tam,(nI,nJ))
906    return tam           #position mask
907
908def Fill2ThetaAzimuthMap(masks,TA,tam,image):
909    'Needs a doc string'
910    Zlim = masks['Thresholds'][1]
911    rings = masks['Rings']
912    arcs = masks['Arcs']
913    TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0]),ma.getdata(TA[2])))    #azimuth, 2-theta, dist
914    tax,tay,tad = np.dsplit(TA,3)    #azimuth, 2-theta, dist**2/d0**2
915    for tth,thick in rings:
916        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.)))
917    for tth,azm,thick in arcs:
918        tamt = ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.))
919        tama = ma.getmask(ma.masked_inside(tax.flatten(),azm[0],azm[1]))
920        tam = ma.mask_or(tam.flatten(),tamt*tama)
921    taz = ma.masked_outside(image.flatten(),int(Zlim[0]),Zlim[1])
922    tabs = np.ones_like(taz)
923    tam = ma.mask_or(tam.flatten(),ma.getmask(taz))
924    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))   #azimuth
925    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))   #2-theta
926    taz = ma.compressed(ma.array(taz.flatten(),mask=tam))   #intensity
927    tad = ma.compressed(ma.array(tad.flatten(),mask=tam))   #dist**2/d0**2
928    tabs = ma.compressed(ma.array(tabs.flatten(),mask=tam)) #ones - later used for absorption corr.
929    return tax,tay,taz,tad,tabs
930
931def MakeUseTA(data,blkSize=128):
932    Nx,Ny = data['size']
933    nXBlks = (Nx-1)//blkSize+1
934    nYBlks = (Ny-1)//blkSize+1
935    useTA = []
936    for iBlk in range(nYBlks):
937        iBeg = iBlk*blkSize
938        iFin = min(iBeg+blkSize,Ny)
939        useTAj = []
940        for jBlk in range(nXBlks):
941            jBeg = jBlk*blkSize
942            jFin = min(jBeg+blkSize,Nx)
943            TA = Make2ThetaAzimuthMap(data,(iBeg,iFin),(jBeg,jFin))          #2-theta & azimuth arrays & create position mask
944            useTAj.append(TA)
945        useTA.append(useTAj)
946    return useTA
947
948def MakeUseMask(data,masks,blkSize=128):
949    Masks = copy.deepcopy(masks)
950    Masks['Points'] = np.array(Masks['Points']).T           #get spots as X,Y,R arrays
951    if np.any(masks['Points']):
952        Masks['Points'][2] = np.square(Masks['Points'][2]/2.)
953    Nx,Ny = data['size']
954    nXBlks = (Nx-1)//blkSize+1
955    nYBlks = (Ny-1)//blkSize+1
956    useMask = []
957    tamp = ma.make_mask_none((1024*1024))       #NB: this array size used in the fortran histogram2d
958    for iBlk in range(nYBlks):
959        iBeg = iBlk*blkSize
960        iFin = min(iBeg+blkSize,Ny)
961        useMaskj = []
962        for jBlk in range(nXBlks):
963            jBeg = jBlk*blkSize
964            jFin = min(jBeg+blkSize,Nx)
965            mask = MakeMaskMap(data,Masks,(iBeg,iFin),(jBeg,jFin),tamp)          #2-theta & azimuth arrays & create position mask
966            useMaskj.append(mask)
967        useMask.append(useMaskj)
968    return useMask
969
970def ImageIntegrate(image,data,masks,blkSize=128,returnN=False,useTA=None,useMask=None):
971    'Integrate an image; called from OnIntegrateAll and OnIntegrate in G2imgGUI'    #for q, log(q) bins need data['binType']
972    import histogram2d as h2d
973    print ('Begin image integration; image range: %d %d'%(np.min(image),np.max(image)))
974    CancelPressed = False
975    LUtth = np.array(data['IOtth'])
976    LRazm = np.array(data['LRazimuth'],dtype=np.float64)
977    numAzms = data['outAzimuths']
978    numChans = (data['outChannels']//4)*4
979    Dazm = (LRazm[1]-LRazm[0])/numAzms
980    if '2-theta' in data.get('binType','2-theta'):
981        lutth = LUtth               
982    elif 'log(q)' in data['binType']:
983        lutth = np.log(4.*np.pi*npsind(LUtth/2.)/data['wavelength'])
984    elif 'q' == data['binType'].lower():
985        lutth = 4.*np.pi*npsind(LUtth/2.)/data['wavelength']
986    dtth = (lutth[1]-lutth[0])/numChans
987    muT = data.get('SampleAbs',[0.0,''])[0]
988    if data['DetDepth'] > 0.5:          #patch - redefine DetDepth
989        data['DetDepth'] /= data['distance']
990    if 'SASD' in data['type']:
991        muT = -np.log(muT)/2.       #Transmission to 1/2 thickness muT
992    Masks = copy.deepcopy(masks)
993    Masks['Points'] = np.array(Masks['Points']).T           #get spots as X,Y,R arrays
994    if np.any(masks['Points']):
995        Masks['Points'][2] = np.square(Masks['Points'][2]/2.)
996    NST = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
997    H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
998    H2 = np.linspace(lutth[0],lutth[1],numChans+1)
999    Nx,Ny = data['size']
1000    nXBlks = (Nx-1)//blkSize+1
1001    nYBlks = (Ny-1)//blkSize+1
1002    tbeg = time.time()
1003    times = [0,0,0,0,0]
1004    tamp = ma.make_mask_none((1024*1024))       #NB: this array size used in the fortran histogram2d
1005    for iBlk in range(nYBlks):
1006        iBeg = iBlk*blkSize
1007        iFin = min(iBeg+blkSize,Ny)
1008        for jBlk in range(nXBlks):
1009            jBeg = jBlk*blkSize
1010            jFin = min(jBeg+blkSize,Nx)
1011            # next is most expensive step!
1012            t0 = time.time()
1013            if useTA:
1014                TA = useTA[iBlk][jBlk]
1015            else:
1016                TA = Make2ThetaAzimuthMap(data,(iBeg,iFin),(jBeg,jFin))           #2-theta & azimuth arrays & create position mask
1017            times[1] += time.time()-t0
1018            t0 = time.time()
1019            if useMask:
1020                tam = useMask[iBlk][jBlk]
1021            else:
1022                tam = MakeMaskMap(data,Masks,(iBeg,iFin),(jBeg,jFin),tamp)
1023            Block = image[iBeg:iFin,jBeg:jFin]
1024            tax,tay,taz,tad,tabs = Fill2ThetaAzimuthMap(Masks,TA,tam,Block)    #and apply masks
1025            times[0] += time.time()-t0
1026            t0 = time.time()
1027            tax = np.where(tax > LRazm[1],tax-360.,tax)                 #put azm inside limits if possible
1028            tax = np.where(tax < LRazm[0],tax+360.,tax)
1029            if data.get('SampleAbs',[0.0,''])[1]:
1030                if 'Cylind' in data['SampleShape']:
1031                    muR = muT*(1.+npsind(tax)**2/2.)/(npcosd(tay))      #adjust for additional thickness off sample normal
1032                    tabs = G2pwd.Absorb(data['SampleShape'],muR,tay)
1033                elif 'Fixed' in data['SampleShape']:    #assumes flat plate sample normal to beam
1034                    tabs = G2pwd.Absorb('Fixed',muT,tay)
1035            if 'log(q)' in data.get('binType',''):
1036                tay = np.log(4.*np.pi*npsind(tay/2.)/data['wavelength'])
1037            elif 'q' == data.get('binType','').lower():
1038                tay = 4.*np.pi*npsind(tay/2.)/data['wavelength']
1039            times[2] += time.time()-t0
1040            t0 = time.time()
1041            taz = np.array((taz*tad/tabs),dtype='float32')
1042            if any([tax.shape[0],tay.shape[0],taz.shape[0]]):
1043                NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,
1044                    numAzms,numChans,LRazm,lutth,Dazm,dtth,NST,H0)
1045            times[3] += time.time()-t0
1046#            print('done block %d %d %d %d %d %d %d %d'%(iBlk,iBeg,iFin,jBlk,jBeg,jFin,np.min(Block),np.max(Block)))
1047    print('End integration loops')
1048    t0 = time.time()
1049#    H2 = np.array([tth for tth in np.linspace(lutth[0],lutth[1],numChans+1)])
1050#    NST = np.array(NST,dtype=np.float32)
1051    #prepare masked arrays of bins with pixels for interpolation setup
1052    H2msk = [ma.array(H2[:-1],mask=np.logical_not(nst)) for nst in NST]
1053    H0msk = [ma.array(np.divide(h0,nst),mask=np.logical_not(nst)) for nst,h0 in zip(NST,H0)]
1054    #make linear interpolators; outside limits give NaN
1055    H0int = [scint.interp1d(h2msk.compressed(),h0msk.compressed(),bounds_error=False) for h0msk,h2msk in zip(H0msk,H2msk)]
1056    #do interpolation on all points - fills in the empty bins; leaves others the same
1057    H0 = np.array([h0int(H2[:-1]) for h0int in H0int])
1058    H0 = np.nan_to_num(H0)
1059    if 'log(q)' in data.get('binType',''):
1060        H2 = 2.*npasind(np.exp(H2)*data['wavelength']/(4.*np.pi))
1061    elif 'q' == data.get('binType','').lower():
1062        H2 = 2.*npasind(H2*data['wavelength']/(4.*np.pi))
1063    if Dazm:       
1064        H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)])
1065    else:
1066        H1 = LRazm
1067    H0 /= npcosd(H2[:-1])           #**2? I don't think so, **1 is right for powders
1068    if 'SASD' in data['type']:
1069        H0 /= npcosd(H2[:-1])           #one more for small angle scattering data?
1070    if data['Oblique'][1]:
1071        H0 /= G2pwd.Oblique(data['Oblique'][0],H2[:-1])
1072    if 'SASD' in data['type'] and data['PolaVal'][1]:
1073        #NB: in G2pwd.Polarization azm is defined from plane of polarization, not image x axis!
1074        H0 /= np.array([G2pwd.Polarization(data['PolaVal'][0],H2[:-1],Azm=azm-90.)[0] for azm in (H1[:-1]+np.diff(H1)/2.)])
1075    times[4] += time.time()-t0
1076    print ('Step times: \n apply masks  %8.3fs xy->th,azm   %8.3fs fill map     %8.3fs \
1077        \n binning      %8.3fs cleanup      %8.3fs'%(times[0],times[1],times[2],times[3],times[4]))
1078    print ("Elapsed time:","%8.3fs"%(time.time()-tbeg))
1079    print ('Integration complete')
1080    if returnN:     #As requested by Steven Weigand
1081        return H0,H1,H2,NST,CancelPressed
1082    else:
1083        return H0,H1,H2,CancelPressed
1084   
1085def MakeStrStaRing(ring,Image,Controls):
1086    ellipse = GetEllipse(ring['Dset'],Controls)
1087    pixSize = Controls['pixelSize']
1088    scalex = 1000./pixSize[0]
1089    scaley = 1000./pixSize[1]
1090    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
1091    if len(Ring):
1092        ring['ImxyObs'] = copy.copy(Ring[:2])
1093        TA = GetTthAzm(Ring[0],Ring[1],Controls)       #convert x,y to tth,azm
1094        TA[0] = Controls['wavelength']/(2.*npsind(TA[0]/2.))      #convert 2th to d
1095        ring['ImtaObs'] = TA
1096        ring['ImtaCalc'] = np.zeros_like(ring['ImtaObs'])
1097        Ring[0] = TA[0]
1098        Ring[1] = TA[1]
1099        return Ring,ring
1100    else:
1101        ring['ImxyObs'] = [[],[]]
1102        ring['ImtaObs'] = [[],[]]
1103        ring['ImtaCalc'] = [[],[]]
1104        return [],[]    #bad ring; no points found
1105   
1106def FitStrSta(Image,StrSta,Controls):
1107    'Needs a doc string'
1108   
1109    StaControls = copy.deepcopy(Controls)
1110    phi = StrSta['Sample phi']
1111    wave = Controls['wavelength']
1112    pixelSize = Controls['pixelSize']
1113    scalex = 1000./pixelSize[0]
1114    scaley = 1000./pixelSize[1]
1115    StaType = StrSta['Type']
1116    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
1117
1118    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
1119        dset = ring['Dset']
1120        Ring,R = MakeStrStaRing(ring,Image,StaControls)
1121        if len(Ring):
1122            ring.update(R)
1123            p0 = ring['Emat']
1124            val,esd,covMat = FitStrain(Ring,p0,dset,wave,phi,StaType)
1125            ring['Emat'] = val
1126            ring['Esig'] = esd
1127            ellipse = FitEllipse(R['ImxyObs'].T)
1128            ringxy,ringazm = makeRing(ring['Dcalc'],ellipse,0,0.,scalex,scaley,Image)
1129            ring['ImxyCalc'] = np.array(ringxy).T[:2]
1130            ringint = np.array([float(Image[int(x*scalex),int(y*scaley)]) for y,x in np.array(ringxy)[:,:2]])
1131            ringint /= np.mean(ringint)
1132            ring['Ivar'] = np.var(ringint)
1133            ring['covMat'] = covMat
1134            print ('Variance in normalized ring intensity: %.3f'%(ring['Ivar']))
1135    CalcStrSta(StrSta,Controls)
1136   
1137def IntStrSta(Image,StrSta,Controls):
1138    StaControls = copy.deepcopy(Controls)
1139    pixelSize = Controls['pixelSize']
1140    scalex = 1000./pixelSize[0]
1141    scaley = 1000./pixelSize[1]
1142    phi = StrSta['Sample phi']
1143    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
1144    RingsAI = []
1145    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
1146        Ring,R = MakeStrStaRing(ring,Image,StaControls)
1147        if len(Ring):
1148            ellipse = FitEllipse(R['ImxyObs'].T)
1149            ringxy,ringazm = makeRing(ring['Dcalc'],ellipse,0,0.,scalex,scaley,Image,5)
1150            ring['ImxyCalc'] = np.array(ringxy).T[:2]
1151            ringint = np.array([float(Image[int(x*scalex),int(y*scaley)]) for y,x in np.array(ringxy)[:,:2]])
1152            ringint /= np.mean(ringint)
1153            print (' %s %.3f %s %.3f %s %d'%('d-spacing',ring['Dcalc'],'sig(MRD):',np.sqrt(np.var(ringint)),'# points:',len(ringint)))
1154            RingsAI.append(np.array(zip(ringazm,ringint)).T)
1155    return RingsAI
1156   
1157def CalcStrSta(StrSta,Controls):
1158
1159    wave = Controls['wavelength']
1160    phi = StrSta['Sample phi']
1161    StaType = StrSta['Type']
1162    for ring in StrSta['d-zero']:
1163        Eij = ring['Emat']
1164        E = [[Eij[0],Eij[1],0],[Eij[1],Eij[2],0],[0,0,0]]
1165        th,azm = ring['ImtaObs']
1166        th0 = np.ones_like(azm)*npasind(wave/(2.*ring['Dset']))
1167        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th0).T/1.e6,axis=2),axis=1)
1168        if StaType == 'True':
1169            ring['ImtaCalc'] = np.array([np.exp(V)*ring['Dset'],azm])
1170        else:
1171            ring['ImtaCalc'] = np.array([(V+1.)*ring['Dset'],azm])
1172        dmin = np.min(ring['ImtaCalc'][0])
1173        dmax = np.max(ring['ImtaCalc'][0])
1174        if ring.get('fixDset',True):
1175            if abs(Eij[0]) < abs(Eij[2]):         #tension
1176                ring['Dcalc'] = dmin+(dmax-dmin)/4.
1177            else:                       #compression
1178                ring['Dcalc'] = dmin+3.*(dmax-dmin)/4.
1179        else:
1180            ring['Dcalc'] = np.mean(ring['ImtaCalc'][0])
1181
1182def calcFij(omg,phi,azm,th):
1183    '''    Uses parameters as defined by Bob He & Kingsley Smith, Adv. in X-Ray Anal. 41, 501 (1997)
1184
1185    :param omg: his omega = sample omega rotation; 0 when incident beam || sample surface,
1186        90 when perp. to sample surface
1187    :param phi: his phi = sample phi rotation; usually = 0, axis rotates with omg.
1188    :param azm: his chi = azimuth around incident beam
1189    :param th:  his theta = theta
1190    '''
1191    a = npsind(th)*npcosd(omg)+npsind(azm)*npcosd(th)*npsind(omg)
1192    b = -npcosd(azm)*npcosd(th)
1193    c = npsind(th)*npsind(omg)-npsind(azm)*npcosd(th)*npcosd(omg)
1194    d = a*npsind(phi)+b*npcosd(phi)
1195    e = a*npcosd(phi)-b*npsind(phi)
1196    Fij = np.array([
1197        [d**2,d*e,c*d],
1198        [d*e,e**2,c*e],
1199        [c*d,c*e,c**2]])
1200    return -Fij
1201
1202def FitStrain(rings,p0,dset,wave,phi,StaType):
1203    'Needs a doc string'
1204    def StrainPrint(ValSig,dset):
1205        print ('Strain tensor for Dset: %.6f'%(dset))
1206        ptlbls = 'names :'
1207        ptstr =  'values:'
1208        sigstr = 'esds  :'
1209        for name,fmt,value,sig in ValSig:
1210            ptlbls += "%s" % (name.rjust(12))
1211            ptstr += fmt % (value)
1212            if sig:
1213                sigstr += fmt % (sig)
1214            else:
1215                sigstr += 12*' '
1216        print (ptlbls)
1217        print (ptstr)
1218        print (sigstr)
1219       
1220    def strainCalc(p,xyd,dset,wave,phi,StaType):
1221        E = np.array([[p[0],p[1],0],[p[1],p[2],0],[0,0,0]])
1222        dspo,azm,dsp = xyd
1223        th = npasind(wave/(2.0*dspo))
1224        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th).T/1.e6,axis=2),axis=1)
1225        if StaType == 'True':
1226            dspc = dset*np.exp(V)
1227        else:
1228            dspc = dset*(V+1.)
1229        return dspo-dspc
1230       
1231    names = ['e11','e12','e22']
1232    fmt = ['%12.2f','%12.2f','%12.2f']
1233    result = leastsq(strainCalc,p0,args=(rings,dset,wave,phi,StaType),full_output=True)
1234    vals = list(result[0])
1235    chisq = np.sum(result[2]['fvec']**2)/(rings.shape[1]-3)     #reduced chi^2 = M/(Nobs-Nvar)
1236    covM = result[1]
1237    covMatrix = covM*chisq
1238    sig = list(np.sqrt(chisq*np.diag(result[1])))
1239    ValSig = zip(names,fmt,vals,sig)
1240    StrainPrint(ValSig,dset)
1241    return vals,sig,covMatrix
1242
1243def FitImageSpots(Image,ImMax,ind,pixSize,nxy):
1244   
1245    def calcMean(nxy,pixSize,img):
1246        gdx,gdy = np.mgrid[0:nxy,0:nxy]
1247        gdx = ma.array((gdx-nxy//2)*pixSize[0]/1000.,mask=~ma.getmaskarray(ImBox))
1248        gdy = ma.array((gdy-nxy//2)*pixSize[1]/1000.,mask=~ma.getmaskarray(ImBox))
1249        posx = ma.sum(gdx)/ma.count(gdx)
1250        posy = ma.sum(gdy)/ma.count(gdy)
1251        return posx,posy
1252   
1253    def calcPeak(values,nxy,pixSize,img):
1254        back,mag,px,py,sig = values
1255        peak = np.zeros([nxy,nxy])+back
1256        nor = 1./(2.*np.pi*sig**2)
1257        gdx,gdy = np.mgrid[0:nxy,0:nxy]
1258        gdx = (gdx-nxy//2)*pixSize[0]/1000.
1259        gdy = (gdy-nxy//2)*pixSize[1]/1000.
1260        arg = (gdx-px)**2+(gdy-py)**2       
1261        peak += mag*nor*np.exp(-arg/(2.*sig**2))
1262        return ma.compressed(img-peak)/np.sqrt(ma.compressed(img))
1263   
1264    def calc2Peak(values,nxy,pixSize,img):
1265        back,mag,px,py,sigx,sigy,rho = values
1266        peak = np.zeros([nxy,nxy])+back
1267        nor = 1./(2.*np.pi*sigx*sigy*np.sqrt(1.-rho**2))
1268        gdx,gdy = np.mgrid[0:nxy,0:nxy]
1269        gdx = (gdx-nxy//2)*pixSize[0]/1000.
1270        gdy = (gdy-nxy//2)*pixSize[1]/1000.
1271        argnor = -1./(2.*(1.-rho**2))
1272        arg = (gdx-px)**2/sigx**2+(gdy-py)**2/sigy**2-2.*rho*(gdx-px)*(gdy-py)/(sigx*sigy)       
1273        peak += mag*nor*np.exp(argnor*arg)
1274        return ma.compressed(img-peak)/np.sqrt(ma.compressed(img))       
1275   
1276    nxy2 = nxy//2
1277    ImBox = Image[ind[1]-nxy2:ind[1]+nxy2+1,ind[0]-nxy2:ind[0]+nxy2+1]
1278    back = np.min(ImBox)
1279    mag = np.sum(ImBox-back)
1280    vals = [back,mag,0.,0.,0.2,0.2,0.]
1281    ImBox = ma.array(ImBox,dtype=float,mask=ImBox>0.75*ImMax)
1282    px = (ind[0]+.5)*pixSize[0]/1000.
1283    py = (ind[1]+.5)*pixSize[1]/1000.
1284    if ma.any(ma.getmaskarray(ImBox)):
1285        vals = calcMean(nxy,pixSize,ImBox)
1286        posx,posy = [px+vals[0],py+vals[1]]
1287        return [posx,posy,6.]
1288    else:
1289        result = leastsq(calc2Peak,vals,args=(nxy,pixSize,ImBox),full_output=True)
1290        vals = result[0]
1291        ratio = vals[4]/vals[5]
1292        if 0.5 < ratio < 2.0 and vals[2] < 2. and vals[3] < 2.:
1293            posx,posy = [px+vals[2],py+vals[3]]
1294            return [posx,posy,min(6.*vals[4],6.)]
1295        else:
1296            return None
1297   
1298def AutoSpotMasks(Image,Masks,Controls):
1299   
1300    print ('auto spot search')
1301    nxy = 15
1302    rollImage = lambda rho,roll: np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1)
1303    pixelSize = Controls['pixelSize']
1304    spotMask = ma.array(Image,mask=(Image<np.mean(Image)))
1305    indices = (-1,0,1)
1306    rolls = np.array([[ix,iy] for ix in indices for iy in indices])
1307    time0 = time.time()
1308    for roll in rolls:
1309        if np.any(roll):        #avoid [0,0]
1310            spotMask = ma.array(spotMask,mask=(spotMask-rollImage(Image,roll)<0.),dtype=float)
1311    mags = spotMask[spotMask.nonzero()]
1312    indx = np.transpose(spotMask.nonzero())
1313    size1 = mags.shape[0]
1314    magind = [[indx[0][0],indx[0][1],mags[0]],]
1315    for ind,mag in list(zip(indx,mags))[1:]:        #remove duplicates
1316#            ind[0],ind[1],I,J = ImageLocalMax(Image,nxy,ind[0],ind[1])
1317        if (magind[-1][0]-ind[0])**2+(magind[-1][1]-ind[1])**2 > 16:
1318            magind.append([ind[0],ind[1],Image[ind[0],ind[1]]]) 
1319    magind = np.array(magind).T
1320    indx = np.array(magind[0:2],dtype=np.int32)
1321    mags = magind[2]
1322    size2 = mags.shape[0]
1323    print ('Initial search done: %d -->%d %.2fs'%(size1,size2,time.time()-time0))
1324    nx,ny = Image.shape
1325    ImMax = np.max(Image)
1326    peaks = []
1327    nxy2 = nxy//2
1328    mult = 0.001
1329    num = 1e6
1330    while num>500:
1331        mult += .0001           
1332        minM = mult*np.max(mags)
1333        num = ma.count(ma.array(mags,mask=mags<=minM))
1334        print('try',mult,minM,num)
1335    minM = mult*np.max(mags)
1336    print ('Find biggest spots:',mult,num,minM)
1337    for i,mag in enumerate(mags):
1338        if mag > minM:
1339            if (nxy2 < indx[0][i] < nx-nxy2-1) and (nxy2 < indx[1][i] < ny-nxy2-1):
1340#                    print ('try:%d %d %d %.2f'%(i,indx[0][i],indx[1][i],mags[i]))
1341                peak = FitImageSpots(Image,ImMax,[indx[1][i],indx[0][i]],pixelSize,nxy)
1342                if peak and not any(np.isnan(np.array(peak))):
1343                    peaks.append(peak)
1344#                    print (' Spot found: %s'%str(peak))
1345    peaks = G2mth.sortArray(G2mth.sortArray(peaks,1),0)
1346    Peaks = [peaks[0],]
1347    for peak in peaks[1:]:
1348        if GetDsp(peak[0],peak[1],Controls) >= 1.:      #toss non-diamond low angle spots
1349            continue
1350        if (peak[0]-Peaks[-1][0])**2+(peak[1]-Peaks[-1][1])**2 > peak[2]*Peaks[-1][2] :
1351            Peaks.append(peak)
1352#            print (' Spot found: %s'%str(peak))
1353    print ('Spots found: %d time %.2fs'%(len(Peaks),time.time()-time0))
1354    Masks['Points'] = Peaks
1355    return None
Note: See TracBrowser for help on using the repository browser.