d3-spring-model/src/interpolation.js

export default function(sampleSet, remainderSet, sampleSubset, distanceFn) {
  var distance = distanceFunction;
  // var distance = calculateEuclideanDistance;
  // console.log("Brute-force");
  for (let node of remainderSet) {
    let nearestSample = undefined,
        minDist = undefined,
        sampleSubsetDistanceCache = [];

    for (let sample of sampleSet) {
      let dist = distanceFn(node, sample);
      if (nearestSample === undefined || dist < minDist) {
        minDist = dist;
        nearestSample = sample;
      }

      let index = sampleSubset.indexOf(sample);
      if (index !== -1) {
        sampleSubsetDistanceCache[index] = dist;
      }
    }

    placeNearToNearestNeighbour(node, nearestSample, minDist, sampleSubset, sampleSubsetDistanceCache);
  }
}


function placeNearToNearestNeighbour(node, minNode, radius, samples, realDistances) {
  let
    dist0 = sumDistError(pointOnCircle(minNode.x, minNode.y, 0, radius), samples, realDistances),
    dist90 = sumDistError(pointOnCircle(minNode.x, minNode.y, 90, radius), samples, realDistances),
    dist180 = sumDistError(pointOnCircle(minNode.x, minNode.y, 180, radius), samples, realDistances),
    dist270 = sumDistError(pointOnCircle(minNode.x, minNode.y, 270, radius), samples, realDistances),
    lowBound = 0.0,
    highBound = 0.0;

  // Determine the closest quadrant
  if (dist0 == dist180) {
    if (dist90 > dist270)
      lowBound = highBound = 270;
    else
      lowBound = highBound = 90;
  } else if (dist90 == dist270) {
    if (dist0 > dist180)
      lowBound = highBound = 180;
    else
      lowBound = highBound = 0;
  } else if (dist0 > dist180) {
    if (dist90 > dist270) {
      lowBound = 180;
      highBound = 270;
    } else {
      lowBound = 90;
      highBound = 180;
    }
  } else {
    if (dist90 > dist270) {
      lowBound = 270;
      highBound = 360;
    } else {
      lowBound = 0;
      highBound = 90;
    }
  }

  let angle = binarySearchMin(lowBound, highBound,
                             function(angle){
                               return sumDistError(pointOnCircle(minNode.x, minNode.y, angle, radius), samples, realDistances);
                             });
  let newPoint = pointOnCircle(minNode.x, minNode.y, angle, radius);

  // console.log(newPoint);
  node.x = newPoint.x;
  node.y = newPoint.y;

  // for (var i = 0; i < 20; i++) {
  //   var forces = sumForcesToSample(node, sample, sampleCache);
  //   // console.log(forces);
  //   node.x += forces.x;
  //   node.y += forces.y;
  // }

}

// With a circle radius r, and center at (h,k),
// Find the coordinate of a point at angle degree
function pointOnCircle(h, k, angle, r) {
  let x = h + r*Math.cos(toRadians(angle));
  let y = k + r*Math.sin(toRadians(angle));

  return {
    x: x,
    y: y
  };
}

function toRadians(degrees) {
  return degrees * (Math.PI / 180);
}

function sumDistError(currentPos, samples, realDistances) {
  let total = 0.0;
  for (let i = 0; i < samples.length; i++) {
    let sample = samples[i];
    let lowDDistance = Math.hypot(sample.x - currentPos.x, sample.y - currentPos.y);
    total += Math.abs(lowDDistance - realDistances[i]);
  }
  return total;
}


function sumForcesToSample(node, sample, sampleCache) {
  var x = 0,
    y = 0,
    // len = 0,
    dist = 0,
    force,
    SPRING_FORCE = 0.7;

  for (var i = 0, unitX, unitY; i < sample.length; i++) {
    var s = sample[i];
    if (s !== node) {
      unitX = s.x - node.x;
      unitY = s.y - node.y;

      // Normalize coordinates
      len = Math.sqrt(unitX * unitX + unitY * unitY);
      unitX /= len;
      unitY /= len;

      console.log(unitX, unitY);

      var realDist = Math.sqrt(unitX * unitX + unitY * unitY);
      var desDist = sampleCache[i];
      dist += realDist - desDist;
      force = (SPRING_FORCE * dist);

      x += unitX * force;
      y += unitY * force;
    }

    x *= (1.0 / sample.length);
    y *= (1.0 / sample.length);

    return {
      x: x,
      y: y
    };
  }
}

function binarySearchMin(lb, hb, fn) {
  while (lb <= hb) {
    if(lb === hb) return lb;
    if(hb-lb == 1) {
      if (fn(lb) >= fn(hb)) return hb;
      else return lb;
    }

    let range = hb-lb;
    let valLowerHalf = fn(lb + range/4);
    let valHigherHalf = fn(lb + range*3/4);

    if (valLowerHalf > valHigherHalf)
      lb = Math.floor((lb + hb) / 2);
    else if (valLowerHalf < valHigherHalf)
      hb = Math.ceil((lb + hb) / 2);
    else {
      lb += Math.floor(range/4);
      hb -= Math.ceil(range/4);
    }
  }
  return -1;
}