summaryrefslogtreecommitdiff
path: root/silx/math/fit/filters.pyx
diff options
context:
space:
mode:
Diffstat (limited to 'silx/math/fit/filters.pyx')
-rw-r--r--silx/math/fit/filters.pyx411
1 files changed, 411 insertions, 0 deletions
diff --git a/silx/math/fit/filters.pyx b/silx/math/fit/filters.pyx
new file mode 100644
index 0000000..1a7aa3b
--- /dev/null
+++ b/silx/math/fit/filters.pyx
@@ -0,0 +1,411 @@
+# coding: utf-8
+#/*##########################################################################
+# Copyright (C) 2016-2017 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`
+
+API documentation:
+-------------------
+
+"""
+
+__authors__ = ["P. Knobel"]
+__license__ = "MIT"
+__date__ = "15/05/2017"
+
+import logging
+import numpy
+
+_logger = logging.getLogger(__name__)
+
+cimport cython
+cimport 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 *Miroslav
+ Morhac et al. Nucl. Instruments and Methods in Physics Research A401
+ (1997) 113-132*.
+
+ The original idea for 1D and the low-statistics-digital-filter (lsdf) come
+ from *C.G. Ryan et al. Nucl. Instruments and Methods in Physics Research
+ B34 (1988) 396-402*.
+
+ :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
+ *Miroslav Morhac et al. Nucl. Instruments and Methods in Physics Research
+ A401 (1997) 113-132.*
+
+ :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 2D described in
+ *Miroslav Morhac et al. Nucl. Instruments and Methods in Physics Research
+ A401 (1997) 113-132.*
+
+ :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)