source: trunk/GSASIIimage.py @ 3653

Last change on this file since 3653 was 3653, checked in by vondreele, 3 years ago

fix image calibration bug - a python 3.6 issue

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