Changeset 29


Ignore:
Timestamp:
Feb 11, 2010 5:46:51 PM (12 years ago)
Author:
vondreel
Message:

next stage of image calibration

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASII.py

    r25 r29  
    398398        dlg = wx.FileDialog(self, 'Choose image file', '.', '', \
    399399            'MAR345 (*.mar3450)|*.mar3450|ADSC Image (*.img)|*.img \
    400             |Perkin-Elmer TIF (*.tif)|*.tif|GE Image sum (*.sum)|*.sum|All files (*.*)|*.*',wx.OPEN)
     400            |Perkin-Elmer TIF (*.tif)|*.tif \
     401            |GE Image sum (*.sum)|*.sum|GE Image avg (*.avg) \
     402            |*.avg|All files (*.*)|*.*',wx.OPEN)
    401403        if self.dirname:
    402404            dlg.SetDirectory(self.dirname)
     
    414416                elif ext == '.mar3450':
    415417                    Comments,Data,Size,Image = G2IO.GetMAR345Data(self.imagefile)
    416                 elif ext == '.sum':
     418                elif ext in ['.sum','.avg']:
    417419                    Comments,Data,Size,Image = G2IO.GetGEsumData(self.imagefile)
    418420                if Comments:
     
    432434                        Data['showLines'] = False
    433435                        Data['ring'] = []
     436                        Data['rings'] = []
    434437                        Data['ellipses'] = []
    435438                        Data['masks'] = []
     
    13031306                    xpos = int(event.xdata)*self.imScale
    13041307                    ypos = int(event.ydata)*self.imScale
    1305                     self.pdplot.canvas.SetToolTipString('%6d'%(self.ImageZ[ypos][xpos]))
     1308                    if (0 <= xpos <= size) and (0 <= ypos <= size):
     1309                        self.pdplot.canvas.SetToolTipString('%6d'%(self.ImageZ[ypos][xpos]))
    13061310
    13071311        def OnImPlotKeyPress(event):
     
    14561460            yring *= scaley
    14571461            ax.text(xring,yring,'+',ha='center',va='center',picker=3)
     1462#        for ring in Data['rings']:
     1463#            for xring,yring in ring:
     1464#                xring *= scalex
     1465#                yring *= scaley
     1466#                ax.text(xring,yring,'+',ha='center',va='center')           
    14581467        for ellipse in Data['ellipses']:
    14591468            cent,phi,[width,height] = ellipse
  • trunk/GSASIIIO.py

    r28 r29  
    399399def GetGEsumData(filename):
    400400    import array as ar
    401     print 'Read GE sum file: ',filename
     401    print 'Read GE sum file: ',filename   
    402402    File = open(filename,'rb')
    403403    size = 2048
     404    if '.sum' in filename:
     405        head = ['GE detector sum data from APS 1-ID',]
     406    if '.avg' in filename:
     407#        image = np.zeros(shape=(size,size),dtype=np.int16)
     408        head = ['GE detector avg data from APS 1-ID',]
    404409    image = np.zeros(shape=(size,size),dtype=np.int32)
    405     head = ['GE detector sum data',]
    406410    row = 0
    407411    pos = 0
    408412    while row < size:
    409413        File.seek(pos)
    410         line = ar.array('f',File.read(4*size))
     414        if '.sum' in filename:
     415            line = ar.array('f',File.read(4*size))
     416            pos += 4*size
     417        elif '.avg' in filename:
     418            line = ar.array('H',File.read(2*size))
     419            pos += 2*size
    411420        image[row] = np.asarray(line)
    412421        row += 1
    413         pos += 4*size
    414     data = {'pixelSize':(200,200),'wavelength':0.10,'distance':250.0,'center':[204.8,204.8]} 
     422    data = {'pixelSize':(200,200),'wavelength':0.15,'distance':250.0,'center':[204.8,204.8]} 
    415423    return head,data,size,image
    416424    File.close()   
  • trunk/GSASIIcomp.py

    r26 r29  
    12151215    if f3/a3 < 0 or f3/b3 < 0:
    12161216        return 0
    1217     sr1 = math.sqrt(f3/a3)
    1218     sr2 = math.sqrt(f3/b3)
    1219     return ell,cent,phi,[sr1,sr2]
     1217    if phi > 0.:       
     1218        sr1 = math.sqrt(f3/a3)
     1219        sr2 = math.sqrt(f3/b3)
     1220    else:
     1221        sr1 = math.sqrt(f3/a3)
     1222        sr2 = math.sqrt(f3/b3)
     1223        phi *= -1.       
     1224    return cent,phi,[sr1,sr2]
    12201225       
    12211226def FitCircle(ring):
     
    12391244def ImageLocalMax(image,w,Xpos,Ypos):
    12401245    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    
    1248 def makeRing(cent,phi,semiradii,scalex,scaley,imScale,image):
     1246    size = len(image)
     1247    if (w < Xpos < size-w) and (w < Ypos < size-w):
     1248        Z = image[Ypos-w:Ypos+w,Xpos-w:Xpos+w]
     1249        Zmax = np.argmax(Z)
     1250        Zmin = np.argmin(Z)
     1251        Xpos += Zmax%w2-w
     1252        Ypos += Zmax/w2-w
     1253        return Xpos,Ypos,np.ravel(Z)[Zmax],np.ravel(Z)[Zmin]
     1254    else:
     1255        return 0,0,0,0
     1256   
     1257def makeRing(ellipse,scalex,scaley,imScale,image):
     1258    cent,phi,semiradii = ellipse
    12491259    cphi = cosd(phi)
    12501260    sphi = sind(phi)
    12511261    ring = []
    1252     for a in range(0,360,4):
     1262    for a in range(0,360):
    12531263        x = semiradii[0]*cosd(a)
    12541264        y = semiradii[1]*sind(a)
    12551265        X = (cphi*x-sphi*y+cent[0])*scalex*imScale
    12561266        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])   
     1267        X,Y,I,J = ImageLocalMax(image,20,X,Y)     
     1268        if I and J:
     1269            X /= scalex*imScale
     1270            Y /= scaley*imScale
     1271            ring.append([X,Y])
     1272    if len(ring) < 75:             #want more than 1/3 of a circle
     1273        return []
    12631274    return ring
    12641275       
    12651276def ImageCalibrate(self,data):
     1277    import copy
    12661278    import ImageCalibrants as calFile
    12671279    print 'image calibrate'
    1268     if not data['calibrant']:
    1269         print 'no calibration material selected'
    1270         return
    1271     Bravais,cell = calFile.Calibrants[data['calibrant']]
    12721280    ring = data['ring']
    1273     data['ellipses'] = ellipses = []
    12741281    scalex = data['scalex']             # = 1000./(pixelSize[0]*self.imScale)
    12751282    scaley = data['scaley']
    12761283    if len(ring) < 5:
    12771284        print 'not enough inner ring points for ellipse'
    1278         return
     1285        return False
     1286       
     1287    #fit start points on inner ring
     1288    data['ellipses'] = []
    12791289    outB = FitEllipse(ring)
    12801290    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()
     1291        ellipse = outB
     1292        print 'start:',ellipse
    12861293    else:
    12871294        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)
     1295       
     1296    #setup 360 points on that ring for "good" fit
     1297    Ring = makeRing(ellipse,scalex,scaley,self.imScale,self.ImageZ)
     1298    if Ring:
     1299        ellipse = FitEllipse(Ring)
     1300    else:
     1301        print '1st ring not sufficiently complete to proceed'
     1302        return False
     1303    print 'inner ring:',ellipse
     1304    data['center'] = ellipse[0]           #not right!! (but useful for now)
     1305#    data['ellipses'].append(ellipse[:])
    12941306    self.PlotImage()
     1307   
     1308    #setup for calibration
     1309    data['rings'] = []
     1310    if not data['calibrant']:
     1311        print 'no calibration material selected'
     1312        return True
     1313    Bravais,cell = calFile.Calibrants[data['calibrant']]
     1314    A = cell2A(cell)
     1315    wave = data['wavelength']
     1316    dist = data['distance']
     1317    refine = data['refine']             #flags for center, wavelength,distance, tilt angle, tilt rotation
     1318    HKL = GenHBravais(0.5,Bravais,A)
     1319    dList = []
     1320    for i,H in enumerate(HKL):
     1321        cent,phi,semiradii = ellipse
     1322        meanRadius = semiradii[0]**2/semiradii[1]
     1323        dsp = 1.0/math.sqrt(calc_rDsq(H,A))
     1324        tth = 2.0*asind(wave/(2.*dsp))
     1325        radius = dist*tand(tth)
     1326        dList.append([dsp,radius])
     1327        semiradii[0] *= radius/meanRadius
     1328        semiradii[1] *= radius/meanRadius       
     1329        Ring = makeRing(ellipse,scalex,scaley,self.imScale,self.ImageZ)
     1330        if Ring:
     1331            data['rings'].append(Ring)
     1332            ellipse = FitEllipse(Ring)
     1333            cosb = (1.-ellipse[2][0]**2/ellipse[2][1]**2)*cosd(tth)**2
     1334            print 'for ring #',i,cosb,(ellipse[2][0]**2/ellipse[2][1])/tand(tth)
     1335            data['ellipses'].append(copy.deepcopy(ellipse))
     1336            self.PlotImage()
     1337        else:
     1338            break
     1339    self.PlotImage()
     1340   
     1341           
     1342   
    12951343       
    12961344   
  • trunk/GSASIIgrid.py

    r27 r29  
    952952    def OnClearCalib(event):
    953953        data['ring'] = []
     954        data['rings'] = []
    954955        data['ellipses'] = []
    955956        self.PlotImage()
  • trunk/ImageCalibrants.py

    r14 r29  
    11Calibrants={
     2'':(0,(0,0,0,0,0,0)),
    23'LaB6  SRM660a':(2,(4.1569162,4.1569162,4.1569162,90,90,90)),
    34'LaB6  SRM660': (2,(4.15695,4.15695,4.15695,90,90,90)),
    45'Si    SRM640c':(0,(5.4311946,5.4311946,5.4311946,90,90,90)),
    5 'CeO2  SRM674b':(2,(5.411651,5.411651,5.411651,90,90,90)),
     6'CeO2  SRM674b':(0,(5.411651,5.411651,5.411651,90,90,90)),
    67'Al2O3 SRM676a':(3,(4.759091,4.759091,12.991779,90,90,120))}
Note: See TracChangeset for help on using the changeset viewer.