Skip to content

Commit abaac4e

Browse files
vondreelevondreele
committed
fix hyperbola drawing, etc.
1 parent 056a72b commit abaac4e

File tree

4 files changed

+49
-157
lines changed

4 files changed

+49
-157
lines changed

‎GSASII/GSASIIimage.py‎

Lines changed: 33 additions & 141 deletions
Original file line numberDiff line numberDiff line change
@@ -42,7 +42,7 @@
4242
npatand = lambda x: 180.*np.arctan(x)/np.pi
4343
npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
4444
nxs = np.newaxis
45-
debug = True
45+
debug = False
4646

4747
def peneCorr(tth,dep,dist):
4848
'Compute empirical position correction due to detector absorption'
@@ -457,94 +457,40 @@ def GetEllipse(dsp,data):
457457
dxy = peneCorr(tth,dep,dist)
458458
return GetEllipse2(tth,dxy,dist,cent,tilt,phi)
459459

460-
# def GetDetectorXY(dsp,azm,data):
461-
# '''Get detector x,y position from d-spacing (dsp), azimuth (azm,deg)
462-
# & image controls dictionary (data) - new version
463-
# it seems to be only used in plotting & wrong
464-
# '''
465-
# def LinePlaneCollision(planeNormal, planePoint, rayDirection, rayPoint, epsilon=1e-6):
466-
467-
# ndotu = planeNormal.dot(rayDirection)
468-
# if ndotu < epsilon:
469-
# return None
470-
# w = rayPoint - planePoint
471-
# si = -planeNormal.dot(w) / ndotu
472-
# Psi = w + si * rayDirection + planePoint
473-
# return Psi
474-
475-
# dist = data['distance']
476-
# cent = data['center']
477-
# phi = data['rotation']-90. #to give rotation of major axis
478-
# T = makeMat(data['tilt'],0) #rotate about X
479-
# R = makeMat(phi,2) #rotate about Z
480-
# MN = np.inner(R,np.inner(R,T))
481-
# iMN= nl.inv(MN)
482-
# tth = 2.0*npasind(data['wavelength']/(2.*dsp))
483-
# vect = np.array([npsind(tth)*npcosd(azm-phi),npsind(tth)*npsind(azm-phi),npcosd(tth)])
484-
# dxyz0 = np.inner(np.array([0.,0.,1.0]),MN) #tilt detector normal
485-
# dxyz0 += np.array([0.,0.,dist]) #translate to distance
486-
# dxyz0 = np.inner(dxyz0,makeMat(data['det2theta'],1).T) #rotate on 2-theta
487-
# dxyz1 = np.inner(np.array([cent[0],cent[1],0.]),MN) #tilt detector cent
488-
# dxyz1 += np.array([0.,0.,dist]) #translate to distance
489-
# dxyz1 = np.inner(dxyz1,makeMat(data['det2theta'],1).T) #rotate on 2-theta
490-
# xyz = LinePlaneCollision(dxyz1,dxyz0,vect,dist*vect)
491-
# if xyz is None:
492-
# return np.zeros(2)
493-
# # return None
494-
# xyz = np.inner(xyz,makeMat(data['det2theta'],1).T)
495-
# xyz -= np.array([0.,0.,dist]) #translate back
496-
# xyz = np.inner(xyz,iMN)
497-
# return np.squeeze(xyz)[:2]+cent
498-
499-
def GetDetectorXY2(dsp,azm,data):
460+
def GetDetectorXY(dsp,azm,data):
500461
'''Get detector x,y position from d-spacing (dsp), azimuth (azm,deg)
501-
& image controls dictionary (data)
502-
it seems to be only used in plotting
462+
& image controls dictionary (data) - new version
503463
'''
504-
elcent,phi,radii = GetEllipse(dsp,data)
505-
phi = data['rotation']-90. #to give rotation of major axis
506-
tilt = data['tilt']
507-
dist = data['distance']
464+
465+
def LinePlaneCollision(planeNormal, planePoint, rayDirection, rayPoint, epsilon=1e-6):
466+
467+
ndotu = planeNormal.dot(rayDirection)
468+
if ndotu < epsilon:
469+
return None
470+
w = rayPoint - planePoint
471+
si = -planeNormal.dot(w) / ndotu
472+
Psi = w + si * rayDirection + planePoint
473+
return Psi
474+
475+
dist = data['distance']/npcosd(data['tilt']) #sample-beam intersection point on detector plane
508476
cent = data['center']
509-
tth = 2.0*asind(data['wavelength']/(2.*dsp))
510-
stth = sind(tth)
511-
cosb = cosd(tilt)
512-
if radii[0] > 0.: #ellipse - correct!
513-
sinb = sind(tilt)
514-
tanb = tand(tilt)
515-
fplus = dist*tanb*stth/(cosb+stth)
516-
fminus = dist*tanb*stth/(cosb-stth)
517-
zdis = (fplus-fminus)/2.
518-
rsqplus = radii[0]**2+radii[1]**2
519-
rsqminus = radii[0]**2-radii[1]**2
520-
R = rsqminus*cosd(2.*azm-2.*phi)+rsqplus
521-
Q = np.sqrt(2.)*radii[0]*radii[1]*np.sqrt(R-2.*zdis**2*sind(azm-phi)**2)
522-
if radii[0] <= 0.:
523-
zdis *= -1.
524-
P = 2.*radii[0]**2*zdis*cosd(azm-phi)
525-
radius = (P+Q)/R
526-
xy = np.array([radius*cosd(azm),radius*sind(azm)])
527-
xy += cent
528-
else: #hyperbola - both branches (one is way off screen!) - not correct
529-
sinb = abs(sind(tilt))
530-
tanb = abs(tand(tilt))
531-
f = dist*tanb*stth/(cosb+stth)
532-
v = dist*(tanb+tand(tth-abs(tilt)))
533-
delt = dist*stth*(1+stth*cosb)/(sinb*cosb*(stth+cosb))
534-
ecc = (v-f)/(delt-v)
535-
R = radii[1]*(ecc**2-1)/(1-ecc*cosd(azm-phi))
536-
if tilt > 0.:
537-
offset = 2.*radii[1]*ecc+f #select other branch
538-
xy = [-R*cosd(azm-phi)-offset,-R*sind(azm-phi)]
539-
else:
540-
offset = -f
541-
xy = [-R*cosd(azm-phi)-offset,2.*R*sind(azm-phi)]
542-
print(azm,ecc,R,offset,cosd(azm-phi),sind(azm-phi),xy)
543-
xy = -np.array([xy[0]*cosd(phi)+xy[1]*sind(phi),xy[0]*sind(phi)-xy[1]*cosd(phi)])
544-
xy += cent
545-
if data['det2theta']:
546-
xy[0] += dist*nptand(data['det2theta']+data['tilt']*npsind(data['rotation']))
547-
return xy
477+
phi = data['rotation']
478+
T = makeMat(data['tilt'],0) #rotate about X
479+
R = makeMat(phi,2) #rotate about Z
480+
MN = np.inner(R,np.inner(R,T))
481+
iMN= nl.inv(MN)
482+
tth = 2.0*npasind(data['wavelength']/(2.*dsp))
483+
vect = np.array([npsind(tth)*npcosd(azm),npsind(tth)*npsind(azm),npcosd(tth)])
484+
dxyzN = np.inner(np.array([0.,0.,1.0]),MN) #tilt detector normal
485+
dxyzO = np.array([0.,0.,dist]) #translate to distance
486+
dorig = np.zeros(3)
487+
xyz = LinePlaneCollision(dxyzN,dxyzO,vect,dorig)
488+
if xyz is None:
489+
return np.zeros(2)
490+
xyz = np.inner(xyz,makeMat(data['det2theta'],1).T)
491+
xyz -= np.array([0.,0.,dist]) #translate back
492+
xyz = np.inner(xyz,iMN)
493+
return np.squeeze(xyz)[:2]+cent
548494

549495
def GetTthAzmDsp2(x,y,data): #expensive
550496
'''Computes a 2theta, etc. from a detector position and calibration constants - checked
@@ -644,28 +590,6 @@ def GetAzm(x,y,data):
644590
else:
645591
return GetTthAzmDsp2(x,y,data)[1]
646592
# these two are used only for integration & finding pixel masks
647-
# def GetTthAzmG2(x,y,data):
648-
# '''Give 2-theta, azimuth & geometric corr. values for detector x,y position;
649-
# calibration info in data - only used in integration for detector 2-theta = 0
650-
# '''
651-
# tilt = data['tilt']
652-
# dist = data['distance']/npcosd(tilt)
653-
# MN = -np.inner(makeMat(data['rotation'],2),makeMat(tilt,0))
654-
# dx = x-data['center'][0]
655-
# dy = y-data['center'][1]
656-
# dz = np.dot(np.dstack([dx.T,dy.T,np.zeros_like(dx.T)]),MN).T[2]
657-
# xyZ = dx**2+dy**2-dz**2
658-
# tth0 = npatand(np.sqrt(xyZ)/(dist-dz))
659-
# dzp = peneCorr(tth0,data['DetDepth'],dist)
660-
# tth = npatan2d(np.sqrt(xyZ),dist-dz+dzp)
661-
# azm = (npatan2d(dy,dx)+data['azmthOff']+720.)%360.
662-
# # G-calculation - use Law of sines
663-
# distm = data['distance']/1000.0
664-
# sinB2 = np.minimum(np.ones_like(tth),(data['distance']*npsind(tth))**2/(dx**2+dy**2))
665-
# #sinB2 = (data['distance']*npsind(tth))**2/(dx**2+dy**2)
666-
# C = 180.-tth-npacosd(np.sqrt(1.- sinB2))
667-
# G = distm**2*sinB2/npsind(C)**2
668-
# return tth,azm,G
669593

