source: trunk/GSASIIimage.py @ 3522

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

fix array overflow problem in histogram2d.for
clean up some stuff in ImageIntegrate?

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