Doug's Compendium of Stuff

Stuff

Main
Where Is Doug?
My Bookmarks
My Resume
Pacbell Outages
Pilot Logbook

Code

Ruby B-Tree
Tournament
Jar Search
DCT6200 Firewire Tuner (For Motorola DCT6200 set top box and Freevo)
Serial Port Tuner (For Motorola DCT2XXX set top box and Freevo)
Gentoo IVTV Ebuild
Java Internet Jukebox

Old

Offsite

My Wife
My Trainer
My Work

Contact

Email:

[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);
}