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