Changeset 2760


Ignore:
Timestamp:
Mar 24, 2017 12:20:38 PM (5 years ago)
Author:
vondreele
Message:

new auto spot mask finder - works for diamond spots (not so well fro gritty patterns)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIimage.py

    r2731 r2760  
    2323import polymask as pm
    2424from scipy.optimize import leastsq
     25import scipy.signal as scsg
     26import scipy.cluster.vq as scvq
    2527import copy
    2628import GSASIIpath
     
    11161118   
    11171119def AutoSpotMasks(Image,Masks,Controls):
    1118     rollImage = lambda rho,roll: np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1)
    11191120    print 'auto spot search'
     1121    maxPks = 500
     1122    dmax = 1.5         #just beyond diamond 111 @ 2.0596
     1123    wave = Controls['wavelength']
     1124    tthMax = 2.0*asind(wave/(2.*dmax))
    11201125    pixelSize = Controls['pixelSize']
    1121     spotMask = ma.array(Image,mask=(Image<np.mean(Image)))
    1122     indices = (-1,0,1)
    1123     rolls = np.array([[ix,iy] for ix in indices for iy in indices])
    1124     for roll in rolls:
    1125         if np.any(roll):        #avoid [0,0]
    1126             spotMask = ma.array(spotMask,mask=(spotMask-rollImage(Image,roll)<0.))
    1127     mags = spotMask[spotMask.nonzero()]
    1128     indx = np.transpose(spotMask.nonzero())
    1129     nx,ny = Image.shape
    1130     jndx = []
    1131     for [ind,mag] in zip(indx,mags):
    1132         if (1 < ind[0] < nx-2) and (1 < ind[1] < ny-2):
    1133             cent = np.zeros((3,3))
    1134             cent[1,1] = mag
    1135             msk = np.array(Image[ind[0]-2:ind[0]+3,ind[1]-2:ind[1]+3])
    1136             msk = ma.array(msk,mask=(msk==mag))
    1137             if mag > 1.5*ma.mean(msk):
    1138                 jndx.append([ind[1]+.5,ind[0]+.5])
    1139     print 'Spots found: ',len(jndx)
    1140     jndx = np.array(jndx)
    1141     peaks = jndx*pixelSize/1000.
    1142     tth = GetTth(peaks.T[0],peaks.T[1],Controls)
    1143     histth,bins = np.histogram(tth,2500)
    1144     for shft in [-1,1]:
    1145         histmsk = ma.array(histth,mask=-(histth-np.roll(histth,shft)<2.))
    1146     histmsk = ma.array(histmsk,mask=(histmsk>5))
    1147     binmsk = np.array(ma.nonzero(histmsk))
    1148     digts = np.digitize(tth,bins)-1
    1149     PeaksList = []
    1150     for digt,peak in zip(digts,peaks):
    1151         if digt in binmsk:
    1152             PeaksList.append(peak)
    1153     #should be able to filter out spotty Bragg rings here
    1154     PeaksList = np.array(PeaksList)
    1155     Points = np.ones((PeaksList.shape[0],3))
    1156     Points[:,:2] = PeaksList
     1126    lap = np.array([[0,1,0],[1,-4,1],[0,1,0]],np.float32)
     1127    ck = scsg.cspline2d(Image.T,24.0)
     1128    spotMask1=-scsg.convolve2d(ck,lap,mode='same',boundary='symm')
     1129    for mul in range(10):
     1130        spotMask = ma.array(Image,mask=(spotMask1<(mul+1)*1.e4))
     1131        indx = np.transpose(spotMask.nonzero())
     1132        if not len(indx):       #smooth image - no sharp points
     1133            break
     1134        pts,distort = scvq.kmeans(np.array(indx,np.float32),maxPks,thresh=1.0)
     1135        peaks = pts*pixelSize/1000.
     1136        tths = GetTth(peaks.T[0],peaks.T[1],Controls)
     1137        masktths = ma.getmask(ma.array(tths,mask=tths>tthMax))
     1138        pts = pts[masktths,:]
     1139        A,B = np.mgrid[0:len(pts),0:len(pts)]   
     1140        M = np.sqrt(np.sum((pts[A]-pts[B])**2,axis=-1))
     1141        MX = ma.array(M,mask=(M<10.))
     1142        Indx = ma.nonzero(ma.getmask(MX))
     1143        print 'Spots found: ',mul,len(pts),distort,len(Indx[0])
     1144        if distort < .3:
     1145            break
     1146    if not len(Indx[0]):
     1147        print 'no auto search spots found'
     1148        return None
     1149    if distort > 2.:
     1150        print 'no big spots found'
     1151        return None
     1152    #use clustered points to get position & spread
     1153   
     1154    Indx = np.sort(Indx,axis=0).T
     1155    Nmax = np.max(Indx)
     1156    nums = [list(set([ind[1] for ind in Indx if ind[1]>=jnd and ind[0] == jnd])) for jnd in range(Nmax+1)]
     1157    delList = []
     1158    for Num in nums:
     1159        delList += list(np.sort(Num))[1:]
     1160    delList = list(set(delList))
     1161    Nums = [num for inum,num in enumerate(nums) if inum not in delList]
     1162    points = []
     1163    esds = []
     1164    for num in Nums:
     1165        points.append(np.mean(pts[num],axis=0))
     1166        esds.append(np.mean(np.var(pts[num],axis=0)))
     1167    peaks = (np.array(points)+np.array([.5,.5]))*pixelSize/1000.
     1168    Points = np.ones((peaks.shape[0],3))*2.
     1169    Points[:,2] += np.array(esds)*pixelSize[0]/1000.
     1170    Points[:,:2] = peaks
    11571171    Masks['Points'] = list(Points)
    11581172    return None
     1173
     1174#    rollImage = lambda rho,roll: np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1)
     1175#    print 'auto spot search'
     1176#    pixelSize = Controls['pixelSize']
     1177#    spotMask = ma.array(Image,mask=(Image<np.mean(Image)))
     1178#    indices = (-1,0,1)
     1179#    rolls = np.array([[ix,iy] for ix in indices for iy in indices])
     1180#    for roll in rolls:
     1181#        if np.any(roll):        #avoid [0,0]
     1182#            spotMask = ma.array(spotMask,mask=(spotMask-rollImage(Image,roll)<0.))
     1183#    mags = spotMask[spotMask.nonzero()]
     1184#    indx = np.transpose(spotMask.nonzero())
     1185#    nx,ny = Image.shape
     1186#    jndx = []
     1187#    for [ind,mag] in zip(indx,mags):
     1188#        if (1 < ind[0] < nx-2) and (1 < ind[1] < ny-2):
     1189#            cent = np.zeros((3,3))
     1190#            cent[1,1] = mag
     1191#            msk = np.array(Image[ind[0]-2:ind[0]+3,ind[1]-2:ind[1]+3])
     1192#            msk = ma.array(msk,mask=(msk==mag))
     1193#            if mag > 1.5*ma.mean(msk):
     1194#                jndx.append([ind[1]+.5,ind[0]+.5])
     1195#    print 'Spots found: ',len(jndx)
     1196#    jndx = np.array(jndx)
     1197#    peaks = jndx*pixelSize/1000.
     1198#    tth = GetTth(peaks.T[0],peaks.T[1],Controls)
     1199#    histth,bins = np.histogram(tth,2500)
     1200#    for shft in [-1,1]:
     1201#        histmsk = ma.array(histth,mask=-(histth-np.roll(histth,shft)<2.))
     1202#    histmsk = ma.array(histmsk,mask=(histmsk>5))
     1203#    binmsk = np.array(ma.nonzero(histmsk))
     1204#    digts = np.digitize(tth,bins)-1
     1205#    PeaksList = []
     1206#    for digt,peak in zip(digts,peaks):
     1207#        if digt in binmsk:
     1208#            PeaksList.append(peak)
     1209#    #should be able to filter out spotty Bragg rings here
     1210#    PeaksList = np.array(PeaksList)
     1211#    Points = np.ones((PeaksList.shape[0],3))
     1212#    Points[:,:2] = PeaksList
     1213#    Masks['Points'] = list(Points)
     1214#    return None
     1215       
    11591216                   
    11601217   
Note: See TracChangeset for help on using the changeset viewer.