670594
def GetTthAzmG(x,y,data):
671595
'''Give 2-theta, azimuth & geometric corr. values for detector x,y position;
@@ -705,31 +629,11 @@ def costth(xyz,d0):
705629
azm = (npatan2d(dxyz[:,:,1],dxyz[:,:,0])+data['azmthOff']+720.)%360.
706630
# new G calculation - parallax & distance (in meters)
707631
G = 1./ctth0*np.sum(dxyz0**2,axis=2)*1.e-6 #parallax*distance^2; convert to m^2 - correct!
708-
# old G-calculation
709-
# x0 = data['distance']*nptand(tilt)
710-
# x0x = x0*npcosd(data['rotation'])
711-
# x0y = x0*npsind(data['rotation'])
712-
# distsq = data['distance']**2
713-
# G = ((dx-x0x)**2+(dy-x0y)**2+distsq)/distsq #for geometric correction = 1/cos(2theta)^2 if tilt=0.
714-
# G *= ctth0
715-
# # G-calculation - use Law of sines - wrong for nonzero det2theta!
716-
# distm = data['distance']/1000.0
717-
# sinB2 = np.minimum(np.ones_like(tth),
718-
# (data['distance']*npsind(tth))**2/(dx**2+dy**2))
719-
# C = 180.-tth-npacosd(np.sqrt(1.- sinB2))
720-
# G = distm**2*sinB2/npsind(C)**2
721632
return tth,azm,G
722633

723634
def meanAzm(a,b):
724635
AZM = lambda a,b: npacosd(0.5*(npsind(2.*b)-npsind(2.*a))/(np.pi*(b-a)/180.))/2.
725636
azm = AZM(a,b)
726-
# quad = int((a+b)/180.)
727-
# if quad == 1:
728-
# azm = 180.-azm
729-
# elif quad == 2:
730-
# azm += 180.
731-
# elif quad == 3:
732-
# azm = 360-azm
733637
return azm
734638

735639
def ImageCompress(image,scale):
@@ -761,7 +665,7 @@ def GetLineScan(image,data):
761665
Tx = np.array([tth for tth in np.linspace(LUtth[0],LUtth[1],numChans+1)])
762666
Ty = np.zeros_like(Tx)
763667
dsp = wave/(2.0*npsind(Tx/2.0))
764-
xy = [GetDetectorXY2(d,azm,data) for d in dsp]
668+
xy = [GetDetectorXY(d,azm,data) for d in dsp]
765669
xy = np.array(xy).T
766670
xy[1] *= scalex
767671
xy[0] *= scaley
@@ -833,10 +737,6 @@ def CalcRings(G2frame,ImageZ,data,masks):
833737
for iH,H in enumerate(HKL):
834738
if debug: print (H)
835739
dsp = H[3]
836-
tth = 2.0*asind(wave/(2.*dsp))
837-
# if tth+abs(data['tilt']) > 90.:
838-
# G2fil.G2Print ('next line is a hyperbola - search stopped')
839-
# break
840740
ellipse = GetEllipse(dsp,data)
841741
if iH not in absent and iH >= skip:
842742
Ring = makeRing(dsp,ellipse,0,-1.,scalex,scaley,ma.array(ImageZ,mask=tam))[0]
@@ -903,18 +803,13 @@ def ImageRecalibrate(G2frame,ImageZ,data,masks,getRingsOnly=False):
903803
'setdist':data.get('setdist',data['distance']),
904804
'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
905805
Found = False
906-
wave = data['wavelength']
907806
frame = masks['Frames']
908807
tam = ma.make_mask_none(ImageZ.shape)
909808
if frame:
910809
tam = ma.mask_or(tam,ma.make_mask(np.abs(polymask(data,frame)-255)))
911810
for iH,H in enumerate(HKL):
912811
if debug: print (H)
913812
dsp = H[3]
914-
tth = 2.0*asind(wave/(2.*dsp))
915-
# if tth+abs(data['tilt']) > 90.:
916-
# G2fil.G2Print ('next line is a hyperbola - search stopped')
917-
# break
918813
ellipse = GetEllipse(dsp,data)
919814
if iH not in absent and iH >= skip:
920815
Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,ma.array(ImageZ,mask=tam))[0]
@@ -1128,9 +1023,6 @@ def ImageCalibrate(G2frame,data):
11281023
for i,H in enumerate(HKL):
11291024
dsp = H[3]
11301025
tth = 2.0*asind(wave/(2.*dsp))
1131-
# if tth+abs(data['tilt']) > 90.:
1132-
# G2fil.G2Print ('next line is a hyperbola - search stopped')
1133-
# break
11341026
if debug: print ('HKLD:'+str(H[:4])+'2-theta: %.4f'%(tth))
11351027
elcent,phi,radii = ellipse = GetEllipse(dsp,data)
11361028
data['ellipses'].append(copy.deepcopy(ellipse+('g',)))

‎GSASII/GSASIIimgGUI.py‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -2060,8 +2060,8 @@ def OnFindPixelMask(event):
20602060
LUtth = np.array(Controls['IOtth'])
20612061
dsp0 = wave/(2.0*sind(LUtth[0]/2.0))
20622062
dsp1 = wave/(2.0*sind(LUtth[1]/2.0))
2063-
x0 = G2img.GetDetectorXY2(dsp0,0.0,Controls)[0]
2064-
x1 = G2img.GetDetectorXY2(dsp1,0.0,Controls)[0]
2063+
x0 = G2img.GetDetectorXY(dsp0,0.0,Controls)[0]
2064+
x1 = G2img.GetDetectorXY(dsp1,0.0,Controls)[0]
20652065
if not np.any(x0) or not np.any(x1):
20662066
raise Exception
20672067
nChans = int(1000*(x1-x0)/Controls['pixelSize'][0])//2
@@ -2126,8 +2126,8 @@ def OnAutoFindPixelMask(event):
21262126
LUtth = np.array(Controls['IOtth'])
21272127
dsp0 = wave/(2.0*sind(LUtth[0]/2.0))
21282128
dsp1 = wave/(2.0*sind(LUtth[1]/2.0))
2129-
x0 = G2img.GetDetectorXY2(dsp0,0.0,Controls)[0]
2130-
x1 = G2img.GetDetectorXY2(dsp1,0.0,Controls)[0]
2129+
x0 = G2img.GetDetectorXY(dsp0,0.0,Controls)[0]
2130+
x1 = G2img.GetDetectorXY(dsp1,0.0,Controls)[0]
21312131
if not np.any(x0) or not np.any(x1):
21322132
raise Exception
21332133
nChans = int(1000*(x1-x0)/Controls['pixelSize'][0])//2

‎GSASII/GSASIIplot.py‎

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -4561,10 +4561,10 @@ def ComputeArc(angI,angO,wave,azm0=0,azm1=362):
45614561
if azm1-azm0 > 180: aR[2] //= 2 # for more than 180 degrees, steps can be 2 deg.
45624562
Azm = np.linspace(*aR)
45634563
for azm in Azm:
4564-
XY = G2img.GetDetectorXY2(Dsp(np.squeeze(angI),wave),azm,Data)
4564+
XY = G2img.GetDetectorXY(Dsp(np.squeeze(angI),wave),azm,Data)
45654565
if np.any(XY):
45664566
xy1.append(XY)
4567-
XY = G2img.GetDetectorXY2(Dsp(np.squeeze(angO),wave),azm,Data)
4567+
XY = G2img.GetDetectorXY(Dsp(np.squeeze(angO),wave),azm,Data)
45684568
if np.any(XY):
45694569
xy2.append(XY)
45704570
return np.array(xy1).T,np.array(xy2).T
@@ -4845,7 +4845,7 @@ def OnDragIntBound(event):
48454845
if ellI:
48464846
xyI = []
48474847
for azm in Azm:
4848-
xy = G2img.GetDetectorXY2(dspI,azm,Data)
4848+
xy = G2img.GetDetectorXY(dspI,azm,Data)
48494849
if np.any(xy):
48504850
xyI.append(xy)
48514851
if len(xyI):
@@ -4854,7 +4854,7 @@ def OnDragIntBound(event):
48544854
if ellO:
48554855
xyO = []
48564856
for azm in Azm:
4857-
xy = G2img.GetDetectorXY2(dspO,azm,Data)
4857+
xy = G2img.GetDetectorXY(dspO,azm,Data)
48584858
if np.any(xy):
48594859
xyO.append(xy)
48604860
if len(xyO):
@@ -4874,9 +4874,9 @@ def OnDragIntBound(event):
48744874
xlim = Plot1.get_xlim()
48754875
azm = Data['linescan'][1]-AzmthOff
48764876
dspI = wave/(2.0*sind(IOtth[0]/2.0))
4877-
xyI = G2img.GetDetectorXY2(dspI,azm,Data)
4877+
xyI = G2img.GetDetectorXY(dspI,azm,Data)
48784878
dspO = wave/(2.0*sind(IOtth[1]/2.0))
4879-
xyO = G2img.GetDetectorXY2(dspO,azm,Data)
4879+
xyO = G2img.GetDetectorXY(dspO,azm,Data)
48804880
pick.set_data([[xyI[0],xyO[0]],[xyI[1],xyO[1]]])
48814881
xy = G2img.GetLineScan(G2frame.ImageZ,Data)
48824882
Plot1.cla()
@@ -5461,7 +5461,7 @@ def OnImRelease(event):
54615461
if ellI:
54625462
xyI = []
54635463
for azm in Azm:
5464-
xy = G2img.GetDetectorXY2(dspI,azm,Data)
5464+
xy = G2img.GetDetectorXY(dspI,azm,Data)
54655465
if np.any(xy):
54665466
xyI.append(xy)
54675467
if len(xyI):
@@ -5472,7 +5472,7 @@ def OnImRelease(event):
54725472
xyO = []
54735473
arcxO = []
54745474
for azm in Azm:
5475-
xy = G2img.GetDetectorXY2(dspO,azm,Data)
5475+
xy = G2img.GetDetectorXY(dspO,azm,Data)
54765476
if np.any(xy):
54775477
xyO.append(xy)
54785478
if len(xyO):
@@ -5496,9 +5496,9 @@ def OnImRelease(event):
54965496
IOtth = Data['IOtth']
54975497
wave = Data['wavelength']
54985498
dspI = wave/(2.0*sind(IOtth[0]/2.0))
5499-
xyI = G2img.GetDetectorXY2(dspI,azm,Data)
5499+
xyI = G2img.GetDetectorXY(dspI,azm,Data)
55005500
dspO = wave/(2.0*sind(IOtth[1]/2.0))
5501-
xyO = G2img.GetDetectorXY2(dspO,azm,Data)
5501+
xyO = G2img.GetDetectorXY(dspO,azm,Data)
55025502
Plot.plot([xyI[0],xyO[0]],[xyI[1],xyO[1]],picker=False,label='linescan')
55035503

55045504
if G2frame.PickId and G2frame.GPXtree.GetItemText(G2frame.PickId) in ['Image Controls',]:

‎GSASII/GSASIIscriptable.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -6868,8 +6868,8 @@ def GeneratePixelMask(self,esdMul=3.0,ttmin=0.,ttmax=180.,
68686868
wave = Controls['wavelength']
68696869
dsp0 = wave/(2.0*sind(LUtth[0]/2.0))
68706870
dsp1 = wave/(2.0*sind(LUtth[1]/2.0))
6871-
x0 = G2img.GetDetectorXY2(dsp0,0.0,Controls)[0]
6872-
x1 = G2img.GetDetectorXY2(dsp1,0.0,Controls)[0]
6871+
x0 = G2img.GetDetectorXY(dsp0,0.0,Controls)[0]
6872+
x1 = G2img.GetDetectorXY(dsp1,0.0,Controls)[0]
68736873
if not np.any(x0) or not np.any(x1):
68746874
raise Exception
68756875
numChans = int(1000*(x1-x0)/Controls['pixelSize'][0])//2

0 commit comments

Comments
 (0)