Convex hull

From Rosetta Code
Convex hull is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Find the points which form a convex hull from a set of arbitrary two dimensional points.

For example, given the points (16,3), (12,17), (0,6), (-4,-6), (16,6), (16,-7), (16,-3), (17,-4), (5,19), (19,-8), (3,16), (12,13), (3,-4), (17,5), (-3,15), (-3,-9), (0,11), (-9,-3), (-4,-2) and (12,10) the convex hull would be (-9,-3), (-3,-9), (19,-8), (17,5), (12,17), (5,19) and (-3,15).

See also



Go[edit]

package main
 
import (
"fmt"
"image"
"sort"
)
 
 
// ConvexHull returns the set of points that define the
// convex hull of p in CCW order starting from the left most.
func (p points) ConvexHull() points {
// From https://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain
// with only minor deviations.
sort.Sort(p)
var h points
 
// Lower hull
for _, pt := range p {
for len(h) >= 2 && !ccw(h[len(h)-2], h[len(h)-1], pt) {
h = h[:len(h)-1]
}
h = append(h, pt)
}
 
// Upper hull
for i, t := len(p)-2, len(h)+1; i >= 0; i-- {
pt := p[i]
for len(h) >= t && !ccw(h[len(h)-2], h[len(h)-1], pt) {
h = h[:len(h)-1]
}
h = append(h, pt)
}
 
return h[:len(h)-1]
}
 
// ccw returns true if the three points make a counter-clockwise turn
func ccw(a, b, c image.Point) bool {
return ((b.X - a.X) * (c.Y - a.Y)) > ((b.Y - a.Y) * (c.X - a.X))
}
 
type points []image.Point
 
func (p points) Len() int { return len(p) }
func (p points) Swap(i, j int) { p[i], p[j] = p[j], p[i] }
func (p points) Less(i, j int) bool {
if p[i].X == p[j].X {
return p[i].Y < p[i].Y
}
return p[i].X < p[j].X
}
 
func main() {
pts := points{
{16, 3}, {12, 17}, {0, 6}, {-4, -6}, {16, 6},
{16, -7}, {16, -3}, {17, -4}, {5, 19}, {19, -8},
{3, 16}, {12, 13}, {3, -4}, {17, 5}, {-3, 15},
{-3, -9}, {0, 11}, {-9, -3}, {-4, -2}, {12, 10},
}
hull := pts.ConvexHull()
fmt.Println("Convex Hull:", hull)
}
Output:
Convex Hull: [(-9,-3) (-3,-9) (19,-8) (17,5) (12,17) (5,19) (-3,15)]

Haskell[edit]

import Data.List (sortBy, groupBy, maximumBy)
import Data.Ord (comparing)
 
(x, y) = ((!! 0), (!! 1))
 
compareFrom
:: (Num a, Ord a)
=> [a] -> [a] -> [a] -> Ordering
compareFrom o l r =
compare ((x l - x o) * (y r - y o)) ((y l - y o) * (x r - x o))
 
distanceFrom
:: Floating a
=> [a] -> [a] -> a
distanceFrom from to = ((x to - x from) ** 2 + (y to - y from) ** 2) ** (1 / 2)
 
convexHull
:: (Floating a, Ord a)
=> [[a]] -> [[a]]
convexHull points =
let o = minimum points
presorted = sortBy (compareFrom o) (filter (/= o) points)
collinears = groupBy (((EQ ==) .) . compareFrom o) presorted
outmost = maximumBy (comparing (distanceFrom o)) <$> collinears
in dropConcavities [o] outmost
 
dropConcavities
:: (Num a, Ord a)
=> [[a]] -> [[a]] -> [[a]]
dropConcavities (left:lefter) (right:righter:rightest) =
case compareFrom left right righter of
LT -> dropConcavities (right : left : lefter) (righter : rightest)
EQ -> dropConcavities (left : lefter) (righter : rightest)
GT -> dropConcavities lefter (left : righter : rightest)
dropConcavities output lastInput = lastInput ++ output
 
main :: IO ()
main =
mapM_ print $
convexHull
[ [16, 3]
, [12, 17]
, [0, 6]
, [-4, -6]
, [16, 6]
, [16, -7]
, [16, -3]
, [17, -4]
, [5, 19]
, [19, -8]
, [3, 16]
, [12, 13]
, [3, -4]
, [17, 5]
, [-3, 15]
, [-3, -9]
, [0, 11]
, [-9, -3]
, [-4, -2]
, [12, 10]
]
Output:
[-3.0,-9.0]
[19.0,-8.0]
[17.0,5.0]
[12.0,17.0]
[5.0,19.0]
[-3.0,15.0]
[-9.0,-3.0]

J[edit]

Restated from the implementation at http://kukuruku.co/hub/funcprog/introduction-to-j-programming-language-2004 which in turn is a translation of http://dr-klm.livejournal.com/42312.html

counterclockwise =: ({. , }. /: 12 o. }. - {.) @ /:~
crossproduct =: 11"_ o. [: (* +)/ }. - {.
removeinner =: #~ 1, 0 > 3 crossproduct\ ], 1:
hull =: [: removeinner^:_ counterclockwise

Example use:

   hull 16j3 12j17 0j6 _4j_6 16j6 16j_7 16j_3 17j_4 5j19 19j_8 3j16 12j13 3j_4 17j5 _3j15 _3j_9 0j11 _9j_3 _4j_2 12j10
_9j_3 _3j_9 19j_8 17j5 12j17 5j19 _3j15

Perl 6[edit]

Works with: Rakudo version 2017.05
Translation of: zkl

Modified the angle sort method as the original could fail if there were multiple points on the same y coordinate as the starting point. Sorts on tangent rather than triangle area. Inexpensive since it still doesn't do any trigonometric math, just calculates the ratio of opposite over adjacent. The original returned the correct answer for the task example, but only by accident. If the points (14,-9), (1,-9) were added to the task example, it wouldn't give a correct answer. Now it does.

class Point {
has Real $.x is rw;
has Real $.y is rw;
method gist { [~] '(', self.x,', ', self.y, ')' };
}
 
sub ccw (Point $a, Point $b, Point $c) {
($b.x - $a.x)*($c.y - $a.y) - ($b.y - $a.y)*($c.x - $a.x);
}
 
sub tangent (Point $a, Point $b) {
my $opp = $b.x - $a.x;
my $adj = $b.y - $a.y;
$adj != 0 ?? $opp / $adj !! Inf
}
 
sub graham-scan (**@coords) {
# sort points by y, secondary sort on x
my @sp = @coords.map( { Point.new( :x($_[0]), :y($_[1]) ) } )
.sort: {.y, .x};
 
# need at least 3 points to make a hull
return @sp if +@sp < 3;
 
# first point on hull is minimum y point
my @h = @sp.shift;
 
# re-sort the points by angle, secondary on x
@sp = @sp.map( { $++ => [tangent(@h[0], $_), $_.x] } )
.sort( {-$_.value[0], $_.value[1] } )
.map: { @sp[$_.key] };
 
# first point of re-sorted list is guaranteed to be on hull
@h.push: @sp.shift;
 
# check through the remaining list making sure that
# there is always a positive angle
for @sp -> $point {
if ccw( |@h.tail(2), $point ) >= 0 {
@h.push: $point;
} else {
@h.pop;
redo;
}
}
@h
}
 
my @hull = graham-scan(
(16, 3), (12,17), ( 0, 6), (-4,-6), (16, 6), (16,-7), (16,-3),
(17,-4), ( 5,19), (19,-8), ( 3,16), (12,13), ( 3,-4), (17, 5),
(-3,15), (-3,-9), ( 0,11), (-9,-3), (-4,-2), (12,10)
);
 
say "Convex Hull ([email protected]} points): ", @hull;
 
@hull = graham-scan(
(16, 3), (12,17), ( 0, 6), (-4,-6), (16, 6), (16,-7), (16,-3),
(17,-4), ( 5,19), (19,-8), ( 3,16), (12,13), ( 3,-4), (17, 5),
(-3,15), (-3,-9), ( 0,11), (-9,-3), (-4,-2), (12,10), (14,-9), (1,-9)
);
 
say "Convex Hull ([email protected]} points): ", @hull;
Output:
Convex Hull (7 points): [(-3, -9) (19, -8) (17, 5) (12, 17) (5, 19) (-3, 15) (-9, -3)]
Convex Hull (9 points): [(-3, -9) (1, -9) (14, -9) (19, -8) (17, 5) (12, 17) (5, 19) (-3, 15) (-9, -3)]

Python[edit]

An approach that uses the shapely library:

from __future__ import print_function
from shapely.geometry import MultiPoint
 
if __name__=="__main__":
pts = MultiPoint([(16,3), (12,17), (0,6), (-4,-6), (16,6), (16,-7), (16,-3), (17,-4), (5,19), (19,-8), (3,16), (12,13), (3,-4), (17,5), (-3,15), (-3,-9), (0,11), (-9,-3), (-4,-2), (12,10)])
print (pts.convex_hull)
Output:
POLYGON ((-3 -9, -9 -3, -3 15, 5 19, 12 17, 17 5, 19 -8, -3 -9))

Racket[edit]

Also an implementation of https://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain (therefore kinda
Translation of: Go
#lang typed/racket
(define-type Point (Pair Real Real))
(define-type Points (Listof Point))
 
(: ⊗ (Point Point Point -> Real))
(define/match (⊗ o a b)
[((cons o.x o.y) (cons a.x a.y) (cons b.x b.y))
(- (* (- a.x o.x) (- b.y o.y)) (* (- a.y o.y) (- b.x o.x)))])
 
(: Point<? (Point Point -> Boolean))
(define (Point<? a b)
(cond [(< (car a) (car b)) #t] [(> (car a) (car b)) #f] [else (< (cdr a) (cdr b))]))
 
;; Input: a list P of points in the plane.
(define (convex-hull [P:unsorted : Points])
 ;; Sort the points of P by x-coordinate (in case of a tie, sort by y-coordinate).
(define P (sort P:unsorted Point<?))
 ;; 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
 ;; TB: U is identical with (reverse P)
(define (upper/lower-hull [P : Points])
(reverse
(for/fold ((rev : Points null))
((P.i (in-list P)))
(let u/l : Points ((rev rev))
(match rev
[(list p-2 p-1 ps ...) #:when (not (positive? (⊗ p-2 P.i p-1))) (u/l (list* p-1 ps))]
[(list ps ...) (cons P.i ps)])))))
 
 ;; Initialize U and L as empty lists.
 ;; The lists will hold the vertices of upper and lower hulls respectively.
(let ((U (upper/lower-hull (reverse P)))
(L (upper/lower-hull P)))
 ;; 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.
(append (drop-right L 1) (drop-right U 1)))) ; Points in the result will be listed in CCW order.)
 
(module+ test
(require typed/rackunit)
(check-equal?
(convex-hull
(list '(16 . 3) '(12 . 17) '(0 . 6) '(-4 . -6) '(16 . 6) '(16 . -7) '(16 . -3) '(17 . -4)
'(5 . 19) '(19 . -8) '(3 . 16) '(12 . 13) '(3 . -4) '(17 . 5) '(-3 . 15) '(-3 . -9)
'(0 . 11) '(-9 . -3) '(-4 . -2) '(12 . 10)))
(list '(-9 . -3) '(-3 . -9) '(19 . -8) '(17 . 5) '(12 . 17) '(5 . 19) '(-3 . 15))))
Output:

silence implies tests pass (and output is as expected)

REXX[edit]

version 1[edit]

/* REXX ---------------------------------------------------------------
* Compute the Convex Hull for a set of points
* Format of the input file:
* (16,3) (12,17) (0,6) (-4,-6) (16,6) (16,-7) (16,-3) (17,-4) (5,19)
* (19,-8) (3,16) (12,13) (3,-4) (17,5) (-3,15) (-3,-9) (0,11) (-9,-3)
* (-4,-2)
*--------------------------------------------------------------------*/

Signal On Novalue
Signal On Syntax
Parse Arg fid
If fid='' Then Do
fid='chullmin.in' /* miscellaneous test data */
fid='chullx.in'
fid='chullt.in'
fid='chulla.in'
fid='chullxx.in'
fid='sq.in'
fid='tri.in'
fid='line.in'
fid='point.in'
fid='chull.in' /* data from task description */
End
g.0debug=''
g.0oid=fn(fid)'.txt'; 'erase' g.0oid
x.=0
yl.=''
Parse Value '1000 -1000' With g.0xmin g.0xmax
Parse Value '1000 -1000' With g.0ymin g.0ymax
/*---------------------------------------------------------------------
* First read the input and store the points' coordinates
* x.0 contains the number of points, x.i contains the x.coordinate
* yl.x contains the y.coordinate(s) of points (x/y)
*--------------------------------------------------------------------*/

Do while lines(fid)>0
l=linein(fid)
Do While l<>''
Parse Var l '(' x ',' y ')' l
Call store x,y
End
End
Call lineout fid
Do i=1 To x.0 /* loop over points */
x=x.i
yl.x=sortv(yl.x) /* sort y-coordinates */
End
Call sho
 
/*---------------------------------------------------------------------
* Now we look for special border points:
* lefthigh and leftlow: leftmost points with higheste and lowest y
* ritehigh and ritelow: rightmost points with higheste and lowest y
* yl.x contains the y.coordinate(s) of points (x/y)
*--------------------------------------------------------------------*/

leftlow=0
lefthigh=0
Do i=1 To x.0
x=x.i
If maxv(yl.x)=g.0ymax Then Do
If lefthigh=0 Then lefthigh=x'/'g.0ymax
ritehigh=x'/'g.0ymax
End
If minv(yl.x)=g.0ymin Then Do
ritelow=x'/'g.0ymin
If leftlow=0 Then leftlow=x'/'g.0ymin
End
End
Call o 'lefthigh='lefthigh
Call o 'ritehigh='ritehigh
Call o 'ritelow ='ritelow
Call o 'leftlow ='leftlow
/*---------------------------------------------------------------------
* Now we look for special border points:
* leftmost_n and leftmost_s: points with lowest x and highest/lowest y
* ritemost_n and ritemost_s: points with largest x and highest/lowest y
* n and s stand foNorth and South, respectively
*--------------------------------------------------------------------*/

x=g.0xmi; leftmost_n=x'/'maxv(yl.x)
x=g.0xmi; leftmost_s=x'/'minv(yl.x)
x=g.0xma; ritemost_n=x'/'maxv(yl.x)
x=g.0xma; ritemost_s=x'/'minv(yl.x)
 
/*---------------------------------------------------------------------
* Now we compute the paths from ritehigh to ritelow (n_end)
* and leftlow to lefthigh (s_end), respectively
*--------------------------------------------------------------------*/

x=g.0xma
n_end=''
Do i=words(yl.x) To 1 By -1
n_end=n_end x'/'word(yl.x,i)
End
Call o 'n_end='n_end
x=g.0xmi
s_end=''
Do i=1 To words(yl.x)
s_end=s_end x'/'word(yl.x,i)
End
Call o 's_end='s_end
 
n_high=''
s_low=''
/*---------------------------------------------------------------------
* Now we compute the upper part of the convex hull (nhull)
*--------------------------------------------------------------------*/

Call o 'leftmost_n='leftmost_n
Call o 'lefthigh ='lefthigh
nhull=leftmost_n
res=mk_nhull(leftmost_n,lefthigh);
nhull=nhull res
Call o 'A nhull='nhull
Do While res<>lefthigh
res=mk_nhull(res,lefthigh); nhull=nhull res
Call o 'B nhull='nhull
End
res=mk_nhull(lefthigh,ritemost_n); nhull=nhull res
Call o 'C nhull='nhull
Do While res<>ritemost_n
res=mk_nhull(res,ritemost_n); nhull=nhull res
Call o 'D nhull='nhull
End
 
nhull=nhull n_end /* attach the right vertical border */
 
/*---------------------------------------------------------------------
* Now we compute the lower part of the convex hull (shull)
*--------------------------------------------------------------------*/

res=mk_shull(ritemost_s,ritelow);
shull=ritemost_s res
Call o 'A shull='shull
Do While res<>ritelow
res=mk_shull(res,ritelow)
shull=shull res
Call o 'B shull='shull
End
res=mk_shull(ritelow,leftmost_s)
shull=shull res
Call o 'C shull='shull
Do While res<>leftmost_s
res=mk_shull(res,leftmost_s);
shull=shull res
Call o 'D shull='shull
End
 
shull=shull s_end
 
chull=nhull shull /* concatenate upper and lower part */
/* eliminate duplicates */
/* too lazy to take care before :-) */
Parse Var chull chullx chull
have.=0
have.chullx=1
Do i=1 By 1 While chull>''
Parse Var chull xy chull
If have.xy=0 Then Do
chullx=chullx xy
have.xy=1
End
End
/* show the result */
Say 'Points of convex hull in clockwise order:'
Say chullx
/**********************************************************************
* steps that were necessary in previous attempts
/*---------------------------------------------------------------------
* Final polish: Insert points that are not yet in chullx but should be
* First on the upper hull going from left to right
*--------------------------------------------------------------------*/

i=1
Do While i<words(chullx)
xya=word(chullx,i)  ; Parse Var xya xa '/' ya
If xa=g.0xmax Then Leave
xyb=word(chullx,i+1); Parse Var xyb xb '/' yb
Do j=1 To x.0
If x.j>xa Then Do
If x.j<xb Then Do
xx=x.j
parse Value kdx(xya,xyb) With k d x
If (k*xx+d)=maxv(yl.xx) Then Do
chullx=subword(chullx,1,i) xx'/'maxv(yl.xx),
subword(chullx,i+1)
i=i+1
End
End
End
Else
i=i+1
End
End
 
Say chullx
 
/*---------------------------------------------------------------------
* Final polish: Insert points that are not yet in chullx but should be
* Then on the lower hull going from right to left
*--------------------------------------------------------------------*/

i=wordpos(ritemost_s,chullx)
Do While i<words(chullx)
xya=word(chullx,i)  ; Parse Var xya xa '/' ya
If xa=g.0xmin Then Leave
xyb=word(chullx,i+1); Parse Var xyb xb '/' yb
Do j=x.0 To 1 By -1
If x.j<xa Then Do
If x.j>xb Then Do
xx=x.j
parse Value kdx(xya,xyb) With k d x
If (k*xx+d)=minv(yl.xx) Then Do
chullx=subword(chullx,1,i) xx'/'minv(yl.xx),
subword(chullx,i+1)
i=i+1
End
End
End
Else
i=i+1
End
End
Say chullx
**********************************************************************/
Call lineout g.0oid
 
Exit
 
store: Procedure Expose x. yl. g.
/*---------------------------------------------------------------------
* arrange the points in ascending order of x (in x.) and,
* for each x in ascending order of y (in yl.x)
* g.0xmin is the smallest x-value, etc.
* g.0xmi is the x-coordinate
* g.0ymin is the smallest y-value, etc.
* g.0ymi is the x-coordinate of such a point
*--------------------------------------------------------------------*/

Parse Arg x,y
Call o 'store' x y
If x<g.0xmin Then Do; g.0xmin=x; g.0xmi=x; End
If x>g.0xmax Then Do; g.0xmax=x; g.0xma=x; End
If y<g.0ymin Then Do; g.0ymin=y; g.0ymi=x; End
If y>g.0ymax Then Do; g.0ymax=y; g.0yma=x; End
Do i=1 To x.0
Select
When x.i>x Then
Leave
When x.i=x Then Do
yl.x=yl.x y
Return
End
Otherwise
Nop
End
End
Do j=x.0 To i By -1
ja=j+1
x.ja=x.j
End
x.i=x
yl.x=y
x.0=x.0+1
Return
 
sho: Procedure Expose x. yl. g.
Do i=1 To x.0
x=x.i
say format(i,2) 'x='format(x,3) 'yl='yl.x
End
Say ''
Return
 
maxv: Procedure Expose g.
Call trace 'O'
Parse Arg l
res=-1000
Do While l<>''
Parse Var l v l
If v>res Then res=v
End
Return res
 
minv: Procedure Expose g.
Call trace 'O'
Parse Arg l
res=1000
Do While l<>''
Parse Var l v l
If v<res Then res=v
End
Return res
 
sortv: Procedure Expose g.
Call trace 'O'
Parse Arg l
res=''
Do Until l=''
v=minv(l)
res=res v
l=remove(v,l)
End
Return space(res)
 
lastword: return word(arg(1),words(arg(1)))
 
kdx: Procedure Expose xy. g.
/*---------------------------------------------------------------------
* Compute slope and y-displacement of a straight line
* that is defined by two points: y=k*x+d
* Specialty; k='*' x=xa if xb=xa
*--------------------------------------------------------------------*/

Call trace 'O'
Parse Arg xya,xyb
Parse Var xya xa '/' ya
Parse Var xyb xb '/' yb
If xa=xb Then
Parse Value '*' '-' xa With k d x
Else Do
k=(yb-ya)/(xb-xa)
d=yb-k*xb
x='*'
End
Return k d x
 
remove:
/*---------------------------------------------------------------------
* Remove a specified element (e) from a given string (s)
*--------------------------------------------------------------------*/

Parse Arg e,s
Parse Var s sa (e) sb
Return space(sa sb)
 
o: Procedure Expose g.
/*---------------------------------------------------------------------
* Write a line to the debug file
*--------------------------------------------------------------------*/

If arg(2)=1 Then say arg(1)
Return lineout(g.0oid,arg(1))
 
is_ok: Procedure Expose x. yl. g. sigl
/*---------------------------------------------------------------------
* Test if a given point (b) is above/on/or below a straight line
* defined by two points (a and c)
*--------------------------------------------------------------------*/

Parse Arg a,b,c,op
Call o 'is_ok' a b c op
Parse Value kdx(a,c) With k d x
Parse Var b x'/'y
If op='U' Then y=maxv(yl.x)
Else y=minv(yl.x)
Call o y x (k*x+d)
If (abs(y-(k*x+d))<1.e-8) Then Return 0
If op='U' Then res=(y<=(k*x+d))
Else res=(y>=(k*x+d))
Return res
 
mk_nhull: Procedure Expose x. yl. g.
/*---------------------------------------------------------------------
* Compute the upper (north) hull between two points (xya and xyb)
* Move x from xyb back to xya until all points within the current
* range (x and xyb) are BELOW the straight line defined xya and x
* Then make x the new starting point
*--------------------------------------------------------------------*/

Parse Arg xya,xyb
Call o 'mk_nhull' xya xyb
If xya=xyb Then Return xya
Parse Var xya xa '/' ya
Parse Var xyb xb '/' yb
iu=0
iv=0
Do xi=1 To x.0
if x.xi>=xa & iu=0 Then iu=xi
if x.xi<=xb Then iv=xi
If x.xi>xb Then Leave
End
Call o iu iv
xu=x.iu
xyu=xu'/'maxv(yl.xu)
Do h=iv To iu+1 By -1 Until good
Call o 'iv='iv,g.0debug
Call o ' h='h,g.0debug
xh=x.h
xyh=xh'/'maxv(yl.xh)
Call o 'Testing' xyu xyh,g.0debug
good=1
Do hh=h-1 To iu+1 By -1 While good
xhh=x.hh
xyhh=xhh'/'maxv(yl.xhh)
Call o 'iu hh iv=' iu hh h,g.0debug
If is_ok(xyu,xyhh,xyh,'U') Then Do
Call o xyhh 'is under' xyu xyh,g.0debug
Nop
End
Else Do
good=0
Call o xyhh 'is above' xyu xyh '-' xyh 'ist nicht gut'
End
End
End
Call o xyh 'is the one'
 
Return xyh
 
p: Return
Say arg(1)
Pull .
Return
 
mk_shull: Procedure Expose x. yl. g.
/*---------------------------------------------------------------------
* Compute the lower (south) hull between two points (xya and xyb)
* Move x from xyb back to xya until all points within the current
* range (x and xyb) are ABOVE the straight line defined xya and x
* Then make x the new starting point
*-----
---------------------------------------------------------------*/

Parse Arg xya,xyb
Call o 'mk_shull' xya xyb
If xya=xyb Then Return xya
Parse Var xya xa '/' ya
Parse Var xyb xb '/' yb
iu=0
iv=0
Do xi=x.0 To 1 By -1
if x.xi<=xa & iu=0 Then iu=xi
if x.xi>=xb Then iv=xi
If x.xi<xb Then Leave
End
Call o iu iv '_' x.iu x.iv
Call o 'mk_shull iv iu' iv iu
xu=x.iu
xyu=xu'/'minv(yl.xu)
good=0
Do h=iv To iu-1 Until good
xh=x.h
xyh=xh'/'minv(yl.xh)
Call o 'Testing' xyu xyh h iu
good=1
Do hh=h+1 To iu-1 While good
Call o 'iu hh h=' iu hh h
xhh=x.hh
xyhh=xhh'/'minv(yl.xhh)
If is_ok(xyu,xyhh,xyh,'O') Then Do
Call o xyhh 'is above' xyu xyh
Nop
End
Else Do
Call o xyhh 'is under' xyu xyh '-' xyh 'ist nicht gut'
good=0
End
End
End
Call o xyh 'is the one'
Return xyh
 
Novalue:
Say 'Novalue raised in line' sigl
Say sourceline(sigl)
Say 'Variable' condition('D')
Signal lookaround
 
Syntax:
Say 'Syntax raised in line' sigl
Say sourceline(sigl)
Say 'rc='rc '('errortext(rc)')'
 
halt:
lookaround:
Say 'You can look around now.'
Trace ?R
Nop
Exit 12
Output:
 1 x= -9 yl=-3
 2 x= -4 yl=-6 -2
 3 x= -3 yl=-9 15
 4 x=  0 yl=6 11
 5 x=  3 yl=-4 16
 6 x=  5 yl=19
 7 x= 12 yl=13 17
 8 x= 16 yl=-7 -3 3 6
 9 x= 17 yl=-4 5
10 x= 19 yl=-8

Points of convex hull in clockwise order:
-9/-3 -3/15 5/19 12/17 17/5 19/-8 -3/-9

version 2[edit]

After learning from https://www.youtube.com/watch?v=wRTGDig3jx8

/* REXX ---------------------------------------------------------------
* Compute the Convex Hull for a set of points
* Format of the input file:
* (16,3) (12,17) (0,6) (-4,-6) (16,6) (16,-7) (16,-3) (17,-4) (5,19)
* (19,-8) (3,16) (12,13) (3,-4) (17,5) (-3,15) (-3,-9) (0,11) (-9,-3)
* (-4,-2)
* Alternate (better) method using slopes
* 1) Compute path from lowest/leftmost to leftmost/lowest
* 2) Compute leftmost vertical border
* 3) Compute path from rightmost/highest to highest/rightmost
* 4) Compute path from highest/rightmost to rightmost/highest
* 5) Compute rightmost vertical border
* 6) Compute path from rightmost/lowest to lowest_leftmost point
*--------------------------------------------------------------------*/

Parse Arg fid
If fid='' Then Do
fid='line.in'
fid='point.in'
fid='chullmin.in' /* miscellaneous test data */
fid='chullxx.in'
fid='chullx.in'
fid='chullt.in'
fid='chulla.in'
fid='sq.in'
fid='tri.in'
fid='z.in'
fid='chull.in' /* data from task description */
End
g.0debug=''
g.0oid=fn(fid)'.txt'; 'erase' g.0oid
x.=0
yl.=''
Parse Value '1000 -1000' With g.0xmin g.0xmax
Parse Value '1000 -1000' With g.0ymin g.0ymax
/*---------------------------------------------------------------------
* First read the input and store the points' coordinates
* x.0 contains the number of points, x.i contains the x.coordinate
* yl.x contains the y.coordinate(s) of points (x/y)
*--------------------------------------------------------------------*/

Do while lines(fid)>0
l=linein(fid)
Do While l<>''
Parse Var l '(' x ',' y ')' l
Call store x,y
End
End
Call lineout fid
g.0xlist=''
Do i=1 To x.0 /* loop over points */
x=x.i
g.0xlist=g.0xlist x
yl.x=sortv(yl.x) /* sort y-coordinates */
End
Call sho
If x.0<3 Then Do
Say 'We need at least three points!'
Exit
End
Call o 'g.0xmin='g.0xmin
Call o 'g.0xmi ='g.0xmi
Call o 'g.0ymin='g.0ymin
Call o 'g.0ymi ='g.0ymi
 
Do i=1 To x.0
x=x.i
If minv(yl.x)=g.0ymin Then Leave
End
lowest_leftmost=i
 
highest_rightmost=0
Do i=1 To x.0
x=x.i
If maxv(yl.x)=g.0ymax Then
highest_rightmost=i
If maxv(yl.x)<g.0ymax Then
If highest_rightmost>0 Then
Leave
End
Call o 'lowest_leftmost='lowest_leftmost
Call o 'highest_rightmost ='highest_rightmost
 
x=x.lowest_leftmost
Call o 'We start at' from x'/'minv(yl.x)
path=x'/'minv(yl.x)
/*---------------------------------------------------------------------
* 1) Compute path from lowest/leftmost to leftmost/lowest
*--------------------------------------------------------------------*/

Call min_path lowest_leftmost,1
/*---------------------------------------------------------------------
* 2) Compute leftmost vertical border
*--------------------------------------------------------------------*/

Do i=2 To words(yl.x)
path=path x'/'word(yl.x,i)
End
cxy=x'/'maxv(yl.x)
/*---------------------------------------------------------------------
* 3) Compute path from rightmost/highest to highest/rightmost
*--------------------------------------------------------------------*/

Call max_path ci,highest_rightmost
/*---------------------------------------------------------------------
* 4) Compute path from highest/rightmost to rightmost/highest
*--------------------------------------------------------------------*/

Call max_path ci,x.0
/*---------------------------------------------------------------------
* 5) Compute rightmost vertical border
*--------------------------------------------------------------------*/

Do i=words(yl.x)-1 To 1 By -1
cxy=x'/'word(yl.x,i)
path=path cxy
End
/*---------------------------------------------------------------------
* 6) Compute path from rightmost/lowest to lowest_leftmost
*--------------------------------------------------------------------*/

Call min_path ci,lowest_leftmost
 
Parse Var path pathx path
have.=0
Do i=1 By 1 While path>''
Parse Var path xy path
If have.xy=0 Then Do
pathx=pathx xy
have.xy=1
End
End
Say 'Points of convex hull in clockwise order:'
Say pathx
Call lineout g.0oid
Exit
 
min_path:
Parse Arg from,tgt
ci=from
cxy=x.ci
Do Until ci=tgt
kmax=-1000
Do i=ci-1 To 1 By sign(tgt-from)
x=x.i
k=k(cxy'/'minv(yl.cxy),x'/'minv(yl.x))
If k>kmax Then Do
kmax=k
ii=i
End
End
ci=ii
cxy=x.ii
path=path cxy'/'minv(yl.cxy)
End
Return
 
max_path:
Parse Arg from,tgt
Do Until ci=tgt
kmax=-1000
Do i=ci+1 To tgt
x=x.i
k=k(cxy,x'/'maxv(yl.x))
If k>kmax Then Do
kmax=k
ii=i
End
End
x=x.ii
cxy=x'/'maxv(yl.x)
path=path cxy
ci=ii
End
Return
 
store: Procedure Expose x. yl. g.
/*---------------------------------------------------------------------
* arrange the points in ascending order of x (in x.) and,
* for each x in ascending order of y (in yl.x)
* g.0xmin is the smallest x-value, etc.
* g.0xmi is the x-coordinate
* g.0ymin is the smallest y-value, etc.
* g.0ymi is the x-coordinate of such a point
*--------------------------------------------------------------------*/

Parse Arg x,y
Call o 'store' x y
If x<g.0xmin Then Do; g.0xmin=x; g.0xmi=x; End
If x>g.0xmax Then Do; g.0xmax=x; g.0xma=x; End
If y<g.0ymin Then Do; g.0ymin=y; g.0ymi=x; End
If y>g.0ymax Then Do; g.0ymax=y; g.0yma=x; End
Do i=1 To x.0
Select
When x.i>x Then
Leave
When x.i=x Then Do
yl.x=yl.x y
Return
End
Otherwise
Nop
End
End
Do j=x.0 To i By -1
ja=j+1
x.ja=x.j
End
x.i=x
yl.x=y
x.0=x.0+1
Return
 
sho: Procedure Expose x. yl. g.
Do i=1 To x.0
x=x.i
say format(i,2) 'x='format(x,3) 'yl='yl.x
End
Say ''
Return
 
maxv: Procedure Expose g.
Call trace 'O'
Parse Arg l
res=-1000
Do While l<>''
Parse Var l v l
If v>res Then res=v
End
Return res
 
minv: Procedure Expose g.
Call trace 'O'
Parse Arg l
res=1000
Do While l<>''
Parse Var l v l
If v<res Then res=v
End
Return res
 
sortv: Procedure Expose g.
Call trace 'O'
Parse Arg l
res=''
Do Until l=''
v=minv(l)
res=res v
l=remove(v,l)
End
Return space(res)
 
lastword: return word(arg(1),words(arg(1)))
 
k: Procedure
/*---------------------------------------------------------------------
* Compute slope of a straight line
* that is defined by two points: y=k*x+d
* Specialty; k='*' x=xa if xb=xa
*--------------------------------------------------------------------*/

Call trace 'O'
Parse Arg xya,xyb
Parse Var xya xa '/' ya
Parse Var xyb xb '/' yb
If xa=xb Then
k='*'
Else
k=(yb-ya)/(xb-xa)
Return k
 
remove:
/*---------------------------------------------------------------------
* Remove a specified element (e) from a given string (s)
*--------------------------------------------------------------------*/

Parse Arg e,s
Parse Var s sa (e) sb
Return space(sa sb)
 
o: Procedure Expose g.
/*---------------------------------------------------------------------
* Write a line to the debug file
*--------------------------------------------------------------------*/

If arg(2)=1 Then say arg(1)
Return lineout(g.0oid,arg(1))
Output:
 1 x= -9 yl=-3
 2 x= -4 yl=-6 -2
 3 x= -3 yl=-9 15
 4 x=  0 yl=6 11
 5 x=  3 yl=-4 16
 6 x=  5 yl=19
 7 x= 12 yl=13 17
 8 x= 16 yl=-7 -3 3 6
 9 x= 17 yl=-4 5
10 x= 19 yl=-8

Points of convex hull in clockwise order:
-3/-9 -9/-3 -3/15 5/19 12/17 17/5 19/-8 -3/-9

Scala[edit]

Scala Implementation to find Convex hull of given points collection. Functional Paradigm followed

 
object convex_hull{
def get_hull(points:List[(Double,Double)], hull:List[(Double,Double)]):List[(Double,Double)] = points match{
case Nil => join_tail(hull,hull.size -1)
case head :: tail => get_hull(tail,reduce(head::hull))
}
def reduce(hull:List[(Double,Double)]):List[(Double,Double)] = hull match{
case p1::p2::p3::rest => {
if(check_point(p1,p2,p3)) hull
else reduce(p1::p3::rest)
}
case _ => hull
}
def check_point(pnt:(Double,Double), p2:(Double,Double),p1:(Double,Double)): Boolean = {
val (x,y) = (pnt._1,pnt._2)
val (x1,y1) = (p1._1,p1._2)
val (x2,y2) = (p2._1,p2._2)
((x-x1)*(y2-y1) - (x2-x1)*(y-y1)) <= 0
}
def m(p1:(Double,Double), p2:(Double,Double)):Double = {
if(p2._1 == p1._1 && p1._2>p2._2) 90
else if(p2._1 == p1._1 && p1._2<p2._2) -90
else if(p1._1<p2._1) 180 - Math.toDegrees(Math.atan(-(p1._2 - p2._2)/(p1._1 - p2._1)))
else Math.toDegrees(Math.atan((p1._2 - p2._2)/(p1._1 - p2._1)))
}
def join_tail(hull:List[(Double,Double)],len:Int):List[(Double,Double)] = {
if(m(hull(len),hull(0)) > m(hull(len-1),hull(0))) join_tail(hull.slice(0,len),len-1)
else hull
}
def main(args:Array[String]){
val points = List[(Double,Double)]((16,3), (12,17), (0,6), (-4,-6), (16,6), (16,-7), (16,-3), (17,-4), (5,19), (19,-8), (3,16), (12,13), (3,-4), (17,5), (-3,15), (-3,-9), (0,11), (-9,-3), (-4,-2), (12,10))
val sorted_points = points.sortWith(m(_,(0.0,0.0)) < m(_,(0.0,0.0)))
println(f"Points:\n" + points + f"\n\nConvex Hull :\n" +get_hull(sorted_points,List[(Double,Double)]()))
}
}
 
Output:
Points:
List((16.0,3.0), (12.0,17.0), (0.0,6.0), (-4.0,-6.0), (16.0,6.0), (16.0,-7.0), (16.0,-3.0), (17.0,-4.0), (5.0,19.0), (19.0,-8.0), (3.0,16.0), (12.0,13.0), (3.0,-4.0), (17.0,5.0), (-3.0,15.0), (-3.0,-9.0), (0.0,11.0), (-9.0,-3.0), (-4.0,-2.0), (12.0,10.0))

Convex Hull :
List((-3.0,-9.0), (-9.0,-3.0), (-3.0,15.0), (5.0,19.0), (12.0,17.0), (17.0,5.0), (19.0,-8.0))

zkl[edit]

// Use Graham Scan to sort points into a convex hull
// https://en.wikipedia.org/wiki/Graham_scan, O(n log n)
// http://www.geeksforgeeks.org/convex-hull-set-2-graham-scan/
// http://geomalgorithms.com/a10-_hull-1.html
fcn grahamScan(points){
N:=points.len();
# find the point with the lowest y-coordinate, x is tie breaker
p0:=points.reduce(fcn([(a,b)]ab,[(x,y)]xy){
if(b<y)ab else if(b==y and a<x)ab else xy });
#sort points by polar angle with p0, ie ccw from p0
points.sort('wrap(p1,p2){ ccw(p0,p1,p2)>0 });
 
# We want points[0] to be a sentinel point that will stop the loop.
points.insert(0,points[-1]);
M:=1; # M will denote the number of points on the convex hull.
foreach i in ([2..N]){
# Find next valid point on convex hull.
while(ccw(points[M-1], points[M], points[i])<=0){
if(M>1) M-=1;
else if(i==N) break; # All points are collinear
else i+=1;
}
points.swap(M+=1,i); # Update M and swap points[i] to the correct place.
}
points[0,M]
}
# 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 twice the signed area of the triangle formed by p1, p2 and p3.
fcn ccw(a,b,c){ // a,b,c are points: (x,y)
((b[0] - a[0])*(c[1] - a[1])) - ((b[1] - a[1])*(c[0] - a[0]))
}
pts:=List( T(16,3), T(12,17), T(0,6), T(-4,-6), T(16,6),
T(16, -7), T(16,-3),T(17,-4), T(5,19), T(19,-8),
T(3,16), T(12,13), T(3,-4), T(17,5), T(-3,15),
T(-3,-9), T(0,11), T(-9,-3), T(-4,-2), T(12,10), )
.apply(fcn(xy){ xy.apply("toFloat") }).copy();
hull:=grahamScan(pts);
println("Convex Hull (%d points): %s".fmt(hull.len(),hull.toString(*)));
Output:
Convex Hull (7 points): L(L(-3,-9),L(19,-8),L(17,5),L(12,17),L(5,19),L(-3,15),L(-9,-3))