Create border from a set of points

[hak_vz]

Here in attachment is a drawing that contains set of terrain points that can be used to develop and test

algorithm I'm looking for. I need algorithm that would create boundary(lwpolyline) that contains edge points.

I would say it's a point set hull algorithm on steroids. I hope attached pictures illustrate what I'm looking for.

Thanks in advance for your code or idea how to tackle this problem.

[marko_ribar]

It was discussed here :

https://www.cadtutor.net/forum/topic/65519-how-to-draw-the-boundaryor-maximum-outline-of-a-point-set... 

I believe you are stucked here... We'll see if there are some more thoughts...

[hak_vz]

Thank you @marko_ribar .

After posting my request, I have spent some time Googling for this particular problem.

This task is part of triangulation problems and is called alpha shape, or concave hull, as is it mentioned in discussion at Cadtutor.

In my works I'm using my autolisp code to connect to Triangle, an excellent and tremendously fast triangulation program written by prof. J. R. Shewchuk at CMU. Program is used for creation of Delaunay triangulations and Voronoi diagrams. Unfortunately it do not have option to create alpha shape, where sets of triangulated (or Voronoi) points are prerequisite for further calculations and sorting. There a some algorithms and academic papers available (Python, CGAL, c++) so I'll try to make autolisp code from it.

[marko_ribar]

Look, ymg wrote triangulation routine and as a side sub he made XSHAPE sub function to be used with triangulation... But even then when you have triangulation, how would you determine what trinagles should be removed and your desired shape become resulting polyline... He opted for removing convex shape in steps and even then result may be unexpected (not like human would do in your posted pictures example)...

For more info search : Triangulation (revisited) topic at theswamp.org...

[hak_vz]

I've seen this post from ymg and have followed the discussion.

I'll have to dig deeper into this subject to see how it's done with available algorithms. Posted samples are looking promising. I have whole triangulation and contours code written (including Catmull-Rom smoothing we've worked on last year) . Now I want to add  option to remove outer triangles instantly. 

[marko_ribar]

This is the closest I could think of without triangulation... It worked on one my example (it has to be tested more deeply...)

(defun c:MR-ConcaveHull ( / *error* car-sort trimlst ucsf ss i p pl loop pp pln lil pea )

  (defun *error* ( m )
    (if pea
      (setvar 'peditaccept pea)
    )
    (if ucsf
      (progn
        (command-s "_.ucs" "_p")
        (command-s "_.zoom" "_p")
      )
    )
    (command-s "_.undo" "_e")
    (if m
      (prompt m)
    )
    (princ)
  )

  (defun car-sort ( lst cmp / rtn )
    (setq rtn (car lst))
    (foreach itm (cdr lst)
      (if (apply cmp (list itm rtn))
        (setq rtn itm)
      )
    )
    rtn
  )

  (defun trimlst ( l n / ll )
    (repeat n
      (setq ll (cons (car l) ll))
      (setq l (cdr l))
    )
    (vl-remove nil (reverse ll))
  )

  (vl-cmdf "_.undo" "_be")
  (if (= 0 (getvar 'worlducs))
    (progn
      (vl-cmdf "_.ucs" "_w")
      (vl-cmdf "_.plan" "")
      (setq ucsf t)
    )
  )
  (prompt "\nSelect points, blocks or circles...")
  (setq ss (ssget '((0 . "POINT,INSERT,CIRCLE"))))
  (if ss
    (progn
      (repeat (setq i (sslength ss))
        (setq p (cdr (assoc 10 (entget (ssname ss (setq i (1- i)))))))
        (setq pl (cons p pl))
      )
      (setq p (car-sort pl '(lambda ( a b ) (if (equal (car a) (car b) 1e-6) (< (cadr a) (cadr b)) (< (car a) (car b))))))
      (setq loop t)
      (while loop
        (if (null pp)
          (setq pp (car-sort (trimlst (vl-sort (vl-remove p pl) '(lambda ( a b ) (< (distance p a) (distance p b)))) 3) '(lambda ( a b ) (< (rem (+ (angle p a) (* 0.5 pi)) (* 2 pi)) (rem (+ (angle p b) (* 0.5 pi)) (* 2 pi))))))
          (progn
            (vl-cmdf "_.ucs" "_3p" "_non" (caar lil) "_non" (cadar lil) "")
            (vl-cmdf "_.ucs" "_z" -90)
            (vl-cmdf "_.ucs" "_m" "_non" (trans (cadar lil) 0 1))
            (setq pln (mapcar '(lambda ( x ) (trans x 0 1)) pl))
            (setq pp (car-sort (trimlst (vl-sort pln '(lambda ( a b ) (< (distance '(0 0 0) a) (distance '(0 0 0) b)))) 3) '(lambda ( a b ) (< (angle '(0 0) a) (angle '(0 0) b)))))
            (setq pp (trans pp 1 0))
            (vl-cmdf "_.ucs" "_p")
            (vl-cmdf "_.ucs" "_p")
            (vl-cmdf "_.ucs" "_p")
          )
        )
        (setq lil (cons (list p pp) lil))
        (setq p pp)
        (setq pl (vl-remove-if '(lambda ( x ) (equal x p 1e-6)) pl))
        (if (equal p (car (last lil)) 1e-6)
          (setq loop nil)
        )
      )
      (setq ss (ssadd))
      (foreach li lil
        (ssadd (entmakex (list '(0 . "LINE") (cons 10 (car li)) (cons 11 (cadr li)))) ss)
      )
      (setq pea (getvar 'peditaccept))
      (setvar 'peditaccept 1)
      (vl-cmdf "_.pedit" "_m" ss "" "_j")
      (while (< 0 (getvar 'cmdactive))
        (vl-cmdf "")
      )
    )
  )
  (*error* nil)
)

HTH.,

Regards, M.R.

[hak_vz]

@marko_ribar  I like this.

It definitively need some more refinements and testing, but it's a good starting point. 

What comes on my mind, after some test is:

- check that segments of boundary not self-intersect

- check for points left outside the boundary

[hak_vz]

I have found this paper, and it seams to me relative easy to implement. It uses combination of nearest neighbor and "Jarvis march" hull algorithm. 

[leeminardi]

@hak_vz The algorithm described in the concave hull paper you reference looks very promising.  I wish the author had included the results using the algorithm for the data shown in Figure 4 for different starting values of k. 

[marko_ribar · Accepted solution]

---

@hak_vz wrote:

@marko_ribar  I like this.

It definitively need some more refinements and testing, but it's a good starting point. 

What comes on my mind, after some test is:

- check that segments of boundary not self-intersect

- check for points left outside the boundary

---

I think I improved it based on your comments...

(defun c:MR-ConcaveHull ( / *error* car-sort trimlst ucsf ss i p pl loop pp pln lil pea p1 p2 li )

  (defun *error* ( m )
    (if pea
      (setvar 'peditaccept pea)
    )
    (if ucsf
      (progn
        (command-s "_.ucs" "_p")
        (command-s "_.zoom" "_p")
      )
    )
    (command-s "_.undo" "_e")
    (if m
      (prompt m)
    )
    (princ)
  )

  (defun car-sort ( lst cmp / rtn )
    (setq rtn (car lst))
    (foreach itm (cdr lst)
      (if (apply cmp (list itm rtn))
        (setq rtn itm)
      )
    )
    rtn
  )

  (defun trimlst ( l n / ll )
    (repeat n
      (setq ll (cons (car l) ll))
      (setq l (cdr l))
    )
    (vl-remove nil (reverse ll))
  )

  (vl-cmdf "_.undo" "_g")
  (if (= 0 (getvar 'worlducs))
    (progn
      (vl-cmdf "_.ucs" "_w")
      (vl-cmdf "_.plan" "")
      (setq ucsf t)
    )
  )
  (prompt "\nSelect points, blocks or circles...")
  (setq ss (ssget '((0 . "POINT,INSERT,CIRCLE"))))
  (if ss
    (progn
      (repeat (setq i (sslength ss))
        (setq p (cdr (assoc 10 (entget (ssname ss (setq i (1- i)))))))
        (setq pl (cons p pl))
      )
      (setq p (car-sort pl '(lambda ( a b ) (if (equal (car a) (car b) 1e-6) (< (cadr a) (cadr b)) (< (car a) (car b))))))
      (setq loop t)
      (while loop
        (if (null pp)
          (setq pp (car-sort (trimlst (vl-sort (vl-remove p pl) '(lambda ( a b ) (< (distance p a) (distance p b)))) 10) '(lambda ( a b ) (< (rem (+ (angle p a) (* 0.5 pi)) (* 2 pi)) (rem (+ (angle p b) (* 0.5 pi)) (* 2 pi))))))
          (progn
            (vl-cmdf "_.ucs" "_3p" "_non" (caar lil) "_non" (cadar lil) "")
            (vl-cmdf "_.ucs" "_z" -90)
            (vl-cmdf "_.ucs" "_m" "_non" (trans (cadar lil) 0 1))
            (setq pln (mapcar '(lambda ( x ) (trans x 0 1)) (vl-remove-if '(lambda ( y ) (vl-some '(lambda ( z / ip ) (and (setq ip (inters p y (car z) (cadr z))) (not (equal ip p 1e-6)) (not (equal ip y 1e-6)) (not (equal ip (car z) 1e-6)) (not (equal ip (cadr z) 1e-6)))) (vl-remove-if '(lambda ( q ) (or (equal p (car q) 1e-6) (equal p (cadr q) 1e-6))) lil))) pl)))
            (setq pln (trimlst (vl-sort pln '(lambda ( a b ) (< (distance '(0 0 0) a) (distance '(0 0 0) b)))) 10))
            (if (and (> (length lil) 3) (vl-some '(lambda ( x ) (equal (car (last lil)) (trans x 1 0) 1e-6)) pln))
              (setq pp (trans (car (last lil)) 0 1))
              (setq pp (car-sort pln '(lambda ( a b ) (< (angle '(0 0) a) (angle '(0 0) b)))))
            )
            (setq pp (trans pp 1 0))
            (vl-cmdf "_.ucs" "_p")
            (vl-cmdf "_.ucs" "_p")
            (vl-cmdf "_.ucs" "_p")
          )
        )
        (setq lil (cons (list p pp) lil))
        (setq p pp)
        (setq pl (vl-remove-if '(lambda ( x ) (equal x p 1e-6)) pl))
        (if (equal p (car (last lil)) 1e-6)
          (setq loop nil)
        )
      )
      (setq pln nil)
      (foreach p pl
        (if (vl-every '(lambda ( x ) (and (not (inters p (polar p 0.0 1e+6) (car x) (cadr x))) (not (inters p (polar p (* 0.5 pi) 1e+6) (car x) (cadr x))))) lil)
          (setq pln (cons p pln))
        )
      )
      (foreach p pln
        (setq p1 (car-sort (apply 'append lil) '(lambda ( a b ) (< (distance p a) (distance p b)))))
        (setq p2 (car-sort (vl-remove-if '(lambda ( x ) (equal x p1 1e-6)) (apply 'append lil)) '(lambda ( a b ) (< (distance p a) (distance p b)))))
        (setq li (vl-some '(lambda ( x ) (if (or (and (equal p1 (car x) 1e-6) (equal p2 (cadr x) 1e-6)) (and (equal p1 (cadr x) 1e-6) (equal p2 (car x) 1e-6))) x)) lil))
        (setq lil (vl-remove li lil))
        (setq lil (cons (list p1 p) lil))
        (setq lil (cons (list p2 p) lil))
      )
      (setq ss (ssadd))
      (foreach li lil
        (ssadd (entmakex (list '(0 . "LINE") (cons 10 (car li)) (cons 11 (cadr li)))) ss)
      )
      (setq pea (getvar 'peditaccept))
      (setvar 'peditaccept 1)
      (vl-cmdf "_.pedit" "_m" ss "" "_j")
      (while (< 0 (getvar 'cmdactive))
        (vl-cmdf "")
      )
    )
  )
  (*error* nil)
)

Regards, M.R.

[marko_ribar]

Just to inform...

Because I had some lacks and had a time (less than 30 min. limit) I edited code... But still, on one my example it breaked with error 2d/3d point : nil...

All in all, this routine is very sensitive and sometimes it works well, and sometimes not...

Now I think that I won't edit it further as I am satisfied and with present coding and probably won't have much nerves to search for errors and debugg it to be bulletproof...

[hak_vz]

@marko_ribar  Excellent work. I really like it.

I've added some changes to your code to enable work with my original set of 3d points.

(defun take (amount lst / ret)(repeat amount (setq ret (cons (car lst) (take (1- amount) (cdr lst))))))

and in  code

(setq p (append (take 2 (cdr (assoc 10 (entget (ssname ss (setq i (1- i))))))) (list 0.0)))

This moves point coordinates to z = 0.

 and

(setvar 'cmdecho 0)

Here is one of my results, it's almost perfect. I'll do some more test.

One suggestion.  For large sets with thousands of points I would let user to create internal polyline boundary to exclude in selection all points inside it, to make code run faster.

I'll select your code as a solution, although there are some improvements needed,  I'm sure you'll add later.

Thank you for your time and effort.

[marko_ribar · Accepted solution]

I improved it further more... Only lack I now see is that it breaks with error if selected more than one point cloud... I've implemented your remarks and added few more things of my own... I think that this is it... Had to add alert instead of prompt when starting routine - info is important to understand and I just decided to put it in alert...

Regards, M.R.

BTW. I noticed that AutoCAD crashed when typing : conc at Command: prompt... This is why I choosed MR-ConcaveHull (just with ConcaveHull - CAD fatal errors after typed Conc...)

(defun c:MR-ConcaveHull ( / *error* unique uniquecoord car-sort trimlst numrayintlil cmde ucsf ss ti i p pl loop pp pln f lil pea p1 p2 li )

  (defun *error* ( m )
    (if ucsf
      (progn
        (command-s "_.ucs" "_p")
        (command-s "_.zoom" "_p")
      )
    )
    (if (= 8 (logand 8 (getvar (quote undoctl))))
      (command-s "_.undo" "_e")
    )
    (if pea
      (setvar (quote peditaccept) pea)
    )
    (if cmde
      (setvar (quote cmdecho) cmde)
    )
    (if m
      (prompt m)
    )
    (princ)
  )

  (defun unique ( lst )
    (if lst
      (cons
        (car lst)
        (unique
          (vl-remove-if
            (function (lambda ( x )
              (equal x (car lst) 1e-6)
            ))
            (cdr lst)
          )
        )
      )
    )
  )

  (defun uniquecoord ( lst a )
    (if lst
      (cons
        (car lst)
        (uniquecoord
          (vl-remove-if
            (function (lambda ( x )
              (cond
                ( (= a 0.0) (equal (cadr x) (cadar lst) 1e-6) )
                ( (= a (* 0.5 pi)) (equal (car x) (caar lst) 1e-6) )
              )
            ))
            (cdr lst)
          )
          a
        )
      )
    )
  )

  (defun car-sort ( lst cmp / rtn )
    (setq rtn (car lst))
    (foreach itm (cdr lst)
      (if (apply cmp (list itm rtn))
        (setq rtn itm)
      )
    )
    rtn
  )

  (defun trimlst ( lst n / ll )
    (repeat n
      (setq ll (cons (car lst) ll))
      (setq lst (cdr lst))
    )
    (vl-remove nil (reverse ll))
  )
;|
  (defun numrayintlil ( p lil a / ip lst )
    (foreach li lil
      (if (setq ip (inters p (polar p a 1e+6) (car li) (cadr li)))
        (setq lst (cons ip lst))
      )
    )
    (length (uniquecoord (unique lst) a))
  )
|;
  (defun numrayintlil ( p lil a / ip l1 l2 l3 )
    (foreach li lil
      (if (setq ip (inters p (polar p a 1e+6) (car li) (cadr li)))
        (setq l1 (cons ip l1))
      )
    )
    (setq l1 (unique l1))
    (foreach li lil
      (if (setq ip (inters (polar p (+ a (* 0.5 pi)) 1e-2) (polar (polar p (+ a (* 0.5 pi)) 1e-2) a 1e+6) (car li) (cadr li)))
        (setq l2 (cons ip l2))
      )
    )
    (setq l2 (unique l2))
    (foreach li lil
      (if (setq ip (inters (polar p (- a (* 0.5 pi)) 1e-2) (polar (polar p (- a (* 0.5 pi)) 1e-2) a 1e+6) (car li) (cadr li)))
        (setq l3 (cons ip l3))
      )
    )
    (setq l3 (unique l3))
    (cond
      ( (and
          (= (length l1) (length l2) (length l3))
          (vl-every
            (function (lambda ( x )
              (vl-some
                (function (lambda ( y )
                  (equal x y 0.1)
                ))
                (append l2 l3)
              )
            ))
            l1
          )
        )
        (length l1)
      )
      ( (and
          (> (length l1) (length l2))
          (> (length l1) (length l3))
        )
        (length l2)
      )
      ( t
        (1+ (length (vl-remove-if (function (lambda ( x ) (vl-some (function (lambda ( y ) (equal x y 0.1))) (append l2 l3)))) l1)))
      )
    )
  )

  (setq cmde (getvar (quote cmdecho)))
  (setvar (quote cmdecho) 0)
  (if (= 8 (logand 8 (getvar (quote undoctl))))
    (vl-cmdf "_.undo" "_e")
  )
  (vl-cmdf "_.undo" "_m")
  (if (= 0 (getvar (quote worlducs)))
    (progn
      (vl-cmdf "_.ucs" "_w")
      (vl-cmdf "_.plan" "")
      (setq ucsf t)
    )
  )
  (alert "Select points, blocks or circles... WARNING : SELECTED ENTITIES MUST BELONG ONE POINT CLOUD GROUP - ISLANDS ARE NOT ALLOWED...")
  (if (setq ss (ssget (list (cons 0 "POINT,INSERT,CIRCLE"))))
    (progn
      (setq ti (car (_vl-times)))
      (repeat (setq i (sslength ss))
        (setq p
          (append
            (mapcar (function +)
              (list 0.0 0.0)
              (cdr (assoc 10 (entget (ssname ss (setq i (1- i))))))
            )
            (list 0.0)
          )
        )
        (setq pl (cons p pl))
      )
      (setq p
        (car-sort pl
          (function (lambda ( a b )
            (if (equal (car a) (car b) 1e-6)
              (< (cadr a) (cadr b))
              (< (car a) (car b))
            )
          ))
        )
      )
      (setq loop t)
      (while loop
        (if (null pp)
          (setq pp
            (car-sort
              (trimlst
                (vl-sort (vl-remove p pl)
                  (function (lambda ( a b )
                    (< (distance p a) (distance p b))
                  ))
                )
                10
              )
              (function (lambda ( a b )
                (if
                  (equal
                    (rem (+ (angle p a) (* 0.5 pi)) (* 2 pi))
                    (rem (+ (angle p b) (* 0.5 pi)) (* 2 pi))
                    1e-6
                  )
                  (< (distance p a) (distance p b))
                  (< (rem (+ (angle p a) (* 0.5 pi)) (* 2 pi)) (rem (+ (angle p b) (* 0.5 pi)) (* 2 pi)))
                )
              ))
            )
          )
          (progn
            (vl-cmdf "_.ucs" "_3p" "_non" (caar lil) "_non" (cadar lil) "")
            (vl-cmdf "_.ucs" "_z" -90.0)
            (vl-cmdf "_.ucs" "_m" "_non" (trans (cadar lil) 0 1))
            (setq pln
              (mapcar
                (function (lambda ( x )
                  (trans x 0 1)
                ))
                (vl-remove-if
                  (function (lambda ( y )
                    (vl-some
                      (function (lambda ( z / ip )
                        (and
                          (setq ip (inters p y (car z) (cadr z)))
                          (not (equal ip p 1e-6))
                          (not (equal ip y 1e-6))
                          (not (equal ip (car z) 1e-6))
                          (not (equal ip (cadr z) 1e-6))
                        )
                      ))
                      (vl-remove-if
                        (function (lambda ( q )
                          (or
                            (equal p (car q) 1e-6)
                            (equal p (cadr q) 1e-6)
                          )
                        ))
                        lil
                      )
                    )
                  ))
                  pl
                )
              )
            )
            (setq pln
              (trimlst
                (vl-sort pln
                  (function (lambda ( a b )
                    (< (distance (list 0.0 0.0 0.0) a) (distance (list 0.0 0.0 0.0) b))
                  ))
                )
                10
              )
            )
            (if
              (and
                (> (length lil) 3)
                (vl-some
                  (function (lambda ( x )
                    (equal (car (last lil)) (trans x 1 0) 1e-6)
                  ))
                  pln
                )
              )
              (progn
                (setq pp
                  (car-sort pln
                    (function (lambda ( a b )
                      (if (equal (angle (list 0.0 0.0) a) (angle (list 0.0 0.0) b) 1e-6)
                        (< (distance (list 0.0 0.0 0.0) a) (distance (list 0.0 0.0 0.0) b))
                        (< (angle (list 0.0 0.0) a) (angle (list 0.0 0.0) b))
                      )
                    ))
                  )
                )
                (setq f t)
              )
              (if (null f)
                (setq pp
                  (car-sort pln
                    (function (lambda ( a b )
                      (if (equal (angle (list 0.0 0.0) a) (angle (list 0.0 0.0) b) 1e-6)
                        (< (distance (list 0.0 0.0 0.0) a) (distance (list 0.0 0.0 0.0) b))
                        (< (angle (list 0.0 0.0) a) (angle (list 0.0 0.0) b))
                      )
                    ))
                  )
                )
                (setq pp (trans (car (last lil)) 0 1))
              )
            )
            (setq pp (trans pp 1 0))
            (vl-cmdf "_.ucs" "_p")
            (vl-cmdf "_.ucs" "_p")
            (vl-cmdf "_.ucs" "_p")
          )
        )
        (setq lil (cons (list p pp) lil))
        (setq p pp)
        (setq pl
          (vl-remove-if
            (function (lambda ( x )
              (equal x p 1e-6)
            ))
            pl
          )
        )
        (if (equal p (car (last lil)) 1e-6)
          (setq loop nil)
        )
      )
      (setq pln nil)
      (foreach p pl
        (if
          (not
            (or
              (= (rem (numrayintlil p lil 0.0) 2) 1)
              (= (rem (numrayintlil p lil (* 0.5 pi)) 2) 1)
            )
          )
          (setq pln (cons p pln))
        )
      )
      (foreach p pln
        (setq p1
          (car-sort
            (apply (function append) lil)
            (function (lambda ( a b )
              (< (distance p a) (distance p b))
            ))
          )
        )
        (setq p2
          (car-sort
            (vl-remove-if
              (function (lambda ( x )
                (equal x p1 1e-6)
              ))
              (apply (function append) lil)
            )
            (function (lambda ( a b )
              (< (distance p a) (distance p b))
            ))
          )
        )
        (setq li
          (vl-some
            (function (lambda ( x )
              (if
                (or
                  (and
                    (equal p1 (car x) 1e-6)
                    (equal p2 (cadr x) 1e-6)
                  )
                  (and
                    (equal p1 (cadr x) 1e-6)
                    (equal p2 (car x) 1e-6)
                  )
                )
                x
              )
            ))
            lil
          )
        )
        (setq lil (vl-remove li lil))
        (setq lil (cons (list p1 p) lil))
        (setq lil (cons (list p2 p) lil))
      )
      (setq ss (ssadd))
      (foreach li lil
        (ssadd
          (entmakex
            (list
              (cons 0 "LINE")
              (cons 10 (car li))
              (cons 11 (cadr li))
            )
          )
          ss
        )
      )
      (setq pea (getvar (quote peditaccept)))
      (setvar (quote peditaccept) 1)
      (vl-cmdf "_.pedit" "_m" ss "" "_j")
      (while (< 0 (getvar (quote cmdactive)))
        (vl-cmdf "")
      )
      (prompt "\nElapsed time : ") (prompt (rtos (- (car (_vl-times)) ti) 2 16)) (prompt " milliseconds...")
    )
  )
  (prompt "\nFor UNDO - type \"UNDO\" - \"Back\" option...")
  (*error* nil)
)

[Sea-Haven]

Two comments about deleting triangles.

Having dealt with TINS for like 40 years removing outside tringles can be done a few different ways, the 1st is a simple "erase long triangles" this removes a massive amount when concave shapes involved. the second is to make a boundary like your trying to do, just offset a tiny bit out and delete all crossing tins, something like Cookiecutter2.lsp method. There is also drag a pline and erase triangles, this is useful for buildings so no contours show through.

Lastly you are re-inventing the wheel the answers are available in commercial software and not just CIV3D, the delete methods I spoke about were available in like 1980's, as suggested it may be a worthwhile function in YMG code to add long triangle check as a 1st step. For me I use Civil Site Design or Stringer reasonably priced software.

[pbejse]

---

@hak_vz wrote:

I have found this paper, and it seams to me relative easy to implement. It uses combination of nearest neighbor and "Jarvis march" hull algorithm. 

---

Interesting read "Algorithm 1: The Concave Hull algorithm. "

Thank you for the link @hak_vz 

[Yasir.Aman]

Hi marko_ribar,
Sir, you saved me a ton of labor with this code of yours. Allow me to say, in fact, you saved my a** today. I don't have enough words to thank you. Just wanted to appreciate your effort and work you put into such pieces of code.
Stay blessed and live long.
~Yasir
(From Pakistan, working on a project in IRAQ)

[ajay_k5NRG4]

IT TAKES LOT OF TIME 

[autoid374ceb4990]

Sea-Haven:

"Lastly you are re-inventing the wheel the answers are available in commercial software and not just CIV3D, the delete methods I spoke about were available in like 1980's, as suggested it may be a worthwhile function in YMG code to add long triangle check as a 1st step. For me I use Civil Site Design or Stringer reasonably priced software."

But that takes all the fun out of writing your own code.  I wrote code to do delaunay triangulations to produce TIN's back in the 80's and my solution to get rid of those long triangles inside the concave areas in the TIN was just to write  a routine to erase the long triangles with sides longer than a user defined distance.  I also had a method of erasing all the triangles outside a closed polyline, but I had to draw the polyline manually.  The code presented in this post would make defining a boundary polyline simple.