Changeset 1163


Ignore:
Timestamp:
Dec 7, 2013 3:00:51 PM (8 years ago)
Author:
vondreele
Message:

ellipse & hyperbola fitting & plotting - much closer
still issues tho.

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIimage.py

    r1159 r1163  
    202202        R1 = (vplus+vminus)/2.                                    #major axis
    203203        zdis = (fplus-fminus)/2.
    204         X = x-x0-zdis*npsind(phi)
     204        X = x-x0-zdis*npsind(phi)       #shift obsd. ellipses to ellipse center= 0,0
    205205        Y = y-y0-zdis*npcosd(phi)
    206         XR = X*npcosd(phi)+Y*npsind(phi)
    207         YR = X*npsind(phi)+Y*npcosd(phi)
    208         return (XR/R0)**2+(YR/R1)**2-1
     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
     209        return M
    209210       
    210211    names = ['dist','det-X','det-Y','tilt','phi','dep','wave']
    … …  
    299300    sinb = sind(tilt)
    300301    d = dist+dxy
    301     if tth+tilt < 90.:      #ellipse
     302    if tth+abs(tilt) < 90.:      #ellipse
    302303        fplus = d*tanb*stth/(cosb+stth)
    303304        fminus = d*tanb*stth/(cosb-stth)
    … …  
    334335    'Needs a doc string'
    335336   
    336     from scipy.optimize import fsolve
    337     def ellip(xy,*args):
    338         azm,phi,R0,R1,A,B = args
    339         cp = cosd(phi)
    340         sp = sind(phi)
    341         x,y = xy
    342         out = []
    343         out.append(y-x*tand(azm))
    344         out.append(R0**2*((x+A)*sp-(y+B)*cp)**2+R1**2*((x+A)*cp+(y+B)*sp)**2-(R0*R1)**2)
    345         return out
    346     def hyperb(xy,*args):
    347         azm,phi,R0,R1,A,B = args
    348         cp = cosd(phi)
    349         sp = sind(phi)
    350         x,y = xy
    351         out = []
    352         out.append(y+x*tand(azm))
    353         out.append(R0**2*((x+A)*sp-(y+B)*cp)**2-R1**2*((x+A)*cp+(y+B)*sp)**2-(R0*R1)**2)
    354         return out               
    355337    elcent,phi,radii = GetEllipse(dsp,data)
     338    phi = data['rotation']-90.          #to give rotation of major axis
     339    tilt = data['tilt']
     340    dist = data['distance']
    356341    cent = data['center']
    357     phi = data['rotation']
    358     wave = data['wavelength']
    359     tilt = data['tilt']
    360     dist = data['distance']/cosd(tilt)
    361     dep = data['DetDepth']
    362     tth = 2.0*asind(wave/(2.*dsp))
    363     dxy = peneCorr(tth,dep)
     342    tth = 2.0*asind(data['wavelength']/(2.*dsp))
    364343    ttth = tand(tth)
    365     radius = (dist+dxy)*ttth
    366344    stth = sind(tth)
     345    ctth = cosd(tth)
     346    sinb = sind(tilt)
    367347    cosb = cosd(tilt)
    368     if tth+abs(tilt) < 90.:
    369         R1 = (dist+dxy)*stth*cosd(tth)*cosb/(cosb**2-stth**2)
    370         R0 = math.sqrt(R1*radius*cosb)
    371         zdis = R1*ttth*tand(tilt)
    372         A = zdis*sind(phi)
    373         B = -zdis*cosd(phi)
    374         xy0 = [radius*cosd(azm),radius*sind(azm)]
    375         xy = fsolve(ellip,xy0,args=(azm,phi,R0,R1,A,B))+cent
    376     else:
    377         R1 = (dist+dxy)*stth*cosd(tth)*cosb/(cosb**2+stth**2)
    378         R0 = math.sqrt(R1*radius*cosb)
    379         zdis = R1*ttth*tand(tilt)
    380         A = zdis*sind(phi)
    381         B = -zdis*cosd(phi)
    382         xy0 = [radius*cosd(azm),radius*sind(azm)]
    383         xy = fsolve(hyperb,xy0,args=(azm,phi,R0,R1,A,B))+cent
     348    tanb = tand(tilt)
     349    if radii[0] > 0.:
     350        fplus = dist*tanb*stth/(cosb+stth)
     351        fminus = dist*tanb*stth/(cosb-stth)
     352        zdis = 0.
     353        if tilt:
     354            zdis = tilt*(fplus-fminus)/(2.*abs(tilt))
     355        rsqplus = radii[0]**2+radii[1]**2
     356        rsqminus = radii[0]**2-radii[1]**2
     357        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))
     360        radius = (P+Q)/R
     361        xy = np.array([radius*cosd(azm),radius*sind(azm)])
     362        xy += cent
     363    else:   #hyperbola - both branches (one is way off screen!)
     364        f = dist*tanb*stth/(cosb+stth)
     365        v = dist*(tanb+tand(tth-tilt))
     366        delt = dist*stth*(1+stth*cosb)/(sinb*cosb*(stth+cosb))
     367        ecc = (v-f)/(delt-v)
     368        R = radii[1]*(ecc**2-1)/(1-ecc*cosd(azm))
     369        xy = [R*cosd(azm)-f,R*sind(azm)]
     370        xy = -np.array([xy[0]*cosd(phi)-xy[1]*sind(phi),xy[0]*sind(phi)+xy[1]*cosd(phi)])
     371        xy += cent
    384372    return xy
    385373   
    … …  
    409397    DY = np.sqrt(dx**2+dy**2-Z**2)
    410398    D = (DX**2+DY**2)/dist**2       #for geometric correction = 1/cos(2theta)^2 if tilt=0.
    411     tth = npatand(DY/DX) #depth corr not correct for tilted detector
     399    tth = npatan2d(DY,DX) #depth corr not correct for tilted detector
    412400    dsp = wave/(2.*npsind(tth/2.))
    413401    azm = (npatan2d(dy,dx)+azmthoff+720.)%360.
    … …  
    637625        print 'WARNING - selected ring was fitted as a circle'
    638626        print ' - if detector was tilted we suggest you skip this ring - WARNING'
     627        tilt = 0.01
    639628#1st estimate of dist: sample to detector normal to plane
    640629    data['distance'] = dist = radii[0]**2/(ttth*radii[1])
    … …  
    642631    zdisp = radii[1]*ttth*tand(tilt)
    643632    zdism = radii[1]*ttth*tand(-tilt)
     633    print 'zdis',zdisp,zdism
    644634#cone axis position; 2 choices. Which is right?
    645635    centp = [elcent[0]+zdisp*sind(phi),elcent[1]+zdisp*cosd(phi)]
    646636    centm = [elcent[0]+zdism*sind(phi),elcent[1]+zdism*cosd(phi)]
     637    print 'cent',centp,centm
    647638#check get same ellipse parms either way
    648639#now do next ring; estimate either way & check sum of Imax/Imin in 3x3 block around each point
    … …  
    650641    tth = 2.0*asind(wave/(2.*dsp))
    651642    ellipsep = GetEllipse2(tth,0.,dist,centp,tilt,phi)
    652     Ringp,delt,sumIp = makeRing(dsp,ellipsep,2,cutoff,scalex,scaley,self.ImageZ)
     643#    data['ellipses'].append(copy.deepcopy(ellipsep+('g',)))
     644    Ringp,delt,sumIp = makeRing(dsp,ellipsep,3,cutoff,scalex,scaley,self.ImageZ)
    653645    if len(Ringp) > 10:
    654646        outEp,errp = FitRing(Ringp,True)
    655647    else:
    656648        errp = 1e6
    657 #    print '+',ellipsep,errp
     649    print '+',sumIp,errp
    658650    ellipsem = GetEllipse2(tth,0.,dist,centm,-tilt,phi)
    659     Ringm,delt,sumIm = makeRing(dsp,ellipsem,2,cutoff,scalex,scaley,self.ImageZ)
     651#    data['ellipses'].append(copy.deepcopy(ellipsem+('r',)))
     652    Ringm,delt,sumIm = makeRing(dsp,ellipsem,3,cutoff,scalex,scaley,self.ImageZ)
    660653    if len(Ringm) > 10:
    661654        outEm,errm = FitRing(Ringm,True)
    662655    else:
    663656        errm = 1e6
    664 #    print '-',ellipsem,errm
    665     if errp < errm:
     657    print '-',sumIm,errm
     658#    if errp < errm:
     659    if sumIp > sumIm:
    666660        data['tilt'] = tilt
    667661        data['center'] = centp

  • trunk/GSASIIimgGUI.py

    r1158 r1163  
    464464       
    465465        def OnLRazim(event):
    466             Lazm =int(G2frame.Lazim.GetValue())
     466            Lazm = int(G2frame.Lazim.GetValue())%360
     467            Razm = int(G2frame.Razim.GetValue())%360
     468            if Lazm > Razm:
     469                Razm += 360
    467470            if data['fullIntegrate']:
    468                G2frame.Razim.SetValue("%6d" % (Lazm+360))
    469             Razm = int(G2frame.Razim.GetValue())
    470             if Lazm > Razm:
    471                 Lazm -= 360
     471                Razm = Lazm+360
     472            G2frame.Lazim.SetValue("%6d" % (Lazm))
     473            G2frame.Razim.SetValue("%6d" % (Razm))
    472474            data['LRazimuth'] = [Lazm,Razm]
    473475            G2plt.PlotExposedImage(G2frame,event=event)

  • trunk/GSASIIplot.py

    r1153 r1163  
    23232323                    elif 'line3' in itemPicked:
    23242324                        Data['LRazimuth'][0] = int(azm)
    2325                         if Data['fullIntegrate']:
    2326                             Data['LRazimuth'][1] = Data['LRazimuth'][0]+360
    23272325                    elif 'line4' in itemPicked and not Data['fullIntegrate']:
    23282326                        Data['LRazimuth'][1] = int(azm)
    2329                        
     2327                   
     2328                    Data['LRazimuth'][0] %= 360
     2329                    Data['LRazimuth'][1] %= 360
    23302330                    if Data['LRazimuth'][0] > Data['LRazimuth'][1]:
    2331                         Data['LRazimuth'][0] -= 360
    2332                        
    2333                     azLim = np.array(Data['LRazimuth'])   
    2334                     if np.any(azLim>360):
    2335                         azLim -= 360
    2336                         Data['LRazimuth'] = list(azLim)
     2331                        Data['LRazimuth'][1] += 360                       
     2332                    if Data['fullIntegrate']:
     2333                        Data['LRazimuth'][1] = Data['LRazimuth'][0]+360
    23372334                       
    23382335                    if  Data['IOtth'][0] > Data['IOtth'][1]:
    … …  
    24942491                xyI = []
    24952492                for azm in Azm:
    2496                     xyI.append(G2img.GetDetectorXY(dspI,azm,Data))      #what about hyperbola?
    2497                 xyI = np.array(xyI)
    2498                 arcxI,arcyI = xyI.T
    2499                 Plot.plot(arcxI,arcyI,picker=3)
     2493                    xy = G2img.GetDetectorXY(dspI,azm,Data)
     2494                    if np.any(xy):
     2495                        xyI.append(xy)
     2496                if len(xyI):
     2497                    xyI = np.array(xyI)
     2498                    arcxI,arcyI = xyI.T
     2499                    Plot.plot(arcxI,arcyI,picker=3)
    25002500            if ellO:
    25012501                xyO = []
    25022502                for azm in Azm:
    2503                     xyO.append(G2img.GetDetectorXY(dspO,azm,Data))      #what about hyperbola?
    2504                 xyO = np.array(xyO)
    2505                 arcxO,arcyO = xyO.T
    2506                 Plot.plot(arcxO,arcyO,picker=3)
     2503                    xy = G2img.GetDetectorXY(dspO,azm,Data)
     2504                    if np.any(xy):
     2505                        xyO.append(xy)
     2506                if len(xyO):
     2507                    xyO = np.array(xyO)
     2508                    arcxO,arcyO = xyO.T               
     2509                    Plot.plot(arcxO,arcyO,picker=3)
    25072510            if ellO and ellI:
    25082511                Plot.plot([arcxI[0],arcxO[0]],[arcyI[0],arcyO[0]],picker=3)
Note: See TracChangeset for help on using the changeset viewer.