source: trunk/GSASIIimage.py @ 3857

Last change on this file since 3857 was 3857, checked in by toby, 4 years ago

Fix scriptable to work w/o wx & mpl; reorg of docs

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