diff options
Diffstat (limited to 'silx/math/fit/filters.pyx')
| -rw-r--r-- | silx/math/fit/filters.pyx | 416 |
1 files changed, 0 insertions, 416 deletions
diff --git a/silx/math/fit/filters.pyx b/silx/math/fit/filters.pyx deleted file mode 100644 index da1f6f5..0000000 --- a/silx/math/fit/filters.pyx +++ /dev/null @@ -1,416 +0,0 @@ -# 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) |