Find Peaks

README.md

This example illustrates how to use the findPeaks API. You can download the library from https://github.com/efekarakus/d3-peaks.

The algorithm is based on "Improved peak detection in mass spectrum by incorporating continuous wavelet transform-based pattern matching" by Pan Du, Warren A. Kibbe and Simon M. Lin. The paper can be found here.

d3-peaks.js

(function (global, factory) {
  typeof exports === 'object' && typeof module !== 'undefined' ? factory(exports) :
  typeof define === 'function' && define.amd ? define(['exports'], factory) :
  (factory((global.d3_peaks = global.d3_peaks || {})));
}(this, function (exports) { 'use strict';

  /**
   * See https://en.wikipedia.org/wiki/Mexican_hat_wavelet
   */
  function ricker() {
    var σ = 1;
    
    var ricker = function(t) {
      var t2 = t*t,
          variance = σ*σ;
      
      var C = 2.0 / ( Math.sqrt(3 * σ) * (Math.pow(Math.PI, 0.25)) );
      var norm = (1.0 - (t2)/(variance));
      var gauss = Math.exp( -(t2) / (2*variance) );
      
      return C*norm*gauss;
    }
    
    ricker.std = function(_) {
      return arguments.length ? (σ = _, ricker) : σ;
    }
    
    /**
     * Range of points to sample from the wavelet. [-reach, reach]
     */
    ricker.reach = function() {
      return 5 * σ;
    }
    
    return ricker;
  };

  function convolve() {
    var kernel = ricker();
    
    /**
     * y[n] = Sum_k{x[k] * h[n-k]}
     * y: output
     * x: input
     * h: smoother
     */
    var convolve = function(signal) {
      var size = signal.length,
          n = -1,
          convolution = new Array(size);
          
      while (++n < size) {
        var y = 0;
        
        var box = boundingBox(n, kernel.reach(), 0, size - 1);
        box.forEach(function(δ) {
          var k = n + δ;
          y += signal[k] * kernel(δ);
        });
        convolution[n] = y;
      }
      
      return convolution;
    };
    
    convolve.kernel = function(_) {
      return arguments.length ? (kernel = _, convolve) : kernel;
    }
    
    function range(reach) {
      reach = +reach;
      var i = -1,
          n = 2*reach + 1,
          range = new Array(n);
      while(++i < n) {
        range[i] = (-reach) + i;
      }
      return range;
    }
    
    function boundingBox(n, reach, lo, hi) {
      for (var i = 1; i <= reach; i++) {
        var left  = n - i,
            right = n + i;
        if (left >= lo && right <= hi) continue;
        return range(i - 1);
      }
      return range(reach);
    }
    
    return convolve;
  };

  function isLocalMaxima(arr, index) {
    var current = arr[index],
        left = arr[index - 1],
        right = arr[index + 1];
        
    if (left !== undefined && right !== undefined) {
      if (current > left && current > right) { return true; }
      else if (current >= left && current > right) { return true; }
      else if (current > left && current >= right) { return true; }
    }
    else if (left !== undefined && current > left) { return true; }
    else if (right !== undefined && current > right) { return true; }
    
    return false;
  }

  /**
   * @param {arr} row in the CWT matrix.
   * @return Array of indices with relative maximas.
   */
  function maximas(arr) {
    var maximas = [];
    var length = arr.length;
    arr.forEach(function(value, index) {
      if (isLocalMaxima(arr, index)) maximas.push({x: index, y: value});
    });
    return maximas;
  };

  function nearestNeighbor(line, maximas, window) {
    var cache = {};
    maximas.forEach(function(d) {
      cache[d.x] = d.y;
    });
    
    var point = line.top();
    for (var i = 0; i <= window; i++) {
      var left = point.x + i;
      var right = point.x - i;
      
      if ( (left in cache) && (right in cache) ) {
        if (cache[left] > cache[right]) {
          return left;
        }
        return right;
      }
      else if (left in cache) {
        return left;
      }
      else if (right in cache) {
        return right;
      }
    }
    return null;
  }

  function percentile(arr, perc) {
    var length = arr.length;
    var index = Math.min(length - 1, Math.ceil(perc * length));
    
    arr.sort(function(a, b) { return a - b; });
    return arr[index];
  }

  function Point(x, y, width) {
    this.x = x;
    this.y = y;
    this.width = width;
    this.snr = undefined;
  }

  Point.prototype.SNR = function(conv) {
    var smoothingFactor = 0.00001;
    var signal = this.y;
    
    var lowerBound = Math.max(0, this.x - this.width);
    var upperBound = Math.min(conv.length, this.x + this.width + 1);
    var neighbors = conv.slice(lowerBound, upperBound);
    var noise = percentile(neighbors, 0.95);
    
    signal += smoothingFactor;
    noise += smoothingFactor;
    this.snr = signal/noise;
    return this.snr;
  }

  Point.prototype.serialize = function() {
    return {index: this.x, width: this.width, snr: this.snr};
  }

  function RidgeLine() {
    this.points = [];
    this.gap = 0;
  }

  /**
   * If the point is valid append it to the ridgeline, and reset the gap.
   * Otherwise, increment the gap and do nothing.
   * 
   * @param {point} Point object.
   */
  RidgeLine.prototype.add = function(point) {
    if (point === null || point === undefined) {
      this.gap += 1;
      return;
    } else {
      this.points.push(point);
      this.gap = 0;
    }
  }

  /**
   * @return {Point} Last point added into the ridgeline.
   */
  RidgeLine.prototype.top = function() {
    return this.points[this.points.length - 1];
  }

  /**
   * @return {number} Length of points on the ridgeline.
   */
  RidgeLine.prototype.length = function() {
    return this.points.length;
  }

  /**
   * @return {boolean} True if the gap in the line is above a threshold. False otherwise.
   */
  RidgeLine.prototype.isDisconnected = function (threshold) {
    return this.gap > threshold;
  }

  /**
   * @param {Array} Smallest scale in the convolution matrix
   */
  RidgeLine.prototype.SNR = function(conv) {
    var maxSnr = Number.NEGATIVE_INFINITY;
    this.points.forEach(function(point) {
      var snr = point.SNR(conv);
      if (snr > maxSnr) maxSnr = snr;
    });
    return maxSnr;
  }

  function findPeaks() {
    var kernel = ricker,
        gapThreshold = 1,
        minLineLength = 1,
        minSNR = 1.0,
        widths = [1];
    
    var findPeaks = function(signal) {
      var M = CWT(signal);
      
      var ridgeLines = initializeRidgeLines(M);
      ridgeLines = connectRidgeLines(M, ridgeLines);
      ridgeLines = filterRidgeLines(signal, ridgeLines);
      
      return peaks(signal, ridgeLines);
    };
    
    /**
     * Smoothing function.
     */
    findPeaks.kernel = function(_) {
      return arguments.length ? (kernel = _, findPeaks) : kernel;
    }
    
    /**
     * Expected widths of the peaks.
     */
    findPeaks.widths = function(_) {
      _.sort(function(a, b) { return a - b; });
      return arguments.length ? (widths = _, findPeaks) : widths;
    }
    
    /**
     * Number of gaps that we allow in the ridge lines.
     */
    findPeaks.gapThreshold = function(_) {
      return arguments.length ? (gapThreshold = _, findPeaks) : gapThreshold;
    }
    
    /**
     * Minimum ridge line length.
     */
    findPeaks.minLineLength = function(_) {
      return arguments.length ? (minLineLength = _, findPeaks) : minLineLength;
    }
    
    /**
     * Minimum signal to noise ratio for the peaks.
     */
    findPeaks.minSNR = function(_) {
      return arguments.length ? (minSNR = _, findPeaks) : minSNR;
    }
    
    var CWT = function(signal) {
      var M = new Array(widths.length);
      widths.forEach(function(width, i) {
        var smoother = kernel()
          .std(width);
        var transform = convolve()
          .kernel(smoother);
        
        var convolution = transform(signal);
        M[i] = convolution;
      });
      return M;
    }
    
    
    var initializeRidgeLines = function(M) {
      var n = widths.length;
      var locals = maximas(M[n - 1], widths[n - 1]);
      var ridgeLines = [];
      locals.forEach(function(d) {
        var point = new Point(d.x, d.y, widths[n - 1]);
        var line = new RidgeLine();
        line.add(point);
        ridgeLines.push(line);
      });
      return ridgeLines;
    }
    
    var connectRidgeLines = function(M, ridgeLines) {
      var n = widths.length;
      for (var row = n - 2; row >= 0; row--) {
        var locals = maximas(M[row], widths[row]);
        var addedLocals = [];
        
        // Find nearest neighbor at next scale and add to the line
        ridgeLines.forEach(function(line, i) {
          var x = nearestNeighbor(line, locals, widths[row]);
          line.add(x === null ? null : new Point(x, M[row][x], widths[row]));
          
          if (x !== null) {
            addedLocals.push(x);
          }
        });
        
        // Remove lines that has exceeded the gap threshold
        ridgeLines = ridgeLines.filter(function(line) {
          return !line.isDisconnected(gapThreshold);
        });
        
        // Add all the unitialized ridge lines
        locals.forEach(function(d) {
          if (addedLocals.indexOf(d.x) !== -1) return;
          
          var point = new Point(d.x, d.y, widths[row]);
          var ridgeLine = new RidgeLine();
          ridgeLine.add(point);
          ridgeLines.push(ridgeLine);
        });
      }
      return ridgeLines;
    }
    
    var filterRidgeLines = function(signal, ridgeLines) {
      var smoother = kernel()
          .std(1.0);
      var transform = convolve()
        .kernel(smoother);
      var convolution = transform(signal);
        
      ridgeLines = ridgeLines.filter(function(line) {
        var snr = line.SNR(convolution);
        return (snr >= minSNR) && (line.length() >= minLineLength);
      });
      return ridgeLines
    }
    
    /**
     * Pick the point on the ridgeline with highest SNR.
     */
    var peaks = function(signal, ridgeLines) {
      var peaks = ridgeLines.map(function(line) {
        var points = line.points;
        var maxValue = Number.NEGATIVE_INFINITY,
            maxPoint = undefined;
        points.forEach(function(point) {
          var y = signal[point.x];
          if (y > maxValue) {
            maxPoint = point;
            maxValue = y;
          }
        });
        return maxPoint.serialize();
      });
      return peaks;
    }
    
    return findPeaks;
  };

  var version = "0.0.1";

  exports.version = version;
  exports.ricker = ricker;
  exports.convolve = convolve;
  exports.findPeaks = findPeaks;

}));

index.html

<!DOCTYPE html>
<meta charset="utf-8">
<title>Find Peaks</title>
<style>
  .axis path,
  .axis line {
    fill: none;
    stroke: black;
    shape-rendering: crispEdges;
  }
  text {
    font-family: sans-serif;
    font-size: 11px;
  }
  
  path {
    fill: none;
    stroke: black;
    stroke-width: 1.5px;
  }
  
  path.cardinal {
    stroke: crimson;
    stroke-width: 3px;
  }
  
  path.area {
    stroke: black;
    stroke-width: 0.5px;
    fill: crimson;
    fill-opacity: 0.2;
  }
  
  circle.peak {
    fill: crimson;
  }
  
</style>
<body>
<script src="//d3js.org/d3.v3.min.js" charset="utf-8"></script>
<script src="d3-peaks.js" charset="utf-8"></script>
<script>
  var width = 920,
      height = 464,
      margin = {top: 10, left: 10, bottom: 20, right: 10};
  
  var ricker = d3_peaks.ricker;
  var findPeaks = d3_peaks.findPeaks()
    .kernel(ricker)
    .gapThreshold(1)
    .minLineLength(2)
    .minSNR(1.0)
    .widths([1,2,3]);
  
  d3.json("signal.json", function(error, json) {
    if (error) console.warn(error);
    
    var signal = json.signal;
    var peaks = findPeaks(signal);
    
    var svg = d3.select("body").append("svg")
      .attr("width", width + margin.left + margin.right)
      .attr("height", height + margin.top + margin.bottom);
    var g = svg.append("g")
      .attr("transform", "translate(" + margin.left + "," + margin.top + ")");
      
    var x = d3.scale.linear()
      .domain ([0, signal.length])
      .range([20, width]);
      
    var y = d3.scale.linear()
      .domain([0, d3.max(signal)])
      .range([height, 0]);
      
    var xAxis = d3.svg.axis()
      .scale(x)
      .ticks(20);
      
    var line = d3.svg.line()
      .x(function(d, i) { return x(i); })
      .y(function(d) { return y(d); });
      
    var area = d3.svg.area()
      .x(function(d) { return x(d.x); })
      .y0(height)
      .y1(function(d) { return y(d.y); });
      
      
    g.append("g")
      .datum(signal)
      .append("path")
      .attr("d", line);
      
    g.append("g")
      .attr("class", "axis")
      .attr("transform", "translate(0," + y(0) + ")")
      .call(xAxis);
      
    g.selectAll(".peak")
      .data(peaks)
      .enter().append("circle")
        .attr("cx", function(peak) { return x(peak.index); })
        .attr("cy", function(peak) { return y(signal[peak.index]); })
        .attr("r", 4)
        .attr("class", "peak");

    line.interpolate("cardinal")
      .x(function(d) { return x(d.index); })
      .y(function(d) { return y(signal[d.index]); });
      
    peaks.sort(function(a, b) { return a.index - b.index; });
    g.append("g")
      .datum(peaks)
      .append("path")
        .attr("d", line)
        .attr("class", "cardinal");

    var areas = [];
    peaks.forEach(function(peak) {
      var width = peak.width;
      var lowerBound = Math.max(0, peak.index - width);
      var upperBound = Math.min(signal.length, peak.index + width + 1);
      
      var a = [];
      for (var i = lowerBound; i < upperBound; i++) {
        a.push({
          x: i,
          y: signal[i]
        });
      }
      areas.push(a);
    });
    
    areas.forEach(function(data) {
      g.append("g")
        .datum(data)
        .append("path")
          .attr("d", area)
          .attr("class", "area");
    })
  })
