Efficient Algorithms to Test Digital Convexity

Digital geometry is the field of mathematics that studies the geometry of points with integer coordinates, also known as lattice points [1]. Convexity is a fundamental concept in digital geometry, as well as in continuous geometry [2]. From a historical perspective, the study of digital convexity dates back to the works of Minkowski [3] and it is the main subject of the mathematical field of geometry of numbers.

While convexity has a unique well stated definition in any linear space, different definitions have been investigated in \(\mathbb {Z}^2\) and \(\mathbb {Z}^3\) [4,5,6,7,8]. In two dimensions, we encounter at least five different approaches, called respectively digital line, triangle, line [4], HV (for Horizontal and Vertical [9]), and Q (for Quadrant [10]) convexities. These definitions were created in order to guarantee that a digital convex set is connected (in terms of the induced grid subgraph), which simplifies several algorithmic problems.

The original definition of digital convexity in the geometry of number does not guarantee connectivity of the grid subgraph, but provides several other important mathematical properties, such as being preserved under certain affine transformations (Fig. 1). The definition is the following. A set of lattice points \(S \subset \mathbb {Z}^d\) is digital convex if \({{\,\mathrm{conv}\,}}(S) \cap \mathbb {Z}^d = S\), where \({{\,\mathrm{conv}\,}}(S)\) denotes the convex hull of S.

Shearing a digital convex set. Example of a set whose connectivity is lost after a linear shear.

Full size image

Herein, we consider the fundamental problem of verifying whether a given set of lattice points is digital convex.

The input of TestConvexity(S) is an unstructured finite lattice set (without repeating elements). Related work considered more structured data in dimension 2, in which S is assumed to be connected. The contour of a connected set S of lattice points is the ordered list of the points of S having a grid neighbor outside S. When S is connected, it is possible to represent S by its contour, either directly as in [11] or encoded as binary word [12]. The algorithms presented in [11, 12] test digital convexity in linear time on the respective input representations.

Our work, however, does not make any assumption on S being connected, or any particular ordering of the input. In this setting, a naive approach to test the digital convexity is:

1. 1.

   Compute the convex hull \({{\,\mathrm{conv}\,}}(S)\) of the n lattice points of S.

2. 2.

   Compute the number \(n'\) of lattice points inside the convex hull of S.

3. 3.

   If \(n=n'\), then S is convex. Otherwise, it is not.

Step 1 consists of computing the convex hull of n points. The field of computational geometry provides a plethora of algorithms to compute the convex hull of a finite set \(S\subset \mathbb {R}^d\) of n points [13]. The fastest algorithms for dimensions 2 and 3 take \(O(n \log n)\) time [14], which matches the lower bound in the algebraic decision tree model of computation [15]. In dimension \(d \le 3\), if we also take into consideration the output size h, i.e. the number of vertices of the convex hull, the fastest algorithms take \(O(n \log h)\) time [16, 17]. Some polytopes with n vertices (e.g., the cyclic polytope) have \(\varTheta (n^{\left\lfloor (d-1)/2\right\rfloor })\) facets. Therefore, any algorithm that outputs this facet description of the convex hull requires \(\varOmega (n^{\left\lfloor (d-1)/2\right\rfloor })\) time. Optimal algorithms to compute the convex hull in dimension \(d \ge 4\) match this lower bound [18].

Step 2 consists of computing the number of lattice points inside a polytope (represented by its vertices), which is a well studied problem. In dimension 2, it can be solved using Pick’s formula [19]. The question has been widely investigated in the framework of the geometry of numbers, from Ehrhart theory [20] to Barvinok’s algorithm [21]. Currently best known algorithms have a complexity of \(O(n^{O(d)})\) for fixed dimension d [22]. As conclusion, the time complexity of this naive approach is at least the one of the computation of the convex hull.

1.1 Results

In Sect. 2, we consider the 2-dimensional version of the problem and show that the convex hull of digital convex sets can be computed in linear time. Our main result is an algorithm for dimension \(d=2\) to solve TestConvexity(S) in \(O(n + h \log r)\) time, where h is the number of edges of the convex hull and r is the diameter of S.

In Sect. 3, we consider the problem in fixed dimension d. We present the first polynomial-time algorithm to test digital convexity, as well as a simpler and more practical algorithm whose running time may not be polynomial in n for certain inputs.

The purpose of this section is to provide an algorithm to test the convexity of a finite lattice \(S \subset \mathbb {Z}^2\) in linear time in n. To this endeavour, we show that the convex hull of a digital convex set S can be computed in linear time. In fact, we show that this linear running time is achieved by the well-known quickhull algorithm [23].

Quickhull is one the many early algorithms to compute the convex hull in dimension 2. Its worst case time is \(O(n^2)\), which makes it generally less attractive than the \(O(n \log n)\) algorithm. However for certain inputs and variations of the algorithm, the average time complexity is reduced to \(O(n \log n)\) or O(n) [13, 24].

The quickhull algorithm starts by initializing a convex polygon in the following manner. First it computes the top-most and bottom-most points of the set. Then it computes the two extreme points in the normal direction of the line supported by the top-most and bottom-most points. Those four points describe a convex polygon that we call a partial hull, which is contained inside the convex hull of S. The points contained in the interior of the partial hull are discarded. Furthermore, horizontal lines and lines parallel to the top-most to bottom-most line passing through these points describe an outlying bounding box in which the convex hull lies (Fig. 2).

Quickhull initialization. Points inside the partial hull (light brown) are discarded. The remaining points are potentially part of the hull. (Color figure online)

Full size image

The algorithm adds vertices to the partial hull until it obtains the actual convex hull. This is done by inserting new vertices in the partial hull one by one. Given an edge of the partial hull, let v denote its outwards normal vector. The algorithm searches for the extreme point in direction v. If this point is already an edge point, then the edge is part of the convex hull. Otherwise, we insert the farthest point found between the two edge vertices, discarding the points that are inside the new partial hull. Throughout this paper, we call a step of the quickhull algorithm the computation of the farthest point of every edge for a given partial hull. When adding new vertices to the partial hull, the region inside the partial hull expands. Points inside that expansion are discarded by quickhull and herein we name this region discarded region. The points that still lie outside the partial hull are preserved, and we call the region within which points might still lies preserved region (Fig. 3).

We show that quickhull steps takes linear time and that at each step half of the remaining input points of the convex hull is discarded. Therefore, as in standard decimation algorithms, the total running time remains linear. In Sect. 2.2, we explain how to use this algorithm to test the digital convexity of any lattice set in linear time in n.

Theorem 1

If the input is a digital convex set of n points, then QuickHull has O(n) time and space complexities.

2.1 Proof of Theorem 1

We prove Theorem 1 with the help of the following lemma.

Lemma 1

The area of the discarded region is larger than the area of the preserved region.

Proof

Consider one step of the algorithm: Let ab be the edge associated to the step. When a was added to the hull, it was as the farthest point in a given direction. Hence, there is no point behind the line orthogonal to this direction going through a. (Fig. 3b). The same can be said for b. Let c be the intersection point of those two lines. Every point that lies within \(\triangle abc\) will be fed to the following steps. At this step, we are looking for the point that is the farthest from the supporting line of ab and outside the partial hull (let that point be d) (Fig. 4). Let e and f be the intersections between the line parallel to ab going through d, and respectively ac and bc. There are no points from S inside the triangle \(\triangle cef\). Adding d to the partial hull creates two other edges to further be treated: one with ad as an edge that will be fed the points inside \(\triangle ade\) and one with bd as the edge that will be fed the points inside \(\triangle bdf\). The triangle \(\triangle abd\) lies within the partial hull, therefore \(\triangle abd\) is the region in which points are discarded (Fig. 4).

Quickhull regions. The preserved region (region in which we look for the next vertex to be added to the partial hull) is a triangle. This stays true when adding new vertices to the hull (as shown here in the bottom right corner). The partial hull (whose interior is shown in light brown) grows at each vertex insertions to the partial hull. The new region added to the partial hull is called discarded region. (Color figure online)

Full size image

Symmetrical regions. The next step of the algorithm will only be fed the points inside the dark brown regions (search regions). Each lattice points inside the light brown region (discarded region) is inside the partial hull and is therefore discarded. Each search region (in dark brown) has a symmetrical region (either through \(c_1\) or \(c_2\)) that lies inside the discarded region. Furthermore, this symmetrical transformation also preserve lattice points. (Color figure online)

Full size image

We established that the preserved lattice points are the lattice points within \(\triangle ade\) and \(\triangle bdf\). Also the discarded lattice points are those within \(\triangle abd\). Let \(\text {c}_{\text {1}}\) be the middle of ad and \(\text {c}_{\text {2}}\) be the middle of bd. As shown in Fig. 4, the symmetrical of \(\triangle ade\) and \(\triangle bdf\) through respectively \(c_1\) and \(c_2\) both lie inside \(\triangle abd\) and do not intersect each other. Hence \(\triangle abd\) is larger in terms of area than \(\triangle aed \cup \triangle bdf\).    \(\square \)

Remark 1

Pick’s formula does not apply here since all vertices of the triangle (namely c in Fig. 4) are not necessarily lattice points.

Remark 2

As there is no direct relation between the area of a triangle and the number of lattice points inside it, this result is not sufficient to conclude that a constant proportion of points are discarded at each step.

Corollary 1

The reflection of lattice points inside \(\triangle aed\) and \(\triangle bdf\) across respectively \(\text {c}_{\text {1}}\) and \(\text {c}_{\text {2}}\) are lattice points.

Proof

The points a, b, d are lattice points so \(\text {c}_{\text {1}}\) and \(\text {c}_{\text {2}}\) (middle of respectively ad and bd) have their coordinates in multiple of half integers. Hence the reflection of a lattice point across \(\text {c}_{\text {1}}\) or \(\text {c}_{\text {2}}\) is a lattice point. Therefore, every lattice point within \(\triangle aed\) has a lattice point reflection across \(\text {c}_{\text {1}}\) within \(\triangle a\text {e}_{\text {1}}d\) and every lattice point within \(\triangle bfd\) has a lattice point reflection across \(\text {c}_{\text {2}}\) within \(\triangle b\text {f}_{\text {2}}d\).

   \(\square \)

Remark 3

This previous result would prove that half the points are discarded at each step if it were not for the lattice points on the diagonals ad and bd.

We will now show that quickhull discards at least half of the remaining points at each step, hence proving Theorem 1

Proof

We established in Corollary 1 that lattice points inside the search regions (\(\triangle aed\) and \(\triangle bfd\)) have symmetrical counterparts inside the discarded region (more precisely inside \(\triangle ae_1d\) and \(\triangle bf_2d\)) (Fig. 4). By preserving each points inside \(\triangle aed\) and \(\triangle bfd\) at each step, we do not have a discarded symmetrical counterpart for the lattice points lying on ad and bd. But we do not need to preserve those points, since ad and bd are at this step edges of the partial hull. Removing lattice points from ad and bd implies that in the following step there will be no lattice points on ab, leaving lattice points on ef without a discarded symmetrical counterparts (Fig. 5).

Lonely points. The lattice points without discarded symmetrical counterparts are shown in red. On the left: if every points inside the triangle is preserved, and in the center: if the points on the edges of the partial hull are discarded. Finally on the right a visualization of what happens if we discard all the farthest points and update the partial hull accordingly. (Color figure online)

Full size image

Let actually discard every points on ef, since they all are equally farthest from ab in the outer direction, they all belong to the hull. Hence we can add the first and last lattice point on ef to the partial hull (Fig. 5). Note that this only takes linear time and does not change the time complexity of each individual step. Hence, at each step of quickhull, for every preserved points there is at least a discarded point. Consequently, the number of operations is proportional to \(n\sum \limits _{i=0}^\infty (\frac{1}{2})^i = 2n\) and quickhull takes linear time for digital convex sets.    \(\square \)

2.2 Determining the Digital Convexity of a Set

We showed in Theorem 1 that the quickhull algorithm computes the convex hull of digital convex sets in linear time thanks to the fact that at each step quickhull discards at least half of the remaining points. By running quickhull on any given set S, and stopping the computation if any step of the algorithm discards less than half of the remaining points, we ensure both that the running time is linear, and that if S is digital convex, quickhull finishes and we get the convex hull of S. If the computation finishes for S, we still need to test its digital convexity. To do so, we use the previously computed convex hull and compute \(|conv(S) \cap \mathbb {Z}^2|\) using Pick’s formula [19]. The set S is digital convex if \(|conv(S) \cap \mathbb {Z}^2| = |S|\). Hence the resulting Algorithm 1.

Theorem 2

Algorithm 1 tests digital convexity of any 2 dimensional set S, and runs in \(O(n + h \log r)\) time, where h is the number of edges of conv(S) and r is the diameter of S.

Proof

As Algorithm 1 runs quickhull, but stops as soon as less than half the remaining points have been removed, the running time of the quickhull part is bounded by the series \(n\sum \limits _{i=0}^\infty (\frac{1}{2})^i = 2n\), and is hence linear. Thanks to Theorem 1 we know that the computation of quickhull will not stop for any digital convex sets. Computing \(|conv(S) \cap \mathbb {Z}^2|\) using Pick’s formula requires the computation of the area of \({{\,\mathrm{conv}\,}}(S)\) and of the number of lattice points lying on its boundary, which requires the computation of a greatest common divisor. Hence this takes \(O(h \log r)\) time where h is the number of edge of conv(S) and r is the diameter of S. As S is digital convex if and only if \(|S| = |conv(S) \cap \mathbb {Z}^2|\), Algorithm 1 effectively tests the digital convexity of a 2 dimensional set in \(O(n + h \log r)\) time.    \(\square \)

We provide two algorithms for verifying the digital convexity in any fixed dimension.

3.1 Naive Algorithm

The naive algorithm mentioned in the Introduction is based on the following equivalence: the set \(S \subset \mathbb {Z}^d\) is digital convex if and only if its cardinality is equal to the cardinality of \(conv(S) \cap \mathbb {Z}^d\). In Step 1, we compute the convex hull of S (in \(O(n \log n + n ^{\lfloor \frac{d}{2}\rfloor })\) time [18]). In Step 2, we need to count the number of integer points inside \({{\,\mathrm{conv}\,}}(S)\). The classical algorithm to achieve this goal is known as Barvinok algorithm [21]. This approach determines only the number of missing points. If we want to enumerate the points, it is possible to do so through a formal computation of the generating functions used in Barvinok algorithm.

Theorem 3

The naive algorithm tests digital convexity in any fixed dimension d and runs in polynomial time.

Proof

Computing the convex hull of any set can be done in \(O(n \log n + n ^{\lfloor \frac{d}{2}\rfloor })\) time [18]). Counting lattice points inside a convex lattice polytope can be done in polynomial time [22]. A direct consequence of the digital convexity definition is that a set \(S \subset \mathbb {Z}^d\) is digital convex if and only if \(|S| = |conv(S) \cap \mathbb {Z}^d|\), hence the naive algorithm tests digital convexity in any fixed dimension d and runs in polynomial time.    \(\square \)

3.2 Alternative Algorithm

This new algorithm computes all integer points in the convex hull of S with a more direct approach. Its principle is to enumerate the points x of a finite lattice set \(S' \subset \mathbb {Z}^d\) surrounding \(conv(S) \cap \mathbb {Z}^d\) (\(conv(S) \cap \mathbb {Z}^d \subset S' \)). In a first variant, we count the number of points of \(S'\) belonging to conv(S). At the end, the set S is convex if and only if \(|conv (S) \cap \mathbb {Z}^d|\) is equal to the cardinality of S. In a second variant, for each point of \(S'\), we test whether it belongs to S and in the negative case, we test whether it belongs to the convex hull of S. If a point of \(S' \setminus S \cap conv(S)\) is found, then S is not convex (Fig. 6).

We define the set \(S'\) as the set of points \(x\in \mathbb {Z}^d\) such that the cube \(x+[-\frac{1}{2}, \frac{1}{2} ]^d \) has a nonempty intersection with the convex hull of S, where \(+\) denotes the Minkowski sum. It can be easily proved that \(S'\) is 2d-connected (the 2d neighbors of a lattice point \(x\in \mathbb {Z}^d\) are the 2d integer points at Euclidean distance 1) and by construction, it contains S. The graph structure induced by the 2d-connectivity on \(S'\) allows to visit all the points of \(S'\) efficiently: for each point \(x\in S'\), we consider its 2d neighbors and test whether they belong to \(S'\). If they do, we add them to the stack of the remaining points of \(S'\). The goal is to test whether a point of \(S' \setminus S\) is in the convex hull of S.

Practical algorithm. A lattice set S, its convex hull and its dilation by a centered cube of side 1. The intersection of \({{\,\mathrm{conv}\,}}(S) + [-\frac{1}{2}, \frac{1}{2}] ^d\) with the lattice is the set \(S'\). It is 2d-connected and contains the convex hull of S. The principle of the algorithm is either to count the points of \(S'\) in \({{\,\mathrm{conv}\,}}(S)\) (variant 1) or to search for a point of \(S' \setminus S\) (blue points) in the convex hull of S (variant 2). (Color figure online)

Full size image

Then the algorithm has two main routines:

• \(\mathtt {InConvexHull}_S\) tests whether a given point \(x\in \mathbb {R}^d\) belongs to the convex hull of S. It is equivalent with testing whether there exists a hyperplane separating x from the points of S. It can be done by linear programming with a worst-case time complexity of O(n) for fixed dimension d [13].

• \(\mathtt {InConvexHull}_{S+ [-\frac{1}{2}, \frac{1}{2}] ^d }\) tests whether a given point x belongs to the convex hull of \(S+ [-\frac{1}{2}, \frac{1}{2}] ^d\). It follows the same principle as \(\mathtt {InConvexHull}_S\) with \(2^d n\) points. The time complexity remains linear in fixed dimension. This routine is used to test whether an integer point belongs to \(S'\).

The algorithm is the following. First, we create a stack T of the points of \(S'\) to visit and initialize it with the set S. For each point x in T, we remove it from the stack T and label it as already visited. Then, we consider its 2d neighbors \(x'\). If \(x'\) belongs to \(S'\) and has not been visited previously, we add it in the stack T. We test whether x belongs to conv(S) and increment the cardinality of \(conv (S) \cap \mathbb {Z}^d\) accordingly (variant 1) or test whether x is in S and conv(S) and return \(\mathtt {S}\) not convex if \(x \in conv(S)\setminus S\) (variant 2).

The running time is strongly dependent on the cardinality of \(S'\). It is \( O(n |S'| )\). If the size of \(S'\) is of the same magnitude as the initial set, the algorithm runs in \(O(n^2)\) time. It is unfortunately not possible to bound \(|S' |\) as a function of n. The ratio \(\frac{|S' |}{|S|}\) can go to infinity. It is easy to build such an example with a set S consisting of only two lattice points, for instance for any \(k \in \mathbb {Z}\) the set \(S=\{(0,0);(1,2k)\}\) induces \(\frac{|S' |}{|S|} \ge k\). A direction of improvement could be to consider a linear transformation of the lattice \(\mathbb {Z}^d\) in order to obtain a more compact lattice set and then a lower ratio \(\frac{|S' |}{|S|}\). LLL algorithm [25] could be useful to achieve this goal in future work.

As in the naive algorithm, a variant of this approach can be easily developed in order to enumerate the missing points.

In this paper, we presented an algorithm to test digital convexity in time linear in n for dimension \(d=2\). In higher dimensions, our running time depends on the complexity of general convex hull algorithms. The questions of whether digital convexity can be tested in linear time in 3 dimensions, or faster than convex hull computation in arbitrary dimensions remain open. A tentative approach consists of changing the lattice base, in order to obtain certain connectivity properties.

We showed that the convex hull of a digital convex set in dimension 2 can be computed in linear time. Can the convex hull of digital convex sets be computed in linear time in dimension 3, or more generally, what is the complexity of convex hull computation of a digital convex set in any fixed dimension? We note that the number of faces of any digital convex set in d dimensions is \(O(V^{(d-1)/(d+1)})\), where V is the volume of the polytope [26, 27]. Therefore, the lower bound of \(\varOmega (n^{\left\lfloor (d-1)/2\right\rfloor })\) for the complexity of the convex hull of arbitrary polytopes does not hold for digital convex sets.