Changeset 26


Ignore:
Timestamp:
Feb 8, 2010 9:52:06 AM (12 years ago)
Author:
vondreel
Message:
 
File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIcomp.py

    r24 r26  
    1919    # this error is non-recoverable, so just quit
    2020    exit()
     21import fitellipse as fte
    2122
    2223# trig functions in degrees
     
    2425asind = lambda x: 180.*math.asin(x)/math.pi
    2526tand = lambda x: math.tan(x*math.pi/180.)
     27atand = lambda x: 180.*math.atan(x)/math.pi
    2628cosd = lambda x: math.cos(x*math.pi/180.)
    2729acosd = lambda x: 180.*math.acos(x)/math.pi
     
    11831185    else:
    11841186        return False,0,0
    1185 
    1186 def ImageCalibrate(data,image,PlotImage):
     1187       
     1188def FitEllipse(ring):
     1189    N = len(ring)
     1190    A = np.zeros(shape=(N,5),order='F',dtype=np.float32)
     1191    B = np.zeros(shape=N,dtype=np.float32)
     1192    X = np.zeros(shape=N,dtype=np.float32)
     1193    Y = np.zeros(shape=N,dtype=np.float32)
     1194    for i,(x,y) in enumerate(ring):
     1195        X[i] = x
     1196        Y[i] = y
     1197    B,Krank,Rnorm = fte.fitqdr(N,A,B,X,Y)
     1198   
     1199    ell = [1.+B[0],B[1],1.-B[0],B[2],B[3],B[4]]
     1200    # ell is [A,B,C,D,E,F] for Ax^2+Bxy+Cy^2+Dx+Ey+F = 0
     1201    det = 4.*ell[0]*ell[2]-ell[1]**2
     1202    if det < 0.:
     1203        print 'hyperbola!'
     1204        return 0
     1205    elif det == 0.:
     1206        print 'parabola!'
     1207        return 0
     1208    cent = [(ell[1]*ell[4]-2.0*ell[2]*ell[3])/det, \
     1209        (ell[1]*ell[3]-2.0*ell[0]*ell[4])/det]
     1210    phi = 0.5*atand(ell[1]/(ell[0]-ell[2]))
     1211   
     1212    a3 = 0.5*(ell[0]+ell[2]+(ell[0]-ell[2])/cosd(2.0*phi))
     1213    b3 = ell[0]+ell[2]-a3
     1214    f3 = ell[0]*cent[0]**2+ell[2]*cent[1]**2+ell[1]*cent[0]*cent[1]-ell[5]
     1215    if f3/a3 < 0 or f3/b3 < 0:
     1216        return 0
     1217    sr1 = math.sqrt(f3/a3)
     1218    sr2 = math.sqrt(f3/b3)
     1219    return ell,cent,phi,[sr1,sr2]
     1220       
     1221def FitCircle(ring):
     1222    N = len(ring)
     1223    A = np.zeros(shape=(N,3),order='F',dtype=np.float32)
     1224    B = np.zeros(shape=N,dtype=np.float32)
     1225    X = np.zeros(shape=N,dtype=np.float32)
     1226    Y = np.zeros(shape=N,dtype=np.float32)
     1227    for i,(x,y) in enumerate(ring):
     1228        X[i] = x
     1229        Y[i] = y
     1230    B,Krank,Rnorm = fte.fitcrc(N,A,B,X,Y)
     1231    cent = (-0.5*B[0],-0.5*B[1])
     1232    R2 = cent[0]**2+cent[1]**2-B[2]
     1233    if R2 > 0.0:
     1234        radius = math.sqrt(R2)
     1235    else:
     1236        return 0
     1237    return cent,radius
     1238   
     1239def ImageLocalMax(image,w,Xpos,Ypos):
     1240    w2 = w*2
     1241    Z = image[Ypos-w:Ypos+w,Xpos-w:Xpos+w]
     1242    Zmax = np.argmax(Z)
     1243    Zmin = np.argmin(Z)
     1244    Xpos += Zmax%w2-w
     1245    Ypos += Zmax/w2-w
     1246    return Xpos,Ypos,np.ravel(Z)[Zmax],np.ravel(Z)[Zmin]
     1247   
     1248def makeRing(cent,phi,semiradii,scalex,scaley,imScale,image):
     1249    cphi = cosd(phi)
     1250    sphi = sind(phi)
     1251    ring = []
     1252    for a in range(0,360,4):
     1253        x = semiradii[0]*cosd(a)
     1254        y = semiradii[1]*sind(a)
     1255        X = (cphi*x-sphi*y+cent[0])*scalex*imScale
     1256        Y = (sphi*x+cphi*y+cent[1])*scaley*imScale
     1257        if (0 < X < len(image)) and (0 < Y < len(image)):       
     1258            X,Y,I,J = ImageLocalMax(image,20,X,Y)
     1259            if I and J:
     1260                X /= scalex*imScale
     1261                Y /= scaley*imScale
     1262                ring.append([X,Y])   
     1263    return ring
     1264       
     1265def ImageCalibrate(self,data):
     1266    import ImageCalibrants as calFile
    11871267    print 'image calibrate'
     1268    if not data['calibrant']:
     1269        print 'no calibration material selected'
     1270        return
     1271    Bravais,cell = calFile.Calibrants[data['calibrant']]
    11881272    ring = data['ring']
     1273    data['ellipses'] = ellipses = []
     1274    scalex = data['scalex']             # = 1000./(pixelSize[0]*self.imScale)
     1275    scaley = data['scaley']
    11891276    if len(ring) < 5:
    11901277        print 'not enough inner ring points for ellipse'
    11911278        return
    1192     A = np.zeros(shape=(5,len(ring)))
    1193     B = np.zeros(shape=len(ring))
    1194     for i,(x,y) in enumerate(ring):
    1195         A[0,i] = x**2-y**2
    1196         A[1,i] = x*y
    1197         A[2,i] = x
    1198         A[3,i] = y
    1199         A[4,i] = 1.
    1200         B[i] = -(x**2+y**2)
     1279    outB = FitEllipse(ring)
     1280    if outB:
     1281        ell,cent,phi,semiradii = outB
     1282        print 'start ellipse coeff.',ell
     1283        print 'start:',cent,phi,semiradii
     1284        data['ellipses'].append([cent,phi,semiradii])
     1285        self.PlotImage()
     1286    else:
     1287        return False
     1288    ring = makeRing(cent,phi,semiradii,scalex,scaley,self.imScale,self.ImageZ)
     1289    ell,cent,phi,semiradii = FitEllipse(ring)
     1290    print '1st ring ellipse coeff.',ell
     1291    print '1st ring:',cent,phi,semiradii
     1292    data['ellipses'].append([cent,phi,semiradii])
     1293    data['center'] = cent           #not right!! (but useful for now)
     1294    self.PlotImage()
    12011295       
    1202     for a in A: print a
    1203     print B
     1296   
     1297   
     1298       
     1299       
     1300    return True
     1301       
     1302       
Note: See TracChangeset for help on using the changeset viewer.