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

[It's Alive!]

gile, your lambda expression to sort the points is awesome!

[pkohut]

last one ... maybe 

There ain't much skin left on that cat. 

Nice clean code and then -

Code: [Select]

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.

[ElpanovEvgeniy]

Sure master Evgeniy can easily improve the code and quick quick than me, but  I am just try to study this funny algorithm. 

The algorithm is as following picture.

An interesting algorithm, similar to I use in many of its programs.
Hopefully, tonight will be free time and I will write the new version of the code.

[gile]

gile, your lambda expression to sort the points is awesome!

Thanks.
Maybe a 'LISP like' way of thinking/writing

Code: [Select]

(vl-sort pts '(lambda (p1 p2) (< 'car p1) (car p2))))

[ElpanovEvgeniy]

new version:

Code: [Select]

(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)
)

[ElpanovEvgeniy]

new version:

recommend testing after compiling:

[gile]

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...

[ElpanovEvgeniy]

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...

I used an algorithm similar to the picture of the messages Chen

[qjchen]

Super cool, Evgeniy, very very quick.

Now about a 500,000 points dwg generated by the following code.

Your code is about only 48s (vlx) or 58s (lsp), my code is about 160s.

sure, Daniel and pkohut and Gile 's  code is much quick. But maybe C can 10 times quick than Lisp?

I need to deep study for your code structure, Thank you

---------------------------------

Just find a bugs, when in a only two points dwg, p1=(60,0,0) p2=(180,0,0), the code seems to get error, it seems that the error happen in the sortting data into xy grid.

----------------------------------

Recently, I buy the book you recommend to me :"Structure and interpretation of computer programs" by Harold Abeleson in MIT.

It is a rather good book for LISP, I am learning now about the recursion and so on.

I read your article about "recursion" in dwg.ru and learn a lot, could you post it here also for theswamp friends.   

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 500000
    (entmakepoint (list (random-n 100000) (random-n 60000) 0) "0")
  )
)


[pkohut]

last one ... maybe  

Times a staggering -

10,000     -   0.0449 s
100,000    -  0.885 s
500,000    -  5.829 s
1,000,000 - 19.376 s

edit: what is the project optimizations? I say the numbers last night, but being so tired they didn't register.

[It's Alive!]

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

[pkohut]

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.

[MP]

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.

[It's Alive!]

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.

the only thing I did was add temps inside the loop.. I thought the compiler would have done this but..

Code: [Select]

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


[pkohut]

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.

This is my reorg'ed version (2).  To meet your requirements (I think) the only change needed just now was to add and use a Distance3DSquared function, everything else is already in place.

Code notes: This version follows Daniel's version and layout to facilitate post mortem comparisions.

Just replace the contents acrxEntryPoint.cpp in the VC++ project at
http://www.theswamp.org/index.php?topic=32874.msg383657#msg383657
and recompile.

Code: [Select]

// 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)