manipulations.py @ bb0b12c

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload

Last change on this file since bb0b12c was bb0b12c, checked in by Jae Cho <jhjcho@…>, 14 years ago
Working on 2d circular and sector average for improvement on pix size effect
Property mode set to `100644`
File size: 43.5 KB

Line
1	"""
2	Data manipulations for 2D data sets.
3	Using the meta data information, various types of averaging
4	are performed in Q-space
5	"""
6
7	"""
8	This software was developed by the University of Tennessee as part of the
9	Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
10	project funded by the US National Science Foundation.
11
12	See the license text in license.txt
13
14	copyright 2008, University of Tennessee
15	"""
16	#TODO: copy the meta data from the 2D object to the resulting 1D object
17
18	from data_info import plottable_2D, Data1D
19	import math
20	import numpy
21
22	def get_q(dx, dy, det_dist, wavelength):
23	"""
24	@param dx: x-distance from beam center [mm]
25	@param dy: y-distance from beam center [mm]
26	@return: q-value at the given position
27	"""
28	# Distance from beam center in the plane of detector
29	plane_dist = math.sqrt(dxdx + dydy)
30	# Half of the scattering angle
31	theta = 0.5*math.atan(plane_dist/det_dist)
32	return (4.0math.pi/wavelength)math.sin(theta)
33
34	def get_q_compo(dx, dy, det_dist, wavelength,compo=None):
35	#This reduces tiny error at very large q.
36	#Implementation of this func is not started yet.<--ToDo
37	if dy==0:
38	if dx>=0:
39	angle_xy=0
40	else:
41	angle_xy=math.pi
42	else:
43	angle_xy=math.atan(dx/dy)
44
45	if compo=="x":
46	out=get_q(dx, dy, det_dist, wavelength)*cos(angle_xy)
47	elif compo=="y":
48	out=get_q(dx, dy, det_dist, wavelength)*sin(angle_xy)
49	else:
50	out=get_q(dx, dy, det_dist, wavelength)
51	return out
52
53	class _Slab(object):
54	"""
55	Compute average I(Q) for a region of interest
56	"""
57	def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0, bin_width=0.001):
58	# Minimum Qx value [A-1]
59	self.x_min = x_min
60	# Maximum Qx value [A-1]
61	self.x_max = x_max
62	# Minimum Qy value [A-1]
63	self.y_min = y_min
64	# Maximum Qy value [A-1]
65	self.y_max = y_max
66	# Bin width (step size) [A-1]
67	self.bin_width = bin_width
68	# If True, I(\|Q\|) will be return, otherwise, negative q-values are allowed
69	self.fold = False
70
71	def __call__(self, data2D): return NotImplemented
72
73	def _avg(self, data2D, maj):
74	"""
75	Compute average I(Q_maj) for a region of interest.
76	The major axis is defined as the axis of Q_maj.
77	The minor axis is the axis that we average over.
78
79	@param data2D: Data2D object
80	@param maj_min: min value on the major axis
81	@return: Data1D object
82	"""
83	if len(data2D.detector) != 1:
84	raise RuntimeError, "_Slab._avg: invalid number of detectors: %g" % len(data2D.detector)
85
86	pixel_width_x = data2D.detector[0].pixel_size.x
87	pixel_width_y = data2D.detector[0].pixel_size.y
88	det_dist = data2D.detector[0].distance
89	wavelength = data2D.source.wavelength
90	center_x = data2D.detector[0].beam_center.x/pixel_width_x
91	center_y = data2D.detector[0].beam_center.y/pixel_width_y
92
93	# Build array of Q intervals
94	if maj=='x':
95	nbins = int(math.ceil((self.x_max-self.x_min)/self.bin_width))
96	qbins = self.bin_width*numpy.arange(nbins)+self.x_min
97	elif maj=='y':
98	nbins = int(math.ceil((self.y_max-self.y_min)/self.bin_width))
99	qbins = self.bin_width*numpy.arange(nbins)+self.y_min
100	else:
101	raise RuntimeError, "_Slab._avg: unrecognized axis %s" % str(maj)
102
103	x = numpy.zeros(nbins)
104	y = numpy.zeros(nbins)
105	err_y = numpy.zeros(nbins)
106	x_counts = numpy.zeros(nbins)
107	y_counts = numpy.zeros(nbins)
108
109	for i in range(numpy.size(data2D.data,1)):
110	# Min and max x-value for the pixel
111	minx = pixel_width_x*(i - center_x)
112	maxx = pixel_width_x*(i+1.0 - center_x)
113
114	qxmin = get_q(minx, 0.0, det_dist, wavelength)
115	qxmax = get_q(maxx, 0.0, det_dist, wavelength)
116
117	# Get the count fraction in x for that pixel
118	frac_min = get_pixel_fraction_square(self.x_min, qxmin, qxmax)
119	frac_max = get_pixel_fraction_square(self.x_max, qxmin, qxmax)
120	frac_x = frac_max - frac_min
121
122	if frac_x == 0:
123	continue
124
125	if maj=='x':
126	dx = pixel_width_x*(i+0.5 - center_x)
127	q_value = get_q(dx, 0.0, det_dist, wavelength)
128	if self.fold==False and dx<0:
129	q_value = -q_value
130	i_q = int(math.ceil((q_value-self.x_min)/self.bin_width)) - 1
131
132	if i_q<0 or i_q>=nbins:
133	continue
134
135	for j in range(numpy.size(data2D.data,0)):
136	# Min and max y-value for the pixel
137	miny = pixel_width_y*(j - center_y)
138	maxy = pixel_width_y*(j+1.0 - center_y)
139
140	qymin = get_q(0.0, miny, det_dist, wavelength)
141	qymax = get_q(0.0, maxy, det_dist, wavelength)
142
143	# Get the count fraction in x for that pixel
144	frac_min = get_pixel_fraction_square(self.y_min, qymin, qymax)
145	frac_max = get_pixel_fraction_square(self.y_max, qymin, qymax)
146	frac_y = frac_max - frac_min
147
148	frac = frac_x * frac_y
149
150	if frac == 0:
151	continue
152
153	if maj=='y':
154	dy = pixel_width_y*(j+0.5 - center_y)
155	q_value = get_q(0.0, dy, det_dist, wavelength)
156	if self.fold==False and dy<0:
157	q_value = -q_value
158	i_q = int(math.ceil((q_value-self.y_min)/self.bin_width)) - 1
159
160	if i_q<0 or i_q>=nbins:
161	continue
162
163	x[i_q] += q_value
164	y[i_q] += frac * data2D.data[j][i]
165	if data2D.err_data == None or data2D.err_data[j][i]==0.0:
166	err_y[i_q] += frac * frac * math.fabs(data2D.data[j][i])
167	else:
168	err_y[i_q] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
169	x_counts[i_q] += 1
170	y_counts[i_q] += frac
171
172	# Average the sums
173	for i in range(nbins):
174	if x_counts[i]>0 and y_counts[i]>0:
175	err_y[i] = math.sqrt(err_y[i])/y_counts[i]
176	x[i] = x[i]/x_counts[i]
177	y[i] = y[i]/y_counts[i]
178
179	return Data1D(x=x, y=y, dy=err_y)
180
181	class SlabY(_Slab):
182	"""
183	Compute average I(Qy) for a region of interest
184	"""
185	def __call__(self, data2D):
186	"""
187	Compute average I(Qy) for a region of interest
188
189	@param data2D: Data2D object
190	@return: Data1D object
191	"""
192	return self._avg(data2D, 'y')
193
194	class SlabX(_Slab):
195	"""
196	Compute average I(Qx) for a region of interest
197	"""
198	def __call__(self, data2D):
199	"""
200	Compute average I(Qx) for a region of interest
201
202	@param data2D: Data2D object
203	@return: Data1D object
204	"""
205	return self._avg(data2D, 'x')
206
207	class Boxsum(object):
208	"""
209	Perform the sum of counts in a 2D region of interest.
210	"""
211	def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0):
212	# Minimum Qx value [A-1]
213	self.x_min = x_min
214	# Maximum Qx value [A-1]
215	self.x_max = x_max
216	# Minimum Qy value [A-1]
217	self.y_min = y_min
218	# Maximum Qy value [A-1]
219	self.y_max = y_max
220
221	def __call__(self, data2D):
222	"""
223	Perform the sum in the region of interest
224
225	@param data2D: Data2D object
226	@return: number of counts, error on number of counts
227	"""
228	y, err_y, y_counts = self._sum(data2D)
229
230	# Average the sums
231	counts = 0 if y_counts==0 else y
232	error = 0 if y_counts==0 else math.sqrt(err_y)
233
234	return counts, error
235
236	def _sum(self, data2D):
237	"""
238	Perform the sum in the region of interest
239	@param data2D: Data2D object
240	@return: number of counts, error on number of counts, number of entries summed
241	"""
242	if len(data2D.detector) != 1:
243	raise RuntimeError, "Circular averaging: invalid number of detectors: %g" % len(data2D.detector)
244
245	pixel_width_x = data2D.detector[0].pixel_size.x
246	pixel_width_y = data2D.detector[0].pixel_size.y
247	det_dist = data2D.detector[0].distance
248	wavelength = data2D.source.wavelength
249	center_x = data2D.detector[0].beam_center.x/pixel_width_x
250	center_y = data2D.detector[0].beam_center.y/pixel_width_y
251
252	y = 0.0
253	err_y = 0.0
254	y_counts = 0.0
255	sign=1
256	for i in range(numpy.size(data2D.data,1)):
257	# Min and max x-value for the pixel
258	minx = pixel_width_x*(i - center_x)
259	maxx = pixel_width_x*(i+1.0 - center_x)
260	if minx>=0:
261	sign=1
262	else:
263	sign=-1
264	qxmin = sign*get_q(minx, 0.0, det_dist, wavelength)
265	if maxx>=0:
266	sign=1
267	else:
268	sign=-1
269	qxmax = sign*get_q(maxx, 0.0, det_dist, wavelength)
270
271	# Get the count fraction in x for that pixel
272	frac_min = get_pixel_fraction_square(self.x_min, qxmin, qxmax)
273	frac_max = get_pixel_fraction_square(self.x_max, qxmin, qxmax)
274	frac_x = frac_max - frac_min
275
276	for j in range(numpy.size(data2D.data,0)):
277	# Min and max y-value for the pixel
278	miny = pixel_width_y*(j - center_y)
279	maxy = pixel_width_y*(j+1.0 - center_y)
280	if miny>=0:
281	sign=1
282	else:
283	sign=-1
284
285	qymin = sign*get_q(0.0, miny, det_dist, wavelength)
286	if maxy>=0:
287	sign=1
288	else:
289	sign=-1
290
291	qymax = sign*get_q(0.0, maxy, det_dist, wavelength)
292
293	# Get the count fraction in x for that pixel
294	frac_min = get_pixel_fraction_square(self.y_min, qymin, qymax)
295	frac_max = get_pixel_fraction_square(self.y_max, qymin, qymax)
296	frac_y = frac_max - frac_min
297
298	frac = frac_x * frac_y
299
300	y += frac * data2D.data[j][i]
301	if data2D.err_data == None or data2D.err_data[j][i]==0.0:
302	err_y += frac * frac * math.fabs(data2D.data[j][i])
303	else:
304	err_y += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
305	y_counts += frac
306
307	return y, err_y, y_counts
308
309	class Boxavg(Boxsum):
310	"""
311	Perform the average of counts in a 2D region of interest.
312	"""
313	def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0):
314	super(Boxavg, self).__init__(x_min=x_min, x_max=x_max, y_min=y_min, y_max=y_max)
315
316	def __call__(self, data2D):
317	"""
318	Perform the sum in the region of interest
319
320	@param data2D: Data2D object
321	@return: average counts, error on average counts
322	"""
323	y, err_y, y_counts = self._sum(data2D)
324
325	# Average the sums
326	counts = 0 if y_counts==0 else y/y_counts
327	error = 0 if y_counts==0 else math.sqrt(err_y)/y_counts
328
329	return counts, error
330
331	def get_pixel_fraction_square(x, xmin, xmax):
332	"""
333	Return the fraction of the length
334	from xmin to x.
335
336	A B
337	+-----------+---------+
338	xmin x xmax
339
340	@param x: x-value
341	@param xmin: minimum x for the length considered
342	@param xmax: minimum x for the length considered
343	@return: (x-xmin)/(xmax-xmin) when xmin < x < xmax
344
345	"""
346	if x<=xmin:
347	return 0.0
348	if x>xmin and x<xmax:
349	return (x-xmin)/(xmax-xmin)
350	else:
351	return 1.0
352
353
354	class CircularAverage(object):
355	"""
356	Perform circular averaging on 2D data
357
358	The data returned is the distribution of counts
359	as a function of Q
360	"""
361	def __init__(self, r_min=0.0, r_max=0.0, bin_width=0.0005):
362	# Minimum radius included in the average [A-1]
363	self.r_min = r_min
364	# Maximum radius included in the average [A-1]
365	self.r_max = r_max
366	# Bin width (step size) [A-1]
367	self.bin_width = bin_width
368
369	def __call__(self, data2D):
370	"""
371	Perform circular averaging on the data
372
373	@param data2D: Data2D object
374	@return: Data1D object
375	"""
376	if len(data2D.detector) != 1:
377	raise RuntimeError, "Circular averaging: invalid number of detectors: %g" % len(data2D.detector)
378
379	pixel_width_x = data2D.detector[0].pixel_size.x
380	pixel_width_y = data2D.detector[0].pixel_size.y
381	det_dist = data2D.detector[0].distance
382	wavelength = data2D.source.wavelength
383	center_x = data2D.detector[0].beam_center.x/pixel_width_x
384	center_y = data2D.detector[0].beam_center.y/pixel_width_y
385
386	# Find out the maximum Q range
387	xwidth = (numpy.size(data2D.data,1))*pixel_width_x
388	dx_max = xwidth - data2D.detector[0].beam_center.x
389	if xwidth-dx_max>dx_max:
390	dx_max = xwidth-dx_max
391
392	ywidth = (numpy.size(data2D.data,0))*pixel_width_y
393	dy_max = ywidth - data2D.detector[0].beam_center.y
394	if ywidth-dy_max>dy_max:
395	dy_max = ywidth-dy_max
396
397	qmax = get_q(dx_max, dy_max, det_dist, wavelength)
398
399	# Build array of Q intervals
400	nbins = int(math.ceil((qmax-self.r_min)/self.bin_width))
401	qbins = self.bin_width*numpy.arange(nbins)+self.r_min
402
403	x = numpy.zeros(nbins)
404	x_counts = numpy.zeros(nbins)
405	y = numpy.zeros(nbins)
406	err_y = numpy.zeros(nbins)
407	y_counts = numpy.zeros(nbins)
408
409	for i in range(numpy.size(data2D.data,1)):
410
411
412	# Min and max x-value for the pixel
413	minx = pixel_width_x*(i - center_x)
414	maxx = pixel_width_x*(i+1.0 - center_x)
415
416	for j in range(numpy.size(data2D.data,0)):
417
418
419	# Min and max y-value for the pixel
420	miny = pixel_width_y*(j - center_y)
421	maxy = pixel_width_y*(j+1.0 - center_y)
422
423	# Calculate the q-value for each corner
424	# q_[x min or max][y min or max]
425	q_00 = get_q(minx, miny, det_dist, wavelength)
426	q_01 = get_q(minx, maxy, det_dist, wavelength)
427	q_10 = get_q(maxx, miny, det_dist, wavelength)
428	q_11 = get_q(maxx, maxy, det_dist, wavelength)
429
430	#number of sub-bin: default = 1 (ie., 1 bin: no sub-bin)
431	nsubbins = 2
432
433	#if is_intercept(self.r_max, q_00, q_01, q_10, q_11) \
434	# or is_intercept(self.r_min, q_00, q_01, q_10, q_11):
435	# # Number of sub-bins near the boundary of ROI in x-direction
436	# nsubbins = 3
437
438	#Sub divide one pixel into nine sub-pixels only for the pixels where the qmax or qmin line crosses over.
439	for x_subbins in range(nsubbins):
440	for y_subbins in range(nsubbins):
441
442	dx = pixel_width_x*(i+(x_subbins)/nsubbins - center_x)
443	dy = pixel_width_y*(j+(y_subbins)/nsubbins - center_y)
444
445	q_value = get_q(dx, dy, det_dist, wavelength)
446	# Min and max x-value for the subpixel
447	minx = pixel_width_x*(i+(-0.5+x_subbins)/nsubbins - center_x)
448	maxx = pixel_width_x*(i+(0.5 +x_subbins)/nsubbins - center_x)
449	# Min and max y-value for the subpixel
450	miny = pixel_width_y*(j+(-0.5+y_subbins)/nsubbins - center_y)
451	maxy = pixel_width_y*(j+(0.5 +y_subbins)/nsubbins - center_y)
452	# Look for intercept between each side of the pixel
453	# and the constant-q ring for qmax
454	# Calculate the q-value for each corner of sub pixel
455	# q_[x min or max][y min or max]
456	q_00 = get_q(minx, miny, det_dist, wavelength)
457	q_01 = get_q(minx, maxy, det_dist, wavelength)
458	q_10 = get_q(maxx, miny, det_dist, wavelength)
459	q_11 = get_q(maxx, maxy, det_dist, wavelength)
460	frac_max = get_pixel_fraction(self.r_max, q_00, q_01, q_10, q_11)/(nsubbins*nsubbins)
461
462	# Look for intercept between each side of the pixel
463	# and the constant-q ring for qmin
464	frac_min = get_pixel_fraction(self.r_min, q_00, q_01, q_10, q_11)/(nsubbins*nsubbins)
465
466	# We are interested in the region between qmin and qmax
467	# therefore the fraction of the surface of the pixel
468	# that we will use to calculate the number of counts to
469	# include is given by:
470	frac = (frac_max - frac_min) #??? Todo: What if frac_max - fracmin <0?
471
472
473	i_q = int(math.floor((q_value-self.r_min)/self.bin_width)) #- 1
474	if q_value > qmax or q_value < self.r_min:
475	continue
476	#print q_value
477
478	x[i_q] += q_value
479	x_counts[i_q] += 1
480
481	y[i_q] += frac * data2D.data[j][i]
482
483	if data2D.err_data == None or data2D.err_data[j][i]==0.0:
484	err_y[i_q] += frac * frac * math.fabs(data2D.data[j][i])
485	else:
486	err_y[i_q] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
487	y_counts[i_q] += frac
488
489	# Average the sums
490	nzero = 0
491
492	for i in range(nbins):
493	#print x_counts[i],x[i],y[i]
494	if y_counts[i]>0 and x_counts[i]>0:
495	err_y[i] = math.sqrt(err_y[i])/y_counts[i]
496	y[i] = y[i]/y_counts[i]
497	x[i] = (qmax-self.r_min)/nbins(1.0i + 0.5)+ self.r_min#x[i]/x_counts[i]
498
499	else:
500	nzero += 1
501	## Get rid of NULL points
502	tx = numpy.zeros(nbins-nzero)
503	ty = numpy.zeros(nbins-nzero)
504	terr_y = numpy.zeros(nbins-nzero)
505	j=0
506	for i in range(nbins):
507	if err_y[i] != 0 and y_counts[i]>0 and x_counts[i]>0:
508	#Make sure wintin range
509	if x[i] <= self.r_max and x[i] >= self.r_min:
510	tx[j] = x[i]
511	ty[j] = y[i]
512	terr_y[j] = err_y[i]
513	j+=1
514	return Data1D(x=tx, y=ty, dy=terr_y)
515
516
517	class Ring(object):
518	"""
519	Defines a ring on a 2D data set.
520	The ring is defined by r_min, r_max, and
521	the position of the center of the ring.
522
523	The data returned is the distribution of counts
524	around the ring as a function of phi.
525
526	"""
527
528	def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0,nbins=20 ):
529	# Minimum radius
530	self.r_min = r_min
531	# Maximum radius
532	self.r_max = r_max
533	# Center of the ring in x
534	self.center_x = center_x
535	# Center of the ring in y
536	self.center_y = center_y
537	# Number of angular bins
538	self.nbins_phi = nbins
539
540	def __call__(self, data2D):
541	"""
542	Apply the ring to the data set.
543	Returns the angular distribution for a given q range
544
545	@param data2D: Data2D object
546	@return: Data1D object
547	"""
548	if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
549	raise RuntimeError, "Ring averaging only take plottable_2D objects"
550
551	data = data2D.data
552	qmin = self.r_min
553	qmax = self.r_max
554
555	if len(data2D.detector) != 1:
556	raise RuntimeError, "Ring averaging: invalid number of detectors: %g" % len(data2D.detector)
557	pixel_width_x = data2D.detector[0].pixel_size.x
558	pixel_width_y = data2D.detector[0].pixel_size.y
559	det_dist = data2D.detector[0].distance
560	wavelength = data2D.source.wavelength
561	#center_x = self.center_x/pixel_width_x
562	#center_y = self.center_y/pixel_width_y
563	center_x = data2D.detector[0].beam_center.x/pixel_width_x
564	center_y = data2D.detector[0].beam_center.y/pixel_width_y
565
566
567	phi_bins = numpy.zeros(self.nbins_phi)
568	phi_counts = numpy.zeros(self.nbins_phi)
569	phi_values = numpy.zeros(self.nbins_phi)
570	phi_err = numpy.zeros(self.nbins_phi)
571
572	for i in range(numpy.size(data,1)):
573	dx = pixel_width_x*(i+0.5 - center_x)
574
575	# Min and max x-value for the pixel
576	minx = pixel_width_x*(i - center_x)
577	maxx = pixel_width_x*(i+1.0 - center_x)
578
579	for j in range(numpy.size(data,0)):
580	dy = pixel_width_y*(j+0.5 - center_y)
581
582	q_value = get_q(dx, dy, det_dist, wavelength)
583
584	# Min and max y-value for the pixel
585	miny = pixel_width_y*(j - center_y)
586	maxy = pixel_width_y*(j+1.0 - center_y)
587
588	# Calculate the q-value for each corner
589	# q_[x min or max][y min or max]
590	q_00 = get_q(minx, miny, det_dist, wavelength)
591	q_01 = get_q(minx, maxy, det_dist, wavelength)
592	q_10 = get_q(maxx, miny, det_dist, wavelength)
593	q_11 = get_q(maxx, maxy, det_dist, wavelength)
594
595	# Look for intercept between each side of the pixel
596	# and the constant-q ring for qmax
597	frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11)
598
599	# Look for intercept between each side of the pixel
600	# and the constant-q ring for qmin
601	frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11)
602
603	# We are interested in the region between qmin and qmax
604	# therefore the fraction of the surface of the pixel
605	# that we will use to calculate the number of counts to
606	# include is given by:
607
608	frac = frac_max - frac_min
609
610	i_phi = int(math.ceil(self.nbins_phi(math.atan2(dy, dx)+math.pi)/(2.0math.pi))) - 1
611
612	phi_bins[i_phi] += frac * data[j][i]
613
614	if data2D.err_data == None or data2D.err_data[j][i]==0.0:
615	phi_err[i_phi] += frac * frac * math.fabs(data2D.data[j][i])
616	else:
617	phi_err[i_phi] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
618	phi_counts[i_phi] += frac
619
620	for i in range(self.nbins_phi):
621	phi_bins[i] = phi_bins[i] / phi_counts[i]
622	phi_err[i] = math.sqrt(phi_err[i]) / phi_counts[i]
623	phi_values[i] = 2.0math.pi/self.nbins_phi(1.0*i + 0.5)-math.pi # move the pi back to -pi <-->+pi
624
625	return Data1D(x=phi_values, y=phi_bins, dy=phi_err)
626
627	def get_pixel_fraction(qmax, q_00, q_01, q_10, q_11):
628	"""
629	Returns the fraction of the pixel defined by
630	the four corners (q_00, q_01, q_10, q_11) that
631	has q < qmax.
632
633	q_01 q_11
634	y=1 +--------------+
635	\| \|
636	\| \|
637	\| \|
638	y=0 +--------------+
639	q_00 q_10
640
641	x=0 x=1
642
643	"""
644
645	# y side for x = minx
646	x_0 = get_intercept(qmax, q_00, q_01)
647	# y side for x = maxx
648	x_1 = get_intercept(qmax, q_10, q_11)
649
650	# x side for y = miny
651	y_0 = get_intercept(qmax, q_00, q_10)
652	# x side for y = maxy
653	y_1 = get_intercept(qmax, q_01, q_11)
654
655	# surface fraction for a 1x1 pixel
656	frac_max = 0
657
658	if x_0 and x_1:
659	frac_max = (x_0+x_1)/2.0
660
661	elif y_0 and y_1:
662	frac_max = (y_0+y_1)/2.0
663
664	elif x_0 and y_0:
665	if q_00 < q_10:
666	frac_max = x_0*y_0/2.0
667	else:
668	frac_max = 1.0-x_0*y_0/2.0
669
670	elif x_0 and y_1:
671	if q_00 < q_10:
672	frac_max = x_0*y_1/2.0
673	else:
674	frac_max = 1.0-x_0*y_1/2.0
675
676	elif x_1 and y_0:
677	if q_00 > q_10:
678	frac_max = x_1*y_0/2.0
679	else:
680	frac_max = 1.0-x_1*y_0/2.0
681
682	elif x_1 and y_1:
683	if q_00 < q_10:
684	frac_max = 1.0 - (1.0-x_1)*(1.0-y_1)/2.0
685	else:
686	frac_max = (1.0-x_1)*(1.0-y_1)/2.0
687
688	# If we make it here, there is no intercept between
689	# this pixel and the constant-q ring. We only need
690	# to know if we have to include it or exclude it.
691	elif (q_00+q_01+q_10+q_11)/4.0 < qmax:
692	frac_max = 1.0
693
694	return frac_max
695
696	def get_intercept(q, q_0, q_1):
697	"""
698	Returns the fraction of the side at which the
699	q-value intercept the pixel, None otherwise.
700	The values returned is the fraction ON THE SIDE
701	OF THE LOWEST Q.
702
703
704
705	A B
706	+-----------+--------+
707	0 1 <--- pixel size
708
709	Q_0 -------- Q ----- Q_1 <--- equivalent Q range
710
711
712	if Q_1 > Q_0, A is returned
713	if Q_1 < Q_0, B is returned
714
715	if Q is outside the range of [Q_0, Q_1], None is returned
716
717	"""
718	if q_1 > q_0:
719	if (q > q_0 and q <= q_1):
720	return (q-q_0)/(q_1 - q_0)
721	else:
722	if (q > q_1 and q <= q_0):
723	return (q-q_1)/(q_0 - q_1)
724	return None
725
726	class _Sector:
727	"""
728	Defines a sector region on a 2D data set.
729	The sector is defined by r_min, r_max, phi_min, phi_max,
730	and the position of the center of the ring
731	where phi_min and phi_max are defined by the right and left lines wrt central line
732	and phi_max could be less than phi_min.
733
734	Phi is defined between 0 and 2pi
735	"""
736	def __init__(self, r_min, r_max, phi_min, phi_max,nbins=20):
737	self.r_min = r_min
738	self.r_max = r_max
739	self.phi_min = phi_min
740	self.phi_max = phi_max
741	self.nbins = nbins
742
743	def _agv(self, data2D, run='phi'):
744	"""
745	Perform sector averaging.
746
747	@param data2D: Data2D object
748	@param run: define the varying parameter ('phi' or 'q')
749	@return: Data1D object
750	"""
751	if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
752	raise RuntimeError, "Ring averaging only take plottable_2D objects"
753
754	data = data2D.data
755	qmax = self.r_max
756	qmin = self.r_min
757
758	if len(data2D.detector) != 1:
759	raise RuntimeError, "Ring averaging: invalid number of detectors: %g" % len(data2D.detector)
760	pixel_width_x = data2D.detector[0].pixel_size.x
761	pixel_width_y = data2D.detector[0].pixel_size.y
762	det_dist = data2D.detector[0].distance
763	wavelength = data2D.source.wavelength
764	center_x = data2D.detector[0].beam_center.x/pixel_width_x
765	center_y = data2D.detector[0].beam_center.y/pixel_width_y
766
767	y = numpy.zeros(self.nbins)
768	y_counts = numpy.zeros(self.nbins)
769	x = numpy.zeros(self.nbins)
770	x_counts = numpy.zeros(self.nbins)
771	y_err = numpy.zeros(self.nbins)
772
773
774	center_x=center_x
775	center_y=center_y
776
777	# This If finds qmax within ROI defined by sector lines
778	temp=0 #to find qmax within ROI or phimax and phimin
779	temp0=1000000
780	temp1=0
781	for i in range(numpy.size(data,1)):
782	for j in range(numpy.size(data,0)):
783	#number of sub-bin: default = 2 (ie., 1 bin: no sub-bin)
784	nsubbins = 2
785
786	#Sub divide one pixel into nine sub-pixels only for the pixels where the qmax or qmin line crosses over.
787	for x_subbins in range(nsubbins):
788	for y_subbins in range(nsubbins):
789
790	dx = pixel_width_x*(i+(x_subbins)/nsubbins - center_x)
791	dy = pixel_width_y*(j+(y_subbins)/nsubbins - center_y)
792
793	q_value = get_q(dx, dy, det_dist, wavelength)
794
795	# Compute phi and check whether it's within the limits
796	phi_value=math.atan2(dy,dx)+math.pi
797	if self.phi_max>2*math.pi:
798	self.phi_max=self.phi_max-2*math.pi
799	if self.phi_min<0:
800	self.phi_max=self.phi_max+2*math.pi
801
802	#In case of two ROI (symmetric major and minor wings)(for 'q2')
803	if run.lower()=='q2':
804	if ((self.phi_max>=0 and self.phi_max<math.pi)and (self.phi_min>=0 and self.phi_min<math.pi)):
805	temp_max=self.phi_max+math.pi
806	temp_min=self.phi_min+math.pi
807	else:
808	temp_max=self.phi_max
809	temp_min=self.phi_min
810
811	if ((temp_max>=math.pi and temp_max<2math.pi)and (temp_min>=math.pi and temp_min<2math.pi)):
812	if (phi_value<temp_min or phi_value>temp_max):
813	if (phi_value<temp_min-math.pi or phi_value>temp_max-math.pi):
814	continue
815	if (self.phi_max<self.phi_min):
816	tmp_max=self.phi_max+math.pi
817	tmp_min=self.phi_min-math.pi
818	else:
819	tmp_max=self.phi_max
820	tmp_min=self.phi_min
821	if (tmp_min<math.pi and tmp_max>math.pi):
822	if((phi_value>tmp_max and phi_value<tmp_min+math.pi)or (phi_value>tmp_max-math.pi and phi_value<tmp_min)):
823	continue
824	#In case of one ROI (major only)(i.e.,for 'q')and do nothing for 'phi'.
825	elif run.lower()=='q':
826	if (self.phi_max>=self.phi_min):
827	if (phi_value<self.phi_min or phi_value>self.phi_max):
828	continue
829	else:
830	if (phi_value<self.phi_min and phi_value>self.phi_max):
831	continue
832	if q_value<qmin or q_value>qmax:
833	continue
834
835	if run.lower()=='phi':
836	if temp1<phi_value:
837	temp1=phi_value
838	if temp0>phi_value:
839	temp0=phi_value
840
841	elif temp<q_value:
842	temp=q_value
843
844	if run.lower()=='phi':
845	self.phi_max=temp1
846	self.phi_min=temp0
847	else:
848	qmax=temp
849	#Beam center is already corrected, but the calculation below assumed it was not.
850	# Thus Beam center shifted back to uncorrected value. ToDo: cleanup the mess.
851	#center_x=center_x+0.5
852	#center_y=center_y+0.5
853	for i in range(numpy.size(data,1)):
854	for j in range(numpy.size(data,0)):
855	#number of sub-bin: default = 1 (ie., 1 bin: no sub-bin)
856	nsubbins = 2
857
858	#if is_intercept(self.r_max, q_00, q_01, q_10, q_11) \
859	# or is_intercept(self.r_min, q_00, q_01, q_10, q_11):
860	# # Number of sub-bins near the boundary of ROI in x-direction
861	# nsubbins = 3
862
863	#Sub divide one pixel into nine sub-pixels only for the pixels where the qmax or qmin line crosses over.
864	for x_subbins in range(nsubbins):
865	for y_subbins in range(nsubbins):
866
867	dx = pixel_width_x*(i+(x_subbins)/nsubbins - center_x)
868	dy = pixel_width_y*(j+(y_subbins)/nsubbins - center_y)
869
870	q_value = get_q(dx, dy, det_dist, wavelength)
871	# Min and max x-value for the subpixel
872	minx = pixel_width_x*(i+(-0.5+x_subbins)/nsubbins - center_x)
873	maxx = pixel_width_x*(i+(0.5 +x_subbins)/nsubbins - center_x)
874	# Min and max y-value for the subpixel
875	miny = pixel_width_y*(j+(-0.5+y_subbins)/nsubbins - center_y)
876	maxy = pixel_width_y*(j+(0.5 +y_subbins)/nsubbins - center_y)
877
878
879	# Calculate the q-value for each corner
880	# q_[x min or max][y min or max]
881	q_00 = get_q(minx, miny, det_dist, wavelength)
882	q_01 = get_q(minx, maxy, det_dist, wavelength)
883	q_10 = get_q(maxx, miny, det_dist, wavelength)
884	q_11 = get_q(maxx, maxy, det_dist, wavelength)
885
886	# Compute phi and check whether it's within the limits
887	phi_value=math.atan2(dy,dx)+math.pi
888	if self.phi_max>2*math.pi:
889	self.phi_max=self.phi_max-2*math.pi
890	if self.phi_min<0:
891	self.phi_max=self.phi_max+2*math.pi
892
893	# Look for intercept between each side of the pixel
894	# and the constant-q ring for qmax
895	frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11)
896
897	# Look for intercept between each side of the pixel
898	# and the constant-q ring for qmin
899	frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11)
900
901	# We are interested in the region between qmin and qmax
902	# therefore the fraction of the surface of the pixel
903	# that we will use to calculate the number of counts to
904	# include is given by:
905
906	frac = frac_max - frac_min
907
908	#In case of two ROI (symmetric major and minor regions)(for 'q2')
909	if run.lower()=='q2':
910	if ((self.phi_max>=0 and self.phi_max<math.pi)and (self.phi_min>=0 and self.phi_min<math.pi)):
911	temp_max=self.phi_max+math.pi
912	temp_min=self.phi_min+math.pi
913	else:
914	temp_max=self.phi_max
915	temp_min=self.phi_min
916
917	if ((temp_max>=math.pi and temp_max<2math.pi)and (temp_min>=math.pi and temp_min<2math.pi)):
918	if (phi_value<temp_min or phi_value>temp_max):
919	if (phi_value<temp_min-math.pi or phi_value>temp_max-math.pi):
920	continue
921	if (self.phi_max<self.phi_min):
922	tmp_max=self.phi_max+math.pi
923	tmp_min=self.phi_min-math.pi
924	else:
925	tmp_max=self.phi_max
926	tmp_min=self.phi_min
927	if (tmp_min<math.pi and tmp_max>math.pi):
928	if((phi_value>tmp_max and phi_value<tmp_min+math.pi)or (phi_value>tmp_max-math.pi and phi_value<tmp_min)):
929	continue
930	#In case of one ROI (major only)(i.e.,for 'q' and 'phi')
931	else:
932	if (self.phi_max>=self.phi_min):
933	if (phi_value<self.phi_min or phi_value>self.phi_max):
934	continue
935
936	# Check which type of averaging we need
937	if run.lower()=='phi':
938	i_bin = int(math.floor(self.nbins*(phi_value-self.phi_min)/(self.phi_max-self.phi_min)))
939	else:
940	# If we don't need this pixel, skip the rest of the work
941	#TODO: an improvement here would be to compute the average
942	# Q for the pixel from the part that is covered by
943	# the ring defined by q_min/q_max rather than the complete
944	# pixel
945	i_bin = int(math.floor(self.nbins*(q_value-qmin)/(qmax-qmin)))
946
947	else:
948	if (phi_value<self.phi_min and phi_value>self.phi_max):
949	continue
950	#print "qmax=",qmax,qmin
951
952	if q_value<qmin or q_value>qmax:
953	continue
954
955	# Check which type of averaging we need
956	if run.lower()=='phi':
957	i_bin = int(math.ceil(self.nbins*(phi_value-self.phi_min)/(self.phi_max-self.phi_min))) - 1
958	else:
959	# If we don't need this pixel, skip the rest of the work
960	#TODO: an improvement here would be to compute the average
961	# Q for the pixel from the part that is covered by
962	# the ring defined by q_min/q_max rather than the complete
963	# pixel
964	i_bin = int(math.ceil(self.nbins*(q_value-qmin)/(qmax-qmin))) - 1
965
966	try:
967	y[i_bin] += frac * data[j][i]
968	except:
969	import sys
970	print sys.exc_value
971	print i_bin, frac
972
973	if data2D.err_data == None or data2D.err_data[j][i]==0.0:
974	y_err[i_bin] += frac * frac * math.fabs(data2D.data[j][i])
975	else:
976	y_err[i_bin] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
977	y_counts[i_bin] += frac
978
979
980	for i in range(self.nbins):
981	y[i] = y[i] / y_counts[i]
982	y_err[i] = math.sqrt(y_err[i]) / y_counts[i]
983	# Check which type of averaging we need
984	if run.lower()=='phi':
985	#Calculate x[i] and put back the origin of angle back to the right hand side (from the left).
986	x[i] = ((self.phi_max-self.phi_min)/self.nbins(1.0i + 0.5)+self.phi_min-2math.pi)180/math.pi
987	if x[i]<0:
988	x[i]=180+x[i]
989	else:
990	x[i] = (qmax-qmin)/self.nbins(1.0i + 0.5)+qmin
991
992	return Data1D(x=x, y=y, dy=y_err)
993
994	class SectorPhi(_Sector):
995	"""
996	Sector average as a function of phi.
997	I(phi) is return and the data is averaged over Q.
998
999	A sector is defined by r_min, r_max, phi_min, phi_max.
1000	The number of bin in phi also has to be defined.
1001	"""
1002	def __call__(self, data2D):
1003	"""
1004	Perform sector average and return I(phi).
1005
1006	@param data2D: Data2D object
1007	@return: Data1D object
1008	"""
1009	return self._agv(data2D, 'phi')
1010
1011	class SectorQ(_Sector):
1012	"""
1013	Sector average as a function of Q for both symatric wings.
1014	I(Q) is return and the data is averaged over phi.
1015
1016	A sector is defined by r_min, r_max, phi_min, phi_max.
1017	r_min, r_max, phi_min, phi_max >0.
1018	The number of bin in Q also has to be defined.
1019	"""
1020	def __call__(self, data2D):
1021	"""
1022	Perform sector average and return I(Q).
1023
1024	@param data2D: Data2D object
1025	@return: Data1D object
1026	"""
1027	return self._agv(data2D, 'q2')
1028	if __name__ == "__main__":
1029
1030	from loader import Loader
1031
1032
1033	d = Loader().load('test/MAR07232_rest.ASC')
1034	#d = Loader().load('test/MP_New.sans')
1035
1036
1037	r = SectorQ(r_min=.000001, r_max=.01, phi_min=0.0, phi_max=2*math.pi)
1038	o = r(d)
1039
1040	s = Ring(r_min=.000001, r_max=.01)
1041	p = s(d)
1042
1043	for i in range(len(o.x)):
1044	print o.x[i], o.y[i], o.dy[i]
1045
1046
1047

Note: See TracBrowser for help on using the repository browser.

SasView

source: sasview/DataLoader/manipulations.py @ bb0b12c

Download in other formats: