Topic: {Challenge} connect each two points whose distance is shorter than given value

[pkohut]

Hey Highflybird,

I've taken the liberty to use your QuickSortY and GetClosestPairs (renamed to DandC, for Divide and Conquer) routines to put another test together.  They've been reworked to get max speed, at least as much as I could squeeze out of 'em.  More speed could probably be had if the routines used only 1 or 2 parameters, that way the parameters would fit in the CPU registers.  That would require setting up a structure to hold all the needed data for the routines and pass it by reference.

Other notable changes - arrays are replaced with std::vector, the quick sort routine replaced with a faster vertex sort, and the collection of points is done in a manner similar to what Dan and I did previously.

These timings are the total run times after entering the input distance and include the time spent collecting and sorting 1,000,000 points using a distance of 70.
AcDbLine creation code disabled
1.464 seconds

AcDbLine creation code enabled
11.147 seconds

VS project file attached.

Code: [Select]

#include "StdAfx.h"
#include <vector>
#include <algorithm>

struct IndexAndSide
{
    int N; // index of point
    bool R; // is right side
};

/////////////////////////////////////////////////////////////////////////////
//* Neat little sort routine, shuffling an array of IndexAndSide objects. *//
/////////////////////////////////////////////////////////////////////////////
void QuickSortY(const AcGePoint3d * pts, IndexAndSide * data, int low, int high)
{

    if(low < high) {
        IndexAndSide tmp = data[low];
        int i = low;
        int j = high;
        while(i < j) {
            while(i < j && pts[data[j].N].y > pts[tmp.N].y)
                --j;
            if(i < j)
                data[i++] = data[j];
            while(i < j && pts[data[i].N].y < pts[tmp.N].y)
                ++i;
            if(i < j)
                data[j--] = data[i];
        }
        data[i] = tmp;

        QuickSortY(pts, data, low, i - 1);
        QuickSortY(pts, data, i + 1, high);
    }
}

///////////////////////////////////////////////////////////////////////
//* Hmm, append a line (from pt1 to pt2) to the block table record. *//
///////////////////////////////////////////////////////////////////////
void AddLine(AcDbBlockTableRecord * pRec, const AcGePoint3d & pt1, const AcGePoint3d & pt2)
{
    AcDbLine *pLine = new AcDbLine(pt1, pt2);
    if(pRec->appendAcDbEntity(pLine) != Acad::eOk)
        delete pLine;
    pLine->close();
}

/////////////////////////////////////////////////////////////////////////////////
//* DandC (Divide and Conquer), this is where the actual work-work begins and *//
//* its smokin'                                                               *//
/////////////////////////////////////////////////////////////////////////////////
void DandC( AcDbBlockTableRecord * pRec, const AcGePoint3d * pts, int low, int high, const double & distance )
{
    int index = high - low + 1;

    if(index < 3) {
        if(index == 2) {
            // only 2 points in the array, are they less than distance appart?
            if( fabs(pts[high].x - pts[low].x) <= distance &&
                fabs(pts[high].y - pts[low].y) <= distance &&
                pts[low].distanceTo(pts[high]) <= distance) {
                    // yep, build a line
                    AddLine(pRec, pts[low], pts[high]);
            }
        }
        return;
    }

    int u = low + index / 2;

    DandC(pRec, pts, low, u - 1, distance);
    DandC(pRec, pts, u, high, distance);

    int k1 = 0, k2 = 0;

    // note pts[u].x is middle x
    for(int i = u - 1; i >= low && pts[i].x >= (pts[u].x - distance); --i, ++k1);

    if(k1 == 0)
        return;

    for(int i = u; i <= high && pts[i].x <= (pts[u].x + distance); ++i, ++k2);

    std::vector<IndexAndSide> pts3(k1 + k2);

    for(int i = 0; i < k1; ++i) {
        pts3[i].N = u - 1 - i;
        pts3[i].R = false;
    }

    for(int i = 0; i < k2; ++i) {
        pts3[i + k1].N = u + i;
        pts3[i + k1].R = true;
    }

    QuickSortY(pts, &pts3.front(), 0, k1 + k2 - 1);

    // look at points in array from 0 to k1 + k2
    for(int i = 0; i < k1 + k2 - 1; ++i) {
        AcGePoint3d pt1 = pts[pts3[i].N];
        // look at points in array from i + 1 to k1 + k2
        for(int j = i + 1; j < k1 + k2; ++j) {
            // compare the points to see if they are within distance of each other
            AcGePoint3d pt2 = pts[pts3[j].N];
            if(fabs(pt2.y - pt1.y) > distance)
                break;
            if(fabs(pt2.x - pt1.x) <= distance && fabs(pt2.y - pt1.y) <= distance &&
                pt1.distanceTo(pt2) <= distance &&
                (pts3[i].R ^ pts3[j].R)) {
                    // yep, build a line
                    AddLine(pRec, pt1, pt2);
            }
        }
    }
}

//////////////////////////////////////////////////////////////////////////
//* Coolect all the points in the block table record reference in pRec *//
//////////////////////////////////////////////////////////////////////////
size_t CollectPoints(AcDbBlockTableRecord * pRec, std::vector<AcGePoint3d> & points)
{
    AcDbBlockTableRecordIterator *pIter;
    if(pRec->newIterator(pIter) == Acad::eOk)
    {
        for(pIter->start(); !pIter->done(); pIter->step())
        {
            AcDbEntity *pEnt = NULL;
            if(pIter->getEntity(pEnt, AcDb::kForRead) == Acad::eOk)
            {
                AcDbPoint *pPoint = AcDbPoint::cast(pEnt);
                if(pPoint)
                    points.push_back(pPoint->position());
                pEnt->close();
            }
        }
        delete pIter;
    }
    return points.size();
}

/////////////////////////////////////////////////////////////////////////////////
//* A very fast vertex sort routine written by Jonathan Richard Shewchuk      *//
//* and is part of his "Triangle, A Two-Dimensional Quality Mesh Generator    *//
//* and Delaunay Triangulator, copyright 2005                                 *//
//////                                                                     //////
//*  This is part of the insertvertex routine starting at line 6507 - 6525    *//
//* of version 1.5 of trinalge.  More specific it from about line 7707 - 7822 *//
/////////////////////////////////////////////////////////////////////////////////
unsigned long randomseed = 0;
size_t randomnation(size_t choices)
{

    randomseed = (randomseed * 13661 + 1508891) % 7140251;
    return randomseed / (7140251 / choices + 1);
}

void VertexSortX(AcGePoint3d * sortarray, int arraysize)
{
    int left, right;
    int pivot;
    double pivotx, pivoty;


    if (arraysize == 2) {
        /* Recursive base case. */
        if ((sortarray[0].x > sortarray[1].x) ||
            ((sortarray[0].x == sortarray[1].x) &&
            (sortarray[0].y > sortarray[1].y))) {
                std::swap(sortarray[1], sortarray[0]);
        }
    }
    /* Choose a random pivot to split the array. */
    pivot = (int) randomnation((unsigned int) arraysize);
    pivotx = sortarray[pivot].x;
    pivoty = sortarray[pivot].y;
    /* Split the array. */
    left = -1;
    right = arraysize;
    while (left < right) {
        /* Search for a vertex whose x-coordinate is too large for the left. */
        do {
            left++;
        } while ((left <= right) && ((sortarray[left].x < pivotx) ||
            ((sortarray[left].x == pivotx) &&
            (sortarray[left].y < pivoty))));
        /* Search for a vertex whose x-coordinate is too small for the right. */
        do {
            right--;
        } while ((left <= right) && ((sortarray[right].x > pivotx) ||
            ((sortarray[right].x == pivotx) &&
            (sortarray[right].y > pivoty))));
        if (left < right) {
            /* Swap the left and right vertices. */
            std::swap(sortarray[left], sortarray[right]);
        }
    }

    if (left > 1) {
        /* Recursively sort the left subset. */
        VertexSortX(sortarray, left);
    }
    if (right < arraysize - 2) {
        /* Recursively sort the right subset. */
        VertexSortX(&sortarray[right + 1], arraysize - right - 1);
    }
}

/////////////////////////////////////////////////////////////////
//* This is where all the work begins getting user input,     *//
//* collecting and sorting the points then finally qualifying *//
//* the points. Oh ya, and the all important timer reporting  *//
/////////////////////////////////////////////////////////////////
void ConnectPoints(void)
{
    double distance = 0.0;

    if(acedGetDist(0,_T("\nGetDist: "),&distance) != RTNORM)
        return;

    const double distancesqrd = pow(distance,2);

    __int64 freq, start, end, diff;
    QueryPerformanceFrequency((LARGE_INTEGER*)&freq);
    QueryPerformanceCounter((LARGE_INTEGER*)&start);

    std::vector<AcGePoint3d> points;
    AcDbBlockTableRecord *pRec;
    AcDbDatabase *pDb = acdbHostApplicationServices()->workingDatabase();

    if(acdbOpenAcDbObject((AcDbObject *&) pRec,pDb->currentSpaceId(), AcDb::kForWrite) == Acad::eOk)
    {           
        CollectPoints(pRec, points);
        VertexSortX(&points.front(), points.size());

        // QueryPerformanceCounter((LARGE_INTEGER*)&end);
        // acutPrintf(_T("\n%g Seconds collecting and sorting points\n"), (double)(end - start) / (double) freq);

        DandC(pRec, &points.front(), 0, points.size(), distance);

        pRec->close();
    }

    QueryPerformanceCounter((LARGE_INTEGER*)&end);
    acutPrintf(_T("\n%g Seconds\n"), (double)(end - start) / (double) freq);
}


[highflyingbird]

Pkonut,you are very good man,I like to talk to you.
you gave me  some very good advice and a big challenge,I am happy  to accept  it.
Today I am very busy, so I won't be able to take a look,tomorrow,I am free.
BTW,I have a msn, highflybird@msn,would you make friend with me ?

[pkohut]

you gave me  some very good advice and a big challenge,I am happy  to accept  it.

Well I just want to say thanks for sharing your research and code, that approach to the problem is very interesting.

The very best way to contact me through email (see profile).

[highflyingbird]

I analyzed my algorithm  carefully a few days ago, I found  it takes O(n*log(n) *log(n) + K) time,This is in the worst situation. and I learned the Algorithm of  KDtree, it take  O (n*sqrt(n)+k) time ,so I think my algorithm is better than the algorithm of KDtree ,  ,maybe  my analyzing is wrong. 
Yesterday, I thought a better way to improver my code, I will realize it  after I finish my work.
pkohut, your email address is hidden.

[pkohut]

I analyzed my algorithm  carefully a few days ago, I found  it takes O(n*log(n) *log(n) + K) time,This is in the worst situation. and I learned the Algorithm of  KDtree, it take  O (n*sqrt(n)+k) time ,so I think my algorithm is better than the algorithm of KDtree ,  ,maybe  my analyzing is wrong. 
Yesterday, I thought a better way to improver my code, I will realize it  after I finish my work.
pkohut, your email address is hidden.

Setting up the KDTree is about as fast as doing a quick sort.  I was working on a remove node function, but kept muffing it up, which I think once implemented would allow it to perform at least as fast as Daniel's code (gut feeling is it would smoke Daniel's entry).  But your algorithm will be faster for this challenge.

Code: [Select]

   while kd-tree has points {
         process root node
         remove root node
    }
This would keep the total iterations to an absolute minimum.  As it is now it must visit work that has already been completed in order to find a leaf that needs to be worked on, which takes time.

But, KD-tree is a spatial container for N-dimensional objects and allows some types of spatial queries, especially if the data to queried is stored on a drive, network, or even in a SQL database.  I use the KDTree code presented here in a few programs , and it's generic enough for easy reuse. What's very nice about the KDTree routine presented is that each leaf node can contain user defined data.

When trying to speed up my entry I used the extra data as a flag indicating the leaf has been previously processed (speedup was only a couple seconds over my best posted speed).

For this challenge I wanted to show a different way to accomplish the work.  Because my best posted speed of 29 seconds (about 20 seconds raw speed) people may think its not very good, but it's pretty much left in its generic state. Adding the remove node will speed it up dramatically, and tweaking the code a bit for this challenge could get a little more speed.

My email address is paulkohut at hotmail dot com

[pkohut]

<...snip...>
Setting up the KDTree is about as fast as doing a quick sort.  I was working on a remove node function, but kept muffing it up, which I think once implemented would allow it to perform at least as fast as Daniel's code (gut feeling is it would smoke Daniel's entry).
<...snip...>
For this challenge I wanted to show a different way to accomplish the work.  Because my best posted speed of 29 seconds (about 20 seconds raw speed)

Erase node implemented and first pass clocks in at 12.24 seconds raw speed. Looks like I'll be eatin' gut and breathing Daniel's smoke  :kewl:

[Kerry]

That statement may be a little difficult for 'bing' to translate (sensibly, anyway)

[pkohut]

That statement may be a little difficult for 'bing' to translate (sensibly, anyway)

Ya, that dawned on me a little while after reading your post. All of a sudden there was an ugly picture in my mind.

Without some major restructuring (beyond that CKdTree is now a template and does conditional compiling) I've got it down to 11.20 seconds RS.

[pkohut]

For the next test I change where the timer is started in your code to more closely match Daniel's original, and I subtracted 1 from your total time as an allowance to type "all" at the select objects prompt and hit return twice.
Your's  - 22.63 seconds
Daniel's - 8.24 seconds
(note: Your's is HighFlyBird's code)

And finally enable the AcDbLine creation.
Your's  - 31.48 seconds.
Daniel's - 18.90 seconds  (oops forgot to cache a value, now down to 16.4)

New numbers for my KDTree entry 
Prior to my final optimization it was running raw speed at 11.20. Could not get it any faster unless I turned on SSE2 and favor speed as a compiler option which then dropped it to 10.83 seconds raw speed for the million point test.  I tried 5 or 6 methods of tracking previous nodes and such and non provided any speed improvement, two methods I thought would work for sure added about 5 seconds to the elapse time (making the starting point in the tree the current node position, then walk backward up the tree till a candidate is found, current method starts at the root node and walks forward).  Which got me thinking, is my KDTree Optimized, in other words balanced and the answer to that was, not in a million freaking years.

So, after implementing a tree optimization routine the final numbers are:
6.15 seconds raw speed, with the 1 million point test. (SSE2 off, favor speed or size: neither, 5.91 with them on)
16.46 seconds adding AcDbLine's to the drawing.

For kicks I turned off the routine that erases previously visited points (to reduce redundant work) and it runs 7.99 seconds RS, prior to that it was in the 15 second range.

I've got to do some code clean up then I'll post. Probably after 9pm PST Tuesday.

2008 ARX attached. Load it and test drawing with points, type "doit" at the command prompt.

[pkohut]

So, after implementing a tree optimization routine the final numbers are:
6.15 seconds raw speed, with the 1 million point test.

While documenting the code found a spot to further optimize by delaying when an expensive set of calcs is done.
Now down to 5.79 seconds. 

There's 14 functions left to document, once finished will release code. Got a busy day today, so shooting for midnight tonight PST.

[It's Alive!]

Paul try this arx (07) and see if there is any difference, note: its not creating the lines
the command is doit

[It's Alive!]

Hey Highflybird,

I've taken the liberty to use your QuickSortY ...

Woa! this one is fast! 

[pkohut]

Paul try this arx (07) and see if there is any difference, note: its not creating the lines
the command is doit

Ya, thats the raw speed test, it reports how many lines would have been created.

[It's Alive!]

Paul try this arx (07) and see if there is any difference, note: its not creating the lines
the command is doit

Ya, thats the raw speed test, it reports how many lines would have been created.

Sorry, I just compiled this.. I wanted to know if there was any difference on your machine.
...This is before I tested the DandC version which crushes mine so there is no need.... 

[highflyingbird]

         
            if(fabs(pt2.x - pt1.x) <= distance && fabs(pt2.y - pt1.y) <= distance &&
                pt1.distanceTo(pt2) <= distance &&
                (pts3.R ^ pts3[j].R)) {
                    // yep, build a line
                    AddLine(pRec, pt1, pt2);
            }

I think these codes can be improved.
           if(fabs(pt2.x - pt1.x) <= distance &&
               //fabs(pt2.y - pt1.y) <= distance &&         //This line can be cut. because
                                                                         // we have "if(fabs(pt2.y - pt1.y) > distance) break; "    
                (pts3.R ^ pts3[j].R)  &&                                  
                pt1.distanceTo(pt2) <= distance ){         //The distance calculation is more expensive.
                    // yep, build a line
                    AddLine(pRec, pt1, pt2);
            }
I found an another way for this topics
that is called "layered range tree"

http://cgi.di.uoa.gr/~vfisikop/compgeom/Layered_Range_trees/doc_tex/report.pdf
use this algorithm ,we can get the time O(n*log(n)) .  I wanted to change it to arx,but I failed.