</script>
</body>

signal.json

{"signal": [ 1, 10, 12, 12, 8, 15, 10, 4, 13, 11, 8, 7, 8, 9, 3, 10, 8, 5, 9, 6, 9, 7, 17, 5, 8, 6, 4, 5, 3, 5, 9, 6, 3, 8, 1, 7, 4, 8, 7, 5, 3, 6, 3, 13, 6, 6, 9, 6, 8, 9, 5, 6, 6, 5, 9, 6, 5, 7, 7, 7, 7, 6, 6, 7, 6, 4, 7, 8, 6, 18, 6, 11, 12, 7, 19, 13, 2, 14, 16, 4, 11, 10, 11, 10, 8, 6, 5, 6, 5, 6, 6, 6, 4, 5, 3, 7, 6, 7, 5, 5, 6, 3, 7, 2, 3, 1, 3, 4, 3, 2, 4, 4, 1, 2, 1, 6, 1, 1, 1, 4, 3, 1, 5, 1, 2, 3, 2, 6, 2, 2, 3, 2, 2, 2, 1, 4, 7, 6, 5, 7, 7, 6, 4, 12, 9, 7, 5, 2, 4, 7, 5, 4, 4, 4, 2, 5, 3, 3, 3, 5, 7, 4, 4, 4, 2, 4, 1, 3, 2, 3, 4, 3, 6, 6, 4, 4, 8, 11, 8, 9, 3, 3, 3, 5, 6, 4, 7, 9, 7, 13, 9, 7, 12, 7, 14, 11, 8, 5, 2, 6, 6, 5, 5, 2, 4, 2, 4, 6, 6, 14, 16, 16, 18, 14, 3, 18, 26, 2, 23, 26, 8, 18, 16, 3, 14, 27, 4, 10, 14, 20, 2, 15, 22, 8, 9, 11, 15, 8, 8, 15]}

thumbnail.png