manipulations.py @ 8dcfb2e

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

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

SasView

source: sasview/DataLoader/manipulations.py @ 8dcfb2e

Download in other formats: