source: trunk/GSASIIimage.py @ 3362

Last change on this file since 3362 was 3362, checked in by vondreele, 5 years ago

fix azimuth offset problems with linescan
add sqrt & log (scanline) options

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 53.3 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASII image calculations: ellipse fitting & image integration       
3########### SVN repository information ###################
4# $Date: 2018-04-27 21:23:53 +0000 (Fri, 27 Apr 2018) $
5# $Author: vondreele $
6# $Revision: 3362 $
7# $URL: trunk/GSASIIimage.py $
8# $Id: GSASIIimage.py 3362 2018-04-27 21:23:53Z 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: 3362 $")
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 copy
631    import ImageCalibrants as calFile
632    print ('Image calibration:')
633    time0 = time.time()
634    ring = data['ring']
635    pixelSize = data['pixelSize']
636    scalex = 1000./pixelSize[0]
637    scaley = 1000./pixelSize[1]
638    pixLimit = data['pixLimit']
639    cutoff = data['cutoff']
640    varyDict = data['varyList']
641    if varyDict['dist'] and varyDict['wave']:
642        print ('ERROR - you can not simultaneously calibrate distance and wavelength')
643        return False
644    if len(ring) < 5:
645        print ('ERROR - not enough inner ring points for ellipse')
646        return False
647       
648    #fit start points on inner ring
649    data['ellipses'] = []
650    data['rings'] = []
651    outE = FitEllipse(ring)
652    fmt  = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f'
653    fmt2 = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f, chi**2: %.3f, Np: %d'
654    if outE:
655        print (fmt%('start ellipse: ',outE[0][0],outE[0][1],outE[1],outE[2][0],outE[2][1]))
656        ellipse = outE
657    else:
658        return False
659       
660    #setup 360 points on that ring for "good" fit
661    data['ellipses'].append(ellipse[:]+('g',))
662    Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0]
663    if Ring:
664        ellipse = FitEllipse(Ring)
665        Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0]    #do again
666        ellipse = FitEllipse(Ring)
667    else:
668        print ('1st ring not sufficiently complete to proceed')
669        return False
670    if debug:
671        print (fmt2%('inner ring:    ',ellipse[0][0],ellipse[0][1],ellipse[1],
672            ellipse[2][0],ellipse[2][1],0.,len(Ring)))     #cent,phi,radii
673    data['ellipses'].append(ellipse[:]+('r',))
674    data['rings'].append(np.array(Ring))
675    G2plt.PlotImage(G2frame,newImage=True)
676   
677#setup for calibration
678    data['rings'] = []
679    if not data['calibrant']:
680        print ('no calibration material selected')
681        return True
682   
683    skip = data['calibskip']
684    dmin = data['calibdmin']
685#generate reflection set
686    Bravais,SGs,Cells = calFile.Calibrants[data['calibrant']][:3]
687    HKL = []
688    for bravais,sg,cell in zip(Bravais,SGs,Cells):
689        A = G2lat.cell2A(cell)
690        if sg:
691            SGData = G2spc.SpcGroup(sg)[1]
692            hkl = G2pwd.getHKLpeak(dmin,SGData,A)
693            HKL += hkl
694        else:
695            hkl = G2lat.GenHBravais(dmin,bravais,A)
696            HKL += 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    Nx,Ny = data['size']
974    nXBlks = (Nx-1)//blkSize+1
975    nYBlks = (Ny-1)//blkSize+1
976    tbeg = time.time()
977    times = [0,0,0,0,0]
978    tamp = ma.make_mask_none((1024*1024))       #NB: this array size used in the fortran histogram2d
979    for iBlk in range(nYBlks):
980        iBeg = iBlk*blkSize
981        iFin = min(iBeg+blkSize,Ny)
982        for jBlk in range(nXBlks):
983            jBeg = jBlk*blkSize
984            jFin = min(jBeg+blkSize,Nx)
985            # next is most expensive step!
986            t0 = time.time()
987            if useTA:
988                TA = useTA[iBlk][jBlk]
989            else:
990                TA = Make2ThetaAzimuthMap(data,(iBeg,iFin),(jBeg,jFin))           #2-theta & azimuth arrays & create position mask
991            times[1] += time.time()-t0
992            t0 = time.time()
993            if useMask:
994                tam = useMask[iBlk][jBlk]
995            else:
996                tam = MakeMaskMap(data,Masks,(iBeg,iFin),(jBeg,jFin),tamp)
997            Block = image[iBeg:iFin,jBeg:jFin]
998            tax,tay,taz,tad,tabs = Fill2ThetaAzimuthMap(Masks,TA,tam,Block)    #and apply masks
999            times[0] += time.time()-t0
1000            t0 = time.time()
1001            tax = np.where(tax > LRazm[1],tax-360.,tax)                 #put azm inside limits if possible
1002            tax = np.where(tax < LRazm[0],tax+360.,tax)
1003            if data.get('SampleAbs',[0.0,''])[1]:
1004                if 'Cylind' in data['SampleShape']:
1005                    muR = muT*(1.+npsind(tax)**2/2.)/(npcosd(tay))      #adjust for additional thickness off sample normal
1006                    tabs = G2pwd.Absorb(data['SampleShape'],muR,tay)
1007                elif 'Fixed' in data['SampleShape']:    #assumes flat plate sample normal to beam
1008                    tabs = G2pwd.Absorb('Fixed',muT,tay)
1009            if 'log(q)' in data.get('binType',''):
1010                tay = np.log(4.*np.pi*npsind(tay/2.)/data['wavelength'])
1011            elif 'q' == data.get('binType','').lower():
1012                tay = 4.*np.pi*npsind(tay/2.)/data['wavelength']
1013            times[2] += time.time()-t0
1014            t0 = time.time()
1015            taz = np.array((taz*tad/tabs),dtype='float32')
1016            if any([tax.shape[0],tay.shape[0],taz.shape[0]]):
1017                NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,
1018                    numAzms,numChans,LRazm,lutth,Dazm,dtth,NST,H0)
1019            times[3] += time.time()-t0
1020            del tax; del tay; del taz; del tad; del tabs
1021    t0 = time.time()
1022    H2 = np.array([tth for tth in np.linspace(lutth[0],lutth[1],numChans+1)])
1023    NST = np.array(NST,dtype=np.float)
1024    #prepare masked arrays of bins with pixels for interpolation setup
1025    H2msk = [ma.array(H2[:-1],mask=np.logical_not(nst)) for nst in NST]
1026    H0msk = [ma.array(np.divide(h0,nst),mask=np.logical_not(nst)) for nst,h0 in zip(NST,H0)]
1027    #make linear interpolators; outside limits give NaN
1028    H0int = [scint.interp1d(h2msk.compressed(),h0msk.compressed(),bounds_error=False) for h0msk,h2msk in zip(H0msk,H2msk)]
1029    #do interpolation on all points - fills in the empty bins; leaves others the same
1030    H0 = np.array([h0int(H2[:-1]) for h0int in H0int])
1031    H0 = np.nan_to_num(H0)
1032    if 'log(q)' in data.get('binType',''):
1033        H2 = 2.*npasind(np.exp(H2)*data['wavelength']/(4.*np.pi))
1034    elif 'q' == data.get('binType','').lower():
1035        H2 = 2.*npasind(H2*data['wavelength']/(4.*np.pi))
1036    if Dazm:       
1037        H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)])
1038    else:
1039        H1 = LRazm
1040    H0 /= npcosd(H2[:-1])           #**2? I don't think so, **1 is right for powders
1041    if 'SASD' in data['type']:
1042        H0 /= npcosd(H2[:-1])           #one more for small angle scattering data?
1043    if data['Oblique'][1]:
1044        H0 /= G2pwd.Oblique(data['Oblique'][0],H2[:-1])
1045    if 'SASD' in data['type'] and data['PolaVal'][1]:
1046        #NB: in G2pwd.Polarization azm is defined from plane of polarization, not image x axis!
1047        H0 /= np.array([G2pwd.Polarization(data['PolaVal'][0],H2[:-1],Azm=azm-90.)[0] for azm in (H1[:-1]+np.diff(H1)/2.)])
1048    times[4] += time.time()-t0
1049    print ('Step times: \n apply masks  %8.3fs xy->th,azm   %8.3fs fill map     %8.3fs \
1050        \n binning      %8.3fs cleanup      %8.3fs'%(times[0],times[1],times[2],times[3],times[4]))
1051    print ("Elapsed time:","%8.3fs"%(time.time()-tbeg))
1052    print ('Integration complete')
1053    if returnN:     #As requested by Steven Weigand
1054        return H0,H1,H2,NST,CancelPressed
1055    else:
1056        return H0,H1,H2,CancelPressed
1057   
1058def MakeStrStaRing(ring,Image,Controls):
1059    ellipse = GetEllipse(ring['Dset'],Controls)
1060    pixSize = Controls['pixelSize']
1061    scalex = 1000./pixSize[0]
1062    scaley = 1000./pixSize[1]
1063    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
1064    if len(Ring):
1065        ring['ImxyObs'] = copy.copy(Ring[:2])
1066        TA = GetTthAzm(Ring[0],Ring[1],Controls)       #convert x,y to tth,azm
1067        TA[0] = Controls['wavelength']/(2.*npsind(TA[0]/2.))      #convert 2th to d
1068        ring['ImtaObs'] = TA
1069        ring['ImtaCalc'] = np.zeros_like(ring['ImtaObs'])
1070        Ring[0] = TA[0]
1071        Ring[1] = TA[1]
1072        return Ring,ring
1073    else:
1074        ring['ImxyObs'] = [[],[]]
1075        ring['ImtaObs'] = [[],[]]
1076        ring['ImtaCalc'] = [[],[]]
1077        return [],[]    #bad ring; no points found
1078   
1079def FitStrSta(Image,StrSta,Controls):
1080    'Needs a doc string'
1081   
1082    StaControls = copy.deepcopy(Controls)
1083    phi = StrSta['Sample phi']
1084    wave = Controls['wavelength']
1085    pixelSize = Controls['pixelSize']
1086    scalex = 1000./pixelSize[0]
1087    scaley = 1000./pixelSize[1]
1088    StaType = StrSta['Type']
1089    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
1090
1091    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
1092        dset = ring['Dset']
1093        Ring,R = MakeStrStaRing(ring,Image,StaControls)
1094        if len(Ring):
1095            ring.update(R)
1096            p0 = ring['Emat']
1097            val,esd,covMat = FitStrain(Ring,p0,dset,wave,phi,StaType)
1098            ring['Emat'] = val
1099            ring['Esig'] = esd
1100            ellipse = FitEllipse(R['ImxyObs'].T)
1101            ringxy,ringazm = makeRing(ring['Dcalc'],ellipse,0,0.,scalex,scaley,Image)
1102            ring['ImxyCalc'] = np.array(ringxy).T[:2]
1103            ringint = np.array([float(Image[int(x*scalex),int(y*scaley)]) for y,x in np.array(ringxy)[:,:2]])
1104            ringint /= np.mean(ringint)
1105            ring['Ivar'] = np.var(ringint)
1106            ring['covMat'] = covMat
1107            print ('Variance in normalized ring intensity: %.3f'%(ring['Ivar']))
1108    CalcStrSta(StrSta,Controls)
1109   
1110def IntStrSta(Image,StrSta,Controls):
1111    StaControls = copy.deepcopy(Controls)
1112    pixelSize = Controls['pixelSize']
1113    scalex = 1000./pixelSize[0]
1114    scaley = 1000./pixelSize[1]
1115    phi = StrSta['Sample phi']
1116    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
1117    RingsAI = []
1118    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
1119        Ring,R = MakeStrStaRing(ring,Image,StaControls)
1120        if len(Ring):
1121            ellipse = FitEllipse(R['ImxyObs'].T)
1122            ringxy,ringazm = makeRing(ring['Dcalc'],ellipse,0,0.,scalex,scaley,Image)
1123            ring['ImxyCalc'] = np.array(ringxy).T[:2]
1124            ringint = np.array([float(Image[int(x*scalex),int(y*scaley)]) for y,x in np.array(ringxy)[:,:2]])
1125            ringint /= np.mean(ringint)
1126            print (' %s %.3f %s %.3f'%('d-spacing',ring['Dcalc'],'sig(MRD):',np.sqrt(np.var(ringint))))
1127            RingsAI.append(np.array(zip(ringazm,ringint)).T)
1128    return RingsAI
1129   
1130def CalcStrSta(StrSta,Controls):
1131
1132    wave = Controls['wavelength']
1133    phi = StrSta['Sample phi']
1134    StaType = StrSta['Type']
1135    for ring in StrSta['d-zero']:
1136        Eij = ring['Emat']
1137        E = [[Eij[0],Eij[1],0],[Eij[1],Eij[2],0],[0,0,0]]
1138        th,azm = ring['ImtaObs']
1139        th0 = np.ones_like(azm)*npasind(wave/(2.*ring['Dset']))
1140        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th0).T/1.e6,axis=2),axis=1)
1141        if StaType == 'True':
1142            ring['ImtaCalc'] = np.array([np.exp(V)*ring['Dset'],azm])
1143        else:
1144            ring['ImtaCalc'] = np.array([(V+1.)*ring['Dset'],azm])
1145        dmin = np.min(ring['ImtaCalc'][0])
1146        dmax = np.max(ring['ImtaCalc'][0])
1147        if ring.get('fixDset',True):
1148            if abs(Eij[0]) < abs(Eij[2]):         #tension
1149                ring['Dcalc'] = dmin+(dmax-dmin)/4.
1150            else:                       #compression
1151                ring['Dcalc'] = dmin+3.*(dmax-dmin)/4.
1152        else:
1153            ring['Dcalc'] = np.mean(ring['ImtaCalc'][0])
1154
1155def calcFij(omg,phi,azm,th):
1156    '''    Uses parameters as defined by Bob He & Kingsley Smith, Adv. in X-Ray Anal. 41, 501 (1997)
1157
1158    :param omg: his omega = sample omega rotation; 0 when incident beam || sample surface,
1159        90 when perp. to sample surface
1160    :param phi: his phi = sample phi rotation; usually = 0, axis rotates with omg.
1161    :param azm: his chi = azimuth around incident beam
1162    :param th:  his theta = theta
1163    '''
1164    a = npsind(th)*npcosd(omg)+npsind(azm)*npcosd(th)*npsind(omg)
1165    b = -npcosd(azm)*npcosd(th)
1166    c = npsind(th)*npsind(omg)-npsind(azm)*npcosd(th)*npcosd(omg)
1167    d = a*npsind(phi)+b*npcosd(phi)
1168    e = a*npcosd(phi)-b*npsind(phi)
1169    Fij = np.array([
1170        [d**2,d*e,c*d],
1171        [d*e,e**2,c*e],
1172        [c*d,c*e,c**2]])
1173    return -Fij
1174
1175def FitStrain(rings,p0,dset,wave,phi,StaType):
1176    'Needs a doc string'
1177    def StrainPrint(ValSig,dset):
1178        print ('Strain tensor for Dset: %.6f'%(dset))
1179        ptlbls = 'names :'
1180        ptstr =  'values:'
1181        sigstr = 'esds  :'
1182        for name,fmt,value,sig in ValSig:
1183            ptlbls += "%s" % (name.rjust(12))
1184            ptstr += fmt % (value)
1185            if sig:
1186                sigstr += fmt % (sig)
1187            else:
1188                sigstr += 12*' '
1189        print (ptlbls)
1190        print (ptstr)
1191        print (sigstr)
1192       
1193    def strainCalc(p,xyd,dset,wave,phi,StaType):
1194        E = np.array([[p[0],p[1],0],[p[1],p[2],0],[0,0,0]])
1195        dspo,azm,dsp = xyd
1196        th = npasind(wave/(2.0*dspo))
1197        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th).T/1.e6,axis=2),axis=1)
1198        if StaType == 'True':
1199            dspc = dset*np.exp(V)
1200        else:
1201            dspc = dset*(V+1.)
1202        return dspo-dspc
1203       
1204    names = ['e11','e12','e22']
1205    fmt = ['%12.2f','%12.2f','%12.2f']
1206    result = leastsq(strainCalc,p0,args=(rings,dset,wave,phi,StaType),full_output=True)
1207    vals = list(result[0])
1208    chisq = np.sum(result[2]['fvec']**2)/(rings.shape[1]-3)     #reduced chi^2 = M/(Nobs-Nvar)
1209    covM = result[1]
1210    covMatrix = covM*chisq
1211    sig = list(np.sqrt(chisq*np.diag(result[1])))
1212    ValSig = zip(names,fmt,vals,sig)
1213    StrainPrint(ValSig,dset)
1214    return vals,sig,covMatrix
1215
1216def FitImageSpots(Image,ImMax,ind,pixSize,nxy):
1217   
1218    def calcMean(nxy,pixSize,img):
1219        gdx,gdy = np.mgrid[0:nxy,0:nxy]
1220        gdx = ma.array((gdx-nxy//2)*pixSize[0]/1000.,mask=~ma.getmaskarray(ImBox))
1221        gdy = ma.array((gdy-nxy//2)*pixSize[1]/1000.,mask=~ma.getmaskarray(ImBox))
1222        posx = ma.sum(gdx)/ma.count(gdx)
1223        posy = ma.sum(gdy)/ma.count(gdy)
1224        return posx,posy
1225   
1226    def calcPeak(values,nxy,pixSize,img):
1227        back,mag,px,py,sig = values
1228        peak = np.zeros([nxy,nxy])+back
1229        nor = 1./(2.*np.pi*sig**2)
1230        gdx,gdy = np.mgrid[0:nxy,0:nxy]
1231        gdx = (gdx-nxy//2)*pixSize[0]/1000.
1232        gdy = (gdy-nxy//2)*pixSize[1]/1000.
1233        arg = (gdx-px)**2+(gdy-py)**2       
1234        peak += mag*nor*np.exp(-arg/(2.*sig**2))
1235        return ma.compressed(img-peak)/np.sqrt(ma.compressed(img))
1236   
1237    def calc2Peak(values,nxy,pixSize,img):
1238        back,mag,px,py,sigx,sigy,rho = values
1239        peak = np.zeros([nxy,nxy])+back
1240        nor = 1./(2.*np.pi*sigx*sigy*np.sqrt(1.-rho**2))
1241        gdx,gdy = np.mgrid[0:nxy,0:nxy]
1242        gdx = (gdx-nxy//2)*pixSize[0]/1000.
1243        gdy = (gdy-nxy//2)*pixSize[1]/1000.
1244        argnor = -1./(2.*(1.-rho**2))
1245        arg = (gdx-px)**2/sigx**2+(gdy-py)**2/sigy**2-2.*rho*(gdx-px)*(gdy-py)/(sigx*sigy)       
1246        peak += mag*nor*np.exp(argnor*arg)
1247        return ma.compressed(img-peak)/np.sqrt(ma.compressed(img))       
1248   
1249    nxy2 = nxy//2
1250    ImBox = Image[ind[1]-nxy2:ind[1]+nxy2+1,ind[0]-nxy2:ind[0]+nxy2+1]
1251    back = np.min(ImBox)
1252    mag = np.sum(ImBox-back)
1253    vals = [back,mag,0.,0.,0.2,0.2,0.]
1254    ImBox = ma.array(ImBox,dtype=float,mask=ImBox>0.75*ImMax)
1255    px = (ind[0]+.5)*pixSize[0]/1000.
1256    py = (ind[1]+.5)*pixSize[1]/1000.
1257    if ma.any(ma.getmaskarray(ImBox)):
1258        vals = calcMean(nxy,pixSize,ImBox)
1259        posx,posy = [px+vals[0],py+vals[1]]
1260        return [posx,posy,6.]
1261    else:
1262        result = leastsq(calc2Peak,vals,args=(nxy,pixSize,ImBox),full_output=True)
1263        vals = result[0]
1264        ratio = vals[4]/vals[5]
1265        if 0.5 < ratio < 2.0 and vals[2] < 2. and vals[3] < 2.:
1266            posx,posy = [px+vals[2],py+vals[3]]
1267            return [posx,posy,min(6.*vals[4],6.)]
1268        else:
1269            return None
1270   
1271def AutoSpotMasks(Image,Masks,Controls):
1272   
1273    print ('auto spot search')
1274    nxy = 15
1275    rollImage = lambda rho,roll: np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1)
1276    pixelSize = Controls['pixelSize']
1277    spotMask = ma.array(Image,mask=(Image<np.mean(Image)))
1278    indices = (-1,0,1)
1279    rolls = np.array([[ix,iy] for ix in indices for iy in indices])
1280    time0 = time.time()
1281    for roll in rolls:
1282        if np.any(roll):        #avoid [0,0]
1283            spotMask = ma.array(spotMask,mask=(spotMask-rollImage(Image,roll)<0.),dtype=float)
1284    mags = spotMask[spotMask.nonzero()]
1285    indx = np.transpose(spotMask.nonzero())
1286    size1 = mags.shape[0]
1287    magind = [[indx[0][0],indx[0][1],mags[0]],]
1288    for ind,mag in list(zip(indx,mags))[1:]:        #remove duplicates
1289#            ind[0],ind[1],I,J = ImageLocalMax(Image,nxy,ind[0],ind[1])
1290        if (magind[-1][0]-ind[0])**2+(magind[-1][1]-ind[1])**2 > 16:
1291            magind.append([ind[0],ind[1],Image[ind[0],ind[1]]]) 
1292    magind = np.array(magind).T
1293    indx = np.array(magind[0:2],dtype=np.int32)
1294    mags = magind[2]
1295    size2 = mags.shape[0]
1296    print ('Initial search done: %d -->%d %.2fs'%(size1,size2,time.time()-time0))
1297    nx,ny = Image.shape
1298    ImMax = np.max(Image)
1299    peaks = []
1300    nxy2 = nxy//2
1301    mult = 0.001
1302    num = 1e6
1303    while num>500:
1304        mult += .0001           
1305        minM = mult*np.max(mags)
1306        num = ma.count(ma.array(mags,mask=mags<=minM))
1307        print('try',mult,minM,num)
1308    minM = mult*np.max(mags)
1309    print ('Find biggest spots:',mult,num,minM)
1310    for i,mag in enumerate(mags):
1311        if mag > minM:
1312            if (nxy2 < indx[0][i] < nx-nxy2-1) and (nxy2 < indx[1][i] < ny-nxy2-1):
1313#                    print ('try:%d %d %d %.2f'%(i,indx[0][i],indx[1][i],mags[i]))
1314                peak = FitImageSpots(Image,ImMax,[indx[1][i],indx[0][i]],pixelSize,nxy)
1315                if peak and not any(np.isnan(np.array(peak))):
1316                    peaks.append(peak)
1317#                    print (' Spot found: %s'%str(peak))
1318    peaks = G2mth.sortArray(G2mth.sortArray(peaks,1),0)
1319    Peaks = [peaks[0],]
1320    for peak in peaks[1:]:
1321        if GetDsp(peak[0],peak[1],Controls) >= 1.:      #toss non-diamond low angle spots
1322            continue
1323        if (peak[0]-Peaks[-1][0])**2+(peak[1]-Peaks[-1][1])**2 > peak[2]*Peaks[-1][2] :
1324            Peaks.append(peak)
1325#            print (' Spot found: %s'%str(peak))
1326    print ('Spots found: %d time %.2fs'%(len(Peaks),time.time()-time0))
1327    Masks['Points'] = Peaks
1328    return None
Note: See TracBrowser for help on using the repository browser.