PureBasic Forums - English

http://www.purebasic.com
https://www.purebasic.fr/english/

Delaunay triangulation

https://www.purebasic.fr/english/viewtopic.php?t=57820

Page 1 of 1

Delaunay triangulation

Posted: Sun Dec 29, 2013 12:53 pm
by 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

Re: Delaunay triangulation

Posted: Sun Dec 29, 2013 5:39 pm
by einander
Nice!
Thanks for sharing, and happy 2014!

Re: Delaunay triangulation

Posted: Sun Dec 29, 2013 5:54 pm
by IdeasVacuum
Now that is a very educational example 8)

Re: Delaunay triangulation

Posted: Mon Dec 30, 2013 10:32 am
by Samuel
Very nice work, Comtois! :D

The links are nice and informative too.

Re: Delaunay triangulation

Posted: Mon Dec 30, 2013 2:11 pm
by 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

Re: Delaunay triangulation

Posted: Wed Jan 08, 2014 11:28 pm
by 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 :)

Re: Delaunay triangulation

Posted: Fri Jan 10, 2014 7:45 pm
by 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.

Re: Delaunay triangulation

Posted: Tue Jul 07, 2015 10:55 am
by 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
All times are UTC+01:00
Page 1 of 1

Powered by phpBB® Forum Software © phpBB Limited