'use strict'; var Lib = require('../../lib'); var constants = require('./constants'); module.exports = function findAllPaths(pathinfo, xtol, ytol) { var cnt, startLoc, i, pi, j; // Default just passes these values through as they were before: xtol = xtol || 0.01; ytol = ytol || 0.01; for(i = 0; i < pathinfo.length; i++) { pi = pathinfo[i]; for(j = 0; j < pi.starts.length; j++) { startLoc = pi.starts[j]; makePath(pi, startLoc, 'edge', xtol, ytol); } cnt = 0; while(Object.keys(pi.crossings).length && cnt < 10000) { cnt++; startLoc = Object.keys(pi.crossings)[0].split(',').map(Number); makePath(pi, startLoc, undefined, xtol, ytol); } if(cnt === 10000) Lib.log('Infinite loop in contour?'); } }; function equalPts(pt1, pt2, xtol, ytol) { return Math.abs(pt1[0] - pt2[0]) < xtol && Math.abs(pt1[1] - pt2[1]) < ytol; } // distance in index units - uses the 3rd and 4th items in points function ptDist(pt1, pt2) { var dx = pt1[2] - pt2[2]; var dy = pt1[3] - pt2[3]; return Math.sqrt(dx * dx + dy * dy); } function makePath(pi, loc, edgeflag, xtol, ytol) { var locStr = loc.join(','); var mi = pi.crossings[locStr]; var marchStep = getStartStep(mi, edgeflag, loc); // start by going backward a half step and finding the crossing point var pts = [getInterpPx(pi, loc, [-marchStep[0], -marchStep[1]])]; var m = pi.z.length; var n = pi.z[0].length; var startLoc = loc.slice(); var startStep = marchStep.slice(); var cnt; // now follow the path for(cnt = 0; cnt < 10000; cnt++) { // just to avoid infinite loops if(mi > 20) { mi = constants.CHOOSESADDLE[mi][(marchStep[0] || marchStep[1]) < 0 ? 0 : 1]; pi.crossings[locStr] = constants.SADDLEREMAINDER[mi]; } else { delete pi.crossings[locStr]; } marchStep = constants.NEWDELTA[mi]; if(!marchStep) { Lib.log('Found bad marching index:', mi, loc, pi.level); break; } // find the crossing a half step forward, and then take the full step pts.push(getInterpPx(pi, loc, marchStep)); loc[0] += marchStep[0]; loc[1] += marchStep[1]; locStr = loc.join(','); // don't include the same point multiple times if(equalPts(pts[pts.length - 1], pts[pts.length - 2], xtol, ytol)) pts.pop(); var atEdge = (marchStep[0] && (loc[0] < 0 || loc[0] > n - 2)) || (marchStep[1] && (loc[1] < 0 || loc[1] > m - 2)); var closedLoop = loc[0] === startLoc[0] && loc[1] === startLoc[1] && marchStep[0] === startStep[0] && marchStep[1] === startStep[1]; // have we completed a loop, or reached an edge? if((closedLoop) || (edgeflag && atEdge)) break; mi = pi.crossings[locStr]; } if(cnt === 10000) { Lib.log('Infinite loop in contour?'); } var closedpath = equalPts(pts[0], pts[pts.length - 1], xtol, ytol); var totaldist = 0; var distThresholdFactor = 0.2 * pi.smoothing; var alldists = []; var cropstart = 0; var distgroup, cnt2, cnt3, newpt, ptcnt, ptavg, thisdist, i, j, edgepathi, edgepathj; /* * Check for points that are too close together (<1/5 the average dist * *in grid index units* (important for log axes and nonuniform grids), * less if less smoothed) and just take the center (or avg of center 2). * This cuts down on funny behavior when a point is very close to a * contour level. */ for(cnt = 1; cnt < pts.length; cnt++) { thisdist = ptDist(pts[cnt], pts[cnt - 1]); totaldist += thisdist; alldists.push(thisdist); } var distThreshold = totaldist / alldists.length * distThresholdFactor; function getpt(i) { return pts[i % pts.length]; } for(cnt = pts.length - 2; cnt >= cropstart; cnt--) { distgroup = alldists[cnt]; if(distgroup < distThreshold) { cnt3 = 0; for(cnt2 = cnt - 1; cnt2 >= cropstart; cnt2--) { if(distgroup + alldists[cnt2] < distThreshold) { distgroup += alldists[cnt2]; } else break; } // closed path with close points wrapping around the boundary? if(closedpath && cnt === pts.length - 2) { for(cnt3 = 0; cnt3 < cnt2; cnt3++) { if(distgroup + alldists[cnt3] < distThreshold) { distgroup += alldists[cnt3]; } else break; } } ptcnt = cnt - cnt2 + cnt3 + 1; ptavg = Math.floor((cnt + cnt2 + cnt3 + 2) / 2); // either endpoint included: keep the endpoint if(!closedpath && cnt === pts.length - 2) newpt = pts[pts.length - 1]; else if(!closedpath && cnt2 === -1) newpt = pts[0]; // odd # of points - just take the central one else if(ptcnt % 2) newpt = getpt(ptavg); // even # of pts - average central two else { newpt = [(getpt(ptavg)[0] + getpt(ptavg + 1)[0]) / 2, (getpt(ptavg)[1] + getpt(ptavg + 1)[1]) / 2]; } pts.splice(cnt2 + 1, cnt - cnt2 + 1, newpt); cnt = cnt2 + 1; if(cnt3) cropstart = cnt3; if(closedpath) { if(cnt === pts.length - 2) pts[cnt3] = pts[pts.length - 1]; else if(cnt === 0) pts[pts.length - 1] = pts[0]; } } } pts.splice(0, cropstart); // done with the index parts - remove them so path generation works right // because it depends on only having [xpx, ypx] for(cnt = 0; cnt < pts.length; cnt++) pts[cnt].length = 2; // don't return single-point paths (ie all points were the same // so they got deleted?) if(pts.length < 2) return; else if(closedpath) { pts.pop(); pi.paths.push(pts); } else { if(!edgeflag) { Lib.log('Unclosed interior contour?', pi.level, startLoc.join(','), pts.join('L')); } // edge path - does it start where an existing edge path ends, or vice versa? var merged = false; for(i = 0; i < pi.edgepaths.length; i++) { edgepathi = pi.edgepaths[i]; if(!merged && equalPts(edgepathi[0], pts[pts.length - 1], xtol, ytol)) { pts.pop(); merged = true; // now does it ALSO meet the end of another (or the same) path? var doublemerged = false; for(j = 0; j < pi.edgepaths.length; j++) { edgepathj = pi.edgepaths[j]; if(equalPts(edgepathj[edgepathj.length - 1], pts[0], xtol, ytol)) { doublemerged = true; pts.shift(); pi.edgepaths.splice(i, 1); if(j === i) { // the path is now closed pi.paths.push(pts.concat(edgepathj)); } else { if(j > i) j--; pi.edgepaths[j] = edgepathj.concat(pts, edgepathi); } break; } } if(!doublemerged) { pi.edgepaths[i] = pts.concat(edgepathi); } } } for(i = 0; i < pi.edgepaths.length; i++) { if(merged) break; edgepathi = pi.edgepaths[i]; if(equalPts(edgepathi[edgepathi.length - 1], pts[0], xtol, ytol)) { pts.shift(); pi.edgepaths[i] = edgepathi.concat(pts); merged = true; } } if(!merged) pi.edgepaths.push(pts); } } // special function to get the marching step of the // first point in the path (leading to loc) function getStartStep(mi, edgeflag, loc) { var dx = 0; var dy = 0; if(mi > 20 && edgeflag) { // these saddles start at +/- x if(mi === 208 || mi === 1114) { // if we're starting at the left side, we must be going right dx = loc[0] === 0 ? 1 : -1; } else { // if we're starting at the bottom, we must be going up dy = loc[1] === 0 ? 1 : -1; } } else if(constants.BOTTOMSTART.indexOf(mi) !== -1) dy = 1; else if(constants.LEFTSTART.indexOf(mi) !== -1) dx = 1; else if(constants.TOPSTART.indexOf(mi) !== -1) dy = -1; else dx = -1; return [dx, dy]; } /* * Find the pixel coordinates of a particular crossing * * @param {object} pi: the pathinfo object at this level * @param {array} loc: the grid index [x, y] of the crossing * @param {array} step: the direction [dx, dy] we're moving on the grid * * @return {array} [xpx, ypx, xi, yi]: the first two are the pixel location, * the next two are the interpolated grid indices, which we use for * distance calculations to delete points that are too close together. * This is important when the grid is nonuniform (and most dramatically when * we're on log axes and include invalid (0 or negative) values. * It's crucial to delete these extra two before turning an array of these * points into a path, because those routines require length-2 points. */ function getInterpPx(pi, loc, step) { var locx = loc[0] + Math.max(step[0], 0); var locy = loc[1] + Math.max(step[1], 0); var zxy = pi.z[locy][locx]; var xa = pi.xaxis; var ya = pi.yaxis; // Interpolate in linear space, then convert to pixel if(step[1]) { var dx = (pi.level - zxy) / (pi.z[locy][locx + 1] - zxy); // Interpolate, but protect against NaN linear values for log axis (dx will equal 1 or 0) var dxl = (dx !== 1 ? (1 - dx) * xa.c2l(pi.x[locx]) : 0) + (dx !== 0 ? dx * xa.c2l(pi.x[locx + 1]) : 0); return [xa.c2p(xa.l2c(dxl), true), ya.c2p(pi.y[locy], true), locx + dx, locy]; } else { var dy = (pi.level - zxy) / (pi.z[locy + 1][locx] - zxy); var dyl = (dy !== 1 ? (1 - dy) * ya.c2l(pi.y[locy]) : 0) + (dy !== 0 ? dy * ya.c2l(pi.y[locy + 1]) : 0); return [xa.c2p(pi.x[locx], true), ya.c2p(ya.l2c(dyl), true), locx, locy + dy]; } }