manipulations.py @ 053c769

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 053c769 was 342a506, checked in by Jae Cho <jhjcho@…>, 14 years ago
Added x_err and more correction on determining q values, some utest are updates but many more to go
Property mode set to `100644`
File size: 37.5 KB

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

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

SasView

source: sasview/DataLoader/manipulations.py @ 053c769

Download in other formats: