manipulations.py @ afb93df

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.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload

Last change on this file since afb93df was b699768, checked in by Piotr Rozyczko <piotr.rozyczko@…>, 9 years ago
Initial commit of the refactored SasCalc? module.
Property mode set to `100644`
File size: 38.4 KB

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

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

SasView

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

Download in other formats: