Changeset 29

Timestamp:

Feb 11, 2010 5:46:51 PM (12 years ago)

Author:

vondreel

Message:

next stage of image calibration

Location:

trunk

Files:

• GSASII.py (modified) (5 diffs)
• GSASIIIO.py (modified) (1 diff)
• GSASIIcomp.py (modified) (2 diffs)
• GSASIIgrid.py (modified) (1 diff)
• ImageCalibrants.py (modified) (1 diff)

trunk/GSASII.py

@@ -398 +398 @@
         dlg = wx.FileDialog(self, 'Choose image file', '.', '', \
             'MAR345 (*.mar3450)|*.mar3450|ADSC Image (*.img)|*.img \
-            |Perkin-Elmer TIF (*.tif)|*.tif|GE Image sum (*.sum)|*.sum|All files (*.*)|*.*',wx.OPEN)
+            |Perkin-Elmer TIF (*.tif)|*.tif \
+            |GE Image sum (*.sum)|*.sum|GE Image avg (*.avg) \
+            |*.avg|All files (*.*)|*.*',wx.OPEN)
         if self.dirname:
             dlg.SetDirectory(self.dirname)
@@ -414 +416 @@
                 elif ext == '.mar3450':
                     Comments,Data,Size,Image = G2IO.GetMAR345Data(self.imagefile)
-                elif ext == '.sum':
+                elif ext in ['.sum','.avg']:
                     Comments,Data,Size,Image = G2IO.GetGEsumData(self.imagefile)
                 if Comments:
@@ -432 +434 @@
                         Data['showLines'] = False
                         Data['ring'] = []
+                        Data['rings'] = []
                         Data['ellipses'] = []
                         Data['masks'] = []
@@ -1303 +1306 @@
                     xpos = int(event.xdata)*self.imScale
                     ypos = int(event.ydata)*self.imScale
-                    self.pdplot.canvas.SetToolTipString('%6d'%(self.ImageZ[ypos][xpos]))
+                    if (0 <= xpos <= size) and (0 <= ypos <= size):
+                        self.pdplot.canvas.SetToolTipString('%6d'%(self.ImageZ[ypos][xpos]))
 
         def OnImPlotKeyPress(event):
@@ -1456 +1460 @@
             yring *= scaley
             ax.text(xring,yring,'+',ha='center',va='center',picker=3)
+#        for ring in Data['rings']:
+#            for xring,yring in ring:
+#                xring *= scalex
+#                yring *= scaley
+#                ax.text(xring,yring,'+',ha='center',va='center')            
         for ellipse in Data['ellipses']:
             cent,phi,[width,height] = ellipse

trunk/GSASIIIO.py

@@ -399 +399 @@
 def GetGEsumData(filename):
     import array as ar
-    print 'Read GE sum file: ',filename
+    print 'Read GE sum file: ',filename    
     File = open(filename,'rb')
     size = 2048
+    if '.sum' in filename:
+        head = ['GE detector sum data from APS 1-ID',]
+    if '.avg' in filename:
+#        image = np.zeros(shape=(size,size),dtype=np.int16)
+        head = ['GE detector avg data from APS 1-ID',]
     image = np.zeros(shape=(size,size),dtype=np.int32)
-    head = ['GE detector sum data',]
     row = 0
     pos = 0
     while row < size:
         File.seek(pos)
-        line = ar.array('f',File.read(4*size))
+        if '.sum' in filename:
+            line = ar.array('f',File.read(4*size))
+            pos += 4*size
+        elif '.avg' in filename:
+            line = ar.array('H',File.read(2*size))
+            pos += 2*size
         image[row] = np.asarray(line)
         row += 1
-        pos += 4*size
-    data = {'pixelSize':(200,200),'wavelength':0.10,'distance':250.0,'center':[204.8,204.8]}  
+    data = {'pixelSize':(200,200),'wavelength':0.15,'distance':250.0,'center':[204.8,204.8]}  
     return head,data,size,image
     File.close()    

trunk/GSASIIcomp.py

@@ -1215 +1215 @@
     if f3/a3 < 0 or f3/b3 < 0:
         return 0
-    sr1 = math.sqrt(f3/a3)
-    sr2 = math.sqrt(f3/b3)
-    return ell,cent,phi,[sr1,sr2]
+    if phi > 0.:        
+        sr1 = math.sqrt(f3/a3)
+        sr2 = math.sqrt(f3/b3)
+    else:
+        sr1 = math.sqrt(f3/a3)
+        sr2 = math.sqrt(f3/b3)
+        phi *= -1.        
+    return cent,phi,[sr1,sr2]
         
 def FitCircle(ring):
@@ -1239 +1244 @@
 def ImageLocalMax(image,w,Xpos,Ypos):
     w2 = w*2
-    Z = image[Ypos-w:Ypos+w,Xpos-w:Xpos+w]
-    Zmax = np.argmax(Z)
-    Zmin = np.argmin(Z)
-    Xpos += Zmax%w2-w
-    Ypos += Zmax/w2-w
-    return Xpos,Ypos,np.ravel(Z)[Zmax],np.ravel(Z)[Zmin]
-    
-def makeRing(cent,phi,semiradii,scalex,scaley,imScale,image):
+    size = len(image)
+    if (w < Xpos < size-w) and (w < Ypos < size-w):
+        Z = image[Ypos-w:Ypos+w,Xpos-w:Xpos+w]
+        Zmax = np.argmax(Z)
+        Zmin = np.argmin(Z)
+        Xpos += Zmax%w2-w
+        Ypos += Zmax/w2-w
+        return Xpos,Ypos,np.ravel(Z)[Zmax],np.ravel(Z)[Zmin]
+    else:
+        return 0,0,0,0
+    
+def makeRing(ellipse,scalex,scaley,imScale,image):
+    cent,phi,semiradii = ellipse
     cphi = cosd(phi)
     sphi = sind(phi)
     ring = []
