contour_plot.js 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455
  1. (function (global, factory) {
  2. typeof exports === 'object' && typeof module !== 'undefined' ? factory(exports) :
  3. typeof define === 'function' && define.amd ? define('contour_plot', ['exports'], factory) :
  4. factory((global.contour_plot = {}));
  5. }(this, function (exports) { 'use strict';
  6. /** finds the zeros of a function, given two starting points (which must
  7. * have opposite signs */
  8. function bisect(f, a, b, parameters) {
  9. parameters = parameters || {};
  10. var maxIterations = parameters.maxIterations || 100,
  11. tolerance = parameters.tolerance || 1e-10,
  12. fA = f(a),
  13. fB = f(b),
  14. delta = b - a;
  15. if (fA * fB > 0) {
  16. throw "Initial bisect points must have opposite signs";
  17. }
  18. if (fA === 0) return a;
  19. if (fB === 0) return b;
  20. for (var i = 0; i < maxIterations; ++i) {
  21. delta /= 2;
  22. var mid = a + delta,
  23. fMid = f(mid);
  24. if (fMid * fA >= 0) {
  25. a = mid;
  26. }
  27. if ((Math.abs(delta) < tolerance) || (fMid === 0)) {
  28. return mid;
  29. }
  30. }
  31. return a + delta;
  32. }
  33. // This file is modified from the d3.geom.contour
  34. // plugin found here https://github.com/d3/d3-plugins/tree/master/geom/contour
  35. /*
  36. Copyright (c) 2012-2015, Michael Bostock
  37. All rights reserved.
  38. Redistribution and use in source and binary forms, with or without
  39. modification, are permitted provided that the following conditions are met:
  40. * Redistributions of source code must retain the above copyright notice, this
  41. list of conditions and the following disclaimer.
  42. * Redistributions in binary form must reproduce the above copyright notice,
  43. this list of conditions and the following disclaimer in the documentation
  44. and/or other materials provided with the distribution.
  45. * The name Michael Bostock may not be used to endorse or promote products
  46. derived from this software without specific prior written permission.
  47. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  48. AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  49. IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
  50. DISCLAIMED. IN NO EVENT SHALL MICHAEL BOSTOCK BE LIABLE FOR ANY DIRECT,
  51. INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
  52. BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  53. DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
  54. OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  55. NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
  56. EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  57. */
  58. function d3_contour(grid, start) {
  59. var s = start || d3_geom_contourStart(grid), // starting point
  60. c = [], // contour polygon
  61. x = s[0], // current x position
  62. y = s[1], // current y position
  63. dx = 0, // next x direction
  64. dy = 0, // next y direction
  65. pdx = NaN, // previous x direction
  66. pdy = NaN, // previous y direction
  67. i = 0;
  68. do {
  69. // determine marching squares index
  70. i = 0;
  71. if (grid(x-1, y-1)) i += 1;
  72. if (grid(x, y-1)) i += 2;
  73. if (grid(x-1, y )) i += 4;
  74. if (grid(x, y )) i += 8;
  75. // determine next direction
  76. if (i === 6) {
  77. dx = pdy === -1 ? -1 : 1;
  78. dy = 0;
  79. } else if (i === 9) {
  80. dx = 0;
  81. dy = pdx === 1 ? -1 : 1;
  82. } else {
  83. dx = d3_geom_contourDx[i];
  84. dy = d3_geom_contourDy[i];
  85. }
  86. // update contour polygon
  87. if (dx != pdx && dy != pdy) {
  88. c.push([x, y]);
  89. pdx = dx;
  90. pdy = dy;
  91. } else {
  92. c.push([x, y]);
  93. }
  94. x += dx;
  95. y += dy;
  96. } while (s[0] != x || s[1] != y);
  97. return c;
  98. }
  99. // lookup tables for marching directions
  100. var d3_geom_contourDx = [1, 0, 1, 1,-1, 0,-1, 1,0, 0,0,0,-1, 0,-1,NaN];
  101. var d3_geom_contourDy = [0,-1, 0, 0, 0,-1, 0, 0,1,-1,1,1, 0,-1, 0,NaN];
  102. function d3_geom_contourStart(grid) {
  103. var x = 0,
  104. y = 0;
  105. // search for a starting point; begin at origin
  106. // and proceed along outward-expanding diagonals
  107. while (true) {
  108. if (grid(x,y)) {
  109. return [x,y];
  110. }
  111. if (x === 0) {
  112. x = y + 1;
  113. y = 0;
  114. } else {
  115. x = x - 1;
  116. y = y + 1;
  117. }
  118. }
  119. }
  120. function isoline(f, value, xScale, yScale) {
  121. var xRange = xScale.range(), yRange = yScale.range();
  122. return function(x, y) {
  123. if ((x < xRange[0]) || (x > xRange[1]) ||
  124. (y < yRange[0]) || (y > yRange[1])) return false;
  125. return f(xScale.invert(x), yScale.invert(y)) < value;
  126. };
  127. }
  128. function smoothPoints(f, points, level, xScale, yScale) {
  129. var xRange = xScale.range(), yRange = yScale.range();
  130. var ySmooth = function(y) {
  131. return f(xScale.invert(x), yScale.invert(y)) - level;
  132. };
  133. var xSmooth = function(x) {
  134. return f(xScale.invert(x), yScale.invert(y)) - level;
  135. };
  136. for (var k = 0; k < points.length; ++k) {
  137. var point = points[k],
  138. x = point[0], y = point[1];
  139. if ((x <= xRange[0]) || (x >= xRange[1]) ||
  140. (y <= yRange[0]) || (y >= yRange[1])) continue;
  141. var currentSmooth = ySmooth(y);
  142. var p = {'maxIterations' : 9};
  143. for (var delta = 0.5; delta <= 3; delta += 0.5) {
  144. if (ySmooth(y - delta) * currentSmooth < 0) {
  145. y = bisect(ySmooth, y, y - delta, p);
  146. } else if (xSmooth(x - delta) * currentSmooth < 0) {
  147. x = bisect(xSmooth, x, x - delta, p);
  148. } else if (ySmooth(y + delta) * currentSmooth < 0) {
  149. y = bisect(ySmooth, y, y + delta, p);
  150. } else if (xSmooth(x + delta) * currentSmooth < 0) {
  151. x = bisect(xSmooth, x, x + delta, p);
  152. } else {
  153. continue;
  154. }
  155. break;
  156. }
  157. point[0] = x;
  158. point[1] = y;
  159. }
  160. }
  161. function getLogLevels(f, xScale, yScale, count) {
  162. var xRange = xScale.range(), yRange = yScale.range();
  163. // figure out min/max values by sampling pointson a grid
  164. var maxValue, minValue, value;
  165. maxValue = minValue = f(xScale.invert(xRange[0]), yScale.invert(yRange[0]));
  166. for (var y = yRange[0]; y < yRange[1]+1; ++y) {
  167. for (var x = xRange[0]; x < xRange[1]+1; ++x) {
  168. value = f(xScale.invert(x),yScale.invert(y));
  169. minValue = Math.min(value, minValue);
  170. maxValue = Math.max(value, maxValue);
  171. }
  172. }
  173. // lets get contour lines on a log scale, keeping
  174. // values on an integer scale (if possible)
  175. var levels = [];
  176. var logRange = Math.log(maxValue - Math.floor(minValue));
  177. var base = Math.ceil(Math.exp(logRange / (count))),
  178. upper = Math.pow(base, Math.ceil(logRange / Math.log(base)));
  179. for (var i = 0; i < count; ++i) {
  180. var current = Math.floor(minValue) + upper;
  181. if (current < minValue) {
  182. break;
  183. }
  184. levels.push(current);
  185. upper /= base;
  186. }
  187. return levels;
  188. }
  189. function getStartingPoint(lineFunc, x, y) {
  190. x = Math.floor(x);
  191. y = Math.floor(y);
  192. var j = 0;
  193. while (true) {
  194. j += 1;
  195. if (!lineFunc(x+j, y)) {
  196. return [x+j, y];
  197. }
  198. if (!lineFunc(x, y+j)) {
  199. return [x, y+j];
  200. }
  201. }
  202. }
  203. function getContours(f, xScale, yScale, count, minima) {
  204. // figure out even distribution in log space of values
  205. var levels = getLogLevels(f, xScale, yScale, count);
  206. // use marching squares algo from d3.geom.contour to build up a series of paths
  207. var ret = [];
  208. for (var i = 0; i < levels.length; ++i) {
  209. var level = levels[i];
  210. var lineFunc = isoline(f, level, xScale, yScale);
  211. var points= [];
  212. if (minima) {
  213. var initialPoints = [];
  214. for (var m = 0; m < minima.length; ++m) {
  215. var initial = getStartingPoint(lineFunc, xScale(minima[m].x), yScale(minima[m].y));
  216. var current = d3_contour(lineFunc, initial);
  217. // don't add points if already seen
  218. var duplicate = false;
  219. for (var j = 0 ; j < current.length; ++j) {
  220. var point = current[j];
  221. for (var k = 0; k < initialPoints.length; ++k) {
  222. var other = initialPoints[k];
  223. if ((point[0] == other[0]) &&
  224. (point[1] == other[1])) {
  225. duplicate = true;
  226. break;
  227. }
  228. }
  229. if (duplicate) break;
  230. }
  231. if (duplicate) continue;
  232. initialPoints.push(initial);
  233. smoothPoints(f, current, level, xScale, yScale);
  234. if (points.length) points.push(null);
  235. points = points.concat(current);
  236. }
  237. } else {
  238. points = d3_contour(lineFunc);
  239. smoothPoints(f, points, level, xScale, yScale);
  240. }
  241. ret.push(points);
  242. }
  243. // return the contours
  244. return {'paths': ret, 'levels': levels};
  245. }
  246. function ContourPlot() {
  247. var drawAxis = false,
  248. f = function (x, y) { return (1 - x) * (1 - x) + 100 * (y - x * x) * ( y - x * x); },
  249. yDomain = [3, -3],
  250. xDomain = [-2, 2],
  251. minima = null,
  252. contourCount = 14,
  253. colourScale = d3.scaleLinear().domain([0, contourCount]).range(["white", d3.schemeCategory10[0]]);
  254. // todo: resolution independent (sample say 200x200)
  255. // todo: handle function with multiple local minima
  256. function chart(selection) {
  257. var width = selection.nodes()[0].offsetWidth,
  258. height = width * 0.75,
  259. padding = (drawAxis) ? 24 : 0,
  260. yScale = d3.scaleLinear()
  261. .range([padding, height - padding])
  262. .domain(yDomain),
  263. xScale = d3.scaleLinear()
  264. .range([padding, width - padding])
  265. .domain(xDomain);
  266. // create tooltip if doesn't exist
  267. d3.select("body").selectAll(".contour_tooltip").data([0]).enter()
  268. .append("div")
  269. .attr("class", "contour_tooltip")
  270. .style("font-size", "12px")
  271. .style("position", "absolute")
  272. .style("text-align", "center")
  273. .style("width", "128px")
  274. .style("height", "32px")
  275. .style("background", "#333")
  276. .style("color", "#ddd")
  277. .style("padding", "0px")
  278. .style("border", "0px")
  279. .style("border-radius", "8px")
  280. .style("opacity", "0");
  281. var tooltip = d3.selectAll(".contour_tooltip");
  282. // create the svg element if it doesn't already exist
  283. selection.selectAll("svg").data([0]).enter().append("svg");
  284. var svg = selection.selectAll("svg").data([0]);
  285. svg.attr("width", width)
  286. .attr("height", height)
  287. .on("mouseover", function() {
  288. tooltip.transition().duration(400).style("opacity", 0.9);
  289. tooltip.style("z-index", "");
  290. })
  291. .on("mousemove", function() {
  292. var point = d3.mouse(this),
  293. x = xScale.invert(point[0]),
  294. y = yScale.invert(point[1]),
  295. fx = f(x, y);
  296. tooltip.style("left", (d3.event.pageX) + "px")
  297. .style("top", (d3.event.pageY - 44) + "px");
  298. tooltip.html("x = " + x.toFixed(2) + " y = " + y.toFixed(2) + "<br>f(x,y) = " + fx.toFixed(2) );
  299. })
  300. .on("mouseout", function() {
  301. tooltip.transition().duration(400).style("opacity", 0);
  302. tooltip.style("z-index", -1);
  303. });
  304. var contours = getContours(f, xScale, yScale, contourCount, minima);
  305. var paths = contours.paths,
  306. levels = contours.levels;
  307. var line = d3.line()
  308. .x(function(d) { return d[0]; })
  309. .y(function(d) { return d[1]; })
  310. .curve(d3.curveLinearClosed)
  311. .defined(function(d) { return d; });
  312. var pathGroup = svg.append("g");
  313. pathGroup.selectAll("path").data(paths).enter()
  314. .append("path")
  315. .attr("d", line)
  316. .style("fill", function(d, i) { return colourScale(i); })
  317. .style("stroke-width", 1.5)
  318. .style("stroke", "white")
  319. .on("mouseover", function() {
  320. d3.select(this).style("stroke-width", "4");
  321. })
  322. .on("mouseout", function() {
  323. d3.select(this).style("stroke-width", "1.5");
  324. });
  325. // draw axii
  326. if (drawAxis) {
  327. var xAxis = d3.axisBottom().scale(xScale),
  328. yAxis = d3.axisLeft().scale(yScale);
  329. svg.append("g")
  330. .attr("class", "axis")
  331. .attr("transform", "translate(0," + (height - 1.0 * padding) + ")")
  332. .call(xAxis);
  333. svg.append("g")
  334. .attr("class", "axis")
  335. .attr("transform", "translate(" + (padding) + ",0)")
  336. .call(yAxis);
  337. }
  338. return {'xScale' : xScale, 'yScale' : yScale, 'svg' : svg};
  339. }
  340. chart.drawAxis = function(_) {
  341. if (!arguments.length) return drawAxis;
  342. drawAxis = _;
  343. return chart;
  344. };
  345. chart.xDomain = function(_) {
  346. if (!arguments.length) return xDomain;
  347. xDomain = _;
  348. return chart;
  349. };
  350. chart.yDomain = function(_) {
  351. if (!arguments.length) return yDomain;
  352. yDomain = _;
  353. return chart;
  354. };
  355. chart.colourScale = function(_) {
  356. if (!arguments.length) return colourScale;
  357. colourScale = _;
  358. return chart;
  359. };
  360. chart.contourCount = function(_) {
  361. if (!arguments.length) return contourCount;
  362. contourCount = _;
  363. return chart;
  364. };
  365. chart.minima = function(_) {
  366. if (!arguments.length) return minima;
  367. minima = _;
  368. return chart;
  369. };
  370. chart.f = function(_) {
  371. if (!arguments.length) return f;
  372. f = _;
  373. return chart;
  374. };
  375. return chart;
  376. }
  377. var version = "0.0.1";
  378. exports.version = version;
  379. exports.isoline = isoline;
  380. exports.smoothPoints = smoothPoints;
  381. exports.getLogLevels = getLogLevels;
  382. exports.getContours = getContours;
  383. exports.ContourPlot = ContourPlot;
  384. }));