manipulations.py @ 55a0dc1

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 55a0dc1 was 12c5b87, checked in by Jae Cho <jhjcho@…>, 14 years ago
fixed a bug in sector average
Property mode set to `100644`
File size: 35.4 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
433	q_data_max = numpy.max(q_data)
434
435	if len(data2D.q_data) == None:
436	msg = "Circular averaging: invalid q_data: %g" % data2D.q_data
437	raise RuntimeError, msg
438
439	# Build array of Q intervals
440	nbins = int(math.ceil((self.r_max - self.r_min) / self.bin_width))
441	qbins = self.bin_width * numpy.arange(nbins) + self.r_min
442
443	x = numpy.zeros(nbins)
444	y = numpy.zeros(nbins)
445	err_y = numpy.zeros(nbins)
446	y_counts = numpy.zeros(nbins)
447
448	for npt in range(len(data)):
449	frac = 0
450
451	# q-value at the pixel (j,i)
452	q_value = q_data[npt]
453	data_n = data[npt]
454
455	## No need to calculate the frac when all data are within range
456	if self.r_min >= self.r_max:
457	raise ValueError, "Limit Error: min > max"
458
459	if self.r_min <= q_value and q_value <= self.r_max:
460	frac = 1
461	if frac == 0:
462	continue
463	i_q = int(math.floor((q_value - self.r_min) / self.bin_width))
464
465	# Take care of the edge case at phi = 2pi.
466	if i_q == nbins:
467	i_q = nbins -1
468	y[i_q] += frac * data_n
469
470	if err_data == None or err_data[npt] == 0.0:
471	if data_n < 0:
472	data_n = -data_n
473	err_y[i_q] += frac * frac * data_n
474	else:
475	err_y[i_q] += frac * frac * err_data[npt] * err_data[npt]
476	y_counts[i_q] += frac
477
478	## x should be the center value of each bins
479	x = qbins + self.bin_width / 2
480
481	# Average the sums
482	for n in range(nbins):
483	if err_y[n] < 0: err_y[n] = -err_y[n]
484	err_y[n] = math.sqrt(err_y[n])
485
486	err_y = err_y / y_counts
487	y = y / y_counts
488	idx = (numpy.isfinite(y)) & (numpy.isfinite(x))
489
490	if not idx.any():
491	msg = "Average Error: No points inside ROI to average..."
492	raise ValueError, msg
493	#elif len(y[idx])!= nbins:
494	# print "resulted",nbins- len(y[idx])
495	#,"empty bin(s) due to tight binning..."
496	return Data1D(x=x[idx], y=y[idx], dy=err_y[idx])
497
498
499	class Ring(object):
500	"""
501	Defines a ring on a 2D data set.
502	The ring is defined by r_min, r_max, and
503	the position of the center of the ring.
504
505	The data returned is the distribution of counts
506	around the ring as a function of phi.
507
508	Phi_min and phi_max should be defined between 0 and 2*pi
509	in anti-clockwise starting from the x- axis on the left-hand side
510	"""
511	#Todo: remove center.
512	def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0, nbins=20):
513	# Minimum radius
514	self.r_min = r_min
515	# Maximum radius
516	self.r_max = r_max
517	# Center of the ring in x
518	self.center_x = center_x
519	# Center of the ring in y
520	self.center_y = center_y
521	# Number of angular bins
522	self.nbins_phi = nbins
523
524	def __call__(self, data2D):
525	"""
526	Apply the ring to the data set.
527	Returns the angular distribution for a given q range
528
529	:param data2D: Data2D object
530
531	:return: Data1D object
532	"""
533	if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
534	raise RuntimeError, "Ring averaging only take plottable_2D objects"
535
536	Pi = math.pi
537
538	# Get data
539	data = data2D.data[numpy.isfinite(data2D.data)]
540	q_data = data2D.q_data[numpy.isfinite(data2D.data)]
541	err_data = data2D.err_data[numpy.isfinite(data2D.data)]
542	qx_data = data2D.qx_data[numpy.isfinite(data2D.data)]
543	qy_data = data2D.qy_data[numpy.isfinite(data2D.data)]
544
545	q_data_max = numpy.max(q_data)
546
547	# Set space for 1d outputs
548	phi_bins = numpy.zeros(self.nbins_phi)
549	phi_counts = numpy.zeros(self.nbins_phi)
550	phi_values = numpy.zeros(self.nbins_phi)
551	phi_err = numpy.zeros(self.nbins_phi)
552
553	for npt in range(len(data)):
554	frac = 0
555	# q-value at the point (npt)
556	q_value = q_data[npt]
557	data_n = data[npt]
558
559	# phi-value at the point (npt)
560	phi_value = math.atan2(qy_data[npt], qx_data[npt]) + Pi
561
562	if self.r_min <= q_value and q_value <= self.r_max:
563	frac = 1
564	if frac == 0:
565	continue
566	# binning
567	i_phi = int(math.floor((self.nbins_phi) * phi_value / (2 * Pi)))
568
569	# Take care of the edge case at phi = 2pi.
570	if i_phi == self.nbins_phi:
571	i_phi = self.nbins_phi - 1
572	phi_bins[i_phi] += frac * data[npt]
573
574	if err_data == None or err_data[npt] == 0.0:
575	if data_n < 0:
576	data_n = -data_n
577	phi_err[i_phi] += frac * frac * math.fabs(data_n)
578	else:
579	phi_err[i_phi] += frac * frac * err_data[npt] * err_data[npt]
580	phi_counts[i_phi] += frac
581
582	for i in range(self.nbins_phi):
583	phi_bins[i] = phi_bins[i] / phi_counts[i]
584	phi_err[i] = math.sqrt(phi_err[i]) / phi_counts[i]
585	phi_values[i] = 2.0 * math.pi / self.nbins_phi * (1.0 * i + 0.5)
586
587	idx = (numpy.isfinite(phi_bins))
588
589	if not idx.any():
590	msg = "Average Error: No points inside ROI to average..."
591	raise ValueError, msg
592	#elif len(phi_bins[idx])!= self.nbins_phi:
593	# print "resulted",self.nbins_phi- len(phi_bins[idx])
594	#,"empty bin(s) due to tight binning..."
595	return Data1D(x=phi_values[idx], y=phi_bins[idx], dy=phi_err[idx])
596
597	def get_pixel_fraction(qmax, q_00, q_01, q_10, q_11):
598	"""
599	Returns the fraction of the pixel defined by
600	the four corners (q_00, q_01, q_10, q_11) that
601	has q < qmax.::
602
603	q_01 q_11
604	y=1 +--------------+
605	\| \|
606	\| \|
607	\| \|
608	y=0 +--------------+
609	q_00 q_10
610
611	x=0 x=1
612
613	"""
614	# y side for x = minx
615	x_0 = get_intercept(qmax, q_00, q_01)
616	# y side for x = maxx
617	x_1 = get_intercept(qmax, q_10, q_11)
618
619	# x side for y = miny
620	y_0 = get_intercept(qmax, q_00, q_10)
621	# x side for y = maxy
622	y_1 = get_intercept(qmax, q_01, q_11)
623
624	# surface fraction for a 1x1 pixel
625	frac_max = 0
626
627	if x_0 and x_1:
628	frac_max = (x_0 + x_1) / 2.0
629	elif y_0 and y_1:
630	frac_max = (y_0 + y_1) / 2.0
631	elif x_0 and y_0:
632	if q_00 < q_10:
633	frac_max = x_0 * y_0 / 2.0
634	else:
635	frac_max = 1.0 - x_0 * y_0 / 2.0
636	elif x_0 and y_1:
637	if q_00 < q_10:
638	frac_max = x_0 * y_1 / 2.0
639	else:
640	frac_max = 1.0 - x_0 * y_1 / 2.0
641	elif x_1 and y_0:
642	if q_00 > q_10:
643	frac_max = x_1 * y_0 / 2.0
644	else:
645	frac_max = 1.0 - x_1 * y_0 / 2.0
646	elif x_1 and y_1:
647	if q_00 < q_10:
648	frac_max = 1.0 - (1.0 - x_1) * (1.0 - y_1) / 2.0
649	else:
650	frac_max = (1.0 - x_1) * (1.0 - y_1) / 2.0
651
652	# If we make it here, there is no intercept between
653	# this pixel and the constant-q ring. We only need
654	# to know if we have to include it or exclude it.
655	elif (q_00 + q_01 + q_10 + q_11)/4.0 < qmax:
656	frac_max = 1.0
657
658	return frac_max
659
660	def get_intercept(q, q_0, q_1):
661	"""
662	Returns the fraction of the side at which the
663	q-value intercept the pixel, None otherwise.
664	The values returned is the fraction ON THE SIDE
665	OF THE LOWEST Q. ::
666
667
668	A B
669	+-----------+--------+ <--- pixel size
670	0 1
671	Q_0 -------- Q ----- Q_1 <--- equivalent Q range
672	if Q_1 > Q_0, A is returned
673	if Q_1 < Q_0, B is returned
674	if Q is outside the range of [Q_0, Q_1], None is returned
675
676	"""
677	if q_1 > q_0:
678	if (q > q_0 and q <= q_1):
679	return (q - q_0)/(q_1 - q_0)
680	else:
681	if (q > q_1 and q <= q_0):
682	return (q - q_1)/(q_0 - q_1)
683	return None
684
685	class _Sector:
686	"""
687	Defines a sector region on a 2D data set.
688	The sector is defined by r_min, r_max, phi_min, phi_max,
689	and the position of the center of the ring
690	where phi_min and phi_max are defined by the right
691	and left lines wrt central line
692	and phi_max could be less than phi_min.
693
694	Phi is defined between 0 and 2*pi in anti-clockwise
695	starting from the x- axis on the left-hand side
696	"""
697	def __init__(self, r_min, r_max, phi_min=0, phi_max=2*math.pi, nbins=20):
698	self.r_min = r_min
699	self.r_max = r_max
700	self.phi_min = phi_min
701	self.phi_max = phi_max
702	self.nbins = nbins
703
704
705	def _agv(self, data2D, run='phi'):
706	"""
707	Perform sector averaging.
708
709	:param data2D: Data2D object
710	:param run: define the varying parameter ('phi' , 'q' , or 'q2')
711
712	:return: Data1D object
713	"""
714	if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
715	raise RuntimeError, "Ring averaging only take plottable_2D objects"
716	Pi = math.pi
717
718	# Get the all data & info
719	data = data2D.data[numpy.isfinite(data2D.data)]
720	q_data = data2D.q_data[numpy.isfinite(data2D.data)]
721	err_data = data2D.err_data[numpy.isfinite(data2D.data)]
722	qx_data = data2D.qx_data[numpy.isfinite(data2D.data)]
723	qy_data = data2D.qy_data[numpy.isfinite(data2D.data)]
724
725	#set space for 1d outputs
726	x = numpy.zeros(self.nbins)
727	y = numpy.zeros(self.nbins)
728	y_err = numpy.zeros(self.nbins)
729	y_counts = numpy.zeros(self.nbins)
730
731	# Get the min and max into the region: 0 <= phi < 2Pi
732	phi_min = flip_phi(self.phi_min)
733	phi_max = flip_phi(self.phi_max)
734
735	q_data_max = numpy.max(q_data)
736
737	for n in range(len(data)):
738	frac = 0
739
740	# q-value at the pixel (j,i)
741	q_value = q_data[n]
742	data_n = data[n]
743
744	# Is pixel within range?
745	is_in = False
746
747	# phi-value of the pixel (j,i)
748	phi_value = math.atan2(qy_data[n], qx_data[n]) + Pi
749
750	## No need to calculate the frac when all data are within range
751	if self.r_min <= q_value and q_value <= self.r_max:
752	frac = 1
753	if frac == 0:
754	continue
755	#In case of two ROIs (symmetric major and minor regions)(for 'q2')
756	if run.lower()=='q2':
757	## For minor sector wing
758	# Calculate the minor wing phis
759	phi_min_minor = flip_phi(phi_min - Pi)
760	phi_max_minor = flip_phi(phi_max - Pi)
761	# Check if phis of the minor ring is within 0 to 2pi
762	if phi_min_minor > phi_max_minor:
763	is_in = (phi_value > phi_min_minor or \
764	phi_value < phi_max_minor)
765	else:
766	is_in = (phi_value > phi_min_minor and \
767	phi_value < phi_max_minor)
768
769	#For all cases(i.e.,for 'q', 'q2', and 'phi')
770	#Find pixels within ROI
771	if phi_min > phi_max:
772	is_in = is_in or (phi_value > phi_min or \
773	phi_value < phi_max)
774	else:
775	is_in = is_in or (phi_value >= phi_min and \
776	phi_value < phi_max)
777
778	if not is_in:
779	frac = 0
780	if frac == 0:
781	continue
782	# Check which type of averaging we need
783	if run.lower() == 'phi':
784	temp_x = (self.nbins) * (phi_value - self.phi_min)
785	temp_y = (self.phi_max - self.phi_min)
786	i_bin = int(math.floor(temp_x / temp_y))
787	else:
788	temp_x = (self.nbins) * (q_value - self.r_min)
789	temp_y = (self.r_max - self.r_min)
790	i_bin = int(math.floor(temp_x / temp_y))
791
792	# Take care of the edge case at phi = 2pi.
793	if i_bin == self.nbins:
794	i_bin = self.nbins - 1
795
796	## Get the total y
797	y[i_bin] += frac * data_n
798
799	if err_data == None or err_data[n] == 0.0:
800	if data_n < 0:
801	data_n = -data_n
802	y_err[i_bin] += frac * frac * data_n
803	else:
804	y_err[i_bin] += frac * frac * err_data[n] * err_data[n]
805	y_counts[i_bin] += frac
806
807	# Organize the results
808	for i in range(self.nbins):
809	y[i] = y[i] / y_counts[i]
810	y_err[i] = math.sqrt(y_err[i]) / y_counts[i]
811
812	# The type of averaging: phi,q2, or q
813	# Calculate x[i]should be at the center of the bin
814	if run.lower()=='phi':
815	x[i] = (self.phi_max - self.phi_min) / self.nbins * \
816	(1.0 * i + 0.5) + self.phi_min
817	else:
818	x[i] = (self.r_max - self.r_min) / self.nbins * \
819	(1.0 * i + 0.5) + self.r_min
820
821	idx = (numpy.isfinite(y) & numpy.isfinite(y_err))
822
823	if not idx.any():
824	msg = "Average Error: No points inside sector of ROI to average..."
825	raise ValueError, msg
826	#elif len(y[idx])!= self.nbins:
827	# print "resulted",self.nbins- len(y[idx]),
828	#"empty bin(s) due to tight binning..."
829	return Data1D(x=x[idx], y=y[idx], dy=y_err[idx])
830
831	class SectorPhi(_Sector):
832	"""
833	Sector average as a function of phi.
834	I(phi) is return and the data is averaged over Q.
835
836	A sector is defined by r_min, r_max, phi_min, phi_max.
837	The number of bin in phi also has to be defined.
838	"""
839	def __call__(self, data2D):
840	"""
841	Perform sector average and return I(phi).
842
843	:param data2D: Data2D object
844	:return: Data1D object
845	"""
846
847	return self._agv(data2D, 'phi')
848
849	class SectorQ(_Sector):
850	"""
851	Sector average as a function of Q for both symatric wings.
852	I(Q) is return and the data is averaged over phi.
853
854	A sector is defined by r_min, r_max, phi_min, phi_max.
855	r_min, r_max, phi_min, phi_max >0.
856	The number of bin in Q also has to be defined.
857	"""
858	def __call__(self, data2D):
859	"""
860	Perform sector average and return I(Q).
861
862	:param data2D: Data2D object
863
864	:return: Data1D object
865	"""
866	return self._agv(data2D, 'q2')
867
868	class Ringcut(object):
869	"""
870	Defines a ring on a 2D data set.
871	The ring is defined by r_min, r_max, and
872	the position of the center of the ring.
873
874	The data returned is the region inside the ring
875
876	Phi_min and phi_max should be defined between 0 and 2*pi
877	in anti-clockwise starting from the x- axis on the left-hand side
878	"""
879	def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0 ):
880	# Minimum radius
881	self.r_min = r_min
882	# Maximum radius
883	self.r_max = r_max
884	# Center of the ring in x
885	self.center_x = center_x
886	# Center of the ring in y
887	self.center_y = center_y
888
889
890	def __call__(self, data2D):
891	"""
892	Apply the ring to the data set.
893	Returns the angular distribution for a given q range
894
895	:param data2D: Data2D object
896
897	:return: index array in the range
898	"""
899	if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
900	raise RuntimeError, "Ring cut only take plottable_2D objects"
901
902	# Get data
903	qx_data = data2D.qx_data
904	qy_data = data2D.qy_data
905	mask = data2D.mask
906	q_data = numpy.sqrt(qx_data * qx_data + qy_data * qy_data)
907	#q_data_max = numpy.max(q_data)
908
909	# check whether or not the data point is inside ROI
910	out = (self.r_min <= q_data) & (self.r_max >= q_data)
911
912	return (out)
913
914
915	class Boxcut(object):
916	"""
917	Find a rectangular 2D region of interest.
918	"""
919	def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0):
920	# Minimum Qx value [A-1]
921	self.x_min = x_min
922	# Maximum Qx value [A-1]
923	self.x_max = x_max
924	# Minimum Qy value [A-1]
925	self.y_min = y_min
926	# Maximum Qy value [A-1]
927	self.y_max = y_max
928
929	def __call__(self, data2D):
930	"""
931	Find a rectangular 2D region of interest.
932
933	:param data2D: Data2D object
934	:return: mask, 1d array (len = len(data))
935	with Trues where the data points are inside ROI, otherwise False
936	"""
937	mask = self._find(data2D)
938
939	return mask
940
941	def _find(self, data2D):
942	"""
943	Find a rectangular 2D region of interest.
944
945	:param data2D: Data2D object
946
947	:return: out, 1d array (length = len(data))
948	with Trues where the data points are inside ROI, otherwise Falses
949	"""
950	if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
951	raise RuntimeError, "Boxcut take only plottable_2D objects"
952	# Get qx_ and qy_data
953	qx_data = data2D.qx_data
954	qy_data = data2D.qy_data
955	mask = data2D.mask
956
957	# check whether or not the data point is inside ROI
958	outx = (self.x_min <= qx_data) & (self.x_max > qx_data)
959	outy = (self.y_min <= qy_data) & (self.y_max > qy_data)
960
961	return (outx & outy)
962
963	class Sectorcut(object):
964	"""
965	Defines a sector (major + minor) region on a 2D data set.
966	The sector is defined by phi_min, phi_max,
967	where phi_min and phi_max are defined by the right
968	and left lines wrt central line.
969
970	Phi_min and phi_max are given in units of radian
971	and (phi_max-phi_min) should not be larger than pi
972	"""
973	def __init__(self, phi_min=0, phi_max=math.pi):
974	self.phi_min = phi_min
975	self.phi_max = phi_max
976
977	def __call__(self, data2D):
978	"""
979	Find a rectangular 2D region of interest.
980
981	:param data2D: Data2D object
982
983	:return: mask, 1d array (len = len(data))
984
985	with Trues where the data points are inside ROI, otherwise False
986	"""
987	mask = self._find(data2D)
988
989	return mask
990
991	def _find(self, data2D):
992	"""
993	Find a rectangular 2D region of interest.
994
995	:param data2D: Data2D object
996
997	:return: out, 1d array (length = len(data))
998
999	with Trues where the data points are inside ROI, otherwise Falses
1000	"""
1001	if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
1002	raise RuntimeError, "Sectorcut take only plottable_2D objects"
1003	Pi = math.pi
1004	# Get data
1005	qx_data = data2D.qx_data
1006	qy_data = data2D.qy_data
1007	phi_data = numpy.zeros(len(qx_data))
1008
1009	# get phi from data
1010	phi_data = numpy.arctan2(qy_data, qx_data)
1011
1012	# Get the min and max into the region: -pi <= phi < Pi
1013	phi_min_major = flip_phi(self.phi_min + Pi) - Pi
1014	phi_max_major = flip_phi(self.phi_max + Pi) - Pi
1015	# check for major sector
1016	if phi_min_major > phi_max_major:
1017	out_major = (phi_min_major <= phi_data) + (phi_max_major > phi_data)
1018	else:
1019	out_major = (phi_min_major <= phi_data) & (phi_max_major > phi_data)
1020
1021	# minor sector
1022	# Get the min and max into the region: -pi <= phi < Pi
1023	phi_min_minor = flip_phi(self.phi_min) - Pi
1024	phi_max_minor = flip_phi(self.phi_max) - Pi
1025
1026	# check for minor sector
1027	if phi_min_minor > phi_max_minor:
1028	out_minor = (phi_min_minor <= phi_data) + \
1029	(phi_max_minor >= phi_data)
1030	else:
1031	out_minor = (phi_min_minor <= phi_data) & \
1032	(phi_max_minor >= phi_data)
1033	out = out_major + out_minor
1034
1035	return out
1036
1037	if __name__ == "__main__":
1038
1039	from loader import Loader
1040
1041
1042	d = Loader().load('test/MAR07232_rest.ASC')
1043	#d = Loader().load('test/MP_New.sans')
1044
1045
1046	r = SectorQ(r_min=.000001, r_max=.01, phi_min=0.0, phi_max=2*math.pi)
1047	o = r(d)
1048
1049	s = Ring(r_min=.000001, r_max=.01)
1050	p = s(d)
1051
1052	for i in range(len(o.x)):
1053	print o.x[i], o.y[i], o.dy[i]
1054
1055
1056

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

SasView

source: sasview/DataLoader/manipulations.py @ 55a0dc1

Download in other formats: