diff options
Diffstat (limited to 'src/silx/math/fit/filters.pyx')
| -rw-r--r-- | src/silx/math/fit/filters.pyx | 416 |
1 files changed, 416 insertions, 0 deletions
diff --git a/src/silx/math/fit/filters.pyx b/src/silx/math/fit/filters.pyx new file mode 100644 index 0000000..da1f6f5 --- /dev/null +++ b/src/silx/math/fit/filters.pyx @@ -0,0 +1,416 @@ +# coding: utf-8 +#/*########################################################################## +# Copyright (C) 2016-2018 European Synchrotron Radiation Facility +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +#############################################################################*/ +"""This module provides background extraction functions and smoothing +functions. These functions are extracted from PyMca module SpecFitFuns. + +Index of background extraction functions: +------------------------------------------ + + - :func:`strip` + - :func:`snip1d` + - :func:`snip2d` + - :func:`snip3d` + +Smoothing functions: +-------------------- + + - :func:`savitsky_golay` + - :func:`smooth1d` + - :func:`smooth2d` + - :func:`smooth3d` + +References: +----------- + +.. [Morhac97] Miroslav Morháč et al. + Background elimination methods for multidimensional coincidence γ-ray spectra. + Nucl. Instruments and Methods in Physics Research A401 (1997) 113-132. + https://doi.org/10.1016/S0168-9002(97)01023-1 + +.. [Ryan88] C.G. Ryan et al. + SNIP, a statistics-sensitive background treatment for the quantitative analysis of PIXE spectra in geoscience applications. + Nucl. Instruments and Methods in Physics Research B34 (1988) 396-402*. + https://doi.org/10.1016/0168-583X(88)90063-8 + +API documentation: +------------------- + +""" + +__authors__ = ["P. Knobel"] +__license__ = "MIT" +__date__ = "15/05/2017" + +import logging +import numpy + +_logger = logging.getLogger(__name__) + +cimport cython +cimport silx.math.fit.filters_wrapper as filters_wrapper + + +def strip(data, w=1, niterations=1000, factor=1.0, anchors=None): + """Extract background from data using the strip algorithm, as explained at + http://pymca.sourceforge.net/stripbackground.html. + + In its simplest implementation it is just as an iterative procedure + depending on two parameters. These parameters are the strip background + width ``w``, and the number of iterations. At each iteration, if the + contents of channel ``i``, ``y(i)``, is above the average of the contents + of the channels at ``w`` channels of distance, ``y(i-w)`` and + ``y(i+w)``, ``y(i)`` is replaced by the average. + At the end of the process we are left with something that resembles a spectrum + in which the peaks have been stripped. + + :param data: Data array + :type data: numpy.ndarray + :param w: Strip width + :param niterations: number of iterations + :param factor: scaling factor applied to the average of ``y(i-w)`` and + ``y(i+w)`` before comparing to ``y(i)`` + :param anchors: Array of anchors, indices of points that will not be + modified during the stripping procedure. + :return: Data with peaks stripped away + """ + cdef: + double[::1] input_c + double[::1] output + long[::1] anchors_c + + if not isinstance(data, numpy.ndarray): + if not hasattr(data, "__len__"): + raise TypeError("data must be a sequence (list, tuple) " + + "or a numpy array") + data_shape = (len(data), ) + else: + data_shape = data.shape + + input_c = numpy.array(data, + copy=True, + dtype=numpy.float64, + order='C').reshape(-1) + + output = numpy.empty(shape=(input_c.size,), + dtype=numpy.float64) + + if anchors is not None and len(anchors): + # numpy.int_ is the same as C long (http://docs.scipy.org/doc/numpy/user/basics.types.html) + anchors_c = numpy.array(anchors, + copy=False, + dtype=numpy.int_, + order='C') + len_anchors = anchors_c.size + else: + # Make a dummy length-1 array, because if I use shape=(0,) I get the error + # IndexError: Out of bounds on buffer access (axis 0) + anchors_c = numpy.empty(shape=(1,), + dtype=numpy.int_) + len_anchors = 0 + + + status = filters_wrapper.strip(&input_c[0], input_c.size, + factor, niterations, w, + &anchors_c[0], len_anchors, &output[0]) + + return numpy.asarray(output).reshape(data_shape) + + +def snip1d(data, snip_width): + """Estimate the baseline (background) of a 1D data vector by clipping peaks. + + Implementation of the algorithm SNIP in 1D is described in [Morhac97]_. + The original idea for 1D and the low-statistics-digital-filter (lsdf) comes + from [Ryan88]_. + + :param data: Data array, preferably 1D and of type *numpy.float64*. + Else, the data array will be flattened and converted to + *dtype=numpy.float64* prior to applying the snip filter. + :type data: numpy.ndarray + :param snip_width: Width of the snip operator, in number of samples. + A sample will be iteratively compared to it's neighbors up to a + distance of ``snip_width`` samples. This parameters has a direct + influence on the speed of the algorithm. + :type width: int + :return: Baseline of the input array, as an array of the same shape. + :rtype: numpy.ndarray + """ + cdef: + double[::1] data_c + + if not isinstance(data, numpy.ndarray): + if not hasattr(data, "__len__"): + raise TypeError("data must be a sequence (list, tuple) " + + "or a numpy array") + data_shape = (len(data), ) + else: + data_shape = data.shape + + data_c = numpy.array(data, + copy=True, + dtype=numpy.float64, + order='C').reshape(-1) + + filters_wrapper.snip1d(&data_c[0], data_c.size, snip_width) + + return numpy.asarray(data_c).reshape(data_shape) + + +def snip2d(data, snip_width): + """Estimate the baseline (background) of a 2D data signal by clipping peaks. + + Implementation of the algorithm SNIP in 2D described in [Morhac97]_. + + :param data: 2D array + :type data: numpy.ndarray + :param width: Width of the snip operator, in number of samples. A wider + snip operator will result in a smoother result (lower frequency peaks + will be clipped), and a longer computation time. + :type width: int + :return: Baseline of the input array, as an array of the same shape. + :rtype: numpy.ndarray + """ + cdef: + double[::1] data_c + + if not isinstance(data, numpy.ndarray): + if not hasattr(data, "__len__") or not hasattr(data[0], "__len__"): + raise TypeError("data must be a 2D sequence (list, tuple) " + + "or a 2D numpy array") + nrows = len(data) + ncolumns = len(data[0]) + data_shape = (len(data), len(data[0])) + + else: + data_shape = data.shape + nrows = data_shape[0] + if len(data_shape) == 2: + ncolumns = data_shape[1] + else: + raise TypeError("data array must be 2-dimensional") + + data_c = numpy.array(data, + copy=True, + dtype=numpy.float64, + order='C').reshape(-1) + + filters_wrapper.snip2d(&data_c[0], nrows, ncolumns, snip_width) + + return numpy.asarray(data_c).reshape(data_shape) + + +def snip3d(data, snip_width): + """Estimate the baseline (background) of a 3D data signal by clipping peaks. + + Implementation of the algorithm SNIP in 3D described in [Morhac97]_. + + :param data: 3D array + :type data: numpy.ndarray + :param width: Width of the snip operator, in number of samples. A wider + snip operator will result in a smoother result (lower frequency peaks + will be clipped), and a longer computation time. + :type width: int + + :return: Baseline of the input array, as an array of the same shape. + :rtype: numpy.ndarray + """ + cdef: + double[::1] data_c + + if not isinstance(data, numpy.ndarray): + if not hasattr(data, "__len__") or not hasattr(data[0], "__len__") or\ + not hasattr(data[0][0], "__len__"): + raise TypeError("data must be a 3D sequence (list, tuple) " + + "or a 3D numpy array") + nx = len(data) + ny = len(data[0]) + nz = len(data[0][0]) + data_shape = (len(data), len(data[0]), len(data[0][0])) + else: + data_shape = data.shape + nrows = data_shape[0] + if len(data_shape) == 3: + nx = data_shape[0] + ny = data_shape[1] + nz = data_shape[2] + else: + raise TypeError("data array must be 3-dimensional") + + data_c = numpy.array(data, + copy=True, + dtype=numpy.float64, + order='C').reshape(-1) + + filters_wrapper.snip3d(&data_c[0], nx, ny, nz, snip_width) + + return numpy.asarray(data_c).reshape(data_shape) + + +def savitsky_golay(data, npoints=5): + """Smooth a curve using a Savitsky-Golay filter. + + :param data: Input data + :type data: 1D numpy array + :param npoints: Size of the smoothing operator in number of samples + Must be between 3 and 100. + :return: Smoothed data + """ + cdef: + double[::1] data_c + double[::1] output + + data_c = numpy.array(data, + dtype=numpy.float64, + order='C').reshape(-1) + + output = numpy.empty(shape=(data_c.size,), + dtype=numpy.float64) + + status = filters_wrapper.SavitskyGolay(&data_c[0], data_c.size, + npoints, &output[0]) + + if status: + _logger.error("Smoothing failed. Check that npoints is greater " + + "than 3 and smaller than 100.") + + return numpy.asarray(output).reshape(data.shape) + + +def smooth1d(data): + """Simple smoothing for 1D data. + + For a data array :math:`y` of length :math:`n`, the smoothed array + :math:`ys` is calculated as a weighted average of neighboring samples: + + :math:`ys_0 = 0.75 y_0 + 0.25 y_1` + + :math:`ys_i = 0.25 (y_{i-1} + 2 y_i + y_{i+1})` for :math:`0 < i < n-1` + + :math:`ys_{n-1} = 0.25 y_{n-2} + 0.75 y_{n-1}` + + + :param data: 1D data array + :type data: numpy.ndarray + :return: Smoothed data + :rtype: numpy.ndarray(dtype=numpy.float64) + """ + cdef: + double[::1] data_c + + if not isinstance(data, numpy.ndarray): + if not hasattr(data, "__len__"): + raise TypeError("data must be a sequence (list, tuple) " + + "or a numpy array") + data_shape = (len(data), ) + else: + data_shape = data.shape + + data_c = numpy.array(data, + copy=True, + dtype=numpy.float64, + order='C').reshape(-1) + + filters_wrapper.smooth1d(&data_c[0], data_c.size) + + return numpy.asarray(data_c).reshape(data_shape) + + +def smooth2d(data): + """Simple smoothing for 2D data: + :func:`smooth1d` is applied succesively along both axis + + :param data: 2D data array + :type data: numpy.ndarray + :return: Smoothed data + :rtype: numpy.ndarray(dtype=numpy.float64) + """ + cdef: + double[::1] data_c + + if not isinstance(data, numpy.ndarray): + if not hasattr(data, "__len__") or not hasattr(data[0], "__len__"): + raise TypeError("data must be a 2D sequence (list, tuple) " + + "or a 2D numpy array") + nrows = len(data) + ncolumns = len(data[0]) + data_shape = (len(data), len(data[0])) + + else: + data_shape = data.shape + nrows = data_shape[0] + if len(data_shape) == 2: + ncolumns = data_shape[1] + else: + raise TypeError("data array must be 2-dimensional") + + data_c = numpy.array(data, + copy=True, + dtype=numpy.float64, + order='C').reshape(-1) + + filters_wrapper.smooth2d(&data_c[0], nrows, ncolumns) + + return numpy.asarray(data_c).reshape(data_shape) + + +def smooth3d(data): + """Simple smoothing for 3D data: + :func:`smooth2d` is applied on each 2D slice of the data volume along all + 3 axis + + :param data: 2D data array + :type data: numpy.ndarray + :return: Smoothed data + :rtype: numpy.ndarray(dtype=numpy.float64) + """ + cdef: + double[::1] data_c + + if not isinstance(data, numpy.ndarray): + if not hasattr(data, "__len__") or not hasattr(data[0], "__len__") or\ + not hasattr(data[0][0], "__len__"): + raise TypeError("data must be a 3D sequence (list, tuple) " + + "or a 3D numpy array") + nx = len(data) + ny = len(data[0]) + nz = len(data[0][0]) + data_shape = (len(data), len(data[0]), len(data[0][0])) + else: + data_shape = data.shape + nrows = data_shape[0] + if len(data_shape) == 3: + nx = data_shape[0] + ny = data_shape[1] + nz = data_shape[2] + else: + raise TypeError("data array must be 3-dimensional") + + data_c = numpy.array(data, + copy=True, + dtype=numpy.float64, + order='C').reshape(-1) + + filters_wrapper.smooth3d(&data_c[0], nx, ny, nz) + + return numpy.asarray(data_c).reshape(data_shape) |