Topic: Triangulation (re-visited)

[ymg]

Rick,

I am afraid that you need the following:

Code - Auto/Visual Lisp: [Select]

1. (setq lt        (cons 6 "Depression")
2.                sc        (cons 48 0.5) ;only this line needed to be changed
3.                majdepcol (cons 62 majdepcolor)
4.                mindepcol (cons 62 mindepcolor)
5.                lgclosed  (length ccont)
6.          )
7.  

For the godep you either remove the (if  (= godep "1") and the progn
or do as you did setting it to "1".

Variable godep will be used to make the depression processing optionnal.

ymg

[d.valkanov]

Here is an implementation of Sloan algorithm. I am still working on it, but but for nor now it works quite fast. It work in 2016. Boundaries, holes and breaklines are implemented. Use SMOOTH command to smooth contours. Point on edge of triangulation is also implemented. There is no command for labeling of contours for now, but it will be ready in couple of days. Profiles, cross sections and alingment tools are under development for now. Any suggestions are welcome.

[ymg]

d.valkanov,

Command: NETLOAD
Cannot load assembly. Error details: System.IO.FileLoadException: Could not
load file or assembly 'file:///C:\Users\ymg\Downloads\CDT2015.dll' or one of
its dependencies. Operation is not supported. (Exception from HRESULT:
0x80131515)
File name: 'file:///C:\Users\ymg\Downloads\CDT2015.dll' --->
System.NotSupportedException: An attempt was made to load an assembly from a
network location which would have caused the assembly to be sandboxed in
previous versions of the .NET Framework. This release of the .NET Framework
does not enable CAS policy by default, so this load may be dangerous. If this
load is not intended to sandbox the assembly, please enable the
loadFromRemoteSources switch. See http://go.microsoft.com/fwlink/?LinkId=155569
for more information.
   at System.Reflection.RuntimeAssembly._nLoad(AssemblyName fileName, String
codeBase, Evidence assemblySecurity, RuntimeAssembly locationHint,
StackCrawlMark& stackMark, IntPtr pPrivHostBinder, Boolean throwOnFileNotFound,
Boolean forIntrospection, Boolean suppressSecurityChecks)
   at System.Reflection.RuntimeAssembly.InternalLoadAssemblyName(AssemblyName
assemblyRef, Evidence assemblySecurity, RuntimeAssembly reqAssembly,
StackCrawlMark& stackMark, IntPtr pPrivHostBinder, Boolean throwOnFileNotFound,
Boolean forIntrospection, Boolean suppressSecurityChecks)
   at System.Reflection.RuntimeAssembly.InternalLoadAssemblyName(AssemblyName
assemblyRef, Evidence assemblySecurity, RuntimeAssembly reqAssembly,
StackCrawlMark& stackMark, Boolean throwOnFileNotFound, Boolean
forIntrospection, Boolean suppressSecurityChecks)
   at System.Reflection.RuntimeAssembly.InternalLoadFrom(String assemblyFile,
Evidence securityEvidence, Byte[] hashValue, AssemblyHashAlgorithm
hashAlgorithm, Boolean forIntrospection, Boolean suppressSecurityChecks,
StackCrawlMark& stackMark)
   at System.Reflection.Assembly.LoadFrom(String assemblyFile)
   at Autodesk.AutoCAD.Runtime.ExtensionLoader.Load(String fileName)
   at loadmgd()

Did not work for me as I am on acad 2012.

Would I've liked to see it run.

I did implement sloan in Autolisp, but it was too slow
to be usable.

ymg

[d.valkanov]

Well, it it takes  about 7 secs. for 1 000 000 points, 23 secs for saving TIN binary file and 87 secs. for contouring. It is actually C# application, but it is developed for 2016. I saw your work and I was  really impressed. It is really great and may be the fastest written in Lisp. I could not find how to implement boundary and holes. There is a command like XSHAPE or something like that, but I am not sure it is for boundary or holes. What is the algorithm in your application? What command should I use for boundary and holes? About smoothing I will suggest Bezier curves. See http://www.antigrain.com/research/bezier_interpolation/
https://en.wikipedia.org/wiki/B%C3%A9zier_curve
It looks much more natural.
Cheers

[d.valkanov]

No way to run on Acad2012. It is written for 2016. I can upload you the source code. CDT to start

[ymg]

d. valkanov,

Isn't it possible to compile it for the earlier versions ?

I understand that it will certainly be much faster,
however at the price of having to recompile at every new versions.

For big triangulation you certainly need to go that route. Either Sloan
or Triangle by Shewchuk.

Still prefer autolisp for my needs.

Xshape is simply to draw a boundary around a point cloud .

For island clean-up you use a closed boundaries (Breaklines), and clean the one
that goes ccw or have no points inside.  Depends on your convention
The outside boundary would need to go clockwise, If you have island
ccw.

ymg

[geobid]

Hi
I have a proposal for a small modification of your procedure:

(defun get_extcont ( / i l r)

(defun LsDiff ( l1 l2)
  (vl-remove-if '(lambda ( x ) (member x l2)) l1)
)
   
   (setq i 0  l nil)
   (foreach n nl
      (if (not (car   n)) (setq l (cons (nth i el) l)))
      (if (not (cadr  n)) (setq l (cons (nth (1+ i) el) l)))
      (if (not (caddr n)) (setq l (cons (nth (+ 2 i) el) l)))
      (setq i (+ 3 i))
   )   
       
   (setq l (reverse l))
   (setq r (reverse (car l)))
   (while (setq a (assoc (car r) (cdr l)))
      (setq r (cons (cadr a) r)
            l (lsDiff l (list a))
      )
   )
   (reverse r)
)


It seems to me that there will be more stable and more resistant to errors.
Of course, the decision is yours.

[ymg]

Geobid,

Your method is about 70% slower.

See below, on a 500 points tin with 984 3dfaces.

(benchmark '((get_extcont) (get_extcont1)))
Benchmarking .............Elapsed milliseconds / relative speed for 1024 iteration(s):

    (GET_EXTCONT)......1934 / 1.69 <fastest>
    (GET_EXTCONT1).....3260 / 1.00 <slowest>
 

However as the list grows the speed penalty gets smaller

For  5000 points and  9980 3DFACES:

(benchmark '((get_extcont) (get_extcont1)))
Benchmarking .........Elapsed milliseconds / relative speed for 64 iteration(s):

    (GET_EXTCONT)......1185 / 1.36 <fastest>
    (GET_EXTCONT1).....1607 / 1.00 <slowest>

One thing I did not check is If the ECO routine by VVA or a the similar one
by Lee Mac would be faster.

Generating the Neighbour List is actually quite lenghty.

ymg

[d.valkanov]

OK. It will be soon. I found ObjectARX 2012, but there some differences in functions for registries, so I have to fix it. May be next week. I am on holiday for now. My question is how to select outer boundary or inner hole as a polyline* or I have to clear triangles manually?

[ymg]

d.valkanov,

Currently, island are defined by closed polylines breaklines that have no points inside.
I use the breakline as a window polygon and do a selection for points or blocks.
These island are cleaned of triangle inside of it.

For the external boundaries we could also do the same, that is if a closed breakline contains
all the points, it would means that we would delete all the triangles which have edges outside
of it.

A small added difficulty is that sometimes a breakline introduce new points which are not
physically there but are part of the point list pl.  We could work around this by keeping a
separate list of the additionnal point or simply how many there are when we creates the point list.

The outside boundary is not currently implemented, but I lean on the above definition.

ymg

[rw2691]

YMG,

Formerly the depression code was hanging or crashing.

In it you evaluate (< zb za) and (> za zb).

What happens with (= zb za) situations?

Could that hang a repeat or while, or mismatch data?

Rick

[d.valkanov]

d.valkanov,

Have a look at this code:

        private void MarkTriangle(int pa, int pb, int tNum, Polygon boundary, RemoveTriangles In_Or_Out)
        {
            if (tNum != -1)
            {
                // Current triangle
                Triangle t = triangles[tNum];

                // Find which edge passes through points pa and pb
                int pc = t.GetThirdVertex(pa, pb);

                // Centre of gravity
                Point mc = ONE_THIRD * (points[pa] + points[pb] + points[pc]);

                // If current triangle is not marked and ...
                if (!t.ToBeDeleted)
                {
                    bool inside; // Point is inside or outside of polygon

                    if (In_Or_Out == RemoveTriangles.Inside)     // If In_Or_Out is true it is a hole in the mesh
                        inside = boundary.InPolygon(mc);         // mc MUST be inside of hole
                    else           // otherwise it is a outer boundary.
                        inside = !boundary.InPolygon(mc);

                    // centre of gravity is outside of boubdary
                    if (inside)
                    {
                        // Mark triangle to be deleted.
                        t.ToBeDeleted = true;

                        // Get opposed triangle of pa
                        int tOpp = t.OpposedTriangle(pa);

                        // Recursively mark neighbour.
                        MarkTriangle(pb, pc, tOpp, boundary, In_Or_Out);

                        // Get opposed triangle of pb
                        tOpp = t.OpposedTriangle(pb);

                        // Recursively mark neighbour.
                        MarkTriangle(pc, pa, tOpp, boundary, In_Or_Out);
                    }

                }
            }
        }

[ymg]

d.valkanov,

Not sure your code accounts for a closed breaklines with point inside
that is simply a breakline which we do not wnt to clean ?

ymg

[ymg]

Rick,

I check for closed polyline that are lower only and change them to depression.

When a polyline does not contain any polyline, we check if there are points in it. if the z of selecting entity is higher we got a sink else it is a peak.

On my side I've never crash in the depression checking so I cannot duplicate
the behaviour you describe.

But maybe the code you have is different from the latest ?

Code - Auto/Visual Lisp: [Select]

1. (if (= godep "1")
2.       (progn
3.          (setq ti (time))
4.  
5.          (or
6.             (tblsearch "ltype" "Depression")
7.             (mk_linetype
8.                "Depression"
9.                "____|____|____|____|____|____|__"
10.                "A,.5,-.005,[\"|\",STANDARD,S=.06,R=0.0,X=-0.01125,Y=-.0725],-.005"
11.             )
12.          )
13.  
14.          (setq lgclosed  (length ccont))
15.  
16.          (if (< lgclosed 10)
17.             (setq div   lgclosed
18.                   count 1
19.             )
20.             (setq div   (/ lgclosed 10)
21.                   count div
22.             )
23.          )
24.          (acet-ui-progress "Depression Processing:"
25.                            (if (< lgclosed 10)
26.                               lgclosed
27.                               10
28.                            )
29.          )
30.  
31.          (while (setq ena (caar ccont))
32.             (setq za      (cadar ccont)
33.                   enta    (entget ena)
34.                   fence   (distinct (caddar ccont))
35.                   aisdepr t
36.                   ccont  (cdr ccont)
37.             )
38.             (if (setq ssb (ssget "_WP" fence '((0 . "LWPOLYLINE"))))
39.                
40.                (repeat (setq i (sslength ssb))
41.                   (setq enb  (ssname ssb (setq i (1- i)))
42.                         entb (entget enb)
43.                         zb   (cdr (assoc 38 entb))
44.                   )
45.                   (if (< zb za)
46.                      (progn
47.                         ;(or first (setq first t ** (chg2depr enta)))
48.                         (chg2depr entb)
49.                         (vl-remove (assoc enb ccont) ccont)
50.                      )
51.                      (setq aisdepr nil)
52.                   )
53.                )
54.  
55.          ;                                                                    ;
56.          ; Else Contour Had No Other Contour Inside, then Need to Check       ;
57.          ; If Contour Contains Only Points or Blocks that Are Lower than      ;
58.          ; Itself, Then It is a Depression Contour Else It is a Peak          ;
59.          ; In a next revision I might add some kind of report about the       ;
60.          ; Volume of the Depression. sink is lowest point in a depression,    ;
61.          ; peak could be used to mark spot elevation on Contour map.          ;
62.          ;                                                                    ;
63.  
64.                (progn
65.  
66.                   (if (setq ssb (ssget "_WP"
67.                                        fence
68.                                        '((-4 . "<OR")
69.                                          (0 . "POINT")
70.                                          (0 . "INSERT")
71.                                          (-4 . "OR>")
72.                                         )
73.                                 )
74.                       )
75.                      (progn
76.                         (setq ps nil)
77.                         (repeat (setq i (sslength ssb))
78.                            (setq enb  (ssname ssb (setq i (1- i)))
79.                                  entb (entget enb)
80.                                  ps   (cons (cdr (assoc 10 entb)) ps)
81.                            )
82.                         )
83.                         (setq ps (vl-sort ps '(lambda (a b) (< (caddr a) (caddr b))))
84.                               zb (caddar ps)
85.                         )
86.                         (if (> za zb)
87.                            (setq sink (car ps))
88.                       ; Else we have a Peak                                   ;
89.                            (setq peak (last ps) aisdepr nil)
90.                         )
91.                      )
92.                   )
93.                )
94.             )
95.             (if aisdepr (chg2depr enta))
96.             (progress)
97.          )
98.          (acet-ui-progress)
99.  
100.          (princ (strcat "\nDEPRESSION " version " - Elapsed time: " (rtos (- (time) ti) 2 3) " secs, " (itoa lgclosed) " Closed Poly\n" ))        
101.       )
102.  

ymg

[d.valkanov]

It is exactly for holes and boundary

private void ClearTriangles(Polygon polygon, RemoveTriangles InOut)
        {
            if (polygon.Count > 0)
            {
                // Start point of boundary
                Point a = polygon.polygonPoints[0];
                int pa = points.IndexOf(a);

                // Find the closest triangle
                int tStart = ChooseStartingTriangle(a);

                // Triangle that contains point pa.
                TriangleEdge edge;
                BarycentricWalk(ref tStart, a); //, out edge);

                int bc = polygon.Count;

                for (int i = 1; i < bc; i++)
                {
                    Point b = polygon.polygonPoints;

                    int pb = points.IndexOf(b);

                    int pc, tOpp;     // Third vertex of triangle

                    Triangle t = triangles[tStart];     // Current triangle

                    // Find triangle containing both points pa and pb.

                    if (!t.Contains(pb))
                    {
                        int r, l;

                        if (pa == t.a)          // pa is vertex a of triangle
                        {
                            r = t.b; l = t.c;   // CCW order
                        }
                        else if (pa == t.b)     // pa is vertex b
                        {
                            r = t.c; l = t.a;   // CCW order
                        }
                        else                    // pa is vertex c
                        {
                            r = t.a; l = t.b;   // CCW order
                        }

                        int tOld = tStart;

                        int tmp;

                        while (!t.Contains(pb))
                        {
                            GetNextCW(triangles[tOld], l, pa, out tStart, out tmp);

                            r = l;
                            t = triangles[tStart];
                            l = tmp;
                            tOld = tStart;
                        }
                        // Now tStart points to the one of triangles passing through pa and pb.
                    }

                    edge = t.GetEdge(pa, pb);

                    // Get third vertex of triangle
                    if (edge == TriangleEdge.ab)
                    {
                        pc = t.c; tOpp = t.ab;
                    }
                    else if (edge == TriangleEdge.bc)
                    {
                        pc = t.a; tOpp = t.bc;
                    }
                    else
                    {
                        pc = t.b; tOpp = t.ca;
                    }

                    if (CommonTools.Orient2D(a, b, points[pc]) == Orientation.CCW)    // pc is on the left side of the boundary
                    {
                        if (tOpp != -1)
                            MarkTriangle(pb, pa, tOpp, polygon, InOut);
                    }
                    else
                        MarkTriangle(pb, pa, tStart, polygon, InOut);

                    // a = b;
                    pa = pb;
                }
            }
        }