diff --git a/clusters/cluster.go b/clusters/cluster.go
index f0be99e..2202f07 100644
--- a/clusters/cluster.go
+++ b/clusters/cluster.go
@@ -46,6 +46,11 @@ func (c *Cluster) Append(point space.Point) {
 	c.PointList = append(c.PointList, point)
 }
 
+// ConvexHull returns the convex hull of the clusters PointList
+func (c *Cluster) ConvexHull() PointList {
+	return c.PointList.ConvexHull()
+}
+
 // Nearest returns the index of the cluster nearest to point
 func (c Clusters) Nearest(point space.Point) int {
 	var ci int
diff --git a/clusters/point_list.go b/clusters/point_list.go
index 97cd94f..52fe7df 100644
--- a/clusters/point_list.go
+++ b/clusters/point_list.go
@@ -2,6 +2,7 @@ package clusters
 
 import (
 	"fmt"
+	"sort"
 
 	"github.com/spatial-go/geoos/planar"
 	"github.com/spatial-go/geoos/space"
@@ -51,3 +52,59 @@ func AverageDistance(point space.Point, points PointList) float64 {
 	}
 	return d / float64(l)
 }
+
+// ConvexHull returns the convex hull of a list of points
+// Implementation of Andrew's Monotone Chain algorithm as specified on https://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain
+func (points PointList) ConvexHull() PointList {
+	// empty slice, single point, and two points are already their own convex hull
+	if len(points) <= 2 {
+		return points
+	}
+
+	// sort points by x coordinates (ties stay in original order)
+	sort.SliceStable(points, func(i, j int) bool {
+		return points[i].X() < points[j].X() || (points[i].X() == points[j].X() && points[i].Y() < points[j].Y())
+	})
+
+	// build lower hull
+	lowerHull := PointList{}
+	for _, p := range points {
+		for len(lowerHull) >= 2 && cross(lowerHull[len(lowerHull)-2], lowerHull[len(lowerHull)-1], p) <= 0 {
+			// pop last point
+			lowerHull = lowerHull[:len(lowerHull)-1]
+		}
+		// add p
+		lowerHull = append(lowerHull, p)
+	}
+	// build upper hull
+	upperHull := PointList{}
+	// iterate in reverse order
+	for i := len(points) - 1; i >= 0; i-- {
+		p := points[i]
+		for len(upperHull) >= 2 && cross(upperHull[len(upperHull)-2], upperHull[len(upperHull)-1], p) <= 0 {
+			// pop last point
+			upperHull = upperHull[:len(upperHull)-1]
+		}
+		// add p
+		upperHull = append(upperHull, p)
+	}
+
+	// concatenate lower and upper hull to build convexHull.
+	// omit the last point of lowerHull as it is the beginning of upperHull
+	// keep last point of upperHull to get a closed polygon shape (last point = first point)
+	hull := append(lowerHull[:len(lowerHull)-1], upperHull...)
+
+	if len(hull) == 3 {
+		// for lines, remove last point (so no closed polygon shape)
+		hull = hull[:len(hull)-1]
+	}
+	return hull
+}
+
+// cross is a helper function for ConvexHull()
+// cross product of OA and OB vectors.
+// returns positive value, if OAB makes a counter-clockwise turn,
+// negative for clockwise turn, and zero if the points are colinear.
+func cross(o, a, b space.Point) float64 {
+	return (a.X()-o.X())*(b.Y()-o.Y()) - (a.Y()-o.Y())*(b.X()-o.X())
+}
diff --git a/clusters/point_list_test.go b/clusters/point_list_test.go
new file mode 100644
index 0000000..5a08ef0
--- /dev/null
+++ b/clusters/point_list_test.go
@@ -0,0 +1,60 @@
+package clusters
+
+import (
+	"math/rand"
+	"reflect"
+	"testing"
+
+	"github.com/spatial-go/geoos/space"
+)
+
+func TestConvexHull(t *testing.T) {
+	// empty list
+	points := PointList{}
+	expected := PointList{}
+	testExample(t, points, expected)
+
+	// single point
+	points = PointList{{1, 1}}
+	expected = PointList{{1, 1}}
+	testExample(t, points, expected)
+
+	// two points
+	points = PointList{{1, 1}, {1, 2}}
+	expected = PointList{{1, 1}, {1, 2}}
+	testExample(t, points, expected)
+
+	// line
+	points = PointList{{1, 1}, {2, 2}, {3, 3}}
+	expected = PointList{{1, 1}, {3, 3}} // intermediate point omitted
+	testExample(t, points, expected)
+
+	// triangle
+	points = PointList{{2, 2}, {1, 2}, {1, 1}}
+	expected = PointList{{1, 1}, {2, 2}, {1, 2}, {1, 1}} // closed polygon shape
+	testExample(t, points, expected)
+
+	// square
+	points = PointList{{1, 1}, {2, 1}, {2, 2}, {1, 2}}
+	expected = PointList{{1, 1}, {2, 1}, {2, 2}, {1, 2}, {1, 1}} // closed polygon-shape
+	testExample(t, points, expected)
+
+	// generate random point cloud with set outer bound
+	size := 50
+	min, max := 0., 90.
+	points = make(PointList, size)
+	for i := 0; i < size; i++ {
+		points[i] = space.Point{min + rand.Float64()*max, min + rand.Float64()*max}
+	}
+	points = append(points, PointList{{min - 1, min - 1}, {min - 1, max + 1}, {max + 1, max + 1}, {max + 1, min - 1}}...)    // square outer boundary in clockwise order
+	expected = PointList{{min - 1, min - 1}, {max + 1, min - 1}, {max + 1, max + 1}, {min - 1, max + 1}, {min - 1, min - 1}} // closed polygon-shape in counter-clockwise order
+	testExample(t, points, expected)
+}
+
+func testExample(t *testing.T, points PointList, expected PointList) {
+	hull := points.ConvexHull()
+	if !reflect.DeepEqual(expected, hull) {
+		t.Errorf("Expected %v but got %v", expected, hull)
+		t.FailNow()
+	}
+}