Changeset 1149


Ignore:
Timestamp:
Nov 24, 2013 2:50:12 PM (8 years ago)
Author:
vondreele
Message:

cleanup of image stuff
remove wx from G2image, limit fitting to ellipses - break on hyperbola

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIimage.py

    r1148 r1149  
    1414Ellipse fitting & image integration
    1515
    16 needs some minor refactoring to remove wx code
    17 actually not easy in this case wx.ProgressDialog needs #of blocks to process when started
    18 not known in G2imageGUI.
    1916'''
    2017
    2118import math
    22 import wx
    2319import time
    2420import numpy as np
     
    510506        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
    511507    Found = False
     508    wave = data['wavelength']
    512509    for H in HKL:
    513510        dsp = H[3]
     511        tth = 2.0*asind(wave/(2.*dsp))
     512        if tth+abs(data['tilt']) > 90.:
     513            print 'next line is a hyperbola - search stopped'
     514            break
    514515        ellipse = GetEllipse(dsp,data)
    515516        Ring,delt,sumI = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
     
    585586    data['rings'] = []
    586587    data['ellipses'] = []
    587     if not data['calibrant']  or 'None' in data['calibrant']:
     588    if not data['calibrant']:
    588589        print 'no calibration material selected'
    589590        return True
     
    609610#1st estimate of tilt; assume ellipse - don't know sign though
    610611    tilt = npasind(np.sqrt(max(0.,1.-(radii[0]/radii[1])**2))*ctth)
     612    if not tilt:
     613        print 'WARNING - selected ring was fitted as a circle'
     614        print ' - if detector was tilted we suggest you skip this ring - WARNING'
    611615#1st estimate of dist: sample to detector normal to plane
    612616    data['distance'] = dist = radii[0]**2/(ttth*radii[1])
     
    623627    ellipsep = GetEllipse2(tth,0.,dist,centp,tilt,phi)
    624628    Ringp,delt,sumIp = makeRing(dsp,ellipsep,2,cutoff,scalex,scaley,self.ImageZ)
    625     outEp,errp = FitRing(Ringp,True)
    626     print '+',ellipsep,errp
     629    if len(Ringp) > 10:
     630        outEp,errp = FitRing(Ringp,True)
     631    else:
     632        errp = 1e6
     633#    print '+',ellipsep,errp
    627634    ellipsem = GetEllipse2(tth,0.,dist,centm,-tilt,phi)
    628635    Ringm,delt,sumIm = makeRing(dsp,ellipsem,2,cutoff,scalex,scaley,self.ImageZ)
    629     outEm,errm = FitRing(Ringm,True)
    630     print '-',ellipsem,errm
     636    if len(Ringm) > 10:
     637        outEm,errm = FitRing(Ringm,True)
     638    else:
     639        errm = 1e6
     640#    print '-',ellipsem,errm
    631641    if errp < errm:
    632642        data['tilt'] = tilt
     
    646656        dsp = H[3]
    647657        tth = 2.0*asind(wave/(2.*dsp))
     658        if tth+abs(data['tilt']) > 90.:
     659            print 'next line is a hyperbola - search stopped'
     660            break
    648661        print 'HKLD:',H[:4],'2-theta: %.4f'%(tth)
    649662        elcent,phi,radii = ellipse = GetEllipse(dsp,data)
     
    746759    return tax,tay,taz
    747760   
    748 def ImageIntegrate(image,data,masks):
     761def ImageIntegrate(image,data,masks,dlg=None):
    749762    'Needs a doc string'
    750763    import histogram2d as h2d
     
    766779    nYBlks = (Ny-1)/blkSize+1
    767780    Nup = nXBlks*nYBlks*3+3
    768     dlg = wx.ProgressDialog("Elapsed time","2D image integration",Nup,
    769         style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE)
    770     try:
    771         t0 = time.time()
    772         Nup = 0
     781    t0 = time.time()
     782    Nup = 0
     783    if dlg:
    773784        dlg.Update(Nup)
    774         for iBlk in range(nYBlks):
    775             iBeg = iBlk*blkSize
    776             iFin = min(iBeg+blkSize,Ny)
    777             for jBlk in range(nXBlks):
    778                 jBeg = jBlk*blkSize
    779                 jFin = min(jBeg+blkSize,Nx)               
    780 #                print 'Process map block:',iBlk,jBlk,' limits:',iBeg,iFin,jBeg,jFin
    781                 TA,tam = Make2ThetaAzimuthMap(data,masks,(iBeg,iFin),(jBeg,jFin))           #2-theta & azimuth arrays & create position mask
    782                
    783                 Nup += 1
     785    for iBlk in range(nYBlks):
     786        iBeg = iBlk*blkSize
     787        iFin = min(iBeg+blkSize,Ny)
     788        for jBlk in range(nXBlks):
     789            jBeg = jBlk*blkSize
     790            jFin = min(jBeg+blkSize,Nx)               
     791#            print 'Process map block:',iBlk,jBlk,' limits:',iBeg,iFin,jBeg,jFin
     792            TA,tam = Make2ThetaAzimuthMap(data,masks,(iBeg,iFin),(jBeg,jFin))           #2-theta & azimuth arrays & create position mask
     793           
     794            Nup += 1
     795            if dlg:
    784796                dlg.Update(Nup)
    785                 Block = image[iBeg:iFin,jBeg:jFin]
    786                 tax,tay,taz = Fill2ThetaAzimuthMap(masks,TA,tam,Block)    #and apply masks
    787                 del TA,tam
    788                 Nup += 1
     797            Block = image[iBeg:iFin,jBeg:jFin]
     798            tax,tay,taz = Fill2ThetaAzimuthMap(masks,TA,tam,Block)    #and apply masks
     799            del TA,tam
     800            Nup += 1
     801            if dlg:
    789802                dlg.Update(Nup)
    790                 tax = np.where(tax > LRazm[1],tax-360.,tax)                 #put azm inside limits if possible
    791                 tax = np.where(tax < LRazm[0],tax+360.,tax)
    792                 if any([tax.shape[0],tay.shape[0],taz.shape[0]]):
    793                     NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,numAzms,numChans,LRazm,LUtth,Dazm,Dtth,NST,H0)
    794 #                print 'block done'
    795                 del tax,tay,taz
    796                 Nup += 1
     803            tax = np.where(tax > LRazm[1],tax-360.,tax)                 #put azm inside limits if possible
     804            tax = np.where(tax < LRazm[0],tax+360.,tax)
     805            if any([tax.shape[0],tay.shape[0],taz.shape[0]]):
     806                NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,numAzms,numChans,LRazm,LUtth,Dazm,Dtth,NST,H0)
     807#            print 'block done'
     808            del tax,tay,taz
     809            Nup += 1
     810            if dlg:
    797811                dlg.Update(Nup)
    798         NST = np.array(NST)
    799         H0 = np.divide(H0,NST)
    800         H0 = np.nan_to_num(H0)
    801         del NST
    802         if Dtth:
    803             H2 = np.array([tth for tth in np.linspace(LUtth[0],LUtth[1],numChans+1)])
    804         else:
    805             H2 = np.array(LUtth)
    806         if Dazm:       
    807             H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)])
    808         else:
    809             H1 = LRazm
    810         H0 /= npcosd(H2[:-1])**2
    811         if data['Oblique'][1]:
    812             H0 /= G2pwd.Oblique(data['Oblique'][0],H2[:-1])
    813         if 'SASD' in data['type'] and data['PolaVal'][1]:
    814             H0 /= np.array([G2pwd.Polarization(data['PolaVal'][0],H2[:-1],Azm=azm)[0] for azm in H1[:-1]+np.diff(H1)/2.])
    815         Nup += 1
     812    NST = np.array(NST)
     813    H0 = np.divide(H0,NST)
     814    H0 = np.nan_to_num(H0)
     815    del NST
     816    if Dtth:
     817        H2 = np.array([tth for tth in np.linspace(LUtth[0],LUtth[1],numChans+1)])
     818    else:
     819        H2 = np.array(LUtth)
     820    if Dazm:       
     821        H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)])
     822    else:
     823        H1 = LRazm
     824    H0 /= npcosd(H2[:-1])**2
     825    if data['Oblique'][1]:
     826        H0 /= G2pwd.Oblique(data['Oblique'][0],H2[:-1])
     827    if 'SASD' in data['type'] and data['PolaVal'][1]:
     828        H0 /= np.array([G2pwd.Polarization(data['PolaVal'][0],H2[:-1],Azm=azm)[0] for azm in H1[:-1]+np.diff(H1)/2.])
     829    Nup += 1
     830    if dlg:
    816831        dlg.Update(Nup)
    817         t1 = time.time()
    818     finally:
    819         dlg.Destroy()
     832    t1 = time.time()
    820833    print 'Integration complete'
    821834    print "Elapsed time:","%8.3f"%(t1-t0), "s"
  • trunk/GSASIIimgGUI.py

    r1148 r1149  
    7373           
    7474    def OnIntegrate(event):
    75         if data['background image'][0]:
    76             maskCopy = copy.deepcopy(masks)
    77             backImg = data['background image'][0]
    78             backScale = data['background image'][1]
    79             id = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, backImg)
    80             Npix,imagefile = G2frame.PatternTree.GetItemPyData(id)
    81             backImage = G2IO.GetImageData(G2frame,imagefile,True)*backScale
    82             sumImage = G2frame.ImageZ+backImage
    83             sumMin = np.min(sumImage)
    84             sumMax = np.max(sumImage)
    85             maskCopy['Thresholds'] = [(sumMin,sumMax),[sumMin,sumMax]]
    86             G2frame.Integrate = G2img.ImageIntegrate(sumImage,data,maskCopy)
    87         else:
    88             G2frame.Integrate = G2img.ImageIntegrate(G2frame.ImageZ,data,masks)
    89 #        G2plt.PlotIntegration(G2frame,newPlot=True)
    90         G2IO.SaveIntegration(G2frame,G2frame.PickId,data)
     75        dlg = wx.ProgressDialog("Elapsed time","2D image integration",Nup,
     76            style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE)
     77        try:
     78            if data['background image'][0]:
     79                maskCopy = copy.deepcopy(masks)
     80                backImg = data['background image'][0]
     81                backScale = data['background image'][1]
     82                id = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, backImg)
     83                Npix,imagefile = G2frame.PatternTree.GetItemPyData(id)
     84                backImage = G2IO.GetImageData(G2frame,imagefile,True)*backScale
     85                sumImage = G2frame.ImageZ+backImage
     86                sumMin = np.min(sumImage)
     87                sumMax = np.max(sumImage)
     88                maskCopy['Thresholds'] = [(sumMin,sumMax),[sumMin,sumMax]]
     89                G2frame.Integrate = G2img.ImageIntegrate(sumImage,data,maskCopy,dlg)
     90            else:
     91                G2frame.Integrate = G2img.ImageIntegrate(G2frame.ImageZ,data,masks,dlg)
     92    #        G2plt.PlotIntegration(G2frame,newPlot=True)
     93            G2IO.SaveIntegration(G2frame,G2frame.PickId,data)
     94        finally:
     95            dlg.Destroy()
    9196        for item in G2frame.MakePDF: item.Enable(True)
    9297       
     
    116121                        ifintegrate,name,id = item
    117122                        if ifintegrate:
    118                             id = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, name)
    119                             Npix,imagefile = G2frame.PatternTree.GetItemPyData(id)
    120                             image = G2IO.GetImageData(G2frame,imagefile,True)
    121                             Id = G2gd.GetPatternTreeItemId(G2frame,id, 'Image Controls')
    122                             Data = G2frame.PatternTree.GetItemPyData(Id)
    123                             backImage = []
    124                             if Data['background image'][0]:
    125                                 backImg = Data['background image'][0]
    126                                 backScale = Data['background image'][1]
    127                                 id = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, backImg)
     123                            dlgp = wx.ProgressDialog("Elapsed time","2D image integration",Nup,
     124                                style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE)
     125                            try:
     126                                id = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, name)
    128127                                Npix,imagefile = G2frame.PatternTree.GetItemPyData(id)
    129                                 backImage = G2IO.GetImageData(G2frame,imagefile,True)*backScale
    130                             try:
    131                                 Masks = G2frame.PatternTree.GetItemPyData(
    132                                     G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Masks'))
    133                             except TypeError:       #missing Masks
    134                                 Imin,Imax = Data['Range']
    135                                 Masks = {'Points':[],'Rings':[],'Arcs':[],'Polygons':[],'Frames':[],'Thresholds':[(Imin,Imax),[Imin,Imax]]}
    136                                 G2frame.PatternTree.SetItemPyData(
    137                                     G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Masks'),Masks)
    138                             if len(backImage):                               
    139                                 G2frame.Integrate = G2img.ImageIntegrate(image+backImage,Data,Masks)
    140                             else:
    141                                 G2frame.Integrate = G2img.ImageIntegrate(image,Data,Masks)
    142 #                            G2plt.PlotIntegration(G2frame,newPlot=True,event=event)
    143                             G2IO.SaveIntegration(G2frame,Id,Data)
     128                                image = G2IO.GetImageData(G2frame,imagefile,True)
     129                                Id = G2gd.GetPatternTreeItemId(G2frame,id, 'Image Controls')
     130                                Data = G2frame.PatternTree.GetItemPyData(Id)
     131                                backImage = []
     132                                if Data['background image'][0]:
     133                                    backImg = Data['background image'][0]
     134                                    backScale = Data['background image'][1]
     135                                    id = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, backImg)
     136                                    Npix,imagefile = G2frame.PatternTree.GetItemPyData(id)
     137                                    backImage = G2IO.GetImageData(G2frame,imagefile,True)*backScale
     138                                try:
     139                                    Masks = G2frame.PatternTree.GetItemPyData(
     140                                        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Masks'))
     141                                except TypeError:       #missing Masks
     142                                    Imin,Imax = Data['Range']
     143                                    Masks = {'Points':[],'Rings':[],'Arcs':[],'Polygons':[],'Frames':[],'Thresholds':[(Imin,Imax),[Imin,Imax]]}
     144                                    G2frame.PatternTree.SetItemPyData(
     145                                        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Masks'),Masks)
     146                                if len(backImage):                               
     147                                    G2frame.Integrate = G2img.ImageIntegrate(image+backImage,Data,Masks,dlgp)
     148                                else:
     149                                    G2frame.Integrate = G2img.ImageIntegrate(image,Data,Masks,dlgp)
     150#                               G2plt.PlotIntegration(G2frame,newPlot=True,event=event)
     151                                G2IO.SaveIntegration(G2frame,Id,Data)
     152                            finally:
     153                                dlgp.Destroy()
    144154            finally:
    145155                dlg.Destroy()
Note: See TracChangeset for help on using the changeset viewer.