Skip to content

Instantly share code, notes, and snippets.

@efekarakus

efekarakus/README.md

Last active March 15, 2016 16:09
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save efekarakus/cc7303456841523f37dd to your computer and use it in GitHub Desktop.
Save efekarakus/cc7303456841523f37dd to your computer and use it in GitHub Desktop.
Download ZIP
Find Peaks
Raw
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.

Raw
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;
}));
Raw
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>
Raw
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]}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment