Delaunay triangulation

David Williams · Post by **David Williams** » Wed 13 Jan 2021, 12:53

Okay, here's a first attempt using a naive brute force approach which takes a whopping 7.6 seconds on my Core-i7 laptop to triangulate only 100 points. The trouble is, I may need to triangulate potentially thousands of points!

Perhaps someone here could kindly come up with something significantly faster (which I can subsequently and shamelessly nick for a game project) without resorting to assembly language?

(I shall resort to C or assembler if I really have to...)

Code: Select all

      REM Delaunay triangulation (slow!)
      REM by DW

      MODE 12 : OFF
      ScrW% = @vdu%!208
      ScrH% = @vdu%!212

      nPts% = 100
      DIM p(nPts%-1,1)

      minDist% = 50 : REM No points can be closer than this

      R% = RND(-123456)

      REM Generate random points (subject to specified minimum distance):
      FOR I% = 0 TO nPts%-1
        IF I% = 0 THEN
          p(I%,0) = ScrW%*RND(1)
          p(I%,1) = ScrH%*RND(1)
        ELSE
          REPEAT
            x = ScrW%*RND(1)
            y = ScrH%*RND(1)
            F% = FALSE
            FOR J% = 0 TO I%-1
              IF SQR( (x - p(J%,0))^2 + (y - p(J%,1))^2 ) < minDist% THEN
                F% = TRUE
                EXIT FOR
              ENDIF
            NEXT J%
          UNTIL NOT F%
          p(I%,0) = x
          p(I%,1) = y
        ENDIF
      NEXT I%

      REM Draw all the points:
      GCOL 7
      FOR I% = 0 TO nPts%-1
        CIRCLE FILL 2*p(I%,0), 2*p(I%,1), 8
      NEXT I%

      REM Do the triangulation:
      GCOL 2
      TIME = 0
      FOR K% = 0 TO nPts%-3
        FOR J% = K%+1 TO nPts%-2
          FOR I% = J%+1 TO nPts%-1
            PROCcalc_circle( p(I%,0), p(I%,1), p(J%,0), p(J%,1), p(K%,0), p(K%,1), cx, cy, r_sq )
            F% = FALSE
            FOR L% = 0 TO nPts%-1
              IF L%<>I% AND L%<>J% AND L%<>K% THEN
                dx = p(L%,0) - cx
                dy = p(L%,1) - cy
                d_sq = dx*dx + dy*dy
                IF d_sq < r_sq THEN
                  F% = TRUE
                  EXIT FOR
                ENDIF
              ENDIF
            NEXT L%
            IF NOT F% THEN
              REM Draw triangle:
              LINE 2*p(I%,0), 2*p(I%,1), 2*p(J%,0), 2*p(J%,1)
              LINE 2*p(J%,0), 2*p(J%,1), 2*p(K%,0), 2*p(K%,1)
              LINE 2*p(K%,0), 2*p(K%,1), 2*p(I%,0), 2*p(I%,1)
            ENDIF
          NEXT I%
        NEXT J%
      NEXT K%

      T% = TIME

      REM Redraw all the points:
      GCOL 15
      FOR I% = 0 TO nPts%-1
        CIRCLE FILL 2*p(I%,0), 2*p(I%,1), 8
      NEXT I%

      PRINT "Finished!"'
      PRINT "That took "; T%/100; " seconds"''
      END

      DEF PROCcalc_circle( x0, y0, x1, y1, x2, y2, RETURN cx, RETURN cy, RETURN r_sq )
      LOCAL m1, m2, n1, n2, midx1, midy1, midx2, midy2, c1, c2
      m1 = (y1 - y0) / (x1 - x0)
      m2 = (y2 - y0) / (x2 - x0)
      n1 = -1 / m1
      n2 = -1 / m2
      midx1 = (x0 + x1)/2
      midy1 = (y0 + y1)/2
      midx2 = (x0 + x2)/2
      midy2 = (y0 + y2)/2
      c1 = midy1 - n1*midx1
      c2 = midy2 - n2*midx2
      cx = (c2 - c1) / (n1 - n2)
      cy = cx*n1 + c1
      r_sq = (cx - x0)^2 + (cy - y0)^2
      ENDPROC

David Williams · Post by **David Williams** » Wed 13 Jan 2021, 13:18

Using a revised PROCcalc_circle subroutine, this version lops about 2 seconds compared with the original version. One alternative approach to computing the circumcircle of a triangle does away with at least 2 divisions, but I don't think it will make much difference in terms of speed.

