manipulations.py @ 8c9e65c

magnetic_scattrelease-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1249

Last change on this file since 8c9e65c was e4e9162, checked in by Paul Kienzle <pkienzle@…>, 6 years ago
update code comments
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	with np.errstate(divide='ignore', invalid='ignore'):
931	y = y/y_counts
932	y_err = np.sqrt(y_err)/y_counts
933	# The type of averaging: phi, q2, or q
934	# Calculate x values at the center of the bin
935	if run.lower() == 'phi':
936	step = (self.phi_max - self.phi_min) / self.nbins
937	x = (np.arange(self.nbins) + 0.5) * step + self.phi_min
938	else:
939	# set q to the average of the q values within each bin
940	x = x/y_counts
941
942	### Alternate algorithm
943	## We take the center of ring area, not radius.
944	## This is more accurate than taking the radial center of ring.
945	#step = (self.r_max - self.r_min) / self.nbins
946	#r_inner = self.r_min + step * np.arange(self.nbins)
947	#x = math.sqrt((r_inner2 + (r_inner + step)2) / 2)
948
949	idx = (np.isfinite(y) & np.isfinite(y_err))
950	if x_err is not None:
951	d_x = x_err[idx] / y_counts[idx]
952	else:
953	d_x = None
954	if not idx.any():
955	msg = "Average Error: No points inside sector of ROI to average..."
956	raise ValueError(msg)
957	# elif len(y[idx])!= self.nbins:
958	# print "resulted",self.nbins- len(y[idx]),
959	# "empty bin(s) due to tight binning..."
960	return Data1D(x=x[idx], y=y[idx], dy=y_err[idx], dx=d_x)
961
962
963	class SectorPhi(_Sector):
964	"""
965	Sector average as a function of phi.
966	I(phi) is return and the data is averaged over Q.
967
968	A sector is defined by r_min, r_max, phi_min, phi_max.
969	The number of bin in phi also has to be defined.
970	"""
971
972	def __call__(self, data2D):
973	"""
974	Perform sector average and return I(phi).
975
976	:param data2D: Data2D object
977	:return: Data1D object
978	"""
979	return self._agv(data2D, 'phi')
980
981
982	class SectorQ(_Sector):
983	"""
984	Sector average as a function of Q for both symatric wings.
985	I(Q) is return and the data is averaged over phi.
986
987	A sector is defined by r_min, r_max, phi_min, phi_max.
988	r_min, r_max, phi_min, phi_max >0.
989	The number of bin in Q also has to be defined.
990	"""
991
992	def __call__(self, data2D):
993	"""
994	Perform sector average and return I(Q).
995
996	:param data2D: Data2D object
997
998	:return: Data1D object
999	"""
1000	return self._agv(data2D, 'q2')
1001
1002	################################################################################
1003
1004	class Ringcut(object):
1005	"""
1006	Defines a ring on a 2D data set.
1007	The ring is defined by r_min, r_max, and
1008	the position of the center of the ring.
1009
1010	The data returned is the region inside the ring
1011
1012	Phi_min and phi_max should be defined between 0 and 2*pi
1013	in anti-clockwise starting from the x- axis on the left-hand side
1014	"""
1015
1016	def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0):
1017	# Minimum radius
1018	self.r_min = r_min
1019	# Maximum radius
1020	self.r_max = r_max
1021	# Center of the ring in x
1022	self.center_x = center_x
1023	# Center of the ring in y
1024	self.center_y = center_y
1025
1026	def __call__(self, data2D):
1027	"""
1028	Apply the ring to the data set.
1029	Returns the angular distribution for a given q range
1030
1031	:param data2D: Data2D object
1032
1033	:return: index array in the range
1034	"""
1035	if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
1036	raise RuntimeError("Ring cut only take plottable_2D objects")
1037
1038	# Get data
1039	qx_data = data2D.qx_data
1040	qy_data = data2D.qy_data
1041	q_data = np.sqrt(qx_data * qx_data + qy_data * qy_data)
1042
1043	# check whether or not the data point is inside ROI
1044	out = (self.r_min <= q_data) & (self.r_max >= q_data)
1045	return out
1046
1047	################################################################################
1048
1049	class Boxcut(object):
1050	"""
1051	Find a rectangular 2D region of interest.
1052	"""
1053
1054	def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0):
1055	# Minimum Qx value [A-1]
1056	self.x_min = x_min
1057	# Maximum Qx value [A-1]
1058	self.x_max = x_max
1059	# Minimum Qy value [A-1]
1060	self.y_min = y_min
1061	# Maximum Qy value [A-1]
1062	self.y_max = y_max
1063
1064	def __call__(self, data2D):
1065	"""
1066	Find a rectangular 2D region of interest.
1067
1068	:param data2D: Data2D object
1069	:return: mask, 1d array (len = len(data))
1070	with Trues where the data points are inside ROI, otherwise False
1071	"""
1072	mask = self._find(data2D)
1073
1074	return mask
1075
1076	def _find(self, data2D):
1077	"""
1078	Find a rectangular 2D region of interest.
1079
1080	:param data2D: Data2D object
1081
1082	:return: out, 1d array (length = len(data))
1083	with Trues where the data points are inside ROI, otherwise Falses
1084	"""
1085	if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
1086	raise RuntimeError("Boxcut take only plottable_2D objects")
1087	# Get qx_ and qy_data
1088	qx_data = data2D.qx_data
1089	qy_data = data2D.qy_data
1090
1091	# check whether or not the data point is inside ROI
1092	outx = (self.x_min <= qx_data) & (self.x_max > qx_data)
1093	outy = (self.y_min <= qy_data) & (self.y_max > qy_data)
1094
1095	return outx & outy
1096
1097	################################################################################
1098
1099	class Sectorcut(object):
1100	"""
1101	Defines a sector (major + minor) region on a 2D data set.
1102	The sector is defined by phi_min, phi_max,
1103	where phi_min and phi_max are defined by the right
1104	and left lines wrt central line.
1105
1106	Phi_min and phi_max are given in units of radian
1107	and (phi_max-phi_min) should not be larger than pi
1108	"""
1109
1110	def __init__(self, phi_min=0, phi_max=math.pi):
1111	self.phi_min = phi_min
1112	self.phi_max = phi_max
1113
1114	def __call__(self, data2D):
1115	"""
1116	Find a rectangular 2D region of interest.
1117
1118	:param data2D: Data2D object
1119
1120	:return: mask, 1d array (len = len(data))
1121
1122	with Trues where the data points are inside ROI, otherwise False
1123	"""
1124	mask = self._find(data2D)
1125
1126	return mask
1127
1128	def _find(self, data2D):
1129	"""
1130	Find a rectangular 2D region of interest.
1131
1132	:param data2D: Data2D object
1133
1134	:return: out, 1d array (length = len(data))
1135
1136	with Trues where the data points are inside ROI, otherwise Falses
1137	"""
1138	if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
1139	raise RuntimeError("Sectorcut take only plottable_2D objects")
1140	Pi = math.pi
1141	# Get data
1142	qx_data = data2D.qx_data
1143	qy_data = data2D.qy_data
1144
1145	# get phi from data
1146	phi_data = np.arctan2(qy_data, qx_data)
1147
1148	# Get the min and max into the region: -pi <= phi < Pi
1149	phi_min_major = flip_phi(self.phi_min + Pi) - Pi
1150	phi_max_major = flip_phi(self.phi_max + Pi) - Pi
1151	# check for major sector
1152	if phi_min_major > phi_max_major:
1153	out_major = (phi_min_major <= phi_data) + \
1154	(phi_max_major > phi_data)
1155	else:
1156	out_major = (phi_min_major <= phi_data) & (
1157	phi_max_major > phi_data)
1158
1159	# minor sector
1160	# Get the min and max into the region: -pi <= phi < Pi
1161	phi_min_minor = flip_phi(self.phi_min) - Pi
1162	phi_max_minor = flip_phi(self.phi_max) - Pi
1163
1164	# check for minor sector
1165	if phi_min_minor > phi_max_minor:
1166	out_minor = (phi_min_minor <= phi_data) + \
1167	(phi_max_minor >= phi_data)
1168	else:
1169	out_minor = (phi_min_minor <= phi_data) & \
1170	(phi_max_minor >= phi_data)
1171	out = out_major + out_minor
1172
1173	return out

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

SasView

source: sasview/src/sas/sascalc/dataloader/manipulations.py @ 8c9e65c

Download in other formats: