Changeset 2411 for trunk/GSASIIimage.py


Ignore:
Timestamp:
Aug 8, 2016 5:09:05 PM (5 years ago)
Author:
vondreele
Message:

fix image plot rescaling problem
add ring graininess determination - can be done sequentially
some plot fixes

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIimage.py

    r2348 r2411  
    2121import numpy.linalg as nl
    2222import numpy.ma as ma
     23import numpy.fft as fft
     24import scipy.signal as signal
    2325import polymask as pm
    2426from scipy.optimize import leastsq
     
    216218    xpix = int(Xpix)            #get reference corner of pixel chosen
    217219    ypix = int(Ypix)
    218     if (w < xpix < sizex-w) and (w < ypix < sizey-w) and image[ypix,xpix]:
    219         Z = image[ypix-w:ypix+w,xpix-w:xpix+w]
    220         Zmax = np.argmax(Z)
    221         Zmin = np.argmin(Z)
     220    if not w:
     221        ZMax = np.sum(image[ypix-1:ypix+1,xpix-1:xpix+1])
     222        return xpix,ypix,ZMax,0.0001
     223    if (w2 < xpix < sizex-w2) and (w2 < ypix < sizey-w2) and image[ypix,xpix]:
     224        ZMax = image[ypix-w:ypix+w,xpix-w:xpix+w]
     225        Zmax = np.argmax(ZMax)
     226        ZMin = image[ypix-w2:ypix+w2,xpix-w2:xpix+w2]
     227        Zmin = np.argmin(ZMin)
    222228        xpix += Zmax%w2-w
    223229        ypix += Zmax/w2-w
    224         return xpix,ypix,np.ravel(Z)[Zmax],max(0.0001,np.ravel(Z)[Zmin])   #avoid neg/zero minimum
     230        return xpix,ypix,np.ravel(ZMax)[Zmax],max(0.0001,np.ravel(ZMin)[Zmin])   #avoid neg/zero minimum
    225231    else:
    226232        return 0,0,0,0     
     
    243249    ring = []
    244250    C = int(ellipseC())         #ring circumference
     251    azm = []
    245252    for i in range(0,C,1):      #step around ring in 1mm increments
    246253        a = 360.*i/C
     
    257264            if [X,Y,dsp] not in ring:           #no duplicates!
    258265                ring.append([X,Y,dsp])
     266                azm.append(a)
    259267    if len(ring) < 10:
    260268        ring = []
    261     return ring
     269        azm = []
     270    return ring,azm
    262271   
    263272def GetEllipse2(tth,dxy,dist,cent,tilt,phi):
     
    554563            break
    555564        ellipse = GetEllipse(dsp,data)
    556         Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,ma.array(G2frame.ImageZ,mask=tam))
     565        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,ma.array(G2frame.ImageZ,mask=tam))[0]
    557566        if Ring:
    558567            if iH >= skip:
     
    623632    #setup 360 points on that ring for "good" fit
    624633    data['ellipses'].append(ellipse[:]+('g',))
    625     Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)
     634    Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0]
    626635    if Ring:
    627636        ellipse = FitEllipse(Ring)
    628         Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)    #do again
     637        Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0]    #do again
    629638        ellipse = FitEllipse(Ring)
    630639    else:
     
    670679        tth = npatan2d(radii[0],dist)
    671680        data['wavelength'] = wave =  2.0*dsp*sind(tth/2.0)
    672     Ring0 = makeRing(dsp,ellipse,3,cutoff,scalex,scaley,G2frame.ImageZ)
     681    Ring0 = makeRing(dsp,ellipse,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
    673682    ttth = nptand(tth)
    674683    stth = npsind(tth)
     
    706715            ellipsep = GetEllipse2(tth,0.,dist,centp,tilt,phi)
    707716            print fmt%('plus ellipse :',ellipsep[0][0],ellipsep[0][1],ellipsep[1],ellipsep[2][0],ellipsep[2][1])
    708             Ringp = makeRing(dsp,ellipsep,3,cutoff,scalex,scaley,G2frame.ImageZ)
     717            Ringp = makeRing(dsp,ellipsep,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
    709718            parmDict = {'dist':dist,'det-X':centp[0],'det-Y':centp[1],
    710719                'tilt':tilt,'phi':phi,'wave':wave,'dep':0.0}       
     
    720729            ellipsem = GetEllipse2(tth,0.,dist,centm,-tilt,phi)
    721730            print fmt%('minus ellipse:',ellipsem[0][0],ellipsem[0][1],ellipsem[1],ellipsem[2][0],ellipsem[2][1])
    722             Ringm = makeRing(dsp,ellipsem,3,cutoff,scalex,scaley,G2frame.ImageZ)
     731            Ringm = makeRing(dsp,ellipsem,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
    723732            if len(Ringm) > 10:
    724733                parmDict['tilt'] *= -1
     
    760769        data['ellipses'].append(copy.deepcopy(ellipse+('g',)))
    761770        if debug:   print fmt%('predicted ellipse:',elcent[0],elcent[1],phi,radii[0],radii[1])
    762         Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)
     771        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0]
    763772        if Ring:
    764773            data['rings'].append(np.array(Ring))
     
    979988    scalex = 1000./pixSize[0]
    980989    scaley = 1000./pixSize[1]
    981     Ring = np.array(makeRing(ring['Dset'],ellipse,ring['pixLimit'],ring['cutoff'],scalex,scaley,Image)).T   #returns x,y,dsp for each point in ring
     990    Ring = np.array(makeRing(ring['Dset'],ellipse,ring['pixLimit'],ring['cutoff'],scalex,scaley,Image)[0]).T   #returns x,y,dsp for each point in ring
    982991    if len(Ring):
    983992        ring['ImxyObs'] = copy.copy(Ring[:2])
     
    10011010    phi = StrSta['Sample phi']
    10021011    wave = Controls['wavelength']
     1012    pixelSize = Controls['pixelSize']
     1013    scalex = 1000./pixelSize[0]
     1014    scaley = 1000./pixelSize[1]
    10031015    StaType = StrSta['Type']
    10041016    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
     
    10131025            ring['Emat'] = val
    10141026            ring['Esig'] = esd
     1027            ellipse = FitEllipse(R['ImxyObs'].T)
     1028            ringxy,ringazm = makeRing(ring['Dcalc'],ellipse,0,0.,scalex,scaley,Image)
     1029            ringint = np.array([float(Image[int(y*scaley),int(x*scalex)]) for x,y in np.array(ringxy)[:,:2]])
     1030            ringint /= np.mean(ringint)
     1031            ring['Ivar'] = np.var(ringint)
     1032            print 'Variance in normalized ring intensity: %.3f'%(ring['Ivar'])
    10151033    CalcStrSta(StrSta,Controls)
     1034   
     1035def IntStrSta(Image,StrSta,Controls):
     1036    StaControls = copy.deepcopy(Controls)
     1037    pixelSize = Controls['pixelSize']
     1038    scalex = 1000./pixelSize[0]
     1039    scaley = 1000./pixelSize[1]
     1040    phi = StrSta['Sample phi']
     1041    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
     1042    RingsAI = []
     1043    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
     1044        Ring,R = MakeStrStaRing(ring,Image,StaControls)
     1045        if len(Ring):
     1046            ellipse = FitEllipse(R['ImxyObs'].T)
     1047            ringxy,ringazm = makeRing(ring['Dcalc'],ellipse,0,0.,scalex,scaley,Image)
     1048            ringint = np.array([float(Image[int(y*scaley),int(x*scalex)]) for x,y in np.array(ringxy)[:,:2]])
     1049            ringint /= np.mean(ringint)
     1050            print 'variance:',ring['Dcalc'],np.var(ringint)
     1051            RingsAI.append(np.array(zip(ringazm,ringint)).T)
     1052#            GSASIIpath.IPyBreak()
     1053    return RingsAI
    10161054   
    10171055def CalcStrSta(StrSta,Controls):
Note: See TracChangeset for help on using the changeset viewer.