23

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)} {
    \fill[black] \point circle[radius=1pt];
  }
\end{tikzpicture}
\end{document}

output

3
  • 1
    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
  • 1
    @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

8 Answers 8

22

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.

enter image description here

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.

18

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

enter image description here

1
18

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,
  @prepare coordinates/.code={\advance\pgfutil@tempcnta1
    \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
      \edef\CH@lista{\CH@lista\CH@element{#1}{#2}{#3}{#4}}\expandafter\CH@addLeftover
    \fi}}
\newcommand*\CHtestforLeftOrRight[6]{%
  \begingroup
    \dimen6=#6 \advance\dimen6-#2 % (#6-#2)
    \dimen3=#3 \advance\dimen3-#1 % (#3-#1)
    \dimen5=#5 \advance\dimen5-#1 % (#5-#1)
    \dimen4=#4 \advance\dimen4-#2 % (#4-#2)
    % 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)
    \advance\pgf@x\pgf@sys@tonumber{\dimen3}\dimen6 % add (#3-#1)(#6-#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@addLeftover#1\@@#2#3#4{\edef\CH@lista{\CH@lista#1}}
\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

enter image description here

enter image description hereenter image description hereenter image description here

16

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
% (example : \shade[left color=red])
%--------------------------
  \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.

enter image description here

14

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{%
  \addto@pscode{%
    \pst@intersectdict
      ] ArrayToPointArray
      ConvexHull PointArrayToArray
      mark exch aload pop
    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

enter image description here

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.

1
  • 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
12

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},
  point = {\pgfcoordinate{convex hull-\the\pgf@plot@mark@count}{##1}\global\advance\pgf@plot@mark@count by1\relax},
  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}

enter image description here

11

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.

enter image description here

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

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}

enter image description here

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .