manipulations.py @ 3a41d4e

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload

Last change on this file since 3a41d4e was 8f12385, checked in by Jae Cho <jhjcho@…>, 14 years ago
included mask in masked circular average in 2D
Property mode set to `100644`
File size: 41.1 KB

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

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

SasView

source: sasview/DataLoader/manipulations.py @ 3a41d4e

Download in other formats: