Convex hull: Difference between revisions
m (→{{header|Haskell}}: (generalised one map to fmap, using <$> and one less pair of parentheses for legibility)) |
Walterpachl (talk | contribs) (add REXX) |
||
Line 237: | Line 237: | ||
silence implies tests pass (and output is as expected) |
silence implies tests pass (and output is as expected) |
||
=={{header|REXX}}== |
|||
<lang rexx>/* 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 |
|||
fid='chullx.in' |
|||
fid='chullxx.in' |
|||
fid='chull.in' |
|||
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 |
|||
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) |
|||
*--------------------------------------------------------------------*/ |
|||
lefthigh=0 |
|||
leftlow=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 |
|||
/*--------------------------------------------------------------------- |
|||
* 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 lefthigh to ritehigh (n_high) |
|||
* and ritelow to leftlow (s_low), respectively |
|||
*--------------------------------------------------------------------*/ |
|||
y=g.0ymax |
|||
n_high='' |
|||
Do xi=1 To x.0 |
|||
x=x.xi |
|||
If maxv(yl.x)=y Then n_high=n_high x'/'y |
|||
End |
|||
Call o 'n_high='n_high |
|||
y=g.0ymin |
|||
s_low='' |
|||
Do xi=x.0 To 1 By -1 |
|||
x=x.xi |
|||
If minv(yl.x)=y Then s_low=s_low x'/'y |
|||
End |
|||
Call o 's_low='s_low |
|||
/*--------------------------------------------------------------------- |
|||
* 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 |
|||
/*--------------------------------------------------------------------- |
|||
* Now we compute the upper part of the convex hull (nhull) |
|||
*--------------------------------------------------------------------*/ |
|||
res=mk_nhull(leftmost_n,lefthigh); nhull=leftmost_n res |
|||
Do until res=lefthigh |
|||
res=mk_nhull(res,lefthigh); nhull=nhull res |
|||
End |
|||
nhull=nhull n_high |
|||
res=mk_nhull(ritehigh,ritemost_n); nhull=nhull res |
|||
Do While res<>ritemost_n |
|||
res=mk_nhull(res,ritemost_n); nhull=nhull res |
|||
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 |
|||
Do While res<>ritelow |
|||
res=mk_shull(res,ritelow); |
|||
shull=shull res |
|||
End |
|||
shull=shull s_low |
|||
res=mk_shull(leftlow,leftmost_s); shull=shull leftlow res |
|||
Do While res<>leftmost_s |
|||
res=mk_shull(res,leftmost_s); |
|||
shull=shull res |
|||
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 |
|||
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 |
|||
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 'x='format(x,2) '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))) |
|||
between: Procedure Expose g. |
|||
Parse Arg a,b,c |
|||
If sign(b-a)=sign(c-b) Then |
|||
Return 1 |
|||
Else |
|||
Return 0 |
|||
kdx: Procedure Expose xy. g. |
|||
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 |
|||
sort_xyl: Procedure |
|||
Parse Arg xyl |
|||
n=0 |
|||
Do While xyl<>'' |
|||
Parse Var xyl xy xyl |
|||
Parse Var xy x '/' y |
|||
Do i=1 To n |
|||
Select |
|||
When x<x.i Then Leave |
|||
When x=x.i & y<y.i Then Leave |
|||
Otherwise Nop |
|||
End |
|||
End |
|||
Do j=n To i By -1 |
|||
ja=j+1 |
|||
Parse Value x.j y.j With x.ja y.ja |
|||
End |
|||
Parse Value x y with x.i y.i |
|||
n=n+1 |
|||
End |
|||
res='' |
|||
Do i=1 To n |
|||
res=res x.i'/'y.i |
|||
End |
|||
Return space(res) |
|||
show: Procedure |
|||
Parse Arg xyl |
|||
Do i=1 To words(xyl) |
|||
Say format(i,2) word(xyl,i) |
|||
End |
|||
Return |
|||
remove: |
|||
Parse Arg e,s |
|||
Parse Var s sa (e) sb |
|||
Return space(sa sb) |
|||
o: Return lineout(g.0oid,arg(1)) |
|||
is_ok: Procedure Expose x. yl. g. |
|||
Parse Arg a,b,c,op |
|||
Parse Value kdx(a,c) With k d x |
|||
Parse Var b x'/'y |
|||
y=maxv(yl.x) |
|||
Call o y x (k*x+d) |
|||
If op='U' Then res=(y<=(k*x+d))|(abs(y-(k*x+d))<1.e-8) |
|||
Else res=(y>=(k*x+d))|(abs(y-(k*x+d))<1.e-8) |
|||
Return res |
|||
mk_nhull: Procedure Expose x. yl. g. |
|||
Parse Arg xya,xyb |
|||
Call o 'mk_nhull' xya xyb |
|||
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 |
|||
Do xi=iu To iv |
|||
x=x.xi; Call o xi format(x,2) '->' yl.x |
|||
End |
|||
xu=x.iu |
|||
xyu=xu'/'maxv(yl.xu) |
|||
Do h=iv To iu+1 By -1 Until good |
|||
xh=x.h |
|||
xyh=xh'/'maxv(yl.xh) |
|||
Call o 'Testing' xyu xyh |
|||
good=1 |
|||
Do hh=iv To iu+1 By -1 While good |
|||
xhh=x.hh |
|||
xyhh=xhh'/'maxv(yl.xhh) |
|||
If is_ok(xyu,xyhh,xyh,'U') Then Do |
|||
Call o xyhh 'is under' xyu xyh |
|||
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 |
|||
mk_shull: Procedure Expose x. yl. g. |
|||
Parse Arg xya,xyb |
|||
Call o 'mk_shull' xya xyb |
|||
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 'iu iv' iu iv |
|||
Do xi=iu To iv By -1 |
|||
x=x.xi; Call o xi format(x,2) '->' yl.x |
|||
End |
|||
xu=x.iu |
|||
xyu=xu'/'minv(yl.xu) |
|||
Do h=iv To iu+1 Until good |
|||
xh=x.h |
|||
xyh=xh'/'minv(yl.xh) |
|||
Call o 'Testing' xyu xyh |
|||
good=1 |
|||
Do hh=iv To iu+1 While good |
|||
xhh=x.hh |
|||
xyhh=xhh'/'minv(yl.xhh) |
|||
Call o hh xyhh |
|||
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: |
|||
If fore() Then Do |
|||
Say 'You can look around now.' |
|||
Trace ?R |
|||
Nop |
|||
End |
|||
Exit 12</lang> |
|||
{{out}} |
|||
<pre>x=-9 yl=-3 |
|||
x=-4 yl=-6 -2 |
|||
x=-3 yl=-9 15 |
|||
x= 0 yl=6 11 |
|||
x= 3 yl=-4 16 |
|||
x= 5 yl=19 |
|||
x=12 yl=13 17 |
|||
x=16 yl=-7 -3 3 6 |
|||
x=17 yl=-4 5 |
|||
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 -9/-3</pre> |
|||
=={{header|zkl}}== |
=={{header|zkl}}== |
Revision as of 20:13, 15 January 2017
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
<lang go>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) }</lang>
- Output:
Convex Hull: [(-9,-3) (-3,-9) (19,-8) (17,5) (12,17) (5,19) (-3,15)]
Haskell
<lang Haskell>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</lang>
- 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
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
<lang J>counterclockwise =: ({. , }. /: 12 o. }. - {.) @ /:~ crossproduct =: 11"_ o. [: (* +)/ }. - {. removeinner =: #~ 1, 0 > 3 crossproduct\ ], 1: hull =: [: removeinner^:_ counterclockwise</lang>
Example use:
<lang J> 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</lang>
Python
An approach that uses the shapely library:
<lang python>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)</lang>
- Output:
POLYGON ((-3 -9, -9 -3, -3 15, 5 19, 12 17, 17 5, 19 -8, -3 -9))
Racket
Also an implementation of https://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain (therefore kinda
<lang racket>#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))))</lang>
- Output:
silence implies tests pass (and output is as expected)
REXX
<lang rexx>/* 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
fid='chullx.in' fid='chullxx.in' fid='chull.in' 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
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)
- --------------------------------------------------------------------*/
lefthigh=0 leftlow=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
/*---------------------------------------------------------------------
- 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 lefthigh to ritehigh (n_high)
- and ritelow to leftlow (s_low), respectively
- --------------------------------------------------------------------*/
y=g.0ymax n_high= Do xi=1 To x.0
x=x.xi If maxv(yl.x)=y Then n_high=n_high x'/'y End
Call o 'n_high='n_high y=g.0ymin s_low= Do xi=x.0 To 1 By -1
x=x.xi If minv(yl.x)=y Then s_low=s_low x'/'y End
Call o 's_low='s_low
/*---------------------------------------------------------------------
- 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
/*---------------------------------------------------------------------
- Now we compute the upper part of the convex hull (nhull)
- --------------------------------------------------------------------*/
res=mk_nhull(leftmost_n,lefthigh); nhull=leftmost_n res Do until res=lefthigh
res=mk_nhull(res,lefthigh); nhull=nhull res End
nhull=nhull n_high res=mk_nhull(ritehigh,ritemost_n); nhull=nhull res Do While res<>ritemost_n
res=mk_nhull(res,ritemost_n); nhull=nhull res 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 Do While res<>ritelow
res=mk_shull(res,ritelow); shull=shull res End
shull=shull s_low res=mk_shull(leftlow,leftmost_s); shull=shull leftlow res Do While res<>leftmost_s
res=mk_shull(res,leftmost_s); shull=shull res 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 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 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 'x='format(x,2) '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)))
between: Procedure Expose g.
Parse Arg a,b,c If sign(b-a)=sign(c-b) Then Return 1 Else Return 0
kdx: Procedure Expose xy. g.
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
sort_xyl: Procedure
Parse Arg xyl n=0 Do While xyl<> Parse Var xyl xy xyl Parse Var xy x '/' y Do i=1 To n Select When x<x.i Then Leave When x=x.i & y<y.i Then Leave Otherwise Nop End End Do j=n To i By -1 ja=j+1 Parse Value x.j y.j With x.ja y.ja End Parse Value x y with x.i y.i n=n+1 End res= Do i=1 To n res=res x.i'/'y.i End Return space(res)
show: Procedure
Parse Arg xyl Do i=1 To words(xyl) Say format(i,2) word(xyl,i) End Return
remove:
Parse Arg e,s Parse Var s sa (e) sb Return space(sa sb)
o: Return lineout(g.0oid,arg(1))
is_ok: Procedure Expose x. yl. g.
Parse Arg a,b,c,op Parse Value kdx(a,c) With k d x Parse Var b x'/'y y=maxv(yl.x) Call o y x (k*x+d) If op='U' Then res=(y<=(k*x+d))|(abs(y-(k*x+d))<1.e-8) Else res=(y>=(k*x+d))|(abs(y-(k*x+d))<1.e-8) Return res
mk_nhull: Procedure Expose x. yl. g.
Parse Arg xya,xyb Call o 'mk_nhull' xya xyb 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 Do xi=iu To iv x=x.xi; Call o xi format(x,2) '->' yl.x End xu=x.iu xyu=xu'/'maxv(yl.xu) Do h=iv To iu+1 By -1 Until good xh=x.h xyh=xh'/'maxv(yl.xh) Call o 'Testing' xyu xyh good=1 Do hh=iv To iu+1 By -1 While good xhh=x.hh xyhh=xhh'/'maxv(yl.xhh) If is_ok(xyu,xyhh,xyh,'U') Then Do Call o xyhh 'is under' xyu xyh 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
mk_shull: Procedure Expose x. yl. g.
Parse Arg xya,xyb Call o 'mk_shull' xya xyb 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 'iu iv' iu iv Do xi=iu To iv By -1 x=x.xi; Call o xi format(x,2) '->' yl.x End xu=x.iu xyu=xu'/'minv(yl.xu) Do h=iv To iu+1 Until good xh=x.h xyh=xh'/'minv(yl.xh) Call o 'Testing' xyu xyh good=1 Do hh=iv To iu+1 While good xhh=x.hh xyhh=xhh'/'minv(yl.xhh) Call o hh xyhh 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:
If fore() Then Do Say 'You can look around now.' Trace ?R Nop End Exit 12</lang>
- Output:
x=-9 yl=-3 x=-4 yl=-6 -2 x=-3 yl=-9 15 x= 0 yl=6 11 x= 3 yl=-4 16 x= 5 yl=19 x=12 yl=13 17 x=16 yl=-7 -3 3 6 x=17 yl=-4 5 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 -9/-3
zkl
<lang zkl>// 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]))
}</lang> <lang zkl>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(*)));</lang>
- 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))