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

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.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since fa4af76 was fa4af76, checked in by Ricardo Ferraz Leal <ricleal@…>, 7 years ago

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