Polygon-Triangle

Auf den ersten Blick scheint das Aufteilen einer Polygonfläche in Dreiecke ein sehr spezielles Programmierproblem darzustellen. Bei näherer Betrachtung eröffnen sich eine Reihe interessanter Anwendungsgebiete. Da sich die so gebildeten Dreiecksflächen, im Gegensatz zum Polygonzug, mathematisch behandeln lassen, kann man z.B. 3D-Routinen vereinfachen, die Fläche eines Polygons berechnen oder einen Button beliebiger Gestalt verwalten.

Die hier vorgestellte Assembler-Routine ist für den Aufruf aus GFA-BASIC 3.xx gedacht. Dank Markus Fritze und Kollegen dürfte der verwendete Assembler (TURBOASS) in jeder Programmsammlung vorhanden sein. Der Quellcode ist zwar gut dokumentiert, aber um eine Konvertierung auf andere Sprachen zu ermöglichen folgt nun die „verständliche Erschlagung“ des Problems.

Bild 1

Das eigentliche Problem bei der Zerlegung des Polygons liegt darin, festzustellen, ob ein zu generierendes Dreieck innerhalb oder außerhalb des Linienzuges liegt. Die Vektorrechnung bietet eine einfache Lösung in Form des vektoriellen Produkts (a x b = c). Der Vektor c steht senkrecht auf den Vektoren a und b, welche ein Rechtssystem bilden. Vereinfacht gesagt wird damit ermittelt, in welche Richtung sich ein Rechtsgewinde (z.B. eine Glühbirne) bewegen würde, wenn man es in die Richtung dreht, in welche sich der Vektor a auf kürzestem Wege auf den Vektor b rotieren läßt (Bild 1). Mit dieser einfachen Methode läßt sich also feststellen, ob die beiden untersuchten Linien (Vektoren) einen konkaven oder konvexen Anteil am Polygon haben. Sollten sie konvex sein, also eine Ausbuchtung bilden, sind es schon mal Kandidaten für ein Dreieck. Dieses Verfahren läßt sich, wie gesagt, nur anwenden, wenn man sich auf einen Umlaufsinn für den Linienzug geeinigt hat (hier wird das übliche Rechtssystem verwendet). Die Größe des Normalvektors c ist also gar nicht von Interesse, lediglich sein Vorzeichen wird benötigt. Die gesuchte z-Komponente des Vektors läßt sich wie folgt berechnen:

mit a=(x10,y10) und b=(x12,y12)
a x b = c = x10y12-x12y10

Sollte c nun negativ sein, steht noch immer nicht fest ob das Dreieck existiert, denn es kann noch von anderen Linien überschnitten werden. Es muß also noch die Bedingung erfüllt sein, daß keiner der Polygoneckpunkte im Dreieck liegt. Dies läßt sich ebenfalls mit der Vektorrechnung überprüfen (Bild 2). Die Linienvektoren, jeweils mit einem positiven Faktor zwischen 0 und 1 multipliziert, beschreiben jeden Punkt auf der von ihnen aufgespannten Raute. Mit der Bedingung, daß die Summe der beiden Faktoren 1 nicht überschreiten darf, begrenzt man die Betrachtung auf die halbe Raute (Dreieck). Die Vektorgleichung lautet mit s und t als Faktoren:

ptx = s * x10 + t * x12
pty = s * y10 + t * y12

Dieses Gleichungssystem mit zwei Unbekannten läßt sich am einfachsten über Determinanten lösen. Die Lösungsdeterminante D ist hierbei angenehmerweise gleich dem Normalvektor c (s.o.).

s = Ds/D mit Ds = ptx * y12 - pty * x12
t = Dt/D mit Dt = pty * x10 - ptx * y10

Da Divisionen zu Rundungsverlusten führen (besonders in der Maschinensprache), sollte man statt der Bedingung s + t <= 1 besser Ds + Dt <= D schreiben.

Das Programm prüft nun mit Hilfe dieser beiden Ansätze für jeweils zwei aufeinanderfolgende Linien, ob sie als Dreiecksseiten aufgenommen werden. Wenn dies der Fall ist, werden die Dreiecksdaten gesichert und der Poygonzug um die beiden Linien verkleinert.

Auf diese Weise wird der Linienzug solange eingekreist, bis er vollständig zerlegt ist - die Routine gibt als Erfolgsbestätigung -1 zurück. Dabei sollte man wissen, daß sich ein Polygon (ohne Überschneidungen) mit n Ecken immer in n-2 Dreiecke zerlegen läßt.

Die beiden zusätzlichen Test-Routinen gestatten es, ein einzelnes Dreieck oder ein ganzes Polygon auf die Inhärenz eines Punktes zu kontrollieren. Das Beispiel-Listing sollte die Verwendung hinreichend erklären.

Bild 2
;************************************************ 
;* polygon_triangle             V1.0 / 28.03.90 *
;* ingo weirauch        (c) 1992 MAXON Computer *
;************************************************ 
;* gfa-basic aufruf:                            *
;* e = C:a_0(L:p_tri,L:p_puf,L:p_x,L:p_y,anz_e) * 
;* f = C:a_1(L:p tri,L:*num,anz_d,xt,yt)        *
;* f = C:a_2(x0,y0,x1,y1,x2,y2,xt,yt)           *
;*                                              *
;* p.tri   = zeiger auf feld form: t&(5,anz d&) * 
;* p_x/p_y = zeiger auf x/y-felder x&(n)/y&(n)  *
;* p_puf   = zeiger auf ram : (anz_e&+3)*4 byte *
;* anz_e   = anzahl der polygonecken - 1 !      *
;* anz_d   = anzahl der dreiecke (anz_e&-3)     *
;* xt/yt   = zu testende koordinaten            *
;* x0-y2   = eckpunkte vom zu testenden dreieck *
;* nun     = zeiger auf card (dreiecksnummer)   *
;* e       = fehler (-1) sonst kein fehler      *
;* f       = -1 = punkt innerhalb des dreiecks  *
;* a_0-a_2 = programmstart (+4 byte / +8 byte)  *
;************************************************

        bra     p_tri           ; polygon triangle
        bra     t_poly          ; punkt im polygon
        bra     t_tri           ; punkt im dreieck
;------------------------------------------------
p_tri:  movem.l 4(sp),a0-a3     ; felder u. puffer
        move.w  20(sp),d7       ; anzahl der ecken
        move.w  d7,d0           ; als zaehler
        movea.l a1,a4           ; pufferadr retten
feld_1: move.w  (a2)+,(a4)+     ; x-koor in puffer
        move.w  (a3)+,(a4)+     ; y-koor in puffer
        dbra    d0,feld_1
        move.l  (a1),(a4)+      ; punkte 0 und 1
        move.l  4(a1),(a4)      ; kopieren
        movea.w d7,a3           ; anzahl der ecken
        adda.w  a3,a3           ; * 4 (als zeiger)
        adda.w  a3,a3
        suba.w  a2,a2           ; ringzaehler init
        subq.w  #2,d7           ; -> dreiecke
        bmi     p_t_e           ; fehler ...
        movea.w d7,a4           ; fehler-zaehler

p_t_lo: lea     0(a1,a2.w),a5   ; datenanfang
        bsr     n_vek           ; vektor-produkt
        tst.l   d2              ; richtung <= 0 ?
        bpi     p_t_0           ; nein .. (konvex)
        move.w  d7,d6           ; retten
        lea     8(a2),a6        ; akt. pos.+2 koor
        bra.s   p_t_2           ; abweisender loop
p_t_1:  bsr     ring_m          ; ringzaehler
        movem.l d0-d7,-(sp)
        movem.w 0(a1,a6.w),d6-d7 ; vergleichspkt 
        sub.w   4(a1,a2.w),d6   ; relativieren
        sub.w   6(a1,a2.w),d7
        bsr     t_t_r           ; punkt testen
        movem.l (sp)+,d0-d7
        bmi.s   p_t_0           ; kein dreieck !
p_t_2:  dbra    d6,p_t_1        ; naechster punkt
        movem.l 0(a1,a2.w),d0-d2; in die 
        movem.l d0-d2,(a0)      ; dreiecksliste
        lea     12(a0),a0       ; neues dreieck
        subq.w  #1,d7           ; 1 punkt weniger
        bmi.s   p_t_e

        movea.w a2,a6           ; akt. pos. retten
        bsr     ring_m          ; pkt. 1 entfernen
        lea     0(a1,a6.w),a5   ; position pkt. 1
        lea     0(a1,a3.w),a6   ; listenende
del_pk: move.l  4(a5),(a5)+     ; verschieben
        cmpa.l  a5,a6           ; ende der liste ?
        bqt.s   del_pk

        subq.w  #4,a3           ; ecken-pointer
        movem.l (a1),d0-d1      ; ende der liste
        movem.l d0-d1,(a6)      ; erneuern
        movea.w d7,a4           ; fehler-zahler
        subq.w  #4,a2           ; alte position

p_t_0:  movea.w a2,a6           ; ringzaehler
        bsr     ring_m          ; erneuern
        movea.w a6,a2
        subq.w  #1,a4           ; fehler wenn < -2
        cmpa.w  #-3,a4
        bge     p_t_lo          ; fehler ?
p_t_e:  move.w  d7,d0           ; zaehler return
        ext.l   d0
        rts
; - - - - - - - - - - - - - - - - - - - - - - - - 
; (a6.w+4) mod a3.w (a3 = akt. anzahl ecken * 4) 
ring_m: addq.w  #4,a6           ; + 1 koordinate
        cmpa.w  a6,a3           ; ecken => a6 ?
        bge.s   ring_e
        suba.w  a6,a6           ; sonst auf 0
ring_e: rts
;------------------------------------------------
t_poly: movem.l 4(sp),a0-a1     ; feldadr und card
        movem.w 12(sp),d5-d7    ; anz d / xt / yt

t_po_l: movem.w d5-d7,-(sp)     ; werte retten
        movem.w (a0)+,d0-d5     ; dreieck laden
        bsr     t_t_d           ; testen
        movem.w (sp)+,d5-d7
        dbmi    d5,t_po_l       ; in die schleife
        sub.w   12(sp),d5       ; rueckrechnen
        neg.w   d5              ; positiv
        move.w  d5,(a1)         ; in die variable
        rts                     ; fertig
;------------------------------------------------
t_tri:  movem.w 4(sp),d0-d7     ; koordinaten-adr.
t_t_d:  sub.w   d2,d6           ; d6= rel. x-koor.
        sub.w   d3,d7           ; d7= rel. y-koor.
        bsr     n_vek0          ; deter. berechnen
t_t_r:  muls    d6,d5           ; 1 unbekannte der
        muls    d7,d4           ; vektor-gl. (s*d)
        sub.l   d4,d5           ; d5 = s*d
        muls    d0,d7           ; 2 unbekannte der
        muls    d1,d6           ; vektor-gl. (t*d)
        sub.l   d6,d7           ; d7 = t*d
        tst.l   d2              ; d negativ ?
        bpl.s   t_t_0           ; Vorzeichen test
        neg.l   d5              ; sonst alle
        neg.l   d7              ; Vorzeichen
        neg.l   d2              ; tauschen
t_t_0:  moveq   #0,d0           ; flag zunaechst 0
        tst.l   d5              ; ausserhalb ?
        bmi.s   t_t_e           ; ja ...
        tst.l   d7              ; ausserhalb ?
        bmi.s   t_t_e           ; ja ...
        add.l   d5,d7           ; <= 1 (hier d)
        cmp.l   d2,d7           ; (dx+dy)/d <= d ?
        bgt.s   t_t_e           ; nein ...
        moveq   #-1,d0          ; im dreieck !
t_t_e:  tst.w   d0              ; flags setzen
        rts
; - - - - - - - - - - - - - - - - - - - - - - - - 
;loesungsmatrix mit den koor. ab (a5.1) berechnen 
n_vek:  movem.w (a5),d0-d5      ; koordinaten
n_vek0: sub.w   d2,d0           ; vektor 0: d0=dx0
        sub.w   d3,d1           ; d1=dy0
        sub.w   d2,d4           ; vektor 1: d4=dx1
        sub.w   d3,d5           ; d5=dy1
        move.w  d0,d2           ; dx0 retten
        move.w  d1,d3           ; dy0 retten
        muls    d5,d2           ; loesungsdeter. d
        muls    d4,d3           ; ( dx+dy <= d )
        sub.l   d3,d2           ; d normalrichtung
        rts 
        end

Listing 1

' **********************************************
' * polygon_triangle - (testprogramm) 29.03.90 *
' *                    in gfa-basic v3.07      *
' * ingo weirauch      (c) 1992 MAXON Computer *
' **********************************************
'
max_pa=200 ! maximale punktzahl
'
INLINE adr%,326
p_tri%=adr% ! polygon -> dreiecke
test_p%=adr%+4 ! polygon ueberwachen
test_tri%=adr%+8 ! ein dreieck testen
'
EVERY 25 GOSUE taste 
DO
    CLS
    GRAPHMODE 3
    PRINT AT(1,1);CHR$(27);"p";SPACE$(80)
    PRINT AT(10,1);"GEBEN SIE EIN POLYGON OHNE "; 
    PRINT "ÜBERSCHNEIDUNGEN IM UHRZEIGERSINN EIN" 
    PRINT AT(1,25);CHR$(27);"w[ESC = ENDE] ";
    PRINT " [LINKE MAUSTASTE = POLYGOPUNKT] ";
    PRINT " [RECHTE MAUSTASTE = POLYGONENDE] ";
    polyline
    anz_d&=anz_e&-2
    ERASE tri&()
    DIM tri&(5,anz_e&)
    a$=SPACE${(anz_e&*3)*4) ! puffer
    tri%=V:tri&(0,0) ! dreiecksliste
    x%=V:x&(0) ! x-koordinaten
    y%=V:y&(0) ! y-koordinaten
    e&=C:p_tri%(L:tri%,L:V:a$,L:x%,L:y%,anz_e&)
    IP e&=-1
        '
        CLR sf
        FOR i&=0 TO anz_da 
            x0&=tri&(0,i&) 
            y0&=tri&(1,i&) 
            x1&=tri&(2,i&) 
            y1&=tri&(3,i&) 
            x2&=tri&(4,i&) 
            y2&=tri&(5,i&)
            a=SQR((x0&-x1&)^2+(y0&-y1&)^2) 
            b=SQR((x1&-x2&)^2+(y1&-y2&)^2) 
            c=SQR((x2&-x0&)^2+(y2&-y0&)^2) 
            s=(a+b+c)*0.5 
            f=SQR(s*(s-a)*(s-b)*(s-c))
            ADD sf,f 
        NEXT i&
        p=sf*100/((320*200)*BSET(0,XBIOS(4)))
        PRINT AT(1,1);" DIE POLYGONFLÄCHE NIMMT ";INT(p*100)/100;" % DER ";
        PRINT "BILDSCHIRMOBERFLÄCHE EIN";SPACES(20) 
        MOUSE ox&,oy&,k& 
        REPEAT
            MOUSE x&,y&,k&
        UNTIL k&=1 OR ABS(x&-ox&)>20 OR ABS(y&-oy&)>20
        '
        PRINT AT(1,1);"ERR: ";e&;SPACE$(70)
        PRINT AT(12,1);"DIE DREIECKSNUMMER LAUTET:" 
        REPEAT
        UNTIL MOUSEK=0 
        nn&=-1
        '
        REPEAT
            MOUSE mx&,my&,k& 
            tri%=V:tri&(0,0)
            IF C:test_p%(L:tri%,L:*na,anz_d&,mx&,my&)
                IF n&<>nn OR anz_d&=0 
                    PRINT AT(39,1);n&;"        ""
                    GRAPKMODE (k&+1)*3 
                    IF nn&>-1
                        draw_tri(nn&)
                    ENDIF 
                    nn&=n& 
                    draw_tri(n&)
                ENDIF 
                nh&=n& 
            ELSE IF nh&<>n& 
                nh&=na
                PRINT AT(39,1);"AUßERHALB"
            ENDIF 
        UNTIL k&=2 
    ELSE 
        CLS
        PRINT AT(30,12);"EINGABEFEHLER !!!!"
        PRINT AT(1,14);"?? ";
        PRINT "FALSCHER UMLAUFSINN / ";
        PRINT "ÜBERSCHNEIDUNGEN ! ";
        PRINT "DOPPELTE PUNKTE (Z.B. START=ENDE) ???"" 
        ~INP(2)
    ENDIF
LOOP
'
PROCEDURE draw_tri(x&)
    ERASE xh&(),yh&()
    DIM xh&(2),yh&(2)
    FOR i&=0 TO 2
        xh&(i&)=tri&(i&*2,x&) 
        yh&(i&)=tri&(i&*2+1,x&)
    NEXT i&
    POLYFILL 3,xh&(),yh&()
RETURN
'
PROCEDURE polyline
    LOCAL n&,m&,z&,key&,t& 
    ERASE x&(),y&()
    DIM x&(max_p&),y&(max_p&)
    CLR z& 
    REPEAT
        MOUSE x&(z&),y&(z&),k& 
    UNTIL k&=1 
    n&=x&(z&) 
    m&=y&(z&)
    REPEAT
        INC z& 
        key&=1 
        REPEAT 
            REPEAT
                MOUSE x&(z&),y&(z&),k& 
                IF x&(z&)<>n& OR y&(z&)<>m& 
                    LINE x&(z&-1),y&(z&-1),n&,m& 
                    LINE x&(z&-1),y&(z&-1),x&(z&),y&(z&) 
                    PRINT AT(1,1);"E: ";z& 
                    SHOWM 
                    n&=x&(z&) 
                    m&=y&(z&)
                ENDIF 
            UNTIL k&<key& 
            key&=3 
        UNTIL k&>0 
    UNTIL z&=max_p& OR k&=2 
    DEC z&
    DRAW TO x&(z&),y&(z&) to x&(0),y&(0) 
    anz_e&=z& 
RETURN
'
PROCEDURE taste 
    IF ASC(INKEY$)=27
        PRINT CHR$(27);"q";CHR$(27);"v”
        GRAPHMODE 1 
        END 
        EDIT 
    ENDIF 
RETURN

Listing 2


Ingo Weirauch
Links

Copyright-Bestimmungen: siehe Über diese Seite