Code: Select all

      REM Delaunay triangulation (slow!)
      REM by DW

      MODE 12 : OFF
      ScrW% = @vdu%!208
      ScrH% = @vdu%!212

      nPts% = 100
      DIM p(nPts%-1,1)

      minDist% = 50 : REM No points can be closer than this

      R% = RND(-123456)

      REM Generate random points (subject to specified minimum distance):
      FOR I% = 0 TO nPts%-1
        IF I% = 0 THEN
          p(I%,0) = ScrW%*RND(1)
          p(I%,1) = ScrH%*RND(1)
        ELSE
          REPEAT
            x = ScrW%*RND(1)
            y = ScrH%*RND(1)
            F% = FALSE
            FOR J% = 0 TO I%-1
              IF SQR( (x - p(J%,0))^2 + (y - p(J%,1))^2 ) < minDist% THEN
                F% = TRUE
                EXIT FOR
              ENDIF
            NEXT J%
          UNTIL NOT F%
          p(I%,0) = x
          p(I%,1) = y
        ENDIF
      NEXT I%

      REM Draw all the points:
      GCOL 7
      FOR I% = 0 TO nPts%-1
        CIRCLE FILL 2*p(I%,0), 2*p(I%,1), 8
      NEXT I%

      REM Do the triangulation:
      GCOL 2
      TIME = 0
      FOR K% = 0 TO nPts%-3
        FOR J% = K%+1 TO nPts%-2
          FOR I% = J%+1 TO nPts%-1
            PROCcalc_circle( p(I%,0), p(I%,1), p(J%,0), p(J%,1), p(K%,0), p(K%,1), cx, cy, r_sq )
            F% = FALSE
            FOR L% = 0 TO nPts%-1
              IF L%<>I% AND L%<>J% AND L%<>K% THEN
                dx = p(L%,0) - cx
                dy = p(L%,1) - cy
                d_sq = dx*dx + dy*dy
                IF d_sq < r_sq THEN
                  F% = TRUE
                  EXIT FOR
                ENDIF
              ENDIF
            NEXT L%
            IF NOT F% THEN
              REM Draw triangle:
              LINE 2*p(I%,0), 2*p(I%,1), 2*p(J%,0), 2*p(J%,1)
              LINE 2*p(J%,0), 2*p(J%,1), 2*p(K%,0), 2*p(K%,1)
              LINE 2*p(K%,0), 2*p(K%,1), 2*p(I%,0), 2*p(I%,1)
            ENDIF
          NEXT I%
        NEXT J%
      NEXT K%

      T% = TIME

      REM Redraw all the points:
      GCOL 15
      FOR I% = 0 TO nPts%-1
        CIRCLE FILL 2*p(I%,0), 2*p(I%,1), 8
      NEXT I%

      PRINT "Finished!"'
      PRINT "That took "; T%/100; " seconds"''
      END

      DEF PROCcalc_circle( x0, y0, x1, y1, x2, y2, RETURN cx, RETURN cy, RETURN r_sq )
      LOCAL n1, n2, mx1, my1, mx2, my2, c1, c2, dx, dy
      n1 = (x1 - x0) / (y1 - y0)
      n2 = (x2 - x0) / (y2 - y0)
      mx1 = x0 + x1
      my1 = y0 + y1
      mx2 = x0 + x2
      my2 = y0 + y2
      c1 = 0.5*(my1 + n1*mx1)
      c2 = 0.5*(my2 + n2*mx2)
      cx = (c1 - c2) / (n1 - n2)
      cy = c1 - cx*n1
      dx = cx - x0
      dy = cy - y0
      r_sq = dx*dx + dy*dy
      ENDPROC

David Williams · Post by **David Williams** » Wed 13 Jan 2021, 16:59

I've made progress by rejecting large triangles (which imply large or huge 'circumcircles'). A pretty obvious speed-up.

DDRM · Post by **DDRM** » Thu 14 Jan 2021, 10:52

I'll have a look at Delaunay triangulation, to see if I can make sense of it....

One immediate thought (which may be completely stupid since I haven't really worked out what is going on!) is, can you use matrix operations in the innermost (L) loop to calculate them all at once, then select the first/best afterwards? I haven't traced the logic: I guess you'll get some 0's for the points that eI,J,K represent? But you can ignore those? Don't see any divide by 0's that will crash it.

Best wishes (and good to hear from you again!),

D

David Williams · Post by **David Williams** » Thu 14 Jan 2021, 16:02

DDRM wrote: ↑Thu 14 Jan 2021, 10:52 I'll have a look at Delaunay triangulation, to see if I can make sense of it....

I've done an integer-only C version, the main optimisation being reducing (drastically) the number of 3-point circumcircles to check, and whilst much faster as you'd expect, the algorithm remains naive (and there's an annoying bug I haven't yet tracked down). The only reason I wanted to triangulate 2D points in the first place is to create fancy, procedurally-generated background images for a vertical scroller I'm working on, although it's far from urgent.

can you use matrix operations in the innermost (L) loop to calculate them all at once

I don't know - it wouldn't surprise me if that was possible, but for my purposes any matrix operations one might resort to must also be compatible with ARM BASIC V (the game is being written to be mostly BB4W/BASIC V agnostic - large chunks of assembler code aside!).

Don't see any divide by 0's that will crash it.

It was sloppy of me to (deliberately) not check for zero (or close-to-zero) divisors, but I did choose to use floats and fractional-valued random coordinates to minimise the risk of division-by-zero errors. Not great software engineering, admittedly.

David Williams · Post by **David Williams** » Thu 14 Jan 2021, 22:28

Delaunay triangulation and Voronoi diagrams are intimately related so, inspired by this pretty animation (one of many) at Shadertoy...

https://www.shadertoy.com/view/ldl3W8

I made this very quickly with BB4W + GfxLib2, purely for a few seconds of entertainment:

http://www.proggies.uk/misc/voronoi2.zip (EXE)

It makes use of a very peculiar GfxLib module called 'PlotZ'.

Source (GfxLib2 needed in order to run it):

Code: Select all

      *ESC OFF

      ON ERROR OSCLI "REFRESH ON" : CLS : ON : REPORT : PRINT " at line ";ERL : END

      MODE 8 : OFF

      ScrW% = @vdu%!208
      ScrH% = @vdu%!212

      nPts% = 50

      INSTALL @lib$+"GFXLIB2" : PROCInitGFXLIB(g{},0)
      INSTALL @lib$+"GFXLIB_modules\PlotZ" : PROCInitModule(0)

      G% = g{}

      PRINT "Please wait a mo..."

      REM Create sprite:
      sprSz% = 512
      sprSz_2% = sprSz% / 2
      spr% = FNmalloc( 4 * sprSz%^2 )
      SYS GFXLIB_SaveAndSetDispVars%, g{}, spr%, sprSz%, sprSz%
      maxRSq% = (sprSz%/2)^2
      FOR Y% = 0 TO sprSz%-1
        FOR X% = 0 TO sprSz%-1
          dx = X% - sprSz_2%
          dy = Y% - sprSz_2%
          r = dx*dx + dy*dy
          IF r < maxRSq% THEN
            f = 1 - r/maxRSq%
            rgba% = FNrgba( 255*f, 255*f*ABSSIN(4*PI*f+0.75), 255, 255*f )
            SYS GFXLIB_PlotPixel%, G%, X%, Y%, rgba%
          ENDIF
        NEXT X%
      NEXT Y%
      SYS GFXLIB_RestoreDispVars%, g{}

      REM Create random positions and velocities:
      DIM p{(nPts%-1) x, y, xv, yv, rgb%}
      FOR I% = 0 TO nPts%-1
        p{(I%)}.x = RND(ScrW%)-1
        p{(I%)}.y = RND(ScrH%)-1
        p{(I%)}.xv = 4*RND(1)*SGN(RND-RND)
        p{(I%)}.yv = 4*RND(1)*SGN(RND-RND)
        p{(I%)}.rgb% = FNrgb( 200+RND(55), 100+RND(50), 32+RND(32) )
      NEXT I%

      REM Animate:
      *REFRESH OFF
      REPEAT
        SYS GFXLIB_Clr%, g{}, 0
        FOR I% = 0 TO nPts%-1
          SYS GFXLIB_PlotZ%, G%, spr%, sprSz%, sprSz%, p{(I%)}.x-0.5*sprSz%, p{(I%)}.y-0.5*sprSz%
        NEXT I%
        FOR I% = 0 TO nPts%-1
          p{(I%)}.x += p{(I%)}.xv
          p{(I%)}.y += p{(I%)}.yv
          IF p{(I%)}.x < 0 OR p{(I%)}.x >= ScrW% THEN p{(I%)}.xv *= -1
          IF p{(I%)}.y < 0 OR p{(I%)}.y >= ScrH% THEN p{(I%)}.yv *= -1
        NEXT I%
        PROCdisplay
      UNTIL FALSE

Anyway, back to other things...

DDRM · Post by **DDRM** » Sat 16 Jan 2021, 15:14

Richard has posted at the groups.io group something that might be relevant/interesting:

https://groups.io/g/bb4w/topic/animated ... 0,79728856

BBC BASIC forum

Delaunay triangulation

Delaunay triangulation

Re: Delaunay triangulation

Re: Delaunay triangulation

Re: Delaunay triangulation

Re: Delaunay triangulation

Re: Delaunay triangulation

Re: Delaunay triangulation