I needed some PHP code that finds the area of the minimum bounding rectangle (MBR) of a polygon. I came across a lot of pieces of useful info scattered across the web, but none of these pieces provided enough of a solution. If you are looking for a good implementation that finds the MBR in JAVA then check out Bart Kiers' code. I started porting Bart's implementation to PHP, and although it is a great implementation, it is very tightly coupled with the nuissances of JAVA and several abstractions that Bart implemented. So, I found Godfried Toussaint's 1983 whitepaper on "Solving geometric problems with the rotating calipers" and the Wikipedia page for "Rotating Calipers" and used them to code a simple and brief implementation in PHP. I hope this implementation proves useful in its current form or as a good template for those wishing to create a similar implementation in the language of their choice.
My code can be found on GitHub under the project name BbsRotatingCalipers.
Another Page in the Brain Book
Monday, November 28, 2011
Sunday, November 27, 2011
Finding the Minimum Bounding Width in PHP Using the Rotating Calipers Algorithm
I had a difficult time finding a simple implementation of the Rotating Calipers algorithm, which has various applications like finding the minimum bounding width of a convex polygon. So, I decided to create my own implementation in PHP. My hope is that this post will help others to create simple implementations in the languages that they desire.
Note: This implementation is only for a convex polygon, but you can adapt it to work with all polygons by first finding the Convex Hull using Westhoff's implementation of the QuickHull algorithm. You can call the following function to calculate the convex hull before calling minBoundingWidth().
And here are some samples of how it is used:
/*
* Copyright (c) 2011, Geoffrey Cox
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY
* OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT
* LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE
* FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR
* OTHERWISE, ARISING FROM, OUT OF OR IN
* CONNECTION WITH THE SOFTWARE OR THE
* USE OR OTHER DEALINGS IN THE SOFTWARE.
*
* Developed by: Geoffrey Cox, http://brainbooksoftware.com
*
*/
class Point
{
public $x;
public $y;
public function __construct($x=0, $y=0)
{
$this->x = $x;
$this->y = $y;
}
public function __toString()
{
return "(" . $this->x . "," . $this->y . ")";
}
}
/**
* Used to make the code more readable
*/
class Vector extends Point {}
/**
* Provides a kind of wrapping so that array indexes are easier to use
*/
function getItem(&$a, $i)
{
$sz = sizeof($a);
if ($i >= $sz)
$i = $i % $sz;
return $a[$i];
}
/**
* Returns the angle in radians between vector v1 and vector v2
*/
function angle($v1, $v2)
{
$n = $v1->x*$v2->x + $v1->y*$v2->y;
$d = sqrt($v1->x*$v1->x + $v1->y*$v1->y)
*sqrt($v2->x*$v2->x + $v2->y*$v2->y);
return acos($n/$d);
}
/**
* Rotates the vector v r radians
*/
function rotate($v, $r)
{
$v2 = new Vector();
$v2->x = $v->x*cos($r) - $v->y*sin($r);
$v2->y = $v->x*sin($r) + $v->y*cos($r);
return $v2;
}
/**
* Shortest distance from point p to the line formed by extending the
* vector v through point t
*/
function distance($p, $t, $v)
{
if ($v->x == 0)
return abs($p->x - $t->x);
$a = $v->y/$v->x;
$c = $t->y - $a*$t->x;
return abs($p->y - $c - $a*$p->x)/sqrt($a*$a + 1);
}
/**
* Calculates the minimum bounding width of a convex polygon in
* O(n) time where n is the number of points in the polygon.
* This function uses the "Rotating Calipers" algorithm. For a
* detailed explanation of this algorithm, see: Toussaint, Godfried
* T. (1983). Solving geometric problems with the rotating
* calipers. Proc. MELECON '83, Athens.
* Special thanks to the wikipedia page that contained pseudocode
* for this algorithm and Bart Kiers for his JAVA implementation.
* Precondition: The points must be in a counterclockwise order
* e.g. echo minBoundingWidth(
* array(new Point(1, 1),
* new Point(2, 0),
* new Point(3, 2)));
*/
function minBoundingWidth($points)
{
$pA = 0; // index of vertex with minimum y-coordinate
$pB = 0; // index of vertex with maximum y-coordinate
for ($i = 1; $i < sizeof($points); $i++) {
if ($points[$i]->y < $points[$pA]->y)
$pA = $i;
if ($points[$i]->y > $points[$pB]->y)
$pB = $i;
}
$rotatedAngle = 0;
$minWidth = null;
$caliperA = new Vector(1, 0); // along the positive x-axis
$caliperB = new Vector(-1, 0); // along the negative x-axis
while ($rotatedAngle < M_PI) {
// Determine the angle between each caliper and the next
// adjacent edge in the polygon
$edgeA = new Vector(
getItem($points, $pA + 1)->x - getItem($points, $pA)->x,
getItem($points, $pA + 1)->y - getItem($points, $pA)->y);
$edgeB = new Vector(
getItem($points, $pB + 1)->x - getItem($points, $pB)->x,
getItem($points, $pB + 1)->y - getItem($points, $pB)->y);
$angleA = angle($edgeA, $caliperA);
$angleB = angle($edgeB, $caliperB);
$width = 0;
// Rotate the calipers by the smaller of these angles
$caliperA = rotate($caliperA, min($angleA, $angleB));
$caliperB = rotate($caliperB, min($angleA, $angleB));
if ($angleA < $angleB) {
$width = distance(getItem($points, $pB),
getItem($points, $pA),
$caliperA);
$pA++;
} else {
$width = distance(getItem($points, $pA),
getItem($points, $pB),
$caliperB);
$pB++;
}
$rotatedAngle = $rotatedAngle + min($angleA, $angleB);
if ($minWidth === null || $width < $minWidth)
$minWidth = $width;
}
return $minWidth;
}
echo minBoundingWidth(array(new Point(0, 0), new Point(1, 0),
new Point(1, 1), new Point(0, 1))) . "\n"; // echos 1
echo minBoundingWidth(array(new Point(1, 1), new Point(2, 0),
new Point(3, 2))) . "\n"; // echos 1.3416407864999
echo minBoundingWidth(array(new Point(0, 0), new Point(2, 0),
new Point(0, 2))) . "\n"; // echos 1.4142135623731
Note: This implementation is only for a convex polygon, but you can adapt it to work with all polygons by first finding the Convex Hull using Westhoff's implementation of the QuickHull algorithm. You can call the following function to calculate the convex hull before calling minBoundingWidth().
require_once(dirname(__FILE__) . "/convex_hull.php");
/**
* A wrapper that gets the convex hull in our desired format using
* Westhoff's function
*/
function convexHullPoints($points)
{
$whPoints = array();
foreach ($points as $p) {
array_push($whPoints, array($p->x, $p->y));
}
$hull = new ConvexHull($whPoints);
$hullPoints = $hull->getHullPoints();
$newPoints = array();
for ($i = sizeof($hullPoints) - 1; $i >= 0; $i--) {
array_push($newPoints, new Point($hullPoints[$i][0],
$hullPoints[$i][1]));
}
return $newPoints;
}
Thursday, November 24, 2011
PHP Function that Calculates Distance Between a Line and a Point
Here is a brief PHP function that calculates the distance between a line and a point. I had trouble finding a good example online so I decided to create my own:
/**
* Copyright (c) 2011, Geoffrey Cox
* You are granted unlimited rights to this code.
* Finds the distance from a point (px, py) to a line that
* passes through points ($x1, $y1) and ($x2, $y2).
* Based on logic found at
* http://www.worsleyschool.net/science/files/...
* linepoint/method5.html
*/
function distanceFromLine($px, $py, $x1, $y1, $x2, $y2)
{
$ad = $x2 - $x1;
if ($ad == 0)
return abs($px - $x1);
$a = ($y2 - $y1)/$ad;
$c = $y1 - $a*$x1;
return abs($py - $c - $a*$px)/sqrt($a*$a + 1);
}
Subscribe to:
Posts (Atom)