Changeset 48


Ignore:
Timestamp:
Apr 12, 2010 1:56:43 PM (13 years ago)
Author:
vondreel
Message:

finished ellipse fitting
finished integration range display

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASII.py

    r45 r48  
    409409                        Data['calibrant'] = ''
    410410                        Data['IOtth'] = [2.0,5.0]
    411                         Data['LRazimuth'] = [-45,45]
     411                        Data['LRazimuth'] = [-135,-45]
    412412                        Data['outChannels'] = 2500
    413413                        Data['fullIntegrate'] = False

  • trunk/GSASIIcomp.py

    r47 r48  
    12061206        return False,0,0
    12071207       
     1208def FitRing(ring):
     1209    Err,parms = FitCircle(ring)
     1210    Err /= len(ring)
     1211    print 'circle error:','%8f'%(Err)
     1212    if Err > 20000.:
     1213        eparms = FitEllipse(ring)
     1214        if eparms:
     1215            parms = eparms
     1216    return parms
     1217       
     1218def FitCircle(ring):
     1219    import numpy.linalg as nl
     1220   
     1221    def makeParmsCircle(B):
     1222        cent = [-B[0]/2,-B[1]/2]
     1223        phi = 0.
     1224        sr1 = sr2 = math.sqrt(cent[0]**2+cent[1]**2-B[2])
     1225        return cent,phi,[sr1,sr2]
     1226       
     1227    ring = np.array(ring)
     1228    x = np.asarray(ring.T[0])
     1229    y = np.asarray(ring.T[1])
     1230   
     1231    M = np.array((x,y,np.ones_like(x)))
     1232    B = np.array(-(x**2+y**2))
     1233    result = nl.lstsq(M.T,B)
     1234    return result[1],makeParmsCircle(result[0])
     1235       
    12081236def FitEllipse(ring):
     1237    import numpy.linalg as nl
     1238           
     1239    def makeParmsEllipse(B):
     1240        det = 4.*(1.-B[0]**2)-B[1]**2
     1241        if det < 0.:
     1242            print 'hyperbola!'
     1243            return 0
     1244        elif det == 0.:
     1245            print 'parabola!'
     1246            return 0
     1247        cent = [(B[1]*B[3]-2.*(1.-B[0])*B[2])/det, \
     1248            (B[1]*B[2]-2.*(1.+B[0])*B[3])/det]
     1249        phi = 0.5*atand(0.5*B[1]/B[0])
     1250       
     1251        a = (1.+B[0])/cosd(2*phi)
     1252        b = 2.-a
     1253        f = (1.+B[0])*cent[0]**2+(1.-B[0])*cent[1]**2+B[1]*cent[0]*cent[1]-B[4]
     1254        if f/a < 0 or f/b < 0:
     1255            return 0
     1256        sr1 = math.sqrt(f/a)
     1257        sr2 = math.sqrt(f/b)
     1258        if sr1 > sr2:
     1259            sr1,sr2 = SwapXY(sr1,sr2)
     1260            phi -= 90.
     1261            if phi < -90.:
     1262                phi += 180.
     1263        return cent,phi,[sr1,sr2]
     1264               
     1265    ring = np.array(ring)
     1266    x = np.asarray(ring.T[0])
     1267    y = np.asarray(ring.T[1])
     1268    M = np.array((x**2-y**2,x*y,x,y,np.ones_like(x)))
     1269    B = np.array(-(x**2+y**2))
     1270    result = nl.lstsq(M.T,B)
     1271    return makeParmsEllipse(result[0])
     1272   
     1273def FitDetector(rings,p0,wave):
    12091274    from scipy.optimize import leastsq
    1210     def ellipseCalc(A,xy):
    1211         x = xy[0]
    1212         y = xy[1]
    1213         return (1.+A[0])*x**2+(1.-A[0])*y**2+A[1]*x*y+A[2]*x+A[3]*y+A[4]
    1214     ring = np.array(ring)
    1215     p0 = np.zeros(shape=5)
    1216     result = leastsq(ellipseCalc,p0,args=(ring.T,))
    1217     B = result[0]
    1218    
    1219     ell = [1.+B[0],B[1],1.-B[0],B[2],B[3],B[4]]
    1220     # ell is [A,B,C,D,E,F] for Ax^2+Bxy+Cy^2+Dx+Ey+F = 0
    1221     det = 4.*ell[0]*ell[2]-ell[1]**2
    1222     if det < 0.:
    1223         print 'hyperbola!'
    1224         return 0
    1225     elif det == 0.:
    1226         print 'parabola!'
    1227         return 0
    1228     cent = [(ell[1]*ell[4]-2.0*ell[2]*ell[3])/det, \
    1229         (ell[1]*ell[3]-2.0*ell[0]*ell[4])/det]
    1230     phi = 0.5*atand(ell[1]/(ell[0]-ell[2]+1e-32))
    1231    
    1232     a3 = 0.5*(ell[0]+ell[2]+(ell[0]-ell[2])/cosd(2.0*phi))
    1233     b3 = ell[0]+ell[2]-a3
    1234     f3 = ell[0]*cent[0]**2+ell[2]*cent[1]**2+ell[1]*cent[0]*cent[1]-ell[5]
    1235     if f3/a3 < 0 or f3/b3 < 0:
    1236         return 0
    1237     sr1 = math.sqrt(f3/a3)
    1238     sr2 = math.sqrt(f3/b3)
    1239     if sr1 > sr2:
    1240         sr1,sr2 = SwapXY(sr1,sr2)
    1241         phi -= 90.
    1242         if phi < -90.:
    1243             phi += 180.
    1244     return cent,phi,[sr1,sr2]
     1275    def ellipseCalc(B,xyd,wave):
     1276        sind = lambda x: np.sin(x*math.pi/180.)
     1277        asind = lambda x: 180.*np.arcsin(x)/math.pi
     1278        cosd = lambda x: np.cos(x*math.pi/180.)
     1279        tand = lambda x: np.tan(x*math.pi/180.)
     1280        x = xyd[0]
     1281        y = xyd[1]
     1282        dsp = xyd[2]
     1283        dist,x0,y0,phi,tilt = B
     1284        tth = 2.0*asind(wave/(2.*dsp))
     1285        ttth = tand(tth)
     1286        radius = dist*ttth
     1287        stth = sind(tth)
     1288        cosb = cosd(tilt)
     1289        R1 = dist*stth*cosd(tth)*cosb/(cosb**2-stth**2)
     1290        R0 = np.sqrt(R1*radius*cosb)
     1291        zdis = R1*ttth*tand(tilt)
     1292        X = x-x0+zdis*sind(phi)
     1293        Y = y-y0-zdis*cosd(phi)
     1294        XR = X*cosd(phi)-Y*sind(phi)
     1295        YR = X*sind(phi)+Y*cosd(phi)
     1296        return (XR/R0)**2+(YR/R1)**2-1
     1297    result = leastsq(ellipseCalc,p0,args=(rings.T,wave))
     1298    return result[0]
     1299
    12451300           
    12461301def ImageLocalMax(image,w,Xpix,Ypix):
    … …  
    12571312        return 0,0,0,0
    12581313   
    1259 def makeRing(ellipse,pix,reject,scalex,scaley,image):
     1314def makeRing(dsp,ellipse,pix,reject,scalex,scaley,image):
    12601315    cent,phi,radii = ellipse
    12611316    cphi = cosd(phi)
    … …  
    12711326            X /= scalex                         #convert to mm
    12721327            Y /= scaley
    1273             ring.append([X,Y])
     1328            ring.append([X,Y,dsp])
    12741329    if len(ring) < 45:             #want more than 1/4 of a circle
    12751330        return []
    … …  
    13021357    stth = sind(tth)
    13031358    ctth = cosd(tth)
    1304     cosb = sind(tilt)
    1305     sinb = math.sqrt(1.-cosb**2)
    1306     sinp = sind(phi)
    1307     cosp = cosd(phi)
     1359    cosb = cosd(tilt)
    13081360    radius = dist*ttth
    1309     radii[1] = dist*stth*ctth*sinb/(sinb**2-stth**2)
     1361    radii[1] = dist*stth*ctth*cosb/(cosb**2-stth**2)
    13101362    if radii[1] > 0:
    1311         radii[0] = math.sqrt(radii[1]*radius*sinb)
    1312         zdis = radii[1]*ttth*cosb/sinb
    1313         elcent = [cent[0]-zdis*sinp,cent[1]+zdis*cosp]
     1363        radii[0] = math.sqrt(radii[1]*radius*cosb)
     1364        zdis = radii[1]*ttth*tand(tilt)
     1365        elcent = [cent[0]-zdis*sind(phi),cent[1]+zdis*cosd(phi)]
    13141366        return elcent,phi,radii
    13151367    else:
    13161368        return False
    1317    
     1369       
     1370def GetDetectorXY(dsp,azm,data):
     1371    from scipy.optimize import fsolve
     1372    def func(xy,*args):
     1373       azm,phi,R0,R1,A,B = args
     1374       cp = cosd(phi)
     1375       sp = sind(phi)
     1376       x,y = xy
     1377       out = []
     1378       out.append(y-x*tand(azm))
     1379       out.append(R0**2*((x+A)*sp-(y+B)*cp)**2+R1**2*((x+A)*cp+(y+B)*sp)**2-(R0*R1)**2)
     1380       return out
     1381    elcent,phi,radii = GetEllipse(dsp,data)
     1382    cent = data['center']
     1383    tilt = data['tilt']
     1384    phi = data['rotation']
     1385    wave = data['wavelength']
     1386    dist = data['distance']
     1387    tth = 2.0*asind(wave/(2.*dsp))
     1388    ttth = tand(tth)
     1389    radius = dist*ttth
     1390    stth = sind(tth)
     1391    cosb = cosd(tilt)
     1392    R1 = dist*stth*cosd(tth)*cosb/(cosb**2-stth**2)
     1393    R0 = math.sqrt(R1*radius*cosb)
     1394    zdis = R1*ttth*tand(tilt)
     1395    A = zdis*sind(phi)
     1396    B = -zdis*cosd(phi)
     1397    xy0 = [radius*cosd(azm),radius*sind(azm)]
     1398    xy = fsolve(func,xy0,args=(azm,phi,R0,R1,A,B))+cent
     1399    return xy
     1400                   
    13181401def GetTthDspAzm(x,y,data):
     1402    #make numpy array compliant
     1403    atand = lambda x: 180.*np.arctan(x)/np.pi
     1404    sind = lambda x: np.sin(x*np.pi/180.)
     1405    atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
    13191406    wave = data['wavelength']
    13201407    dist = data['distance']
    … …  
    13271414    X = np.sum(X*makeMat(-phi,2),axis=1)
    13281415    X = np.sum(X*makeMat(-tilt,0),axis=1)
    1329     tth = atand(math.sqrt(dx**2+dy**2-X[2]**2)/(dist-X[2]))
     1416    tth = atand(np.sqrt(dx**2+dy**2-X[2]**2)/(dist-X[2]))
    13301417    dsp = wave/(2.*sind(tth/2.))
    13311418    azm = atan2d(dy,dx)   
    … …  
    13591446    #fit start points on inner ring
    13601447    data['ellipses'] = []
    1361     outB = FitEllipse(ring)
    1362     if outB:
    1363         ellipse = outB
    1364         print 'start:',ellipse
     1448    outE = FitRing(ring)
     1449    if outE:
     1450        print 'start ellipse:',outE
     1451        ellipse = outE
    13651452    else:
    13661453        return False
    13671454       
    13681455    #setup 180 points on that ring for "good" fit
    1369     Ring = makeRing(ellipse,20,cutoff,scalex,scaley,self.ImageZ)
     1456    Ring = makeRing(1.0,ellipse,20,cutoff,scalex,scaley,self.ImageZ)
    13701457    if Ring:
    1371         ellipse = FitEllipse(Ring)
    1372         Ring = makeRing(ellipse,20,cutoff,scalex,scaley,self.ImageZ)    #do again
    1373         ellipse = FitEllipse(Ring)
     1458        ellipse = FitRing(Ring)
     1459        Ring = makeRing(1.0,ellipse,20,cutoff,scalex,scaley,self.ImageZ)    #do again
     1460        ellipse = FitRing(Ring)
    13741461    else:
    13751462        print '1st ring not sufficiently complete to proceed'
    … …  
    14001487    radius = dist*tand(tth)
    14011488    zdis,cosB = calcZdisCosB(radius,tth,radii)
    1402     sinp = sind(ellipse[1])
    1403     cosp = cosd(ellipse[1])
    14041489    cent1 = []
    14051490    cent2 = []
    … …  
    14281513        elcent = [cent[0]+zsinp,cent[1]-zcosp]
    14291514        ratio = radii[1]/radii[0]
    1430         Ring = makeRing(ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
     1515        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
    14311516        if Ring:
    14321517            numZ = len(Ring)
    1433             data['rings'].append(np.column_stack((np.array(Ring),dsp*np.ones(numZ))))
    1434             elcent,phi,radii = ellipse = FitEllipse(Ring)
     1518            data['rings'].append(np.array(Ring))
     1519            ellipse = FitRing(Ring)
     1520            elcent,phi,radii = ellipse               
    14351521            if abs(phi) > 45. and phi < 0.:
    14361522                phi += 180.
    … …  
    14991585    phi = data['rotation'] = phiSum/Zsum
    15001586    rings = np.concatenate((data['rings']),axis=0)
    1501     print wave,dist,cent,phi,tilt
    1502    
    1503    
     1587    p0 = [dist,cent[0],cent[1],phi,tilt]
     1588    result = FitDetector(rings,p0,wave)
     1589    data['distance'] = result[0]
     1590    data['center'] = result[1:3]
     1591    data['rotation'] = np.mod(result[3],360.0)
     1592    data['tilt'] = result[4]
    15041593    N = len(data['ellipses'])
     1594    data['ellipses'] = []           #clear away individual ellipse fits
    15051595    for H in HKL[:N]:
    15061596        ellipse = GetEllipse(H[3],data)
    … …  
    15091599       
    15101600    return True
    1511    
    1512    
     1601       
    15131602def test():
    15141603    cell = [5,5,5,90,90,90]

  • trunk/GSASIIgrid.py

    r43 r48  
    11641164        if data['fullIntegrate']:
    11651165            data['fullIntegrate'] = False
    1166             data['LRazimuth'] = [-45,45]
    1167             self.Lazim.SetValue("%6d" % (-45))           
    1168             self.Razim.SetValue("%6d" % (45))           
     1166            self.Lazim.SetEditable(True)           
     1167            self.Razim.SetEditable(True)           
    11691168        else:
    11701169            data['fullIntegrate'] = True
    1171             data['LRazimuth'] = [0,360]
    1172             self.Lazim.SetValue("%6d" % (0))           
    1173             self.Razim.SetValue("%6d" % (360))           
     1170            self.Lazim.SetEditable(False)           
     1171            self.Razim.SetEditable(False)           
    11741172        G2plt.PlotImage(self)
    11751173       

  • trunk/GSASIIplot.py

    r47 r48  
    186186                    ypos = event.ydata-ycent
    187187                    if 'line3' in  str(item) or 'line4' in str(item) and not Data['fullIntegrate']:
    188                         ang = int(atan2d(-ypos,xpos))
     188                        ang = int(atan2d(xpos,ypos))
    189189                        imgPage.canvas.SetToolTipString('%6d deg'%(ang))
    190                     elif 'inner' in  str(item.get_gid()) or 'outer' in str(item.get_gid()):
     190                    elif 'line1' in  str(item) or 'line2' in str(item):
    191191                        tth = G2cmp.GetTth(event.xdata,event.ydata,Data)
    192192                        imgPage.canvas.SetToolTipString('%8.3fdeg'%(tth))                           
    … …  
    335335        IOtth = Data['IOtth']
    336336        wave = Data['wavelength']
    337         dsp = wave/(2.0*sind(IOtth[0]/2.0))
    338         ell0 = G2cmp.GetEllipse(dsp,Data)           #=False if dsp didn't yield an ellipse (ugh! a parabola or a hyperbola)
    339         dsp = wave/(2.0*sind(IOtth[1]/2.0))
    340         ell1 = G2cmp.GetEllipse(dsp,Data)           #Ditto & more likely for outer ellipse
     337        dspI = wave/(2.0*sind(IOtth[0]/2.0))
     338        ellI = G2cmp.GetEllipse(dspI,Data)           #=False if dsp didn't yield an ellipse (ugh! a parabola or a hyperbola)
     339        dspO = wave/(2.0*sind(IOtth[1]/2.0))
     340        ellO = G2cmp.GetEllipse(dspO,Data)           #Ditto & more likely for outer ellipse
    341341        if Data['fullIntegrate']:
    342             arcxI = arcyI = np.array(range(0,361))
    343             arcxO = arcyO = np.array(range(0,361))
     342            Azm = np.array(range(0,361))
    344343        else:
    345             arcxI = arcyI = np.array(range(LRAzim[0],LRAzim[1]+1))
    346             arcxO = arcyO = np.array(range(LRAzim[0],LRAzim[1]+1))
    347         if ell0:
    348             cent0,phi0,radii0 = ell0
    349             radii0 = np.sum(G2cmp.makeMat(phi0,2).T*np.array(radii0+[0.,]),axis=1)
    350             arcxI = np.cos(arcxI*math.pi/180.)*radii0[0]+cent0[0]
    351             arcyI = np.sin(arcyI*math.pi/180.)*radii0[1]+cent0[1]
     344            Azm = np.array(range(LRAzim[0],LRAzim[1]+1))
     345        if ellI:
     346            xyI = []
     347            for azm in Azm:
     348                xyI.append(G2cmp.GetDetectorXY(dspI,azm,Data))
     349            xyI = np.array(xyI)
     350            arcxI,arcyI = xyI.T
    352351            imgPlot.plot(arcxI,arcyI,picker=3)
    353         if ell1:
    354             cent1,phi1,radii1 = ell1
    355             radii1 = np.sum(G2cmp.makeMat(phi1,2).T*np.array(radii1+[0.,]),axis=1)
    356             arcxO = np.cos(arcxO*math.pi/180.)*radii1[0]+cent1[0]
    357             arcyO = np.sin(arcyO*math.pi/180.)*radii1[1]+cent1[1]
    358             imgPlot.plot(arcxO,arcyO,picker=3)       
    359         if ell0 and ell1 and not Data['fullIntegrate']:
     352        if ellO:
     353            xyO = []
     354            for azm in Azm:
     355                xyO.append(G2cmp.GetDetectorXY(dspO,azm,Data))
     356            xyO = np.array(xyO)
     357            arcxO,arcyO = xyO.T
     358            imgPlot.plot(arcxO,arcyO,picker=3)
     359        if ellO and ellI and not Data['fullIntegrate']:
    360360            imgPlot.plot([arcxI[0],arcxO[0]],[arcyI[0],arcyO[0]],picker=3)
    361361            imgPlot.plot([arcxI[-1],arcxO[-1]],[arcyI[-1],arcyO[-1]],picker=3)
    … …  
    363363        imgPlot.text(xring,yring,'+',color='b',ha='center',va='center',picker=3)
    364364    if Data['setRings']:
    365         rings = Data['rings']
    366         for xring,yring,tth in rings:
     365        rings = np.concatenate((Data['rings']),axis=0)
     366        for xring,yring,dsp in rings:
    367367            imgPlot.text(xring,yring,'+',ha='center',va='center')           
    368368    for ellipse in Data['ellipses']:
Note: See TracChangeset for help on using the changeset viewer.