Changeset 1165 for trunk/GSASIIimage.py


Ignore:
Timestamp:
Dec 14, 2013 8:54:05 AM (9 years ago)
Author:
vondreele
Message:

image calibration much better now - works on transposed data

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIimage.py

    r1164 r1165  
    8181    return np.roll(np.roll(M,Axis,axis=0),Axis,axis=1)
    8282                   
    83 def FitRing(ring,delta):
    84     'Needs a doc string'
    85     parms = []
    86     if delta:
    87         err,parms = FitEllipse(ring)
    88         errc,parmsc = FitCircle(ring)
    89         errc = errc[0]/(len(ring)*parmsc[2][0]**2)
    90         if not parms or errc < .1:
    91             parms = parmsc
    92     else:
    93         err,parms = FitCircle(ring)
    94     return parms,err
    95        
    96 def FitCircle(ring):
    97     'Needs a doc string'
    98    
    99     def makeParmsCircle(B):
    100         cent = [-B[0]/2,-B[1]/2]
    101         phi = 0.
    102         sr1 = sr2 = math.sqrt(cent[0]**2+cent[1]**2-B[2])
    103         return cent,phi,[sr1,sr2]
    104        
    105     ring = np.array(ring)
    106     x = np.asarray(ring.T[0])
    107     y = np.asarray(ring.T[1])
    108    
    109     M = np.array((x,y,np.ones_like(x)))
    110     B = np.array(-(x**2+y**2))
    111     result = nl.lstsq(M.T,B)
    112     return result[1],makeParmsCircle(result[0])
    113        
    114 def FitEllipse(ring):
    115     'Needs a doc string'
    116            
    117     def makeParmsEllipse(B):
    118         det = 4.*(1.-B[0]**2)-B[1]**2
    119         if det < 0.:
    120             print 'hyperbola!'
    121             print B
    122             return 0
    123         elif det == 0.:
    124             print 'parabola!'
    125             print B
    126             return 0
    127         cent = [(B[1]*B[3]-2.*(1.-B[0])*B[2])/det, \
    128             (B[1]*B[2]-2.*(1.+B[0])*B[3])/det]
    129         phi = 0.5*atand(0.5*B[1]/B[0])
    130        
    131         a = (1.+B[0])/cosd(2*phi)
    132         b = 2.-a
    133         f = (1.+B[0])*cent[0]**2+(1.-B[0])*cent[1]**2+B[1]*cent[0]*cent[1]-B[4]
    134         if f/a < 0 or f/b < 0:
    135             return 0
    136         sr1 = math.sqrt(f/a)
    137         sr2 = math.sqrt(f/b)
    138         if sr1 > sr2:
    139             sr1,sr2 = sr2,sr1
    140             phi -= 90.
    141             if phi < -90.:
    142                 phi += 180.
    143         return cent,phi,[sr1,sr2]
    144                
    145     ring = np.array(ring)
    146     x = np.asarray(ring.T[0])
    147     y = np.asarray(ring.T[1])
    148     M = np.array((x**2-y**2,x*y,x,y,np.ones_like(x)))
    149     B = np.array(-(x**2+y**2))
    150     bb,err = fel.fellipse(len(x),x,y,1.E-7)
    151 #    print nl.lstsq(M.T,B)[0]
    152 #    print bb
    153     return err,makeParmsEllipse(bb)
    154    
     83def FitEllipse(xy):
     84   
     85    def ellipse_center(p):
     86        ''' gives ellipse center coordinates
     87        '''
     88        b,c,d,f,g,a = p[1]/2, p[2], p[3]/2, p[4]/2, p[5], p[0]
     89        num = b*b-a*c
     90        x0=(c*d-b*f)/num
     91        y0=(a*f-b*d)/num
     92        return np.array([x0,y0])
     93   
     94    def ellipse_angle_of_rotation( p ):
     95        ''' gives rotation of ellipse major axis from x-axis
     96        range will be -90 to 90 deg
     97        '''
     98        b,c,d,f,g,a = p[1]/2, p[2], p[3]/2, p[4]/2, p[5], p[0]
     99        return 0.5*npatand(2*b/(a-c))
     100   
     101    def ellipse_axis_length( p ):
     102        ''' gives ellipse radii in [minor,major] order
     103        '''
     104        b,c,d,f,g,a = p[1]/2, p[2], p[3]/2, p[4]/2, p[5], p[0]
     105        up = 2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g)
     106        down1=(b*b-a*c)*( (c-a)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
     107        down2=(b*b-a*c)*( (a-c)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
     108        res1=np.sqrt(up/down1)
     109        res2=np.sqrt(up/down2)
     110        return np.array([ res2,res1])
     111   
     112    xy = np.array(xy)
     113    x = np.asarray(xy.T[0])[:,np.newaxis]
     114    y = np.asarray(xy.T[1])[:,np.newaxis]
     115    D =  np.hstack((x*x, x*y, y*y, x, y, np.ones_like(x)))
     116    S = np.dot(D.T,D)
     117    C = np.zeros([6,6])
     118    C[0,2] = C[2,0] = 2; C[1,1] = -1
     119    E, V =  nl.eig(np.dot(nl.inv(S), C))
     120    n = np.argmax(np.abs(E))
     121    a = V[:,n]
     122    cent = ellipse_center(a)
     123    phi = ellipse_angle_of_rotation(a)
     124    radii = ellipse_axis_length(a)
     125    phi += 90.
     126    if radii[0] > radii[1]:
     127        radii = [radii[1],radii[0]]
     128        phi -= 90.
     129    return cent,phi,radii
     130
    155131def FitDetector(rings,varyList,parmDict,Print=True):
    156132    'Needs a doc string'
    … …  
    163139        for name,fmt,value,sig in ValSig:
    164140            ptlbls += "%s" % (name.rjust(12))
    165             if name == 'rotate':
    166                 ptstr += fmt % (value-90.)      #kluge to get rotation from vertical - see GSASIIimgGUI
     141            if name == 'phi':
     142                ptstr += fmt % (value%360.)
    167143            else:
    168144                ptstr += fmt % (value)
    … …  
    179155        wave = parmDict['wave']
    180156        if 'dep' in varyList:
    181             dist,x0,y0,tilt,phi,dep = B[:6]
     157            dist,x0,y0,tilt,chi,dep = B[:6]
    182158        else:
    183             dist,x0,y0,tilt,phi = B[:5]
     159            dist,x0,y0,tilt,chi = B[:5]
    184160            dep = parmDict['dep']
    185161        if 'wave' in varyList:
    186162            wave = B[-1]
     163        phi = chi-90.               #get rotation of major axis from tilt axis
    187164        tth = 2.0*npasind(wave/(2.*dsp))
    188165        dxy = peneCorr(tth,dep)
    … …  
    202179        R1 = (vplus+vminus)/2.                                    #major axis
    203180        zdis = (fplus-fminus)/2.
    204         X = x-x0-zdis*npsind(phi)       #shift obsd. ellipses to ellipse center= 0,0
    205         Y = y-y0-zdis*npcosd(phi)
    206         XR = X*npcosd(phi)+Y*npsind(phi)    #rotate by phi to put major axis along x
    207         YR = -X*npsind(phi)+Y*npcosd(phi)
    208         M = (XR/R0)**2+(YR/R1)**2-1         #=0 for all points exactly on ellipse
     181        Robs = np.sqrt((x-x0)**2+(y-y0)**2)
     182        phi0 = npatan2d(y-y0,x-x0)
     183        rsqplus = R0**2+R1**2
     184        rsqminus = R0**2-R1**2
     185        R = rsqminus*npcosd(2.*phi0-2.*phi)+rsqplus
     186        Q = np.sqrt(2.)*R0*R1*np.sqrt(R-2.*zdis**2*npsind(phi0-phi)**2)
     187        P = 2.*R0**2*zdis*npcosd(phi0-phi)
     188        Rcalc = (P+Q)/R
     189        M = Robs-Rcalc
    209190        return M
    210191       
    … …  
    266247        x = radii[0]*cosd(a)
    267248        y = radii[1]*sind(a)
    268         X = (cphi*x-sphi*y+cent[0])*scalex      #convert mm to pixels
    269         Y = (sphi*x+cphi*y+cent[1])*scaley
     249        X = (cphi*x+sphi*y+cent[0])*scalex      #convert mm to pixels
     250        Y = (-sphi*x+cphi*y+cent[1])*scaley
    270251        X,Y,I,J = ImageLocalMax(image,pix,X,Y)
    271252        if I and J and I/J > reject:
    … …  
    279260            if [X,Y,dsp] not in ring:
    280261                ring.append([X,Y,dsp])
    281     delt = amax-amin
    282262    if len(ring) < 10:             #want more than 10 deg
    283         return [],(delt > 90)
    284     return ring,(delt > 90)
     263        ring = []
     264    return ring
    285265   
    286266def GetEllipse2(tth,dxy,dist,cent,tilt,phi):
    … …  
    316296        radii[1] = eps*(delt-f)/(eps**2-1.)                             #major axis
    317297        zdis = f+radii[1]*eps
    318     elcent = [cent[0]+zdis*sind(phi),cent[1]+zdis*cosd(phi)]
     298#NB: zdis is || to major axis & phi is rotation of minor axis
     299#thus shift from beam to ellipse center is [Z*sin(phi),-Z*cos(phi)]
     300    elcent = [cent[0]+zdis*sind(phi),cent[1]-zdis*cosd(phi)]
    319301    return elcent,phi,radii
    320302   
    … …  
    350332        fplus = dist*tanb*stth/(cosb+stth)
    351333        fminus = dist*tanb*stth/(cosb-stth)
    352         zdis = 0.
    353         if tilt:
    354             zdis = tilt*(fplus-fminus)/(2.*abs(tilt))
     334        zdis = (fplus-fminus)/2.
    355335        rsqplus = radii[0]**2+radii[1]**2
    356336        rsqminus = radii[0]**2-radii[1]**2
    357337        R = rsqminus*cosd(2.*azm-2.*phi)+rsqplus
    358         Q = np.sqrt(2.)*radii[0]*radii[1]*np.sqrt(R-2.*zdis**2*sind(azm)**2)
    359         P = zdis*(rsqminus*cosd(azm-2.*phi)+rsqplus*cosd(azm))
     338        Q = np.sqrt(2.)*radii[0]*radii[1]*np.sqrt(R-2.*zdis**2*sind(azm-phi)**2)
     339        P = 2.*radii[0]**2*zdis*cosd(azm-phi)
    360340        radius = (P+Q)/R
    361341        xy = np.array([radius*cosd(azm),radius*sind(azm)])
    … …  
    363343    else:   #hyperbola - both branches (one is way off screen!)
    364344        f = dist*tanb*stth/(cosb+stth)
    365         v = dist*(tanb+tand(tth-tilt))
     345        v = dist*(tanb+tand(tth-abs(tilt)))
    366346        delt = dist*stth*(1+stth*cosb)/(sinb*cosb*(stth+cosb))
    367347        ecc = (v-f)/(delt-v)
    368348        R = radii[1]*(ecc**2-1)/(1-ecc*cosd(azm))
    369         xy = [R*cosd(azm)-f,R*sind(azm)]
     349        if tilt > 0.:
     350            offset = 2.*radii[1]*ecc+f      #select other branch
     351            xy = [-R*cosd(azm)-offset,R*sind(azm)]
     352        else:
     353            offset = -f
     354            xy = [-R*cosd(azm)-offset,-R*sind(azm)]
    370355        xy = -np.array([xy[0]*cosd(phi)-xy[1]*sind(phi),xy[0]*sind(phi)+xy[1]*cosd(phi)])
    371356        xy += cent
    … …  
    525510            break
    526511        ellipse = GetEllipse(dsp,data)
    527         Ring,delt = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,ma.array(self.ImageZ,mask=tam))
     512        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,ma.array(self.ImageZ,mask=tam))
    528513        if Ring:
    529514            if iH >= skip:
    … …  
    572557    #fit start points on inner ring
    573558    data['ellipses'] = []
    574     outE,err = FitRing(ring,True)
     559    data['rings'] = []
     560    outE = FitEllipse(ring)
    575561    fmt  = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f'
    576562    fmt2 = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f, chi^2: %.3f, N: %d'
    … …  
    582568       
    583569    #setup 360 points on that ring for "good" fit
    584     Ring,delt = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
     570    data['ellipses'].append(ellipse[:]+('g',))
     571    Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
    585572    if Ring:
    586         ellipse,err = FitRing(Ring,delt)
    587         Ring,delt = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)    #do again
    588         ellipse,err = FitRing(Ring,delt)
     573        ellipse = FitEllipse(ring)
     574        Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)    #do again
     575        ellipse = FitEllipse(ring)
    589576    else:
    590577        print '1st ring not sufficiently complete to proceed'
    591578        return False
    592     print fmt2%('inner ring:    ',ellipse[0][0],ellipse[0][1],ellipse[1],ellipse[2][0],ellipse[2][1],err,len(Ring))     #cent,phi,radii
     579    print fmt2%('inner ring:    ',ellipse[0][0],ellipse[0][1],ellipse[1],ellipse[2][0],ellipse[2][1],0.,len(Ring))     #cent,phi,radii
    593580    data['ellipses'].append(ellipse[:]+('r',))
    594581    data['rings'].append(np.array(Ring))
    595582    G2plt.PlotImage(self,newImage=True)
    596583   
    597     #setup for calibration
     584#setup for calibration
    598585    data['rings'] = []
    599     data['ellipses'] = []
    600586    if not data['calibrant']:
    601587        print 'no calibration material selected'
    … …  
    617603    dsp = HKL[0][3]
    618604    tth = 2.0*asind(wave/(2.*dsp))
    619     Ring0 = makeRing(dsp,ellipse,3,cutoff,scalex,scaley,self.ImageZ)[0]
     605    Ring0 = makeRing(dsp,ellipse,3,cutoff,scalex,scaley,self.ImageZ)
    620606    ttth = nptand(tth)
    621607    stth = npsind(tth)
    … …  
    626612        print 'WARNING - selected ring was fitted as a circle'
    627613        print ' - if detector was tilted we suggest you skip this ring - WARNING'
    628         tilt = 0.01
    629614#1st estimate of dist: sample to detector normal to plane
    630615    data['distance'] = dist = radii[0]**2/(ttth*radii[1])
    … …  
    632617    zdisp = radii[1]*ttth*tand(tilt)
    633618    zdism = radii[1]*ttth*tand(-tilt)
    634 #cone axis position; 2 choices. Which is right?
    635     centp = [elcent[0]+zdisp*sind(phi),elcent[1]+zdisp*cosd(phi)]
    636     centm = [elcent[0]+zdism*sind(phi),elcent[1]+zdism*cosd(phi)]
     619#cone axis position; 2 choices. Which is right?     
     620#NB: zdisp is || to major axis & phi is rotation of minor axis
     621#thus shift from beam to ellipse center is [Z*sin(phi),-Z*cos(phi)]
     622    centp = [elcent[0]+zdisp*sind(phi),elcent[1]-zdisp*cosd(phi)]
     623    centm = [elcent[0]+zdism*sind(phi),elcent[1]-zdism*cosd(phi)]
    637624#check get same ellipse parms either way
    638625#now do next ring; estimate either way & do a FitDetector each way; best fit is correct one
    … …  
    640627    tth = 2.0*asind(wave/(2.*dsp))
    641628    ellipsep = GetEllipse2(tth,0.,dist,centp,tilt,phi)
    642     ellipsem = GetEllipse2(tth,0.,dist,centm,-tilt,phi)
    643     Ringp = makeRing(dsp,ellipsep,3,cutoff,scalex,scaley,self.ImageZ)[0]
    644     Ringm = makeRing(dsp,ellipsem,3,cutoff,scalex,scaley,self.ImageZ)[0]
     629    print fmt%('plus ellipse :',ellipsep[0][0],ellipsep[0][1],ellipsep[1],ellipsep[2][0],ellipsep[2][1])
     630    Ringp = makeRing(dsp,ellipsep,3,cutoff,scalex,scaley,self.ImageZ)
    645631    parmDict = {'dist':dist,'det-X':centp[0],'det-Y':centp[1],
    646632        'tilt':tilt,'phi':phi,'wave':wave,'dep':0.0}
    647633    varyList = ['dist','det-X','det-Y','tilt','phi']
    648634    if len(Ringp) > 10:
    649         chip,chisqp = FitDetector(np.array(Ring0+Ringp),varyList,parmDict,False)
     635        chip,chisqp = FitDetector(np.array(Ring0+Ringp),varyList,parmDict,True)
     636        tiltp = parmDict['tilt']
     637        phip = parmDict['phi']
     638        centp = [parmDict['det-X'],parmDict['det-Y']]
    650639    else:
    651640        chip = 1e6
     641    ellipsem = GetEllipse2(tth,0.,dist,centm,-tilt,phi)
     642    print fmt%('minus ellipse:',ellipsem[0][0],ellipsem[0][1],ellipsem[1],ellipsem[2][0],ellipsem[2][1])
     643    Ringm = makeRing(dsp,ellipsem,3,cutoff,scalex,scaley,self.ImageZ)
    652644    if len(Ringm) > 10:
    653645        parmDict['tilt'] *= -1
    654         chim,chisqm = FitDetector(np.array(Ring0+Ringm),varyList,parmDict,False)
     646        chim,chisqm = FitDetector(np.array(Ring0+Ringm),varyList,parmDict,True)
     647        tiltm = parmDict['tilt']
     648        phim = parmDict['phi']
     649        centm = [parmDict['det-X'],parmDict['det-Y']]
    655650    else:
    656651        chim = 1e6
    657652    if chip < chim:
    658         data['tilt'] = tilt
     653        data['tilt'] = tiltp
    659654        data['center'] = centp
     655        data['rotation'] = phip
    660656    else:
    661         data['tilt'] = -tilt
     657        data['tilt'] = tiltm
    662658        data['center'] = centm
    663     data['rotation'] = phi
     659        data['rotation'] = phim
     660    data['ellipses'].append(ellipsep[:]+('b',))
     661    data['rings'].append(np.array(Ringp))
     662    data['ellipses'].append(ellipsem[:]+('r',))
     663    data['rings'].append(np.array(Ringm))
     664    G2plt.PlotImage(self,newImage=True)
    664665    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
    665666        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
    … …  
    668669        varyList.append('dep')
    669670    data['rings'] = []
     671    data['ellipses'] = []
    670672    for i,H in enumerate(HKL):
    671673        dsp = H[3]
    … …  
    678680        data['ellipses'].append(copy.deepcopy(ellipse+('g',)))
    679681        print fmt%('predicted ellipse:',elcent[0],elcent[1],phi,radii[0],radii[1])
    680         Ring,delt = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
     682        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
    681683        if Ring:
    682684            data['rings'].append(np.array(Ring))
    … …  
    686688                data['distance'] = parmDict['dist']
    687689                data['center'] = [parmDict['det-X'],parmDict['det-Y']]
    688                 data['rotation'] = np.mod(parmDict['phi']+720.,360.0)
     690                data['rotation'] = parmDict['phi']
    689691                data['tilt'] = parmDict['tilt']
    690                 if data['rotation'] >180.:
    691                     data['rotation'] -= 180.
    692                     data['tilt'] *= -1.
    693692                data['DetDepth'] = parmDict['dep']
    694693                elcent,phi,radii = ellipse = GetEllipse(dsp,data)
    695694                print fmt2%('fitted ellipse:   ',elcent[0],elcent[1],phi,radii[0],radii[1],chisq,len(rings))
    696695            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
     696#            G2plt.PlotImage(self,newImage=True)
    697697        else:
    698698            print 'insufficient number of points in this ellipse to fit'
    … …  
    707707        data['distance'] = parmDict['dist']
    708708        data['center'] = [parmDict['det-X'],parmDict['det-Y']]
    709         data['rotation'] = np.mod(parmDict['phi']+720.,360.0)
     709        data['rotation'] = parmDict['phi']
    710710        data['tilt'] = parmDict['tilt']
    711         if data['rotation'] >180.:
    712             data['rotation'] -= 180.
    713             data['tilt'] *= -1.
    714711        data['DetDepth'] = parmDict['dep']
    715712    for H in HKL[:N]:
    … …  
    863860    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
    864861        ellipse = GetEllipse(ring['Dset'],StaControls)
    865         Ring,delt = makeRing(ring['Dset'],ellipse,ring['pixLimit'],ring['cutoff'],scalex,scaley,Image)
     862        Ring = makeRing(ring['Dset'],ellipse,ring['pixLimit'],ring['cutoff'],scalex,scaley,Image)
    866863        Ring = np.array(Ring).T
    867864        ring['ImxyObs'] = np.array(Ring[:2])      #need to apply masks to this to eliminate bad points
Note: See TracChangeset for help on using the changeset viewer.