Introduction

Good C++ computational geometry libraries to date have been hideously over-designed and
incorporate usage patterns that in most cases require extensive redesigns and rewrites
of code in order to functionally integrate within an existing project.

Sometimes a lightweight portable solution that has a bit of error is deemed to
be more appropriate and reasonable. However these libraries even though being
more than able to cater for such a requirement still burden the end user with
undue code clutter, very steep learning curves and in some cases unnecessary
overheads.

The solution to such a situation is to simplify use, implementation and application.
This can be achieved by reducing the number of contact points between the computational
geometry back-end and the developer-application combo. But at the same time giving
full control of the computations and ensuing folding processes to the user to the
extent where the user can decide between using a general solution for a particular
problem or a more specialized solution.

For example, say you have two line segments and you want to know if they are intersecting,
one could use a general solution, but say if you had prior knowledge that the line segments
were either always going to be vertical or horizontal, this would allow one to use a more
efficient method to obtain the same result. Typically (but not always) a generalized result
is less efficient than a specialised result for the specialised case - by virtue of the
fact that the generalised result has to take into account the 1001 other possible
scenarios.

A possible solution to the above mentioned problem is Wykobi.
Wykobi is an efficient, robust and simple to use multi-platform 2D/3D computational geometry
library. Wykobi provides a concise, predictable, and deterministic interface for geometric
primitives and complex geometric routines using and conforming to the ISO/IEC 14882:2003 C++
language specification.

The design and structure of Wykobi lends itself to easy and seamless integration into
projects of any scale that require a robust yet efficient 2D/3D computational geometry
back-end.

Wykobi as a library can be used to efficiently and seamlessly solve complex geometric problems
such as collision and proximity detection, efficient spatial queries and geometric constructions
used in areas as diverse as gaming, computer aided design and manufacture, electronic design
and geographic information systems - just to name a few.

Wykobi provides a series of primitive geometric structures for use within the various
algorithms of interest such as intersections, distances, inclusions and clipping operations.

The Wykobi Data Structures

The Point Type

Basic point types, which are zero dimensional entities that exist
in either 2D, 3D or n-dimensions.

The Ray Type

Ray type, A directed half-infinite line or half-line. A ray has an origin
point and a vector that describes the direction in which all the points
that are members of the set of points that make up the ray exist upon.

The Polygon Type

Using the code

There are many different things that can be done with the Wykobi Computational
Geometry Library. The following are some of the slightly more interesting capabilities...

Calculating A Convex Hull

The convex hull of a set of points, is the subset of points from the original
set that comprise a convex shaped polygon or polytope which bounds all the points
in the original set.

Many different techniques exist for calculating the convex hull of a set of points.
Various methods such as the Melkman algorithm rely on special properties of the
points. Complexities for calculating the convex hull range from naive algorithms
which have a complexity of O(N3) to more specialised algorithms such Graham scan
and Melkman that have complexities of O(nlogn) and O(n) respectively.

Calculating A Minimum Bounding Ball

Given a set of n k-dimensional points, the minimum bounding ball is the smallest circle,
sphere or hypersphere that contains all the points. This problem is sometimes called the
smallest enclosing circle or the smallest enclosing disk where by the points in contention
must all be coplanar to each other.

Randomized Algorithm

The randomized algorithm is a stable algorithm which is used to solve
the minimum bounding ball problem for 2D with a space and time complexity O(n).

Ritter Algorithm

An approximation algorithm devised by Jack Ritter [ritter 1990]. It has a complexity
of O(n), can be easily extended to higher dimensions yet does not guarantee an optimal
minimum bounding ball, just something very close.

Note:

All the 2D minimum bounding ball algorithms have been extended to perform a
convex hull filter operation before calculating the bounding ball. Even though
obtaining the convex hull is not of linear complexity, the resulting points from
the hull guarantee a somewhat better result with regards to the optimal minimum
bounding ball when considering the ritter algorithm, When considering the naive
algorithm there is a large linear scaling down of computing time though not of
the complexity, and finally when considering the random algorithm if the point
set is large enough, preprocessing the set by computing its convex hull decreases
(through a constant multiplier to its complexity )the computing time that is
incurred during the randomized rotation process carried out in each step - essentially
in all cases the less points there are the more efficient the algorithms become.

Sutherland Hodgman Polygon Clipping

Clipping one object or more precisely a polygon or polytope against another is essentially the process of
computing the intersecting area or volume between the pair of objects. Depending on the structural nature
of the objects such as convexity and disjointness, the resulting clipped object may itself be disjoint or
may contain islands and other interesting properties.

The Sutherland Hodgman polygon clipping algorithm is a simplified clipping algorithm with the constraint
that the clip boundary be convex where as the other object may be a concave non-self intersecting polygon.

Calculating The Axis Projection Descriptor

In the book flatland and subsequent flatterland, the flatlanders would query an
object's boundary to determine its identity, similar objects would have
similar boundaries. The projection of a 2D object onto various axises produces a
view of that object on those axises. The combinations of views are somewhat unique
to that object and using various normalization methods and difference metrics can
be used to define, in a somewhat invariant to rotation and scaling manner, how
different or how similar one 2D object is from another. The descriptor works best
with convex shapes. When comparing concave shapes, due to the possibility that edges
may occlude others, a definitive differencing metric is difficult to realize.
The concepts used in calculating the descriptor are very similar to the concepts
used in the separating axis theorem.

The 3D form of the axis projection descriptor involves projecting a voluminous object
on various planes. The areas of the projections as in the 2D sense are normalized
producing a histogram covering the planes. A differencing metric such as the bhattacharya
distance can then be used to efficiently determine relative equivalency between two
or more objects.

Simple Examples

The Orientation Predicate

In the world of computational geometry there are many predicates that form
the basis of some of the most complex calculations known. One of the most
important predicates is known as the orientation predicate. The question this
predicate tries to answer is that, given two points P0 and P1
that compose a directed line formed by the vector V (P1 - P0),
and a third point P2, on what side of the directed line does the
point P2 reside?

The concept of a directed line leads to the definition of three unique spaces.
The first two represent the set of points to the left and right of the directed
line whereas the third is the set of points that comprise the directed line.
Another way to view the problem is to assume one is standing on P0
and looking towards P1, the question as defined before is on what
side does P2 reside?

The diagram below depicts the point P2 as being on the left-hand side of
the defined directed line.

The diagram below depicts all three possibilities of the orientation predicate.
Wykobi provides an implementation of the orientation routine with
overloads for line-segments and lines (directed-line representation).
The routine returns one of either values: RightHandSide,
LeftHandSide or CollinearOrientation

Note: The real underlying computation that occurs within the orientation
predicate is called the signed area of a triangle or the 3-point determinant.
The orientation predicate can be extended to 3D where by the testing object
is a plane rather than a directed-line constructed from three points, and
where by a fourth point is queried as being either above, below or coplanar
to the defined plane. In 3D the signed area of the triangle is equivalent to
the signed area of the tetrahedron. There has been extensive research done on
how best to compute the orientation predicate. The sign of the area is the most
important aspect of the computation and getting that right using finite precision
arithmetic seems to be an ongoing research topic.

Segment Intersection Via The Orientation Predicate

The orientation predicate can be used to solve some very interesting problems. Below we have
two line-segments, S0 and S1, comprised of point pairs (P0,P1)
and (Q0,Q1) respectively. The question we would like to solve is whether
or not the two line-segments intersect with each other.

By applying the orientation rule, with the directed-line constructed from S1
as D0 = Q0 - Q1, and using the end points from
S0 P0 and P1, we discover a very interesting
property, that is if S0 intersects D0 then the orientations
for the end points of S0 will not be equivalent.

As demonstrated in the diagram below in the event the given line-segments do not
intersect the computed orientations to the directed-lines D0 and
D1 will be equivalent.

There is a caveat to this rule and that is if the either one or both of the points are collinear
to their respective test directed-line and exist within the axis-aligned bounding box defined by
the end-points used to construct said directed-line, then the line-segments intersect.

Further more the situation below shows the need for the above described operation to
occur for both segments. An orientation test for the end-points of S0 against
the directed-line D0 constructed from S1 and an orientation
test for the end-points of S1 against the directed-line D1
constructed from S0.

Wykobi provides a routine that does the above namely simple_intersect, The
algorithm's pesduo-code is as follows:

if (orientation(Q<sub>0</sub>,Q<sub>1</sub>,P<sub>0</sub>) to orientation(Q<sub>0</sub>,Q<sub>1</sub>,P<sub>1</sub>) is either not equal or collinear) and
(orientation(P<sub>0</sub>,P<sub>1</sub>,Q<sub>0</sub>) to orientation(P<sub>0</sub>,P<sub>1</sub>,Q<sub>1</sub>) is either not equal or collinear) then
begin
return segments intersect.
end

Note: There are some issues relating to this method, initially the algorithm only returns the
boolean state of intersection between the provided segments, not the intersection point itself.
Also the issues relating to numerical stability of the orientation predicate may cause this routine
to provide inconsistent results when given the same input but in different orders - this only
occurs when very large values and small values are mixed and is not necessarily the algorithm's
problem but more so due to the use of finite precision numbers which are generally the main source
of problems when performing numerical computations - after buggy code of course.

Point In Convex Polygon Via The Orientation Predicate

Following on from the previous theme of using the orientation predicate, we have another
problem at hand, given a convex polygon and point, determine if the point is within the
polygon.

As before another interesting property is discovered when using the orientation predicate,
and that is if one were to traverse the edges of the convex polygon in a consistent order
computing the orientation between the given point and the directed-line constructed from the
edge at-hand, then one could say the point exists within the convex polygon if and only if
all the orientations are equivalent.

The diagram below demonstrates that edges |P1,P0|, |P0,P4|,
|P4,P3| and |P2,P1| all have the same orientation
but the edges |P4,P3| and |P3,P2| have differing orientations
to the point Q0. From this it is determined that the point Q0 is outside
the convex polygon.

Again as before there is a caveat to this rule and that is if any of the orientations are collinear
then that orientation can be considered the same as the others if and only if the point is within
the axis aligned bounding box of the edge at-hand. This principle can also be used to easily
determine if a polygon is convex or not. Wykobi provides a set of routines that use the above
mentioned principle in various ways, they include point_in_convex_polygon,
convex_quadix and is_convex_polygon.

The Point Of Reflection

Ray-tracing processes primarily concern themselves with shooting out rays from a point
of view and tracing the rays through a scene as they intersect, diffuse and reflect
various surfaces within the scene. Another problem that arises from this area is the
definition of points of reflection upon surfaces that satisfy a constraint, for example
a billiard table simulation may require the solution to determining the point on the
side of a table where the white ball must hit so as to reflect off-of and hit a target
ball. Another question that could be asked is if the target point is visible from the
source point via a reflection from an object within the scene.

In the diagram below we have two points of interest a source and destination point,
blue and green respectively, and a mirror-like object represented as a black
line-segment. Lets assume this object provides totally elastic collisions and
energy-loss free reflections from object or light interactions.

The question then is where abouts on the reflective line-segment should a ray
from the blue source object be targeted at so as to reflect (angle of incidence
is equal to angle of reflection) and hit the green target object?

It turns out there is a very simple geometric construction that provides a
closed-form solution to this problem. The first step is to extend perpendiculars
from the target and source points onto the line extension of the reflective
line-segment. Or in other words determine the closest points on the line extension
from the target and source points.

The next step is to extend line-segments from the target and source points to the opposite point's
perpendicular extensions. These two line segments should result in an intersection point. At this
point we extend a perpendicular from the intersection point onto the line extension of the reflective
line-segment.

If the intersection point's perpendicular extension is upon the line-segment
(not the line extension), then that point is said to be the point of reflection, the
red point in the diagram below. Casting a ray from the blue point towards the red
point will cause the ray to reflect at that point towards to green point.

Wykobi provides a routine specifically for this kind of computation point_of_reflection

Note: It is assumed that both the target and source points are of the same orientation with regards
to the reflective line-segment. Furthermore there are some very interesting triangle similarity
properties of the above mentioned construction. Can you determine what they are?

Bisecting An Angle (The Other Method...)

Previously it was shown that computing the bisector of an angle with Wykobi can be done very
easily. The basis of the computation was to exploit certain angle and side similarities to
deduce the bisector. Now we have a similar geometric construction problem but with the following
restrictions, Assume you are given an angle constructed from two lines as shown below and only
using a compass (not a protractor) and a non-etched ruler you are to determine the bisector for
the specified angle |ABC|.

We initially begin by drawing a circle centered at B with sufficient radius. We shall call the intersection
points invwards of the angle |ABC| with lines |AB| and |CB| and the circle, N0 and N1
respectively. Furthermore we draw two circles with centers at N0 and N1 with radii such
that they intersect each other at least at one other point.

In the diagram below it so happens that the two circles centered at N0 and N1 intersect
each other at two locations I0 and I1. From this we determine the line that passes
through B, I0 and I1 is the line which bisects the angle |ABC|.

The above mentioned computation can be easily achieved using Wykobi as follows:

Smoothing Of Sharp Corners

There are many reasons why one would need to smooth-out or unsharpen the corner of
an object. Below we'll describe a very simple method to non-uniformly add superfluous
edges around a corner.

Initially we project the end points of each edge outwards by length K in a
direction perpendicular to the edge being processed. The new points will generate
one new edge.

A third edge is generated by using the projected points that arise from the corner
(or common source), this edge is inserted between the two previously generated edges.

The process is repeated over again with the new edges generating a new
set of edges that are then fed back into the algorithm, resulting in a
some-what smoothed-out corner. The following is a solution using Wykobi
to the above mentioned process:

The diagram below demonstrates how a triangle converted to a polygon type has had
its sharp corners smoothed-out. The size of the new boundary being larger is purely
for demonstration purposes. In practice the new boundary would be nearly superimposed
on top of the original.

The above is a very brief overview of only some of the computational geometry algorithms
and processes that are available in the Wykobi C++ computational geometry library. There
is plenty more out there...

Some Notes On Usage

1. Type

The above examples demonstrate the generality/genericity of the Wykobi library
routines with regards to numerical type. However this may be somewhat misleading
as not all types can provide the necessary precision required to obtain satisfactory
results from the routines. Consequently one must approach problems with at least
some information relating to bounds and required precision and efficiency. This is
a problem that a library can never solve but rather provide the end developer the
tools and options by which they can make the necessary decisions to solve their
problem.

2. Robustness

Wykobi's routines make assumptions about the validity of types being passed to them.
Typically these assumptions are manifest by the lack of assertions and type degeneracy
checks within the routines themselves. This is done so as to provide the most optimal
implementation of the routine without causing the routine to fail, and to leave the
details of type validation to the end user as they see fit. Theoretically each of the
routines could verify object degeneracy (e.g: does the triangle have 3 unique points),
then type value validity (e.g: does the value lie within some plausible range) but the
unnecessary overhead one must endure would make using the routines quite inefficient.
As an example consider what the circumcircle of a triangle that has all 3 of its points
being collinear would look like, how would you write the routine to be robust, when would
you need to have a robust routine like that?

2. Correctness

Typical usage patterns involve chaining the output of one routine as the input of
another so on and so forth. Not knowing the exact nature of the computation will
lead to an aggregation of errors that might result in the final outcome being highly
erroneous and subsequently unusable. An example of this is as follows, assume you
have an arm of length x with one end statically positioned at the origin, requests
for rotations of the arm come through, in degree form, +1, -13.5, +290 etc.

This particular form of computation has a number of errors that will result in the arm's
final position of its other end not being accurately calculated if the computation
is carried out on every request. The errors that can occur range from the conversion
of degrees to radian and the way trigonometric identities are calculated in the
processor and of course the application of the rotation transform to the arm itself.
A simple solution to this accumulating error would be to sum-up a certain amount of
rotation requests, and then to apply the final sum modulo 360 to the routine. This
would significantly reduce the total error and as a result would provide some degree
of correctness for the routine.

Wykobi .NET Demo

A very simple WinForms based application is available from the Wykobi download
site which demonstrates the above code pieces. It should be noted there is a massive speed hit (roughly 75-80%) when working under the managed
environment as when compared to native compilation (YMMV).

History

The Wykobi C++ computational geometry library is currently at version 0.0.4. The most recent
versions and updates can be obtained from the official Wykobi site www.wykobi.com

The code available with this article has a VCL graphics file named 'wykobi_graphics_vcl.hpp' for use with the IDE C++Builder from Embarcadero Technology. Are there any instructions to compile these files into a C++Builder project? Some of the files are the VC file type .inl which I am not sure how to handle.

Based on the problem described, it seems to require Boolean operations over polylines - with some kind of intelligent alpha-hull for cleanup.

For the 2D side of things, all the building block routines are present in Wykobi. The general operations however are not - I would suggest approximating all your objects into polylines then using GPC for the Boolean subtraction from your test object. As for alpha-hull for cleaning up the result - that is pretty much beyond the scope of Wykobi.

Thank you for the great article. Is there any possibility to use the library for partitioning set of points into clusters? Cluser would mean a subset of points, which distance to center of cluster is lower than distance to center of any other cluster.

<layer> Is there any possibility to use the library for partitioning set of points into clusters? Cluser would mean a subset of points, which distance to center of cluster is lower than distance to center of any other cluster.

Sounds like you're looking for something similar to k-means or Canopy clustering. At the moment Wykobi doesn't have those capabilities, it only has some of the utility functions that would be required when performing some computations.

That said in the past I've used OpenCVs[^] k-means implementation which was quite good.