GNU Astronomy Utilities


Next: , Previous: , Up: Gnuastro library   [Contents][Index]


12.3.17 Polygons (polygon.h)

Polygons are commonly necessary in image processing. for example, in Crop they are used for cutting out non-rectangular regions of a image (see Crop), and in Warp, for mapping different pixel grids over each other (see Warp).

Polygons come in two classes: convex and concave (or generally, non-convex!), see below for a demonstration. Convex polygons are those where all inner angles are less than 180 degrees. By contrast, a convex polygon is one where an inner angle may be more than 180 degress.

            Concave Polygon        Convex Polygon

             D --------C          D------------- C
              \        |        E /              |
               \E      |          \              |
               /       |           \             |
              A--------B             A ----------B

In all the functions here the vertices (and points) are defined as an array. So a polygon with 4 vertices will be identified with an array of 8 elements with the first two elements keeping the 2D coordinates of the first vertice and so on.

Macro: GAL_POLYGON_MAX_CORNERS

The largest number of vertices a polygon can have in this library.

Macro: GAL_POLYGON_ROUND_ERR

We have to consider floating point round-off errors when dealing with polygons. for example, we will take A as the maximum of A and B when A>B-GAL_POLYGON_ROUND_ERR.

Function:
void
gal_polygon_vertices_sort_convex (double *in, size_t n, size_t *ordinds)

We have a simple polygon (that can result from projection, so its edges do not collide or it does not have holes) and we want to order its corners in an anticlockwise fashion. This is necessary for clipping it and finding its area later. The input vertices can have practically any order.

The input (in) is an array containing the coordinates (two values) of each vertice. n is the number of corners. So in should have 2*n elements. The output (ordinds) is an array with n elements specifying the indices in order. This array must have been allocated before calling this function. The indexes are output for more generic usage, for example, in a homographic transform (necessary in warping an image, see Linear warping basics), the necessary order of vertices is the same for all the pixels. In other words, only the positions of the vertices change, not the way they need to be ordered. Therefore, this function would only be necessary once.

As a summary, the input is unchanged, only n values will be put in the ordinds array. Such that calling the input coordinates in the following fashion will give an anti-clockwise order when there are 4 vertices:

1st vertice: in[ordinds[0]*2], in[ordinds[0]*2+1]
2nd vertice: in[ordinds[1]*2], in[ordinds[1]*2+1]
3rd vertice: in[ordinds[2]*2], in[ordinds[2]*2+1]
4th vertice: in[ordinds[3]*2], in[ordinds[3]*2+1]

The implementation of this is very similar to the Graham scan in finding the Convex Hull. However, in projection we will never have a concave polygon (the left condition below, where this algorithm will get to E before D), we will always have a convex polygon (right case) or E will not exist! This is because we are always going to be calculating the area of the overlap between a quadrilateral and the pixel grid or the quadrilateral itself.

The GAL_POLYGON_MAX_CORNERS macro is defined so there will be no need to allocate these temporary arrays separately. Since we are dealing with pixels, the polygon cannot really have too many vertices.

Function:
int
gal_polygon_is_convex (double *v, size_t n)

Returns 1 if the polygon is convex with vertices defined by v and 0 if it is a concave polygon. Note that the vertices of the polygon should be sorted in an anti-clockwise manner.

Function:
double
gal_polygon_area (double *v, size_t n)

Find the area of a polygon with vertices defined in v. v points to an array of doubles which keep the positions of the vertices such that v[0] and v[1] are the positions of the first vertice to be considered.

Function:
int
gal_polygon_is_inside (double *v, double *p, size_t n)

Returns 0 if point p in inside a polygon, either convex or concave. The vertices of the polygon are defined by v and 0 otherwise, they have to be ordered in an anti-clockwise manner. This function uses the winding number algorithm, to check the points. Note that this is a generic function (working on both concave and convex polygons, so if you know before-hand that your polygon is convex, it is much more efficient to use gal_polygon_is_inside_convex.

Function:
int
gal_polygon_is_inside_convex (double *v, double *p, size_t n)

Return 1 if the point p is within the polygon whose vertices are defined by v. The polygon is assumed to be convex, for a more generic function that deals with concave and convex polygons, see gal_polygon_is_inside. Note that the vertices of the polygon have to be sorted in an anti-clock-wise manner.

Function:
int
gal_polygon_ppropin (double *v, double *p, size_t n)

Similar to gal_polygon_is_inside_convex, except that if the point p is on one of the edges of a polygon, this will return 0.

Function:
int
gal_polygon_is_counterclockwise (double *v, size_t n)

Returns 1 if the sorted polygon has a counter-clockwise orientation and 0 otherwise. This function uses the concept of “winding”, which defines the relative order in which the vertices of a polygon are listed to determine the orientation of vertices. For complex polygons (where edges, or sides, intersect), the most significant orientation is returned. In a complex polygon, when the alternative windings are equal (for example, an 8-shape) it will return 1 (as if it was counter-clockwise). Note that the polygon vertices have to be sorted before calling this function.

Function:
int
gal_polygon_to_counterclockwise (double *v, size_t n)

Arrange the vertices of the sorted polygon in place, to be in a counter-clockwise direction. If the input polygon already has a counter-clockwise direction it will not touch the input. The return value is 1 on successful execution. This function is just a wrapper over gal_polygon_is_counterclockwise, and will reverse the order of the vertices when necessary.

Function:
void
gal_polygon_clip (double *s, size_t n, double *c, size_t m, double *o, size_t *numcrn)

Clip (find the overlap of) two polygons. This function uses the Sutherland-Hodgman polygon clipping algorithm. Note that the vertices of both polygons have to be sorted in an anti-clock-wise manner.

The Pseudocode from Wikipedia:

List outputList = subjectPolygon;
for (Edge clipEdge in clipPolygon) do
  List inputList = outputList;
  outputList.clear();
  Point S = inputList.last;
  for (Point E in inputList) do
     if (E inside clipEdge) then
        if (S not inside clipEdge) then
           outputList.add(ComputeIntersection(S,E,clipEdge));
        end if
        outputList.add(E);
     else if (S inside clipEdge) then
        outputList.add(ComputeIntersection(S,E,clipEdge));
     end if
     S = E;
  done
done

The difference is that we are not using lists, but arrays to keep polygon vertices. The two polygons are called Subject s and Clip c with n and m vertices respectively. The output is stored in o and the number of elements in the output are stored in what *numcrn (for number of corners) points to.

Function:
void
gal_polygon_vertices_sort (double *vertices, size_t n, size_t *ordinds)

Sort the indices of the un-ordered vertices array to a counter-clockwise polygon in the already allocated space of ordinds. It is assumed that there are n vertices, and thus that vertices contains 2*n elements where the two coordinates of the first vertice occupy the first two elements of the array and so on.

The polygon can be both concave and convex (see the start of this section). However, note that for concave polygons there is no unique sort from an un-ordered set of vertices. So after this function you may want to use gal_polygon_is_convex and print a warning to check the output if the polygon was concave.

Note that the contents of the vertices array are left untouched by this function. If you want to write the ordered vertice coordinates in another array with the same size, you can use a loop like this:

for(i=0;i<n;++i)
{
  ordered[i*2  ] = vertices[ ordinds[i]*2    ];
  ordered[i*2+1] = vertices[ ordinds[i]*2 + 1];
}

In this algorithm, we find the rightmost and leftmost points (based on their x-coordinate) and use the diagonal vector between those points to group the points in arrays based on their position with respect to this vector. For anticlockwise sorting, all the points below the vector are sorted by their ascending x-coordinates and points above the vector are sorted in decreasing order using qsort. Finally, both these arrays are merged together to get the final sorted array of points, from which the points are indexed into the ordinds using linear search.


Next: Qsort functions (qsort.h), Previous: Bounding box (box.h), Up: Gnuastro library   [Contents][Index]