Delaunay triangulation

[Comtois]

Button left = add point
Button right = Triangulation
Space = init
escape = quit

Code: Select all

;Button left = add point
;Button right = Triangulation
;Space = init
;escape = quit

;http://www.ai.univ-paris8.fr/~audibert/ens/11-Triangulation.pdf
;
;http://www.geom.uiuc.edu/~samuelp/del_project.html

InitSprite()
InitKeyboard()
InitMouse()
OpenScreen(1024, 768, 32, "Delaunay")

#Max = 2500
Global NombrePoints 

Global Dim Points.point(#Max)
Global Dim voisins(#Max,#Max)
Global Dim longueur(#Max)
Global Dim envconv(#Max)

Declare AfficheSouris()
Declare AffichePoints()
Declare troppres(i)
Declare plusprochevoisin(i)
Declare voisinsuivantdroite(i, j)
Declare voisinsuivantgauche(i, j)
Declare Triangulation()
Declare AfficheTriangles()

;- Main loop
Repeat
  
  ClearScreen(0)
  
  ExamineKeyboard()
  ExamineMouse()
  
  If MouseButton(#PB_MouseButton_Left)
    If mem = 0 And NombrePoints < #Max
      mem = 1
      Points(NombrePoints)\x = MouseX()
      Points(NombrePoints)\y = MouseY()
      NombrePoints + 1
    EndIf 
  Else
    mem = 0  
  EndIf     
  
  If MouseButton(#PB_MouseButton_Right)
    Triangulation()
  EndIf  
  
  If KeyboardReleased(#PB_Key_Space)
    NombrePoints = 0
    Global Dim Points.point(#Max)
    Global Dim voisins(#Max,#Max)
    Global Dim longueur(#Max)
    Global Dim envconv(#Max)
  EndIf  
  
  StartDrawing(ScreenOutput())
    AfficheSouris()
    AfficheTriangles()    
    AffichePoints()
  StopDrawing()
    
  FlipBuffers()
Until KeyboardPushed(#PB_Key_Escape)  

End

Procedure AfficheSouris()
  d = 10
  Line(MouseX(),MouseY(), d,  1, RGB(255, 255, 255))
  Line(MouseX(),MouseY(),-d,  1, RGB(255, 255, 255))
  Line(MouseX(),MouseY(), 1,  d, RGB(255, 255, 255))  
  Line(MouseX(),MouseY(), 1, -d, RGB(255, 255, 255))  
EndProcedure

Procedure AffichePoints()
  For i=0 To NombrePoints-1
    Circle(Points(i)\x, Points(i)\y, 2,RGB(255, 0, 0))
  Next
EndProcedure

Procedure plusprochevoisin(i)
  distmin=900000000
  For j=0 To NombrePoints-1 
    If i<>j
      dist2=(Points(i)\x-Points(j)\x)*(Points(i)\x-Points(j)\x)+(Points(i)\y-Points(j)\y)*(Points(i)\y-Points(j)\y)
      If (dist2<distmin) 
        distmin=dist2
        jmin=j
      EndIf
    EndIf
  Next   
  ProcedureReturn jmin
EndProcedure 

Procedure voisinsuivantdroite(i, j)
  numerovoisin=-1
  Protected.d prodscal,longki2,longki,longkj2,longkj,coskij,cosmin,det
  
  cosmin=1
  For k=0 To NombrePoints-1
    If k<>i And k<>j
      v1x=Points(j)\x-Points(i)\x
      v1y=Points(j)\y-Points(i)\y
      v2x=Points(k)\x-Points(i)\x
      v2y=Points(k)\y-Points(i)\y
      det=v1x*v2y-v1y*v2x
      If det<0 ; on cherche un point à droite 
        prodscal=(Points(i)\x-Points(k)\x)*(Points(j)\x-Points(k)\x)+
                 (Points(i)\y-Points(k)\y)*(Points(j)\y-Points(k)\y)
        longki2 =(Points(i)\x-Points(k)\x)*(Points(i)\x-Points(k)\x)+
                 (Points(i)\y-Points(k)\y)*(Points(i)\y-Points(k)\y)
        longkj2 =(Points(j)\x-Points(k)\x)*(Points(j)\x-Points(k)\x)+
                 (Points(j)\y-Points(k)\y)*(Points(j)\y-Points(k)\y)
        longki=Sqr(longki2)
        longkj=Sqr(longkj2)
        coskij=prodscal/(longki*longkj); on veut le cosinus le plus petit possible 
        If (coskij<cosmin) 
          cosmin=coskij 
          numerovoisin=k
        EndIf
      EndIf
    EndIf  
  Next
  
  ProcedureReturn numerovoisin
EndProcedure

Procedure voisinsuivantgauche(i, j)
  numerovoisin=-1
  Protected.d prodscal,longki2,longki,longkj2,longkj,coskij,cosmin,det
  
  cosmin=1
  For k=0 To NombrePoints-1
    If k<>i And k<>j
      v1x=Points(j)\x-Points(i)\x
      v1y=Points(j)\y-Points(i)\y
      v2x=Points(k)\x-Points(i)\x
      v2y=Points(k)\y-Points(i)\y
      det=v1x*v2y-v1y*v2x
      If det>0 ; on cherche un point à gauche 
        prodscal=(Points(i)\x-Points(k)\x)*(Points(j)\x-Points(k)\x)+
                 (Points(i)\y-Points(k)\y)*(Points(j)\y-Points(k)\y)
        longki2=(Points(i)\x-Points(k)\x)*(Points(i)\x-Points(k)\x)+
                (Points(i)\y-Points(k)\y)*(Points(i)\y-Points(k)\y)
        longkj2=(Points(j)\x-Points(k)\x)*(Points(j)\x-Points(k)\x)+
                (Points(j)\y-Points(k)\y)*(Points(j)\y-Points(k)\y)
        longki=Sqr(longki2)
        longkj=Sqr(longkj2)
        coskij=prodscal/(longki*longkj); on veut le cosinus le plus petit possible 
        If (coskij<cosmin) 
          cosmin=coskij 
          numerovoisin=k
        EndIf
      EndIf
    EndIf  
  Next
  
  ProcedureReturn numerovoisin
EndProcedure

Procedure Triangulation()
  For i=0 To NombrePoints-1 ; on prend chaque site 
    ppv=plusprochevoisin(i)
    voisins(i,0)=ppv
    longueur(i)=1
    kk=0
    Repeat 
      vsd=voisinsuivantdroite(i,voisins(i,kk))
      longueur(i) + 1
      voisins(i,kk+1)=vsd
      kk+1
    Until (vsd=ppv Or vsd=-1)
    If vsd=-1
      envconv(i)=1
      Repeat 
        vsg=voisinsuivantgauche(i,voisins(i,0))
        If vsg<>-1
          For kkk=longueur(i)-1 To 0 Step -1
            voisins(i,kkk+1)=voisins(i,kkk)
          Next
          voisins(i,0)=vsg
          longueur(i)+1
        EndIf
      Until(vsg=-1)
    EndIf
  Next
EndProcedure

Procedure AfficheTriangles()
  Color = RGB(0, 0, 255)
  For i=0 To NombrePoints-1 ; dessin des triangles autour de i 
    For ii=0 To longueur(i)-2
      If voisins(i,ii+1)<>-1
        LineXY(Points(i)\x,Points(i)\y,Points(voisins(i,ii))\x,Points(voisins(i,ii))\y, Color)
        LineXY(Points(i)\x,Points(i)\y,Points(voisins(i,ii+1))\x,Points(voisins(i,ii+1))\y,Color)
        LineXY(Points(voisins(i,ii))\x,Points(voisins(i,ii))\y,Points(voisins(i,ii+1))\x,Points(voisins(i,ii+1))\y,Color)
      EndIf
    Next
  Next  
EndProcedure

[einander]

Nice!
Thanks for sharing, and happy 2014!

[IdeasVacuum]

Now that is a very educational example

[Samuel]

Very nice work, Comtois!

The links are nice and informative too.

[applePi]

certainly great example. i will read what is Delaunay triangulation, and will see if we can turn such triangles to a mesh
Happy new year 2014

[Comtois]

google translation :

Note that this algorithm is a constructive proof of the existence and uniqueness of the Delaunay triangulation.
There are other more efficient algorithms. That we present here has the great disadvantage of repeatedly giving the same triangles, as these are in the ranges of each of its three vertices.
Its advantage is to obtain for the price the same convex envelope, and especially to allow easily to deduce the Voronoi diagram, using the Delaunay triangles adjacent surrounding each site (point), the algorithm gives us.

So yes you are right. And sorry for french naming

[Comtois]

In this tutorial, the author uses two methods to calculate the Delaunay triangles.
- Divide and conquer
- Flip algorithm

Sorry this is still in French, but the language used is C.
The code is not too difficult to adapt in PureBasic, unfortunately I do not have time to do yet.

[applePi]

converting such Delaunay triangulation to a mesh we get a funny mesh like Mosaic, picasso will like it, i have added a procedure makeMesh to the Comtois code above to call after finishing the design by pressing 'Z' key to save the mesh as (plane.mesh). it seems hard to texture. and some triangles are duplicated

Code: Select all

Procedure makeMesh()
  CreateMesh(#mesh, #PB_Mesh_TriangleList, #PB_Mesh_Dynamic)
  tu.f:tv.f
  Color = RGB(0, 0, 255)
  For i=0 To NombrePoints-1 ; dessin des triangles autour de i 
   
    For ii=0 To longueur(i)-2
      
      If voisins(i,ii+1)<>-1
        
      MeshVertexPosition(0.02*Points(i)\x,0.02*Points(i)\y, 0)
      MeshVertexNormal(1,1,1)
      MeshVertexColor(RGB(255,255,255))
      MeshVertexTextureCoordinate(tu, tv)
      
      MeshVertexPosition(0.02*Points(voisins(i,ii))\x,0.02*Points(voisins(i,ii))\y,0)
      MeshVertexNormal(1,1,1)
      MeshVertexColor(RGB(255,255,255))
      MeshVertexTextureCoordinate(tu, tv)
      
      MeshVertexPosition(0.02*Points(voisins(i,ii+1))\x,0.02*Points(voisins(i,ii+1))\y, 0)
      MeshVertexNormal(1,1,1)
      MeshVertexColor(RGB(255,255,255))
      MeshVertexTextureCoordinate(tu, tv)
      tu.f+1/longueur(i)
      
      
      MeshFace(ind+2, ind+1, ind)
        
        ind + 3 ; important
      EndIf
     

    Next
    tu=0
    tv+1/NombrePoints
    
  Next
  
  
  FinishMesh(#True)
  
  UpdateMeshBoundingBox(#mesh)
  
  SaveMesh(#mesh, "plane.mesh")  
  
  
EndProcedure

for your convenience download the windowed version of the program
+code to view the 'plane.mesh' with a ball falling over it
+code to load the plane.mesh and jumping to its triangles one by one. the three colored spheres show the direction of every triangle (red > green > yellow, anticlockwise direction). you can see here that some triangles are duplicated
also attached a specimen of plane.mesh and converting it by deepMesh to 3ds + obj + x version (loading modern *.mesh file in deepMesh is possible if you apply some editing to the file: look here http://purebasic.fr/english/viewtopic.p ... 45#p466340
i have added the 3ds, x, obj files if someone are able to texture it.
http://www.2shared.com/file/xDNvdIhG/De ... ation.html (click on the smaller download button at the bottom and not the big one)
Edit: added ShowCursor_(1) so no need for the crossed lines cursor
http://www.2shared.com/file/AL6AY24R/De ... tion2.html