Click here to Skip to main content
Click here to Skip to main content
Go to top

Computational Geometry, C++ and Wykobi

, 24 Sep 2010
Rate this:
Please Sign up or sign in to vote.
A brief introduction in computational geometry algorithms using Wykobi and C++

Wykobi Polygon Clipping - Copyright Arash Partow

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.

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 {};

The Line Type

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

Wykobi Line - Copyright Arash Partow

template <typename T, unsigned int Dimension>
class line : public geometric_entity {};

The Segment (Line-Segment) Type

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

Wykobi Segment - Copyright Arash Partow

template <typename T, unsigned int Dimension>
class segment : public geometric_entity {};

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.

Wykobi Ray - Copyright Arash Partow

template <typename T, unsigned int Dimension>
class ray : public geometric_entity {};

The Triangle Type

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

Wykobi Triangle - Copyright Arash Partow

template <typename T, unsigned int Dimension>
class triangle : public geometric_entity {};

The Rectangle Type

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.

Wykobi Rectangle - Copyright Arash Partow

template<typename T>

class rectangle : public geometric_entity {};

The Quadix (Quadrilateral) Type

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.

Wykobi Quadix - Copyright Arash Partow

template <typename T, unsigned int Dimension>

class quadix : public geometric_entity {};

The Polygon Type

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

Wykobi Polygon - Copyright Arash Partow

template<typename T, unsigned int Dimension>

class polygon : public geometric_entity {};

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.

Graham Scan (Complexity O(nlogn))

Wykobi Graham Scan Convex Hull - Copyright Arash Partow

const std::size_t max_points = 100000;
std::vector<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));

Jarvis March (Complexity O(nh))

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

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_jarvis_march<wykobi::point2d<T>>(point_list.begin(),
                                                                point_list.end(),
                                                                std::back_inserter(convex_hull));

Melkman (Complexity O(n))

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 Melkman Convex Hull - Copyright Arash Partow

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);

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.

Wykobi Minimum Bounding Ball - Copyright Arash Partow

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

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);

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.

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);

Naive Algorithm (O(N4))

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

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);

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.

Concave Polygon Clipped Against A Rectangle

Wykobi Sutherland Hodgman Rectangle Polygon Clip - Copyright Arash Partow

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);

Concave Polygon Clipped Against A 2D Triangle

Wykobi Sutherland Hodgman Triangle Polygon Clip - Copyright Arash Partow

wykobi::triangle<T,2> 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);

Concave Polygon Clipped Against A 2D Quadix

Wykobi Sutherland Hodgman Quadix Polygon Clip - Copyright Arash Partow

wykobi::quadix<T,2> 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);

Concave Polygon Clipped Against A 2D Convex Polygon

Wykobi Sutherland Hodgman Polygon Polygon Clip - Copyright Arash Partow

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::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);

Cohen-Sutherland Line Segment Clipping

Line Segments Clipped Against An 2D Axis Aligned Bounding Box

Wykobi Cohen-Sutherland Line Segment Clipped Agains A Rectangle - Copyright Arash Partow

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);
   }
}

Line Segments Clipped Against A 2D Triangle

Wykobi Cohen-Sutherland Line Segment Clipped Agains A Triangle - Copyright Arash Partow

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);
   }
}

Line Segments Clipped Against A 2D Quadix

Wykobi Cohen-Sutherland Line Segment Clipped Agains A Quadix - Copyright Arash Partow

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,1000.0,1000.0,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);
   }
}

Line Segments Clipped Against A Circle

Wykobi Cohen-Sutherland Line Segment Clipped Agains A Circle - Copyright Arash Partow

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);
   }
}

Rectangles Clipped Against A Rectangle

Wykobi Rectangles Clipped Against a Rectangle - Copyright Arash Partow

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);
   }
}

Group Based Pairwise Intersections

Segment To Segment Intersections

Wykobi Segment to Segment Pairwise Intersections - Copyright Arash Partow

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));

Circle To Circle (Disk) Intersections

Wykobi Circle to Circle Pairwise Intersections - Copyright Arash Partow

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));

Simple Polygon Triangulation (Ear-Clipping Algorithm)

Wykobi Simple Polygon Triangulation - Copyright Arash Partow

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 Simple Convex Polygon Triangulation - Copyright Arash Partow

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 Simple Convex Polygon Triangulation - Copyright Arash Partow

wykobi::polygon<T,2> polygon;
generate_polygon_type1<T>(width,height,polygon);  // generate simple complex polygon

std::vector<wykobi::triangle<T,2>> triangle_list;
wykobi::algorithm::polygon_triangulate< wykobi::point2d<T> >(polygon,std::back_inserter(triangle_list));

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.

Wykobi Axis Projection Descriptor - Copyright Arash Partow

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.


Beziers And Splines

Random 2D Quadratic Beziers

Wykobi Random Quadratic Beziers - Copyright Arash Partow

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);
}

Random 2D Cubic Beziers

Wykobi Random Cubic Beziers - Copyright Arash Partow

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);
}

Pairwise Segment To Quadratic Bezier Intersections

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.

Wykobi Pairwise Segment To Quadratic Bezier Intersections - Copyright Arash Partow

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);
   }
}

Pairwise Segment To Cubic Bezier Intersections

Wykobi Pairwise Segment To Cubic Bezier Intersections - Copyright Arash Partow

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::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[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);
   }
}

Random Points

Random Points In 2D AABB

Wykobi Random Points In 2D AABB - Copyright Arash Partow

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));

Random Points In 2D Triangle

Wykobi Random Points In 2D Triangle - Copyright Arash Partow

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));

Random Points In 2D Quadix

Wykobi Random Points In 2D Quadix - Copyright Arash Partow

const std::size_t max_points = 1000;
wykobi::quadix<T,2> quadix;
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));

Random Points In Circle

Wykobi Random Points In 2D Circle - Copyright Arash Partow

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));

Constructions, Triangles And Points Of Interest

Vertex Bisector

Wykobi Vertex Bisector - Copyright Arash Partow

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);

Circle Tangent Line Segments

Wykobi Point To Circle Tangents - Copyright Arash Partow

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);

Triangle Circumcircle And Inscribed Circle

Wykobi Circumcircle And Inscribed Circle - Copyright Arash Partow

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);

Construction Of Triangle's Excentral Triangle And Excircles

Wykobi Excentral Triangle - Copyright Arash Partow

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);

Calculation Of The Torricelli Point (Fermat Point)

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 Torricelli Point - Copyright Arash Partow

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);

Closest Point On Segment From External Points

Wykobi Closest Point On Segment From External Points - Copyright Arash Partow

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);
}

Closest Point On Triangle From External Points

Wykobi Closest Point On Triangle From External Points - Copyright Arash Partow

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);
}

Closest Point On Rectangle From External Points

Wykobi Closest Point On Rectangle From External Points - Copyright Arash Partow

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);
}

Closest Point On Quadix From External Points

Wykobi Closest Point On Quadix From External Points - Copyright Arash Partow

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);
}

Closest Point On Polygon From External Points

Wykobi Closest Point On Polygon From External Points - Copyright Arash Partow

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);
}

Closest Point On Circle From External Points

Wykobi Closest Point On Circle From External Points - Copyright Arash Partow

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);
}

Closest Point On Quadratic Bezier From External Points

Wykobi Closest Point On Quadratic Bezier From External Points - Copyright Arash Partow

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);
}

Closest Point On Cubic Bezier From External Points

Wykobi Closest Point On Cubic Bezier From External Points - Copyright Arash Partow

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);
}

Closest Point On Circle From External Segments

Wykobi Closest Point On Circle From An External Segment - Copyright Arash Partow

const unsigned int max_segments = 10;
std::vector<wykobi::segment<T,2> > segment_list;
generate_random_segments<double>(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);
   }
}

Closest Point On Circle From Another Circle

Wykobi Closest Point On Circle From Another Circle - Copyright Arash Partow

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);


Invert A Circle Across Another Circle

Wykobi Invert A Circle Across Another Circle - Copyright Arash Partow

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);


Mirroring Of Objects About An Arbitary Axis

Wykobi Mirror Objects About An Arbitary Axis - Copyright Arash Partow Wykobi Mirror Objects About An Arbitary Axis - Copyright Arash Partow

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);


Other Various Triangle Constructions

wykobi::triangle<T,2> result_triangle;
result_triangle = wykobi::create_morley_triangle(triangle);
result_triangle = wykobi::create_cevian_triangle(triangle,point);
result_triangle = wykobi::create_anticevian_triangle(triangle,point);
result_triangle = wykobi::create_anticomplementary_triangle(triangle);
result_triangle = wykobi::create_inner_napoleon_triangle(triangle);
result_triangle = wykobi::create_outer_napoleon_triangle(triangle);
result_triangle = wykobi::create_inner_vecten_triangle(triangle);
result_triangle = wykobi::create_outer_vecten_triangle(triangle);
result_triangle = wykobi::create_medial_triangle(triangle);
result_triangle = wykobi::create_contact_triangle(triangle);
result_triangle = wykobi::create_symmedial_triangle(triangle,point);
result_triangle = wykobi::create_orthic_triangle(triangle);
result_triangle = wykobi::create_pedal_triangle(point, triangle);
result_triangle = wykobi::create_antipedal_triangle(point,triangle);
result_triangle = wykobi::create_excentral_triangle(triangle);
result_triangle = wykobi::create_incentral_triangle(triangle);
result_triangle = wykobi::create_extouch_triangle(triangle);
result_triangle = wykobi::create_feuerbach_triangle(triangle);

wykobi::circle<T> result_circle;
result_circle = wykobi::inscribed_circle(triangle);
result_circle = wykobi::circumecircle(triangle);
result_circle = wykobi::nine_point_circle(triangle);

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?

Wykobi Orientation Predicate 01 - Copyright Arash Partow

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.

Wykobi Orientation Predicate 02 - Copyright Arash Partow

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 Orientation Predicate 03 - Copyright Arash Partow
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. 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.

Wykobi Segment Intersection Via The Orientation Predicate 01 - Copyright Arash Partow

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.

Wykobi Segment Intersection Via The Orientation Predicate 02 - Copyright Arash Partow

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.

Wykobi Segment Intersection Via The Orientation Predicate 03 - Copyright Arash Partow

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.

Wykobi Segment Intersection Via The Orientation Predicate 04 - Copyright Arash Partow

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 Segment Intersection Via The Orientation Predicate 05 - Copyright Arash Partow

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.

Wykobi Point In Convex Polygon Via The Orientation Predicate 01 - Copyright Arash Partow

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.

Wykobi Point In Convex Polygon Via The Orientation Predicate 02 - Copyright Arash Partow

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.

Wykobi Point In Convex Polygon Via The Orientation Predicate 03 - Copyright Arash Partow

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.

Wykobi Point Of Reflection 01 - Copyright Arash Partow

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.

Wykobi Point Of Reflection 02 - Copyright Arash Partow

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 Point Of Reflection 03 - Copyright Arash Partow

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?


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

Wykobi Bisecting Angle 01 - Copyright Arash Partow

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.

Wykobi Bisecting Angle 02 - Copyright Arash Partow

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

Wykobi Bisecting Angle 03 - Copyright Arash Partow

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 (2 != n.size()) { 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);

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.

Wykobi Smoothing Sharp Corner 01 - Copyright Arash Partow

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.

Wykobi Smoothing Sharp Corner 02 - Copyright Arash Partow

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.

Wykobi Smoothing Sharp Corner 03 - Copyright Arash Partow

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.

Wykobi Smoothing Sharp Corner 04 - Copyright Arash Partow

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.

Wykobi Correctness In Rotations - Copyright Arash Partow

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

Wykobi .NET Demo - Copyright Arash Partow

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

License

This article, along with any associated source code and files, is licensed under The GNU General Public License (GPLv3)

Share

About the Author

Arash Partow

Heard Island And Mcdonald Islands Heard Island And Mcdonald Islands
I'm a computer scientist by ambition and an artisan engineer by trade.

You may also be interested in...

Comments and Discussions

 
GeneralVisual Studio 2010 Pinmembermike taj.10-Mar-11 5:34 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

| Advertise | Privacy | Mobile
Web02 | 2.8.140916.1 | Last Updated 25 Sep 2010
Article Copyright 2007 by Arash Partow
Everything else Copyright © CodeProject, 1999-2014
Terms of Service
Layout: fixed | fluid