Merge pull request #3214 from tschaub/raster

Pixel manipulation with raster sources.

examples/raster.css

@@ -0,0 +1,31 @@
+.rel {
+  position: relative
+}
+
+#plot {
+  pointer-events: none;
+  position: absolute;
+  bottom: 10px;
+  left: 10px;
+}
+
+.bar {
+  pointer-events: auto;
+  fill: #AFAFB9;
+}
+
+.bar.selected {
+  fill: green;
+}
+
+.tip {
+  position: absolute;
+  background: black;
+  color: white;
+  padding: 6px;
+  font-size: 12px;
+  border-radius: 4px;
+  margin-bottom: 10px;
+  display: none;
+  opacity: 0;
+}

examples/raster.html

@@ -0,0 +1,31 @@
+---
+template: example.html
+title: Raster Source
+shortdesc: Demonstrates pixelwise operations with a raster source.
+docs: >
+  <p>
+    This example uses a <code>ol.source.Raster</code> to generate data
+    based on another source.  The raster source accepts any number of
+    input sources (tile or image based) and runs a pipeline of
+    operations on the input pixels.  The return from the final
+    operation is used as the data for the output source.
+  </p>
+  <p>
+    In this case, a single tiled source of imagery is used as input.
+    For each pixel, the Vegetaion Greenness Index
+    (<a href="http://www.tandfonline.com/doi/abs/10.1080/10106040108542184#.Vb90ITBViko">VGI</a>)
+    is calculated from the input pixels.  A second operation colors
+    those pixels based on a threshold value (values above the
+    threshold are green and those below are transparent).
+  </p>
+tags: "raster, pixel"
+resources:
+  - http://d3js.org/d3.v3.min.js
+  - raster.css
+---
+<div class="row-fluid">
+  <div class="span12 rel">
+    <div id="map" class="map"></div>
+    <div id="plot"></div>
+  </div>
+</div>

examples/raster.js

@@ -0,0 +1,200 @@
+// NOCOMPILE
+// this example uses d3 for which we don't have an externs file.
+goog.require('ol.Map');
+goog.require('ol.View');
+goog.require('ol.layer.Image');
+goog.require('ol.layer.Tile');
+goog.require('ol.source.BingMaps');
+goog.require('ol.source.Raster');
+
+var minVgi = 0;
+var maxVgi = 0.25;
+var bins = 10;
+
+
+/**
+ * Calculate the Vegetation Greenness Index (VGI) from an input pixel.  This
+ * is a rough estimate assuming that pixel values correspond to reflectance.
+ * @param {ol.raster.Pixel} pixel An array of [R, G, B, A] values.
+ * @return {number} The VGI value for the given pixel.
+ */
+function vgi(pixel) {
+  var r = pixel[0] / 255;
+  var g = pixel[1] / 255;
+  var b = pixel[2] / 255;
+  return (2 * g - r - b) / (2 * g + r + b);
+}
+
+
+/**
+ * Summarize values for a histogram.
+ * @param {numver} value A VGI value.
+ * @param {Object} counts An object for keeping track of VGI counts.
+ */
+function summarize(value, counts) {
+  var min = counts.min;
+  var max = counts.max;
+  var num = counts.values.length;
+  if (value < min) {
+    // do nothing
+  } else if (value >= max) {
+    counts.values[num - 1] += 1;
+  } else {
+    var index = Math.floor((value - min) / counts.delta);
+    counts.values[index] += 1;
+  }
+}
+
+
+/**
+ * Use aerial imagery as the input data for the raster source.
+ */
+var bing = new ol.source.BingMaps({
+  key: 'Ak-dzM4wZjSqTlzveKz5u0d4IQ4bRzVI309GxmkgSVr1ewS6iPSrOvOKhA-CJlm3',
+  imagerySet: 'Aerial'
+});
+
+
+/**
+ * Create a raster source where pixels with VGI values above a threshold will
+ * be colored green.
+ */
+var raster = new ol.source.Raster({
+  sources: [bing],
+  operation: function(pixels, data) {
+    var pixel = pixels[0];
+    var value = vgi(pixel);
+    summarize(value, data.counts);
+    if (value >= data.threshold) {
+      pixel[0] = 0;
+      pixel[1] = 255;
+      pixel[2] = 0;
+      pixel[3] = 128;
+    } else {
+      pixel[3] = 0;
+    }
+    return pixel;
+  },
+  lib: {
+    vgi: vgi,
+    summarize: summarize
+  }
+});
+raster.set('threshold', 0.1);
+
+function createCounts(min, max, num) {
+  var values = new Array(num);
+  for (var i = 0; i < num; ++i) {
+    values[i] = 0;
+  }
+  return {
+    min: min,
+    max: max,
+    values: values,
+    delta: (max - min) / num
+  };
+}
+
+raster.on('beforeoperations', function(event) {
+  event.data.counts = createCounts(minVgi, maxVgi, bins);
+  event.data.threshold = raster.get('threshold');
+});
+
+raster.on('afteroperations', function(event) {
+  schedulePlot(event.resolution, event.data.counts, event.data.threshold);
+});
+
+var map = new ol.Map({
+  layers: [
+    new ol.layer.Tile({
+      source: bing
+    }),
+    new ol.layer.Image({
+      source: raster
+    })
+  ],
+  target: 'map',
+  view: new ol.View({
+    center: [-9651695, 4937351],
+    zoom: 13,
+    minZoom: 12,
+    maxZoom: 19
+  })
+});
+
+
+var timer = null;
+function schedulePlot(resolution, counts, threshold) {
+  if (timer) {
+    clearTimeout(timer);
+    timer = null;
+  }
+  timer = setTimeout(plot.bind(null, resolution, counts, threshold), 1000 / 60);
+}
+
+var barWidth = 15;
+var plotHeight = 150;
+var chart = d3.select('#plot').append('svg')
+    .attr('width', barWidth * bins)
+    .attr('height', plotHeight);
+
+var chartRect = chart[0][0].getBoundingClientRect();
+
+var tip = d3.select(document.body).append('div')
+    .attr('class', 'tip');
+
+function plot(resolution, counts, threshold) {
+  var yScale = d3.scale.linear()
+      .domain([0, d3.max(counts.values)])
+      .range([0, plotHeight]);
+
+  var bar = chart.selectAll('rect').data(counts.values);
+
+  bar.enter().append('rect');
+
+  bar.attr('class', function(count, index) {
+    var value = counts.min + (index * counts.delta);
+    return 'bar' + (value >= threshold ? ' selected' : '');
+  })
+  .attr('width', barWidth - 2);
+
+  bar.transition().attr('transform', function(value, index) {
+    return 'translate(' + (index * barWidth) + ', ' +
+        (plotHeight - yScale(value)) + ')';
+  })
+  .attr('height', yScale);
+
+  bar.on('mousemove', function(count, index) {
+    var threshold = counts.min + (index * counts.delta);
+    if (raster.get('threshold') !== threshold) {
+      raster.set('threshold', threshold);
+      raster.changed();
+    }
+  });
+
+  bar.on('mouseover', function(count, index) {
+    var area = 0;
+    for (var i = counts.values.length - 1; i >= index; --i) {
+      area += resolution * resolution * counts.values[i];
+    }
+    tip.html(message(counts.min + (index * counts.delta), area));
+    tip.style('display', 'block');
+    tip.transition().style({
+      left: (chartRect.left + (index * barWidth) + (barWidth / 2)) + 'px',
+      top: (d3.event.y - 60) + 'px',
+      opacity: 1
+    });
+  });
+
+  bar.on('mouseout', function() {
+    tip.transition().style('opacity', 0).each('end', function() {
+      tip.style('display', 'none');
+    });
+  });
+
+}
+
+function message(value, area) {
+  var acres = (area / 4046.86).toFixed(0).replace(/\B(?=(\d{3})+(?!\d))/g, ',');
+  return acres + ' acres at<br>' + value.toFixed(2) + ' VGI or above';
+}

examples/region-growing.css

@@ -0,0 +1,4 @@
+table.controls td {
+  min-width: 110px;
+  padding: 2px 5px;
+}

examples/region-growing.html

@@ -0,0 +1,38 @@
+---
+template: example.html
+title: Region Growing
+shortdesc: Grow a region from a seed pixel
+docs: >
+  <p>Click a region on the map. The computed region will be red.</p>
+  <p>
+    This example uses a <code>ol.source.Raster</code> to generate data
+    based on another source.  The raster source accepts any number of
+    input sources (tile or image based) and runs a pipeline of
+    operations on the input data.  The return from the final
+    operation is used as the data for the output source.
+  </p>
+  <p>
+    In this case, a single tiled source of imagery data is used as input.
+    The region is calculated in a single "image" operation using the "seed"
+    pixel provided by the user clicking on the map. The "threshold" value
+    determines whether a given contiguous pixel belongs to the "region" - the
+    difference between a candidate pixel's RGB values and the seed values must
+    be below the threshold.
+  </p>
+  <p>
+    This example also shows how an additional function can be made available
+    to the operation.
+  </p>
+tags: "raster, region growing"
+---
+<div class="row-fluid">
+  <div class="span12">
+    <div id="map" class="map" style="cursor: pointer"></div>
+    <table class="controls">
+      <tr>
+        <td>Threshold: <span id="threshold-value"></span></td>
+        <td><input id="threshold" type="range" min="1" max="50" value="20"></td>
+      </tr>
+    </table>
+  </div>
+</div>

examples/region-growing.js

@@ -0,0 +1,132 @@
+// NOCOMPILE
+goog.require('ol.Map');
+goog.require('ol.View');
+goog.require('ol.layer.Image');
+goog.require('ol.layer.Tile');
+goog.require('ol.proj');
+goog.require('ol.source.BingMaps');
+goog.require('ol.source.Raster');
+
+function growRegion(inputs, data) {
+  var image = inputs[0];
+  var seed = data.pixel;
+  var delta = parseInt(data.delta);
+  if (!seed) {
+    return image;
+  }
+
+  seed = seed.map(Math.round);
+  var width = image.width;
+  var height = image.height;
+  var inputData = image.data;
+  var outputData = new Uint8ClampedArray(inputData);
+  var seedIdx = (seed[1] * width + seed[0]) * 4;
+  var seedR = inputData[seedIdx];
+  var seedG = inputData[seedIdx + 1];
+  var seedB = inputData[seedIdx + 2];
+  var edge = [seed];
+  while (edge.length) {
+    var newedge = [];
+    for (var i = 0, ii = edge.length; i < ii; i++) {
+      // As noted in the Raster source constructor, this function is provided
+      // using the `lib` option. Other functions will NOT be visible unless
+      // provided using the `lib` option.
+      var next = nextEdges(edge[i]);
+      for (var j = 0, jj = next.length; j < jj; j++) {
+        var s = next[j][0], t = next[j][1];
+        if (s >= 0 && s < width && t >= 0 && t < height) {
+          var ci = (t * width + s) * 4;
+          var cr = inputData[ci];
+          var cg = inputData[ci + 1];
+          var cb = inputData[ci + 2];
+          var ca = inputData[ci + 3];
+          // if alpha is zero, carry on
+          if (ca === 0) {
+            continue;
+          }
+          if (Math.abs(seedR - cr) < delta && Math.abs(seedG - cg) < delta &&
+              Math.abs(seedB - cb) < delta) {
+            outputData[ci] = 255;
+            outputData[ci + 1] = 0;
+            outputData[ci + 2] = 0;
+            outputData[ci + 3] = 255;
+            newedge.push([s, t]);
+          }
+          // mark as visited
+          inputData[ci + 3] = 0;
+        }
+      }
+    }
+    edge = newedge;
+  }
+  return new ImageData(outputData, width, height);
+}
+
+function next4Edges(edge) {
+  var x = edge[0], y = edge[1];
+  return [
+    [x + 1, y],
+    [x - 1, y],
+    [x, y + 1],
+    [x, y - 1]
+  ];
+}
+
+var key = 'Ak-dzM4wZjSqTlzveKz5u0d4IQ4bRzVI309GxmkgSVr1ewS6iPSrOvOKhA-CJlm3';
+
+var imagery = new ol.layer.Tile({
+  source: new ol.source.BingMaps({key: key, imagerySet: 'Aerial'})
+});
+
+var raster = new ol.source.Raster({
+  sources: [imagery.getSource()],
+  operationType: 'image',
+  operation: growRegion,
+  // Functions in the `lib` object will be available to the operation run in
+  // the web worker.
+  lib: {
+    nextEdges: next4Edges
+  }
+});
+
+var rasterImage = new ol.layer.Image({
+  opacity: 0.7,
+  source: raster
+});
+
+var map = new ol.Map({
+  layers: [imagery, rasterImage],
+  target: 'map',
+  view: new ol.View({
+    center: ol.proj.fromLonLat([-119.07, 47.65]),
+    zoom: 11
+  })
+});
+
+var coordinate;
+
+map.on('click', function(event) {
+  coordinate = event.coordinate;
+  raster.changed();
+});
+
+raster.on('beforeoperations', function(event) {
+  // the event.data object will be passed to operations
+  var data = event.data;
+  data.delta = thresholdControl.value;
+  if (coordinate) {
+    data.pixel = map.getPixelFromCoordinate(coordinate);
+  }
+});
+
+var thresholdControl = document.getElementById('threshold');
+
+function updateControlValue() {
+  document.getElementById('threshold-value').innerText = thresholdControl.value;
+}
+updateControlValue();
+
+thresholdControl.addEventListener('input', function() {
+  updateControlValue();
+  raster.changed();
+});

examples/shaded-relief.css

@@ -0,0 +1,4 @@
+table.controls td {
+  text-align: center;
+  padding: 2px 5px;
+}

examples/shaded-relief.html

@@ -0,0 +1,46 @@
+---
+template: example.html
+title: Shaded Relief
+shortdesc: Calculate shaded relief from elevation data
+docs: >
+  <p>
+    This example uses a <code>ol.source.Raster</code> to generate data
+    based on another source.  The raster source accepts any number of
+    input sources (tile or image based) and runs a pipeline of
+    operations on the input data.  The return from the final
+    operation is used as the data for the output source.
+  </p>
+  <p>
+    In this case, a single tiled source of elevation data is used as input.
+    The shaded relief is calculated in a single "image" operation.  By setting
+    <code>operationType: 'image'</code> on the raster source, operations are
+    called with an <code>ImageData</code> object for each of the input sources.
+    Operations are also called with a general purpose <code>data</code> object.
+    In this example, the sun elevation and azimuth data from the inputs above
+    are assigned to this <code>data</code> object and accessed in the shading
+    operation.  The shading operation returns an array of <code>ImageData</code>
+    objects.  When the raster source is used by an image layer, the first
+    <code>ImageData</code> object returned by the last operation in the pipeline
+    is used for rendering.
+  </p>
+tags: "raster, shaded relief"
+---
+<div class="row-fluid">
+  <div class="span12">
+    <div id="map" class="map"></div>
+    <table class="controls">
+      <tr>
+        <td>vertical exaggeration: <span id="vertOut"></span>x</td>
+        <td><input id="vert" type="range" min="1" max="5" value="1"/></td>
+      </tr>
+      <tr>
+        <td>sun elevation: <span id="sunElOut"></span>°</td>
+        <td><input id="sunEl" type="range" min="0" max="90" value="45"/></td>
+      </tr>
+      <tr>
+        <td>sun azimuth: <span id="sunAzOut"></span>°</td>
+        <td><input id="sunAz" type="range" min="0" max="360" value="45"/></td>
+      </tr>
+    </table>
+  </div>
+</div>

examples/shaded-relief.js

@@ -0,0 +1,158 @@
+// NOCOMPILE
+goog.require('ol.Map');
+goog.require('ol.View');
+goog.require('ol.layer.Image');
+goog.require('ol.layer.Tile');
+goog.require('ol.source.Raster');
+goog.require('ol.source.TileJSON');
+goog.require('ol.source.XYZ');
+
+
+/**
+ * Generates a shaded relief image given elevation data.  Uses a 3x3
+ * neighborhood for determining slope and aspect.
+ * @param {Array.<ImageData>} inputs Array of input images.
+ * @param {Object} data Data added in the "beforeoperations" event.
+ * @return {Array.<ImageData>} Output images (only the first is rendered).
+ */
+function shade(inputs, data) {
+  var elevationImage = inputs[0];
+  var width = elevationImage.width;
+  var height = elevationImage.height;
+  var elevationData = elevationImage.data;
+  var shadeData = new Uint8ClampedArray(elevationData.length);
+  var dp = data.resolution * 2;
+  var maxX = width - 1;
+  var maxY = height - 1;
+  var pixel = [0, 0, 0, 0];
+  var twoPi = 2 * Math.PI;
+  var halfPi = Math.PI / 2;
+  var sunEl = Math.PI * data.sunEl / 180;
+  var sunAz = Math.PI * data.sunAz / 180;
+  var cosSunEl = Math.cos(sunEl);
+  var sinSunEl = Math.sin(sunEl);
+  var pixelX, pixelY, x0, x1, y0, y1, offset,
+      z0, z1, dzdx, dzdy, slope, aspect, cosIncidence, scaled;
+  for (pixelY = 0; pixelY <= maxY; ++pixelY) {
+    y0 = pixelY === 0 ? 0 : pixelY - 1;
+    y1 = pixelY === maxY ? maxY : pixelY + 1;
+    for (pixelX = 0; pixelX <= maxX; ++pixelX) {
+      x0 = pixelX === 0 ? 0 : pixelX - 1;
+      x1 = pixelX === maxX ? maxX : pixelX + 1;
+
+      // determine elevation for (x0, pixelY)
+      offset = (pixelY * width + x0) * 4;
+      pixel[0] = elevationData[offset];
+      pixel[1] = elevationData[offset + 1];
+      pixel[2] = elevationData[offset + 2];
+      pixel[3] = elevationData[offset + 3];
+      z0 = data.vert * (pixel[0] + pixel[1] * 2 + pixel[2] * 3);
+
+      // determine elevation for (x1, pixelY)
+      offset = (pixelY * width + x1) * 4;
+      pixel[0] = elevationData[offset];
+      pixel[1] = elevationData[offset + 1];
+      pixel[2] = elevationData[offset + 2];
+      pixel[3] = elevationData[offset + 3];
+      z1 = data.vert * (pixel[0] + pixel[1] * 2 + pixel[2] * 3);
+
+      dzdx = (z1 - z0) / dp;
+
+      // determine elevation for (pixelX, y0)
+      offset = (y0 * width + pixelX) * 4;
+      pixel[0] = elevationData[offset];
+      pixel[1] = elevationData[offset + 1];
+      pixel[2] = elevationData[offset + 2];
+      pixel[3] = elevationData[offset + 3];
+      z0 = data.vert * (pixel[0] + pixel[1] * 2 + pixel[2] * 3);
+
+      // determine elevation for (pixelX, y1)
+      offset = (y1 * width + pixelX) * 4;
+      pixel[0] = elevationData[offset];
+      pixel[1] = elevationData[offset + 1];
+      pixel[2] = elevationData[offset + 2];
+      pixel[3] = elevationData[offset + 3];
+      z1 = data.vert * (pixel[0] + pixel[1] * 2 + pixel[2] * 3);
+
+      dzdy = (z1 - z0) / dp;
+
+      slope = Math.atan(Math.sqrt(dzdx * dzdx + dzdy * dzdy));
+
+      aspect = Math.atan2(dzdy, -dzdx);
+      if (aspect < 0) {
+        aspect = halfPi - aspect;
+      } else if (aspect > halfPi) {
+        aspect = twoPi - aspect + halfPi;
+      } else {
+        aspect = halfPi - aspect;
+      }
+
+      cosIncidence = sinSunEl * Math.cos(slope) +
+          cosSunEl * Math.sin(slope) * Math.cos(sunAz - aspect);
+
+      offset = (pixelY * width + pixelX) * 4;
+      scaled = 255 * cosIncidence;
+      shadeData[offset] = scaled;
+      shadeData[offset + 1] = scaled;
+      shadeData[offset + 2] = scaled;
+      shadeData[offset + 3] = elevationData[offset + 3];
+    }
+  }
+
+  return new ImageData(shadeData, width, height);
+}
+
+var elevation = new ol.source.XYZ({
+  url: 'https://{a-d}.tiles.mapbox.com/v3/aj.sf-dem/{z}/{x}/{y}.png',
+  crossOrigin: 'anonymous'
+});
+
+var raster = new ol.source.Raster({
+  sources: [elevation],
+  operationType: 'image',
+  operation: shade
+});
+
+var map = new ol.Map({
+  target: 'map',
+  layers: [
+    new ol.layer.Tile({
+      source: new ol.source.TileJSON({
+        url: 'http://api.tiles.mapbox.com/v3/tschaub.miapgppd.jsonp'
+      })
+    }),
+    new ol.layer.Image({
+      opacity: 0.3,
+      source: raster
+    })
+  ],
+  view: new ol.View({
+    extent: [-13675026, 4439648, -13580856, 4580292],
+    center: [-13615645, 4497969],
+    minZoom: 10,
+    maxZoom: 16,
+    zoom: 13
+  })
+});
+
+var controlIds = ['vert', 'sunEl', 'sunAz'];
+var controls = {};
+controlIds.forEach(function(id) {
+  var control = document.getElementById(id);
+  var output = document.getElementById(id + 'Out');
+  control.addEventListener('input', function() {
+    output.innerText = control.value;
+    raster.changed();
+  });
+  output.innerText = control.value;
+  controls[id] = control;
+});
+
+raster.on('beforeoperations', function(event) {
+  // the event.data object will be passed to operations
+  var data = event.data;
+  data.resolution = event.resolution;
+  for (var id in controls) {
+    data[id] = Number(controls[id].value);
+  }
+});