source: trunk/GSASIIimage.py @ 4021

Last change on this file since 4021 was 4021, checked in by toby, 2 years ago

implement filter for screen messages; start to replace print() with G2Print(); scripting changes for PDF

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