[back]

var util = require('util');

Point = function(x,y) {
  return {
    x: x,
    y: y,
    normalLeft: function() {
      return Point(this.y, -this.x);
    },
    center: function(other_point) {
      return Point((this.x + other_point.x)/2.0, (this.y + other_point.y)/2.0);
    },
    distance: function(other_point) {
      return Math.pow(Math.pow(this.x - other_point.x, 2.0) + Math.pow(this.y + other_point.y, 2.0), 0.5);
    },
    polar: function() {
      var theta = 0.0;
      if (closeTo(this.x, 0.0)) {
        if (this.y > 0) {
          theta = 90.0;
        } else if (this.y < 0) {
          theta = 270.0;
        }
      } else {
        theta = rad2deg(Math.atan(this.y/this.x));
        if (this.x < 0.0) {
          theta += 180.0;
        }
        if (theta < 0.0) {
          theta += 360.0;
        }
      }
      return theta;
    },
    angleBetween: function(p1, p2) {
      vect_p1 = difference(p1,this);
      vect_p2 = difference(p2,this);
      theta = vect_p1.polar() - vect_p2.polar();
      if (theta < 0.0) {
        theta += 360.0;
      }
      if (theta > 180.0) {
        theta = 360.0 - theta;
      }
      return theta;
    },
    equalTo: function(other_point) {
      return closeTo(this.x, other_point.x) && closeTo(this.y, other_point.y);
    },
    to_s: function() {
      return "(" + this.x + "," + this.y + ")";
    },
    compareTo: function(other_point) {
      //console.log("Comparing self " + this.to_s() + " to " + other_point.to_s());
      if (closeTo(this.x, other_point.x)) {
        if (closeTo(this.y, other_point.y)) {
          return 0;
        }
        return (this.y - other_point.y);
      } else {
        //console.log("x coords are different, returning " + (this.x - other_point.x));
        return (this.x - other_point.x);
      }
    },
    scale: function(factor) {
      return Point(this.x * factor, this.y * factor);
    }
  };
}

function sum(p1, p2) {
  return Point(p1.x + p2.x, p1.y + p2.y);
}

function difference(p1, p2) {
  return Point(p1.x - p2.x, p1.y - p2.y);
}

function randomPoint() {
  var radius = Math.random();
  var theta = 360 * Math.random();
  return Point(radius * Math.cos(theta), radius * Math.sin(theta));
}

function closeTo(x, y) {
  return Math.abs(x-y) < 1e-7;
}

function rad2deg(a) {
  return (a * 180.0) / Math.PI;
}

function infLineIntersection(base1, v1, base2, v2) {
  if (closeTo(v1.x, 0.0) && closeTo(v2.y, 0.0)) {
    return Point(base1.x, base2.y);
  }
  if (closeTo(v1.y, 0.0) && closeTo(v2.x, 0.0)) {
    return Point(base2.x, base1.y);
  }

  var m1 = !closeTo(v1.x, 0.0) ? v1.y / v1.x : 0.0
  var m2 = !closeTo(v2.x, 0.0) ? v2.y / v2.x : 0.0

  if (closeTo(m1, m2)) {
    return null;
  }

  var c1 = base1.y - m1 * base1.x;
  var c2 = base2.y = m2 * base2.x;
  if (closeTo(v1.x, 0.0)) {
    return Point(base1.x, base1.x * m2 + c2);
  }
  if (closeTo(v2.x, 0.0)) {
    return Point(base2.x, base2.x * m1 + c1);
  }
  var ix = (c2 - c1) / (m1-m2);
  var iy = m1 * ix + c1;
  return Point(ix,iy);
}

Circle = function() {
  return {
    center: null,
    radius: null,
    contains: function(point) {
      return point.distance(this.center) < this.radius || closeTo(point.distance(this.center) - this.radius, 0.0);
    }
  };
}

