Wykobi - Computational Geometry Tutorial

Wykobi is a lightweight and simple to use C++ Computational Geometry Library. The library focuses primarily on 2D and 3D based geometric problems, though it does have support for some N-D versions of those same problems. The following contains a more extensive listing of all the available features



Wykobi as a library can be used to efficiently and seamlessly solve 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 following is a list of the supported primitives:

Point (2D, 3D, ND)

Line

Segment

Ray

Triangle

Rectangle

Quadix (Quadrilateral)

Polygon

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

template <typename T = Float> class point2d : public geometric_entity {}; template <typename T = Float> class point3d : public geometric_entity {}; template <typename T = Float, std::size_t Dimension> class pointnd : public geometric_entity {};

Line type, which is a 1 dimensional entity of infinite length that is described by two points within its present dimension.

template <typename T = Float, std::size_t Dimension> class line : public geometric_entity {};

Segment type, similar to the line type, but is of finite length bounded by the two points which describe it within its present dimension.

template <typename T = Float, std::size_t Dimension> class segment : public geometric_entity {};

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.

template <typename T = Float, std::size_t Dimension> class ray : public geometric_entity {};

Triangle type, A geometric primitive that is comprised of 3 unique points, which produce 3 unique edges.

template <typename T = Float, std::size_t Dimension> class triangle : public geometric_entity {};

Rectangle type, An axis aligned 4 sided geometric primitive, described by two bounding points in 2D. A rectangle's form in 3D and higher dimensions is a box.

template <typename T = Float, std::size_t Dimension> class rectangle : public geometric_entity {};

Quadix type, A convex quadrilateral or polygon that comprises of 4 unique points which produce 4 unique edges. In the 3D and higher dimensions sense all 4 points have to be coplanar.

template <typename T = Float, std::size_t Dimension> class quadix : public geometric_entity {};

Polygon type, A set of closed sequentially connected coplanar points.

template <typename T = Float, std::size_t Dimension> class polygon : public geometric_entity {};

The polygon type in Wykobi supports only a limited set of possible polygon types. As such the following polygon traits and features are not supported:

Disjoint areas

Holes

Islands

Overlapping areas

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...

The convex hull of a set of points, is the subset of points from the original set that comprise the smallest possible 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.

const std::size_t max_points = 100000; std::vector<wykobi::point2d<T>> point_list; point_list.reserve(max_points); wykobi::generate_random_points<T> ( 0.0, 0.0, width, height, max_points, std::back_inserter(point_list) ); wykobi::polygon<T,2> convex_hull; wykobi::algorithm::convex_hull_graham_scan<wykobi::point2d<T>> ( point_list.begin(), point_list.end(), std::back_inserter(convex_hull) );

The Jarvis march algorithm is also known as the gift-wrapping algorithm. It can be naturally extended to higher dimensions.

wykobi::polygon<T,2> polygon; generate_polygon_type2<T>(width,height,polygon); std::vector<wykobi::point2d<T>> convex_hull; wykobi::algorithm::convex_hull_jarvis_march<wykobi::point2d<T>> ( polygon.begin(), polygon.end(), std::back_inserter(convex_hull) ); wykobi::polygon<T,2> convex_hull_polygon = wykobi::make_polygon<T>(convex_hull);

The Melkman algorithm achieves a complexity of O(n) by assuming that the points in the set are ordered such that they represent a concave non-selfintersecting polygon or polyline.

wykobi::polygon<T,2> polygon; generate_polygon_type2<T>(width,height,polygon); std::vector<wykobi::point2d<T>> convex_hull; wykobi::algorithm::convex_hull_melkman<wykobi::point2d<T>> ( polygon.begin(), polygon.end(), std::back_inserter(convex_hull) ); wykobi::polygon<T,2> convex_hull_polygon = wykobi::make_polygon<T>(convex_hull);

Given a set of N k-dimensional points, the minimal 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.

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

const std::size_t max_points = 100000; std::vector<wykobi::point2d<T>> point_list; point_list.reserve(max_points); wykobi::generate_random_points<T> ( 0.0, 0.0, width, height, max_points, std::back_inserter(point_list) ); wykobi::circle<T> minimum_bounding_ball; wykobi::algorithm::randomized_minimum_bounding_ball<wykobi::point2d<T>> ( point_list.begin(), point_list.end(), minimum_bounding_ball );

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 minimal bounding ball, but rather a very close approximate.

const std::size_t max_points = 100000; std::vector<wykobi::point2d<T>> point_list; point_list.reserve(max_points); wykobi::generate_random_points<T> ( 0.0, 0.0, width, height, max_points, std::back_inserter(point_list) ); wykobi::circle<T> minimum_bounding_ball; wykobi::algorithm::ritter_minimum_bounding_ball<wykobi::point2d<T>> ( point_list.begin(), point_list.end(), minimum_bounding_ball );

const std::size_t max_points = 100000; std::vector<wykobi::point2d<T>> point_list; point_list.reserve(max_points); wykobi::generate_random_points<T> ( 0.0, 0.0, width, height, max_points, std::back_inserter(point_list) ); wykobi::circle<T> minimum_bounding_ball; wykobi::algorithm::naive_minimum_bounding_ball<wykobi::point2d<T>> ( point_list.begin(), point_list.end(), minimum_bounding_ball );

Note: All the 2D minimal bounding ball algorithms have been extended to perform a convex hull filter operation before calculating the bounding ball.



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 minimal 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.

wykobi::algorithm::randomized_minimum_bounding_ball_with_ch_filter<wykobi::point2d<T>> ( point_list.begin(), point_list.end(), minimum_bounding_ball ); wykobi::algorithm::ritter_minimum_bounding_ball_with_ch_filter<wykobi::point2d<T>> ( point_list.begin(), point_list.end(), minimum_bounding_ball ); wykobi::algorithm::naive_minimum_bounding_ball_with_ch_filter<wykobi::point2d<T>> ( point_list.begin(), point_list.end(), minimum_bounding_ball );

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.

wykobi::rectangle<T> clip_boundry; wykobi::generate_random_object<T> (0.0, 0.0, width - 1.0, height - 1.0, clip_boundry); wykobi::polygon<T,2> polygon; generate_polygon_type2<T>(width,height,polygon); wykobi::polygon<T,2> clipped_polygon; wykobi::algorithm::sutherland_hodgman_polygon_clipper<wykobi::point2d<T>> (clip_boundry, polygon, clipped_polygon);

wykobi::triangle<T> clip_boundry; wykobi::generate_random_object<T> (0.0, 0.0, width - 1.0, height - 1.0, clip_boundry); wykobi::polygon<T,2> polygon; generate_polygon_type2<T>(width,height,polygon); wykobi::polygon<T,2> clipped_polygon; wykobi::algorithm::sutherland_hodgman_polygon_clipper<wykobi::point2d<T>> (clip_boundry, polygon, clipped_polygon);

wykobi::quadix<T> clip_boundry; wykobi::generate_random_object<T> (0.0, 0.0, width - 1.0, height - 1.0, clip_boundry); wykobi::polygon<T,2> polygon; generate_polygon_type2<T>(width,height,polygon); wykobi::polygon<T,2> clipped_polygon; wykobi::algorithm::sutherland_hodgman_polygon_clipper<wykobi::point2d<T>> (clip_boundry, polygon, clipped_polygon);

wykobi::circle<T> circle; wykobi::generate_random_object<T>(0.0, 0.0, width, height, circle); wykobi::polygon<T,2> clip_boundry = wykobi::make_polygon<T>(circle,12); wykobi::generate_random_object<T> (0.0, 0.0, width - 1.0, height - 1.0, circle); wykobi::polygon<T,2> polygon; generate_polygon_type2<T>(width,height,polygon); wykobi::polygon<T,2> clipped_polygon; wykobi::algorithm::sutherland_hodgman_polygon_clipper<wykobi::point2d<T>> (clip_boundry, polygon, clipped_polygon);

const std::size_t max_segments = 100; std::vector<wykobi::segment<T,2>> segment_list; for (std::size_t i = 0; i < max_segments; ++i) { wykobi::segment<T,2> tmp_segment; wykobi::generate_random_object(0.0, 0.0, width, height, tmp_segment); segment_list.push_back(tmp_segment); } wykobi::rectangle<T> rectangle; wykobi::generate_random_object<T>(0.0, 0.0, width, height, rectangle); std::vector<wykobi::segment<T,2>> clipped_segment_list; for (std::size_t i = 0; i < segment_list.size(); ++i) { wykobi::segment<T,2> clipped_segment; if (wykobi::clip(segment_list[i], rectangle, clipped_segment)) { clipped_segment_list.push_back(clipped_segment); } }

const std::size_t max_segments = 100; std::vector<wykobi::segment<T,2>> segment_list; for (std::size_t i = 0; i < max_segments; ++i) { wykobi::segment<T,2> tmp_segment; wykobi::generate_random_object(0.0, 0.0, width, height, tmp_segment); segment_list.push_back(tmp_segment); } wykobi::triangle<T,2> triangle; wykobi::generate_random_object<T>(0.0, 0.0, width, height, triangle); std::vector<wykobi::segment<T,2>> clipped_segment_list; for (std::size_t i = 0; i < segment_list.size(); ++i) { wykobi::segment<T,2> clipped_segment; if (wykobi::clip(segment_list[i], triangle, clipped_segment)) { clipped_segment_list.push_back(clipped_segment); } }

const std::size_t max_segments = 100; std::vector<wykobi::segment<T,2>> segment_list; for (std::size_t i = 0; i < max_segments; ++i) { wykobi::segment<T,2> tmp_segment; wykobi::generate_random_object(0.0, 0.0, width, height, tmp_segment); segment_list.push_back(tmp_segment); } wykobi::quadix<T,2> quadix; wykobi::generate_random_object<T>(0.0, 0.0, width, height, quadix); std::vector<wykobi::segment<T,2>> clipped_segment_list; for (std::size_t i = 0; i < segment_list.size(); ++i) { wykobi::segment<T,2> clipped_segment; if (wykobi::clip(segment_list[i], quadix, clipped_segment)) { clipped_segment_list.push_back(clipped_segment); } }

const std::size_t max_segments = 100; std::vector<wykobi::segment<T,2>> segment_list; for (std::size_t i = 0; i < max_segments; ++i) { wykobi::segment<T,2> tmp_segment; wykobi::generate_random_object(0.0, 0.0, width, height, tmp_segment); segment_list.push_back(tmp_segment); } wykobi::circle<T> circle; wykobi::generate_random_object<T>(0.0, 0.0, width, height, circle); std::vector<wykobi::segment<T,2>> clipped_segment_list; for (std::size_t i = 0; i < segment_list.size(); ++i) { wykobi::segment<T,2> clipped_segment; if (wykobi::clip(segment_list[i], circle, clipped_segment)) { clipped_segment_list.push_back(clipped_segment); } }

const std::size_t max_rectangles = 10; std::vector<wykobi::rectangle<T>> rectangle_list; for (std::size_t i = 0; i < max_rectangles; ++i) { wykobi::rectangle<T> tmp_rectangle; wykobi::generate_random_object(0.0,0.0,width,height,tmp_rectangle); rectangle_list.push_back(tmp_rectangle); } wykobi::rectangle<T> clip_rectangle; wykobi::generate_random_object<T>(0.0, 0.0, 1000.0, 1000.0, clip_rectangle); std::vector<wykobi::rectangle<T>> clipped_rectangle_list; for (std::size_t i = 0; i < rectangle_list.size(); ++i) { wykobi::rectangle<T> clipped_rectangle; if (wykobi::clip(rectangle_list[i], clip_rectangle, clipped_rectangle)) { clipped_rectangle_list.push_back(clipped_rectangle); } }

const std::size_t max_segments = 100; std::vector<wykobi::segment<T,2>> segment_list; for (std::size_t i = 0; i < max_segments; ++i) { wykobi::segment<T,2> tmp_segment; wykobi::generate_random_object(0.0, 0.0, width, height, tmp_segment); segment_list.push_back(tmp_segment); } std::vector<wykobi::point2d<T>> intersection_list; wykobi::algorithm::naive_group_intersections<segment<T,2>> ( segment_list.begin(), segment_list.end(), std::back_inserter(intersection_list) );

const std::size_t max_circles = 100; std::vector<wykobi::circle<T>> circle_list; for (std::size_t i = 0; i < max_circles; ++i) { wykobi::circle<T> tmp_circle; wykobi::generate_random_object(0.0, 0.0, width, height, tmp_circle); circle_list.push_back(tmp_circle); } std::vector<wykobi::point2d<T>> intersection_list; wykobi::algorithm::naive_group_intersections<circle<T>> ( circle_list.begin(), circle_list.end(), std::back_inserter(intersection_list) );

There are three primary modes of circle-to-circle boundary intersections. These range from no intersection points up to two intersection points. Lets assume D is the distance between the centers C 0 and C 1 , then the modes can be described in terms of inequalities between the radii of circles (R 0 , R 1 ) and D.

D > (R 0 + R 1 ) - No intersections

- No intersections D = (R 0 + R 1 ) - Single intersection

- Single intersection D < (R 0 + R 1 ) - Two intersections

- Two intersections D < |R 0 - R 1 | - No intersections Complete overlap

The following diagram illustrates the various intersection modes:

The following diagram is a breakdown of all the points and segments of interest when determining the intersection points between two circles where the distance between the centers is less than the sum and greater than the difference of the radii:

In order to determine the intersection points I 0 and I 1 , we need to first determine the location of point P. Which can be calculated in terms of a unit vector v along the path C 0 to C 1 scaled by the length d 0 originating from C 0 . Using the identities on the left, we then derive the value of d 0 using the process described on the right:

Now that we have a computable form for the length d 0 , we can proceed to derive the point P and from there determine the intersection points I 0 and I 1 as being the composition of P and the perpendicular vector to v scaled by the length k:

wykobi::polygon<T,2> polygon; polygon.push_back(wykobi::make_point<T>( 25.0, 191.0)); polygon.push_back(wykobi::make_point<T>( 55.0, 191.0)); polygon.push_back(wykobi::make_point<T>( 52.0, 146.0)); polygon.push_back(wykobi::make_point<T>( 98.0, 134.0)); polygon.push_back(wykobi::make_point<T>(137.0, 200.0)); polygon.push_back(wykobi::make_point<T>(157.0, 163.0)); polygon.push_back(wykobi::make_point<T>(251.0, 188.0)); polygon.push_back(wykobi::make_point<T>(151.0, 138.0)); polygon.push_back(wykobi::make_point<T>(164.0, 116.0)); polygon.push_back(wykobi::make_point<T>(125.0, 141.0)); polygon.push_back(wykobi::make_point<T>( 78.0, 99.0)); polygon.push_back(wykobi::make_point<T>( 29.0, 139.0)); std::vector<wykobi::triangle<T,2>> triangle_list; wykobi::algorithm::polygon_triangulate<wykobi::point2d<T>> (polygon, std::back_inserter(triangle_list));

wykobi::polygon<T,2> polygon = wykobi::make_polygon(wykobi::make_circle<T>(0.0, 0.0, 100.0),10); std::vector<wykobi::triangle<T,2>> triangle_list; wykobi::algorithm::polygon_triangulate<wykobi::point2d<T>> (polygon, std::back_inserter(triangle_list));

wykobi::polygon<T,2> polygon; generate_polygon_type1<T>(width, height, polygon); // generate simple complex polygon std::vector triangle_list; wykobi::algorithm::polygon_triangulate<wykobi::point2d<T>> (polygon, std::back_inserter(triangle_list));

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.

wykobi::quadix<T,2> quadix; wykobi::generate_random_object<T>(0.0, 0.0, width, height, quadix); wykobi::polygon<T,2> polygon = wykobi::make_polygon(quadix); std::vector<T> descriptor; wykobi::algorithm::generate_axis_projection_descriptor<T> (wykobi::polygon, std::back_inserter(descriptor));

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.

const std::size_t bezier_count = 15; const std::size_t resolution = 1000; std::vector<wykobi::quadratic_bezier<T,2>> bezier_list; for (std::size_t i = 0; i < bezier_count; ++i) { wykobi::quadratic_bezier<T,2> bezier; bezier[0] = wykobi::generate_random_point<T>(width, height); bezier[1] = wykobi::generate_random_point<T>(width, height); bezier[2] = wykobi::generate_random_point<T>(width, height); bezier_list.push_back(bezier) } for (std::size_t i = 0; i < bezier_list.size(); ++i) { std::vector<point2d<T>> point_list; wykobi::generate_bezier(bezier_list[i], resolution, point_list); draw_polyline(point_list); }

const std::size_t bezier_count = 15; const std::size_t resolution = 1000; std::vector<wykobi::cubic_bezier<T,2>> bezier_list; for (std::size_t i = 0; i < bezier_count; ++i) { wykobi::cubic_bezier<T,2> bezier; bezier[0] = wykobi::generate_random_point<T>(width, height); bezier[1] = wykobi::generate_random_point<T>(width, height); bezier[2] = wykobi::generate_random_point<T>(width, height); bezier[3] = wykobi::generate_random_point<T>(width, height); bezier_list.push_back(bezier) } for (std::size_t i = 0; i < bezier_list.size(); ++i) { std::vector<point2d<T>> point_list; wykobi::generate_bezier(bezier_list[i], resolution, point_list); draw_polyline(point_list); }

Note: The current method uses an iterative approximation approach. The correct method is to use a polynomial root solver to find the quadratic or cubic polynomial roots for every dimension.

const std::size_t bezier_count = 20; const std::size_t segment_count = 10; const std::size_t resolution = 1000; std::vector<wykobi::quadratic_bezier<T,2>> bezier_list; std::vector<wykobi::segment<T,2>> segment_list; for (std::size_t i = 0; i < bezier_count; ++i) { wykobi::quadratic_bezier<T,2> bezier; bezier[0] = wykobi::generate_random_point<T>(width, height); bezier[1] = wykobi::generate_random_point<T>(width, height); bezier[2] = wykobi::generate_random_point<T>(width, height); bezier_list.push_back(bezier) } for (std::size_t i = 0; i < segment_count; ++i) { wykobi::segment<T,2>; segment; wykobi::generate_random_object<T>(0.0, 0.0, width, height, segment); segment_list.push_back(segment); } std::vector<wykobi::point2d<T>> intersection_point_list; for (std::size_t i = 0; i < bezier_list.size(); ++i) { for (std::size_t j = 0; j < segment_list.size(); ++j) { wykobi::intersection_point (segment_list[j], bezier_list[i], intersection_point_list); } }

const std::size_t bezier_count = 20; const std::size_t segment_count = 10; const std::size_t resolution = 1000; std::vector<wykobi::cubic_bezier<T,2>> bezier_list; std::vector<wykobi::segment<T,2>>segment_list; for (std::size_t i = 0; i < bezier_count; ++i) { wykobi::cubic_bezier<T,2> bezier; bezier[0] = wykobi::generate_random_point<T>(width, height); bezier[1] = wykobi::generate_random_point<T>(width, height); bezier[2] = wykobi::generate_random_point<T>(width, height); bezier[3] = wykobi::generate_random_point<T>(width, height); bezier_list.push_back(bezier) } for (std::size_t i = 0; i < segment_count; ++i) { wykobi::segment<T,2> segment; wykobi::generate_random_object<T>(0.0, 0.0, width, height, segment); segment_list.push_back(segment); } std::vector<wykobi::point2d<T>> intersection_point_list; for (std::size_t i = 0; i < bezier_list.size(); ++i) { for (std::size_t j = 0; j < segment_list.size(); ++j) { wykobi::intersection_point (segment_list[j], bezier_list[i], intersection_point_list); } }

The routines in this section will attempt to generate uniformly distributed points within the boundries of the following 2D objects:

Rectangle

Triangle

Quadix

Circle

const std::size_t max_points = 1000; wykobi::rectangle<T> rectangle; wykobi::generate_random_object<T>(0.0, 0.0, width, height, rectangle); std::vector<wykobi::point2d<T>> point_list; point_list.reserve(max_points); wykobi::generate_random_points(rectangle, max_points, std::back_inserter(point_list));

const std::size_t max_points = 1000; wykobi::triangle<T,2> triangle; wykobi::generate_random_object<T>(0.0, 0.0, width, height, triangle); std::vector<wykobi::point2d<T>> point_list; point_list.reserve(max_points); wykobi::generate_random_points(triangle, max_points, std::back_inserter(point_list));

const std::size_t max_points = 1000; wykobi::quadix<T,2> triangle; wykobi::generate_random_object<T>(0.0, 0.0, width, height, quadix); std::vector<wykobi::point2d<T>> point_list; point_list.reserve(max_points); wykobi::generate_random_points(quadix, max_points, std::back_inserter(point_list));

const std::size_t max_points = 1000; wykobi::circle<T> circle; wykobi::generate_random_object<T>(0.0, 0.0, width, height, circle); std::vector<wykobi::point2d<T>> point_list; point_list.reserve(max_points); wykobi::generate_random_points(circle, max_points, std::back_inserter(point_list));

wykobi::point2d<T> point_a = wykobi::make_point(...,...); wykobi::point2d<T> point_b = wykobi::make_point(...,...); wykobi::point2d<T> point_c = wykobi::make_point(...,...); wykobi::line<T,2> bisector_line = wykobi::create_line_from_bisector(point_a, point_b, point_c);

The following defines the tangent lines (or line-segments) between a given circle C and an external point P. Furthermore this calculation is the basis of the Inner and Outer tangent lines routines.

wykobi::circle<T> circle = wykobi::make_circle<T>(0.0, 0.0, 100.0); wykobi::point2d<T> external_point = wykobi::make_point<T>(-1000.0, 1000.0); wykobi::point2d<T> tangent_point1; wykobi::point2d<T> tangent_point2; wykobi::circle_tangent_points(circle, external_point, tangent_point1, tangent_point2);

Note: The length of the line segments defined between point P and the tangent points T 0 and T 1 are equal, let's call this length K. The value of K is equal to: sqrt(r2 + |C - P|2). As such the problem of finding the tangent points on the circle C simplifies to that of finding the points of intersection between two circles centered at points C and P, with radii of R and K respectively.

wykobi::circle<T> circle0 = wykobi::make_circle<T>(......); wykobi::circle<T> circle1 = wykobi::make_circle<T>(......); std::vector<wykobi::line<T,2>> tangent_lines; wykobi::circle_internal_tangent_lines(circle0,circle1,tangent_lines);

wykobi::circle<T> circle0 = wykobi::make_circle<T>(......); wykobi::circle<T> circle1 = wykobi::make_circle<T>(......); std::vector<wykobi::line<T,2>> tangent_lines; wykobi::circle_outer_tangent_lines(circle0,circle1,tangent_lines);

wykobi::triangle<T,2> triangle; wykobi::generate_random_object<T>(0.0, 0.0, width, height, triangle); wykobi::circle<T> circumcircle = wykobi::circumcircle(triangle); wykobi::circle<T> inscribed_circle = wykobi::inscribed_circle(triangle);

wykobi::triangle<T,2> triangle; wykobi::generate_random_object<T>(0.0, 0.0, width, height, triangle); wykobi::circle<T> excircle1 = wykobi::excircle(triangle,0); wykobi::circle<T> excircle2 = wykobi::excircle(triangle,1); wykobi::circle<T> excircle3 = wykobi::excircle(triangle,2); wykobi::triangle<T,2> excentral_triangle = wykobi::create_excentral_triangle(triangle);

The Torricelli point, also known as the fermat point, is the point within the triangle constructed from 3 unique points that minimizes the total distance from each of the 3 points to itself.

wykobi::triangle<T,2> triangle; wykobi::generate_random_object<T>(0.0, 0.0, width, height, triangle); wykobi::point2d<T> torricelli_point = wykobi::torricelli_point(triangle);

const std::size_t max_points = 50; wykobi::segment<T,2> segment; wykobi::generate_random_object<T> (0.0, 0.0, width - 1.0, height - 1.0, segment); std::vector<wykobi::point2d<T>> point_list; point_list.reserve(max_points); wykobi::generate_random_points ( 0.0, 0.0, width - 5.0, height - 5.0, max_points, std::back_inserter(point_list) ); graphic.draw(segment); for (std::size_t i = 0; i < point_list.size(); ++i) { wykobi::point2d<T> closest_point = wykobi::closest_point_on_segment_from_point(segment, point_list[i]); if (wykobi::distance(closest_point,point_list[i]) > T(1.0)) { graphic.draw(wykobi::make_segment(closest_point,point_list[i])); } graphic.draw_circle(point_list[i],3); graphic.draw_circle(closest_point,2); }

const std::size_t max_points = 50; wykobi::triangle<T,2> triangle; wykobi::generate_random_object<T> (0.0, 0.0, width - 1.0, height - 1.0, triangle); std::vector<wykobi::point2d<T>> point_list; point_list.reserve(max_points); wykobi::generate_random_points ( 0.0, 0.0, width - 5.0, height - 5.0, max_points, std::back_inserter(point_list) ); graphic.draw(triangle); for (std::size_t i = 0; i < point_list.size(); ++i) { wykobi::point2d<T> closest_point = wykobi::closest_point_on_triangle_from_point(triangle, point_list[i]); if (wykobi::distance(closest_point,point_list[i]) > T(1.0)) { graphic.draw(wykobi::make_segment(closest_point,point_list[i])); } graphic.draw_circle(point_list[i],3); graphic.draw_circle(closest_point,2); }

const std::size_t max_points = 50; wykobi::rectangle<T> rectangle; wykobi::generate_random_object<T> (0.0, 0.0, width - 1.0, height - 1.0, rectangle); std::vector<wykobi::point2d<T>> point_list; point_list.reserve(max_points); wykobi::generate_random_points ( 5.0, 5.0, width - 1.0, height - 1.0, max_points, std::back_inserter(point_list) ); graphic.draw(rectangle); for (std::size_t i = 0; i < point_list.size(); ++i) { wykobi::point2d<T> closest_point = wykobi::closest_point_on_rectangle_from_point(rectangle, point_list[i]); if (wykobi::distance(closest_point,point_list[i]) > T(1.0)) { graphic.draw(wykobi::make_segment(closest_point,point_list[i])); } graphic.draw_circle(point_list[i],3); graphic.draw_circle(closest_point,2); }

const std::size_t max_points = 50; wykobi::quadix<T,2> quadix; wykobi::generate_random_object<T> (0.0, 0.0, width - 1.0, height - 1.0, quadix); std::vector<wykobi::point2d<T>> point_list; point_list.reserve(max_points); wykobi::generate_random_points ( 5.0, 5.0, width - 5.0, height - 5.0, max_points, std::back_inserter(point_list) ); graphic.draw(quadix); for (std::size_t i = 0; i < point_list.size(); ++i) { wykobi::point2d<T> closest_point = wykobi::closest_point_on_quadix_from_point(quadix, point_list[i]); if (wykobi::distance(closest_point,point_list[i]) > T(1.0)) { graphic.draw(wykobi::make_segment(closest_point, point_list[i])); } graphic.draw_circle(point_list[i],3); graphic.draw_circle(closest_point,2); }

const std::size_t max_points = 50; wykobi::polygon<T,2> polygon; generate_polygon_type2<T>(graphic.width(), graphic.height(), polygon); std::vector<wykobi::point2d<T>> point_list; point_list.reserve(max_points); wykobi::generate_random_points ( 0.0, 0.0, width - 5.0, height - 5.0, max_points, std::back_inserter(point_list) ); graphic.draw(segment); for (std::size_t i = 0; i < point_list.size(); ++i) { wykobi::point2d<T> closest_point = wykobi::closest_point_on_polygon_from_point(polygon, point_list[i]); if (wykobi::distance(closest_point,point_list[i]) > T(1.0)) { graphic.draw(wykobi::make_segment(closest_point, point_list[i])); } graphic.draw_circle(point_list[i],3); graphic.draw_circle(closest_point,2); }

const std::size_t max_points = 50; wykobi::circle<T> circle; wykobi::generate_random_object<T> (0.0, 0.0, width - 1.0, height - 1.0, circle); std::vector<wykobi::point2d<T>> point_list; point_list.reserve(max_points); wykobi::generate_random_points ( 5.0, 5.0, width - 5.0, height - 5.0, max_points, std::back_inserter(point_list) ); graphic.draw(circle); for (std::size_t i = 0; i < point_list.size(); ++i) { wykobi::point2d<T> closest_point = wykobi::closest_point_on_circle_from_point(circle, point_list[i]); if (wykobi::distance(closest_point,point_list[i]) > T(1.0)) { graphic.draw(wykobi::make_segment(closest_point, point_list[i])); } graphic.draw_circle(point_list[i],3); graphic.draw_circle(closest_point,2); }

const std::size_t max_points = 50; wykobi::quadratic_bezier<T,2> bezier; bezier[0] = wykobi::generate_random_point<T>(width - 1.0, height - 1.0); bezier[1] = wykobi::generate_random_point<T>(width - 1.0, height - 1.0); bezier[2] = wykobi::generate_random_point<T>(width - 1.0, height - 1.0); std::vector<wykobi::point2d<T>> point_list; point_list.reserve(max_points); wykobi::generate_random_points ( 0.0, 0.0, width - 1.0, height - 1.0, max_points, std::back_inserter(point_list) ); graphic.draw(bezier,100); for (std::size_t i = 0; i < point_list.size(); ++i) { wykobi::point2d<T> closest_point = wykobi::closest_point_on_bezier_from_point(bezier, point_list[i]); if (wykobi::distance(closest_point,point_list[i]) > T(1.0)) { graphic.draw(wykobi::make_segment(closest_point, point_list[i])); } graphic.draw_circle(point_list[i],3); graphic.draw_circle(closest_point,2); }

const std::size_t max_points = 50; wykobi::cubic_bezier<T,2> bezier; bezier[0] = wykobi::generate_random_point<T>(width - 1.0, height - 1.0); bezier[1] = wykobi::generate_random_point<T>(width - 1.0, height - 1.0); bezier[2] = wykobi::generate_random_point<T>(width - 1.0, height - 1.0); bezier[3] = wykobi::generate_random_point<T>(width - 1.0, height - 1.0); std::vector<wykobi::point2d<T>> point_list; point_list.reserve(max_points); wykobi::generate_random_points ( 0.0, 0.0, width - 1.0, height - 1.0, max_points, std::back_inserter(point_list) ); graphic.draw(bezier,100); for (std::size_t i = 0; i < point_list.size(); ++i) { wykobi::point2d<T> closest_point = wykobi::closest_point_on_bezier_from_point(bezier, point_list[i]); if (wykobi::distance(closest_point,point_list[i]) > T(1.0)) { graphic.draw(wykobi::make_segment(closest_point,point_list[i])); } graphic.draw_circle(point_list[i],3); graphic.draw_circle(closest_point,2); }

const unsigned int max_segments = 10; std::vector<wykobi::segment<T,2>> segment_list; generate_random_segments ( 0.0, 0.0, width - 10.0, height - 10.0, max_segments, std::back_inserter(segment_list) ); wykobi::circle<T> circle; wykobi::generate_random_object<T>(0.0,0.0,width - 1.0,height - 1.0,circle); graphic.draw(circle); for (std::size_t i = 0; i < segment_list.size(); ++i) { if (!wykobi::intersect(segment_list[i],circle)) { wykobi::point2d<T> closest_point_on_segment = wykobi::closest_point_on_segment_from_point (segment_list[i], wykobi::make_point(circle)); wykobi::point2d<T> closest_point_on_circle = wykobi::closest_point_on_circle_from_point (circle, closest_point_on_segment); graphic.draw(segment_list[i]); graphic.draw( wykobi::make_segment (closest_point_on_segment, closest_point_on_circle)); graphic.draw_circle(closest_point_on_segment,4); graphic.draw_circle(closest_point_on_circle ,4); } else { graphic.draw(segment_list[i]); graphic.draw_circle(segment_list[i][0],3); graphic.draw_circle(segment_list[i][1],3); } }

wykobi::circle<T> circle_a = wykobi::make_circle<T>(..., ..., ...); wykobi::circle<T> circle_b = wykobi::make_circle<T>(..., ..., ...); wykobi::point2d<T> closest_point_on_circle_a = wykobi::closest_point_on_circle_from_circle(circle_a, circle_b); wykobi::point2d<T> closest_point_on_circle_b = wykobi::closest_point_on_circle_from_circle(circle_b, circle_a);

wykobi::circle<T> circle_a = wykobi::make_circle<T>(0.0, 0.0, 120.0); wykobi::circle<T> circle_b = wykobi::make_circle<T>(180.0, 140.0, 60.0); wykobi::circle<T> circle_b_inverted = wykobi::invert_circle_across_circle(circle_b, circle_a);

wykobi::line<T,2> mirror_axis = wykobi::make_line<T>(...); wykobi::triangle<T,2> triangle = wykobi::make_triangle<T>(...); wykobi::circle<T> circle = wykobi::make_circle<T>(...); wykobi::triangle<T,2> mirrored_triangle = wykobi::mirror(triangle,mirror_axis); wykobi::circle<T> mirrored_circle = wykobi::mirror(circle,mirror_axis);

wykobi::create_morley_triangle

wykobi::create_cevian_triangle

wykobi::create_anticevian_triangle

wykobi::create_anticomplementary_triangle

wykobi::create_inner_napoleon_triangle

wykobi::create_outer_napoleon_triangle

wykobi::create_inner_vecten_triangle

wykobi::create_outer_vecten_triangle

wykobi::create_medial_triangle

wykobi::create_contact_triangle

wykobi::create_symmedial_triangle

wykobi::create_orthic_triangle

wykobi::create_pedal_triangle

wykobi::create_antipedal_triangle

wykobi::create_excentral_triangle

wykobi::create_incentral_triangle

wykobi::create_extouch_triangle

wykobi::create_feuerbach_triangle

wykobi::inscribed_circle

wykobi::circumecircle

wykobi::nine_point_circle

In the world of computational geometry there are several predicates that form the basis of some of the most complex calculations known. One of the most important predicates is called the orientation predicate.



The question this predicate tries to answer is: Given two points P 0 and P 1 that compose a directed line formed by the vector V (P 1 - P 0 ), and a third point P 2 , on what side of the directed line does the point P 2 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 P 0 and looking towards P 1 , the question as defined before is on what side does P2 reside?



The diagram below depicts the point P 2 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

wykobi::point2d<T> point_0 = wykobi::make_point<T>(...); wykobi::point2d<T> point_1 = wykobi::make_point<T>(...); wykobi::point2d<T> point_2 = wykobi::make_point<T>(...); switch(wykobi::orientation(point_0, point_1, point_2)) { case wykobi::RightHandSide : ... break; case wykobi::LeftHandSide : ... break; case wykobi::CollinearOrientation : ... break; }

Note: The real underlying computation that occurs within the orientation predicate is called the signed area of a triangle or the 3-point determinant. Specifically this is twice the area of the triangle defined by the points P 0 , P 1 , P 2 , or the area of the parallelogram created by mirroring the point P 2 about the axis defined by the vector (P 1 - P 0 )

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.

The orientation predicate can be used to solve some very interesting problems. Below we have two line-segments, S 0 and S 1 , comprised of point pairs (P 0 ,P 1 ) and (Q 0 ,Q 1 ) respectively. The question we would like to answer is whether or not the two line-segments intersect with each other.

By applying the orientation rule, with the directed-line constructed from S 1 as D 0 = Q 0 - Q 1 , and using the end points P 0 and P 1 of the line-segment S 0 , we discover a very interesting property, that is if S 0 intersects D 0 then the orientations for the end points of S 0 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 D 0 and D 1 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 S 0 against the directed-line D 0 constructed from S 1 and an orientation test for the end-points of S 1 against the directed-line D 1 constructed from S 0 .

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

if (orientation(Q 0 ,Q 1 ,P 0 ) to orientation(Q 0 ,Q 1 ,P 1 ) is either not equal or collinear) and (orientation(P 0 ,P 1 ,Q 0 ) to orientation(P 0 ,P 1 ,Q 1 ) 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.

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 |P 1 ,P 0 |, |P 0 ,P 4 |, |P 4 ,P 3 | and |P 2 ,P 1 | all have the same orientation but the edges |P 4 ,P 3 | and |P 3 ,P 2 | have differing orientations to the point Q 0 . From this it is determined that the point Q 0 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.

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

wykobi::segment<T,2> reflection_object = wykobi::make_line<T>(...); wykobi::point2d<T> blue_point = wykobi::make_point<T>(...); wykobi::point2d<T> green_point = wykobi::make_point<T>(...); wykobi::point2d<T> reflection_point; if (point_of_reflection(reflection_object, blue_point, green_point, reflection_point)) { ... }

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?

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 inwards of the angle |ABC| with lines |AB| and |CB| and the circle, N 0 and N 1 respectively. Furthermore we draw two circles with centers at N 0 and N 1 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 N 0 and N 1 intersect each other at two locations I 0 and I 1 . From this we determine the line that passes through B, I 0 and I 1 is the line which bisects the angle |ABC|.

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

wykobi::point2d<T> A = wykobi::make_point<T>(...,...); wykobi::point2d<T> B = wykobi::make_point<T>(...,...); wykobi::point2d<T> C = wykobi::make_point<T>(...,...); wykobi::segment<T,2> AB = wykobi::make_segment<T>(A,B); wykobi::segment<T,2> CB = wykobi::make_segment<T>(C,B); T B_radius = std::min(wykobi::distance(AB), wykobi::distance(CB)) / T(2.0); std::vector<wykobi::point2d<T>> n; wykobi::intersection_point(AB,wykobi::make_circle(B,B_radius),n); wykobi::intersection_point(CB,wykobi::make_circle(B,B_radius),n); if (n.size() != 2) { return; //error - not the right number of intersections } wykobi::point2d<T> i0; wykobi::point2d<T> i1; T N_radius = wykobi::distance(n[0], n[1]) * T(3.0/4.0); wykobi::intersection_point ( wykobi::make_circle(n[0],N_radius), wykobi::make_circle(n[1],N_radius), i0,i1 ); wykobi::line<T,2> bisector1 = wykobi::make_line(B,i0); wykobi::line<T,2> bisector2 = wykobi::make_line(B,i1);

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 (depicted in yellow).

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:

template <typename T> class smooth_corners { public: smooth_corners(const T& expansion_length, const wykobi::polygon<T,2> input_convex_polygon, wykobi::polygon<T,2>& output_convex_polygon) { wykobi::segment<T,2> edge = wykobi::edge(input_polygon,0); T inverter = T(std::abs(expansion_length)); wykobi::point2d<T> test_point = wykobi::segment_mid_point(edge) + wykobi::normalize(perpendicular(edge[1] - edge[0])); if (point_in_polygon(test_point, input_convex_polygon)) { inverter = T(-expansion_length); } for (std::size_t i = 0; i < input_convex_polygon.size(); ++i) { wykobi::segment<T,2> edge = wykobi::edge(input_convex_polygon,i); wykobi::vector2d<T> v = wykobi::normalize(wykobi::perpendicular(edge[1] - edge[0])) * inverter; output_convex_polygon.push_back(edge[0] + v); output_convex_polygon.push_back(edge[1] + v); } } };

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...

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.

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?

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.

A simple .NET WinForms based application is available from the Wykobi downloads page that demonstrates with graphical visualizations all of the above mentioned computational geometry routines and processes.