Source code for ccpn.core.lib.PeakPickers.PeakPicker1D

"""
Simple 1D PeakPicker; for testing only
"""
#=========================================================================================
# Licence, Reference and Credits
#=========================================================================================
__copyright__ = "Copyright (C) CCPN project (https://www.ccpn.ac.uk) 2014 - 2022"
__credits__ = ("Ed Brooksbank, Joanna Fox, Victoria A Higman, Luca Mureddu, Eliza Płoskoń",
               "Timothy J Ragan, Brian O Smith, Gary S Thompson & Geerten W Vuister")
__licence__ = ("CCPN licence. See http://www.ccpn.ac.uk/v3-software/downloads/license",
               )
__reference__ = ("Skinner, S.P., Fogh, R.H., Boucher, W., Ragan, T.J., Mureddu, L.G., & Vuister, G.W.",
                 "CcpNmr AnalysisAssign: a flexible platform for integrated NMR analysis",
                 "J.Biomol.Nmr (2016), 66, 111-124, http://doi.org/10.1007/s10858-016-0060-y"
                )
#=========================================================================================
# Last code modification
#=========================================================================================
__modifiedBy__ = "$modifiedBy: Luca Mureddu $"
__dateModified__ = "$dateModified: 2022-04-04 15:19:15 +0100 (Mon, April 04, 2022) $"
__version__ = "$Revision: 3.1.0 $"
#=========================================================================================
# Created
#=========================================================================================
__author__ = "$Author: Luca Mureddu $"
__date__ = "$Date: 2022-02-02 14:08:56 +0000 (Wed, February 02, 2022) $"
#=========================================================================================
# Start of code
#=========================================================================================

from ccpn.core.lib.PeakPickers.PeakPickerABC import PeakPickerABC, SimplePeak
from numba import jit
import numpy as np
# from ccpn.framework.Application import getApplication

@jit(nopython=True, nogil=True)
def _find1DMaxima(y, x, positiveThreshold, negativeThreshold=None, findNegative=False):
    """
    from https://gist.github.com/endolith/250860#file-readme-md which was translated from
    http://billauer.co.il/peakdet.html Eli Billauer, 3.4.05.
    Explicitly not copyrighted and any uses allowed.
    """
    maxtab = []
    mintab = []
    mn, mx = np.Inf, -np.Inf
    mnpos, mxpos = np.NaN, np.NaN
    lookformax = True
    if negativeThreshold is None: negativeThreshold = 0

    for i in np.arange(len(y)):
        this = y[i]
        if this > mx:
            mx = this
            mxpos = x[i]
        if this < mn:
            mn = this
            mnpos = x[i]
        if lookformax:
            if not findNegative:  # just positives
                this = abs(this)
            if this < mx - positiveThreshold:
                maxtab.append((float(mxpos), float(mx)))
                mn = this
                mnpos = x[i]
                lookformax = False
        else:
            if this > mn + positiveThreshold:
                mintab.append((float(mnpos), float(mn)))
                mx = this
                mxpos = x[i]
                lookformax = True
    filteredNeg = []
    for p in mintab:
        pos, height = p
        if height <= negativeThreshold:
            filteredNeg.append(p)
    filtered = []
    for p in maxtab:
        pos, height = p
        if height >= positiveThreshold:
            filtered.append(p)
    return filtered, filteredNeg

[docs]class PeakPicker1D(PeakPickerABC): """A peak picker based on Eli Billauer, 3.4.05. algorithm (see _findMaxima function). """ peakPickerType = "PeakPicker1D" onlyFor1D = True def __init__(self, spectrum): super().__init__(spectrum=spectrum) self.noise = None application = spectrum.project.application self._doNegativePeaks = application.preferences.general.negativePeakPick1D def _setThresholds(self): # first make sure the noiseLevel is set for the spectrum. Don't change this. if not self.spectrum.noiseLevel: self.spectrum.noiseLevel = self.spectrum.estimateNoise() if not self.spectrum.negativeNoiseLevel: self.spectrum.negativeNoiseLevel = -self.spectrum.noiseLevel if not self.positiveThreshold: self.positiveThreshold = self.spectrum.noiseLevel if not self.negativeThreshold: self.negativeThreshold = self.spectrum.negativeNoiseLevel # don't pick below the noise level. if self.positiveThreshold <= self.spectrum.noiseLevel: self.positiveThreshold = self.spectrum.noiseLevel if abs(self.negativeThreshold) <= abs(self.spectrum.negativeNoiseLevel): self.negativeThreshold = self.spectrum.negativeNoiseLevel def _isHeightWithinIntesityLimits(self, height): """ check if value is within the intensity limits. This is a different check from the noise Thresholds. Used when picking in a Gui Spectrum Display to make sure is picking only within the drawn box. """ if self._intensityLimits is None or len(self._intensityLimits)==0: self._intensityLimits = (np.inf, -np.inf) value = min(self._intensityLimits) < height < max(self._intensityLimits) return value def _isPositionWithinLimits(self, pointValue): withinLimits = [] ppmValue = self.spectrum.point2ppm(pointValue, axisCode=self.spectrum.axisCodes[0]) excludePpmRegions = self._excludePpmRegions.get(self.spectrum.axisCodes[0], [[0,0],]) # Default ER [0,0] for limits in excludePpmRegions: if len(limits)>0: value = min(limits) < ppmValue < max(limits) withinLimits.append(not value) return all(withinLimits)
[docs] def findPeaks(self, data): peaks = [] x = np.arange(int(self.spectrum.referencePoints[0]),len(data)) self._setThresholds() maxValues, minValues = _find1DMaxima(y=data, x=x, positiveThreshold=self.positiveThreshold, negativeThreshold=self.negativeThreshold, findNegative=self._doNegativePeaks) for position, height in maxValues: if self._isHeightWithinIntesityLimits(height) and self._isPositionWithinLimits(position): pk = SimplePeak(points=(float(position),), height=float(height)) peaks.append(pk) if self._doNegativePeaks: for position, height in minValues: if self._isHeightWithinIntesityLimits(height) and self._isPositionWithinLimits(position): pk = SimplePeak(points=(float(position),), height=float(height)) peaks.append(pk) return peaks
PeakPicker1D._registerPeakPicker()