Solver = function() {
  return {
    mec: Circle(),
    points: [],
    solve: function(a) {
      a.sort(function(p1,p2) { return p1.compareTo(p2); });
      for (var i in a) {
        this.addPoint(a[i]);
      }
      return this.mec;
    },
    direction: function(p0, p1, p2) {
      return (p0.x-p1.x)*(p2.y-p1.y) - (p2.x-p1.x)*(p0.y-p1.y);
    },
    halfHull: function(left, right, input, factor) {
      input = input.slice(0);
      input.push(right);
      var half = [];
      half.push(left);
      for (var i in input) {
        half.push(input[i]);
        //console.log("HALF: " + util.inspect(half) + " LENGTH: " + half.length);
        while (half.length >= 3) {
          var len = half.length;
          var dir = factor * this.direction(half[len-3], half[len-1], half[len-2]);
          if (dir <= 0) {
            half.splice(len-2,1)
          } else {
            break;
          }
        }
      }
      return half;
    },
    convexHull: function() {
      var left = this.points[0];
      var a = this.points.slice(1);
      var right = a.pop();
      var upper = [];
      var lower = [];
      for (var i in a) {
        var p = a[i];
        var dir = this.direction(left, right, p);
        if (dir < 0) {
          upper.push(p);
        } else {
          lower.push(p);
        }
      }
      var upper_hull = this.halfHull(left, right, upper, -1);
      var lower_hull = this.halfHull(left, right, lower, 1);
      var combined = upper_hull.concat(lower_hull);
      var ret = [];
      for (var i in combined) {
        var e = combined[i];
        var has = false;
        for (var j in ret) {
          if (ret[j].equalTo(e)) {
            has = true;
            break;
          }
        }
        if (!has) {
          ret.push(e);
        }
      }
      return ret;
    },
    addPoint: function(point) {
      this.points.push(point);
      if (this.mec.center) {
        if (!this.mec.contains(point)) {
          this.iterate();
        }
      } else {
        this.iterate();
      }
    },
    iterate: function() {
      if (this.points.length == 1) {
        this.mec.center = this.points[0];
        this.mec.radius = 0;
      } else if (this.points.length == 2) {
        this.mec.center = this.points[0].center(this.points[1]);
        this.mec.radius = this.points[0].distance(this.mec.center);
      } else {
        var a = this.convexHull();
        var point_a = a[0];
        var point_c = a[1];
        while (true) {
          var point_b = null;
          var best_theta = 180.0;
          for (var i in this.points) {
            var point = this.points[i];
            if (point != point_a && point != point_c) {
              var theta_abc = point.angleBetween(point_a, point_c);
              if (theta_abc < best_theta) {
                point_b = point;
                best_theta = theta_abc;
              }
            }
          }
          if (best_theta >= 90.0 || point_b == null) {
            this.mec.center = point_a.center(point_c);
            this.mec.radius = point_a.distance(this.mec.center);
            return this.mec;
          }
          var ang_bca = point_c.angleBetween(point_b, point_a);
          var ang_cab = point_a.angleBetween(point_c, point_b);
          if (ang_bca > 90.0) {
            point_c = point_b;
          } else if (ang_cab <= 90.0) {
            break;
          } else {
            point_a = point_b;
          }
        }
        var ch1 = difference(point_b, point_a).scale(0.5);
        var ch2 = difference(point_c, point_a).scale(0.5);
        var n1 = ch1.normalLeft();
        var n2 = ch2.normalLeft();
        ch1 = sum(point_a, ch1);
        ch2 = sum(point_a, ch2);
        this.mec.center = infLineIntersection(ch1, n1, ch2, n2);
        this.mec.radius = this.mec.center.distance(point_a);
      }
      return this.mec;
    }
  };
}

var s = Solver();
console.log("Solver: " + util.inspect(s));

for (var iter = 0; iter < 5; iter++) {
  var start = new Date().getTime();
  var n = Math.pow(10,iter);
  for (var j = 0; j < 25; j++) {
    var a = [];
    for (var i = 0; i < n; i++) {
      a.push(randomPoint());
    }
    var circle = s.solve(a);
  }
  var stop = new Date().getTime();
  var delta = (stop - start) / 1000.0;
  console.log("Points: " + n + ": " + delta);
}

Doug's Compendium of Stuff

Doug Seifert: A Compendium of Stuff. Random thoughts, ramblings, code and other stuff by an experienced software engineer and architect specializing in ruby, rails and Java web and server side development.