Mreženje v 2D - Voronoi

Abstrakt

Crowd problems in the cumputer graphics to construct solvations on net final elements. At the most time this elements is a triangles but that triangles have corrected define a normal and that triangles are enough adaptabled edgeing conditional. At beforehand fixed net of knot point can appear a problems with stability of badly generated net. This work is represented a metod which with use Voronoi diagrams and its dual Delauney triangulation

Izvleček

Množica problemov v računalniški grafiki gradi rešitve na mreži končnih elementov. Največkrat so ti elementi trikotniki, saj imajo točno definirano normalo in so dovolj prilagodljivi robnim pogojem. Pri vnaprej podani mreži vozliščnih točk se lahko pojavijo težave s stabilnostjo slabo generirane mreže. V delu je prikazana metoda z uporabo Voronoi diagramov in njegove dualnosti Delauneyeva triangulacije.

Uvod

Za množico točk izdelajte program, ki bo naredil optimalno trikotnižko mrežo v 2D z uporabo Voronoi diagramov in njegovega duala - Delaunay triangulacija. Z uporabo knjižnice PHIGS je potrebno izrisati generirano mrežo trikotnikov.

Kazalo vsebine

Voronoi diagrami
Definicija Voronoi diagrama
Delaunayeva triangulacija
Odnos med Delaunayevo triangulacijo in Voronoi diagramom
Aplikacije povezane z Voronoi diagrami
Principi izvedbe Delaunayeve triangulacije
Princip tretje točke
Računanje polmera največjega praznega kroga
Predstavitev programskega dela naloge
Literatura

Voronoi diagrami

Voronoi diagrami obravnavajo sosedske odnose med točkami niza N: katera točka je kateri najbližja, katera je od katere najbolj oddaljena, itd.

Definicija Voronoi diagrama

Naj bo P = {p1, p2, ..., pn} niz točk na ravnini. Voronoi diagram V(pi) ravnino razdeli tako, da je diagram V(pi) sestavljen iz vseh točk, ki so vsaj tako blizu točki pi, kot so blizu katerikoli drugi točki pk :

Naj bosta na ravnini samo dve točki p1 in p2. Naj bo premica B(p1, p2) = B12, ki pravokotno razdeli segment p1p2 na dva enaka dela. Potem je vsaka točka x na premici B12 enako oddaljena od točke p1, kakor tudi od točke p2. Velja: |p1x| = |p2x| (Slika 1).


Slika 1 Dve točki na ravnini: |p1x| = |p2x|.

Če ležijo na ravnini tri točke (p1, p2, p3), te definirajo trikotnik. Razpolovnice stranic trikornika B12, B23, B31 se sekajo v točki, ki je središče očrtanega kroga (Slika 2). Ta točka ni nujno v notranjosti trikotnika.


Slika 2 Tri točke na ravnini: razpolovnice se sekajo v središču trikotniku očrtanega kroga.

Iz zgoraj omenjenih primerov je razvidno, da imajo razpolovnice Bij pomembno vlogo pri proučevanju odnosov med točkami na ravnini. Naj bo H(pi, pj) zaprta polravnina, ki je omejena z Bij in vsebuje točko pi. Potem se lahko H(pi, pj) razume kot niz vseh točk, ki so bliže točki pi kot točki pj. Kot že omenjeno, je V(pi) niz točk, ki so bliže točki pi, kot katerikoli drugi točki ali z drugimi besedami, to so točke, ki so bliže pi kot p1, bliže pi kot p2, bliže pi kot p3, itd. Upoštevaje navedeno je mogoče napisati enačbo za Vornoi diagram V(pi):

pri čemer je potrebno upoštevati presek preko vseh i in j, za katere velja . Iz enačbe sledi lastnost Voronoi diagrama: območja Voronoi diagrama so konveksna, ker nastanejo kot presečišča večih polravnin. Kadar so ta območja omejena, so to konveksni poligoni. Robovi se imenujejo Voronoi robovi in temena Voronoi temena. Katerakoli točka na Voronoi robu ima dve najbližji točki niza točk na ravnini, Voronoi teme pa ima najmanj tri najbližje take točke.
V primeru štirih točk na ravnini, ki oblikujejo vogale kvadrata, je to Voronoi diagram kot je prikazano na Sliki 3a. Voronoi teme je v tem primeru četrte stopnje. Če se eno točko nekoliko premakne (Slika 3b), postane diagram normalen. Prvi primer je izrojen zaradi centričnega položaja štirih točk na ravnini.


