Spectrum noise filter
//Implementation of Fastie's macro for identifying emission or absorption lines | |
//see: https://github.com/jywarren/spectral-workbench/issues/184 | |
/* | |
Description: | |
To objectively determine the wavelengths of emission or absorption lines, | |
you have to separate signal from noise in a graph of wavelength intensities. | |
Especially when looking for peaks or valleys that disrupt a continuous (black body radiation) spectrum, | |
the low frequency (maybe > 20 nm) information becomes noise. | |
To eliminate it, fit a cubic spline or other easily adjustable curve to each color channel separately. | |
Compute the deviations of each channel from its curve and average the three new curves into a single curve, | |
ignoring the parts of each curve with original values near zero | |
(e.g., so the red curve does not contribute to the mean in the blue part of the spectrum). | |
Fit a new smoothing curve to that one and set a threshold -- | |
any deviation from the smoothing curve greater than x is a signal | |
(emission peak or absorption valley). | |
Notes: | |
The parameters use pixel units, so they will need to be tuned to | |
individual spectrums | |
(because different spectrums have different pixel to spectrumlength ratios) | |
*/ | |
setup: function() { | |
//SMOOTHING_RANGE is the pixel width the the blur filter (i.e. running average) | |
//that is used to smooth the curves. | |
//It should be odd or else filter will be biased to one side. | |
var SMOOTHING_RANGE = 15; | |
//This sets how high up to draw the emission absorbtion curve. | |
var EA_CURVE_Y_OFFSET = 50; | |
//In order to make the peaks of the ea curve clearer, it is amplified by this factor. | |
var EA_CURVE_SCALE_FACTOR = 20; | |
//This set the first row that gets sampled from the spectrum as a percent of the image height. | |
var IMAGE_SAMPLE_PERCENT_LOCATION = 0.5; | |
//Instead of just sampling one row from the spectrum we can average a few rows | |
//to reduce noice. This sets how many we use. | |
var BLUR_ROWS = 5; | |
//When we do peak detection, this sets how far apart peaks have to be. | |
//Only the maximum peak within the range will be highlighted. | |
var PEAK_RANGE = 40; | |
//When doing peak detection you can also filter out peaks below the given magnitude. | |
var PEAK_THRESH = 10; | |
//The gradient parameters are used to find ends of rgb components. | |
//This is done by finding the left-most and right-most points where the slope | |
//is greater than MAX_GRADIENT. | |
//GRADIENT_KERNEL_WIDTH essentially smooths the curve for the slope detector | |
//so that it doesn't detect noise as points with a steep slope. | |
//The width must be odd or else the slope will be biased. | |
//The formula ensures that it is odd. | |
var MAX_GRADIENT = 2; | |
var GRADIENT_KERNEL_WIDTH = Math.floor(SMOOTHING_RANGE / 2) * 2 + 1; | |
//Apply the given kernel filter to graph_data at the given index and return its value. | |
//If the filter runs off the edge of the data it will treat those values as 0s, | |
//unless repeatEdgePixels is on, in which case they will be the value of the last data. | |
var applyFilter = function(kernel, graph_data, index, repeatEdgePixels){ | |
var kernelCenter = Math.floor(kernel.length / 2); | |
var sum = 0; | |
$.each(kernel, function(index2, kvalue) { | |
var loc = index + index2 - kernelCenter; | |
if(loc < 0 || loc >= graph_data.length) { | |
if(repeatEdgePixels) { | |
sum += graph_data[Math.min(Math.max(0, loc), graph_data.length - 1)][1] * kvalue; | |
} | |
} else { | |
sum += graph_data[loc][1] * kvalue; | |
} | |
}); | |
//console.log('idx: '+index+', '+sum); | |
return sum; | |
}; | |
var buildBlurKernel = function(width){ | |
var mykernel = []; | |
var i = 0; | |
while(i<width){ | |
i++; | |
mykernel.push(1/width); | |
} | |
return mykernel; | |
}; | |
var buildGradientKernel = function(width){ | |
if(width % 2 !== 1) { | |
alert("Gradient kernel should have an odd width"); | |
} | |
var mykernel = []; | |
var i = 0; | |
var center = Math.floor(width/2); | |
while(i<width){ | |
i++; | |
if(i === center) { | |
mykernel.push(0); | |
} else{ | |
mykernel.push(1/center * (i - center) / Math.abs(i - center)); | |
} | |
} | |
return mykernel; | |
}; | |
//1d blurred gradient filter | |
var gradientKernel = buildGradientKernel(GRADIENT_KERNEL_WIDTH); | |
//Find where the real data starts by looking for the first point where the average slope | |
//surpasses MAX_GRADIENT | |
var findStart = function(graph_data){ | |
//doesn't take into account the wavelengths because the data is | |
//taken from evenly spaced pixels | |
for(var index = 0; index<graph_data.length; index++) { | |
if(Math.abs(applyFilter(gradientKernel, graph_data, index, true)) > MAX_GRADIENT) { | |
//console.log('idx: '+index+', '+graph_data[index]); | |
return index; | |
} | |
} | |
console.error("Failed to find start of color component."); | |
return 0; | |
}; | |
//copy-pasta of findStart, but in reverse. | |
var findEnd = function(graph_data){ | |
//doesn't take into account the wavelengths because the data is | |
//taken from evenly spaced pixels | |
for(var index = graph_data.length - 1; index>=0; index--) { | |
if(Math.abs(applyFilter(gradientKernel, graph_data, index, true)) > MAX_GRADIENT) { | |
//console.log('idx: '+index+', '+graph_data[index]); | |
return index; | |
} | |
} | |
console.error("Failed to find end of color component."); | |
return graph_data.length - 1; | |
}; | |
//Generate a smooth version of the graph_data curve | |
//by applying a blur convolution filter | |
var generateSmoothingCurve = function(graph_data){ | |
var mykernel = buildBlurKernel(SMOOTHING_RANGE); | |
//The jQuery map function is really wierd. | |
//The callback args are reversed from $.each, and if you return a list, it gets flattened. | |
//Including a good functional programming library would be nice | |
//because the native js map/reduce functions are not cross browser. | |
return $.map(graph_data, function(point, index){ | |
return [[point[0], applyFilter(mykernel, graph_data, index, true)]]; | |
}); | |
}; | |
//curveA - curveB | |
var curveDiff = function(curveA, curveB){ | |
return $.map(curveA, function(pointA, index){ | |
var pointB = curveB[index]; | |
console.assert(pointA[0] === pointB[0]); | |
return [[pointA[0], pointA[1] - pointB[1]]]; | |
}); | |
}; | |
//Create a curve as long as baseCurve with | |
//a uniform slope of the given value. | |
var createHCurve = function(baseCurve, value){ | |
return $.map(baseCurve, function(point){ | |
return [[point[0], value]] | |
}); | |
}; | |
//Apply the given kernel to the entire curve. | |
var convolve = function(mykernel, curve, repeatEdgePixels){ | |
return $.map(curve, function(point, index){ | |
return [[point[0], applyFilter(mykernel, curve, index, repeatEdgePixels)]] | |
}); | |
} | |
//Find the maximum and minimum points on the curve, | |
//filtering out those that are to close together or too close to the middle. | |
var findPeaks = function(curve, repeatEdgePixels){ | |
var peaks = []; | |
$.each(curve, function(index, pointA){ | |
if(Math.abs(pointA[1] - EA_CURVE_Y_OFFSET) < PEAK_THRESH) { | |
return; | |
} | |
var isMax = true; | |
var isMin = true; | |
$.each(curve.slice(index - (PEAK_RANGE / 2), index + (PEAK_RANGE / 2)), function(index2, pointB){ | |
if(pointA === pointB) return; | |
if(pointA[1] < pointB[1]) { | |
isMax = false; | |
} else if(pointA[1] > pointB[1]){ | |
isMin = false; | |
} else { | |
isMax = isMin = false; | |
} | |
}); | |
if(isMax || isMin) peaks.push(index); | |
}); | |
return peaks; | |
} | |
//vertically blur the image | |
var orig_img = $('#image')[0]; | |
var canvas = $("<canvas>")[0]; | |
var img = new Image(); | |
img.src = orig_img.src; | |
canvas.width = img.width; | |
canvas.height = img.height; | |
var ctx = canvas.getContext("2d"); | |
ctx.drawImage(img,0,0); | |
var pixels = ctx.getImageData(0, IMAGE_SAMPLE_PERCENT_LOCATION * img.height,canvas.width,BLUR_ROWS).data; | |
//reset data | |
$W.data = $W.data.slice(0,1); | |
// now, sum pixels & cheat: just overwrite existing calibrated (or uncalibrated) pixel rows | |
$.each($W.spectrum['lines'],function(index,line) { | |
// clear out existing data | |
line['r'] = 0; | |
line['g'] = 0; | |
line['b'] = 0; | |
// for each row of data, add new data | |
for (var i=0;i<BLUR_ROWS;i++) { | |
line['r'] += pixels[i*canvas.width*4+index*4]; | |
line['g'] += pixels[i*canvas.width*4+index*4+1]; | |
line['b'] += pixels[i*canvas.width*4+index*4+2]; | |
} | |
// divide to get average per channel | |
line['r'] = line['r']/BLUR_ROWS; | |
line['g'] = line['g']/BLUR_ROWS; | |
line['b'] = line['b']/BLUR_ROWS; | |
// average colors to get overall | |
line['average'] = (line['r'] + line['g'] + line['b'])/3; | |
}); | |
var annotatedCurves = []; | |
$.each(['r','g','b'],function(index,channel) { | |
// graph smoothed data: | |
var graph_data = []; | |
$.each($W.spectrum.lines,function(index,line) { | |
if (line.wavelength == null) { | |
line.wavelength = index | |
} | |
graph_data.push([line['wavelength'], line[channel]/2.55]) | |
}); | |
annotatedCurves.push({ | |
label: channel, | |
data: graph_data, | |
gradient: convolve([-1, 1], graph_data, true) | |
}); | |
}); | |
var filterCurves = function(annotatedCurves){ | |
var THRESH = 1; | |
$.each(annotatedCurves[0].data, function(index, point){ | |
var A = annotatedCurves[0].gradient[index][1]; | |
var B = annotatedCurves[1].gradient[index][1]; | |
var C = annotatedCurves[2].gradient[index][1]; | |
if(Math.abs(A - B) + Math.abs(A - C) + Math.abs(B - C) < THRESH) { | |
annotatedCurves[0].gradient[index][1] -= (A + B + C) / 3; | |
annotatedCurves[1].gradient[index][1] -= (A + B + C) / 3; | |
annotatedCurves[2].gradient[index][1] -= (A + B + C) / 3; | |
} | |
}); | |
$.each(annotatedCurves, function(index, aCurve){ | |
$.each(aCurve.gradient, function(index, gradient){ | |
var prevValue = index === 0 ? 0 : aCurve.data[index - 1][1]; | |
aCurve.data[index][1] = prevValue + gradient[1]; | |
}); | |
}); | |
}; | |
var combineCurves = function(annotatedCurves){ | |
//make a curve at 0 | |
var result = createHCurve(annotatedCurves[0].diffFromSmooth, 0); | |
$.each(annotatedCurves, function(index, aCurve){ | |
$.each(aCurve.diffFromSmooth, function(index2, point){ | |
var closestBound = Math.abs(index2 - aCurve.bound[0]) < Math.abs(aCurve.bound[1] < index2) ? aCurve.bound[0] : aCurve.bound[1]; | |
if(index2 < aCurve.bound[0]){ | |
//Make the cure fade in/out when it goes out of bounds so it doesn't create discontinuities. | |
result[index2][1] += point[1] / annotatedCurves.length / Math.abs(index2 - aCurve.bound[0]); | |
} else if(aCurve.bound[1] < index2) { | |
//Make the cure fade in/out when it goes out of bounds so it doesn't create discontinuities. | |
result[index2][1] += point[1] / annotatedCurves.length / Math.abs(index2 - aCurve.bound[1]); | |
} else { | |
result[index2][1] += point[1] / annotatedCurves.length; | |
} | |
}); | |
}); | |
return result; | |
}; | |
filterCurves(annotatedCurves); | |
$W.data.push({ | |
label: 'combinedCurve', | |
data: createHCurve(annotatedCurves[0].data, -EA_CURVE_Y_OFFSET) | |
}); | |
//colors of each curve are recorded in the colors | |
//array at the index of their corresponding curve in $W.data | |
flotoptions.colors = ["#0ff", "#ff0"]; | |
$.each(annotatedCurves, function(idx, aCurve){ | |
var color; | |
if(aCurve.label === 'r'){ | |
color = "#f00"; | |
} else if(aCurve.label === 'g'){ | |
color = "#0f0"; | |
} else if(aCurve.label === 'b'){ | |
color = "#00f"; | |
} | |
flotoptions.colors[idx+2] = color; | |
$W.data.push(aCurve); | |
}); | |
$W.plot = $.plot($("#graph"), $W.data, flotoptions); | |
}, | |
draw: function() { | |
// code to run every frame | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment