Fast and improved 2D Convex Hull algorithm and its implementation in O(n log h)

Source code is available at: GitHub and binaries (executable) here: Download OuelletOnlineConvexHull.zip

Note: I got better printing results selecting landscape instead of portraits.

Introduction

This article is the 2^nd of a series of 3

Little request

After reading this article, if you think this algorithm is good enough to be in Wikipedia – Convex hull algorithms, I would be grateful to add a link to Liu and Chen article (or any of the 2 articles I wrote, this one and/or A Convex Hull Algorithm and its implementation in O(n log h)). But please be sure to read this section first: Appendix B – My Wikipedia experience.

Reason for writing back a new article about almost the same algorithm

This article is the continuation of the first one: A Convex Hull Algorithm and its implementation in O(n log h). New article motivations are:

1. The initial reason was to show a newly developed implementation of Ouellet algorithm done in C++. According to provided performance tests, the new implementation is around 4 times faster than Chan implementation in C (also provided).
2. To add and show another implementation of Ouellet algorithm using an AVL tree as a container for Convex Hull points instead of an array based container. There is 2 implementations based on AVL tree: a simple one and a more efficient one with many added optimization.
3. To explain what optimizations had been added to the AVL version to get better performance (Ouellet AVL v2 vs Ouellet AVL).
4. To share my workbench tool to easily compare 2D algorithms using points like the Convex Hull algorithm or any other algorithm that take points as input and produce points as output. For example, the workbench can also be used to test the smallest enclosing circle algorithms (code provided).
5. To show and document few improvements over my first iteration of my algorithm which makes my C# Convex Hull implementation now faster than Pat Morin implementation of Chan algorithm written in C (at least with all "normal" random point generators that I used for testing).
6. To show many real life convex hull algorithm implementation and compare their performance.
7. To show the fastest implementation of a generator of all permutations possibilities of a set of points (code provided).
8. To share other thoughts/discoveries (Multithreading, Wikipedia, Smallest Enclosing Circle, …).

What is a Convex Hull

Definition

According to Wikipedia: "In mathematics, the Convex Hull or convex envelope of a set X of points in the Euclidean plane or Euclidean space is the smallest convex set that contains X. For instance, when X is a bounded subset of the plane, the Convex Hull may be visualized as the shape enclosed by a rubber band stretched around X."

Diagram 1: Example, in red of a convex null (from the provided Workbench)

Applications of Convex Hull in 2D and in 3D

Convex Hull algorithm is a fundamental algorithm in computation geometry, on which are many algorithms in computation geometry based. Also there are a lot of applications that use Convex Hull algorithm.

The Convex Hull in used in many areas where the path surrounding the space taken by all points become a valuable information. There is some example:

• To determine the impedance zone of electrical public utility simulations of their network (IEEE). This is the reason why I met the convex hull few years ago.
• It could serve for input to other algorithms, to improve their performance. For example, a common Smallest Enclosing Circle algorithm can have its performance greatly improved by pre-processing points by a fast Convex Hull algorithm implementation. Difference in performance can be shown in provided benchmark.
• Convex hulls can be used in image processing.
• Used in GIS system to create a feature class containing polygons which represent a specified minimum bounding geometry enclosing each input feature or each group of input features.
• And many others ... 11 applications on Quora

The finding of the Convex Hull is also largely used and studied in universities.

Performance and complexity of Convex Hull algorithms

There are many Convex Hull algorithms with different complexities. Wikipedia - Convex hull algorithms list many of them. Ouellet Convex Hull (or Liu and Chen) algorithm is not in Wikipedia. For more information about why, read Appendix B – Wikipedia at the end of this paper. A good indicator of algorithm performance is usually related to its complexity. It is recommended to choose an algorithm that has the fastest complexity. In the case of convex hull, Results for Speed Test confirms that rule of thumb. Currently, the best known complexity order for a Convex Hull algorithm is O(n log h) where:

• “n” is the number of all source points
• “h” is the set of result points, among all source points, that delimits the path of the Convex Hull.

According to tests and personal experience, in general use cases, “h” is a lot smaller than “n”.  For example, in my tests for a random set of 20 000 000 points in a circle, the Convex Hull is usually made of 200 to 600 points for regular random generators (circle or throw away). Considering the fact that it exists algorithm where the complexity is either: O(n²), O(n log n) and O(n log h). It is quite easy to realize that having a complexity based on “log h” instead of “log n” would improve the performance drastically. But it is also possible to have 2 algorithms in the same complexity where one is really faster than the other. That is the main reason of this article. But as you will see and could test yourself with provided code, results from general cases suggest that Ouellet algorithm is the fastest one and by a pretty good margin. It also has other advantages depending of your needs. When comparing, you should also consider the language the algorithm is implemented. In actual test results, you have to take into account that C# has the reputation to be slower than C. My personal test seems to be in accordance to the fact that C# is quite slower than C++, at least for the implementation of this algorithm.

Ouellet Convex Hull Algorithm

How it works

The best detailed description of the inner working of Ouellet algorithm could be found at

• My first code project article about it: A Convex Hull Algorithm and its implementation in O(n log h)
• Liu and Chen original Scientific Paper: A new algorithm for computing the convex hull of a planar point set.

Additional information can be found later in section: Differences of Ouellet VS. Liu and Chen Convex Hull Algorithm

Algorithm main idea

This is a quick Summary. More information can be found at A convex Hull Algorithm and its implementation in O(n log h)

The algorithm has 3 main steps that occur in sequence:

• Step 1 - Identify virtual quadrant boundary
  • We found the boundary (left, top, right, bottom) of all input points (first pass on all input points).
  • We define 4 virtual quadrant based on boundary points (Q1 is based on top and most right points, its root is based on x of the top point and y of most right point). Same things apply for each quadrant.
  • Each Hull Points will be kept per quadrant (virtual quadrant) and will always be ordered. Note: in all supplied implementations all points are stored counterclockwise for each quadrant.
• Step 2 - Find Hull Points
  • For each input points, we find if it is part of a quadrant and insert as Convex Hull point if applied.
    • Using dichotomy, we search where an input point should be inserted if apply.
    • We insert the input point if it is part of the Convex Hull. If it is the case, we also remove any neighbors for which the new Convex Hull point invalidated them.
• Step 3 – Merging Quadrant results
  • Merge the 4 quadrant Hull Points results, removing any duplicate at quadrant boundaries

Differences of Ouellet VS. Liu and Chen Convex Hull Algorithm

Before stepping into Convex Hull optimizations, it should be clear Liu and Chen Convex Hull algorithm and Ouellet algorithm are based on the same principle: virtual quadrant, at least according to what I understood from the article of Liu and Chen: A new algorithm for computing the convex hull of a planar point set. The article of Liu and Chen has no code or implementation detail and/or link to any of those. Because of that, I decided to rename my first implementation code: Liu and Chen. Many optimizations have been added since then which improve performance. Included benchmark results show those differences.

According to the date of publication: 2007-07-01, Liu and Chen are the first to discover the principle of using virtual quadrant. Optimizations is what could be seen as differences between "Liu and Chen" and "Ouellet".

Optimizations

Following optimizations are optional optimizations that could be applied to one or more of the Ouellet convex hull algorithm. Some are only used in the Ouellet CPP version while most of them were applied in the Ouellet AVL v2. None of them are included in Liu and Chen implementation.

Slope VS Cross Product

In the case of the Ouellet Convex Hull, to insert a new point as a potential convex hull point, it requires to know if it is to the left or right of 2 neighbor points where it should be inserted. It could be done in 2 ways:

Method for determining the side of a point from 2 other ones	Slope	Cross Product
Formula
Description	By comparing 2 slopes. For A, B, C points where A and B are points each side of new point to compare C. Compare slope of (AB) with slope of (AC). Note: The slope (AB) can be easily cached in order to save some calculation next time.	By using cross product directly. In previous formula for A, B, C: if (d>0) then C is to the left if (d=0) then C is on the same line if (d<0) then C is to the right.
Required calculation	2 x (2 subtractions and 1 division) = 4 subtractions and 2 divisions **	5 subtractions and 2 multiplications
Final result in the algorithm	Slightly slower	Slightly faster

Results in algorithm

Usage results in algorithm was a little gain in performance for Cross Product. Just enough to beat Chan algorithm in some cases. But it also had the effect of simplifying the selection of potential convex hull candidate. The cross product could be used exactly the same for all quadrants which enables three things:

• Simplification and better code architecture (more generic) by raising point selection on the base class of each quadrant.
• Simplification by removing the cache management of slopes into point structure
• Smaller footprint by removing the cached slope from the point structure. Can use the already defined framework point structure.

Multiplication VS. Division

Benchmark results :

According to an included benchmark test, the multiplication is 1.14 faster than the division. To know more about testing condition, see Hardware used for testing. Difference in speed between the 2 partly explain why Cross Product is a little faster. Many other variables should be considered like CPU pipelining, memory proximity, and more.

Quick Discard (QD)

The quick discard was the first difference that exists between "Liu and Chen" algorithm and Ouellet algorithm. It is a simple check, added to Ouellet algorithm, which is done at each iteration of the dichotomy search, which in most cases, would quickly stop the search because the point will be rejected without going any further into the dichotomy search. The check to know if we can reject the currently evaluated point does not require any calculation, only one simple comparaison is required against the each point of the dicotomy travel. This advantage comes from the fact that we are using quadrants, and also, because quadrant Convex Hull points are always sorted. In general, the time to do this simple test will save time to perform the all the dichotomy steps because a point could be discarded quickly in first steps. Quick discard is made possible because dichotomy is merging to its insertion point. If we take into account that "h" (count of Convex Hull points) is, in general, a lot less than "n" (count of input points), it should save time in most cases. In other words, there is a great amount of point that should be discarded and the chance of discarding one quickly is high.

Tests have been done between Liu and Chen and Ouellet with the 2 most valuable (closer to real life) random generators: throw away and circles. Liu and Chen has no QD while Ouellet has QD. In bad case (throw away point random generator), QD is just a little bit slower (negligible) and pretty stable for any amount of input points. In good case (circle point random generator), QD is quicker by 1.2 times for 100 000 points and performance increase slowly with the number of points. It could be worse without QD, but so slightly that it is a nice addition and should improve performance in general usage. The same test could be done easily with the provided benchmark.

Opposite quadrant check (OQC)

In Ouellet algorithm, due to the way quadrant are defined, 2 adjacent quadrants can’t have any intersection zone. Odd and even quadrant are mutually exclusive while 2 opposite quadrant can share a zone.

Point quadrant belonging is calculated based on quadrant root point, in first versions of the algorithm. Opposite quadrant can have an intersection zone (see example below). Because of this, the algorithm should do opposite quadrant check (OQC) for any point in the intersection. OQC come from the necessity to check a point on its opposite quadrant although it has been identified as being part of a specific quadrant. That situation only occurs in some specific point distribution but should be taken into account.

The following example shows a situation where the distribution force the algorithm to do opposite quadrant check. In that example, P1 is part of 2 quadrants (Q1 and Q3) at the same time, according to each quadrant root point. This situation occurs mainly when the distribution of point is a narrow diagonal:

Here we have Virtual Q1 in pink and Q3 in green, with each quadrant root point as diamond and same color as its quadrant. As you can see, P1, the yellow point, which should be evaluated to see if it is part of the Convex Hull, is included in Q1 (pink zone) and Q3 (green zone). But if we would evaluate P1 for Q1 and discard it without verifying it against Q3, then we would have discarded a good candidate. To be sure to not miss any quadrant check, the lazy way would be to check each point against each quadrant. But it could easily be optimized by verifying only against an opposite quadrant when a point has been identified as being part of a quadrant. In most algorithm implementation, each point is checked against all possible quadrant in order to be sure to not miss any point in any quadrant. But in Ouellet CPP and Ouellet AVL v2, some optimization has been added (see specific information about each implementation of Ouellet Convex Hull). We can easily see that odd and even quadrants are mutually exclusive. In other words, if a point is part of Q1 or Q3, it can’t be part of Q2 or Q4 and vice-versa.

Opposite quadrant check bypass (OQCB)

Another optimization is as soon as a point is estimated to be a Convex Hull point (not just being part of a quadrant) in any dichotomy step, it is not required to check it against any other quadrant, either the opposite one. It can’t be part of convex hull points in 2 quadrants simultaneously. The check for an opposite quadrant is only mandatory for a point that has been identified has being part of a quadrant without being identified as being a potential hull point. Only quadrant limit points can be in 2 quadrants at the same time which would be merged into only one point at the 3^rd steps of the algorithm. That optimization is only done in Ouellet CPP and Ouellet AVL v2.

Disjoint quadrant check (DQC)

PCZ

First and last points are always defined counterclockwise per quadrant. PCZ is the "potential candidate zone" and is defined by the first and last quadrant points. It is the zone to the right of the first and last point of a quadrant. This zone is always unique per quadrant and define the zone of potential candidates to be added. Any other place outside that zone can’t contain any valid hull point for this quadrant.

DQC Using PCZ

There is a nice optimization that can also be done which is only used in Ouellet AVL v2. It is to verify at start if all quadrants are disjoint. To verify that, we only have to check if the root point of a quadrant is outside of its opposite quadrant PCZ. If it is the case for all quadrants, then if an input point lies into one quadrant, we can skip any other quadrant check, either its opposite one, always.

That optimization would require 2 things:

1. A verification to do before the 2^nd step of the algorithm to determine if quadrants are disjoint.
2. Different code for the 2^nd step of the algorithm depending if quadrants are disjoint or not:
   • If disjoint. use DQZ
   • If not disjoint, use original algorithm code

     

     Q3 is part of Q1 PCZ (yellow). Quadrants are not disjoint because at least one root point (diamond) is in its opposite quadrant PCZ

     

     All root points (diamond) are not part of any of their opposite quadrant PCZ (yellow). DQC optimization can be used.

PCZ point belonging (PCZB)

PCZ is explained at the beginning of disjoint quadrant check. Usage of the quadrant PCZ to determine current point quadrant belonging.

To quickly keep/discard potential candidate point, we can start by checking if it is part of the target quadrant PCZ by using Cross Product instead of using the quadrant root point. Doing cross product has a calculation has a cost but using it could prevent us to do many dichotomy iteration. It also has the advantage that being in PCZ of a quadrant, a point can’t be a hull point of any other quadrant. Ouellet AVL v2 uses this optimization.

Previous quadrant (PQ)

Another optimization only done in "Ouellet AVL v2", is to always start quadrant check with the last one where a point had been added, then cycle with remaining quadrant in order. That optimization could be useful in some situations where consecutive input point to evaluate have a high risk to be neighbors together, and hence being part of the same quadrant. Doing so will also enable us to skip unnecessary other quadrant checks and improve performance. Although that optimization does not bring any advantage in random point placement, it would not add any downside either. It’s an optimization that could be nice in some cases and have no impact in all others.

Algorithms

Algorithms Comparison

General description

*In place mean that point are moved directly into its source container affecting original order. If you can’t accept that, a copy of source points should be done prior to call the algorithm.

Complexity

* Given for AVL implementation. Exact: n+(h*2)+(h*2*sizeof(pointer))

Algorithm implementation information

This is a list of implementation details in order to better understand any differences.

* The current implementation has a very little bug which appears rarely and mainly when there is a small number of points. The bug could be shown by the workbench pretty easily but it is random (I don’t know the reason why). Code provided is made by Pat Morin from University of Carleton.
 

Ouellet Convex Hull Algorithm Implementation comparison

The Code

The benchmark is made with Visual Studio 2017. The framework is WPF with some MVVM. Algorithms are all made in C# except few which are clearly stated which are made of C/C++, see Algorithm implementation information for more information. All tests have been made with code optimized, compiled for "Release".

Note: The only missing algorithm that is in O (n log h) is "The ultimate planar convex hull algorithm". Because there is calculation to find the median, there is a very high probability to be slower than any other actual O(n log h) algorithms. If I would have found any implementation, I would have included it here.

Implementation choices

Array VS. List

A "List" class is a C# collection which uses an array as its underlying container. Using a "List" instead of an array should have similar performance. Tests confirm a very small increase in performance for managing an array directly. The difference is so small that it is hard to justify the lost in clarity using an array. Both collections have been used in different implementation and can be compared together.

Array vs Tree

Array based containers are used in all Ouellet (and Liu and Chen) implementations except for the Ouellet AVL version. Ouellet AVL and Ouellet AVL v2 uses an AVL tree to store potential candidate instead of an array based container. Using an array based container implies a manual dichotomy management. While the tree, by its nature, implement a dichotomy internally. Using dichotomy to get appropriate insertion point ensure good performance and is the main key to stay in O(n log h) . According to my tests and one real usage in life, I suspect that "h" should stay under 1000 in most cases which would not make any difference on using an array or a tree. But having a tree would be safer for cases where "h" could be very large (~ more than 500 000 points on my tests). Otherwise, the items shifting that should occur in array based container for insertion/deletion would become too much prominent and affect performance too much.

The advantages of arrays are (compared to a tree structure):

1. Data is contiguous together, achieving greater CPU cache hit.
2. In most cases, only one call to the heap allocator is required to get an array large enough to accommodate all hull points. That is based on many tests using random point generators as sources of input points as described in benchmark here. Actual implementation reserve 1000 Points at start to reduce heap allocator requests. For 1 000 000 points, results almost always stay under 800 Hull Points with all normal random generators (not "Arc" one).
3. Can use indexing, which enable direct access to point. In fact, direct access would enable to achieve O(n) performance in a multithreaded algorithm but some tests showed slower results than actual implementation. Slower performance appears to be caused by multithreaded locking mechanism. See Convex Hull in O(n) for a little bit more information on a quick test that I did which did not work as expected. Also because "h" is very low in most cases, the addition of multithreading mechanism which applies on "n", should be really low in order to not become more significant than the time related to "h". The O(n) test is not included in provided code.

The major drawback of using an array over a tree structure is:

When the amount of "h" become large, bigger than ~ 500 000 (Ouellet C# implementation), the array container solution become expensive for 2 reasons:

1. The main one is because we have to shift points every time we insert or delete a candidate.
2. Also, we have to allocate a new array each time we reach the reserved space boundary and copy all data to the new array.

The tree solution (in code using an AVL tree) is newer. I mean using an AVL tree instead of an array is the latest major implementation change made. AVL tree implementation appears to be a great choice in order to get an algorithm that is directly dependent on its output size, not higher. That would assure a more constant performance in all cases. But it should not be required in general use cases because "h" (hull points) is a lots smaller than "n" (source points) in all tested cases. The reason why the C++ implementation is not using a tree is because I just recently realized the negative impact of array when "h" become too large and I’m missing time to code another implementation in C++ .

There is 2 most popular tree management algorithms for balanced tree: AVL and "Red-Black" tree. I chose AVL over "Red-Black" tree for 2 reasons:

1. Laziness, it appears to me to be easier to implement AVL over Red-Black tree.
2. According to my tests, because "h" stay extremely small in regards to "n", there should be a lot more read than insertion into the tree which should favor AVL tree over Red-Black tree to get best performance.

Convex Hull in O(n)

I attempted to create an algorithm in O(n) instead of O(n log h). I think I it is possible with the help of threads, array as containers and good design. At least 2 threads are necessary to achieve the goal. One thread that filter potential candidate quickly and another one which insert potential candidate in the result Convex Hull container, an array.

• Thread 1 - Filter candidate in O(1) per point. That thread reject bad candidates and put good one into a stack. To get approximate insertion place into the array of already found candidate convex hull points, the algorithm should use linear positioning between first and last points. That should not give exact location of insertion (linearization of an arc) but should be enough to quickly reject most points. This way, there is no more dichotomy, which corresponds to removing the "log h". That thread does its job in O (n).
• Thread 2- Try to insert potential convex hull point filtered by thread one. That is done in O (~h++ log h). Insertion should occur less and less over time, when approaching n points and/or approaching the solution. That should give enough time to empty the stack of point to try to insert.

If thread 1 finishes always before thread 2 or very close to that (constant time), then we could say we have a Convex Hull in O(n).

The preferred container between the 2 threads is a stack because most recently evaluated points have a higher risk to be better convex hull candidate.

In fact it was working but the time to do the job was twice as slow as Ouellet Single threaded version. Please note that in most cases, the thread 1 finished with no point in the stack of potential candidates, which should indicate a success of being in O(n), at least for general cases.

Also, the algorithm alternate between 2 copies of immutable array in order to stay coherent and not slowdown the thread which filter candidate. That was not a real success for few estimated reasons:

• In general case, it is twice as slow as any other of my implementations, probably due to locking mechanism required to access the shared stack of filtered candidate. Although I used a "SpinLock" instead of a standard "Lock", which should be faster, I feel that there were too much blocking situation which added delays. I should find a not blocking way to insert potential candidates.
• It is highly dependent on a good random distribution of points. It could be pretty bad otherwise because filtering could have pretty bad discard hit, at least the way I implemented my algorithm. I would have to reworks everything in order to make it more efficient and less dependent on the data…if ever possible.
• I fixed an amount of thread per Quadrant which would sleep instead of helping other quadrant to complete their task when their initial Quadrant/job were completed. After seeing results, I decided to stop dreaming about O(n) for real world usage. Because of simplicity and speed of each step in the regular algorithm, adding a locking mechanism slow things down too much. Perhaps in many years, when thousands of millions of points would be few points, and the amount of thread would be higher, then perhaps it would be much regarded and useful.
• I think that in short, the slowness problem comes from the small "h" which in general stay under 1000 which is a tree level of 8-10. Doing 10 iterations is so fast that it competes well against any lock mechanism.

As a side note: Due to my lack of knowledge, I’m not sure if I have the right to say that an algorithm is in "O(n)" if it requires many threads to achieve that performance? I have to look for that also. Everything about finding Hull Points in O(n) would have to be proven and would have been really fun to do but time and money rules the world and I’m part of it. Because of bad performance and because I think I should improve the algorithm, I preferred to not publish that part.

C# vs C++

Although I did not implement C++ and C# version the same, it is easy to see that C++ flexibility, or machine language proximity, has an edge in terms of performance over C#. But there are few things to note about C++ implementation:

• The code is awful.
• It has one very big function with repetitive code.
• There is many "goto".

The C++ code is hard to maintain. C# implementations are pretty much easier to read and easy to follow. There is also more generic functions in base class, and less repetitive code. There is still some room for improvement but most peoples should understand more easily C# implementation over C++ one.

Every choice made in C++ were intentional in order to beat "Chan". Base on test results, in most general cases, at least based on my 4 general generators (circles, 5 circles, rectangles, throw away), for one million of points, the Ouellet CPP implementation is around 4 times faster than Chan. Results come from my benchmark and are available at: Results for "Deep performance test".

Multithreading

According to benchmark results, it is easy to see that multithreaded version is more than twice as fast as the mono thread one. Because of the nature of Ouellet Hull algorithm, it does make it a good candidate for an easy implementation of multithreading, at least for the 2 first steps of the 3. The 2 first steps are dependent on "n" while the 3^rd is only dependent on "h". Of the 3 steps, the first could easily be fully multithreaded, the second could be implemented easily on 4 totally independent threads (no synchronous mechanism required) and the 3^rd step would be a little harder to multithreaded but it is not affected by "n" only by "h".

Although multithreading usage is not as important as the complexity (Big O) of an algorithm, it could help get nice boost in performance. In the case of Ouellet algorithm, it is very easy to add multithreaded and kept threads fully independent which it is an additional bonus.

C# - Safe code VS. Unmanaged code

Unsafe code usage instead of safe code would have probably resulted in better performance. But unsafe code requires a lot of artifacts. Decision has been made that potential speed gain of unsafe code would not suffice to counterparts the difficulty it would add to read and maintain the code.

Container usage remark

2 different types of container had been used:

• Array (or list which uses array under its hood)
• AVL tree

Both containers have pretty similar speed in general cases. I started with a "List" because it was simpler to me. But there is one main drawback with array derived container: when "h" (the count of hull point) become too big, its performance becomes quickly accursed. The number of points where array become a real problem of performance is around ~100000 - ~500000. That’s why I decided to create a new version with AVL tree for those pretty rare cases. By its nature, an AVL tree is the perfect container for Ouellet algorithm in single thread. It is not affected by "h" like array is. The main difference is

• Array : Dicotomy is done externally. Insertion and deletion require moving the objects close the index implied (close because of invalidated neighbors).
• AVL Tree: Dicotomy come from the internal way a tree is done. We need to rebalance the tree in order to stay in O(n log h).

When the number of result points (hull points) become too big, it is obvious that AVL tree win easily.

C# Array.Copy vs MSVCRT memmove vs MSVCRT memcpy (for immutable)

If you have to copy array of data quickly and you know that type of inner array object is of simple type (or struct of simple type), then you can use MSVCRT memcpy (not overlapping copy) or memmove (overlapping copy). The function MSVCRT - memcpy if the fastest but, according to MSDN, it can’t be used for overlapping data. Otherwise Array.Copy and memmove are identical with a very consistent but negligible time advantage for memmove. To use memcpy, we require to copy the entire array each time (same as using immutable array).

Please note, in a language with a garbage collector, moving any object (not simple type or struct containing simple type) require to pin that object in memory in order for them to not being moved by the garbage collector when being moved in memory by unmanaged code. An array is an object. An array itself should be pinned when used from non-GC aware code.

C# Buffer.Copy

I can’t use Buffer.Copy to move points. Buffer.Copy does not support struct either if they are made of simple types only. See exception below:

If anybody has a clean explanation of why it is like that, I would really appreciate to know it.

C# sizeof(T)

C# sizeof(T) can’t be used on a generic type because the size of object can’t be known at compile time. I wonder if it would not be possible to know the size of the object at compile time, for example for simple type and struct of simple type? Wouldn’t be possible to have a "where" clause that would apply on T which specifies "SimpleType"? If feasible, I think it would be a nice addition the C# language to have the possibility to use "sizeof" on T.

Compiling/running the code

The code has been done using Visual Studio 2013 and had been ported to Visual Studio 2015 and then to Visual Studio 2017. It currently uses some features of VS. 2017. It is advised to use that version. If you can’t do so, you will probably have some trouble with C++ code due to its specific binding to a version of MSVCRT.dll which come with Visual Studio 2017.

If you have an exception saying BadImageFormat, you should ensure to run in 64 bits. The problem is a wrong image CPU target (64 vs 32 bits). The code in C++ is made in 64 bits only. If you remove the C++ code, then you can build with "AnyCPU" setting.

Algorithm Benchmark

All benchmark results have been made using the code provided. Everything is testable and verifiable.

Random generators

The benchmark has five different type of point generator:

Hardware used for testing

Just as additional information, this is my machine specs:

Results for "Speed Test"

Most results presented are for algorithms in O(n log h) to better differentiate performance of fastest ones. Most tests were done with 2 generators, “Circle” and “Throw away” which should be closer to real life usage. Any other combinations of algorithms, random generators can be easily tested with the code provided. You can also add your own generator or algorithm if you prefer.

Results for implementation made in C/C++ are taken inside native language, not in C# code to prevent to take results conversion time in account. That should give a better idea of real time taken by each algorithm/implementation.

A linear regression directly based on result points is also added in order to see the tendency and if it is linear or not.

Conclusion about "Speed Test" results

Observations

1. Advantages of O(n log h) is obvious.
2. Monotone chain algorithm is the slowest of evaluated algorithm implementation
3. Although MiConvexHull and Heap Algorithms are both in O(n log n), Heap is pretty much slower.
4. Speed depends on source points, their placement and order. Speed results could vary depending on point arrangement (random generator).
5. All convex hull algorithms in O(n log h) are faster but good differences exist between them.
6. Optimizations made in Ouellet, in regards to Liu and Chen, bring better performance.
7. Array based container implementations are dependent on very large "h". The problem appears when "h" reach about ~500 000 points.
8. Chan is more independent to the type of random generator than Ouellet. This difference could be seen when comparing results for different random generators.
9. Ouellet AVL and Ouellet AVL v2 are more affected by "h" than Chan. It can be seen by comparing results of "Arc generator" with any other generator. It should be caused by heap allocator delay and/or rebalancing the tree.
10. Multithreading is a real advantage.
11. C++ is a lot faster than C#. Although not being the exact same algorithm, the difference is lots higher than expected. I think it could be attributed to extensive memory read/write which could have less indirection in C++ and also due to array boundary checks which happen automatically in C# for safety but not in C++.
12. C# Single thread vs Chan is different depending on the type of generator. While for "Circle" point generator, C# is faster than Chan in C. When using the "Throw away" point generator, Chan wins easily.
13. Ouellet CPP is 4x faster than Chan where both are made in the same language (C/CPP). That is for a number of points higher than 100000.
14. Immutable array copy is affected by a lot of points a lot quicker than when array is used with overlapping copy (normal usage, like into a list).
15. Although Ouellet CPP clearly wins in normal usage, it performs badly when "h" become too large due to array based usage for storing hull points.
16. The difference between Ouellet AVL vs Chan is pretty big for the arc generator (where n=h). But this difference is reducing when the amount of point augment (become huge). See Results for "Deep performance test".
17. For Ouellet, an AVL tree seems to be the best container to ensure good performance in all cases. Although a little slower than an array in general cases, it has a more contant behavior when "h" become very large.

Warning about the "Arc" Point Generator

Some tests were done using the "Arc" random generator. Output from that generator is the worst case data generator that could be used to feed the Ouellet algorithm when it is based on an array container (array or list). Points generated by the "Arc" point generator will all be part of the convex hull, no exception. All points form an arc which is the border of an ellipse in the 4^th quadrant. This should not happen in real life but it was fun to test here what could be a worst case for Ouellet algorithm and how it could affect some of its implementations.

Derived Conclusions

Based on results, we can extrapolate that it would be possible to get a Ouellet AVL v2 "CPP" implementation that could be faster than Ouellet CPP. That implementation wouldn’t have the "Array" container effect that affects all Ouellet implementation based on array. But due to the usage of the heap allocator for each insertion of potential hull point candidate, the speed gain will probably not be as good as the one between Ouellet ST and Ouellet CPP.

Results for "Deep performance test"

This test mainly shows a comparison of algorithm speed. Results are divided in two sections:

1. The first table show raw speed results
2. The second table shows the same information as a ratio relative to the fastest algorithm on the same row

How "Deep performance test" are achieving

• For each Point count (column: Points)
  • Do 10 times each test
    • Create a random number of points (ptsCount)
    • Execute each algorithm calculating its time taken
  • For each algorithm, do an average based on the 10 tests (to have a better normalized average of each algorithm)

An advantage of this type of test is that it is done on a specific amount of point but for each quantity, an average is made over 10 different sets of random points for which each algorithm is tested against the same set of points. Calculating an average on 10 tests helps to have more reliable results that better reflect reality. It lower the variations that could occur on performance of algorithm calculation.

Some sample results:

About "Arc" random point generator, please read: Warning about the "Arc" point generator

Please keep in mind that most algorithms are written in C# which is a language that has a reputation to be quite slower than C (Chan and Ouellet CPP are in C).

What’s left

There are many things left:

• Test the Ouellet Avl with a Red-Black Tree instead of an Avl tree and compare performance
• Implement AVL versions in C++ to find the language advantage and see if the difference is the same as for Ouellet ST
• Add some or all Optimizations in all versions
• Do a multithreaded version in C++, like the 2 done in C#, with an AVL Tree
• Improve the multithreaded version in order to use all threads in all cases on each pass.
• Find better implementation choices and make more testing on an algorithm that run in O(n)
• Adapt the algorithm to do 3D and/or any dimensions by making it more generic. It should be feasible. As an idea, we can use bit values as quadrant.
• Combine any of the preceding to get the best algorithm implementation ever

Who has enough courage, time, money and knowledge to do any of them or all?

Appendix A - History of Ouellet Convex Hull

Appendix B – Wikipedia

On Wikipedia history page of the "Convex hull algorithms", you can see a trial, by myself, to add a reference to my first Code Project article and Liu and Chen article. Both attempt has been removed by David Eppstein. Up to now, there is no reference to this algorithm in Wikipedia which makes it harder to discover.

Appendix C - WPF Chart Control

In the included code, the graphic control used is OxyPlot. This control works fine and it is simple to use. But its performance is on the slow side. In fact, I never saw any WPF only graph control that could handle many thousands of points smoothly. There is one control that I use at work that deliver proper performance and it has a nice API too: "LightningChart" from Arction. I also think that SciChart make a really good control too. LightningChart uses DirectX under the hood and I think SciChart is also using DirectX to ensure acceptable performance. Both have C# control and are very easy to use.

Arction offer exceptional support and SciChart seems to offer something similar.

My experience of C# Chart control:

Appendix D - Excel

Just as a side note, I used ClosedXML as the library to create my excel report. I really liked it. It was easy to learn and use and it’s free.

In our company, we are using SpreadSheetGear which I prefer but it cost money.

Appendix E - Smallest Enclosing Circle

Also called the minimum enclosing circle, the smallest enclosing circle is, as its name implies, the smallest circle of a set "S" of points that include all "S" points in it. A good definition could be found in Wikipedia: Smallest-circle problem. I added the code from Rod Stephens web page: Find a minimal bounding circle of a set of points in C# as an available algorithm. The smallest enclosing circle takes a set of points as its source and produce a result compose of a point, the center, and a length, the ray of the circle. By adding a small, and easy to do algorithm, we can transform the center point and the ray into a circle. In fact the circle is represented as an array of consecutive points forming that circle. The count of points should be enough to give the illusion of a circle. We then can use that algorithm with the workbench exactly the same as we use for the convex hull. The actual implementation code from Rod Stephens works well, but it is quite slow. I also added another available algorithm which is based on Rod algorithm. The difference is that "S" points are pre-filtered with Ouellet convex hull algorithm. The result of the convex hull is then used as a source for Rod algorithm. The algorithm become a lot faster because it drastically reduce the number of source points. This is another advantage of having an efficient convex hull algorithm. Please note that Rod algorithm seems to be in O(n²). You can easily test performance with the provided test workbench where I included both implementation. I recommend to stay under ~1000 points because of the complexity of the original Rod algorithm. See a quick sample of a benchmark test results between the 2 algorithms:

Appendix F - Generic implementation of the fastest algorithm generating all possible permutations of a set of items.

As a side note. I think I have made the quickest C# implementation for generating all possible permutations of a set of items. You can take a look at my answer at StackOverflow question generating permutations of a set (most efficiently) to have an idea of it. I did it in order to make some test on the Convex Hull. All the code is included in the provided source code of this article.

Main advantages over few algorithms found are:

• Heap's algorithm (Single swap per permutation)
• No multiplication (like some implementations seen on the web)
• Inline swap
• Generic
• No unsafe code
• In place (very low memory usage)
• No modulo (only first bit compare)

Appendix G – The benchmark

The benchmark code has been made with Visual Studio 2017 in C#. It uses ".Net framework" 4.5.2 and newer version. All the source code for the benchmark and each algorithm implementation is provided. The binary could also be downloaded. The language is C# for all but 2 Convex Hull implementation, Chan done by Pat Morin and Ouellet CPP done by me.

It is not fully MVVM but uses some features of it. I had to decide to stop somewhere and it is like it is now. I feel it is working enough to be usable and easy to add your own algorithm.

How to Add Your Own Algorithm

Please look into the project "ConvexHullWorkbench" for the class "AlgorithmManager". Just add a line: "_algorithms.Add(..." at the end of the constructor.

Screenshots

Acknowledgment

• Thanks to Omar Saad, my project manager, for explaining to me what a convex Hull is and left me the time to test my idea, write articles and code.
• Thanks to Paul Labbé, a co-worker, for showing me Cross Product. It was one of the best improvements on my algorithm.
• Thanks to Jacques Bherer and Chuma Francis Mugombozi, co-workers, to revise my article. My article would have been a killer one if I would have followed each and every of their recommendations! But I'm too lazy...

References

• Wikipedia, Convex hull algorithms
• Liu, G. & Chen, C. J. Zhejiang Univ. - Sci. A (2007) 8: 1210. A new algorithm for computing the convex hull of a planar point set (https://rd.springer.com/article/10.1631%2Fjzus.2007.A1210)
• Eric Ouellet, IREQ, A Convex Hull Algorithm and its implementation in O(n log h)
• Wikipedia, Smallest-circle problem
• Rod Stephens, C# Helper, 2014-08-13, Find a minimal bounding circle of a set of points in C#
• Wikipedia, Heap's algorithm
• SimpleVar, Stack Overflow question/answer, Generating permutations of a set (most efficiently)
• Qwertie (David Piepgrass), Monotone Chain implementation at GitHub and 2D Convex hull in C#: 40 lines of code at Loyc.net
• Ceometric, .Net library geometry solutions
• futura-sciences.com forums for Tryss answer about Cross product
• bitlush, Efficient AVL Tree in C#

Version

Date of version submited. There could be a delay (day(s) ?) before the article become online. 

1. 2017-10-13, Initial version of this article.
2. 2017-10-18, Corrected some algorithm complexity. Re-arrange convex hull algorithm/complexity tables in order to make them easier to read. Added some information. Minor corrections. 
3. 2017-10-28, Updated binaries executable files to inlcude MSVCRT. This is to ensure 'C' code works fine on any Windows.
4. 2017-11-17, Minor syntax corrections or few hyperlinks added.
5. 2017-11-30, Minor syntax corrections and added some clarifications. Some updated graphics with Delaunay/Voronoi information. Added 1 graphic (conclusion with circle generator).
6. 2017-12-16, Minor corrections in text.
7. 2018-01-30, Added support for 'Online' Convex Hull (dynamic add) in O (log h) per point. Also added support for Peek (IsHullPoint) in order to know if a point would be a hull point before inserting it. All of that only in code for the moment, included in project OuelletConvexHullAvl2Online in GitHub. REmoved some dead code. Slight improvement (standardization). Added a Graham Scan modified version implementation , from Rod Stephens.
8. 2018-03-01, Update executable to reflect recent changes. Add table or related articles.

Thank you for taking time to let me know what you think about this article ...