Constrained Delauney Triangulation

[SioulisChris]

I read into the following thread at AutoDesk forum:
https://forums.autodesk.com/t5/visual-lisp-autolisp-and-general/lsp-for-topo-survey-boundary/td-p/14013405.
At message 11, user autoid374ceb4990 asks about "breaklines". What he actually talks about is "Constrained Delauney Triangulation - CDT". This is a standard feature of Softdesk (and Land Development and so on) products.

This algorithm would fit very nicely with the Concave Hull, Daniel has already implemented. The generated boundary could be used as a constraint to the triangulation. Also other polylines could be used as internal islands and breaklines.

A C++ project on the subject is https://github.com/artem-ogre/CDT. Could it be used by PyRX?

[CEXT-Dan]

Thanks, looks like a cool library, I’ll have a look. Q. it looks like it already has python bindings, have you tried it in CAD?

It’s an interesting topic, honestly, a bit over my head. I had added these in hopes it would be useful.
Ge.Delaunator (triangles)
Ge.Point3dTree (KDTree)
Ge.Point3dArray.concaveHull

I don’t know if you have seen the Delaunator sample, its pretty powerful in that it gives you triangles and halfedges, so it’s possible to do rudimentary constraints if you already have a polyline boundary(s). Just eliminate the triangle if the edge crosses the polyline, or the circumcircle of the triangle falls in an unwanted area.

In the meantime, I will play around with it in C++.

[SioulisChris]

"Just eliminate the triangle if the edge crosses the polyline, or the circumcircle of the triangle falls in an unwanted area."
Delaunator is GREAT. I use concaveHull a lot. The benefit of the constraints is exactly that: We define the constraints prior to Delaney surface generation. A(n old school) surveyor detects linear features in the field and takes notes. At the CAD he connects the dots and creates polylines for the linear features. These linear features (ponds, walls, curbs, concave hull) are used for the surface definition, together with the point groups. If done properly, there is no need for surface modification after surface generation. No edge deletion or flipping is required, a real time saver.

The biggest disadvantage of CDT is that it deals with 2D points. Your approach with Point3D is the proper one.

[CEXT-Dan]

Okay, I’ll see if I can get it to compile in PyRx

[CEXT-Dan]

I got it to work, caveat, you will need to validate for Intersecting constraint edges, if I add IntersectingConstraintEdges::TryResolve, CDT will add points, invalidating the indexes. In this case, we would no longer be able to map to 3d points.

I don’t know how much of the API you need, I only added a utility method

    #Ge.CTD.triangulate

    def triangulate(
        points: Collection[PyGe.Point3d],
        edges: Collection[tuple[int, int]],
        remove_holes: bool = True,
        /,
    ) -> list: ...

sample:

import traceback
from pyrx import Ap, Ax, Db, Ed, Ge, Rx, Gi

def selectPoly(pntmap: dict[Ge.Point3d, int]):
    ps, id1, _ = Ed.Editor.entSel("\Select Polyline", Db.Polyline.desc())
    if ps != Ed.PromptStatus.eOk:
        raise RuntimeError("oof:")
    pl = Db.Polyline(id1)
    plpnts = pl.toPoint3dArray()
    edgedata = []
    for p in plpnts:
        edgedata.append(pntmap[p])
    return edgedata

@Ap.Command()
def doit():
    try:

        ps, ss = Ed.Editor.select([(0, "POINT")])
        if ps != Ed.PromptStatus.eOk:
            raise RuntimeError("oof:")

        edgedata = []
        pnts = Ge.Point3dArray([Db.Point(id).position() for id in ss.objectIds()])
        
        #use a dict for fast index lookup
        pntmap = {}
        for i, p in enumerate(pnts):
            pntmap[p] = i
        
        result = selectPoly(pntmap)
        edgedata.extend(list(zip(result, result[1:] + [result[0]])))

        result = selectPoly(pntmap)
        edgedata.extend(list(zip(result, result[1:] + [result[0]])))

        faces = []
        for _a, _b, _c in Ge.CTD.triangulate(pnts, edgedata, True):
            faces.append(Db.Face(pnts[_a], pnts[_b], pnts[_c]))

        db = Db.curDb()
        db.addToModelspace(faces)

    except Exception:
        print(traceback.format_exc())

[CEXT-Dan]

fixed that bad triangle above, CDT will use CAD's built in tolerance (AcGeContext::gTol.equalPoint())

same test as above but 3d:

[CEXT-Dan]

Commit 363df60 preview build, you can grab the new ARX file, and probably need PyGe.pyi

[SioulisChris]

Great job Daniel. It's alive!!

Upper part with false argument / Lower part with true.

The initial algorithm/library needs two upgrades:

• Not all Breaklines in Civil/Survey are to be treated as closed. A surveyor expects the following result.

• The API should allow to input both types of breaklines.

[CEXT-Dan]

Ok I’ll play around with this. I see you have played with the python version of this, were you able to get what you wanted from this in 2d? if so it would be very helpful if you have a sample. I’m not familiar with this library at all, so I’m just guessing.

Note: I don’t know if we can use (Conforming Delaunay triangulation), the algorithm generates new points, we would use the Z mapping, but maybe there’s a way to compute this

[CEXT-Dan]

made some progress, generating a new build

[CEXT-Dan]

Hi I committed a new build to the repository; I also added a new sample

I changed the bool to an enum, and triangulate returns a new set of points along with the triangle indices. Apparently flag CDTOpts.KConformToEdges adds new Steiner Points, which is 2d, I have to compute the Z value, hopefully its close.

The enum is bitwise

class CDTOpts(_BoostPythonEnum):
    kNone: ClassVar[Self]  # 0
    KEraseSuperTriangle: ClassVar[Self]  # 1
    kEraseOuterTriangles: ClassVar[Self]  # 2
    kEraseOuterTrianglesAndHoles: ClassVar[Self]  # 4
    KInsertEdges: ClassVar[Self]  # 8
    KConformToEdges: ClassVar[Self]  # 16

[CEXT-Dan]

Hi, while the last release has the a beta implementation, I’ve since refactored the code a bit.
Mainly for performance in edge searching to compute the Z. I’ll probably make another release next week.
Meanwhile, I’ll try to come up with some sort of unit test, scavenge some from CDT and add in Z values