# Compute and draw a convex hull

How can I use TikZ (+ LuaTeX), PSTricks, MetaPost, Asymptote, etc. to draw the convex hull of a set of points specified in no particular order?

Here's an example where the convex hull was computed by hand:

\documentclass[png]{standalone}
\usepackage{tikz}
\begin{document}
\begin{tikzpicture}
% doesn't work:
%\draw[mark=*] convex hull[points={(0,0),(1,1),(2,2),(0,1),(2,0)}];

%should produce something similar to the following:
\draw (0,0) -- (0,1) -- (2,2) -- (2,0) -- cycle;
\foreach \point in {(0,0),(1,1),(2,2),(0,1),(2,0)} {
}
\end{tikzpicture}
\end{document}


• Use one of the well known algorithms for finding the convex hull with the implementation using Lua. Then pass those points to pgfplots to draw the hull and all of the points. Most of the algorithms should implement in nlogn time. A link for one of them "QuickHull" is westhoffswelt.de/blog/2009/10/21/…. Jun 12, 2015 at 3:46
• Gift wrapping is (much) easier to code than quickhull, and plenty fast if you have few points. And you better have few points, or else your figure will be too crowded. Wikipedia has nice psuedocode for it.
– Mark
Jun 12, 2015 at 4:19
• @Mark: A convex hull of many points need not be crowded if the points are not explicitly plotted (which was done mostly for illustration). One possible application is to draw a "random" convex polygon by choosing a bunch of random points inside a circle (or oval, rectangle, ...) and then taking their convex hull. Another application is to approximate the convex hull of a smooth curve by taking a large number of closely spaced points on the curve. Jun 12, 2015 at 8:46

Robert Sedgewick's Algorithms in C has a whole chapter on convex hulls; here is the algorithm that he calls "package wrapping" implemented in Metapost.

prologues := 3;
outputtemplate := "%j%c.eps";

% following Sedgwick, "Algorithms in C", p.364
% make the first M points the hull of the first N points
vardef wrap(expr N) =
save    theta, eta, tx, ty, M, m, a;
numeric theta, eta, tx, ty, M, m, a;
% find the point with the minimum y-coordinate
m = 1;
for i=2 upto N:
if y[i] < y[m]: m:=i; fi
endfor
x[N+1] = x[m];
y[N+1] = y[m];
% work round the points finding the least turning angle
% and swapping the hull points to the front of the list of points
theta := 0;
M = 1;
forever:
tx := x[M]; x[M] := x[m]; x[m] := tx;
ty := y[M]; y[M] := y[m]; y[m] := ty;
m := N+1;
eta := theta;
theta := 360;
for i=M+1 upto N+1:
a := if z[i]=z[M]: 0 else: angle (z[i]-z[M]) mod 360 fi;
if eta < a: if a < theta:
m := i;
theta := a;
fi fi
endfor
exitif (m=N+1);
M := M+1;
exitif (M>N);
endfor
% return the number of points in the hull
% z1, z2, .. z[M] are now the hull
M
enddef;

beginfig(1);

N = 42;
for i=1 upto N:
z[i] = (80 normaldeviate, 80 normaldeviate);
dotlabel.top(decimal i, z[i]);
endfor

m = wrap(N);
draw for i=1 upto m: z[i] -- endfor cycle withcolor .67 red;

for i=1 upto N:
label.bot(decimal i, z[i]) withcolor blue;
endfor

endfig;
end.


In the original you pass an array of points to the wrap function, but passing arrays is a bit cumbersome in MP, so I've just used the normal global arrays x[] and y[] for the points. Note that the algorithm is destructive, in that it rearranges the order of the points so that the points on the hull are at the beginning of the array. I've tried to make this obvious by including before and after labels in the output.

This implementation has running time proportional to N^2 in the worst case, but is reasonably quick on my machine with up to 1200 points. You could speed it up with the interior elimination techniques that Sedgewick discusses in the chapter already mentioned.

I implemented a convex hull generator in asymptote a while back. It uses the gift wrapping algorithm. As the hull is being generated, the function improves performance by eliminating points that are already inside the hull.

path convexHull(pair[] in_pset)
{
pair[] pset = copy(in_pset);
if (pset.length == 0) { path hull; return hull; }
{ // remove duplicate points
int indexDelete = 1;
while (indexDelete > 0)
{
indexDelete = -1;
for (int i = 1; i < pset.length; ++i)
{
for (int j = 0; j < i; ++j)
{
if (pset[i] == pset[j])
{
indexDelete = i;
break;
}
}
if (indexDelete > 0) { break; }
}
if (indexDelete > 0) { pset.delete(indexDelete); }
}
}

path hull;
{ // add point at min y (and min x if tie) to hull, delete point from pset
int minIndex = 0;

for (int i = 1; i < pset.length; ++i)
{
if (pset[i].y < pset[minIndex].y ||
(pset[i].y == pset[minIndex].y && pset[i].x < pset[minIndex].x))
{
minIndex = i;
}
}
hull = pset[minIndex];
pset.delete(minIndex);
}

while (pset.length > 0)
{
{ // add next point to hull
real minAngle = 361.0;
int minAngleIndex = 0;
for (int i = 0; i < pset.length; ++i)
{
real angle = degrees(pset[i] - relpoint(hull, 1.0), false);
if (angle < minAngle)
{
minAngle = angle;
minAngleIndex = i;
}
}
hull = hull--pset[minAngleIndex];
pset.delete(minAngleIndex);
}

{ // remove points interior to current hull from pset
path tempHull = hull--cycle;
int[] deleteIndeces;
for (int i = pset.length - 1; i > -1; --i)
{
if (inside(tempHull, pset[i])) { deleteIndeces.push(i); }
}
for (int i = 0; i < deleteIndeces.length; ++i)
{
pset.delete(deleteIndeces[i]);
}
}
}
return hull--cycle;
}


Now the following code uses this function to draw a convex hull of some random points.

unitsize(1inch);

pair[] pset;
for (int i = 0; i < 100; ++i)
{
pair p = scale(2)*slant(0.5)*(unitrand(), unitrand());
dot(p, 3+red);
pset.push(p);
}

draw(convexHull(pset));


TikZ and TeX and something called Graham Scan.

The macro \CH does all the stuff but the drawing.

You can give it a set of coordinates with the coordinates key which will create named coordinates with the prefix ConvexHullPoint- and saves the number of the coordinate that lies on the hull in \outerPoints, all other are stored in \innerPoints:

\CH[coordinates={(1,1),(2,2),(1,2),(3,3),(4,2),(2,3),(3,2)}]
\path plot[mark=*, samples at=\innerPoints] (ConvexHullPoint-\x);
\draw plot[mark=*, samples at=\outerPoints] (ConvexHullPoint-\x) --cycle;


You can also prepare a set of coordinates beforehand. Name them from <name>-1 to <name>-<n> where <n> is the last and total number of points. This means, that they could also be nodes (of any shape), however, the .center anchor will be used for any calculations:

\foreach \i in {1,...,10} \path (10*rnd,10*rnd) coordinate[label=\tiny\i] (chp-\i);
\CH[total=10]% name=chp
\path plot[mark=*, mark options=blue,  samples at=\innerPoints] (chp-\x);
\draw plot[mark=*, mark options=green, samples at=\outerPoints] (chp-\x) --cycle;


The \CH macro takes one optional argument in the form of a key-value list. The /ch key tree offers these four value keys:

• name, the “base” that all points have in common
• total, the total number <n>,
• outer macro, the macro in which the points on the hull are stored in,
• inner macro, the macro for all other points (inside the hull), and
• coordinates as mentioned above.

## Code

\documentclass[tikz]{standalone}
\usetikzlibrary{backgrounds}
\makeatletter
\newcommand*\chset{\pgfqkeys{/ch}}
\chset{name/.initial=chp, total/.initial=4, outer macro/.initial=\outerPoints,
inner macro/.initial=\innerPoints,
\pgfcoordinate{ConvexHullPoint-\the\pgfutil@tempcnta}
{\tikz@scan@one@point\pgfutil@firstofone#1\relax}},
coordinates/.code={\pgfutil@tempcnta=0
\pgfkeysalso{
/ch/@prepare coordinates/.list={#1},
/ch/name=ConvexHullPoint,
/ch/total/.expanded=\the\pgfutil@tempcnta}}}
\newcommand*\chvof[1]{\pgfkeysvalueof{/ch/#1}}
\newcommand*\CH[1][]{%
\begingroup\chset{#1}%
%% Get the lowest left point
% \CH@Ai stores ID, \CH@Axy stores x, y, \CH@Apoint expands to PGF-point
\def\CH@Ai{0}\pgf@ya=16000pt \pgf@xa=16000pt
\pgfmathloop
\pgf@process{\pgfpointanchor{\chvof{name}-\pgfmathcounter}{center}}%
\ifdim\pgf@y<\pgf@ya
\let\CH@Ai\pgfmathcounter \pgf@xa=\pgf@x \pgf@ya=\pgf@y
\else
\ifdim\pgf@y=\pgf@ya
\ifdim\pgf@x<\pgf@xa
\let\CH@Ai\pgfmathcounter \pgf@xa=\pgf@x \pgf@ya=\pgf@y
\fi
\fi
\fi
\ifnum\pgfmathcounter<\chvof{total}\relax
\repeatpgfmathloop
\edef\CH@Axy{{\the\pgf@xa}{\the\pgf@ya}}%
\edef\CH@Apoint{\noexpand\pgfqpoint\CH@Axy}%
%% Build list of points sorted after angle from lowest left point
% \CH@list will contain stack of (ID, angle, x, y) in TeX groups
\let\CH@list\pgfutil@empty
\pgfmathloop
\ifnum\pgfmathcounter=\CH@Ai\else
\pgfextract@process\CH@p{\pgfpointanchor{\chvof{name}-\pgfmathcounter}{center}}
\edef\pgf@tempa{{\the\pgf@x}{\the\pgf@y}}%
\pgfmathanglebetweenpoints{\CH@Apoint}{\CH@p}%
\edef\CH@element{{\pgfmathcounter}{\pgfmathresult}}%
\let\CH@angle\pgfmathresult
\edef\CH@element{\CH@element\pgf@tempa}%
\ifx\CH@list\pgfutil@empty
\let\CH@list\CH@element
\else
\let\CH@lista\pgfutil@empty
\expandafter\CH@sortin\CH@list\@@{}{}{}%
\let\CH@list\CH@lista
\fi
\fi
\ifnum\pgfmathcounter<\chvof{total}\relax
\repeatpgfmathloop
%% Drop points on the inner side.
% This tests if point[i] is on the right of line through point[i-1] and point[i+1]
% \CH@listb will contain list of outer points (reverse stack)
% \CH@listc will contain list of inner points
\edef\CH@listb{{\CH@Ai}{}\CH@Axy}%
\let\CH@listc\pgfutil@gobble
\expandafter\CH@store\CH@list\CH@stop\CH@Ti\CH@Txy\CH@list
\pgfmathloop
\expandafter\CH@store\CH@list\CH@stop\CH@Bi\CH@Bxy\CH@list
\edef\pgf@marshall{\noexpand\CHtestforLeftOrRight\CH@Axy\CH@Bxy\CH@Txy}%
\pgf@marshall
%      \errmessage{\CH@Ai, \CH@Ti, \CH@Bi; \pgfmathresult; \CH@Axy, \CH@Txy, \CH@Bxy}%
\ifnum\pgfmathresult=-1 % to the right
% woho, add point[i] to the outer list and push everything one down
\edef\CH@listb{{\CH@Ti}{}\CH@Txy\CH@listb}%
\let\CH@Ai\CH@Ti
\let\CH@Axy\CH@Txy
\let\CH@Ti\CH@Bi
\let\CH@Txy\CH@Bxy
\else % otherwise
% ugh, insert point[i+1] back into the source list
% so that it will be point[i+1] in the next iteration and push everything one up
\edef\CH@listc{\CH@listc,\CH@Ti}%
\edef\CH@list{{\CH@Bi}{}\CH@Bxy\CH@list}%
\expandafter\CH@testforfirst\CH@listb\CH@stop\CH@Ti\CH@listb
\expandafter\CH@store\CH@listb\CH@stop\CH@Ti\CH@Txy\CH@listb
\expandafter\CH@store\CH@listb\CH@stop\CH@Ai\CH@Axy\CH@listb
\edef\CH@listb{{\CH@Ai}{}\CH@Axy\CH@listb}%
\fi
\ifx\CH@list\pgfutil@empty % Before we finish, add the last entry
\edef\CH@listb{{\CH@Bi}{}\CH@Bxy\CH@listb}%
\else
\repeatpgfmathloop
%% Alright lets pull only the IDs from \CH@listb and add to "outer" in reverse order
\pgfkeysgetvalue{/ch/outer macro}\CH@outer \pgfkeysgetvalue{/ch/inner macro}\CH@inner
\expandafter\let\CH@outer\pgfutil@empty
\pgfmathloop
\expandafter\CH@store\CH@listb\CH@stop\CH@Ai\CH@Axy\CH@listb
\expandafter\ifx\CH@outer\pgfutil@empty
\expandafter\edef\CH@outer{\CH@Ai}%
\else\expandafter\edef\CH@outer{\CH@Ai,\CH@outer}\fi
\ifx\CH@listb\pgfutil@empty\else
\repeatpgfmathloop
% get "outer" and \CH@listc (in the form of "inner") outside the group
\ifx\CH@listc\pgfutil@gobble\let\CH@listc\pgfutil@empty\fi
\xdef\pgf@marshall{\def\expandafter\noexpand\CH@outer{\CH@outer}%
\def\expandafter\noexpand\CH@inner{\CH@listc}}%
\endgroup
\pgf@marshall}
\newcommand*\CH@sortin[4]{%
\ifx\@@#1%
\edef\CH@lista{\CH@lista\CH@element}%
\expandafter\pgfutil@gobble
\else
\expandafter\pgfutil@firstofone
\fi{%
\ifdim#2pt<\CH@angle pt
\edef\CH@lista{\CH@lista{#1}{#2}{#3}{#4}}\expandafter\CH@sortin
\else
\fi}}
\newcommand*\CHtestforLeftOrRight[6]{%
\begingroup
% numbers too big, scale everything down
\dimen6=.1\dimen6 \dimen3=.1\dimen3
\dimen5=.1\dimen5 \dimen4=-.1\dimen4
\pgf@x=\pgf@sys@tonumber{\dimen5}\dimen4        % - (#5-#1)(#4-#2)
\pgfmath@returnone\pgf@x\endgroup
%  \pgfmathparse{(#3-#1)(#6-#2)-(#5-#1)(#4-#2)}%
\ifdim\pgfmathresult pt<0pt \def\pgfmathresult{-1}%
\else\ifdim\pgfmathresult pt>0pt \def\pgfmathresult{1}%
\else\def\pgfmathresult{0}\fi\fi}
\def\CH@store#1#2#3#4#5\CH@stop#6#7#8{\edef#6{#1}\edef#7{{#3}{#4}}\edef#8{#5}}
\def\CH@testforfirst#1#2#3#4#5\CH@stop#6#7{\ifnum#1=#6 \edef#7{#5}\fi}

\makeatother\tikzset{mark=*, mark size=1pt}

\begin{document}
%
\foreach \n in {4,...,10}{\pgfmathsetseed{249860}%
\begin{tikzpicture}
\useasboundingbox (2,-0.1) -- (10,9.5);
\foreach \i in {1,...,\n} \path (10*rnd,10*rnd) coordinate[label=\tiny\i] (chp-\i);
\CH[total=\n]
\path plot[mark options=blue,  samples at=\innerPoints] (chp-\x);
\draw plot[mark options=green, samples at=\outerPoints] (chp-\x) --cycle;
\end{tikzpicture}}
%
\tikzset{every picture/.append style=gridded}
\begin{tikzpicture}
\CH[coordinates={(0,0),(1,1),(2,2),(0,1),(2,0)}]
\path plot[samples at=\innerPoints] (ConvexHullPoint-\x);
\draw plot[samples at=\outerPoints] (ConvexHullPoint-\x) --cycle;
\end{tikzpicture}
%
\begin{tikzpicture}
\CH[coordinates={(1,1),(2,2),(1,2),(3,3),(4,2),(2,3),(3,2)}]
\path plot[samples at=\innerPoints] (ConvexHullPoint-\x);
\draw plot[samples at=\outerPoints] (ConvexHullPoint-\x) --cycle;
\end{tikzpicture}
%
\begin{tikzpicture}
\CH[coordinates={(3,0),(4,1),(5,2),(3,1),(5,3),(6,0),(4.5,-1),(5,4),(3,2),(3.2,1.7)}]
\path plot[samples at=\innerPoints] (ConvexHullPoint-\x);
\draw plot[samples at=\outerPoints] (ConvexHullPoint-\x) --cycle;
\end{tikzpicture}
\end{document}


## Output

Here is a TikZ kind of solution.

The idea is to obtain cliped region (approximately) equal to the convex hull. For this I rotate the points (the precision is of 1°) and clip the bounding box.

\documentclass[tikz,border=5]{standalone}
\usetikzlibrary{calc}
\tikzset{
point/.style={insert path={node[scale=3,#1]{.}}}
}

\newcommand\inconvexhull[2]{
%--------------------------
% #1 is a set of points
% #2 is the command to execute in the clipped convex hull
%--------------------------
\edef\pts{#1} % the points

% get the bounding box of rotated points
\foreach \i in {0,...,90}{
\begin{scope}[local bounding box=B\i]
\foreach \p in \pts
\path ([rotate=\i]\p);
\end{scope}
}
\begin{scope}
% clip rotated bounding boxes B0, ...,B90
\foreach \i in {0,...,90}
\clip ([rotate=-\i]B\i.south east) -- ([rotate=-\i]B\i.north east)
-- ([rotate=-\i]B\i.north west)  -- ([rotate=-\i]B\i.south west) -- cycle;
% execute command #2 on B0 rectangle
\foreach \i in {0} #2
(B\i.south west) rectangle (B\i.north east);
\end{scope}
}

\begin{document}
% TEST
\begin{tikzpicture}
% generate 10 random points, draw them and stock them in \randpoints
\edef\randpoints{}
\foreach[count=\n, evaluate={\sep=\n>1?",":"";}] \i in {1,...,10} {
\path ({random(0:360)}:{random(1:4)}) coordinate (R\i)[point=red];
\xdef\randpoints{\randpoints\sep R\i}; % add random point R\i to \randpoints
}

% fill the convex hull
\inconvexhull{\randpoints}{\fill[red,opacity=.1]}
\end{tikzpicture}
\end{document}


Note : We can use math library to implement one of the multiple algorithms to calculate the convex hull, but I like this tricky solution using multiple rotated bounding boxes.

The pst-intersect package does such calculations internally (using the Postscript procedures from http://www.math.ubc.ca/~cass/graphics/text/www/). They can be wrapped inside an own macro \convexhull as follows:

\documentclass[margin=12pt,pstricks]{standalone}
\usepackage{pst-intersect}
\makeatletter
\def\convexhull{\pst@object{convexhull}}
\def\convexhull@i{%
\begin@ClosedObj
\pst@getcoors[\convexhull@ii%
}
\def\convexhull@ii{%
\pst@intersectdict
] ArrayToPointArray
ConvexHull PointArrayToArray
end
}%
\pspolygon@ii
}%
\makeatother
\begin{document}
\begin{pspicture}[showgrid](5,4)
\convexhull[fillcolor=red!20, fillstyle=solid](1,1)(2,2)(1,2)(3,3)(4,2)(2,3)(3,2)
\psdots(1,1)(2,2)(1,2)(3,3)(4,2)(2,3)(3,2)
\end{pspicture}
\end{document}


The result is

At the moment, \convexhull cannot yet plot the original points itself, but that must be done with a separate call to \psdots, but that shouldn't be difficult to implement when I have some spare time.

• Plotting the original points and the convex hull in one line was suggested only because it's the sort of thing that happens in TikZ syntax; it's not a requirement for a good answer. Jun 12, 2015 at 8:49

Another late addition to the TikZ class. It defines a new plot handler so that you can use the typical TikZ syntax as you have in the question. Hence the syntax goes like

\draw plot[convex hull,mark=*]  coordinates {(1,1)(2,2)(1,2)(3,3)(4,2)(2,3)(3,2)};


or

\draw plot[convex hull,mark=*]  plot(\x,{0.05*exp(\x)});


That being said, for some reason, if I use the random functions, the seed gets run twice (I guess) and the marks don't match the hull. I don't know why yet.

Another point that needs attention is that if the points are too close to each other, TikZ doesn't compute the angles correctly. Here you can see that it got stuck between the left upper two points and over printed. But if you reduce the number of samples it is OK again.

Long story short, it might not be robust for every case.

\documentclass[tikz]{standalone}
\makeatletter
\pgfkeyssetvalue{/tikz/convex hull}{\let\tikz@plot@handler=\pgfplothandlerconvexhull}
\pgfdeclareplothandler{\pgfplothandlerconvexhull}{}{
start = {\pgf@plot@mark@count = 1},
end   = {% Find the top point
\pgfmathtruncatemacro\pgf@ch@total{\the\pgf@plot@mark@count-1}%
\pgf@plot@mark@count = 0\pgfmathtruncatemacro\xe{int(\pgf@ch@total-1)}%
\foreach\x in{1,...,\xe}{%
\def\pgf@ch@invalid@flag{0}\pgfmathtruncatemacro\xx{int(\x+1)}%
\foreach\y in{\xx,...,\pgf@ch@total}{%
\pgfmathanglebetweenpoints{\pgfpointanchor{convex hull-\x}{center}}%
{\pgfpointanchor{convex hull-\y}{center}}%
%\typeout{current test \x - \y : angle \pgfmathresult}
\ifdim 180pt >\pgfmathresult pt\relax\xdef\pgf@ch@invalid@flag{1}\breakforeach\fi%
}%
\ifnum\pgf@ch@invalid@flag<1\relax\xdef\pgf@ch@first@{\x}\breakforeach\else%
\xdef\pgf@ch@first@{\pgf@ch@total}\fi%
}% Found. Now look at upperright then switch to that and rotate...loop
\edef\pgf@ch@current@{\pgf@ch@first@}\def\pgf@ch@current@angle{0}%
\edef\pgf@ch@collected{\pgf@ch@first@}\def\pgf@ch@temp@max@{0}%
\foreach\x in {1,...,\pgf@ch@total}{
\pgftransformshift{\pgfpointdiff{\pgfpointorigin}%
{\pgfpointanchor{convex hull-\pgf@ch@current@}{center}}}%
\pgftransformrotate{\pgf@ch@current@angle-360}%
\pgfmathloop%
\ifnum\pgfmathcounter=\pgf@ch@current@\relax\else%
\pgfmathanglebetweenpoints{\pgfpointanchor{convex hull-\pgf@ch@current@}{center}}%
{\pgfpointanchor{convex hull-\pgfmathcounter}{center}}%
\edef\pgf@ch@current@angle{\pgfmathresult}%
\ifdim\pgf@ch@current@angle pt>\pgf@ch@temp@max@ pt\relax%
\edef\pgf@ch@temp@max@{\pgf@ch@current@angle}%
\edef\pgf@ch@temp@current@{\pgfmathcounter}%
%\typeout{new max for \pgf@ch@current@ with \the\c@pgf@countd: angle \pgf@ch@current@angle}
\fi%
\fi%
\ifnum\pgfmathcounter<\pgf@ch@total\repeatpgfmathloop%
\xdef\pgf@ch@current@angle{\pgf@ch@temp@max@}%
\xdef\pgf@ch@current@{\pgf@ch@temp@current@}%
\expandafter\xdef\expandafter\pgf@ch@collected\expandafter{\pgf@ch@collected,\pgf@ch@current@}%
\ifnum\pgf@ch@current@=\pgf@ch@first@\relax\breakforeach\fi%
%\typeout{\pgf@ch@collected}
}%
\foreach\x[count=\xi] in\pgf@ch@collected{%
\ifnum\xi=1\relax\pgfpathmoveto{\pgfpointanchor{convex hull-\x}{center}}\else%
\pgfpathlineto{\pgfpointanchor{convex hull-\x}{center}}\fi%
}\pgfpathclose%
}
}
\makeatother
\begin{document}
\begin{tikzpicture}
\draw plot[convex hull,mark=*,samples=25,domain=1:5,mark size={1pt}]  (\x,{sin(\x r)});
%coordinates
%{(3,0)(4,1)(5,2)(3,1)(5,3)(6,0)(4.5,-1)(5,4)(3,2)(3.2,1.7)};
%{(1,1)(2,2)(1,2)(3,3)(4,2)(2,3)(3,2)};
%{(0,0)(1,1)(2,2)(0,1)(2,0)};
%{(3.5,0)(4,1)(5,2)(3.2,1)(5,3)(6,0)(4.5,-1)(5,4)(3,2)(3.2,1.7)};
\end{tikzpicture}
\end{document}


Here's a LaTeX3 and TikZ implementation of Graham's Scan algorithm:

\documentclass[tikz]{standalone}
\usepackage{xparse,l3sort}

\ExplSyntaxOn

\seq_new:N \g_convexhull_input_seq
\seq_new:N \g_convexhull_hull_seq
\int_new:N \g_convexhull_k_int
\int_new:N \g_convexhull_l_int
\bool_new:N \g_convexhull_stop_bool

\cs_generate_variant:Nn \seq_put_right:Nn { Nf }
\cs_generate_variant:Nn \seq_push:Nn  { Nf }
\cs_generate_variant:Nn \tl_set:Nn { Nf }

\NewDocumentCommand{\convexhull}{ m } {

\seq_clear:N \g_convexhull_input_seq
\seq_clear:N \g_convexhull_hull_seq
\seq_set_split:Nnn \g_convexhull_input_seq { ; } { #1 }

% preprocess so that points[1] has lowest y-coordinate; break ties by x-coordinate
\seq_sort:Nn \g_convexhull_input_seq {
\y_compare:nnnTF {##1} {>} {##2} {\sort_reversed:} {\sort_ordered:}
}
\seq_set_eq:NN \l_tmpa_seq \g_convexhull_input_seq
\seq_set_eq:NN \l_tmpb_seq \c_empty_seq
\seq_map_inline:Nn \g_convexhull_input_seq {
\fp_compare:nNnTF {\convexhull_y:f {##1}} {=} {\convexhull_y:f {\seq_item:Nn \g_convexhull_input_seq{1}}}
{ \seq_pop_left:NN \l_tmpa_seq \l_tmpa_tl
\seq_put_right:Nn \l_tmpb_seq {##1}
} {  }
}

% sort by polar angle with respect to base point points[1],
\seq_sort:Nn \l_tmpa_seq {
\angle_compare:nnnTF {##1} {>} {##2} {\sort_reversed:} {\sort_ordered:}
}
\seq_sort:Nn \l_tmpb_seq {
\x_compare:nnnTF {##1} {>} {##2} {\sort_reversed:} {\sort_ordered:}
}
\seq_concat:NNN \g_convexhull_input_seq \l_tmpb_seq \l_tmpa_seq

% p[1] is first extreme point
\seq_set_eq:NN \g_convexhull_hull_seq \c_empty_seq
\seq_push:Nn \g_convexhull_hull_seq { \seq_item:Nn \g_convexhull_input_seq { 1 } }

% find index k of first point not equal to points[1]
\bool_set_false:N \g_convexhull_stop_bool
\int_step_variable:nnnNn { 2 } { 1 } { \seq_count:N \g_convexhull_input_seq } \l_tmpa_tl {
\str_if_eq_x:nnTF
{ \seq_item:Nn \g_convexhull_input_seq { 1 } }
{ \seq_item:Nn \g_convexhull_input_seq { \l_tmpa_tl } }
{  }
{ \bool_if:NTF \g_convexhull_stop_bool
{  }
{
\int_set:Nn \g_convexhull_k_int \l_tmpa_tl
\bool_set_true:N \g_convexhull_stop_bool
}
}
}

% find index l of first point not collinear with points[1] and points[k]
\bool_set_false:N \g_convexhull_stop_bool
\int_step_variable:nnnNn { \g_convexhull_k_int + 1 } { 1 } { \seq_count:N \g_convexhull_input_seq } \l_tmpa_tl {
\fp_compare:nNnTF
{ \convexhull_ccw:fff { \seq_item:Nn \g_convexhull_input_seq { 1 } } { \seq_item:Nn \g_convexhull_input_seq { \g_convexhull_k_int } } { \seq_item:Nn \g_convexhull_input_seq { \l_tmpa_tl } } }
{ = }
{ 0 }
{  }
{ \bool_if:NTF \g_convexhull_stop_bool
{  }
{
\int_set:Nn \g_convexhull_l_int \l_tmpa_tl
\bool_set_true:N \g_convexhull_stop_bool
}
}
}
\seq_push:Nn \g_convexhull_hull_seq { \seq_item:Nn \g_convexhull_input_seq { \g_convexhull_l_int - 1 } }

% Graham scan
\int_step_variable:nnnNn { \g_convexhull_l_int } { 1 } { \seq_count:N \g_convexhull_input_seq } \l_tmpa_tl {
\seq_pop:NN \g_convexhull_hull_seq \g_tmpa_tl
\seq_get:NN \g_convexhull_hull_seq \g_tmpb_tl
\fp_while_do:nNnn
{ \convexhull_ccw:fff {\g_tmpb_tl}{\g_tmpa_tl}{\seq_item:Nn \g_convexhull_input_seq{\l_tmpa_tl}} }
{ < }
{ 0.000001 }
{
\seq_pop:NN \g_convexhull_hull_seq \g_tmpa_tl
\seq_get:NN \g_convexhull_hull_seq \g_tmpb_tl
}
\seq_push:Nf \g_convexhull_hull_seq { \g_tmpa_tl }
\seq_push:Nf \g_convexhull_hull_seq { \seq_item:Nn \g_convexhull_input_seq { \l_tmpa_tl } }
}

% draw points
\int_step_variable:nnnNn { 1 } { 1 } { \seq_count:N \g_convexhull_input_seq } \l_tmpa_tl {
\node at \seq_item:Nn \g_convexhull_input_seq {\l_tmpa_tl} {$\bullet$};
}

% draw convex hull
\draw \seq_item:Nn \g_convexhull_hull_seq {1} -- \seq_item:Nn \g_convexhull_hull_seq {\seq_count:N \g_convexhull_hull_seq};
\int_step_variable:nnnNn { 1 } { 1 } { \seq_count:N \g_convexhull_hull_seq - 1 } \l_tmpa_tl {
\int_set:Nn \g_convexhull_l_int \l_tmpa_tl
\int_incr:N \g_convexhull_l_int
\tl_set:Nf \l_tmpb_tl {\int_use:N \g_convexhull_l_int}
\tl_set:Nf \l_tmpa_tl {\seq_item:Nn \g_convexhull_hull_seq {\l_tmpa_tl}}
\tl_set:Nf \l_tmpb_tl {\seq_item:Nn \g_convexhull_hull_seq {\l_tmpb_tl}}
\draw \l_tmpa_tl -- \l_tmpb_tl;
}
}

\prg_new_conditional:Npnn \y_compare:nnn #1 #2 #3 {TF} {
\fp_compare:nNnTF {\convexhull_y:f {#1}} {#2} {\convexhull_y:f {#3}} {\prg_return_true:} {\prg_return_false:}
}

\prg_new_conditional:Npnn \angle_compare:nnn #1 #2 #3 {TF} {
\fp_compare:nNnTF {\convexhull_polar_angle:ff {\seq_item:Nn \g_convexhull_input_seq{1}} {#1}} {#2} {\convexhull_polar_angle:ff {\seq_item:Nn \g_convexhull_input_seq{1}} {#3}} {\prg_return_true:} {\prg_return_false:}
}

\prg_new_conditional:Npnn \x_compare:nnn #1 #2 #3 {TF} {
\fp_compare:nNnTF {\convexhull_x:ff {\seq_item:Nn \g_convexhull_input_seq{1}} {#1}} {#2} {\convexhull_x:ff {\seq_item:Nn \g_convexhull_input_seq{1}} {#3}} {\prg_return_true:} {\prg_return_false:}
}

\cs_new:Npn \convexhull_ccw:p (#1,#2) (#3,#4) (#5,#6) {     \fp_eval:n { (#3-#1)*(#6-#2)-(#4-#2)*(#5-#1) } }
\cs_new:Npn \convexhull_ccw:nnn #1 #2 #3 { \convexhull_ccw:p #1 #2 #3 }
\cs_generate_variant:Nn \convexhull_ccw:nnn { fff }

\cs_new:Npn \convexhull_polar_angle:p (#1,#2) (#3,#4) {     \fp_eval:n {-1 * (#3-#1) / (#4-#2) } }
\cs_new:Npn \convexhull_polar_angle:nn #1 #2 { \convexhull_polar_angle:p #1 #2 }
\cs_generate_variant:Nn \convexhull_polar_angle:nn { ff }

\cs_new:Npn \convexhull_y:p (#1,#2) { #2 }
\cs_new:Npn \convexhull_y:n #1 { \convexhull_y:p #1 }
\cs_generate_variant:Nn \convexhull_y:n { f }

\cs_new:Npn \convexhull_x:p (#1,#2) (#3,#4) { \fp_eval:n {#3-#1} }
\cs_new:Npn \convexhull_x:nn #1 #2 {    \convexhull_x:p #1 #2 }
\cs_generate_variant:Nn \convexhull_x:nn { ff }

\ExplSyntaxOff

\begin{document}
\begin{tikzpicture}[scale=2]
\convexhull{(0,0);(1,1);(2,2);(0,1);(2,0)}
\convexhull{(3,0);(4,1);(5,2);(3,1);(5,3);(6,0);(4.5,-1);(5,4);(3,2);(3.2,1.7)}
\end{tikzpicture}
\end{document}


There are probably still quite a lot of mistakes in the code as it is my first attempt at LaTeX3. Anyway, it was a lot of fun and at first glance it seems to work fine.

• \sort_ordered: and \sort_reversed: are deprecated as of 2018-12-31; the right functions are \sort_return_same: and \sort_return_swapped: Dec 27, 2018 at 22:47

Using R-knitr with LaTeX to generate convex hulls.

\documentclass{article}
%% Reference:
% https://chitchatr.wordpress.com/2011/12/30/convex-hull-around-scatter-plot-in-r/
% for the R code which is now linked to LaTeX with knitr
\begin{document}

<<>>=
### Plotting function to plot convex hulls
### Filename: Plot_ConvexHull.R
### Notes:
############################################################################

# INPUTS:
# xcoords: x-coordinates of point data
# ycoords: y-coordinates of point data
# lcolor: line color

# OUTPUTS:
# convex hull around data points in a particular color (specified by lcolor)

# FUNCTION:
Plot_ConvexHull<-function(xcoord, ycoord, lcolor){
hpts <- chull(x = xcoord, y = ycoord)
hpts <- c(hpts, hpts[1])
lines(xcoord[hpts], ycoord[hpts], col = lcolor)
}
# END OF FUNCTION
@

<<>>=
# Create 3 sets of random data to plot convex hull around
x1 <- rnorm(100, 0.8, 0.3)
y1 <- rnorm(100, 0.8, 0.3)

x2 <- rnorm(100, 0.2, 0.3)
y2 <- rnorm(100, 0.2, 0.3)

x3 <- rnorm(100, 1.4, 0.3)
y3 <- rnorm(100, 1.4, 0.3)

# get max and min of all x and y data for nice plotting
xrange <- range(c(x1, x2, x3))
yrange <- range(c(y1, y2, y3))

# Plot it up!
par(tck = 0.02, mgp = c(1.7, 0.3, 0))
plot(x1, y1, type = "p", pch = 1, col = "black", xlim = c(xrange), ylim = c(yrange))
Plot_ConvexHull(xcoord = x1, ycoord = y1, lcolor = "black")
points(x2, y2, type = "p", pch = 1, col = "green")
Plot_ConvexHull(xcoord = x2, ycoord = y2, lcolor = "green")
points(x3, y3, type = "p", pch = 1, col = "magenta")
Plot_ConvexHull(xcoord = x3, ycoord = y3, lcolor = "magenta")
@
\end{document}