Slika 3 Izrojen (a) in normalen (b) Voronoi diagram.

Voronoi diagram V(pi) za i = 20 točk na ravnini bi izgledal kot na Sliki 4.

Slika 4 Voronoi diagram za i = 20 točk na ravnini.

Delaunayeva triangulacija

Delaunayeva triangulacija nastane takrat, kadar se med seboj povežejo točke na ravnini, ki imajo skupen Voronoi rob. Delaunayeva triangulacija D(P) in Voronoi diagram V(P) predstavljata dvojnost; isti podatki so predstavljeni na dva različna načina.
Na Sliki 4 je prikazan Voronoi diagram za dvajset točk na ravnini. Na Sliki 5 je prikazana Delaunayeva triangulacija za teh istih dvajset točk. Na Sliki 6 sta prikazana Delaunayeva triangulacija in Voronoi diagram skupaj.


Slika 5Delaunayeva triangulacija za iste točke na ravnini kot na Sliki 4.


Slika 6Delaunayeva triangulacija in Voronoi diagram: Sliki 4 in 5 skupaj.

Odnos med Delaunayevo triangulacijo in Voronoi diagramom

Za razumevanje obeh, v računski geometrijij pomembnih struktur D(P) in V(P), je potrebno poznati odnos med njima. Naj bo P nespremenljiv niz točk na ravnini.

Za Delaunayevo triangulacijo veljajo lastnosti:

D1. D(P) je ravnočrtna dvojnost V(P), po definiciji.
D2. D(P) je triangulacija, če ne obstajajo štiri centrične točke na ravnini: vsaka površina je trikotnik. Ploskve D(P) so Delaunayevi trikotniki.
D3. Vsak trikotnik grafa D(P) ustreza temenu grafa V(P).
D4. Vsak rob grafa D(P) ustreza robu grafa V(P).
D5. Vsak vozel grafa D(P) ustreza območju grafa V(P).
D6. Meja D(P) je konveksna lupina
D7. Notranjost trikotnikov D(P) ne vsebuje nobenih točk niza P.

Za Voronoi diagram veljajo lastnosti:

V1. Vsako Voronoi območje V(pi) je konveksno.
V2. V(pi) je neomejen, če je pi na konveksni lupini niza točk P.
V3. Če je Voronoi teme na stičišču V(p1), V(p2) in V(p3), potem je v središču kroga C(v), ki ga določajo točke p1, p2 in p3.
V4. C(v) je očrtan krog pri Delaunayevi triangulaciji, ki ustreza v.
V5. V notranjosti C(v) ni nobene točke niza P.
V6. Če je pj najbližji pi, potem je (pj,pi) rob triagulacije D(P).
V7. Če obstaja krog skozi točki pj in pi, ki ne vsebuje nobenih drugih točk niza P, potem je (pj,pi) rob triangulacije D(P).

Aplikacije povezane z Voronoi diagrami

S konstrukcijo Voronoi diagrama je povezanih več aplikacij: najbližji sosed, maksimiziranje najmanjših kotov pri triangulaciji, največji prazen krog, najmanjše razvejišče, problem trgovskega potnika, itd.

Principi izvedbe Delaunayeve triangulacije

Rekonstrukcija površine je triangulacija, ki temelji na pravilih Delaunayeve triangulacije na ravnini. Optimalni trikotniki se določajo s pomočjo lokalnega postopka. Poznamo naslednje lokalne postopke: največji prazen krog, pregledovanje ravnine, princip tretje točke, itd.

Princip tretje točke

V našem primeru, za določitev optimalne trikotnižke mreže, smo uporabili princip tretje točke.

Začetek pregledovanja je tako rob, ki ga določata najnižja, najbolj desna točka in točka, ki je prva najbližja po naraščujočem kotu glede na os x (Slika 7).


Slika 7 Začetek pregledovanja niza točk N.

