manipulations.py @ c6f95bb

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

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

SasView

source: sasview/DataLoader/manipulations.py @ c6f95bb

Download in other formats: