Ramer-Douglas-Peucker algorithm

Ramer-Douglas-Peucker is a great algorithm for reducing the number of samples in a given trace, and also for keeping its general shape (as much as he can).

I have a trace with 32.000.000 samples, and I want to compute a reduced version in order to be able to plot it into my winform application.

The first time I used this implementation, time execution was about ~3mn on 32M samples, and the reduced version has 450K samples with Tolerance = 12, and the shape was amazingly preserved.

I tried to improve it, and the time execution is now about ~14sec.

• Getting rid of Math.Pow made me win 121sec...

• Using unsafe code + pointers --> -17sec

• Splitting the work into different threads --> -28sec (on a Core 2 Duo)

Someone could maybe tell me if there is a way to go faster?

```
public class RDPWorker
{
public int start;
public double Tolerance;
public volatile float[] samples;
public List samplesToKeep;

public void DoWork()
{
int h = samples.Length / 2;
float[] half = new float[h];
unsafe
{
fixed (float* ptr = samples, ptrH = half)
{
float* _half = ptrH;
for (int i = start; i pointIndexsToKeep = new List();

pointIndexsToKeep.Add(firstPoint);
pointIndexsToKeep.Add(lastPoint);

int h = points.Length / Environment.ProcessorCount;
List workers = new List();
List threads = new List();
int cpu = 0;

while (cpu ();
_w.Tolerance = Tolerance;
_w.start = h * cpu;
_w.samples = points;
workers.Add(_w);
++cpu;
}

Stopwatch sw = new Stopwatch();
sw.Start();
foreach (RDPWorker worker in workers)
{
Thread t = new Thread(worker.DoWork);
t.IsBackground = true;
threads.Add(t);
t.Start();
}

foreach (Thread thread in threads)
thread.Join()

Note that the algorithm only ever compares distances never using their actual value.

This means that you can speed up your algorithm significantly by directly comparing the squares of the distances instead. This works because the sqrt function is monotonic and thus sqrt(a) < sqrt(b) if and only if a < b.

This allows to change your DouglasPeuckerReduction() function to this:

internal static void DouglasPeuckerReduction(float[] points, int firstPoint, int lastPoint, Double toleranceSquared, ref List pointIndexsToKeep)
{
    float maxDistanceSquared = 0, tmp = 0, areaSquared = 0, X = 0, Y = 0, bottomSquared = 0, distanceSquared = 0;
    int indexFarthest = 0;

    unsafe
    {
        fixed (float* samples = points)
        {
            for (int i = firstPoint; i  maxDistanceSquared)
                {
                    maxDistanceSquared = distanceSquared;
                    indexFarthest = i;
                }
            }
        }
    }

    if (maxDistanceSquared > toleranceSquared && indexFarthest != 0)
    {
        //Add the largest point that exceeds the tolerance
        DouglasPeuckerReduction(points, firstPoint, indexFarthest, toleranceSquared, ref pointIndexsToKeep);
        pointIndexsToKeep.Add(indexFarthest);
        DouglasPeuckerReduction(points, indexFarthest, lastPoint, toleranceSquared, ref pointIndexsToKeep);
    }
}

and call it with the tolerance squared and you'll get the same result in almost half the time. Note that your special Sqrt function is not used anymore and can now be removed.

Also note that I changed the order in which the pointsToKeep list is populated in the last three lines. This adds the indices in order and allows you to remove the sorting step in line 83 which should also save you some time.

StackExchange Code Review Q#29002, answer score: 6