(vl-load-com)
;;stdlib
(defun std-sslist (ss / n lst)
(if (eq 'pickset (type ss))
(repeat (setq n (fix (sslength ss))) ; fixed
(setq lst (cons (ssname ss (setq n (1- n))) lst))
)
)
)
;;make layer or turn on current layer
(defun mylayer (str)
(if (tblsearch "LAYER" str)
(command "._-Layer" "_set" str "")
(command "._-Layer" "_Make" str "")
)
)
;;
;;entmake line
(defun entmakeline (pt1 pt2 layer)
(entmake (list (cons 0 "LINE") ;***
(cons 6 "BYLAYER") (cons 8 layer) (cons 10 pt1) ;***
(cons 11 pt2) ;***
(cons 39 0.0) (cons 62 256) (cons 210 (list 0.0 0.0 1.0))
)
)
);;
;;get the point coordinate
(defun mygetpoint (lst / plst)
(foreach x lst
(setq plst (append
plst
(list (cdr (assoc 10 (entget x))))
)
)
)
plst
)
;;Main
(defun c:gap ( / c i j len len1 limdis pj plst x)
(mylayer "temp")
(setq limdis 70)
(setq c (ssget '((0 . "POINT"))))
(setq c (std-sslist c))
(setq plst (mygetpoint c))
(setq i 0
len (length plst)
)
(foreach x plst
(princ "\n i = ") (princ i);Erase_DV
(setq j (1+ i)
len1 (- len j)
)
(repeat len1
(setq pj (nth j plst))
(if (< (distance x pj) limdis)
(entmakeline x pj "temp")
)
(setq j (1+ j))
)
(setq i (1+ i))
)
)
(defun c:test (/ d e n out pt ss)
(and
(setq n -1)
(setq d 70.)
(setq ss (ssget '((0 . "point"))))
(while (setq e (ssname ss (setq n (1+ n))))
(setq out (cons (cdr (assoc 10 (entget e))) out))
)
(setq out (vl-sort out (function (lambda (x1 x2) (< (car x1) (car x2))))))
(while (setq pt (car out))
(mapcar (function
(lambda (x) (entmakex (list '(0 . "LINE") (cons 10 pt) (cons 11 x))))
)
(vl-remove-if
(function (lambda (x) (>= (distance pt x) d)))
(setq out (cdr out))
)
)
)
)
(princ)
)
(defun PointGap (gap / GetPoints Pts tmp x y)
(defun GetPoints (ss / e i l)
(setq i -1)
(if ss
(while (setq e (ssname ss (setq i (1+ i))))
(setq l (cons (cdr (assoc 10 (entget e))) l))))
l)
(if (setq Pts (GetPoints (ssget "_X" '((0 . "POINT")))))
(while (setq x (car Pts))
(setq tmp (cdr Pts))
(while (setq y (car tmp))
(if (< (distance x y) gap)
(entmakex (list (cons 0 "LINE")
(cons 10 x) (cons 11 y))))
(setq tmp (cdr tmp)))
(setq Pts (cdr Pts))))
(princ))
(_connect (_GetPnts) 70.))
(_connect2 (_GetPnts) 70.))(defun _GetPnts (/ ss i)
(setq i -1)
(if (setq ss (ssget "_X" '((0 . "POINT"))))
(while (setq e (ssname ss (setq i (1+ i))))
(setq l (cons (cdr (assoc 10 (entget e))) l))
) ;_ while
) ;_ if
) ;_ defun
(defun _Line (a b) (entmakex (list '(0 . "LINE") (cons 10 b) (cons 11 a))))
(defun _connect (lst d / lst)
(foreach x lst
(foreach i (setq lst (vl-remove x lst))
(and (< (distance x i) d)
(_Line x i)
) ;_ and
) ;_ foreach
) ;_ foreach
) ;_ defun
(defun _connect2 (lst d / lst)
(while (setq a (car lst))
(foreach x (setq lst (cdr lst))
(and (< (distance a x) d)
(_Line a x)
) ;_ and
) ;_ foreach
) ;_ while
) ;_ defun
(defun foo (dist / n ss ent lst pt)
(if (setq n -1
ss (ssget "_X" '((0 . "POINT")))
)
(while (setq ent (ssname ss (setq n (1+ n))))
(setq lst (cons (cdr (assoc 10 (entget ent))) lst))
)
)
(while lst
(setq pt (car lst)
lst (cdr lst)
)
(foreach p lst
(if (<= (distance pt p) dist)
(entmake (list '(0 . "LINE") (cons 10 pt) (cons 11 p)))
)
)
)
)
Hi,
Here's my way, quite similar to those posted before.Code: [Select]<=
using System;
using Autodesk.AutoCAD.ApplicationServices;
using Autodesk.AutoCAD.DatabaseServices;
using Autodesk.AutoCAD.EditorInput;
using Autodesk.AutoCAD.Geometry;
using Autodesk.AutoCAD.Runtime;
namespace ChallengeConnectPoints
{
public class Connect
{
[CommandMethod("Test")]
public void Test()
{
Editor ed = Application.DocumentManager.MdiActiveDocument.Editor;
PromptDistanceOptions pdo = new PromptDistanceOptions("\nConnection distance: ");
pdo.AllowNegative = false;
pdo.AllowZero = false;
PromptDoubleResult pdr = ed.GetDistance(pdo);
if (pdr.Status == PromptStatus.OK)
ConnectPoints(pdr.Value);
}
private void ConnectPoints(double d)
{
Document doc = Application.DocumentManager.MdiActiveDocument;
Database db = doc.Database;
Editor ed = doc.Editor;
TypedValue[] tv = new TypedValue[1] { new TypedValue(0, "POINT") };
PromptSelectionResult psr = ed.SelectAll(new SelectionFilter(tv));
if (psr.Status != PromptStatus.OK)
return;
using (Transaction tr = db.TransactionManager.StartTransaction())
{
Point3d[] pts = Array.ConvertAll<ObjectId, Point3d>(
psr.Value.GetObjectIds(), x => ((DBPoint)x.GetObject(OpenMode.ForRead)).Position);
int cnt = pts.Length;
BlockTableRecord btr =
(BlockTableRecord)tr.GetObject(db.CurrentSpaceId, OpenMode.ForWrite);
for (int i = 0; i < cnt; i++)
{
for (int j = i + 1; j < cnt; j++)
{
if (pts[i].DistanceTo(pts[j]) <= d)
{
Line l = new Line(pts[i], pts[j]);
btr.AppendEntity(l);
tr.AddNewlyCreatedDBObject(l, true);
}
}
}
tr.Commit();
}
}
}
}
Hi,
A C# solution
Much more faster than LISP. With the sample drawing and a 70.0 length:
LISP = about 45 seconds
C# = about 2 seconds
Maybe ARX shoulde be faster...
Have an hour or so of coding left to do.Man, one really pays for performance. I'm not knocking C/Arx, just taken aback at the amount of time.
Have an hour or so of coding left to do.Man, one really pays for performance. I'm not knocking C/Arx, just taken aback at the amount of time.
:lol:Have an hour or so of coding left to do.Man, one really pays for performance. I'm not knocking C/Arx, just taken aback at the amount of time.
Your right, and this is even repurposing code I've already written. The solution though is non-trivial, but screams. Tested with a million points and it's under a second. Still coding though so back to work :lmao:
:lol:Have an hour or so of coding left to do.Man, one really pays for performance. I'm not knocking C/Arx, just taken aback at the amount of time.
Your right, and this is even repurposing code I've already written. The solution though is non-trivial, but screams. Tested with a million points and it's under a second. Still coding though so back to work :lmao:
Can't wait to see it!!
static void ArxDist_doit(void)
{
double distance = 0.0;
acedGetDist(0,ACRX_T("\nGetDist: "),&distance);
__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)
{
AcDbBlockTableRecordIterator *pIter = NULL;
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;
}
size_t len = points.size();
for (size_t i = 0; i < len; i++)
{
for (size_t j = i + 1; j < len; j++)
{
if ( (points[i].distanceTo(points[j])) <= distance)
{
AcDbLine *pLine = new AcDbLine(points[i],points[j]);
if(pRec->appendAcDbEntity(pLine) != Acad::eOk)
delete pLine;
pLine->close();
}
}
}
pRec->close();
}
QueryPerformanceCounter((LARGE_INTEGER*)&end);
acutPrintf(_T("\n%g Seconds\n"), (double)(end - start) / (double) freq);
}
Thank you all :), Your algorithms are all faster than me.
In my PC, (E7500 CPU, 4G Mem)
Ronjonp 's code, about 34.812967s
Lee Mac's code, about 40.733018s
Gile's code 42.111982s
alanjt's code is about 50s
Mine first code is very slow. :-D, (could you tell me which cause the very slow)
Not very accurate test.
And I hope to see pkohut's super code :)
I am thinking about the "Divide and conquer algorithm", but I need some time to sort the point in some manner, and not so sure about whether it can be improved, because the sort also need some times.
brute force with arx
doit
GetDist: 70
1.036883 Seconds
arx for 07
brute force with arx
doit
GetDist: 70
1.036883 Seconds
arx for 07
That's incredible Daniel! I have got to put Lisp down and learn C.
connectdotsYou will them be prompted for the maximum segment length to connect, and to select the points to operate on.CKdTree kdTree;
kdTree.Create(3);
kdTree.Insert(pt3d.x, pt3d.y, pt3d.z, NULL);RESULT_SET * pResAllPoints = kdTree.NearestRange(0.0, 0.0, 0.0, HUGE_VAL);
...
double * pPos = pResAllPoints->riter->item->pos;
RESULT_SET * pRes = kdTree.NearestRange(pPos[0], pPos[1], 0.0, dRadius);Have an hour or so of coding left to do.Man, one really pays for performance. I'm not knocking C/Arx, just taken aback at the amount of time.
Your right, and this is even repurposing code I've already written. The solution though is non-trivial, but screams. Tested with a million points and it's under a second. Still coding though so back to work :lmao:
(vl-load-com)
;;;start of Time, from theswamp
(defun startTimer ()
(setq time (getvar "DATE"))
)
;;;end of Time, from theswamp
(defun endTimer (/ seconds totaltime)
(setq time (- (getvar "DATE") time)
seconds (* 86400.0 (- time (fix time)))
)
(setq totaltime (rtos seconds 2 6))
(princ (strcat "\nTimed:" (rtos seconds 2 6)))
)
;;entmake line
(defun entmakeline (pt1 pt2 layer)
(entmake (list (cons 0 "LINE") ;***
(cons 6 "BYLAYER") (cons 8 layer) (cons 10 pt1) ;***
(cons 11 pt2) ;***
(cons 39 0.0) (cons 62 256) (cons 210 (list 0.0 0.0 1.0))
)
)
)
;;The idea of gile
(defun getpointscor ( / ent n plst ss)
(if (setq n -1
ss (ssget '((0 . "POINT")))
)
(while (setq ent (ssname ss (setq n (1+ n))))
(setq plst (cons (cdr (assoc 10 (entget ent))) plst))
)
)
)
(defun dividex (plst dis / a b div end index len lst1 plstmp res start xend
xstart
)
(setq start (car plst)
xstart (car start)
end (last plst)
len (- (car end) (car start))
div (+ (fix (/ len dis)) 1)
plstmp plst
res nil
)
(repeat div
(grdraw (list xstart 0 0) (list xstart 10000 0) 8)
(setq xend (+ xstart dis))
(setq index T
lst1 nil
)
(while (and
plstmp
index
)
(if (< (setq a (car plstmp)
b (car a)
)
xend
)
(setq lst1 (append
lst1
(list a)
)
plstmp (cdr plstmp)
)
(setq index nil)
)
)
(setq res (append
res
(list lst1)
)
)
(setq xstart (+ xstart dis))
)
res
)
;;;draw short line in each x region
;;;The idea of Leemac and gile
(defun drawshortline (plst limdis / a b tmp)
(while (setq a (car plst))
(setq plst (cdr plst)
tmp plst
)
(while (setq b (car tmp))
(setq tmp (cdr tmp))
(if (< (distance a b) limdis)
(entmakeline a b "temp")
)
)
)
)
;;;draw short line in each two Neighbor x region
(defun drawshortline12 (plst1 plst2 limdis / a b x)
(while (setq a (car plst1))
(foreach x plst2
(if (< (distance a x) limdis)
(entmakeline a x "temp")
)
)
(setq plst1 (cdr plst1))
)
)
;;Connect short line between each two points
;;by qjchen@gmail.com
;; Thanks to ronjonp, Lee Mac, and Gile
;;Version 1.0
(defun c:gap ( / a b limdis plst plst1) ;(mylayer "temp")
(starttimer)
(setq limdis 70)
;;The idea of gile
(setq plst (getpointscor))
;;The idea of ronjonp
(setq plst (vl-sort plst (function (lambda (x1 x2)(< (car x1) (car x2))))))
;;Divde points list to several region
;;such like ((0, contain x=0~70),(1, contain x=70~140)........)
(setq plst1 (dividex plst limdis))
;;Draw shortline in each x region ,such like (0, contain x=0~70))
(foreach x plst1
(drawshortline x limdis)
)
;;Draw shortline between two Neighbor x region, such between (0, contain x=0~70),(1, contain x=70~140))
(while (setq a (car plst1) b (cadr plst1) plst1 (cdr plst1))
(drawshortline12 a b limdis)
)
(endtimer)
)
brute force with arx
doit
GetDist: 70
1.036883 Seconds
arx for 07
That's incredible Daniel! I have got to put Lisp down and learn C.
:)
Lisp can't be 'put down' it just refuses to die :-D
I agree about the idea of learning C, C++. Great if we can make the time and if the wet stuff is still receptive.
(defun el (d / A L LL TI X)
;;(el 70.)
(setq ti (car (_VL-TIMES))
a -1
ll (ssget "_x" '((0 . "point")))
)
(while (setq x (ssname ll (setq a (1+ a))))
(setq l (cons (cdr (assoc 10 (entget x))) l))
)
(setq l (vl-sort l (function (lambda (a b) (<= (car a) (car b))))))
(while (setq a (car l)
l (cdr l)
)
(setq x (car a)
ll l
)
(while (and ll (<= (- (caar ll) x) d))
(if (and (equal (cdar ll) (cdr a) d) (<= (distance a (car ll)) d))
(entmakex (list '(0 . "LINE") (cons 10 a) (cons 11 (car ll))))
)
(setq ll (cdr ll))
)
)
(princ (strcat "\n " (rtos (/ (- (car (_VL-TIMES)) ti) 1000.) 2 4) "sec."))
(princ)
)
using System;
using System.Collections.Generic;
using System.Text;
using System.Threading;
using System.Threading.Tasks;
using System.Diagnostics;
using Autodesk.AutoCAD.GraphicsInterface;
using Autodesk.AutoCAD.Runtime;
using Autodesk.AutoCAD.Windows;
using Autodesk.AutoCAD.Geometry;
using Autodesk.AutoCAD.EditorInput;
using Autodesk.AutoCAD.LayerManager;
using Autodesk.AutoCAD.ApplicationServices;
using Autodesk.AutoCAD.DatabaseServices;
using Autodesk.AutoCAD.Colors;
using System.Threading.Tasks;
using AcAp = Autodesk.AutoCAD.ApplicationServices;
using AcEd = Autodesk.AutoCAD.EditorInput;
using AcGe = Autodesk.AutoCAD.Geometry;
using AcRx = Autodesk.AutoCAD.Runtime;
using AcDb = Autodesk.AutoCAD.DatabaseServices;
using AcWd = Autodesk.AutoCAD.Windows;
namespace Test
{
public class ExecMethod
{
[CommandMethod("doit")]
public static void MyCommand()
{
Stopwatch stopWatch = new Stopwatch();
stopWatch.Start();
SynchronizedCollection<Line>
lines = new SynchronizedCollection<Line>();
List<Point3d> points = getPoints();
int cnt = points.Count;
Parallel.For(0, cnt, delegate(int i)
{
Point3d pnt = points[i];
for (int j = i + 1; j < cnt; j++)
{
if (pnt.DistanceTo(points[j]) <= 70.0)
{
lines.Add(new Line(pnt, points[j]));
}
}
});
addtoDb(lines);
stopWatch.Stop();
TimeSpan ts = stopWatch.Elapsed;
AcAp.Application.DocumentManager.MdiActiveDocument.
Editor.WriteMessage("\n{0} Seconds", ts.TotalMilliseconds);
}
static void addtoDb(SynchronizedCollection<Line> lines)
{
List<Point3d> points = new List<Point3d>();
Database database = HostApplicationServices.WorkingDatabase;
using (Transaction transaction = database.TransactionManager.StartTransaction())
{
BlockTableRecord record =
transaction.GetObject(
database.CurrentSpaceId, OpenMode.ForWrite) as BlockTableRecord;
foreach (Line line in lines)
{
record.AppendEntity(line);
transaction.AddNewlyCreatedDBObject(line, true);
}
transaction.Commit();
}
}
static List<Point3d> getPoints()
{
List<Point3d> points = new List<Point3d>();
Database database = HostApplicationServices.WorkingDatabase;
using (Transaction transaction = database.TransactionManager.StartTransaction())
{
BlockTableRecord record =
transaction.GetObject(
database.CurrentSpaceId, OpenMode.ForRead) as BlockTableRecord;
foreach (ObjectId id in record)
{
DBPoint point = transaction.GetObject(id, OpenMode.ForRead) as DBPoint;
if (point != null)
{
points.Add(point.Position);
}
}
transaction.Commit();
}
return points;
}
}
}
Command: doit:lol: :lol:
497.0208 Milliseconds 8-)
my variant:Code: [Select](defun el (d / A L LL TI X)
;;(el 70.)
(setq ti (car (_VL-TIMES))
a -1
ll (ssget "_x" '((0 . "point")))
)
(while (setq x (ssname ll (setq a (1+ a))))
(setq l (cons (cdr (assoc 10 (entget x))) l))
)
(setq l (vl-sort l (function (lambda (a b) (<= (car a) (car b))))))
(while (setq a (car l)
l (cdr l)
)
(setq x (car a)
ll l
)
(while (and ll (<= (- (caar ll) x) d))
(if (and (equal (cdar ll) (cdr a) d) (<= (distance a (car ll)) d))
(entmakex (list '(0 . "LINE") (cons 10 a) (cons 11 (car ll))))
)
(setq ll (cdr ll))
)
)
(princ (strcat "\n " (rtos (/ (- (car (_VL-TIMES)) ti) 1000.) 2 4) "sec."))
(princ)
)
my variant:Code: [Select](defun el (d / A L LL TI X)
;;(el 70.)
(setq ti (car (_VL-TIMES))
a -1
ll (ssget "_x" '((0 . "point")))
)
(while (setq x (ssname ll (setq a (1+ a))))
(setq l (cons (cdr (assoc 10 (entget x))) l))
)
(setq l (vl-sort l (function (lambda (a b) (<= (car a) (car b))))))
(while (setq a (car l)
l (cdr l)
)
(setq x (car a)
ll l
)
(while (and ll (<= (- (caar ll) x) d))
(if (and (equal (cdar ll) (cdr a) d) (<= (distance a (car ll)) d))
(entmakex (list '(0 . "LINE") (cons 10 a) (cons 11 (car ll))))
)
(setq ll (cdr ll))
)
)
(princ (strcat "\n " (rtos (/ (- (car (_VL-TIMES)) ti) 1000.) 2 4) "sec."))
(princ)
)
;;from internet
(defun random ()
(setq seed (if seed
(rem (+ (* seed 15625.7) 0.21137152) 1)
0.3171943
)
)
)
(defun random-n (n)
(* n (random))
)
(defun entmakepoint (pt layer)
(entmake (list (cons 0 "POINT")
(cons 8 layer);***
(cons 6 "BYLAYER")
(cons 10 pt) ;***
(cons 39 0.0)
(cons 50 0.0)
(cons 62 256)
(cons 210 (list 0.0 0.0 1.0))
)
)
)
;;generate random points
(defun c:test1 ()
(repeat 100000
(entmakepoint (list (random-n 50000) (random-n 30000) 0) "0")
)
)
using System;
using System.Collections.Generic;
using Autodesk.AutoCAD.ApplicationServices;
using Autodesk.AutoCAD.DatabaseServices;
using Autodesk.AutoCAD.EditorInput;
using Autodesk.AutoCAD.Geometry;
using Autodesk.AutoCAD.Runtime;
namespace ChallengeConnectPoints
{
public class Connect
{
[CommandMethod("Test")]
public void Test()
{
Editor ed = Application.DocumentManager.MdiActiveDocument.Editor;
PromptDistanceOptions pdo = new PromptDistanceOptions("\nConnection distance: ");
pdo.AllowNegative = false;
pdo.AllowZero = false;
PromptDoubleResult pdr = ed.GetDistance(pdo);
if (pdr.Status == PromptStatus.OK)
ConnectPoints(pdr.Value);
}
private void ConnectPoints(double d)
{
DateTime t0 = DateTime.Now;
Document doc = Application.DocumentManager.MdiActiveDocument;
Database db = doc.Database;
Editor ed = doc.Editor;
List<Point3d> pts = new List<Point3d>();
using (Transaction tr = db.TransactionManager.StartTransaction())
{
BlockTableRecord btr =
(BlockTableRecord)tr.GetObject(db.CurrentSpaceId, OpenMode.ForWrite);
foreach (ObjectId id in btr)
{
DBPoint pt = tr.GetObject(id, OpenMode.ForRead) as DBPoint;
if (pt != null)
pts.Add(pt.Position);
}
pts.Sort((p1,p2) => p1.X.CompareTo(p2.X));
int cnt = pts.Count;
for (int i = 0; i < cnt; i++)
{
Point3d p1 = pts[i];
double x1 = p1.X;
for (int j = i+1; j < cnt; j++)
{
Point3d p2 = pts[j];
if ((p2.X - x1) > d)
break;
if (p1.DistanceTo(p2) <= d)
{
Line l = new Line(p1, pts[j]);
btr.AppendEntity(l);
tr.AddNewlyCreatedDBObject(l, true);
}
}
}
tr.Commit();
}
DateTime t1 = DateTime.Now;
TimeSpan ts = t1 - t0;
ed.WriteMessage("\nEllapsed milliseconds : {0}", ts.TotalMilliseconds);
}
}
}
Intel Core 2 Q6600 @ 2.4GHz
3.25 GB Ram
Tested Ac2008
Gile's
10,000 pts - 328 ms - 0.328s
100,000 pts - 4859 ms - 4.86s
500,000 pts - 26015 ms - 26.02s
1,000,000 pts - 84500 ms - 84.5s
Paul's
10,000 pts - 78 ms - 0.078s
100,000 pts - 1297 ms - 1.30s
500,000 pts - 14766 ms - 14.77s
1,000,000 pts - 34625 ms - 34.63s
Daniel's
Can't recreate because I don't have MSVS 2010/.Net 4
also note his reported 497 ms is on his smoking hot i7 machine
Paul, Try this ARX , the command is doit... using Evgeniy's optimization
Times using Evgeniy's optimization, .NET 4, running in parallel.
Just a note, but the 50% of the time is spent adding the lines to the DB as this procedure cannot run in parallel
10,000 - 106.7622 Milliseconds
100,000 - 1127.0502 Milliseconds
1,000,000 - 18485.5115 Milliseconds
Times using Evgeniy's optimization, .NET 4, running in parallel.
Just a note, but the 50% of the time is spent adding the lines to the DB as this procedure cannot run in parallel
10,000 - 106.7622 Milliseconds
100,000 - 1127.0502 Milliseconds
1,000,000 - 18485.5115 Milliseconds
Typo ?
Paul, Try this ARX , the command is doit... using Evgeniy's optimization
Nice numbers! Code please :-D
10,000 - 0.0637 s
100,000 - 1.7528 s
500,000 - 10.215 s
1,000,000 - 36.993 s
static void ArxDist_doit(void)
{
double distance = 0.0;
if(acedGetDist(0,ACRX_T("\nGetDist: "),&distance) != RTNORM)
return;
__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)
{
AcDbBlockTableRecordIterator *pIter = NULL;
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;
}
std::sort(points.begin(),points.end(),psort);
size_t len = points.size();
for (size_t i = 0; i < len; i++)
{
for (size_t j = i + 1; j < len; j++)
{
if ((points[j].x - points[i].x) > distance)
{
break;
}
else if ( (points[i].distanceTo(points[j])) <= distance)
{
AcDbLine *pLine = new AcDbLine(points[i],points[j]);
if(pRec->appendAcDbEntity(pLine) != Acad::eOk)
delete pLine;
pLine->close();
}
}
}
pRec->close();
}
QueryPerformanceCounter((LARGE_INTEGER*)&end);
acutPrintf(_T("\n%g Seconds\n"), (double)(end - start) / (double) freq);
}
(vl-load-com)
;;entmake line
(defun entmakeline (pt1 pt2 layer)
(entmake (list (cons 0 "LINE") ;***
(cons 6 "BYLAYER") (cons 8 layer) (cons 10 pt1) ;***
(cons 11 pt2) ;***
(cons 39 0.0) (cons 62 256) (cons 210 (list 0.0 0.0 1.0))
)
)
)
;;The idea of gile
(defun getpointscor (/ ent n plst ss)
(if (setq n -1
ss (ssget '((0 . "POINT")))
)
(while (setq pubtime (car (_VL-TIMES)) ent (ssname ss (setq n (1+ n))))
(setq plst (cons (cdr (assoc 10 (entget ent))) plst))
)
)
)
;;the bounding box of the points set
(defun getboundingboxpointset (a)
(list (apply 'min (mapcar '(lambda (x) (car x)) a))
(apply 'max (mapcar '(lambda (x) (car x)) a))
(apply 'min (mapcar '(lambda (x) (cadr x)) a))
(apply 'max (mapcar '(lambda (x) (cadr x)) a))
)
)
;;;divide the points set according to a certain x or y or z space
(defun dividexyz (plst dis vmin vmax stringxyz / a b div end index lst1 res start)
(setq start vmin
div (+ (fix (/ (- vmax vmin) dis)) 1)
res nil
)
(repeat div
; (grdraw (list start 0 0) (list start 10000 0) 8)
(setq end (+ start dis)
index T
lst1 nil
)
(while (and plst index)
(if (< (setq a (car plst)
b (cond ((= stringxyz "x") (car a))
((= stringxyz "y") (cadr a))
((= stringxyz "z") (caddr a))
(T nil)
)
)
end
)
(setq lst1 (append
lst1
(list a)
)
plst (cdr plst)
)
(setq index nil)
)
)
(setq res (append
res
(list lst1)
)
)
(setq start (+ start dis))
)
res
)
;;;draw short line in each x region
;;;The idea of Leemac and gile
(defun drawshortline (plst limdis / a b tmp)
(if plst
(while (setq a (car plst))
(setq plst (cdr plst)
tmp plst
)
(while (setq b (car tmp))
(setq tmp (cdr tmp))
(if (< (distance a b) limdis)
(entmakeline a b "temp")
)
)
)
)
)
;;;draw short line in each two Neighbor x region
(defun drawshortline12 (plst1 plst2 limdis / a b x)
(if (and
plst1
plst2
)
(while (setq a (car plst1))
(foreach x plst2
(if (< (distance a x) limdis)
(entmakeline a x "temp")
)
)
(setq plst1 (cdr plst1))
)
)
)
;;draw line between two pointsets and manner
(defun dl12(a b limdis stringabc / ab lst1 lst2)
(cond ((= stringabc "90") (setq lst1 (reverse (cdr (reverse a))) lst2 (cdr b)))
((= stringabc "0") (setq lst1 a lst2 b))
((= stringabc "-45") (setq lst1 (cdr a) lst2 (reverse (cdr (reverse b)))))
((= stringabc "45") (setq lst1 (reverse (cdr (reverse a))) lst2 (cdr b)))
(T nil)
)
(setq ab (mapcar '(lambda (x y) (list x y)) lst1 lst2))
(foreach x ab
(if x
(drawshortline12 (car x) (cadr x) limdis)
)
)
)
;;draw line between the own pointsets
(defun dl1(pl1 limdis)
(foreach x pl1
(drawshortline x limdis)
)
)
;;Connect short line between each two points
;;by qjchen@gmail.com
;; Thanks to ronjonp, Lee Mac, and Gile and alanjt and Evgeniy
;;Version 1.0
(defun c:gap ( / a b bbox limdis plst plst1 plst2) ;(mylayer "temp")
;(starttimer)
(setq limdis 70 ) ;;The idea of gile
(setq plst (getpointscor)) ;;The idea of ronjonp
;;getting the boundingbox of all the points set
(setq bbox (getboundingboxpointset plst) )
;;;sort the point according to the x coordinate
(setq plst (vl-sort plst (function (lambda (x1 x2)
(< (car x1) (car x2))
)
)
)
) ;;Divde points list to several region according to the x value=>plst1
;;such like ((lst 0, contain x=0~70),(lst 1, contain x=70~140)........)
(setq plst1 (dividexyz plst limdis (nth 0 bbox) (nth 1 bbox) "x" )
plst2 nil
)
;;Divde each plst1 to several region according to the y value=>plst2
(foreach x plst1
(if x
(setq x (vl-sort x (function (lambda (x1 x2)
(< (cadr x1) (cadr x2))
)
)
)
plst2 (append
plst2
(list (dividexyz x limdis (nth 2 bbox) (nth 3 bbox) "y"))
)
)
(setq plst2 (append
plst2
(list nil)
)
)
)
)
;;Draw shortline in each x region,such like (0, contain x=0~70))
(while (setq a (car plst2)
b (cadr plst2)
plst2 (cdr plst2)
)
(if a
(progn
(dl1 a limdis)
(dl12 a a limdis "90")
(if b
(progn
(dl12 a b limdis "0")
(dl12 a b limdis "-45")
(dl12 a b limdis "45")
)
(progn
)
)
) ;else continue
)
)
(dl1 a limdis)
(dl12 a a limdis "90")
(princ (strcat "\n " (rtos (/ (- (car (_VL-TIMES)) pubtime) 1000.) 2 4) "sec."))
;(endtimer)
)
Paul, Try this ARX , the command is doit... using Evgeniy's optimization
Nice numbers! Code please :-D
10,000 - 0.0637 s
100,000 - 1.7528 s
500,000 - 10.215 s
1,000,000 - 36.993 s
static inline bool psort(const AcGePoint3d &a, const AcGePoint3d &b)
{
if(a.x == b.x)
return a.y < b.y;
return a.x < b.x;
}
static inline double distsqrd(const AcGePoint3d &a, const AcGePoint3d &b)
{
double dx = a.x-b.x;
double dy = a.y-b.y;
double dz = a.z-b.z;
return dx*dx+dy*dy+dz*dz;
}
static void ArxDist_doit(void)
{
double distance = 0.0;
if(acedGetDist(0,ACRX_T("\nGetDist: "),&distance) != RTNORM)
return;
const double distancesqrd = pow(distance,2) + 1.192092896e-07F;
__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)
{
AcDbBlockTableRecordIterator *pIter = NULL;
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;
}
std::sort(points.begin(),points.end(),psort);
size_t len = points.size();
for (size_t i = 0; i < len; i++)
{
for (size_t j = i + 1; j < len; j++)
{
if ((points[j].x - points[i].x) > distance)
{
break;
}
else if ( distsqrd(points[i],points[j]) < distancesqrd)
{
AcDbLine *pLine = new AcDbLine(points[i],points[j]);
if(pRec->appendAcDbEntity(pLine) != Acad::eOk)
delete pLine;
pLine->close();
}
}
}
pRec->close();
}
QueryPerformanceCounter((LARGE_INTEGER*)&end);
acutPrintf(_T("\n%g Seconds\n"), (double)(end - start) / (double) freq);
}
last one ... maybe :-D
const double distancesqrd = pow(distance,2) + 1.192092896e-07F;Eww! The fabs function is like 1 clock cycle, so your better to use it and an epsilon value of 1e-8 for your comparisons. Otherwise NICE.Sure master Evgeniy can easily improve the code and quick quick than me, but I am just try to study this funny algorithm. :-P
The algorithm is as following picture.
gile, your lambda expression to sort the points is awesome!
(vl-sort pts '(lambda (p1 p2) (< 'car p1) (car p2))))
(defun el-1 (l)
(mapcar (function (lambda (a) (cons (fix (/ (car a) d)) (cons (fix (/ (cadr a) d)) a))))
l
)
)
(defun el-2 (l / a x)
(setq x (caar l)
y (cadar l)
a (list (list (cdar l)))
l (cdr l)
)
(while l
(while (and l (= x (caar l)))
(setq a (cons (cons (cdar l) (car a)) (cdr a))
l (cdr l)
)
)
(setq y (min y (caaar a)))
(cond (l
(repeat (- x (caar l) 1) (setq a (cons (list (setq x (1- x))) a)))
(setq x (1- x)
a (cons (list (cdar l)) a)
)
)
)
)
a
)
(defun el-3 (l y / a)
(while l
(repeat (- (caar l) y)
(setq a (cons nil a)
y (1+ y)
)
)
(setq a (cons (list (cdar l)) a)
l (cdr l)
)
(while (and l (= y (caar l)))
(setq a (cons (cons (cdar l) (car a)) (cdr a))
l (cdr l)
)
)
(setq y (1+ y))
)
(reverse a)
)
(defun el-4 (l / A LL X)
(while (setq a (car l)
l (cdr l)
)
(foreach b l
(if (< 0. (distance a b) d)
(entmakex (list '(0 . "LINE") (cons 10 a) (cons 11 b)))
)
)
)
)
(defun el-5 (l la / A)
(foreach a l
(foreach b la
(if (< 0. (distance a b) d)
(entmakex (list '(0 . "LINE") (cons 10 a) (cons 11 b)))
)
)
)
)
(defun c:el (/ D E I L S TI Y)
(setq ti (car (_VL-TIMES))
i -1
d 70.
s nil
s (ssget "_x" '((0 . "point")))
l nil
)
(while (setq e (ssname s (setq i (1+ i)))) (setq l (cons (cdr (assoc 10 (entget e))) l)))
(setq l (mapcar (function (lambda (a) (el-3 a y)))
(el-2 (vl-sort (el-1 l)
(function (lambda (a b)
(if (= (car a) (car b))
(>= (cadr a) (cadr b))
(> (car a) (car b))
)
)
)
)
)
)
)
(mapcar
(function (lambda (a) (mapcar (function el-4) a) (mapcar (function el-5) a (cdr a))))
l
)
(mapcar (function (lambda (a b)
(mapcar (function el-5) a (mapcar (function append) b (cdr b)))
(mapcar (function el-5) a b)
)
)
l
(cdr l)
)
(princ (strcat "\n " (rtos (/ (- (car (_VL-TIMES)) ti) 1000.) 2 4) "sec."))
(princ)
)
new version:
Awesome Evgeniy (as usual) it's so fast compared to the first LISP codes in this thread. :kewl:
I'll have to study your algorithm...
;;from internet
(defun random ()
(setq seed (if seed
(rem (+ (* seed 15625.7) 0.21137152) 1)
0.3171943
)
)
)
(defun random-n (n)
(* n (random))
)
(defun entmakepoint (pt layer)
(entmake (list (cons 0 "POINT")
(cons 8 layer);***
(cons 6 "BYLAYER")
(cons 10 pt) ;***
(cons 39 0.0)
(cons 50 0.0)
(cons 62 256)
(cons 210 (list 0.0 0.0 1.0))
)
)
)
;;generate random points
(defun c:test1 ()
(repeat 500000
(entmakepoint (list (random-n 100000) (random-n 60000) 0) "0")
)
)
last one ... maybe :-D
10,000 - 0.0449 s
100,000 - 0.885 s
500,000 - 5.829 s
1,000,000 - 19.376 s
I was able to cut another 25% of the time using a bit of caching. only has an effect on larger sets
@ 1,000,000
I was able to cut another 25% of the time using a bit of caching. only has an effect on larger sets
@ 1,000,000
Yep, 15.92 seconds over here. Nice.
I reorg'ed mine earlier (prior to seeing these monster numbers) and was pumped that I shaved 5 seconds so it runs about 29.5 s. While writing a message to brag, fired up your code so I could get some current numbers and BAM 19 freaking seconds, and now 15 seconds. Tremendous job...goes for everyone else too.
...
for (size_t i = 0; i < len; i++)
{
const AcGePoint3d &a = points[i];
for (size_t j = i + 1; j < len; j++)
{
const AcGePoint3d &b = points[j];
if ((b.x - a.x) > distance)
{
break;
}
else if ( distsqrd(a,b) < distancesqrd)
{
AcDbLine *pLine = new AcDbLine(a,b);
if(pRec->appendAcDbEntity(pLine) != Acad::eOk)
delete pLine;
pLine->close();
}
}
}
...
Just to throw another challenge to you keeners ... fill a volume with random 3D points rather than an area with 2D points, connecting them with 3D lines wherever the 3D distance is less than or equal to a specified value. (http://www.theswamp.org/screens/mp/wave.gif)
// acrxEntryPoint.cpp
// Copyright 2010 Paul Kohut
//
//-----------------------------------------------------------------------------
//----- acrxEntryPoint.h
//-----------------------------------------------------------------------------
#include "StdAfx.h"
#include "resource.h"
#include "KdTree.h"
#include "LapseTime.h"
#include <list>
#include <algorithm>
//-----------------------------------------------------------------------------
#define szRDS _RXST("rpk")
struct PT3D{
double x, y, z;
};
inline double Distance2DSquared(PT3D * pt1, PT3D * pt2)
{
return ((pt2->x - pt1->x) * (pt2->x - pt1->x) + (pt2->y - pt1->y) * (pt2->y - pt1->y));
}
inline double Distance3DSquared(PT3D * pt1, PT3D * pt2)
{
return ((pt2->x - pt1->x) * (pt2->x - pt1->x) + (pt2->y - pt1->y) * (pt2->y - pt1->y) + (pt2->z - pt1->z) * (pt2->z - pt1->z));
}
void ConnectDotsWithin(AcDbBlockTableRecord * pRcd, CKdTree & kdTree, double dRadius)
{
RESULT_SET * pResAllPoints = kdTree.NearestRange(0.0, 0.0, 0.0, HUGE_VAL);
while(pResAllPoints->riter) {
RESULT_SET * pRes = kdTree.NearestRange(pResAllPoints->riter->item->pos[0],
pResAllPoints->riter->item->pos[1], pResAllPoints->riter->item->pos[2], dRadius);
while(pRes->riter && pRes->size > 1) {
if(!pRes->riter->item->data && Distance3DSquared((PT3D *) pResAllPoints->riter->item->pos, (PT3D *) pRes->riter->item->pos) > 0.0) {
AcDbLine * pLine = new AcDbLine(asPnt3d(pResAllPoints->riter->item->pos), asPnt3d(pRes->riter->item->pos));
if(pRcd->appendAcDbEntity(pLine) != Acad::eOk) {
delete pLine;
} else {
pLine->close();
}
}
pRes->riter = pRes->riter->next;
}
kdTree.ResultFree(pRes);
pResAllPoints->riter->item->data = (void *) 1;
pResAllPoints->riter = pResAllPoints->riter->next;
}
kdTree.ResultFree(pResAllPoints);
}
void ConnectDotsOptimized(void)
{
double dRadius = 0.0;
if(acedGetDist(0, _T("\nMax segment length: "), &dRadius) != RTNORM)
return;
CKdTree kdTree;
kdTree.Create(3);
CLapseTime timer;
AcDbDatabase * pDb = acdbHostApplicationServices()->workingDatabase();
AcDbBlockTableRecord * pRec;
if(acdbOpenAcDbObject((AcDbObject *&) pRec, pDb->currentSpaceId(), AcDb::kForWrite) == Acad::eOk) {
AcDbBlockTableRecordIterator * pIter = NULL;
if(pRec->newIterator(pIter) == Acad::eOk) {
AcGePoint3d pt;
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) {
pt = pPoint->position();
kdTree.Insert(pt.x, pt.y, pt.z, NULL);
}
pEnt->close();
}
}
delete pIter;
}
ConnectDotsWithin(pRec, kdTree, dRadius);
pRec->close();
}
acutPrintf(_T("\n%.3f milliseconds."), timer.LapseTimeMS());
}
//////////////////////////////////////////////////////////////////////////
// Begin ARX Initialization
//////////////////////////////////////////////////////////////////////////
//-----------------------------------------------------------------------------
//----- ObjectARX EntryPoint
class CConnectDotsApp: public AcRxArxApp {
public:
CConnectDotsApp (): AcRxArxApp () {}
virtual AcRx::AppRetCode On_kInitAppMsg (void *pkt){
// TODO: Load dependencies here
// You *must* call On_kInitAppMsg here
AcRx::AppRetCode retCode =AcRxArxApp::On_kInitAppMsg (pkt) ;
// TODO: Add your initialization code here
return (retCode) ;
}
virtual AcRx::AppRetCode On_kUnloadAppMsg (void *pkt){
// TODO: Add your code here
// You *must* call On_kUnloadAppMsg here
AcRx::AppRetCode retCode =AcRxArxApp::On_kUnloadAppMsg (pkt) ;
// TODO: Unload dependencies here
return (retCode) ;
}
virtual void RegisterServerComponents (){
}
// - rpkConnectDots.ConnectDots command (do not rename)
static void rpkConnectDotsConnectDots(void)
{
ConnectDotsOptimized();
}
} ;
//-----------------------------------------------------------------------------
IMPLEMENT_ARX_ENTRYPOINT(CConnectDotsApp)
ACED_ARXCOMMAND_ENTRY_AUTO(CConnectDotsApp, rpkConnectDots, ConnectDots, ConnectDots, ACRX_CMD_TRANSPARENT, NULL)
I was able to cut another 25% of the time using a bit of caching. only has an effect on larger sets
@ 1,000,000
Single core
No temps 19.304 s
Temps - 15.685 s
Multi core use
No temps - 14.704 s (about 65% of 4 core processor utilization)
With temps as object - 15.790 s (about 65% 4 C PU)
With temps as reference type - 20.770 s (about 50% 4 C PU)
With temps as pointer - 20.074 s (about 50% 4 C PU)
#pragma omp parallel
{
#pragma omp for
for (long i = 0; i < len; i++)
{
for (long j = i + 1; j < len; j++)
{
if ((points[j].x - points[i].x) > distance)
{
break;
}
else if ( distsqrd(points[i],points[j]) < distancesqrd)
{
#pragma omp critical
{
AcDbLine *pLine = new AcDbLine(points[i],points[j]);
if(pRec->appendAcDbEntity(pLine) != Acad::eOk)
delete pLine;
pLine->close();
}
}
}
}
}
#include <omp.h>
...
static void ArxDist_doit(void)
{
double distance = 0.0;
if(acedGetDist(0,_T("\nGetDist: "),&distance) != RTNORM)
return;
const double distancesqrd = pow(distance,2) + 1.192092896e-07F;
__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)
{
AcDbBlockTableRecordIterator *pIter = NULL;
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;
}
std::sort(points.begin(),points.end(),psort);
size_t len = points.size();
#if 1
// No temp variables
#pragma omp parallel
{
#pragma omp for
for (long i = 0; i < len; i++)
{
for (long j = i + 1; j < len; j++)
{
if ((points[j].x - points[i].x) > distance)
{
break;
}
else if ( distsqrd(points[i],points[j]) < distancesqrd)
{
#pragma omp critical
{
AcDbLine *pLine = new AcDbLine(points[i],points[j]);
if(pRec->appendAcDbEntity(pLine) != Acad::eOk)
delete pLine;
pLine->close();
}
}
}
}
}
#else
// With Temp variables as object
#pragma omp parallel
{
#pragma omp for
for (long i = 0; i < len; i++)
{
const AcGePoint3d a = points[i];
for (long j = i + 1; j < len; j++)
{
const AcGePoint3d b = points[j];
if ((b.x - a.x) > distance)
{
break;
}
else if ( distsqrd(a,b) < distancesqrd)
{
#pragma omp critical
{
AcDbLine *pLine = new AcDbLine(a,b);
if(pRec->appendAcDbEntity(pLine) != Acad::eOk)
delete pLine;
pLine->close();
}
}
}
}
}
#endif
pRec->close();
}
QueryPerformanceCounter((LARGE_INTEGER*)&end);
acutPrintf(_T("\n%g Seconds\n"), (double)(end - start) / (double) freq);
}
} ;
(defun el-1 (l / A)
(while (setq a (car l)
l (cdr l)
) ;_ setq
(foreach b l
(if (< 0. (distance a b) d)
(entmakex (list '(0 . "LINE") (cons 10 a) (cons 11 b)))
) ;_ if
) ;_ while
) ;_ while
) ;_ defun
(defun el-2 (l la)
(foreach a l
(foreach b la
(if (< 0. (distance a b) d)
(entmakex (list '(0 . "LINE") (cons 10 a) (cons 11 b)))
) ;_ if
) ;_ foreach
) ;_ foreach
) ;_ defun
(defun c:el (/ A AA B C D E L S TI X Y YY)
;;(repeat 3(c:el))
(setq ti (car (_VL-TIMES))
a -1
d 70.
s (ssget "_x" '((0 . "point")))
y (mapcar (function (lambda (a) (atoi (rtos (/ a d) 2 0)))) (getvar "EXTMIN"))
x (if (minusp (car y))
(- (car y))
0
) ;_ if
y (if (minusp (cadr y))
(- (cadr y))
0
) ;_ if
yy y
) ;_ setq
(while (setq e (ssname s (setq a (1+ a))))
(setq e (cdr (assoc 10 (entget e)))
e (list (car e) (cadr e))
l (cons (cons (fix (+ x (/ (car e) d))) (cons (fix (+ y (/ (cadr e) d))) e)) l)
) ;_ setq
) ;_ while
(setq l (vl-sort l (function (lambda (a b) (<= (car a) (car b)))))
a (list (cdar l))
l (cdr l)
) ;_ setq
(while (and l (= x (caar l)))
(setq a (cons (cdar l) a)
l (cdr l)
) ;_ setq
) ;_ while
(setq a (vl-sort a (function (lambda (a b) (<= (car a) (car b)))))
x (1+ x)
) ;_ setq
(while a
(repeat (- (caar a) y)
(setq c (cons nil c)
y (1+ y)
) ;_ setq
) ;_ repeat
(while (and a (= y (caar a)))
(while (and a (= y (caar a)))
(setq b (cons (cdar a) b)
a (cdr a)
) ;_ setq
) ;_ while
(el-1 b)
(el-2 b (car c))
(setq c (cons b c)
b nil
y (1+ y)
) ;_ setq
) ;_ while
) ;_ while
(setq c (reverse c))
(while l
(setq aa c
a (list (cdar l))
y yy
l (cdr l)
b nil
c nil
) ;_ setq
(while (and l (= x (caar l)))
(setq a (cons (cdar l) a)
l (cdr l)
) ;_ setq
) ;_ while
(setq a (vl-sort a (function (lambda (a b) (<= (car a) (car b)))))
x (1+ x)
) ;_ setq
(while a
(repeat (- (caar a) y)
(setq c (cons nil c)
y (1+ y)
) ;_ setq
) ;_ repeat
(while (and a (= y (caar a)))
(while (and a (= y (caar a)))
(setq b (cons (cdar a) b)
a (cdr a)
) ;_ setq
) ;_ while
(el-1 b)
(el-2 b (car c))
(setq c (cons b c)
b nil
y (1+ y)
) ;_ setq
) ;_ while
) ;_ while
(setq c (reverse c))
(mapcar (function el-2) c (mapcar (function append) aa (cdr aa)))
(mapcar (function el-2) c aa)
) ;_ while
(princ (strcat "\n " (rtos (/ (- (car (_VL-TIMES)) ti) 1000.) 2 4)))
;;; (vl-cmdf "_erase" (ssget "_x" '((0 . "line"))) "")
(princ)
) ;_ defun
;;; the main routine
;;; 主程序
(defun C:test (/ TOL sl to ss Pts T0 pp)
;;set or input Tolerance
;; 输入容差
(if (not (setq TOL (getdist "\nPlease input the tolerance:")))
(setq TOL 70)
)
;;也可以用其他方式取得点集
;;Get the points
(setq sl '((0 . "POINT")))
(setq ss (ssget sl))
(setq Pts (getpt ss))
;;按照X值排序
;;sort by X value
(setq t0 (getvar "TDUSRTIMER"))
(setq Pts (sortx Pts))
(princ "\nThe time for sorting these points:")
(princ (* (- (getvar "TDUSRTIMER") t0) 86400))
(princ "s.")
;;函数用时估算,以了解函数性能
;;how long it will take
(setq t0 (getvar "TDUSRTIMER"))
(setq pp (f2 Pts TOL))
(princ "\nThe time for finding the pairs of points:")
(princ (* (- (getvar "TDUSRTIMER") t0) 86400))
(princ "s.")
(if pp
;;画最短距离的点对集的连线,可能有多条
;;draw the line between these pairs.
(foreach n pp
(entmake
(list
'(0 . "LINE")
(cons 10 (car n))
(cons 11 (cadr n))
(cons 62 1)
)
)
)
)
(princ)
)
;;; 取点函数,其中i为点的编号
;;; get points.
(defun getpt (ss / i l a b c)
(setq i 0)
(if ss
(repeat (sslength ss)
(setq a (ssname ss i))
(setq b (entget a))
(setq c (cdr (assoc 10 b)))
(setq l (cons c l))
(setq i (1+ i))
)
)
(reverse l)
)
;;; 从0到k的点集
;;; the elments of points from 0 to k
(defun Array (Pts k / L1 L2)
(setq L1 Pts)
(repeat k
(setq L2 (cons (car L1) L2))
(setq L1 (cdr L1))
)
(list (reverse L2) L1)
)
;;; 对X排序
;;; sort points by X value.
(defun sortX (pts)
(vl-sort pts (function (lambda (e1 e2) (< (car e1) (car e2))))
)
)
;;; 在带形区域查找
;;; search these points in the Tolerance
(defun searchX (Pts x1 x2 / pp)
(foreach n Pts
(if (and (>= (car n) x1)
(<= (car n) x2)
)
(setq pp (cons n pp))
)
)
(reverse pp)
)
;;; 用递归解决
;;; use recursion to solve this problem (divide-and-conquer)
(defun f2 (Pts TOL / l m n a b pt1 pt2 ppp Pts0 Pts1 Pts2 Pts3 Pts4
midptx mind1 mind2 qqq)
(setq l (length Pts))
(cond
((= l 1) nil)
((= l 2)
;; 两点还用说
;; judge directly if two points
(setq pt1 (car Pts))
(setq pt2 (cadr Pts))
(if (and (equal (car pt1) (car pt2) TOL)
(equal (cadr pt1) (cadr pt2) TOL)
(< (distance pt1 pt2) TOL)
)
(list Pts)
)
)
((> l 2)
(setq n (/ l 2)
m (- l n)
)
;;分治
;;divide
(setq Pts0 (Array Pts n))
(setq Pts1 (car Pts0))
(setq Pts2 (cadr Pts0))
;;递归左边
;;recurse left
(setq mind1 (f2 Pts1 TOL))
;;递归右边
;;recurse right
(setq mind2 (f2 Pts2 TOL))
;;合并左右集合
;;merge them
(setq ppp (append mind1 mind2))
;;对带形区域的合并
;;consider the points in the Tolerance
(setq midptx (caar Pts2))
(setq a (- midptx TOL)
b (+ midptx TOL)
)
(if (setq Pts3 (searchX Pts1 a midptx))
(if (setq Pts4 (searchx Pts2 midptx b))
(foreach n Pts3
(foreach e Pts4
(if (equal (cadr n) (cadr e) TOL)
(if (<= (distance n e) TOL)
(setq ppp (cons (list n e) ppp))
)
)
)
)
)
)
ppp
)
)
)
I improved my program,so now it's faster ,and I made an arx routine. it runs 10 times faster than VLX.
I has compiled the Arx code .see the attachments, maybe it's faster than pkohut's.
I learned a lot from pkohut's. KD-tree,is advanced and hard to understand.
Any ideas, thanks.
Acad::ErrorStatus CreatLine(const AcGePoint3d & ptStart,const AcGePoint3d & ptEnd
AcDbBlockTableRecord *pBlockTableRecord)
{
AcDbLine *pLine = new AcDbLine(ptStart, ptEnd);
es=pBlockTableRecord->appendAcDbEntity(pLine);
pLine->close();
// pBlockTableRecord->close(); close this later
return es;
}
Your CreatLine function is really heavy, you ought to move all the code needed to get the AcDbBlockTableRecord out to the top level function and pass a pointer.Daniel,you are right.Actually ,I realized this question when I begun to optimize my code,but in that time , I was busy, so I didn't post my modified code.
...
;;; 用递归解决
;;; use recursion to solve this problem (divide-and-conquer)
(defun f2 (Pts TOL / l m n pt1 pt2 Pts1 Pts2 Pts3 Pts4 Mx x1 x2)
(setq l (length Pts))
(cond
((= l 1) nil)
((= l 2)
;; 两点还用说
;; judge directly if 2 points
(setq pt1 (car Pts))
(setq pt2 (cadr Pts))
(if (and (equal (car pt1) (car pt2) TOL)
(equal (cadr pt1) (cadr pt2) TOL)
(<= (distance pt1 pt2) TOL)
)
(setq pairs (cons pt1 pairs) ;add the start point into pairs;
pairs (cons pt2 pairs) ;add the end point into pairs;
)
)
)
((> l 2)
(setq n (/ l 2)
m (- l n)
)
(setq Pts1 (Array Pts n)) ;Divide
(setq Pts2 (cadr Pts1)) ;R points
(setq Pts1 (car Pts1)) ;L points
(f2 (reverse Pts1) TOL) ;Recurse Left
(f2 Pts2 TOL) ;Recurse Right
(setq Mx (caar Pts2)) ;Midpoint
(setq x1 (- Mx tol)) ;L points in the TOL
(setq x2 (+ Mx tol)) ;R points in the TOL
(while (and pts2 (<= (caar pts2) x2)) ;gather the points in the Right
(setq Pts3 (cons (cons (car pts2) 1) Pts3)) ;1 means in the Right
(setq pts2 (cdr pts2))
)
(setq pts3 (reverse pts3))
(while (and pts1 (>= (caar pts1) x1)) ;gather the points in the Left
(setq Pts3 (cons (cons (car pts1) 0) Pts3)) ;0 means in the Left
(setq pts1 (cdr pts1))
)
(setq Pts3 (sorty pts3 tol)) ;sort these points by Y value
(setq pts4 (cdr pts3)) ;consider these points one by one
(if pts4
(while pts3
(setq pt1 (caar pts3))
(setq pt2 (caar pts4))
(while (and pts4 (<= (- (cadr pt2) (cadr pt1)) tol))
(if (and (/= (cdar pts3) (cdar pts4)) ;different side;
(equal (car pt1) (car pt2) tol) ;x value in tol;
(<= (distance pt1 pt2) tol) ;distance <= tol;
)
(setq pairs (cons pt1 pairs) ;add the start point into pairs;
pairs (cons pt2 pairs) ;add the end point into pairs;
)
)
(setq pts4 (cdr pts4))
(setq pt2 (caar pts4))
)
(setq pts3 (cdr pts3))
(setq pts4 (cdr pts3))
)
)
Pairs
)
)
)Ok ,this is the main code. I considered the duplicate points. So this program can be used for removing duplicate elements, or finding those closed points,etc. these kinds of finding will get easier and quicklier. I have no idea for analysing the time complexity of this algorithm. maybe it's in O(n*log(n)),because I tested it many times, found that it's faster than pkohut's. I think his algorithm is in O(n*log(n)).//-----------------------------------------------------------------------------
//----- acrxEntryPoint.h
//-----------------------------------------------------------------------------
#include "StdAfx.h"
#include "resource.h"
#include "tchar.h"
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
//-----------------------------------------------------------------------------
#define szRDS _RXST("")
//-----------------------------------------------------------------------------
//----- ObjectARX EntryPoint
class CCreatRandomPointsApp : public AcRxArxApp {
public:
CCreatRandomPointsApp () : AcRxArxApp () {}
virtual AcRx::AppRetCode On_kInitAppMsg (void *pkt) {
// TODO: Load dependencies here
// You *must* call On_kInitAppMsg here
AcRx::AppRetCode retCode =AcRxArxApp::On_kInitAppMsg (pkt) ;
// TODO: Add your initialization code here
acutPrintf(_T("\nThe command is:rand."));
return (retCode) ;
}
virtual AcRx::AppRetCode On_kUnloadAppMsg (void *pkt) {
// TODO: Add your code here
// You *must* call On_kUnloadAppMsg here
AcRx::AppRetCode retCode =AcRxArxApp::On_kUnloadAppMsg (pkt) ;
// TODO: Unload dependencies here
return (retCode) ;
}
virtual void RegisterServerComponents () {
}
// - CreatRandomPoints._Rand command (do not rename)
static void CreatRandomPoints_Rand(void)
{
// Add your code for command CreatRandomPoints._Rand here
// 输入范围
ads_point pt1,pt2;
int ret=acedGetPoint(NULL,_T("\nPlease drag the First coner of region: "),pt1);
if (ret!= RTNORM)
{
acutPrintf(_T("\nInvalid!"));
return;
}
ret=acedGetCorner(pt1,_T("\nPlease drag the region: "),pt2) ;
if (ret!= RTNORM)
{
acutPrintf(_T("\nInvalid!"));
return;
}
//转化坐标系
acdbUcs2Wcs(pt1,pt1,0);
acdbUcs2Wcs(pt2,pt2,0);
//防止点过于密集
double deltaX = (pt2[X]-pt1[X]);
double deltaY = (pt2[Y]-pt1[Y]);
if (abs(deltaX) < 1e-4 || abs(deltaY) < 1e-4 )
{
acutPrintf(_T("\nYour region is too small,Please drag it again."));
return;
}
deltaX = deltaX / 32767.0;
deltaY = deltaY / 32767.0;
//输入数量
double Number;
acedInitGet(7,NULL);
ret = acedGetReal(_T("\nPlease Enter the Number of Points:"),&Number);
if (ret!= RTNORM)
{
acutPrintf(_T("\nInvalid input!"));
return;
}
//初始化种子
srand((unsigned int)time(NULL));
//获得块表
Acad::ErrorStatus es;
AcDbBlockTable * pBlockTable;
es=acdbHostApplicationServices()->workingDatabase()->getBlockTable(pBlockTable, AcDb::kForRead);
if(es!=Acad::eOk)
return ;
pBlockTable->close();
//获得当前空间
resbuf tilemode;
AcDbBlockTableRecord * pBlockTableRecord;
acedGetVar(_T("TILEMODE"),&tilemode);
if(tilemode.resval.rint )
es=pBlockTable->getAt(ACDB_MODEL_SPACE, pBlockTableRecord,AcDb::kForWrite);
else
es=pBlockTable->getAt(ACDB_PAPER_SPACE, pBlockTableRecord,AcDb::kForWrite);
if (es!=Acad::eOk)
return;
AcDbPoint * pPoint;
AcGePoint3d pt;
double Cx,Cy;
Cx = pt1[X];
Cy = pt1[Y];
pt.z = 0;
//画点
int n = (int) Number;
for (int i = 0;i<n;i++)
{
pt.x = Cx + deltaX * rand();
pt.y = Cy + deltaY * rand();
pPoint = new AcDbPoint(pt);
es=pBlockTableRecord->appendAcDbEntity(pPoint);
pPoint->close();
}
pBlockTableRecord->close();
}
} ;
//-----------------------------------------------------------------------------
IMPLEMENT_ARX_ENTRYPOINT(CCreatRandomPointsApp)
ACED_ARXCOMMAND_ENTRY_AUTO(CCreatRandomPointsApp, CreatRandomPoints, _Rand, Rand, ACRX_CMD_TRANSPARENT, NULL)
Now I attached the new code.
BOOL PointArray::equal(const double & v1, const double & v2, const double & e)
{
return (fabs(v2 - v1) <= e);
}
strikeout this code.
BOOL PointArray::equal const point2d & pt1, const point2d & pt2, const double & e) const
{
return distance(pt1, pt2) <= e); return;
}
// line is OK, checking if x and y are in range before doing the more expensive distance calculation.
if (equal (pt1.x, pt2.x,TOL) && equal (pt1.y ,pt2.y,TOL) && distance(pt1,pt2) <= TOL)
strike this code too
if (distance(pt1,pt2) <= TOL)
#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);
}
you gave me some very good advice and a big challenge,I am happy to accept it.
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.
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.<...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)
That statement may be a little difficult for 'bing' to translate :) (sensibly, anyway)
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)
So, after implementing a tree optimization routine the final numbers are:
6.15 seconds raw speed, with the 1 million point test.
Hey Highflybird,
I've taken the liberty to use your QuickSortY ...
Paul try this arx (07) and see if there is any difference, note: its not creating the lines
the command is doit
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.
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);
}
"pts[i]" ---the "[ i ]",it 's dispeared. it can't be displayed? :-)
Quote
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);
}
Based on the code snippet you show, I don't think so. The original is better as it does a quick Manhattan Distance check before doing the more expensive Pythagorean calc.in my test, yep,the distance check is faster,but just a little,nearly the same speed.I don't know the why.
Based on the code snippet you show, I don't think so. The original is better as it does a quick Manhattan Distance check before doing the more expensive Pythagorean calc.in my test, yep,the distance check is faster,but just a little,nearly the same speed.I don't know the why.
in my old opinion,it's slower.
Floating point subtraction and addition are faster operations than multiplication and division. fabs, abs, etc take a single machine clock cycle to perform.
(Pts3[ i ].R ^ Pts3[ j ].R) is a logistic calculation, distanceto includes substracts,additions,and multiplication and squrare root calculation.Ok,this is another topic.pkohut, thanks your replay.
Floating point subtraction and addition are faster operations than multiplication and division. fabs, abs, etc take a single machine clock cycle to perform.
butCode: [Select](Pts3[ i ].R ^ Pts3[ j ].R)is a logistic calculation, distanceto includes substracts,additions,and multiplication and squrare root calculation.Ok,this is another topic.pkohut, thanks your replay.
if((pts3[i].R ^ pts3[j].R) &&
fabs(pt2.x - pt1.x) <= distance && fabs(pt2.y - pt1.y) <= distance &&
pt1.distanceTo(pt2) <= distance) {
// yep, build a line
AddLine(pRec, pt1, pt2);
}
Floating point subtraction and addition are faster operations than multiplication and division. fabs, abs, etc take a single machine clock cycle to perform.
butCode: [Select](Pts3[ i ].R ^ Pts3[ j ].R)is a logistic calculation, distanceto includes substracts,additions,and multiplication and squrare root calculation.Ok,this is another topic.pkohut, thanks your replay.
This is better, check if point 3 i and j are on the right side first, then do the calculations.Code: [Select]
if((pts3[i].R ^ pts3[j].R) &&
fabs(pt2.x - pt1.x) <= distance && fabs(pt2.y - pt1.y) <= distance &&
pt1.distanceTo(pt2) <= distance) {
// yep, build a line
AddLine(pRec, pt1, pt2);
}
if(pts3[i].R ^ pts3[j].R)
if(fabs(pt2.x - pt1.x) <= distance)
if(fabs(pt2.y - pt1.y) <= distance)
if(pt1.distanceTo(pt2) <= distance) {
// yep, build a line
AddLine(pRec, pt1, pt2);
}
The first 3 compares can all be combined into a single logical & statement which will get rid of 2 if's, and we come out ahead in clock cycles. So the above if statement can be rewritten to -if( (fabs(pt2.x - pt1.x) <= distance) & (fabs(pt2.y - pt1.y) <= distance) &
(pts3[i].R ^ pts3[j].R) &&
(pt1.distanceTo(pt2) <= distance))
struct IndexAndSide
{
int N; // index of point
bool R; // is right side
void Set(int n, bool r) { N = n; R = r; }
};
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) {
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].Set(u - 1 - i, false);
//pts3[i].N = u - 1 - i;
//pts3[i].R = false;
}
for(int i = 0; i < k2; ++i) {
pts3[i + k1].Set(u + i, true);
//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))
if( (fabs(pt2.x - pt1.x) <= distance) & (fabs(pt2.y - pt1.y) <= distance) &
(pts3[i].R ^ pts3[j].R) &&
(pt1.distanceTo(pt2) <= distance))
{
// yep, build a line
AddLine(pRec, pt1, pt2);
}
}
}
}
int __fastcall FindNearest(KDNODE<T, DIM> * pNode, RESULT_NODE<T, DIM> *pList)
{
if(!pNode)
return 0;
int added_res = 0; // track number of nodes are added to the pList
// get the Manhattan distance from pNode to the queried position. If Manhattan distance is <= query range
// the go ahead and do more expensive KDTree distance check.
[color=purple]if( fabs(pNode->pos[0] - m_dQueryPos[0]) <=m_dQueryRange && fabs(pNode->pos[1] - m_dQueryPos[1]) <= m_dQueryRange) {[/color]
// get the KDTree distance from pNode to the queried position
double dist_sq = 0.0;
for(int i = 0; i < DIM; i++) {
dist_sq += SQ(pNode->pos[i] - m_dQueryPos[i]);
}
// Only add nodes that are within the query range
if(dist_sq <= SQ(m_dQueryRange)) {
if(RlistInsert(pList, pNode, dist_sq) < 0) {
// error
return -1;
}
added_res = 1;
}
}
// Get the KDTree direction distance from the queryPosition to pNode,
// which determines which leaf node to traverse down next.
double dx = m_dQueryPos[pNode->dir] - pNode->pos[pNode->dir];
[color=red]int ret = FindNearest(dx <= 0.0 ? pNode->left : pNode->right, pList);
if(ret >= 0 && fabs(dx) < m_dQueryRange) {[/color]
added_res += ret;
[color=red]ret = FindNearest(dx <= 0.0 ? pNode->right : pNode->left, pList);[/color]
}
if(ret < 0) {
return -1;
}
added_res += ret;
return added_res;
}
[color=red]if(ret >= 0 && fabs(dx) < m_dQueryRange)[/color]
This basically saysif(ret >= 0)
(if(fabs(dx) < m_dQueryRange) {
...
}
Changing the && operator to & makes the statement a boolean check from a single if -[color=green]if(ret >= 0 & fabs(dx) < m_dQueryRange)[/color]
Run benchmarks to determine overall effect. In this case we got a positive speed boost. Notice the purple highlighted line in the original code. Making to same change to it cause a negative impact on speed. After lots of testing I'm guessing that it is a memory cache miss that caused the slow down, so the line won't be changed.int ret = FindNearest([color=red]dx <= 0.0 ? pNode->left : pNode->right[/color], pList);For those not familiar with C or Java syntax this may look strange, it's called a Ternary operator, and basically we put a big ol if/then/else statement right in the function call parameter. It reads "if dx is less than or equal to 0.0 then use pNode->left else use pNode->right as the first parameter to FindNearest".[color=navy]if/then/else type statement[/color]
Compare dx <= 0.0
Jump to Label1 if false
Use pNode->left
Jump to Label2
Label1:
Use pNode->right
Label2:
[color=navy]if/then type statement[/color]
Compate dx <= 0.0
Jump to Label3 if false
Use pNode->left;
Label3:
See the 2 jumps in the if/then/else statement, for each pass through the if statement one of them must be executed, and those add up. Depending on the condition we will use either left or right. It would be nice if we could use the simpler if/then style statement and still have the proper left or right value set. By preloading a variable we can. The idea being, preload a variable with left or right (we'll use right), then do the check and change the preloaded variable if necessary. Here it is in code -[color=green] KDNODE<T, DIM> * pNode1 = pNode->right; // preload variable
if(dx <= 0.0) {
pNode1 = pNode->left; // oops, needed left
}
int ret = FindNearest(pNode1, pList);[/color]
Test the changes. Nice boost in speed. So when you have 2 conditions and the variables are simple, then preloading may be beneficial (watch out for memory cache misses).[color=red]ret = FindNearest(dx <= 0.0 ? pNode->right : pNode->left, pList);[/color]
... becomes ...[color=purple] pNode1 = pNode->left;
if(dx <= 0.0)
pNode1 = pNode->right;
ret = FindNearest(pNode1, pList);[/color]
Test the changes, oh ya, very nice speed boost. However, -[color=green] pNode1 = pNode->right;
if(dx > 0.0)
pNode1 = pNode->left;
ret = FindNearest(pNode1, pList);[/color]
Gives a massive speed boost. And for a final speed boost to 18% remove the nCount++ from the ConnectDotsWithin function as it causes memory cache misses.[color=green] int __fastcall FindNearest(KDNODE<T, DIM> * pNode, RESULT_NODE<T, DIM> *pList)
{
if(!pNode)
return 0;
int added_res = 0; // track number of nodes are added to the pList
// get the Manhattan distance from pNode to the queried position. If Manhattan distance is <= query range
// the go ahead and do more expensive KDTree distance check.
if( fabs(pNode->pos[0] - m_dQueryPos[0]) <= m_dQueryRange &&
fabs(pNode->pos[1] - m_dQueryPos[1]) <= m_dQueryRange) {
// get the KDTree distance from pNode to the queried position
double dist_sq = 0.0;
for(int i = 0; i < DIM; i++) {
dist_sq += SQ(pNode->pos[i] - m_dQueryPos[i]);
}
// Only add nodes that are within the query range
if(dist_sq <= SQ(m_dQueryRange)) {
if(RlistInsert(pList, pNode, dist_sq) < 0) {
// error
return -1;
}
added_res = 1;
}
}
// Get the KDTree direction distance from the queryPosition to pNode,
// which determines which leaf node to traverse down next.
double dx = m_dQueryPos[pNode->dir] - pNode->pos[pNode->dir];
KDNODE<T, DIM> * pNode1 = pNode->right;
if(dx <= 0.0) {
pNode1 = pNode->left;
}
int ret = FindNearest(pNode1, pList);
if(ret >= 0 & fabs(dx) < m_dQueryRange) {
added_res += ret;
// ret = FindNearest(dx <= 0.0 ? pNode->right : pNode->left, pList);
pNode1 = pNode->right;
if(dx > 0.0)
pNode1 = pNode->left;
ret = FindNearest(pNode1, pList);
}
if(ret < 0) {
return -1;
}
added_res += ret;
return added_res;
}[/color]
int __declspec(nothrow) __fastcall FindNearest(KDNODE<T, DIM> * pNode, RESULT_NODE<T, DIM> *pList)
{
// This function should never be called if pNode is NULL,
// otherwise program will fault
int added_res = 0; // track number of nodes are added to the pList
// get the Manhattan distance from pNode to the queried position. If Manhattan distance is <= query range
// the go ahead and do more expensive KDTree distance check.
if( fabs(pNode->pos[0] - m_dQueryPos[0]) <= m_dQueryRange &&
fabs(pNode->pos[1] - m_dQueryPos[1]) <= m_dQueryRange) {
// get the KDTree distance from pNode to the queried position
double dist_sq = 0.0;
for(int i = 0; i < DIM; i++) {
dist_sq += SQ(pNode->pos[i] - m_dQueryPos[i]);
}
// Only add nodes that are within the query range
if(dist_sq <= SQ(m_dQueryRange)) {
if(RlistInsert(pList, pNode, dist_sq) < 0) {
// error
return -1;
}
added_res = 1;
}
}
// Get the KDTree direction distance from the queryPosition to pNode,
// which determines which leaf node to traverse down next.
double dx = m_dQueryPos[pNode->dir] - pNode->pos[pNode->dir];
// preload pNode1
KDNODE<T, DIM> * pNode1 = pNode->right;
if(dx <= 0.0) {
pNode1 = pNode->left;
}
// preload ret
int ret = 0;
if(pNode1)
ret = FindNearest(pNode1, pList);
if(ret < 0)
return -1;
if(fabs(dx) < m_dQueryRange) {
added_res += ret;
ret = 0; // preload ret
pNode1 = pNode->right; // preload pNode1
if(dx > 0.0)
pNode1 = pNode->left;
if(pNode1)
if( (ret = FindNearest(pNode1, pList)) < 0)
return -1;
}
added_res += ret;
return added_res;
}
Your attitude and persistance should be a shining example to us all Paul.
I admire you very much. 12 hours,assembler code, you work very hard!Your attitude and persistance should be a shining example to us all Paul.Thank you.
The back story might change your view though. I spent about 12 hours writing assembler code to replace the FindNearest function, but could not beat VC's generated code. It ain't like the old Z80 and 6502 days.
edit: what resurrected my interest is the SSE2 intrinsics and running computations in parallel. Come to find out MSVC 2008 doesn't do a very good job with these, which lead me to the back story. Read a MS blog that MSVC 2010 generates much better code from SSE2 intrinsics.
template<class T>
static int __fastcall ManhattenDistance2dArray( T * pos1, T * pos2,
const T & range)
{
return (fabs(pos1[0] - pos2[0]) <= range) &&
(fabs(pos1[1] - pos2[1]) <= range);
}
// Replace the if statement for the Manhattan distance calc in FindNearest with this -
if(ManhattenDistance2dArray<T>(pNode->pos, m_dQueryPos, m_dQueryRange)) {
// add these to the CKdTreeT class
// define some typedefs for use withing the __asm sections. MSVC complains
// about the brackets in the templates.
typedef KDNODE<T, DIM> tKdNode;
typedef CKdTreeT<T, DIM, ORDERED> tKtree;
FORCEINLINE int __fastcall ManhattenDistance2dArray_2(const double * v1) const
{
__asm {
mov eax, offset absMask
movapd xmm4, [eax] // load mask data to r0 and r1 of xmm4
mov eax, dword ptr [ecx] // ecx is pointer to "this" object
lea eax, [ecx + tKtree.m_dQueryPos] // this + CKdTreeT::m_dQueryPos
movapd xmm1, [eax] // load queryPos to r0 and r1 of xmm1
lea eax, [ecx + tKtree.m_dQueryRange] // m_dQueryRange is a single
movsd xmm2, [eax] // double, so have to first
// stuff it into r0 of xmm1.
// Need to stuff it into r1
// too, but do something else
// ...
mov eax, [edx]tKdNode.pos // edx points to pNode
movapd xmm0, [eax] // load to r0 and r1 of xmm0
shufpd xmm2, xmm2, 0 // ... we're back with xmm2. copy r0 to r1
subpd xmm0, xmm1 // subtrack xmm0{r0:r1} - xmm1{r0:r1}
// results go to xmm0
andpd xmm0, xmm4 // make 'em positive with absMask
cmpltpd xmm0, xmm2 // compare the results of the subtraction
// m_dQueryRange. r0 and r1 will be -1 for true
// and 0 for false
movmskpd eax, xmm0 // pack the compare results to an int. They
// will be in 2's complements, so later bitwise
// checks can be applied to see which was true
// or false.
// eax is the return value.
};
}
///////////////////////////////////
// Substitute this if statement for the Manhattan distance calc in FindNearest
if(ManhattenDistance2dArray_2(pNode->pos)) {
Wow... a lot of info in this thread, most of it way over my head :-o
Cool to follow the progression and helps to open my mind a bit to how to look at code differently.
Thanks
double DoHeavyCalc(const double & d1, const double & d2, const double & d3, const double & d4) { ... } Thinking being that on a 32 bit machine it is more efficient to pass a 4 byte reference and make it const to let the compiler do any optimizations it can. For giggles I changed the function signatures todouble DoHeavyCalc(double d1, double d2, double d3, double d4) { ... } which turns out to be slightly more efficient. Running with the 1 million point shaved about 0.025 seconds, but in my new quest to 4, I'll take it.
So far it's just been simple refactoring and I'm learning a great deal. :-)
So far it's just been simple refactoring and I'm learning a great deal. :-)
That's good to hear, I really admire your persistence
Functions -
int __fastcall CKdTreeT::FindNearest(KDNODE<T, DIM> * pNode, RESULT_NODE<T, DIM> *pList);
KDNODE<T, DIM> * __fastcall CKdTreeT::FindMinimum(const KDNODE<T, DIM> * pNode, int axis, int dir) const;
KDNODE<T, DIM> * __fastcall CKdTreeT::EraseNode(KDNODE<T, DIM> * pNode, const T * pos, int dir) const;
struct KDNODE
Header file
MemoryPool.h (and few functions in KdTreeT that use the memory pool)
if(!pNode->right && !pNode->left) { ... }As explained in an earlier post. This is doing 2 if compares; if pNode->right is NULL AND if pNode->left is NULL. This can be simplified (in terms of the computer) toif(!pNode->right & !pNode->left) { ... }; 1257 : if( !pNode->right & !pNode->left) {
000ba 33 c9 xor ecx, ecx
000bc 85 c0 test eax, eax
000be 0f 94 c1 sete cl
000c1 33 db xor ebx, ebx
000c3 85 d2 test edx, edx
000c5 0f 94 c3 sete bl
000c8 85 cb test ecx, ebx
000ca 74 0a je SHORT $LN5@EraseNode
if(!( !!pNode->right & !!pNode->left)) { ... }; 1257 : if(!( !!pNode->right & !!pNode->left)) {
000aa 8b 5c 24 1c mov ebx, DWORD PTR tv271[esp+16]
000ae 85 c3 test eax, ebx
000b0 75 0a jne SHORT $LN5@EraseNode
That's rockin'!Now you need to make a version for X64. doubles are stuck into the XMM registers by default, all you need to do is pack them :evil:
Now you need to make a version for X64. doubles are stuck into the XMM registers by default, all you need to do is pack them :evil:
When you compile applications as 64-bit, the calculations get more complicated. A 64-bit program uses 64-bit pointers, and its instructions are slightly larger, so the memory requirement is slightly increased. This can cause a slight drop in performance. On the other hand, having twice as many registers and having the ability to do 64-bit integer calculations in a single instruction will often more than compensate. The net result is that a 64-bit application might run slightly slower than the same application compiled as 32-bit, but it will often run slightly faster.
It looks like Divide and Conquer has a good performance with a small distance.
and it can be converted to a LISP code.but It's very difficult for KDtree algorithm.
10K points => 12,718 lines
Elapsed milliseconds:
Get points: 35
Connect points: 21
Draw lines: 56
Total: 112
100k points => 128,026 lines
Elapsed milliseconds:
Get points: 254
Connect points: 272
Draw lines: 633
Total: 1159
1M points => 1,280,741 lines
Elapsed milliseconds:
Get points: 2453
Connect points: 2792
Draw lines: 11748
Total: 16993