source: trunk/GSASIIimage.py @ 2627

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

remove progress bar stuff from image Integrate
put progress bar around integrate all
put wx.BusyCursor? around individual integrations & autointegrate

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 45.7 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASII image calculations: ellipse fitting & image integration       
3########### SVN repository information ###################
4# $Date: 2017-01-12 19:22:19 +0000 (Thu, 12 Jan 2017) $
5# $Author: vondreele $
6# $Revision: 2627 $
7# $URL: trunk/GSASIIimage.py $
8# $Id: GSASIIimage.py 2627 2017-01-12 19:22:19Z vondreele $
9########### SVN repository information ###################
10'''
11*GSASIIimage: Image calc module*
12================================
13
14Ellipse fitting & image integration
15
16'''
17
18import math
19import time
20import numpy as np
21import numpy.linalg as nl
22import numpy.ma as ma
23import polymask as pm
24from scipy.optimize import leastsq
25import scipy.signal as scsg
26import copy
27import GSASIIpath
28GSASIIpath.SetVersionNumber("$Revision: 2627 $")
29import GSASIIplot as G2plt
30import GSASIIlattice as G2lat
31import GSASIIpwd as G2pwd
32import GSASIIspc as G2spc
33import GSASIImath as G2mth
34
35# trig functions in degrees
36sind = lambda x: math.sin(x*math.pi/180.)
37asind = lambda x: 180.*math.asin(x)/math.pi
38tand = lambda x: math.tan(x*math.pi/180.)
39atand = lambda x: 180.*math.atan(x)/math.pi
40atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
41cosd = lambda x: math.cos(x*math.pi/180.)
42acosd = lambda x: 180.*math.acos(x)/math.pi
43rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
44#numpy versions
45npsind = lambda x: np.sin(x*np.pi/180.)
46npasind = lambda x: 180.*np.arcsin(x)/np.pi
47npcosd = lambda x: np.cos(x*np.pi/180.)
48npacosd = lambda x: 180.*np.arccos(x)/np.pi
49nptand = lambda x: np.tan(x*np.pi/180.)
50npatand = lambda x: 180.*np.arctan(x)/np.pi
51npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
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,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))         #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['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.2f','%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
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-1:ypix+1,xpix-1:xpix+1])
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)        #major axis
252        y = radii[0]*sind(a)
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    dxy = peneCorr(tth,dep,tilt)
317    dist = data['distance']
318    return GetEllipse2(tth,dxy,dist,cent,tilt,phi)
319       
320def GetDetectorXY(dsp,azm,data):
321    'Needs a doc string'
322   
323    elcent,phi,radii = GetEllipse(dsp,data)
324    phi = data['rotation']-90.          #to give rotation of major axis
325    tilt = data['tilt']
326    dist = data['distance']
327    cent = data['center']
328    tth = 2.0*asind(data['wavelength']/(2.*dsp))
329    stth = sind(tth)
330    cosb = cosd(tilt)
331    if radii[0] > 0.:
332        sinb = sind(tilt)
333        tanb = tand(tilt)
334        fplus = dist*tanb*stth/(cosb+stth)
335        fminus = dist*tanb*stth/(cosb-stth)
336        zdis = (fplus-fminus)/2.
337        rsqplus = radii[0]**2+radii[1]**2
338        rsqminus = radii[0]**2-radii[1]**2
339        R = rsqminus*cosd(2.*azm-2.*phi)+rsqplus
340        Q = np.sqrt(2.)*radii[0]*radii[1]*np.sqrt(R-2.*zdis**2*sind(azm-phi)**2)
341        P = 2.*radii[0]**2*zdis*cosd(azm-phi)
342        radius = (P+Q)/R
343        xy = np.array([radius*cosd(azm),radius*sind(azm)])
344        xy += cent
345    else:   #hyperbola - both branches (one is way off screen!)
346        sinb = abs(sind(tilt))
347        tanb = abs(tand(tilt))
348        f = dist*tanb*stth/(cosb+stth)
349        v = dist*(tanb+tand(tth-abs(tilt)))
350        delt = dist*stth*(1+stth*cosb)/(sinb*cosb*(stth+cosb))
351        ecc = (v-f)/(delt-v)
352        R = radii[1]*(ecc**2-1)/(1-ecc*cosd(azm))
353        if tilt > 0.:
354            offset = 2.*radii[1]*ecc+f      #select other branch
355            xy = [-R*cosd(azm)-offset,R*sind(azm)]
356        else:
357            offset = -f
358            xy = [-R*cosd(azm)-offset,-R*sind(azm)]
359        xy = -np.array([xy[0]*cosd(phi)+xy[1]*sind(phi),xy[0]*sind(phi)-xy[1]*cosd(phi)])
360        xy += cent
361    return xy
362   
363def GetDetXYfromThAzm(Th,Azm,data):
364    'Needs a doc string'
365    dsp = data['wavelength']/(2.0*npsind(Th))   
366    return GetDetectorXY(dsp,Azm,data)
367                   
368def GetTthAzmDsp(x,y,data): #expensive
369    'Needs a doc string - checked OK for ellipses & hyperbola'
370    wave = data['wavelength']
371    cent = data['center']
372    tilt = data['tilt']
373    dist = data['distance']/cosd(tilt)
374    x0 = dist*tand(tilt)
375    phi = data['rotation']
376    dep = data['DetDepth']
377    azmthoff = data['azmthOff']
378    dx = np.array(x-cent[0],dtype=np.float32)
379    dy = np.array(y-cent[1],dtype=np.float32)
380    D = ((dx-x0)**2+dy**2+dist**2)      #sample to pixel distance
381    X = np.array(([dx,dy,np.zeros_like(dx)]),dtype=np.float32).T
382    X = np.dot(X,makeMat(phi,2))
383    Z = np.dot(X,makeMat(tilt,0)).T[2]
384    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z))
385    dxy = peneCorr(tth,dep,tilt,npatan2d(dy,dx))
386    DX = dist-Z+dxy
387    DY = np.sqrt(dx**2+dy**2-Z**2)
388    tth = npatan2d(DY,DX) 
389    dsp = wave/(2.*npsind(tth/2.))
390    azm = (npatan2d(dy,dx)+azmthoff+720.)%360.
391    G = D/dist**2       #for geometric correction = 1/cos(2theta)^2 if tilt=0.
392    return np.array([tth,azm,G,dsp])
393   
394def GetTth(x,y,data):
395    'Give 2-theta value for detector x,y position; calibration info in data'
396    return GetTthAzmDsp(x,y,data)[0]
397   
398def GetTthAzm(x,y,data):
399    'Give 2-theta, azimuth values for detector x,y position; calibration info in data'
400    return GetTthAzmDsp(x,y,data)[0:2]
401   
402def GetTthAzmG(x,y,data):
403    '''Give 2-theta, azimuth & geometric corr. values for detector x,y position;
404     calibration info in data - only used in integration
405    '''
406    'Needs a doc string - checked OK for ellipses & hyperbola'
407    tilt = data['tilt']
408    dist = data['distance']/npcosd(tilt)
409    x0 = data['distance']*nptand(tilt)
410    MN = -np.inner(makeMat(data['rotation'],2),makeMat(tilt,0))
411    distsq = data['distance']**2
412    dx = x-data['center'][0]
413    dy = y-data['center'][1]
414    G = ((dx-x0)**2+dy**2+distsq)/distsq       #for geometric correction = 1/cos(2theta)^2 if tilt=0.
415    X = np.dstack([dx.T,dy.T,np.zeros_like(dx.T)])
416    Z = np.dot(X,MN).T[2]
417    xyZ = dx**2+dy**2-Z**2   
418    tth = npatand(np.sqrt(xyZ)/(dist-Z))
419    dxy = peneCorr(tth,data['DetDepth'],tilt,npatan2d(dy,dx))
420    tth = npatan2d(np.sqrt(xyZ),dist-Z+dxy) 
421    azm = (npatan2d(dy,dx)+data['azmthOff']+720.)%360.
422    return tth,azm,G
423
424def GetDsp(x,y,data):
425    'Give d-spacing value for detector x,y position; calibration info in data'
426    return GetTthAzmDsp(x,y,data)[3]
427       
428def GetAzm(x,y,data):
429    'Give azimuth value for detector x,y position; calibration info in data'
430    return GetTthAzmDsp(x,y,data)[1]
431   
432def meanAzm(a,b):
433    AZM = lambda a,b: npacosd(0.5*(npsind(2.*b)-npsind(2.*a))/(np.pi*(b-a)/180.))/2.
434    azm = AZM(a,b)
435#    quad = int((a+b)/180.)
436#    if quad == 1:
437#        azm = 180.-azm
438#    elif quad == 2:
439#        azm += 180.
440#    elif quad == 3:
441#        azm = 360-azm
442    return azm     
443       
444def ImageCompress(image,scale):
445    ''' Reduces size of image by selecting every n'th point
446    param: image array: original image
447    param: scale int: intervsl between selected points
448    returns: array: reduced size image
449    '''
450    if scale == 1:
451        return image
452    else:
453        return image[::scale,::scale]
454       
455def checkEllipse(Zsum,distSum,xSum,ySum,dist,x,y):
456    'Needs a doc string'
457    avg = np.array([distSum/Zsum,xSum/Zsum,ySum/Zsum])
458    curr = np.array([dist,x,y])
459    return abs(avg-curr)/avg < .02
460
461def EdgeFinder(image,data):
462    '''this makes list of all x,y where I>edgeMin suitable for an ellipse search?
463    Not currently used but might be useful in future?
464    '''
465    import numpy.ma as ma
466    Nx,Ny = data['size']
467    pixelSize = data['pixelSize']
468    edgemin = data['edgemin']
469    scalex = pixelSize[0]/1000.
470    scaley = pixelSize[1]/1000.   
471    tay,tax = np.mgrid[0:Nx,0:Ny]
472    tax = np.asfarray(tax*scalex,dtype=np.float32)
473    tay = np.asfarray(tay*scaley,dtype=np.float32)
474    tam = ma.getmask(ma.masked_less(image.flatten(),edgemin))
475    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))
476    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))
477    return zip(tax,tay)
478   
479def MakeFrameMask(data,frame):
480    pixelSize = data['pixelSize']
481    scalex = pixelSize[0]/1000.
482    scaley = pixelSize[1]/1000.
483    blkSize = 512
484    Nx,Ny = data['size']
485    nXBlks = (Nx-1)/blkSize+1
486    nYBlks = (Ny-1)/blkSize+1
487    tam = ma.make_mask_none(data['size'])
488    for iBlk in range(nXBlks):
489        iBeg = iBlk*blkSize
490        iFin = min(iBeg+blkSize,Nx)
491        for jBlk in range(nYBlks):
492            jBeg = jBlk*blkSize
493            jFin = min(jBeg+blkSize,Ny)               
494            nI = iFin-iBeg
495            nJ = jFin-jBeg
496            tax,tay = np.mgrid[iBeg+0.5:iFin+.5,jBeg+.5:jFin+.5]         #bin centers not corners
497            tax = np.asfarray(tax*scalex,dtype=np.float32)
498            tay = np.asfarray(tay*scaley,dtype=np.float32)
499            tamp = ma.make_mask_none((1024*1024))
500            tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
501                tay.flatten(),len(frame),frame,tamp)[:nI*nJ])-True  #switch to exclude around frame
502            if tamp.shape:
503                tamp = np.reshape(tamp[:nI*nJ],(nI,nJ))
504                tam[iBeg:iFin,jBeg:jFin] = ma.mask_or(tamp[0:nI,0:nJ],tam[iBeg:iFin,jBeg:jFin])
505            else:
506                tam[iBeg:iFin,jBeg:jFin] = True
507    return tam.T
508   
509def ImageRecalibrate(G2frame,data,masks):
510    '''Called to repeat the calibration on an image, usually called after
511    calibration is done initially to improve the fit.
512    '''
513    import ImageCalibrants as calFile
514    print 'Image recalibration:'
515    time0 = time.time()
516    pixelSize = data['pixelSize']
517    scalex = 1000./pixelSize[0]
518    scaley = 1000./pixelSize[1]
519    pixLimit = data['pixLimit']
520    cutoff = data['cutoff']
521    data['rings'] = []
522    data['ellipses'] = []
523    if not data['calibrant']:
524        print 'no calibration material selected'
525        return True   
526    skip = data['calibskip']
527    dmin = data['calibdmin']
528    Bravais,SGs,Cells = calFile.Calibrants[data['calibrant']][:3]
529    HKL = []
530    for bravais,sg,cell in zip(Bravais,SGs,Cells):
531        A = G2lat.cell2A(cell)
532        if sg:
533            SGData = G2spc.SpcGroup(sg)[1]
534            hkl = G2pwd.getHKLpeak(dmin,SGData,A)
535            HKL += hkl
536        else:
537            hkl = G2lat.GenHBravais(dmin,bravais,A)
538            HKL += hkl
539    HKL = G2lat.sortHKLd(HKL,True,False)
540    varyList = [item for item in data['varyList'] if data['varyList'][item]]
541    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
542        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
543    Found = False
544    wave = data['wavelength']
545    frame = masks['Frames']
546    tam = ma.make_mask_none(G2frame.ImageZ.shape)
547    if frame:
548        tam = ma.mask_or(tam,MakeFrameMask(data,frame))
549    for iH,H in enumerate(HKL):
550        if debug:   print H
551        dsp = H[3]
552        tth = 2.0*asind(wave/(2.*dsp))
553        if tth+abs(data['tilt']) > 90.:
554            print 'next line is a hyperbola - search stopped'
555            break
556        ellipse = GetEllipse(dsp,data)
557        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,ma.array(G2frame.ImageZ,mask=tam))[0]
558        if Ring:
559            if iH >= skip:
560                data['rings'].append(np.array(Ring))
561            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
562            Found = True
563        elif not Found:         #skipping inner rings, keep looking until ring found
564            continue
565        else:                   #no more rings beyond edge of detector
566            data['ellipses'].append([])
567            continue
568    if not data['rings']:
569        print 'no rings found; try lower Min ring I/Ib'
570        return True   
571       
572    rings = np.concatenate((data['rings']),axis=0)
573    chisq = FitDetector(rings,varyList,parmDict)
574    data['wavelength'] = parmDict['wave']
575    data['distance'] = parmDict['dist']
576    data['center'] = [parmDict['det-X'],parmDict['det-Y']]
577    data['rotation'] = np.mod(parmDict['phi'],360.0)
578    data['tilt'] = parmDict['tilt']
579    data['DetDepth'] = parmDict['dep']
580    data['chisq'] = chisq
581    N = len(data['ellipses'])
582    data['ellipses'] = []           #clear away individual ellipse fits
583    for H in HKL[:N]:
584        ellipse = GetEllipse(H[3],data)
585        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))   
586    print 'calibration time = %.3f'%(time.time()-time0)
587    G2plt.PlotImage(G2frame,newImage=True)       
588    return True
589           
590def ImageCalibrate(G2frame,data):
591    '''Called to perform an initial image calibration after points have been
592    selected for the inner ring.
593    '''
594    import copy
595    import ImageCalibrants as calFile
596    print 'Image calibration:'
597    time0 = time.time()
598    ring = data['ring']
599    pixelSize = data['pixelSize']
600    scalex = 1000./pixelSize[0]
601    scaley = 1000./pixelSize[1]
602    pixLimit = data['pixLimit']
603    cutoff = data['cutoff']
604    varyDict = data['varyList']
605    if varyDict['dist'] and varyDict['wave']:
606        print 'ERROR - you can not simultaneously calibrate distance and wavelength'
607        return False
608    if len(ring) < 5:
609        print 'ERROR - not enough inner ring points for ellipse'
610        return False
611       
612    #fit start points on inner ring
613    data['ellipses'] = []
614    data['rings'] = []
615    outE = FitEllipse(ring)
616    fmt  = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f'
617    fmt2 = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f, chi**2: %.3f, Np: %d'
618    if outE:
619        print fmt%('start ellipse: ',outE[0][0],outE[0][1],outE[1],outE[2][0],outE[2][1])
620        ellipse = outE
621    else:
622        return False
623       
624    #setup 360 points on that ring for "good" fit
625    data['ellipses'].append(ellipse[:]+('g',))
626    Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0]
627    if Ring:
628        ellipse = FitEllipse(Ring)
629        Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0]    #do again
630        ellipse = FitEllipse(Ring)
631    else:
632        print '1st ring not sufficiently complete to proceed'
633        return False
634    if debug:
635        print fmt2%('inner ring:    ',ellipse[0][0],ellipse[0][1],ellipse[1],
636            ellipse[2][0],ellipse[2][1],0.,len(Ring))     #cent,phi,radii
637    data['ellipses'].append(ellipse[:]+('r',))
638    data['rings'].append(np.array(Ring))
639    G2plt.PlotImage(G2frame,newImage=True)
640   
641#setup for calibration
642    data['rings'] = []
643    if not data['calibrant']:
644        print 'no calibration material selected'
645        return True
646   
647    skip = data['calibskip']
648    dmin = data['calibdmin']
649#generate reflection set
650    Bravais,SGs,Cells = calFile.Calibrants[data['calibrant']][:3]
651    HKL = []
652    for bravais,sg,cell in zip(Bravais,SGs,Cells):
653        A = G2lat.cell2A(cell)
654        if sg:
655            SGData = G2spc.SpcGroup(sg)[1]
656            hkl = G2pwd.getHKLpeak(dmin,SGData,A)
657            HKL += hkl
658        else:
659            hkl = G2lat.GenHBravais(dmin,bravais,A)
660            HKL += hkl
661    HKL = G2lat.sortHKLd(HKL,True,False)[skip:]
662#set up 1st ring
663    elcent,phi,radii = ellipse              #from fit of 1st ring
664    dsp = HKL[0][3]
665    print '1st ring: try %.4f'%(dsp)
666    if varyDict['dist']:
667        wave = data['wavelength']
668        tth = 2.0*asind(wave/(2.*dsp))
669    else:   #varyDict['wave']!
670        dist = data['distance']
671        tth = npatan2d(radii[0],dist)
672        data['wavelength'] = wave =  2.0*dsp*sind(tth/2.0)
673    Ring0 = makeRing(dsp,ellipse,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
674    ttth = nptand(tth)
675    ctth = npcosd(tth)
676#1st estimate of tilt; assume ellipse - don't know sign though
677    if varyDict['tilt']:
678        tilt = npasind(np.sqrt(max(0.,1.-(radii[0]/radii[1])**2))*ctth)
679        if not tilt:
680            print 'WARNING - selected ring was fitted as a circle'
681            print ' - if detector was tilted we suggest you skip this ring - WARNING'
682    else:
683        tilt = data['tilt']
684#1st estimate of dist: sample to detector normal to plane
685    if varyDict['dist']:
686        data['distance'] = dist = radii[0]**2/(ttth*radii[1])
687    else:
688        dist = data['distance']
689    if varyDict['tilt']:
690#ellipse to cone axis (x-ray beam); 2 choices depending on sign of tilt
691        zdisp = radii[1]*ttth*tand(tilt)
692        zdism = radii[1]*ttth*tand(-tilt)
693#cone axis position; 2 choices. Which is right?     
694#NB: zdisp is || to major axis & phi is rotation of minor axis
695#thus shift from beam to ellipse center is [Z*sin(phi),-Z*cos(phi)]
696        centp = [elcent[0]+zdisp*sind(phi),elcent[1]-zdisp*cosd(phi)]
697        centm = [elcent[0]+zdism*sind(phi),elcent[1]-zdism*cosd(phi)]
698#check get same ellipse parms either way
699#now do next ring; estimate either way & do a FitDetector each way; best fit is correct one
700        fail = True
701        i2 = 1
702        while fail:
703            dsp = HKL[i2][3]
704            print '2nd ring: try %.4f'%(dsp)
705            tth = 2.0*asind(wave/(2.*dsp))
706            ellipsep = GetEllipse2(tth,0.,dist,centp,tilt,phi)
707            print fmt%('plus ellipse :',ellipsep[0][0],ellipsep[0][1],ellipsep[1],ellipsep[2][0],ellipsep[2][1])
708            Ringp = makeRing(dsp,ellipsep,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
709            parmDict = {'dist':dist,'det-X':centp[0],'det-Y':centp[1],
710                'tilt':tilt,'phi':phi,'wave':wave,'dep':0.0}       
711            varyList = [item for item in varyDict if varyDict[item]]
712            if len(Ringp) > 10:
713                chip = FitDetector(np.array(Ring0+Ringp),varyList,parmDict,True)
714                tiltp = parmDict['tilt']
715                phip = parmDict['phi']
716                centp = [parmDict['det-X'],parmDict['det-Y']]
717                fail = False
718            else:
719                chip = 1e6
720            ellipsem = GetEllipse2(tth,0.,dist,centm,-tilt,phi)
721            print fmt%('minus ellipse:',ellipsem[0][0],ellipsem[0][1],ellipsem[1],ellipsem[2][0],ellipsem[2][1])
722            Ringm = makeRing(dsp,ellipsem,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
723            if len(Ringm) > 10:
724                parmDict['tilt'] *= -1
725                chim = FitDetector(np.array(Ring0+Ringm),varyList,parmDict,True)
726                tiltm = parmDict['tilt']
727                phim = parmDict['phi']
728                centm = [parmDict['det-X'],parmDict['det-Y']]
729                fail = False
730            else:
731                chim = 1e6
732            if fail:
733                i2 += 1
734        if chip < chim:
735            data['tilt'] = tiltp
736            data['center'] = centp
737            data['rotation'] = phip
738        else:
739            data['tilt'] = tiltm
740            data['center'] = centm
741            data['rotation'] = phim
742        data['ellipses'].append(ellipsep[:]+('b',))
743        data['rings'].append(np.array(Ringp))
744        data['ellipses'].append(ellipsem[:]+('r',))
745        data['rings'].append(np.array(Ringm))
746        G2plt.PlotImage(G2frame,newImage=True)
747    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
748        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
749    varyList = [item for item in varyDict if varyDict[item]]
750    data['rings'] = []
751    data['ellipses'] = []
752    for i,H in enumerate(HKL):
753        dsp = H[3]
754        tth = 2.0*asind(wave/(2.*dsp))
755        if tth+abs(data['tilt']) > 90.:
756            print 'next line is a hyperbola - search stopped'
757            break
758        if debug:   print 'HKLD:',H[:4],'2-theta: %.4f'%(tth)
759        elcent,phi,radii = ellipse = GetEllipse(dsp,data)
760        data['ellipses'].append(copy.deepcopy(ellipse+('g',)))
761        if debug:   print fmt%('predicted ellipse:',elcent[0],elcent[1],phi,radii[0],radii[1])
762        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0]
763        if Ring:
764            data['rings'].append(np.array(Ring))
765            rings = np.concatenate((data['rings']),axis=0)
766            if i:
767                chisq = FitDetector(rings,varyList,parmDict,False)
768                data['distance'] = parmDict['dist']
769                data['center'] = [parmDict['det-X'],parmDict['det-Y']]
770                data['rotation'] = parmDict['phi']
771                data['tilt'] = parmDict['tilt']
772                data['DetDepth'] = parmDict['dep']
773                data['chisq'] = chisq
774                elcent,phi,radii = ellipse = GetEllipse(dsp,data)
775                if debug:   print fmt2%('fitted ellipse:   ',elcent[0],elcent[1],phi,radii[0],radii[1],chisq,len(rings))
776            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
777#            G2plt.PlotImage(G2frame,newImage=True)
778        else:
779            if debug:   print 'insufficient number of points in this ellipse to fit'
780#            break
781    G2plt.PlotImage(G2frame,newImage=True)
782    fullSize = len(G2frame.ImageZ)/scalex
783    if 2*radii[1] < .9*fullSize:
784        print 'Are all usable rings (>25% visible) used? Try reducing Min ring I/Ib'
785    N = len(data['ellipses'])
786    if N > 2:
787        FitDetector(rings,varyList,parmDict)
788        data['wavelength'] = parmDict['wave']
789        data['distance'] = parmDict['dist']
790        data['center'] = [parmDict['det-X'],parmDict['det-Y']]
791        data['rotation'] = parmDict['phi']
792        data['tilt'] = parmDict['tilt']
793        data['DetDepth'] = parmDict['dep']
794    for H in HKL[:N]:
795        ellipse = GetEllipse(H[3],data)
796        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
797    print 'calibration time = %.3f'%(time.time()-time0)
798    G2plt.PlotImage(G2frame,newImage=True)       
799    return True
800   
801def Make2ThetaAzimuthMap(data,masks,iLim,jLim,times): #most expensive part of integration!
802    'Needs a doc string'
803    #transforms 2D image from x,y space to 2-theta,azimuth space based on detector orientation
804    pixelSize = data['pixelSize']
805    scalex = pixelSize[0]/1000.
806    scaley = pixelSize[1]/1000.
807   
808    tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5]         #bin centers not corners
809    tax = np.asfarray(tax*scalex,dtype=np.float32)
810    tay = np.asfarray(tay*scaley,dtype=np.float32)
811    nI = iLim[1]-iLim[0]
812    nJ = jLim[1]-jLim[0]
813    t0 = time.time()
814    #make position masks here
815    frame = masks['Frames']
816    tam = ma.make_mask_none((nI,nJ))
817    if frame:
818        tamp = ma.make_mask_none((1024*1024))
819        tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
820            tay.flatten(),len(frame),frame,tamp)[:nI*nJ])-True  #switch to exclude around frame
821        tam = ma.mask_or(tam.flatten(),tamp)
822    polygons = masks['Polygons']
823    for polygon in polygons:
824        if polygon:
825            tamp = ma.make_mask_none((1024*1024))
826            tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
827                tay.flatten(),len(polygon),polygon,tamp)[:nI*nJ])
828            tam = ma.mask_or(tam.flatten(),tamp)
829    if tam.shape: tam = np.reshape(tam,(nI,nJ))
830    spots = masks['Points']
831    for X,Y,diam in spots:
832        tamp = ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,(diam/2.)**2))
833        tam = ma.mask_or(tam,tamp)
834    times[0] += time.time()-t0
835    t0 = time.time()
836    TA = np.array(GetTthAzmG(tax,tay,data))     #includes geom. corr. as dist**2/d0**2 - most expensive step
837    times[1] += time.time()-t0
838    TA[1] = np.where(TA[1]<0,TA[1]+360,TA[1])
839    return np.array(TA),tam           #2-theta, azimuth & geom. corr. arrays & position mask
840
841def Fill2ThetaAzimuthMap(masks,TA,tam,image):
842    'Needs a doc string'
843    Zlim = masks['Thresholds'][1]
844    rings = masks['Rings']
845    arcs = masks['Arcs']
846    TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0]),ma.getdata(TA[2])))    #azimuth, 2-theta, dist
847    tax,tay,tad = np.dsplit(TA,3)    #azimuth, 2-theta, dist**2/d0**2
848    for tth,thick in rings:
849        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.)))
850    for tth,azm,thick in arcs:
851        tamt = ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.))
852        tama = ma.getmask(ma.masked_inside(tax.flatten(),azm[0],azm[1]))
853        tam = ma.mask_or(tam.flatten(),tamt*tama)
854    taz = ma.masked_outside(image.flatten(),int(Zlim[0]),Zlim[1])
855    tabs = np.ones_like(taz)
856    tam = ma.mask_or(tam.flatten(),ma.getmask(taz))
857    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))   #azimuth
858    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))   #2-theta
859    taz = ma.compressed(ma.array(taz.flatten(),mask=tam))   #intensity
860    tad = ma.compressed(ma.array(tad.flatten(),mask=tam))   #dist**2/d0**2
861    tabs = ma.compressed(ma.array(tabs.flatten(),mask=tam)) #ones - later used for absorption corr.
862    return tax,tay,taz,tad,tabs
863   
864def ImageIntegrate(image,data,masks,blkSize=128,returnN=False):
865    'Integrate an image; called from OnIntegrateAll and OnIntegrate in G2imgGUI'    #for q, log(q) bins need data['binType']
866    import histogram2d as h2d
867    print 'Begin image integration'
868    CancelPressed = False
869    LUtth = np.array(data['IOtth'])
870    LRazm = np.array(data['LRazimuth'],dtype=np.float64)
871    numAzms = data['outAzimuths']
872    numChans = data['outChannels']
873    Dazm = (LRazm[1]-LRazm[0])/numAzms
874    if '2-theta' in data.get('binType','2-theta'):
875        lutth = LUtth               
876    elif 'Q' == data['binType']:
877        lutth = 4.*np.pi*npsind(LUtth/2.)/data['wavelength']
878    elif 'log(q)' in data['binType']:
879        lutth = np.log(4.*np.pi*npsind(LUtth/2.)/data['wavelength'])
880    dtth = (lutth[1]-lutth[0])/numChans
881    muT = data.get('SampleAbs',[0.0,''])[0]
882    if 'SASD' in data['type']:
883        muT = -np.log(muT)/2.       #Transmission to 1/2 thickness muT
884    NST = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
885    H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
886    Nx,Ny = data['size']
887    nXBlks = (Nx-1)/blkSize+1
888    nYBlks = (Ny-1)/blkSize+1
889    Nup = nXBlks*nYBlks*3+3
890    tbeg = time.time()
891    Nup = 0
892    times = [0,0,0,0,0]
893    for iBlk in range(nYBlks):
894        iBeg = iBlk*blkSize
895        iFin = min(iBeg+blkSize,Ny)
896        for jBlk in range(nXBlks):
897            jBeg = jBlk*blkSize
898            jFin = min(jBeg+blkSize,Nx)
899            # next is most expensive step!
900            TA,tam = Make2ThetaAzimuthMap(data,masks,(iBeg,iFin),(jBeg,jFin),times)           #2-theta & azimuth arrays & create position mask
901            Nup += 1
902            Block = image[iBeg:iFin,jBeg:jFin]
903            t0 = time.time()
904            tax,tay,taz,tad,tabs = Fill2ThetaAzimuthMap(masks,TA,tam,Block)    #and apply masks
905            del TA; del tam
906            times[2] += time.time()-t0
907            Nup += 1
908            tax = np.where(tax > LRazm[1],tax-360.,tax)                 #put azm inside limits if possible
909            tax = np.where(tax < LRazm[0],tax+360.,tax)
910            if data.get('SampleAbs',[0.0,''])[1]:
911                if 'Cylind' in data['SampleShape']:
912                    muR = muT*(1.+npsind(tax)**2/2.)/(npcosd(tay))      #adjust for additional thickness off sample normal
913                    tabs = G2pwd.Absorb(data['SampleShape'],muR,tay)
914                elif 'Fixed' in data['SampleShape']:    #assumes flat plate sample normal to beam
915                    tabs = G2pwd.Absorb('Fixed',muT,tay)
916            if 'log(q)' in data.get('binType',''):
917                tay = np.log(4.*np.pi*npsind(tay/2.)/data['wavelength'])
918            elif 'Q' == data.get('binType',''):
919                tay = 4.*np.pi*npsind(tay/2.)/data['wavelength']
920            t0 = time.time()
921            taz = np.array((taz*tad/tabs),dtype='float32')
922            if any([tax.shape[0],tay.shape[0],taz.shape[0]]):
923                NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,
924                    numAzms,numChans,LRazm,lutth,Dazm,dtth,NST,H0)
925            times[3] += time.time()-t0
926            Nup += 1
927            del tax; del tay; del taz; del tad; del tabs
928    t0 = time.time()
929    NST = np.array(NST,dtype=np.float)
930    H0 = np.divide(H0,NST)
931    H0 = np.nan_to_num(H0)
932    H2 = np.array([tth for tth in np.linspace(lutth[0],lutth[1],numChans+1)])
933    if 'log(q)' in data.get('binType',''):
934        H2 = 2.*npasind(np.exp(H2)*data['wavelength']/(4.*np.pi))
935    elif 'Q' == data.get('binType',''):
936        H2 = 2.*npasind(H2*data['wavelength']/(4.*np.pi))
937    if Dazm:       
938        H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)])
939    else:
940        H1 = LRazm
941    H0 /= npcosd(H2[:-1])           #**2? I don't think so, **1 is right for powders
942    if 'SASD' in data['type']:
943        H0 /= npcosd(H2[:-1])           #one more for small angle scattering data?
944    if data['Oblique'][1]:
945        H0 /= G2pwd.Oblique(data['Oblique'][0],H2[:-1])
946    if 'SASD' in data['type'] and data['PolaVal'][1]:
947        #NB: in G2pwd.Polarization azm is defined from plane of polarization, not image x axis!
948        H0 /= np.array([G2pwd.Polarization(data['PolaVal'][0],H2[:-1],Azm=azm-90.)[0] for azm in (H1[:-1]+np.diff(H1)/2.)])
949    Nup += 1
950    times[4] += time.time()-t0
951    print 'Step times: \n apply masks  %8.3fs xy->th,azm   %8.3fs fill map     %8.3fs \
952        \n binning      %8.3fs cleanup      %8.3fs'%(times[0],times[1],times[2],times[3],times[4])
953    print "Elapsed time:","%8.3fs"%(time.time()-tbeg)
954    print 'Integration complete'
955    if returnN:     #As requested by Steven Weigand
956        return H0,H1,H2,NST,CancelPressed
957    else:
958        return H0,H1,H2,CancelPressed
959   
960def MakeStrStaRing(ring,Image,Controls):
961    ellipse = GetEllipse(ring['Dset'],Controls)
962    pixSize = Controls['pixelSize']
963    scalex = 1000./pixSize[0]
964    scaley = 1000./pixSize[1]
965    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
966    if len(Ring):
967        ring['ImxyObs'] = copy.copy(Ring[:2])
968        TA = GetTthAzm(Ring[0],Ring[1],Controls)       #convert x,y to tth,azm
969        TA[0] = Controls['wavelength']/(2.*npsind(TA[0]/2.))      #convert 2th to d
970        ring['ImtaObs'] = TA
971        ring['ImtaCalc'] = np.zeros_like(ring['ImtaObs'])
972        Ring[0] = TA[0]
973        Ring[1] = TA[1]
974        return Ring,ring
975    else:
976        ring['ImxyObs'] = [[],[]]
977        ring['ImtaObs'] = [[],[]]
978        ring['ImtaCalc'] = [[],[]]
979        return [],[]    #bad ring; no points found
980   
981def FitStrSta(Image,StrSta,Controls):
982    'Needs a doc string'
983   
984    StaControls = copy.deepcopy(Controls)
985    phi = StrSta['Sample phi']
986    wave = Controls['wavelength']
987    pixelSize = Controls['pixelSize']
988    scalex = 1000./pixelSize[0]
989    scaley = 1000./pixelSize[1]
990    StaType = StrSta['Type']
991    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
992
993    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
994        dset = ring['Dset']
995        Ring,R = MakeStrStaRing(ring,Image,StaControls)
996        if len(Ring):
997            ring.update(R)
998            p0 = ring['Emat']
999            val,esd = FitStrain(Ring,p0,dset,wave,phi,StaType)
1000            ring['Emat'] = val
1001            ring['Esig'] = esd
1002            ellipse = FitEllipse(R['ImxyObs'].T)
1003            ringxy,ringazm = makeRing(ring['Dcalc'],ellipse,0,0.,scalex,scaley,Image)
1004            ringint = np.array([float(Image[int(y*scaley),int(x*scalex)]) for x,y in np.array(ringxy)[:,:2]])
1005            ringint /= np.mean(ringint)
1006            ring['Ivar'] = np.var(ringint)
1007            print 'Variance in normalized ring intensity: %.3f'%(ring['Ivar'])
1008    CalcStrSta(StrSta,Controls)
1009   
1010def IntStrSta(Image,StrSta,Controls):
1011    StaControls = copy.deepcopy(Controls)
1012    pixelSize = Controls['pixelSize']
1013    scalex = 1000./pixelSize[0]
1014    scaley = 1000./pixelSize[1]
1015    phi = StrSta['Sample phi']
1016    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
1017    RingsAI = []
1018    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
1019        Ring,R = MakeStrStaRing(ring,Image,StaControls)
1020        if len(Ring):
1021            ellipse = FitEllipse(R['ImxyObs'].T)
1022            ringxy,ringazm = makeRing(ring['Dcalc'],ellipse,0,0.,scalex,scaley,Image)
1023            ringint = np.array([float(Image[int(y*scaley),int(x*scalex)]) for x,y in np.array(ringxy)[:,:2]])
1024            ringint /= np.mean(ringint)
1025            print ' %s %.3f %s %.3f'%('d-spacing',ring['Dcalc'],'var(MRD):',np.var(ringint))
1026            RingsAI.append(np.array(zip(ringazm,ringint)).T)
1027#            GSASIIpath.IPyBreak()
1028    return RingsAI
1029   
1030def CalcStrSta(StrSta,Controls):
1031
1032    wave = Controls['wavelength']
1033    phi = StrSta['Sample phi']
1034    StaType = StrSta['Type']
1035    for ring in StrSta['d-zero']:
1036        Eij = ring['Emat']
1037        E = [[Eij[0],Eij[1],0],[Eij[1],Eij[2],0],[0,0,0]]
1038        th,azm = ring['ImtaObs']
1039        th0 = np.ones_like(azm)*npasind(wave/(2.*ring['Dset']))
1040        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th0).T/1.e6,axis=2),axis=1)
1041        if StaType == 'True':
1042            ring['ImtaCalc'] = np.array([np.exp(V)*ring['Dset'],azm])
1043        else:
1044            ring['ImtaCalc'] = np.array([(V+1.)*ring['Dset'],azm])
1045        dmin = np.min(ring['ImtaCalc'][0])
1046        dmax = np.max(ring['ImtaCalc'][0])
1047        if abs(Eij[0]) < abs(Eij[2]):         #tension
1048            ring['Dcalc'] = dmin+(dmax-dmin)/4.
1049        else:                       #compression
1050            ring['Dcalc'] = dmin+3.*(dmax-dmin)/4.
1051#        ring['Dcalc'] = np.mean(ring['ImtaCalc'][0])
1052
1053def calcFij(omg,phi,azm,th):
1054    '''    Uses parameters as defined by Bob He & Kingsley Smith, Adv. in X-Ray Anal. 41, 501 (1997)
1055
1056    :param omg: his omega = sample omega rotation; 0 when incident beam || sample surface,
1057        90 when perp. to sample surface
1058    :param phi: his phi = sample phi rotation; usually = 0, axis rotates with omg.
1059    :param azm: his chi = azimuth around incident beam
1060    :param th:  his theta = theta
1061    '''
1062    a = npsind(th)*npcosd(omg)+npsind(azm)*npcosd(th)*npsind(omg)
1063    b = -npcosd(azm)*npcosd(th)
1064    c = npsind(th)*npsind(omg)-npsind(azm)*npcosd(th)*npcosd(omg)
1065    d = a*npsind(phi)+b*npcosd(phi)
1066    e = a*npcosd(phi)-b*npsind(phi)
1067    Fij = np.array([
1068        [d**2,d*e,c*d],
1069        [d*e,e**2,c*e],
1070        [c*d,c*e,c**2]])
1071    return -Fij
1072
1073def FitStrain(rings,p0,dset,wave,phi,StaType):
1074    'Needs a doc string'
1075    def StrainPrint(ValSig,dset):
1076        print 'Strain tensor for Dset: %.6f'%(dset)
1077        ptlbls = 'names :'
1078        ptstr =  'values:'
1079        sigstr = 'esds  :'
1080        for name,fmt,value,sig in ValSig:
1081            ptlbls += "%s" % (name.rjust(12))
1082            ptstr += fmt % (value)
1083            if sig:
1084                sigstr += fmt % (sig)
1085            else:
1086                sigstr += 12*' '
1087        print ptlbls
1088        print ptstr
1089        print sigstr
1090       
1091    def strainCalc(p,xyd,dset,wave,phi,StaType):
1092        E = np.array([[p[0],p[1],0],[p[1],p[2],0],[0,0,0]])
1093        dspo,azm,dsp = xyd
1094        th = npasind(wave/(2.0*dspo))
1095        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th).T/1.e6,axis=2),axis=1)
1096        if StaType == 'True':
1097            dspc = dset*np.exp(V)
1098        else:
1099            dspc = dset*(V+1.)
1100        return dspo-dspc
1101       
1102    names = ['e11','e12','e22']
1103    fmt = ['%12.2f','%12.2f','%12.2f']
1104    result = leastsq(strainCalc,p0,args=(rings,dset,wave,phi,StaType),full_output=True)
1105    vals = list(result[0])
1106    chisq = np.sum(result[2]['fvec']**2)/(rings.shape[1]-3)     #reduced chi^2 = M/(Nobs-Nvar)
1107    sig = list(np.sqrt(chisq*np.diag(result[1])))
1108    ValSig = zip(names,fmt,vals,sig)
1109    StrainPrint(ValSig,dset)
1110    return vals,sig
1111   
1112def AutoSpotMasks(Image,Masks,Controls):
1113    rollImage = lambda rho,roll: np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1)
1114    print 'auto spot search'
1115    pixelSize = Controls['pixelSize']
1116    spotMask = ma.array(Image,mask=(Image<np.mean(Image)))
1117    indices = (-1,0,1)
1118    rolls = np.array([[ix,iy] for ix in indices for iy in indices])
1119    for roll in rolls:
1120        if np.any(roll):        #avoid [0,0]
1121            spotMask = ma.array(spotMask,mask=(spotMask-rollImage(Image,roll)<0.))
1122    mags = spotMask[spotMask.nonzero()]
1123    indx = np.transpose(spotMask.nonzero())
1124    nx,ny = Image.shape
1125    jndx = []
1126    for [ind,mag] in zip(indx,mags):
1127        if (1 < ind[0] < nx-2) and (1 < ind[1] < ny-2):
1128            cent = np.zeros((3,3))
1129            cent[1,1] = mag
1130            msk = np.array(Image[ind[0]-2:ind[0]+3,ind[1]-2:ind[1]+3])
1131            msk = ma.array(msk,mask=(msk==mag))
1132            if mag > 1.5*ma.mean(msk):
1133                jndx.append([ind[1]+.5,ind[0]+.5])
1134    print 'Spots found: ',len(jndx)
1135    jndx = np.array(jndx)
1136    peaks = jndx*pixelSize/1000.
1137    tth = GetTth(peaks.T[0],peaks.T[1],Controls)
1138    histth,bins = np.histogram(tth,2500)
1139    for shft in [-1,1]:
1140        histmsk = ma.array(histth,mask=-(histth-np.roll(histth,shft)<2.))
1141    histmsk = ma.array(histmsk,mask=(histmsk>5))
1142    binmsk = np.array(ma.nonzero(histmsk))
1143    digts = np.digitize(tth,bins)-1
1144    PeaksList = []
1145    for digt,peak in zip(digts,peaks):
1146        if digt in binmsk:
1147            PeaksList.append(peak)
1148    #should be able to filter out spotty Bragg rings here
1149    PeaksList = np.array(PeaksList)
1150#    Peakarray = np.vstack((tth,peaks.T)).T
1151#    Peakarray = np.array(G2mth.sortArray(Peakarray,0))  #now in 2theta
1152#    if peaks.shape[0] > 100:
1153#        txt = 'More than 100 spots found: %d. Are rings spotty?'%(len(jndx))
1154#        return txt
1155    Points = np.ones((PeaksList.shape[0],3))
1156#    Points[:,:2] = Peakarray[:,1:]
1157    Points[:,:2] = PeaksList
1158    Masks['Points'] = list(Points)
1159    return None
1160                   
1161   
1162       
Note: See TracBrowser for help on using the repository browser.