Tretja točka je določljiva glede na položaj središča očrtanega kroga. Izmed vseh točk, ki so levo od aktivnega roba rg, je potrebno izbrati tisto, katere odsek med središčem očrtanega kroga in razpoloviščem aktivnega roba je najmanjši, oz. največji, če je središče očrtanega kroga na desni strani aktivnega roba. Na Sliki 7 je to točka tt. Orientiranost robov trikotnika je pomembna pri iskanju naslednjih točk. Naslednji pregledovalni rob za Sliko 7 bi bil (r-tt, tt-g). Za vsak rob posebej je potrebno poiskati tretjo točko in tako naprej, dokler ni trianguliran celoten točk N.

Računanje polmera največjega praznega kroga

Določevanje tretje točke se skrči na določitev središča očrtanega kroga. Izračun se poenostavi, če se vse tri točke prestavi v drug koordinatni sistem (Slika 8)tako, da leži rob rg na osi x in je r(0,0).


Slika 8 Lokalni koordinatni sistem.

Predstavitev programskega dela naloge

Naloga je izvedena v računalniškem programu Fortran. Sestavljena je iz dveh programov.

Prvi program, z imenom Princip tretje točke bere podatke iz datoteke Pod.txt in nato po teoriji tretje točke določi robove optimalne trikotniške mreže. Optimalna trikotnižka mreža je konveksna. Sestavljena je optimalnih trikotnikov, kjer naj bi veljalo pravilo maksimiziranja najmanjšega kota v trikotniku. To pravilo ni uporabljeno direktno, temveč po metodi principu tretje točke, kjer določimo središče očrtanega kroga in središče aktivnega roba (ta postopek ja že opisan v poglavju princip tretje točke).

V tem programu imamo dva podprograma. Podprogram tretja točka določuje tretjo točko aktivnega roba ter izpisuje te točke v datoteko Phigs.txt.
Drugi podprogram, z imenom Krog, določuje središče očrtanega kroga trikotniku.

Drugi program Phigs.for prebere točke robov trikotnikov iz datoteke Phigs.txt ter jih s pomočjo grafične knjižnice PHIGS izriše v optimalno trikotniško mrežo.

program Princip tretje tocke     
      
      real x(1000),y(1000),p(1000,1000),alfa(1000)                       
      real pxd,pyd,pxz,pyz

      open(3,file='pod.txt') 
      open(4,file='phigs.txt')
      open(5,file='brez.txt')
      
        do i=1,1000
          read(3,*,end=10)x(i),y(i)
            p(i,1)=x(i)
            p(i,2)=y(i)
              sttock=i                        
          write(*,*)p(i,1),p(i,2)    
        enddo
10    close(3) 
************************************************************************** 
c     Dolocitev zacetne tocke (pxz,pyz).         
              
        pyz=y(1) 
        pxz=x(1)
      do i=1,sttock-1
         if (p(i+1,2).le.pyz)then
              pxz=p(i+1,1)
              pyz=p(i+1,2)
           if(p(i+1,2).eq.pyz)then
              if(p(i+1,1).lt.pxz)then
                 pxz=p(i+1,1)
                 pyz=p(i+1,2)
              else
                 pxz=pxz
                 pyz=pyz
              endif
           endif  
         endif
      enddo
         write(*,*)'zacetna tocka je',pxz,pyz
***********************************************************         
c     Dolocitev druge tocke prvega trikotnika (pxd,pyd).
         alf=180                           
      do i=1,sttock
         if(x(i).eq.pxz.and.y(i).eq.pyz)then
            go to 20
         endif     
         if(x(i).ge.pxz)then  
            alfa(i)=atan((y(i)-pyz)/(x(i)-pxz+0.000001)) 
         else
            alfa(i)=180+atan((y(i)-pyz)/(x(i)-pxz+0.000001))
         endif
            if(alfa(i).lt.alf)then
               alf=alfa(i) 
               pxd=x(i)
               pyd=y(i)
            else            
               alf=alf
               pxd=pxd
               pyd=pyd  
20          endif
      enddo
         write(*,*)'druga tocka je',pxd,pyd      
**********************************************************    
c     pxz,pyz=zacetna tocka
c     pyd,pyd=druga tocka
         write(4,*)pxz,pyz 
         write(4,*)pxd,pyd

         call tretja_t(pyz,pyd,pyt,pxz,pxd,pxt,sttock,x,y,pycrz)

      do j=1,2*sttock
         xz=pxz
         yz=pyz
         xt=pxt
         yt=pyt
         xd=pxd
         yd=pyd
        
         pxz=pxz1
         pyz=pyz1
         pxt=pxt1
         pyt=pyt1
         pxd=pxd1
         pyd=pyd1
***********************************************************
             pxz1=xz
             pyz1=yz
             pxd1=xt 
             pyd1=yt 

         call tretja_t(pyz1,pyd1,pyt1,pxz1,pxd1,pxt1,sttock,x,y,
     1                pycrz)
         if (pycrz.eq.1001)then
            write(*,*)'ni leve zunaj 1'
            goto 70
         endif       
            write(5,*)'nekaj 1 t d z',xt,yt,xd,yd,xz,yz
            write(5,*)'nekaj 1 pt pd pz',pxt1,pyt1,pxd1,pyd1,pxz1,
     1      pyz1,pycrz          
***********************************************************     
70          pxz=xt
            pyz=yt
            pxd=xd
            pyd=yd   
         call tretja_t(pyz,pyd,pyt,pxz,pxd,pxt,sttock,x,y,pycrz)
         if (pycrz.eq.1001)then
            write(*,*)'ni leve zunaj'
               pxz=pxz1
               pyz=pyz1
               pxt=pxt1
               pyt=pyt1
               pxd=pxd1
               pyd=pyd1
         endif       
            write(5,*)'nekaj 2 t d z',xt,yt,xd,yd,xz,yz      
            write(5,*)'nekaj 2 pt pd pz',pxt,pyt,pxd,pyd,pxz,pyz,
     1      pycrz      
      enddo
      stop
100      write(*,*)'nekaj je narobe 100' 
      close(4)
      end
            
**********************************************************     
      SUBROUTINE tretja_t(pyz,pyd,pyt,pxz,pxd,pxt,sttock,x,y,pycrz)       
      
*     subrotina nam izracuna:
*    - kot rot. koord. sistema
*    - doloci lego premice - orientiranost
*    - rotira koord. sistem in vse tocke glede na nov koord. sistem 
*    - glede na levost tock, ocrta krog, sredisce, IZBERE 3. TOCKO

      real d,e,fi,kpremice,x(1000),y(1000),pycrz
      parameter(pi=3.14159265359) 

c     Izracun kota premice zacetnega aktivnega roba (fi)   

         kpremice=(pyz-pyd+0.000001)/(pxz-pxd+0.000001) 
         fi=atan(kpremice)*180/pi 
         
c     Dolocitev lege premice 

      if(pxz.ge.pxd.and.pyz.le.pyd)then
         fi=180+fi
      elseif(pxz.le.pxd.and.pyz.ge.pyd)then
         fi=360+fi  
      elseif(pxz.gt.pxd.and.pyz.gt.pyd)then
         fi=180+fi                 
      elseif(pxz.lt.pxd.and.pyz.lt.pyd)then
         fi=fi    
      endif
         write(*,*)'***************************'
***********************************************************
c     Rotacija koordinatnega sistema

         d=cos(fi*pi/180)
         e=sin(fi*pi/180)
            pxdr=(pxd-pxz)*d+(pyd-pyz)*e
            pydr=-((pxd-pxz)*e)+(pyd-pyz)*d 
         pxzr=(pxz-pxz)*d+(pyz-pyz)*e
         pyzr=-((pxz-pxz)*e)+(pyz-pyz)*d
         
c     Postavitev tock v nov koordinatni sistem. 

         pycrz=1001
      do i=1,sttock
         pxt=x(i)                          ! pxt,pyt=tretja tocka
         pyt=y(i)
            pxtr=(pxt-pxz)*d+(pyt-pyz)*e
            pytr=-(pxt-pxz)*e+(pyt-pyz)*d
         if(pytr.gt.0.001)then
               call krog(pxzr,pyzr,pxdr,pydr,pxtr,pytr,pxcr,pycr)
            if(pycr.lt.pycrz.and.pycr.ne.0.000)then
               pycrz=pycr                             !min rot.y = pycrz
               stevec=i
            endif
         endif 
      enddo  

      if (pycrz.eq.1001)then
        write(*,*)'ni leve od znotraj'
        goto 50
      endif
      
      write(*,*)'min rot. y',pycrz,x(stevec),y(stevec),stevec 
      pxt=x(stevec)
      pyt=y(stevec)      
***********************************************************
*     zapis za PHIGS     

      write(4,*)pxz,pyz
      write(4,*)pxt,pyt
      write(4,*)pxt,pyt
      write(4,*)pxd,pyd
50    return 
      end

***********************************************************
      SUBROUTINE krog(x1,y1,x2,y2,x3,y3,xc,yc)
      
      pr(xs,ys,x1,y1,x2,y2)=(ys-y1)*(x2-x1)-(xs-x1)*(y2-y1)
      
      fp(v1,v2,w1,w2)=((v2*v2-v1*v1)+(w2*w2-w1+w1))
      f(v1,v2,w1,w2)=fp(v1,v2,w1,w2)/(2.0*(v2-v1+0.000001))
      g(v1,v2,w1,w2)=(w2-w1)/(v2-v1+0.000001)
      
      ay(x1,x2,x3,y1,y2,y3)=(f(x1,x3,y1,y3)-f(x1,x2,y1,y2))/
     1                     (g(x1,x3,y1,y3)-g(x1,x2,y1,y2)+0.000001)  
     
      ax(x1,x2,x3,y1,y2,y3)=((x3*x3-x1*x1)+(y3-y1)*
     2      ((y3+y1)-2*ay(x1,x2,x3,y1,y2,y3)))/(2.0*(x3-x1+0.000001))   

      if(pr(x1,y1,x2,y2,x3,y3).eq.0.0) then
         else
             if(x1.ne.x3) then
                xc = ax(x1,x2,x3,y1,y2,y3)
                yc = ay(x1,x2,x3,y1,y2,y3) 
             else
                xc=ax(x2,x1,x3,y2,y1,y3)
                yc=ay(x2,x1,x3,y2,y1,y3) 
             endif
             rk=sqrt((x3-xc)**2+(y3-yc)**2)
          endif 
      return    
      end
                
***********************************************************



program phigs 
        
      include 'phigsdef.f'
      real WindowLimits(4)/0.0,15.0,0.0,15.0/
      real ViewportLimits(4)/0.2,0.8,0.2,0.8/
      real ClipLimits(4)  /0.0,1.0,0.0,1.0/
      real ViewMappingMatrix(3,3)
      integer wkid, ErrorReturn, i,j,k,l
      
      real qx(1000),qy(1000) 

         open(10,file='phigs.txt') 
            do i=1,1000
               read(10,*,end=20)qx(i),qy(i)
                  k=i
            enddo
20       write(*,*)'stevilo ogljisc =',k
      close(10) 

      call popph(2,0)  
         wkid=2
      call popwk(wkid," ",WK151280)
      call pevmm(WindowLimits, ViewportLimits, ErrorReturn, 
     *           ViewMappingMatrix)  
     
      do 50 i=1,3
         do 60 j=1,3
            write (*,*)ViewMappingMatrix(i,j)
60       continue
50    continue   
         
      call psplci(2)
      call pslwsc(5)
      call psln(plsoli)
      call pswkw(wkid,0.0,1.0,0.0,1.0)
      call pschh(0.005) 
        
      call pswkv(wkid,0.0,0.15,0.0,0.12)

******* Izris daljic kot polyline  ********************************** 

      call ppl(k,qx,qy)
         pause
      
      call pclwk(wkid)

c       Zapre PHIGS (zakljuci z graficnim nacinom dela)
c       ( PhigsCLosePHigs )

      call pclph()      
         stop
      end

Literatura

[1] Aleš Drolc, Rekonstrukcija površin, Fakulteta za strojništvo,Diplomska naloga,Ljubljana 1997.

[2] L. Kos, J. Duhovnik, Fakulteta za strojništvo, Ljubljana 1997.

[3] K. Mulmuley,Computational Geomtry: an Introduction Trough Rendomized Algorithms, Prentice -Hall, Inc.,1994

[4] Joseph O'Rourke,Computational Geomtry in C, Cambridge Universty Press, 1993

Vaše mnenje lahko pošljete na:

Andrej Marlin
University of Ljubljana
Faculty of mechanical engineering
CAD lab - LECAD
Askerceva 6
1000 Ljubljana
Slovenia

Telephone: +386 61 1771-436