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

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 6c497b7 was 6c497b7, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

Merge branch 'master' into log_binning

  • Property mode set to 100644
File size: 36.9 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    x_max = data2D.dqx_data[data2D.q_data[z_max]]
202    x_min = data2D.dqx_data[data2D.q_data[z_min]]
203    y_max = data2D.dqy_data[data2D.q_data[z_max]]
204    y_min = data2D.dqy_data[data2D.q_data[z_min]]
205    # Find qdx at q = 0
206    dq_overlap_x = (x_min * z_max - x_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 = (y_min * z_max - y_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 = y_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################################################################################
790
791class _Sector(object):
792    """
793    Defines a sector region on a 2D data set.
794    The sector is defined by r_min, r_max, phi_min, phi_max,
795    and the position of the center of the ring
796    where phi_min and phi_max are defined by the right
797    and left lines wrt central line
798    and phi_max could be less than phi_min.
799
800    Phi is defined between 0 and 2*pi in anti-clockwise
801    starting from the x- axis on the left-hand side
802    """
803
804    def __init__(self, r_min, r_max, phi_min=0, phi_max=2 * math.pi, nbins=20, base = None):
805        self.r_min = r_min
806        self.r_max = r_max
807        self.phi_min = phi_min
808        self.phi_max = phi_max
809        self.nbins = nbins
810        self.base = base
811
812    def _agv(self, data2D, run='phi'):
813        """
814        Perform sector averaging.
815
816        :param data2D: Data2D object
817        :param run:  define the varying parameter ('phi' , 'q' , or 'q2')
818
819        :return: Data1D object
820        """
821        if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
822            raise RuntimeError("Ring averaging only take plottable_2D objects")
823
824        # Get the all data & info
825        data = data2D.data[np.isfinite(data2D.data)]
826        q_data = data2D.q_data[np.isfinite(data2D.data)]
827        err_data = data2D.err_data[np.isfinite(data2D.data)]
828        qx_data = data2D.qx_data[np.isfinite(data2D.data)]
829        qy_data = data2D.qy_data[np.isfinite(data2D.data)]
830
831        dq_data = None
832        if data2D.dqx_data is not None and data2D.dqy_data is not None:
833            dq_data = get_dq_data(data2D)
834
835        # set space for 1d outputs
836        x = np.zeros(self.nbins)
837        y = np.zeros(self.nbins)
838        y_err = np.zeros(self.nbins)
839        x_err = np.zeros(self.nbins)
840        y_counts = np.zeros(self.nbins) # Cycle counts (for the mean)
841
842        # Get the min and max into the region: 0 <= phi < 2Pi
843        phi_min = flip_phi(self.phi_min)
844        phi_max = flip_phi(self.phi_max)
845
846        #  binning object
847        if run.lower() == 'phi':
848            binning = Binning(self.phi_min, self.phi_max, self.nbins, self.base)
849        else:
850            binning = Binning(self.r_min, self.r_max, self.nbins, self.base)
851
852        for n in range(len(data)):
853
854            # q-value at the pixel (j,i)
855            q_value = q_data[n]
856            data_n = data[n]
857
858            # Is pixel within range?
859            is_in = False
860
861            # phi-value of the pixel (j,i)
862            phi_value = math.atan2(qy_data[n], qx_data[n]) + math.pi
863
864            # No need to calculate: data outside of the radius
865            if self.r_min > q_value or q_value > self.r_max:
866                continue
867
868            # In case of two ROIs (symmetric major and minor regions)(for 'q2')
869            if run.lower() == 'q2':
870                # For minor sector wing
871                # Calculate the minor wing phis
872                phi_min_minor = flip_phi(phi_min - math.pi)
873                phi_max_minor = flip_phi(phi_max - math.pi)
874                # Check if phis of the minor ring is within 0 to 2pi
875                if phi_min_minor > phi_max_minor:
876                    is_in = (phi_value > phi_min_minor or
877                             phi_value < phi_max_minor)
878                else:
879                    is_in = (phi_value > phi_min_minor and
880                             phi_value < phi_max_minor)
881
882            # For all cases(i.e.,for 'q', 'q2', and 'phi')
883            # Find pixels within ROI
884            if phi_min > phi_max:
885                is_in = is_in or (phi_value > phi_min or
886                                  phi_value < phi_max)
887            else:
888                is_in = is_in or (phi_value >= phi_min and
889                                  phi_value < phi_max)
890
891            # data oustide of the phi range
892            if not is_in:
893                continue
894
895            # Get the binning index
896            if run.lower() == 'phi':
897                i_bin = binning.get_bin_index(phi_value)
898            else:
899                i_bin = binning.get_bin_index(q_value)
900
901            # Take care of the edge case at phi = 2pi.
902            if i_bin == self.nbins:
903                i_bin = self.nbins - 1
904
905            # Get the total y
906            y[i_bin] += data_n
907            x[i_bin] += q_value
908            if err_data[n] is None or err_data[n] == 0.0:
909                if data_n < 0:
910                    data_n = -data_n
911                y_err[i_bin] += data_n
912            else:
913                y_err[i_bin] += err_data[n]**2
914
915            if dq_data is not None:
916                # To be consistent with dq calculation in 1d reduction,
917                # we need just the averages (not quadratures) because
918                # it should not depend on the number of the q points
919                # in the qr bins.
920                x_err[i_bin] += dq_data[n]
921            else:
922                x_err = None
923            y_counts[i_bin] += 1
924
925        # Organize the results
926        for i in range(self.nbins):
927            y[i] = y[i] / y_counts[i]
928            y_err[i] = math.sqrt(y_err[i]) / y_counts[i]
929
930            # The type of averaging: phi,q2, or q
931            # Calculate x[i]should be at the center of the bin
932            if run.lower() == 'phi':
933                x[i] = (self.phi_max - self.phi_min) / self.nbins * \
934                    (1.0 * i + 0.5) + self.phi_min
935            else:
936                # We take the center of ring area, not radius.
937                # This is more accurate than taking the radial center of ring.
938                #delta_r = (self.r_max - self.r_min) / self.nbins
939                #r_inner = self.r_min + delta_r * i
940                #r_outer = r_inner + delta_r
941                #x[i] = math.sqrt((r_inner * r_inner + r_outer * r_outer) / 2)
942                x[i] = x[i] / y_counts[i]
943        y_err[y_err == 0] = np.average(y_err)
944        idx = (np.isfinite(y) & np.isfinite(y_err))
945        if x_err is not None:
946            d_x = x_err[idx] / y_counts[idx]
947        else:
948            d_x = None
949        if not idx.any():
950            msg = "Average Error: No points inside sector of ROI to average..."
951            raise ValueError(msg)
952        # elif len(y[idx])!= self.nbins:
953        #    print "resulted",self.nbins- len(y[idx]),
954        #"empty bin(s) due to tight binning..."
955        return Data1D(x=x[idx], y=y[idx], dy=y_err[idx], dx=d_x)
956
957
958class SectorPhi(_Sector):
959    """
960    Sector average as a function of phi.
961    I(phi) is return and the data is averaged over Q.
962
963    A sector is defined by r_min, r_max, phi_min, phi_max.
964    The number of bin in phi also has to be defined.
965    """
966
967    def __call__(self, data2D):
968        """
969        Perform sector average and return I(phi).
970
971        :param data2D: Data2D object
972        :return: Data1D object
973        """
974        return self._agv(data2D, 'phi')
975
976
977class SectorQ(_Sector):
978    """
979    Sector average as a function of Q for both symatric wings.
980    I(Q) is return and the data is averaged over phi.
981
982    A sector is defined by r_min, r_max, phi_min, phi_max.
983    r_min, r_max, phi_min, phi_max >0.
984    The number of bin in Q also has to be defined.
985    """
986
987    def __call__(self, data2D):
988        """
989        Perform sector average and return I(Q).
990
991        :param data2D: Data2D object
992
993        :return: Data1D object
994        """
995        return self._agv(data2D, 'q2')
996
997################################################################################
998
999class Ringcut(object):
1000    """
1001    Defines a ring on a 2D data set.
1002    The ring is defined by r_min, r_max, and
1003    the position of the center of the ring.
1004
1005    The data returned is the region inside the ring
1006
1007    Phi_min and phi_max should be defined between 0 and 2*pi
1008    in anti-clockwise starting from the x- axis on the left-hand side
1009    """
1010
1011    def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0):
1012        # Minimum radius
1013        self.r_min = r_min
1014        # Maximum radius
1015        self.r_max = r_max
1016        # Center of the ring in x
1017        self.center_x = center_x
1018        # Center of the ring in y
1019        self.center_y = center_y
1020
1021    def __call__(self, data2D):
1022        """
1023        Apply the ring to the data set.
1024        Returns the angular distribution for a given q range
1025
1026        :param data2D: Data2D object
1027
1028        :return: index array in the range
1029        """
1030        if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
1031            raise RuntimeError("Ring cut only take plottable_2D objects")
1032
1033        # Get data
1034        qx_data = data2D.qx_data
1035        qy_data = data2D.qy_data
1036        q_data = np.sqrt(qx_data * qx_data + qy_data * qy_data)
1037
1038        # check whether or not the data point is inside ROI
1039        out = (self.r_min <= q_data) & (self.r_max >= q_data)
1040        return out
1041
1042################################################################################
1043
1044class Boxcut(object):
1045    """
1046    Find a rectangular 2D region of interest.
1047    """
1048
1049    def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0):
1050        # Minimum Qx value [A-1]
1051        self.x_min = x_min
1052        # Maximum Qx value [A-1]
1053        self.x_max = x_max
1054        # Minimum Qy value [A-1]
1055        self.y_min = y_min
1056        # Maximum Qy value [A-1]
1057        self.y_max = y_max
1058
1059    def __call__(self, data2D):
1060        """
1061       Find a rectangular 2D region of interest.
1062
1063       :param data2D: Data2D object
1064       :return: mask, 1d array (len = len(data))
1065           with Trues where the data points are inside ROI, otherwise False
1066        """
1067        mask = self._find(data2D)
1068
1069        return mask
1070
1071    def _find(self, data2D):
1072        """
1073        Find a rectangular 2D region of interest.
1074
1075        :param data2D: Data2D object
1076
1077        :return: out, 1d array (length = len(data))
1078           with Trues where the data points are inside ROI, otherwise Falses
1079        """
1080        if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
1081            raise RuntimeError("Boxcut take only plottable_2D objects")
1082        # Get qx_ and qy_data
1083        qx_data = data2D.qx_data
1084        qy_data = data2D.qy_data
1085
1086        # check whether or not the data point is inside ROI
1087        outx = (self.x_min <= qx_data) & (self.x_max > qx_data)
1088        outy = (self.y_min <= qy_data) & (self.y_max > qy_data)
1089
1090        return outx & outy
1091
1092################################################################################
1093
1094class Sectorcut(object):
1095    """
1096    Defines a sector (major + minor) region on a 2D data set.
1097    The sector is defined by phi_min, phi_max,
1098    where phi_min and phi_max are defined by the right
1099    and left lines wrt central line.
1100
1101    Phi_min and phi_max are given in units of radian
1102    and (phi_max-phi_min) should not be larger than pi
1103    """
1104
1105    def __init__(self, phi_min=0, phi_max=math.pi):
1106        self.phi_min = phi_min
1107        self.phi_max = phi_max
1108
1109    def __call__(self, data2D):
1110        """
1111        Find a rectangular 2D region of interest.
1112
1113        :param data2D: Data2D object
1114
1115        :return: mask, 1d array (len = len(data))
1116
1117        with Trues where the data points are inside ROI, otherwise False
1118        """
1119        mask = self._find(data2D)
1120
1121        return mask
1122
1123    def _find(self, data2D):
1124        """
1125        Find a rectangular 2D region of interest.
1126
1127        :param data2D: Data2D object
1128
1129        :return: out, 1d array (length = len(data))
1130
1131        with Trues where the data points are inside ROI, otherwise Falses
1132        """
1133        if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
1134            raise RuntimeError("Sectorcut take only plottable_2D objects")
1135        Pi = math.pi
1136        # Get data
1137        qx_data = data2D.qx_data
1138        qy_data = data2D.qy_data
1139
1140        # get phi from data
1141        phi_data = np.arctan2(qy_data, qx_data)
1142
1143        # Get the min and max into the region: -pi <= phi < Pi
1144        phi_min_major = flip_phi(self.phi_min + Pi) - Pi
1145        phi_max_major = flip_phi(self.phi_max + Pi) - Pi
1146        # check for major sector
1147        if phi_min_major > phi_max_major:
1148            out_major = (phi_min_major <= phi_data) + \
1149                (phi_max_major > phi_data)
1150        else:
1151            out_major = (phi_min_major <= phi_data) & (
1152                phi_max_major > phi_data)
1153
1154        # minor sector
1155        # Get the min and max into the region: -pi <= phi < Pi
1156        phi_min_minor = flip_phi(self.phi_min) - Pi
1157        phi_max_minor = flip_phi(self.phi_max) - Pi
1158
1159        # check for minor sector
1160        if phi_min_minor > phi_max_minor:
1161            out_minor = (phi_min_minor <= phi_data) + \
1162                (phi_max_minor >= phi_data)
1163        else:
1164            out_minor = (phi_min_minor <= phi_data) & \
1165                (phi_max_minor >= phi_data)
1166        out = out_major + out_minor
1167
1168        return out
Note: See TracBrowser for help on using the repository browser.