算法实现/几何/凸包/单调链
安德鲁的单调链凸包算法在 时间内构建了一组二维点的 凸包。
它首先按字典顺序对点进行排序（首先按x坐标，如果相等，则按y坐标），然后在 时间内构建点的上包络和下包络。
上包络是凸包的一部分，从上面可以看到。它从最右边的点到最左边的点按逆时针顺序运行。下包络是凸包的剩余部分。
Input: a list P of points in the plane.
Precondition: There must be at least 3 points.
Sort the points of P by x-coordinate (in case of a tie, sort by y-coordinate).
Initialize U and L as empty lists.
The lists will hold the vertices of upper and lower hulls respectively.
for i = 1, 2, ..., n:
while L contains at least two points and the sequence of last two points
of L and the point P[i] does not make a counter-clockwise turn:
remove the last point from L
append P[i] to L
for i = n, n-1, ..., 1:
while U contains at least two points and the sequence of last two points
of U and the point P[i] does not make a counter-clockwise turn:
remove the last point from U
append P[i] to U
Remove the last point of each list (it's the same as the first point of the other list).
Concatenate L and U to obtain the convex hull of P.
Points in the result will be listed in counter-clockwise order.
请注意，此实现适用于排序的输入点。该算法基于数组。
FUNCTION Cross(v1,v2,v3)
! IMPLICIT VARIABLE HANDLING
IMPLICIT NONE
!-----------------------------------------------
! INPUT VARIABLES
REAL,INTENT(IN) :: v1(2) !< input vector 1
REAL,INTENT(IN) :: v2(2) !< input vector 2
REAL,INTENT(IN) :: v3(2) !< input vector 3
!-----------------------------------------------
! OUTPUT VARIABLES
REAL :: Cross !< cross product
!-----------------------------------------------
! LOCAL VARIABLES
!===============================================
Cross=(v2(1)-v1(1))*(v3(2)-v1(2))-(v2(2)-v1(2))*(v3(1)-v1(1))
END FUNCTION Cross
SUBROUTINE ConvHull(nPoints,Points,nHull,Hull)
! IMPLICIT VARIABLE HANDLING
IMPLICIT NONE
!------------------------------------------------
! INPUT VARIABLES
INTEGER,INTENT(IN) :: nPoints
REAL,INTENT(IN) :: Points(2,0:nPoints-1)
!------------------------------------------------
! OUTPUT VARIABLES
INTEGER,INTENT(OUT) :: nHull
! NOTE: allocate Hull always one point greater than Points, because we save the first value twice
REAL,INTENT(OUT) :: Hull(2,0:nPoints)
!------------------------------------------------
! LOCAL VARIABLES
REAL :: Lower(2,0:nPoints-1)
REAL :: Upper(2,0:nPoints-1)
INTEGER :: i,iLower,iUpper
!================================================
IF(nPoints.LE.1)THEN
Hull = Points
nHull = nPoints
ELSE
iLower = 0
Lower = -HUGE(1.)
DO i=0,nPoints-1
DO WHILE(iLower.GE.2.AND.Cross(Lower(:,iLower-2),Lower(:,iLower-1),Points(:,i)).LE.0.)
Lower(:,iLower) = -HUGE(1.)
iLower = iLower - 1
END DO
Lower(:,iLower) = Points(:,i)
iLower = iLower + 1
END DO
iUpper = 0
Upper = HUGE(1.)
DO i=nPoints-1,0,-1
DO WHILE(iUpper.GE.2.AND.Cross(Upper(:,iUpper-2),Upper(:,iUpper-1),Points(:,i)).LE.0.)
Upper(:,iUpper) = HUGE(1.)
iUpper = iUpper - 1
END DO
Upper(:,iUpper) = Points(:,i)
iUpper = iUpper + 1
END DO
iLower = iLower-1
iUpper = iUpper-1
nHull = iLower+iUpper+1
! NOTE: Initialize Hull with zeros
Hull = 0.
! NOTE: save values in Hull
Hull(:,0 :iLower -1) = Lower(:,0:iLower-1)
Hull(:,iLower:iLower+iUpper-1) = Upper(:,0:iUpper-1)
! NOTE: save first value twice
Hull(:, iLower+iUpper ) = Hull (:,0 )
END IF
END SUBROUTINE
function cross(a, b, o) {
return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])
}
/**
* @param points An array of [X, Y] coordinates
*/
function convexHull(points) {
points.sort(function(a, b) {
return a[0] == b[0] ? a[1] - b[1] : a[0] - b[0];
});
var lower = [];
for (var i = 0; i < points.length; i++) {
while (lower.length >= 2 && cross(lower[lower.length - 2], lower[lower.length - 1], points[i]) <= 0) {
lower.pop();
}
lower.push(points[i]);
}
var upper = [];
for (var i = points.length - 1; i >= 0; i--) {
while (upper.length >= 2 && cross(upper[upper.length - 2], upper[upper.length - 1], points[i]) <= 0) {
upper.pop();
}
upper.push(points[i]);
}
upper.pop();
lower.pop();
return lower.concat(upper);
}
import java.io.BufferedReader;
import java.io.FileReader;
import java.io.IOException;
import java.util.Arrays;
import java.util.StringTokenizer;
class Point implements Comparable<Point> {
int x, y;
public int compareTo(Point p) {
if (this.x == p.x) {
return this.y - p.y;
} else {
return this.x - p.x;
}
}
public String toString() {
return "(" + x + "," + y + ")";
}
}
public class ConvexHull {
public static long cross(Point O, Point A, Point B) {
return (A.x - O.x) * (long) (B.y - O.y) - (A.y - O.y) * (long) (B.x - O.x);
}
public static Point[] convex_hull(Point[] P) {
if (P.length > 1) {
int n = P.length, k = 0;
Point[] H = new Point[2 * n];
Arrays.sort(P);
// Build lower hull
for (int i = 0; i < n; ++i) {
while (k >= 2 && cross(H[k - 2], H[k - 1], P[i]) <= 0)
k--;
H[k++] = P[i];
}
// Build upper hull
for (int i = n - 2, t = k + 1; i >= 0; i--) {
while (k >= t && cross(H[k - 2], H[k - 1], P[i]) <= 0)
k--;
H[k++] = P[i];
}
if (k > 1) {
H = Arrays.copyOfRange(H, 0, k - 1); // remove non-hull vertices after k; remove k - 1 which is a duplicate
}
return H;
} else if (P.length <= 1) {
return P;
} else {
return null;
}
}
public static void main(String[] args) throws IOException {
BufferedReader f = new BufferedReader(new FileReader("hull.in")); // "hull.in" Input Sample => size x y x y x y x y
StringTokenizer st = new StringTokenizer(f.readLine());
Point[] p = new Point[Integer.parseInt(st.nextToken())];
for (int i = 0; i < p.length; i++) {
p[i] = new Point();
p[i].x = Integer.parseInt(st.nextToken()); // Read X coordinate
p[i].y = Integer.parseInt(st.nextToken()); // Read y coordinate
}
Point[] hull = convex_hull(p).clone();
for (int i = 0; i < hull.length; i++) {
if (hull[i] != null)
System.out.print(hull[i]);
}
}
}
def convex_hull(points):
"""Computes the convex hull of a set of 2D points.
Input: an iterable sequence of (x, y) pairs representing the points.
Output: a list of vertices of the convex hull in counter-clockwise order,
starting from the vertex with the lexicographically smallest coordinates.
Implements Andrew's monotone chain algorithm. O(n log n) complexity.
"""
# Sort the points lexicographically (tuples are compared lexicographically).
# Remove duplicates to detect the case we have just one unique point.
points = sorted(set(points))
# Boring case: no points or a single point, possibly repeated multiple times.
if len(points) <= 1:
return points
# 2D cross product of OA and OB vectors, i.e. z-component of their 3D cross product.
# Returns a positive value, if OAB makes a counter-clockwise turn,
# negative for clockwise turn, and zero if the points are collinear.
def cross(o, a, b):
return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])
# Build lower hull
lower = []
for p in points:
while len(lower) >= 2 and cross(lower[-2], lower[-1], p) <= 0:
lower.pop()
lower.append(p)
# Build upper hull
upper = []
for p in reversed(points):
while len(upper) >= 2 and cross(upper[-2], upper[-1], p) <= 0:
upper.pop()
upper.append(p)
# Concatenation of the lower and upper hulls gives the convex hull.
# Last point of each list is omitted because it is repeated at the beginning of the other list.
return lower[:-1] + upper[:-1]
# Example: convex hull of a 10-by-10 grid.
assert convex_hull([(i//10, i%10) for i in range(100)]) == [(0, 0), (9, 0), (9, 9), (0, 9)]
# the python code converted to ruby syntax
def convex_hull(points)
points.sort!.uniq!
return points if points.length <= 3
def cross(o, a, b)
(a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])
end
lower = Array.new
points.each{|p|
while lower.length > 1 and cross(lower[-2], lower[-1], p) <= 0 do lower.pop end
lower.push(p)
}
upper = Array.new
points.reverse_each{|p|
while upper.length > 1 and cross(upper[-2], upper[-1], p) <= 0 do upper.pop end
upper.push(p)
}
return lower[0...-1] + upper[0...-1]
end
fail unless convex_hull((0..9).to_a.repeated_permutation(2).to_a) == [[0, 0], [9, 0], [9, 9], [0, 9]]
// Implementation of Andrew's monotone chain 2D convex hull algorithm.
// Asymptotic complexity: O(n log n).
// Practical performance: 0.5-1.0 seconds for n=1000000 on a 1GHz machine.
#include <algorithm>
#include <vector>
using namespace std;
typedef double coord_t; // coordinate type
typedef double coord2_t; // must be big enough to hold 2*max(|coordinate|)^2
struct Point {
coord_t x, y;
bool operator <(const Point &p) const {
return x < p.x || (x == p.x && y < p.y);
}
};
// 3D cross product of OA and OB vectors, (i.e z-component of their "2D" cross product, but remember that it is not defined in "2D").
// Returns a positive value, if OAB makes a counter-clockwise turn,
// negative for clockwise turn, and zero if the points are collinear.
coord2_t cross(const Point &O, const Point &A, const Point &B)
{
return (A.x - O.x) * (B.y - O.y) - (A.y - O.y) * (B.x - O.x);
}
// Returns a list of points on the convex hull in counter-clockwise order.
// Note: the last point in the returned list is the same as the first one.
vector<Point> convex_hull(vector<Point> P)
{
size_t n = P.size(), k = 0;
if (n <= 3) return P;
vector<Point> H(2*n);
// Sort points lexicographically
sort(P.begin(), P.end());
// Build lower hull
for (size_t i = 0; i < n; ++i) {
while (k >= 2 && cross(H[k-2], H[k-1], P[i]) <= 0) k--;
H[k++] = P[i];
}
// Build upper hull
for (size_t i = n-1, t = k+1; i > 0; --i) {
while (k >= t && cross(H[k-2], H[k-1], P[i-1]) <= 0) k--;
H[k++] = P[i-1];
}
H.resize(k-1);
return H;
}
use Sort::Key::Radix qw(nkeysort); # use radix sort for an O(n) algorithm
sub convex_hull {
return @_ if @_ < 2;
my @p = nkeysort { $_->[0] } @_;
my (@u, @l);
my $i = 0;
while ($i < @p) {
my $iu = my $il = $i;
my ($x, $yu) = @{$p[$i]};
my $yl = $yu;
# search for upper and lower Y for the current X
while (++$i < @p and $p[$i][0] == $x) {
my $y = $p[$i][1];
if ($y < $yl) {
$il = $i;
$yl = $y;
}
elsif ($y > $yu) {
$iu = $i;
$yu = $y;
}
}
while (@l >= 2) {
my ($ox, $oy) = @{$l[-2]};
last if ($l[-1][1] - $oy) * ($x - $ox) < ($yl - $oy) * ($l[-1][0] - $ox);
pop @l;
}
push @l, $p[$il];
while (@u >= 2) {
my ($ox, $oy) = @{$u[-2]};
last if ($u[-1][1] - $oy) * ($x - $ox) > ($yu - $oy) * ($u[-1][0] - $ox);
pop @u;
}
push @u, $p[$iu];
}
# remove points from the upper hull extremes when they are already
# on the lower hull:
shift @u if $u[0][1] == $l[0][1];
pop @u if @u and $u[-1][1] == $l[-1][1];
return (@l, reverse @u);
}
C 代码来自 Math::ConvexHull::MonotoneChain Perl 模块。请注意，此实现适用于排序的输入点。否则，它与上面的 C++ 实现非常相似。
typedef struct {
double x;
double y;
} point_t;
typedef point_t* point_ptr_t;
/* Three points are a counter-clockwise turn if ccw > 0, clockwise if
* ccw < 0, and collinear if ccw = 0 because ccw is a determinant that
* gives the signed area of the triangle formed by p1, p2 and p3.
*/
static double
ccw(point_t* p1, point_t* p2, point_t* p3)
{
return (p2->x - p1->x)*(p3->y - p1->y) - (p2->y - p1->y)*(p3->x - p1->x);
}
/* Returns a list of points on the convex hull in counter-clockwise order.
* Note: the last point in the returned list is the same as the first one.
*/
void
convex_hull(point_t* points, ssize_t npoints, point_ptr_t** out_hull, ssize_t* out_hullsize)
{
point_ptr_t* hull;
ssize_t i, t, k = 0;
hull = *out_hull;
/* lower hull */
for (i = 0; i < npoints; ++i) {
while (k >= 2 && ccw(hull[k-2], hull[k-1], &points[i]) <= 0) --k;
hull[k++] = &points[i];
}
/* upper hull */
for (i = npoints-2, t = k+1; i >= 0; --i) {
while (k >= t && ccw(hull[k-2], hull[k-1], &points[i]) <= 0) --k;
hull[k++] = &points[i];
}
*out_hull = hull;
*out_hullsize = k;
}
import Data.List (sort)
-- Coordinate type
type R = Double
-- Vector / point type
type R2 = (R, R)
-- Checks if it's shortest to rotate from the OA to the OB vector in a clockwise
-- direction.
clockwise :: R2 -> R2 -> R2 -> Bool
clockwise o a b = (a `sub` o) `cross` (b `sub` o) <= 0
-- 2D cross product.
cross :: R2 -> R2 -> R
cross (x1, y1) (x2, y2) = x1 * y2 - x2 * y1
-- Subtract two vectors.
sub :: R2 -> R2 -> R2
sub (x1, y1) (x2, y2) = (x1 - x2, y1 - y2)
-- Implements the monotone chain algorithm
convexHull :: [R2] -> [R2]
convexHull [] = []
convexHull [p] = [p]
convexHull points = lower ++ upper
where
sorted = sort points
lower = chain sorted
upper = chain (reverse sorted)
chain :: [R2] -> [R2]
chain = go []
where
-- The first parameter accumulates a monotone chain where the most recently
-- added element is at the front of the list.
go :: [R2] -> [R2] -> [R2]
go acc@(r1:r2:rs) (x:xs) =
if clockwise r2 r1 x
-- Made a clockwise turn - remove the most recent part of the chain.
then go (r2:rs) (x:xs)
-- Made a counter-clockwise turn - append to the chain.
else go (x:acc) xs
-- If there's only one point in the chain, just add the next visited point.
go acc (x:xs) = go (x:acc) xs
-- No more points to consume - finished! Note: the reverse here causes the
-- result to be consistent with the other examples (a ccw hull), but
-- removing that and using (upper ++ lower) above will make it cw.
go acc [] = reverse $ tail acc
main :: IO ()
main =
if result == expected
then putStrLn "convexHull worked!"
else fail $ "convexHull broken, got " ++ show result
where
result = convexHull [(x, y) | x <- [0..9], y <- [0..9]]
expected = [(0, 0), (9, 0), (9, 9), (0, 9)]
defmodule ConvexHull do
@type coordinate :: {number, number}
@type point_list :: [coordinate]
@spec find(point_list) :: point_list
def find(points) do
sorted = points |> Enum.sort
left = sorted |> do_half_calc
right = sorted |> Enum.reverse |> do_half_calc
[left, right] |> Enum.concat
end
@spec do_half_calc(point_list) :: point_list
defp do_half_calc(points) do
points
|> Enum.reduce([], &add_to_convex_list/2)
|> tl
|> Enum.reverse
end
@spec perp_prod(coordinate, coordinate, coordinate) :: number
defp perp_prod({x0, y0}, {x1, y1}, {x2, y2}) do
(x1 - x0) * (y2 - y0) - (y1 - y0) * (x2 - x0)
end
@spec add_to_convex_list(coordinate, point_list) :: point_list
defp add_to_convex_list(p2, list) do
{:ok, new_tail} =
list
|> stream_tails
|> Stream.drop_while(fn tail -> oa_left_of_ob(p2, tail) end)
|> Enum.fetch(0)
[p2 | new_tail]
end
@spec oa_left_of_ob(coordinate, point_list) :: number
defp oa_left_of_ob(b, [a, o | _]), do: perp_prod(o, a, b) <= 0
defp oa_left_of_ob(_, _), do: false
@spec stream_tails(point_list) :: Enum.t(point_list)
defp stream_tails(list) do
tails = Stream.unfold(list, fn [_ | t] -> {t, t} end)
Stream.concat([list], tails)
end
end
ExUnit.start
defmodule ConvexHullStart do
use ExUnit.Case, async: true
@points for x <- 0..9, y <- 0..9, do: {x, y}
test "finding the convex hull for all points (0,0) .. (9,9)" do
assert ConvexHull.find(@points) == [{0, 0}, {9, 0}, {9, 9}, {0, 9}]
end
end
(defn convex-hull
"The function convex-hull accepts a collection of vectors containing pairs of numbers
which define the [x y] bounds of the polygon,
e.g. [[0 0] [0 1] [0 2] [0 3]
[1 0] [1 1] [1 2] [1 3]
[2 0] [2 1] [2 2] [2 3] [2 4] [2 5] [2 6] [2 7] [2 8] [2 9]
[3 0] [3 1] [3 2] [3 3]
[4 0] [4 1] [4 2] [4 3] [4 4]
[5 1]]
When run on the above collection of points the output will be
[[0 0] [0 3] [2 9] [4 4] [5 1] [4 0]]
The input collection of vectors can be either a vector ( e.g. [[0 0] [1 1] [2 2]]) or a quoted sequence
(e.g. '([0 0] [1 1] [2 2]) or (quote ([0 0] [1 1] [2 2])).
The output collection will always be a vector of vectors (e.g. [[0 0] [0 9] [9 9] [9 0]]) which
define the bounds of the convex hull."
[c]
(let [cc (sort-by identity c)
x #(first %) ; returns the 'x' coordinate from an [x y] vector
y #(second %) ; returns the 'y' coordinate from an [x y] vector
ccw #(- (* (- (x %2) (x %1)) (- (y %3) (y %1))) (* (- (y %2) (y %1)) (- (x %3) (x %1))))
half-hull #(loop [h [] ; returns the upper half-hull of the collection of vectors
c %]
(if (not (empty? c))
(if (and (> (count h) 1) (<= 0 (ccw (h (- (count h) 2)) (h (- (count h) 1)) (first c))))
(recur (vec (butlast h)) c)
(recur (conj h (first c)) (rest c)))
h))
upper-hull (butlast (half-hull cc))
lower-hull (butlast (half-hull (reverse cc)))]
(vec (concat upper-hull lower-hull ))))
import scala.collection.mutable.ArrayBuffer
def convexHull(points: ArrayBuffer[(Int, Int)]): ArrayBuffer[(Int, Int)] = {
// 2D cross product of OA and OB vectors, i.e. z-component of their 3D cross product.
// Returns a positive value, if OAB makes a counter-clockwise turn,
// negative for clockwise turn, and zero if the points are collinear.
def cross(o: (Int, Int), a: (Int, Int), b: (Int, Int)): Int = {
(a._1 - o._1) * (b._2 - o._2) - (a._2 - o._2) * (b._1 - o._1)
}
val distinctPoints = points.distinct
// No sorting needed if there are less than 2 unique points.
if(distinctPoints.length < 2) {
return points
} else {
val sortedPoints = distinctPoints.sorted
// Build the lower hull
val lower = ArrayBuffer[(Int, Int)]()
for(i <- sortedPoints){
while(lower.length >= 2 && cross(lower(lower.length - 2), lower(lower.length - 1) , i) <= 0){
lower -= lower.last
}
lower += i
}
// Build the upper hull
val upper = ArrayBuffer[(Int, Int)]()
for(i <- sortedPoints.reverse){
while(upper.size >= 2 && cross(upper(upper.length - 2), upper(upper.length - 1) , i) <= 0){
upper -= upper.last
}
upper += i
}
// Last point of each list is omitted because it is repeated at the beginning of the other list.
lower -= lower.last
upper -= upper.last
// Concatenation of the lower and upper hulls gives the convex hull
return lower ++= upper
}
}
$coords = [
[
'x' => <point_0_x>,
'y' => <point_0_y>,
],
[
'x' => <point_1_x>,
'y' => <point_1_y>,
],
[
'x' => <point_n_x>,
'y' => <point_n_y>,
],
];
function orientation(array $coord1, array $coord2, array $coord3)
{
return ($coord2['x'] - $coord1['x']) * ($coord3['y'] - $coord1['y']) - ($coord2['y'] - $coord1['y']) * ($coord3['x'] - $coord1['x']);
}
usort($coords, function ($a, $b) { return $a['x'] == $b['x'] ? $a['y'] - $b['y'] : $a['x'] - $b['x']; });
$upper_hulls = $lower_hulls = [];
foreach ($coords as $coord) {
while (count($last_two = array_slice($lower_hulls, -2)) >= 2 && orientation($last_two[0], $last_two[1], $coord) <= 0) {
array_pop($lower_hulls);
}
$lower_hulls[] = $coord;
}
foreach (array_reverse($coords) as $coord) {
while (count($last_two = array_slice($upper_hulls, -2)) >= 2 && orientation($last_two[0], $last_two[1], $coord) <= 0) {
array_pop($upper_hulls);
}
$upper_hulls[] = $coord;
}
array_pop($upper_hulls);
array_pop($lower_hulls);
return array_merge($upper_hulls, $lower_hulls);
import Foundation
public func cross(_ o: CGPoint, _ a: CGPoint, _ b: CGPoint) -> CGFloat {
let lhs = (a.x - o.x) * (b.y - o.y)
let rhs = (a.y - o.y) * (b.x - o.x)
return lhs - rhs
}
/// Calculate and return the convex hull of a given sequence of points.
///
/// - Remark: Implements Andrew’s monotone chain convex hull algorithm.
///
/// - Complexity: O(*n* log *n*), where *n* is the count of `points`.
///
/// - Parameter points: A sequence of `CGPoint` elements.
///
/// - Returns: An array containing the convex hull of `points`, ordered
/// lexicographically from the smallest coordinates to the largest,
/// turning counterclockwise.
///
public func convexHull<Source>(_ points: Source) -> [CGPoint]
where Source : Collection,
Source.Element == CGPoint
{
// Exit early if there aren’t enough points to work with.
guard points.count > 1 else { return Array(points) }
// Create storage for the lower and upper hulls.
var lower = [CGPoint]()
var upper = [CGPoint]()
// Sort points in lexicographical order.
let points = points.sorted { a, b in
a.x < b.x || a.x == b.x && a.y < b.y
}
// Construct the lower hull.
for point in points {
while lower.count >= 2 {
let a = lower[lower.count - 2]
let b = lower[lower.count - 1]
if cross(a, b, point) > 0 { break }
lower.removeLast()
}
lower.append(point)
}
// Construct the upper hull.
for point in points.lazy.reversed() {
while upper.count >= 2 {
let a = upper[upper.count - 2]
let b = upper[upper.count - 1]
if cross(a, b, point) > 0 { break }
upper.removeLast()
}
upper.append(point)
}
// Remove each array’s last point, as it’s the same as the first point
// in the opposite array, respectively.
lower.removeLast()
upper.removeLast()
// Join the arrays to form the convex hull.
return lower + upper
}
#lang racket
; Monotone Chain method
; function calculates convex hull from
; unsorted association list of points P
(define (convex-hull P)
;
; given a set of points P, this function will
; recursively build an association list of points H
; that define half of the convex hull
(define (half-hull P H passed-H)
; if there are no further points to add
; report the half-hull except for the last point
(if (zero? (length P))
(cdr H)
; if H is empty (the initial state)
; add the first two points from P
(if (zero? (length H))
(let ([H~ (list (cadr P) (car P))])
(half-hull (cddr P) H~ (list)))
; if we have run out of line segments in H to test
; include the candidate point and
; move on to the next point in P
(if (= 1 (length H))
(let ([H~ (append (cons (car P) passed-H) H)])
(half-hull (cdr P) H~ (list)))
; otherwise test the next line segment in H (p1 p2)
; against the candidate point p~
(let* (
[p1 (cadr H)]
[p2 (car H)]
[p~ (car P)]
[x1 (x p1)]
[y1 (y p1)]
[x2 (x p2)]
[y2 (y p2)]
[x~ (x p~)]
[y~ (y p~)])
; if it passes, move on to the next segment
; if it fails, delete the line segment
(if
(> (* (- y2 y1) (- x~ x1)) (* (- x2 x1) (- y~ y1)))
(half-hull P (cdr H) (append passed-H (list p2)))
(half-hull P (cdr H) passed-H)))))))
;
; body of convex-hull
(if
(= 1 (length P))
P
(let ([P~ (sort (sort P > #:key cdr) > #:key car)])
(append
(half-hull P~ (list) (list))
(half-hull (reverse P~) (list) (list))))))
function convex_hull(points)
# Implements Andrew's monotone chain algorithm
# Input: points - vector of tuples (x,y)
# Ouput: the subset of points that define the convex hull
# not enough points
if length(points) <= 1
return copy(points)
end
# sort the points by x and then by y
points = sort(points)
# function for calculating the cross product of vectors OA and OB
cross(o, a, b) = (a[1] - o[1]) * (b[2] - o[2]) - (a[2] - o[2]) * (b[1] - o[1])
# build lower hull
lower = eltype(points)[]
for p in points
while length(lower) >= 2 && cross(lower[end-1], lower[end], p) <= 0
pop!(lower)
end
push!(lower,p)
end
# build upper hull
upper = eltype(points)[]
for i in length(points):-1:1
p = points[i]
while length(upper) >= 2 && cross(upper[end-1], upper[end], p) <= 0
pop!(upper)
end
push!(upper,p)
end
# concatenates lower hull to upper hull to obtain the convex hull
vcat(lower[1:end-1], upper[1:end-1])
end
@assert convex_hull(map(x -> (div(x,10), mod(x,10)), 0:99)) == [(0, 0), (9, 0), (9, 9), (0, 9)]
program MonotoneChain;
uses crt;
type TPoint = record
x,y:Real;
end;
TFunc = function (a,b:TPoint):integer;
TArray = array of TPoint;
function signedAreaOfParallelogram(O,A,B:TPoint):Real;
begin
signedAreaOfParallelogram := (A.x - O.x)*(B.y - O.y) - (A.y - O.y)*(B.x - O.x);
end;
function monotoneChain(P:TArray;n:integer;var H:TArray):integer;
var k,i,t:integer;
begin
if n > 3 then
begin
k := 0;
for i := 0 to n - 1 do
begin
while (k >= 2)and(signedAreaOfParallelogram(H[k-2],H[k-1],P[i]) <= 0)do
k := k - 1;
H[k] := P[i];
k := k + 1;
end;
t := k + 1;
for i := n - 1 downto 1 do
begin
while(k >= t)and(signedAreaOfParallelogram(H[k-2],H[k-1],P[i-1]) <= 0)do
k := k - 1;
H[k] := P[i-1];
k := k + 1;
end;
monotoneChain := k - 1;
end
else
monotoneChain := n;
end;
function comparePoints(a,b:TPoint):integer;
var t:integer;
begin
t := 1;
if(a.x < b.x)or((a.x = b.x)and(a.y < b.y))then
t := -1;
if(a.x = b.x)and(a.y = b.y)then
t := 0;
comparePoints := t;
end;
procedure quicksort(var A:TArray;l,r:integer;cmp:TFunc);
var i,j:integer;
x,w:TPoint;
begin
i := l;
j := r;
x := A[(l+r)div 2];
repeat
while cmp(A[i],x) < 0 do i := i + 1;
while cmp(x,A[j]) < 0 do j := j - 1;
if i <= j then
begin
w := A[i];
A[i] := A[j];
A[j] := w;
i := i + 1;
j := j - 1;
end
until i > j;
if l < j then quicksort(A,l,j,cmp);
if i < r then quicksort(A,i,r,cmp)
end;
var esc:char;
k,m,n:integer;
P,H:TArray;
BEGIN
clrscr;
repeat
repeat
writeln('How many points you want to read');
readln(n);
until n >= 0;
SetLength(P,n);
SetLength(H,n);
for k := 0 to n-1 do
begin
write('P[',k,']=');
readln(P[k].x,P[k].y);
end;
writeln;
if n > 0 then
quicksort(P,0,n-1,@comparePoints);
m := monotoneChain(P,n,H);
writeln('Array of sorted points');
for k := 0 to n-1 do
write('(',P[k].x:1:6,',',P[k].y:1:6,') ');
writeln;
writeln;
writeln('Points on the hull');
for k := 0 to m-1 do
write('(',H[k].x:1:6,',',H[k].y:1:6,') ');
writeln;
writeln;
esc := readkey;
until esc = #27;
END.
- 德·贝格，范·克雷维尔德，奥弗马斯，施瓦茨科夫。计算几何：算法与应用。第二版，施普林格出版社。 ISBN 3540656200.
- 丹·桑德。 二维点集或多边形的凸包。（webcite 存档副本）
- A. M. 安德鲁，“二维凸包的另一种高效算法”，信息处理快报 9, 216-219 (1979)。