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_flat (double *v, size_t n) ¶

Find the area of a polygon with vertices defined in v on a euclidian (flat) coordinate system. 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:
double
gal_polygon_area_sky (double *v, size_t n) ¶

Find the area of a polygon with vertices defined in v on a celestial coordinate system. This is a coordinate system where the first coordinate goes from 0 to 360 (increasing to the right), while the second coordinate ranges from -90 to +90 (on the poles). 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.

This function uses an approximation to account for the curvature of the sky and the different nature of spherical coordinates with respect to the flat coordinate system. Bug 64617 has been defined in Gnuastro to address this problem. Please check that bug in case it has been fixed. Until this bug is fixed, here are some tips:

• Subtract the RA and Dec of all the vertice coordinates from a constant so the center of the polygon falls on (RA, Dec) of (180,0). The sphere has a similar nature everywhere on it, so shifting the polygon vertices will not change its area; this also removes issues with the RA=0 or RA=360 coordinate and decrease issues caused by RA depending on declination.
• These approximations should not cause any statistically significant error on normal (less than a few degrees) scales. But it won’t hurt to do a small sanity check for your particular usage scenario.
• Any help (even in the mathematics of the problem; not necessary programming) would be appreciated (we didn’t have time to derive the necessary equations), so if you have some background in this and can prepare the mathematical description of the problem, please get in touch.

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.