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.Ord
import Data.List
 
(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
 
example
:: (Floating a, Ord a)
=> [[a]]
example =
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]
]
 
main :: IO ()
main = mapM_ print example
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

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))