Changeset 1184


Ignore:
Timestamp:
Jan 7, 2014 3:45:30 PM (9 years ago)
Author:
vondreele
Message:

First working version of strain refinement from image data

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIgrid.py

    r1183 r1184  
    37553755        G2imG.UpdateStressStrain(G2frame,data)
    37563756        G2plt.PlotImage(G2frame)
     3757        G2plt.PlotStrain(G2frame,data,newPlot=True)
    37573758    elif G2frame.PatternTree.GetItemText(item) == 'HKL Plot Controls':
    37583759        G2frame.PickId = item

  • trunk/GSASIIimage.py

    r1180 r1184  
    236236        amb = radii[1]-radii[0]
    237237        return np.pi*apb*(1+3*(amb/apb)**2/(10+np.sqrt(4-3*(amb/apb)**2)))
     238       
    238239    cent,phi,radii = ellipse
    239     cphi = cosd(phi)
    240     sphi = sind(phi)
     240    cphi = cosd(phi-90.)        #convert to major axis rotation
     241    sphi = sind(phi-90.)
    241242    ring = []
    242     sumI = 0
    243     amin = 0
    244     amax = 360
    245     C = int(ellipseC())
    246     for i in range(0,C,1):
     243    C = int(ellipseC())         #ring circumference
     244    for i in range(0,C,1):      #step around ring in 1mm increments
    247245        a = 360.*i/C
    248         x = radii[0]*cosd(a)
    249         y = radii[1]*sind(a)
     246        x = radii[1]*cosd(a)        #major axis
     247        y = radii[0]*sind(a)
    250248        X = (cphi*x-sphi*y+cent[0])*scalex      #convert mm to pixels
    251249        Y = (sphi*x+cphi*y+cent[1])*scaley
    252250        X,Y,I,J = ImageLocalMax(image,pix,X,Y)
    253251        if I and J and I/J > reject:
    254             sumI += I/J
    255252            X += .5                             #set to center of pixel
    256253            Y += .5
    257             X /= scalex                         #convert to mm
     254            X /= scalex                         #convert back to mm
    258255            Y /= scaley
    259             amin = min(amin,a)
    260             amax = max(amax,a)
    261             if [X,Y,dsp] not in ring:
     256            if [X,Y,dsp] not in ring:           #no duplicates!
    262257                ring.append([X,Y,dsp])
    263     if len(ring) < 10:             #want more than 10 deg
     258    if len(ring) < 10:
    264259        ring = []
    265260    return ring
    … …  
    393388    azm = (npatan2d(dy,dx)+azmthoff+720.)%360.
    394389    G = D/data['distance']**2       #for geometric correction = 1/cos(2theta)^2 if tilt=0.
    395     return tth,azm,G,dsp
     390    return np.array([tth,azm,G,dsp])
    396391   
    397392def GetTth(x,y,data):
    … …  
    404399   
    405400def GetTthAzmD(x,y,data):
    406     '''Give 2-theta, azimuth & geometric correction values for detector x,y position;
     401    '''Give 2-theta, azimuth & d-spacing values for detector x,y position;
    407402     calibration info in data
    408403    '''
    … …  
    858853    return H0,H1,H2
    859854   
    860 def FitStrSta(Image,StrSta,Controls,Masks):
     855def MakeStrStaRing(ring,Image,Controls):
     856    ellipse = GetEllipse(ring['Dset'],Controls)
     857    pixSize = Controls['pixelSize']
     858    scalex = 1000./pixSize[0]
     859    scaley = 1000./pixSize[1]
     860    Ring = np.array(makeRing(ring['Dset'],ellipse,ring['pixLimit'],ring['cutoff'],scalex,scaley,Image)).T   #returns x,y,dsp for each point in ring
     861    ring['ImxyObs'] = Ring[:2]
     862    TA = GetTthAzm(Ring[0],Ring[1],Controls)       #convert x,y to tth,azm
     863    TA[0] = Controls['wavelength']/(2.*npsind(TA[0]/2.))      #convert 2th to d
     864    ring['ImtaObs'] = TA
     865    ring['ImtaCalc'] = np.zeros_like(ring['ImtaObs'])
     866    Ring[0] = TA[0]
     867    Ring[1] = TA[1]
     868    return Ring,ring
     869   
     870def FitStrSta(Image,StrSta,Controls):
    861871    'Needs a doc string'
    862872   
    … …  
    865875    wave = Controls['wavelength']
    866876    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
    867     pixSize = StaControls['pixelSize']
    868     scalex = 1000./pixSize[0]
    869     scaley = 1000./pixSize[1]
    870877    rings = []
    871878
    872879    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
    873         ellipse = GetEllipse(ring['Dset'],StaControls)
    874         Ring = np.array(makeRing(ring['Dset'],ellipse,ring['pixLimit'],ring['cutoff'],scalex,scaley,Image)).T
    875         ring['ImxyObs'] = np.array(Ring[:2])      #need to apply masks to this to eliminate bad points
    876         ring['ImtaObs'] = GetTthAzm(Ring[0],Ring[1],StaControls)       #convert x,y to tth,azm
    877         ring['ImtaCalc'] = np.zeros_like(ring['ImtaObs'])
    878         Ring[0] /= 2.                                           #convert to theta
     880        Ring,R = MakeStrStaRing(ring,Image,StaControls)
     881        ring.update(R)
    879882        if len(rings):
    880883            rings = np.concatenate((rings,Ring),axis=1)
    881884        else:
    882             rings = np.array(Ring)
     885            rings = np.array(Ring)     
    883886    E = StrSta['strain']
    884     p0 = [E[0][0],E[1][1],E[2][2],E[0][1],E[0][2],E[1][2]]
    885     E,dth = FitStrain(rings,p0,wave,phi)
     887    p0 = [E[0][0],E[0][1],E[1][1]]
     888    E = FitStrain(rings,p0,wave,phi)
     889    StrSta['strain'] = E
     890    CalcStrSta(StrSta,Controls)
     891   
     892def CalcStrSta(StrSta,Controls):
     893
     894    wave = Controls['wavelength']
     895    phi = StrSta['Sample phi']
     896    E = StrSta['strain']
    886897    for ring in StrSta['d-zero']:
    887898        th,azm = ring['ImtaObs']
    888         th0 = npasind(wave/(2.*ring['Dset']))
    889         V = np.sum(np.sum(E*calcFij(90.,phi,azm,th0).T,axis=2),axis=1)
    890         dth = 180.*V/(np.pi*10e6) #in degrees & microstrain units
    891         ring['ImtaCalc'] = np.array([(th0-dth)/2,azm])
    892     StrSta['strain'] = E
     899        th0 = np.ones_like(azm)*npasind(wave/(2.*ring['Dset']))
     900        V = 1.+np.sum(np.sum(E*calcFij(90.,phi,azm,th0).T/1.e6,axis=2),axis=1)
     901        ring['ImtaCalc'] = np.array([V*ring['Dset'],azm])
    893902
    894903def calcFij(omg,phi,azm,th):
    … …  
    933942       
    934943    def strainCalc(p,xyd,wave,phi):
    935 #        E = np.array([[p[0],p[3],p[4]],[p[3],p[1],p[5]],[p[4],p[5],p[2]]])
    936         E = np.array([[p[0],0,0],[0,p[1],0],[0,0,0]])
    937         th,azm,dsp = xyd
    938         th0 = npasind(wave/(2.*dsp))
    939         V = np.sum(np.sum(E*calcFij(90.,phi,azm,th).T,axis=2),axis=1)
    940         dth = 180.*V/(np.pi*10e6) #in degrees & microstrain units
    941         th0 += dth
    942         return th-th0
    943        
    944     names = ['e11','e22','e33','e12','e13','e23']
     944        E = np.array([[p[0],p[1],0],[p[1],p[2],0],[0,0,0]])
     945        dspo,azm,dsp = xyd
     946        th = npasind(wave/(2.0*dspo))
     947        V = 1.+np.sum(np.sum(E*calcFij(90.,phi,azm,th).T/1.e6,axis=2),axis=1)
     948        dspc = dsp*V
     949        return dspo-dspc
     950       
     951    names = ['e11','e12','e22','e33','e13','e23']
    945952    fmt = ['%12.2f','%12.2f','%12.2f','%12.2f','%12.2f','%12.5f']
    946     p1 = [1,-1]  
     953    p1 = [-20000.,0,5000.] 
    947954    result = leastsq(strainCalc,p1,args=(rings,wave,phi),full_output=True)
    948955    vals = list(result[0])
    … …  
    951958    ValSig = zip(names,fmt,vals,sig)
    952959    StrainPrint(ValSig)
    953 #    return np.array([[vals[0],vals[3],vals[4]],[vals[3],vals[1],vals[5]],[vals[4],vals[5],vals[2]]])
    954     return np.array([[vals[0],0,0],[0,vals[1],0],[0,0,0]]),result[2]['fvec']
    955    
    956        
     960    return np.array([[vals[0],vals[1],0],[vals[1],vals[2],0],[0,0,0]])
     961   
     962       

  • trunk/GSASIIimgGUI.py

    r1180 r1184  
    13211321   
    13221322    def OnAppendDzero(event):
    1323         data['d-zero'].append({'Dset':1.0,'Dcalc':0.0,'pixLimit':10,'cutoff':10.0,'ImxyObs':[[],[]],'ImxyCalc':[[],[]]})
     1323        data['d-zero'].append({'Dset':1.0,'Dcalc':0.0,'pixLimit':10,'cutoff':10.0,
     1324            'ImxyObs':[[],[]],'ImtaObs':[[],[]],'ImtaCalc':[[],[]]})
    13241325        UpdateStressStrain(G2frame,data)
    13251326           
    … …  
    13831384                    for key in keys2:
    13841385                        File.write("'"+key+"':"+':'+str(data2[key])+',')
    1385                     File.write("'ImxyObs':[[],[]],'Imxycalc':[[],[]]},\n")
     1386                    File.write("'ImxyObs':[[],[]],'ImtaObs':[[],[]],'Imtacalc':[[],[]]},\n")
    13861387                File.write('\t]\n}')
    13871388                File.close()
    … …  
    13941395        Controls = G2frame.PatternTree.GetItemPyData(
    13951396            G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Image Controls'))
    1396         G2img.FitStrSta(G2frame.ImageZ,data,Controls,Masks)
     1397        G2img.FitStrSta(G2frame.ImageZ,data,Controls)
    13971398        UpdateStressStrain(G2frame,data)
    13981399        G2plt.PlotExposedImage(G2frame,event=event)
    1399         G2plt.PlotStrain(G2frame,data,newPlot=False,type='Strain')
     1400        G2plt.PlotStrain(G2frame,data,newPlot=False)
    14001401       
    14011402    def SamSizer():
    … …  
    14521453            Obj.SetValue("%.5f"%(value))
    14531454            data['d-zero'][Indx[Obj.GetId()]]['Dset'] = value
     1455            Ring,R = G2img.MakeStrStaRing(data['d-zero'][Indx[Obj.GetId()]],G2frame.ImageZ,Controls)
     1456            data['d-zero'][Indx[Obj.GetId()]].update(R)
     1457            UpdateStressStrain(G2frame,data)
     1458            G2plt.PlotExposedImage(G2frame,event=event)
     1459            G2plt.PlotStrain(G2frame,data,newPlot=False)
    14541460           
    14551461        def OnDeleteDzero(event):
    … …  
    14571463            del(data['d-zero'][delIndx.index(Obj)])
    14581464            UpdateStressStrain(G2frame,data)
     1465            G2plt.PlotExposedImage(G2frame,event=event)
     1466            G2plt.PlotStrain(G2frame,data,newPlot=False)
    14591467       
    14601468        def OnCutOff(event):
    … …  
    14661474            Obj.SetValue("%.1f"%(value))
    14671475            data['d-zero'][Indx[Obj.GetId()]]['cutoff'] = value
     1476            G2plt.PlotExposedImage(G2frame,event=event)
     1477            G2plt.PlotStrain(G2frame,data,newPlot=False)
    14681478       
    14691479        def OnPixLimit(event):
    14701480            Obj = event.GetEventObject()
    14711481            data['d-zero'][Indx[Obj.GetId()]]['pixLimit'] = int(Obj.GetValue())
     1482            G2plt.PlotExposedImage(G2frame,event=event)
     1483            G2plt.PlotStrain(G2frame,data,newPlot=False)
    14721484           
    14731485        Indx = {}
    … …  
    15111523    def StrainSizer():
    15121524       
     1525        def OnTensor(event):
     1526            Obj = event.GetEventObject()
     1527            try:
     1528                value = float(Obj.GetValue())
     1529            except ValueError:
     1530                value = 1.0
     1531            i,j = Indx[Obj.GetId()]
     1532            data['strain'][i][j] = value
     1533            if i != j:
     1534                data['strain'][j][i] = value
     1535            UpdateStressStrain(G2frame,data)
     1536            G2img.CalcStrSta(data,Controls)
     1537            G2plt.PlotStrain(G2frame,data,newPlot=False)
     1538
     1539        Indx = {}
    15131540        strainSizer = wx.BoxSizer(wx.VERTICAL)
    15141541        strainSizer.Add(wx.StaticText(G2frame.dataDisplay,-1,label=(' Strain tensor:')),
    … …  
    15191546            for j in range(3):
    15201547                tensorSizer.Add(wx.StaticText(G2frame.dataDisplay,-1,label=names[i][j]),0,wx.ALIGN_CENTER_VERTICAL)
    1521                 tensorElem = wx.TextCtrl(G2frame.dataDisplay,-1,value='%.2f'%(data['strain'][i][j]),style=wx.TE_READONLY)
    1522                 tensorElem.SetBackgroundColour(VERY_LIGHT_GREY)
     1548#                tensorElem = wx.TextCtrl(G2frame.dataDisplay,-1,value='%.2f'%(data['strain'][i][j]),style=wx.TE_READONLY)
     1549#                tensorElem.SetBackgroundColour(VERY_LIGHT_GREY)
     1550                tensorElem = wx.TextCtrl(G2frame.dataDisplay,-1,value='%.2f'%(data['strain'][i][j]),style=wx.TE_PROCESS_ENTER)
     1551                tensorElem.Bind(wx.EVT_TEXT_ENTER,OnTensor)
     1552                tensorElem.Bind(wx.EVT_KILL_FOCUS,OnTensor)
     1553                Indx[tensorElem.GetId()] = [i,j]
    15231554                tensorSizer.Add(tensorElem,0,wx.ALIGN_CENTER_VERTICAL)
    15241555        strainSizer.Add(tensorSizer)
    … …  
    15271558    if G2frame.dataDisplay:
    15281559        G2frame.dataDisplay.Destroy()
     1560    Controls = G2frame.PatternTree.GetItemPyData(
     1561        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Image Controls'))       
    15291562    G2gd.SetDataMenuBar(G2frame,G2frame.dataFrame.StrStaMenu)
    15301563    G2frame.dataFrame.Bind(wx.EVT_MENU, OnAppendDzero, id=G2gd.wxID_APPENDDZERO)

  • trunk/GSASIIplot.py

    r1180 r1184  
    12551255################################################################################
    12561256           
    1257 def PlotStrain(G2frame,data,newPlot=False,type=''):
     1257def PlotStrain(G2frame,data,newPlot=False):
    12581258    '''plot of strain data, used for diagnostic purposes
    12591259    '''
    … …  
    12641264            Page.canvas.SetCursor(wx.CROSS_CURSOR)
    12651265            try:
    1266                 G2frame.G2plotNB.status.SetStatusText('X =%9.3f %s =%9.3f'%(xpos,type,ypos),1)                   
     1266                G2frame.G2plotNB.status.SetStatusText('d-spacing =%9.5f Azimuth =%9.3f'%(xpos,ypos),1)                   
    12671267            except TypeError:
    1268                 G2frame.G2plotNB.status.SetStatusText('Select '+type+' pattern first',1)
     1268                G2frame.G2plotNB.status.SetStatusText('Select Strain pattern first',1)
    12691269
    12701270    try:
    1271         plotNum = G2frame.G2plotNB.plotList.index(type)
     1271        plotNum = G2frame.G2plotNB.plotList.index('Strain')
    12721272        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
    12731273        if not newPlot:
    … …  
    12781278    except ValueError:
    12791279        newPlot = True
    1280         Plot = G2frame.G2plotNB.addMpl(type).gca()
    1281         plotNum = G2frame.G2plotNB.plotList.index(type)
     1280        Plot = G2frame.G2plotNB.addMpl('Strain').gca()
     1281        plotNum = G2frame.G2plotNB.plotList.index('Strain')
    12821282        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
    12831283        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
    … …  
    12861286    Page.SetFocus()
    12871287    G2frame.G2plotNB.status.DestroyChildren()
    1288     Plot.set_title(type)
    1289     Plot.set_xlabel(r'$\mathsf{2\theta}$',fontsize=14)
     1288    Plot.set_title('Strain')
     1289    Plot.set_xlabel(r'd-spacing',fontsize=14)
    12901290    Plot.set_ylabel(r'Azimuth',fontsize=14)
    12911291    colors=['b','g','r','c','m','k']
    … …  
    25982598            for N,ring in enumerate(StrSta['d-zero']):
    25992599                xring,yring = ring['ImxyObs']
    2600                 Plot.plot(xring,yring,colors[N%6]+'+')
    2601                 for xring,yring in np.array(ring['ImxyCalc']).T:
    2602                     Plot.add_artist(Polygon(ring['ImxyCalc'].T,ec='b',fc='none'))
    2603                     Plot.plot(xring,yring)
     2600                Plot.plot(xring,yring,colors[N%6]+'o')
    26042601        #masks - mask lines numbered after integration limit lines
    26052602        spots = Masks['Points']
Note: See TracChangeset for help on using the changeset viewer.