-    for a in range(0,360,4):
+    for a in range(0,360):
         x = semiradii[0]*cosd(a)
         y = semiradii[1]*sind(a)
         X = (cphi*x-sphi*y+cent[0])*scalex*imScale
         Y = (sphi*x+cphi*y+cent[1])*scaley*imScale
-        if (0 < X < len(image)) and (0 < Y < len(image)):        
-            X,Y,I,J = ImageLocalMax(image,20,X,Y)
-            if I and J:
-                X /= scalex*imScale
-                Y /= scaley*imScale
-                ring.append([X,Y])    
+        X,Y,I,J = ImageLocalMax(image,20,X,Y)      
+        if I and J:
+            X /= scalex*imScale
+            Y /= scaley*imScale
+            ring.append([X,Y])
+    if len(ring) < 75:             #want more than 1/3 of a circle
+        return []
     return ring
         
 def ImageCalibrate(self,data):
+    import copy
     import ImageCalibrants as calFile
     print 'image calibrate'
-    if not data['calibrant']:
-        print 'no calibration material selected'
-        return
-    Bravais,cell = calFile.Calibrants[data['calibrant']]
     ring = data['ring']
-    data['ellipses'] = ellipses = []
     scalex = data['scalex']             # = 1000./(pixelSize[0]*self.imScale)
     scaley = data['scaley']
     if len(ring) < 5:
         print 'not enough inner ring points for ellipse'
-        return
+        return False
+        
+    #fit start points on inner ring
+    data['ellipses'] = []
     outB = FitEllipse(ring)
     if outB:
-        ell,cent,phi,semiradii = outB
-        print 'start ellipse coeff.',ell
-        print 'start:',cent,phi,semiradii
-        data['ellipses'].append([cent,phi,semiradii])
-        self.PlotImage()
+        ellipse = outB
+        print 'start:',ellipse
     else:
         return False
-    ring = makeRing(cent,phi,semiradii,scalex,scaley,self.imScale,self.ImageZ)
-    ell,cent,phi,semiradii = FitEllipse(ring)
-    print '1st ring ellipse coeff.',ell
-    print '1st ring:',cent,phi,semiradii
-    data['ellipses'].append([cent,phi,semiradii])
-    data['center'] = cent           #not right!! (but useful for now)
+        
+    #setup 360 points on that ring for "good" fit
+    Ring = makeRing(ellipse,scalex,scaley,self.imScale,self.ImageZ)
+    if Ring:
+        ellipse = FitEllipse(Ring)
+    else:
+        print '1st ring not sufficiently complete to proceed'
+        return False
+    print 'inner ring:',ellipse
+    data['center'] = ellipse[0]           #not right!! (but useful for now)
+#    data['ellipses'].append(ellipse[:])
     self.PlotImage()
+    
+    #setup for calibration
+    data['rings'] = []
+    if not data['calibrant']:
+        print 'no calibration material selected'
+        return True
+    Bravais,cell = calFile.Calibrants[data['calibrant']]
+    A = cell2A(cell)
+    wave = data['wavelength']
+    dist = data['distance']
+    refine = data['refine']             #flags for center, wavelength,distance, tilt angle, tilt rotation
+    HKL = GenHBravais(0.5,Bravais,A)
+    dList = []
+    for i,H in enumerate(HKL):
+        cent,phi,semiradii = ellipse
+        meanRadius = semiradii[0]**2/semiradii[1]
+        dsp = 1.0/math.sqrt(calc_rDsq(H,A))
+        tth = 2.0*asind(wave/(2.*dsp))
+        radius = dist*tand(tth)
+        dList.append([dsp,radius])
+        semiradii[0] *= radius/meanRadius
+        semiradii[1] *= radius/meanRadius        
+        Ring = makeRing(ellipse,scalex,scaley,self.imScale,self.ImageZ)
+        if Ring:
+            data['rings'].append(Ring)
+            ellipse = FitEllipse(Ring)
+            cosb = (1.-ellipse[2][0]**2/ellipse[2][1]**2)*cosd(tth)**2
+            print 'for ring #',i,cosb,(ellipse[2][0]**2/ellipse[2][1])/tand(tth)
+            data['ellipses'].append(copy.deepcopy(ellipse))
+            self.PlotImage()
+        else:
+            break
+    self.PlotImage()
+    
+            
+    
         
     

trunk/GSASIIgrid.py

@@ -952 +952 @@
     def OnClearCalib(event):
         data['ring'] = []
+        data['rings'] = []
         data['ellipses'] = []
         self.PlotImage()

trunk/ImageCalibrants.py

@@ -1 +1 @@
 Calibrants={
+'':(0,(0,0,0,0,0,0)),
 'LaB6  SRM660a':(2,(4.1569162,4.1569162,4.1569162,90,90,90)),
 'LaB6  SRM660': (2,(4.15695,4.15695,4.15695,90,90,90)),
 'Si    SRM640c':(0,(5.4311946,5.4311946,5.4311946,90,90,90)),
-'CeO2  SRM674b':(2,(5.411651,5.411651,5.411651,90,90,90)),
+'CeO2  SRM674b':(0,(5.411651,5.411651,5.411651,90,90,90)),
 'Al2O3 SRM676a':(3,(4.759091,4.759091,12.991779,90,90,120))}