Changeset 47


Ignore:
Timestamp:
Apr 5, 2010 11:17:25 AM (12 years ago)
Author:
vondreel
Message:

eliminate fitellipse.for
fitting with scipy.optimize.leastsq
some mods to facilitate this - numpy arrays needed

Location:
trunk
Files:
1 deleted
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIcomp.py

    r46 r47  
    3131    # this error is non-recoverable, so just quit
    3232    exit()
    33 import fitellipse as fte
    3433import GSASIIplot as G2plt
    3534
     
    12081207       
    12091208def FitEllipse(ring):
    1210     N = len(ring)
    1211     A = np.zeros(shape=(N,5),order='F',dtype=np.float32)
    1212     B = np.zeros(shape=N,dtype=np.float32)
    1213     X = np.zeros(shape=N,dtype=np.float32)
    1214     Y = np.zeros(shape=N,dtype=np.float32)
    1215     for i,(x,y) in enumerate(ring):
    1216         X[i] = x
    1217         Y[i] = y
    1218     B,Krank,Rnorm = fte.fitqdr(N,A,B,X,Y)
     1209    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]
    12191218   
    12201219    ell = [1.+B[0],B[1],1.-B[0],B[2],B[3],B[4]]
     
    12441243            phi += 180.
    12451244    return cent,phi,[sr1,sr2]
    1246        
    1247 def FitCircle(ring):
    1248     N = len(ring)
    1249     A = np.zeros(shape=(N,3),order='F',dtype=np.float32)
    1250     B = np.zeros(shape=N,dtype=np.float32)
    1251     X = np.zeros(shape=N,dtype=np.float32)
    1252     Y = np.zeros(shape=N,dtype=np.float32)
    1253     for i,(x,y) in enumerate(ring):
    1254         X[i] = x
    1255         Y[i] = y
    1256     B,Krank,Rnorm = fte.fitcrc(N,A,B,X,Y)
    1257     cent = (-0.5*B[0],-0.5*B[1])
    1258     R2 = cent[0]**2+cent[1]**2-B[2]
    1259     if R2 > 0.0:
    1260         radius = math.sqrt(R2)
    1261     else:
    1262         return 0
    1263     return cent,radius
    1264    
     1245           
    12651246def ImageLocalMax(image,w,Xpix,Ypix):
    12661247    w2 = w*2
     
    13571338    return GetTthDspAzm(x,y,data)[1]
    13581339       
    1359 def GetDetector(data):
    1360     return
    1361 
    13621340def ImageCompress(image,scale):
    13631341    if scale == 1:
     
    14531431        if Ring:
    14541432            numZ = len(Ring)
    1455             data['rings'].append(np.array(Ring))
     1433            data['rings'].append(np.column_stack((np.array(Ring),dsp*np.ones(numZ))))
    14561434            elcent,phi,radii = ellipse = FitEllipse(Ring)
    14571435            if abs(phi) > 45. and phi < 0.:
     
    15121490    t2 = d2/len2
    15131491    if len2 > len1:
    1514         print 'len2 > len1'
    15151492        if -135. < atan2d(t2[1],t2[0]) < 45.:
    15161493            Zsign = -1
    15171494    else:
    1518         print 'len2 < len1'
    15191495        if -135. < atan2d(t1[1],t1[0]) < 45.:
    15201496            Zsign = -1
     
    15221498    tilt = data['tilt'] = Zsign*tiltSum/Zsum
    15231499    phi = data['rotation'] = phiSum/Zsum
     1500    rings = np.concatenate((data['rings']),axis=0)
     1501    print wave,dist,cent,phi,tilt
     1502   
     1503   
    15241504    N = len(data['ellipses'])
    15251505    for H in HKL[:N]:
  • trunk/GSASIIplot.py

    r40 r47  
    363363        imgPlot.text(xring,yring,'+',color='b',ha='center',va='center',picker=3)
    364364    if Data['setRings']:
    365         for ring in Data['rings']:
    366             for xring,yring in ring:
    367                 imgPlot.text(xring,yring,'+',ha='center',va='center')           
     365        rings = Data['rings']
     366        for xring,yring,tth in rings:
     367            imgPlot.text(xring,yring,'+',ha='center',va='center')           
    368368    for ellipse in Data['ellipses']:
    369369        cent,phi,[width,height],col = ellipse
Note: See TracChangeset for help on using the changeset viewer.