Convex hull: Difference between revisions

From Rosetta Code
Content added Content deleted
m (→‎REXX version 2: use subroutines to shorten the program)
(→‎{{header|Perl 6}}: Add a Perl6 example)
Line 167: Line 167:
<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
<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>
_9j_3 _3j_9 19j_8 17j5 12j17 5j19 _3j15</lang>

=={{header|Perl 6}}==
{{works with|Rakudo|2017.05}}
{{trans|zkl}}

<lang perl6>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 graham-scan (**@coords) {
my @sp = @coords.map( { Point.new( :x($_[0]), :y($_[1]) ) } ).sort( {.y, .x} );
my @h = @sp.shift;
my $X = @h[0].clone;
$X.y -= 1;
@sp = @sp.map( { $++ => ccw($X, @h[0], $_) } ).sort( *.value ).map( { @sp[$_.key] } );
@h.push: @sp.shift;
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 ({+@hull} points): ", @hull;</lang>
{{out}}
<pre>Convex Hull (7 points): [(-3, -9) (19, -8) (17, 5) (12, 17) (5, 19) (-3, 15) (-9, -3)]</pre>


=={{header|Python}}==
=={{header|Python}}==

Revision as of 01:15, 9 June 2017

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

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

Perl 6

Works with: Rakudo version 2017.05
Translation of: zkl

<lang perl6>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 graham-scan (**@coords) { 

   my @sp = @coords.map( { Point.new( :x($_[0]), :y($_[1]) ) } ).sort( {.y, .x} );
   my @h  = @sp.shift;
   my $X  = @h[0].clone;
   $X.y  -= 1;
   @sp    = @sp.map( { $++ => ccw($X, @h[0], $_) } ).sort( *.value ).map( { @sp[$_.key] } );
   @h.push: @sp.shift;
   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 ({+@hull} points): ", @hull;</lang>

Output:
Convex Hull (7 points): [(-3, -9) (19, -8) (17, 5) (12, 17) (5, 19) (-3, 15) (-9, -3)]

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

Translation of: Go

<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

version 1

<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

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</lang>
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

After learning from https://www.youtube.com/watch?v=wRTGDig3jx8 <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)
  • 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))</lang>
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

Scala Implementation to find Convex hull of given points collection. Functional Paradigm followed <lang Scala> 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)]()))
   }

} </lang>

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

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

}

  1. Three points are a counter-clockwise turn if ccw > 0, clockwise if
  2. ccw < 0, and collinear if ccw = 0 because ccw is a determinant that
  3. 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))
Retrieved from "https://rosettacode.org/wiki/Convex_hull?oldid=258593"