manipulations.py @ 41d6187

Last change on this file since 41d6187 was 574adc7, checked in by Paul Kienzle <pkienzle@…>, 7 years ago
convert sascalc to python 2/3 syntax
Property mode set to `100644`
File size: 37.1 KB

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

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

SasView

source: sasview/src/sas/sascalc/dataloader/manipulations.py @ 41d6187

Download in other formats: