Spectrum noise filter

filter.js

//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
 
}