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

magnetic_scattrelease-4.2.2ticket-1009ticket-1249
Last change on this file since cb64d86 was 574adc7, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

convert sascalc to python 2/3 syntax

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