Author Topic: {Challenge} connect each two points whose distance is shorter than given value  (Read 58147 times)

0 Members and 4 Guests are viewing this topic.

qjchen

  • Bull Frog
  • Posts: 285
  • Best wishes to all
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)
)

Evgeniy, you are always so cool.   "(<= (- (caar ll) x) d)" is very effective

In a 100000 points dwg generate by this code

Code: [Select]
;;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")
  )
)

your code is about 22s, My revision is about 28s, but you are in a so short code. I will try my best to make my code more quick. :)

Thank you very much.

Daniel, thanks for your c# code, it is rather quick. ~~, Thank you.

Pkohut, thanks for teaching me about the kd-tree algorithm, I will study it carefully and I hope I can also finish such code in Lisp.
« Last Edit: April 03, 2010, 11:14:32 AM by yuanqiu »
http://qjchen.mjtd.com
My blog http://chenqj.blogspot.com (Chinese, can be translate into English)

gile

  • Water Moccasin
  • Posts: 2233
  • Marseille, France
A new C# variant (using .NET Framework 2.0 and inspired by Evgeniy LISP code)

Code: [Select]
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);
        }
    }
}
« Last Edit: April 03, 2010, 12:14:47 PM by gile »
Speaking English as a French Frog

pkohut

  • Guest
Wow, wake up and what do I see but blazing speed. Nice coding Gile, Daniel, and master Evgeniy.  Now I gotta put timers in mine, test your codes, and see how your guys code scale to 100,000 and 500,000 points. Results should be interesting. But first, coffee and lots of it.

edit: btw, can someone post a .Net project of theirs. txs
« Last Edit: April 03, 2010, 01:41:09 PM by pkohut »

pkohut

  • Guest
Added time code to mine (project and ARX attached), and ran Giles for comparison.
Did tests for 10,000, 100,000, 500,000, and 1,000,000 points

Test Results
Quote
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

It's Alive!

  • BricsCAD
  • Needs a day job
  • Posts: 6956
  • AKA Daniel
Paul, Try this ARX , the command is doit... using Evgeniy's optimization
 

It's Alive!

  • BricsCAD
  • Needs a day job
  • Posts: 6956
  • AKA Daniel
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

pkohut

  • Guest
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


pkohut

  • Guest
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 ?

It's Alive!

  • BricsCAD
  • Needs a day job
  • Posts: 6956
  • AKA Daniel
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 ?

no

It's Alive!

  • BricsCAD
  • Needs a day job
  • Posts: 6956
  • AKA Daniel
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




Code: [Select]
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);
  }


It's Alive!

  • BricsCAD
  • Needs a day job
  • Posts: 6956
  • AKA Daniel

qjchen

  • Bull Frog
  • Posts: 285
  • Best wishes to all
Thanks for Daniel and pkohut, your c code all very fast.

But I am not familiar with c, so I only have to do this in Lisp.

Now I use xy grid to sort the point. It seems that now it is quick a little than Evgeiny's. But my code is toooooooo long. Sorry, I need to improve my coding skill. But at least, I think this algorithm is effective.  The next step, I want to study pkohut's kd-tree algorithm,  it seems so amazing.

Code: [Select]
(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)
)





Now In my PC, it is about 0.79s for 10,000 points, and 160s for 500,000 points.
Evgeniy's code is about 1.6s for  10,000 points, and 288s for 500,000 points.

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.



« Last Edit: April 04, 2010, 01:13:28 AM by yuanqiu »
http://qjchen.mjtd.com
My blog http://chenqj.blogspot.com (Chinese, can be translate into English)

pkohut

  • Guest
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

Nice code! For balls out speed I think this is the winner, and for a smart brute force routine this kicks azzstroids. You could squeeze even more by enabling OpenMP and surround the work loop with the proper #pragma omp statement, there for getting the most out of available processors.  Noticed that somewhere along the line the times converged between this one and mine, probably in the 800,000 range.

It's Alive!

  • BricsCAD
  • Needs a day job
  • Posts: 6956
  • AKA Daniel
wait not done yet...
ARX
1,000,000
14.569079 Seconds

It's Alive!

  • BricsCAD
  • Needs a day job
  • Posts: 6956
  • AKA Daniel
last one ... maybe  :-D

Code: [Select]
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);
  }