Ideen für Umlaufzahl-Berechnung?

[Johann L. (gjlayde)]

Hi,

für meine Snake on Scope brauche ich an einer Stelle die Umlaufzahl 
einer Kurve um einen Punkt.

Die Kurve ist eine Jordan-Kurve, und es ist zu bestimmen, ob der Punkt 
im Innern der Kurve liegt oder ausserhalb.

Zur Erinnerung: Eine Jordankurve teilt die Ebene in ein Inneres und ein 
Äusseres, sie ist geschlossen, und sie schneidet sich nicht selbst (also 
nicht 8-förmig oder so).

Die Elemente der Kurve sind Geradenstücke in den Koordinatenrichtungen.
Die Kurve ist weder konvex noch sternförmig.

Momentan mache ich das wie im Anhang in snake_umlauf zu sehen: Für alle 
Mauerelemente wird der Winkel bestimmt, unter dem das Element von dem zu 
testenden Punkt aus zu sehen ist.
Diese Winkel werden für alle Segmente aufaddiert.

Für Punkte, die nicht auf der Kurve liegen (was der hier Fall ist), ist 
die so bestimmt Umlaufzahl ein ganzzahliges Vielfaches des Vollwinkels. 
Ein Punkt liegt im Innern wenn die Umlaufzahl 360° ist.

Irgendwei hab ich den Eindruck, ich mach das viel zu kompliziert, aber 
selbst nach einigem Überlegen ist mir kein einfacherer Weg eingefallen, 
das zu programmieren. Es funktioniert zwar, aber die 
Arcustangens-Tabelle nimmt einiges an Platz weg.

Allein diese Umlaufzahl-Bestimmung brauch an die 500 Bytes Flash!

Codetechnisch muss das ganze natürlich -- wie immer -- schnell sein und 
soll wenig Platz brauchen.

Es läuft auf einem ATmega168 @ 24MHz.
Per 37.5 kHz brauch ich einen neuen Pixel für die Röhre (also alle 640 
Ticks), was nicht sooo viel Zeit ist.

Also: Ideen sind gesucht, wie und ob man die Umlaufzahl einfacher 
berechnen kann.

Hier nochmal das wesentliche Codefragment: Es ist über alle 
Kurvensegmente aufzurufen, welche jeweils von P0 bis P1 reichen.

#include "atan2tab.h"

// (A-B)/4 berechnen und so runden, dass das Vorzeichen gleich bleibt.
// ./4 damit die Differenz in 8 Bit (signed) passt
int8_t subr (uint8_t a, uint8_t b)
{
    return (a < b)
        ? (a-b) >> 2
        : -((b-a) >> 2);
}

void ...
{
    // P0: Beginn des Mauerelements
    uint8_t x0 = ...
    uint8_t y0 = ...

    // P1: Ende des Mauerelements
    uint8_t x1 = ...
    uint8_t y1 = ...

    // Richtungsvektoren vom Punkt, fuer den der Umlauf
    // zu bestimmen ist (s->food) zu den Endpunkten des
    // Mauersegments
    // P2 = P0 - Food
    // P3 = P1 - Food
    int8_t x2 = subr (x0, s->food.x);
    int8_t y2 = subr (y0, s->food.y);
    int8_t x3 = subr (x1, s->food.x);
    int8_t y3 = subr (y1, s->food.y);

    // Element so hindrehen, daß P3 auf (K,0) mit K >= 0
    // abgebildet wird. Matrix ist M=((x3,y3),(y3,-x3))
    // P4 = M*P2
    // M*P3 zeigt in positive x-Richtung
    int16_t x4 = x3*x2 + y3*y2;
    int16_t y4 = y3*x2 - x3*y2;

    // Auf 8 Bits (signed) zusammenschieben
    // P4 >>= ?
    while (1)
    {
        if (x4 < 128 && x4 > -128)
            if (y4 < 128 && y4 > -128)
                break;
        x4 >>= 1;
        y4 >>= 1;
    }

    // Auf 4 Bits (signed) zusammenschieben
    // P5 = P4 >> ?
    int8_t x5 = x4;
    int8_t y5 = y4;

    while (1)
    {
        if (x5 < 8 && x5 > -8)
            if (y5 < 8 && y5 > -8)
                break;
        x5 >>= 1;
        y5 >>= 1;
    }

    // Winkel von P5 gegen die positive X-Achse bestimmen
    // Das ist der Winkel, unter dem (P0,P1) von Food aus
    // erscheint.

    // Index in atan2-Tabelle errechnen
    uint8_t x6 = x5+7;
    uint8_t y6 = y5+7;
    uint8_t index = (uint8_t) (x6 + 15*y6);

    // Winkel aus Tabelle lesen (Vollwinkel ~ 100)
    int8_t phi = pgm_read_byte (& atan2tab[index]);

    s->wall_umlauf += phi;
}

Grüß, Johann

[Johann L. (gjlayde)]

Und noch ein Bild.

Die Kurve, für die so eine Berechnung zu machen ist, sieht aus wie die 
gestrichelte Kurve auf dem Bild, der Rand des Spielfeldes.

Die Kurve kann natürlich auch komplizierter aussehen.

[Karl H. (kbuchegg)]

Es gibt mehrere Möglichkeiten für einen Punkt/Polygon Test.
Und ja eine davon ist die Methode des Umlaufwinkels.

Allerdings muss der Winkel nicht sehr exakt berechnet werden. Es reicht 
völlig aus, wenn er so ungefähr hinkommt. Die Summe kann ja nur 2 Werte 
haben: 0 oder 2*PI (oder Vielfache davon).

Es reicht völlig aus, den Gesamt-Winkel mit 4 Bit Genauigkeit zu 
bestimmen und für jeden Einzelwinkel eine 2 Bit Genauigkeit zu haben. 
Effektiv läuft es damit darauf hinaus, Quadranten abzuzählen :-)
Damit ist deine Tangens-Tabelle nicht notwendig.

Anbei eine Implementierung dieses Verfahrens von Kevin Weiler aus dem 
Buch 'Graphics Gems IV'. Der Code darf beliebig weiter verwendet werden.

(Anpassen an deine Datenstruktur musst du ihn noch selber)

[Karl H. (kbuchegg)]

Und das Header File dazu

[Karl H. (kbuchegg)]

Anbei noch eine Testsuite für verschiedene Point in Polygon Verfahren 
aus dem gleichen Buch. Diesmal von Eric Haines. Auch hier wieder: Der 
Code kann beliebig verwendet werden.

[Johann L. (gjlayde)]

Hi, danke!

Also kommt man am Berechnen der Umluuf- resp. Windungszahl nicht 
vorbei...

Meinen Algorithmus hab ich etwas verbessert und Symmetrien von atan 
ausgenutzt, so daß es jetzt nur noch 150 Bytes insgesamt braucht.

Weilers Algorithmus sieht sehr interessant aus!
Den muss ich mir mal näher anschauen. Die unangenehme Division kann man 
rauserweitern, dabei geht es ja nur um einen Vorzeichentest. Und 16-Bit 
Arithmetik reicht hoffentlich, ansonsten wird das zu komplex.

Johann

[yalu (Gast)]

> Also kommt man am Berechnen der Umluuf- resp. Windungszahl nicht
> vorbei...

Also geht es dir nicht primär darum, den Umlaufwinkel zu bestimmen,
sondern tatsächlich herauszufinden, ob sich ein Punkt innmerhalb eines
gegebenen Polygons befindet? Der Thread-Titel und dein Initialpost haben
mich da etwas verunsichert :)

Eine andere Möglichkeit des Punkt-im-Polygon-Tests besteht darin, dass
man abzählt, wie oft eine in dem zu testenden Punkt beginnende
Halbgerade die Polygonlinie durchquert. Ist die Anzahl ungerade, liegt
der Punkt im Innern.

Das Verfahren könnte in deinem Fall ganz geschickt sein, da alle
Polygonseiten Achsenparallel liegen, so dass die Überprüfung der
Durchquerung relativ leicht ist, wenn man die Halbgerade ebenfalls
achsenparallel wählt. Ob es allerdings schneller ist als das von Karl
Heinz genannte Verfahren, weiß ich nicht.

[Karl H. (kbuchegg)]

yalu wrote:
> Das Verfahren könnte in deinem Fall ganz geschickt sein, da alle
> Polygonseiten Achsenparallel liegen,

Ja das hab ich mir auch überlegt, hab den Vorschlag aber dann wieder 
verworfen, weil ich nicht sicher war, ob er die Achsparallelität nicht 
irgendwann aufgeben möchte.
Ansonsten wäre das Schnittpunkte zählen speziell in diesem Fall sicher 
die beste Lösung, weil man so gut wie nichts rechnen muss und letztendes 
alles auf ein paar Vergleiche rausläuft.

Allerdings gibt es da ein paar unangnehmen Sonderfälle, die Beachtung 
brauchen


(vom Punkt ausgehend, nach rechts die Schnittpunkte zählen

  +------------------+
  |                  |
  |    *  +------+   |
  |       |      |   |
  +-------+      +---+


3 'Schnittpunkte'   -> Punkt ist drinnen, wie es sein soll



                   +-----+
                   |     |
       *     +-----+     |
             |           |
             +-----------+

ebenfalls 3 Schnittpunkte, aber diesmal ist der Punkt draussen


Das Problem besteht in der Behandlung eines 'Schnittpunktes', wenn der 
SChnittpunkt mit dem Endpunkt einer Polygonseite identisch ist.
Ich hab das meist so gelöst, dass in so einem Fall der Schnittpunkt nur 
dann zählt, wenn er mit dem unteren (also dem mit dem kleineren y Wert) 
der beiden Endpunkte identisch ist und im anderen Fall nicht.

[yalu (Gast)]

Habe mal etwas reingehackt, weil gerade Sonntag ist :)

#include <stdio.h>
#include <stdint.h>

typedef struct {
  uint8_t x, y;
} point_t, poly_t[];

// Polygondaten: Es ist nur jeder zweite Punkt eingetragen. Ein Punkt ist mit
// seinem Nachfolger durch ein L-Stück verbunden, dessen erster Schnekel
// horizontal und dessen zweiter vertiakl verläuft.

// Beispielpolygon
poly_t poly = {
  {1, 1}, {3, 2}, {4, 1}, {6, 3}, {5, 4},
  {6, 6}, {4, 5}, {3, 6}, {1, 4}, {2, 3}
};

uint8_t isinside(poly_t poly, uint8_t n, point_t point) {
  uint8_t yold, i, cnt=0;

  yold = poly[n-1].y;
  for(i=0; i<n; i++) {
    if(poly[i].x>point.x && ((yold<=point.y) ^ (poly[i].y<=point.y)))
      cnt++;
    yold = poly[i].y;
  }
  return cnt%2;
}

int main(void) {
  uint8_t ins;

  // Beispielaufruf
  ins = isinside(poly, sizeof poly / sizeof poly[0], (point_t){3, 4});
  return 0;
}

Die Funktion isinside liefert eine 1 zurück, wenn der Punkt innerhalb
und eine 0, wenn er außerhalb des Polygons liegt. Linien, die das
Polygon nach links oder unten begrenzen, gehören zum Polygon dazu,
Linien, die es nach oben oder rechts begrenzen, nicht. Wenn ich den
Eingangspost richtig verstanden habe, liegt der Punkt aber nie auf dem
Polygonumfang, so dass dies keine Rolle spielt. Wenn doch, müsste man
die if-Abfrage etwas erweitern, um zu erreichen, dass der Umfang
entweder komplett zum Polygon dazugehört oder komplett nicht.

Das Polygon wird, da die Seiten achsenparallel sind, durch die Hälfte
seiner Eckpunkte vollständig beschrieben (s. Kommentar im Programm).

Die Funktion braucht 74 Bytes Programmspeicher (AVR-GCC 4.2.4) und
sollte auch einigermaßen schnell sein, da pro Polygonpunkt nur einige
wenige Billigoperation ausgeführt werden.

Ich habe sie an folgendem Polygon getestet, bei das Ergebnis für alle
Punkte x, y in [0, 7] gestimmt hat:

6_   ___   ___
5_  |   |_|   |
4_  |_       _|
3_   _|     |_
2_  |    _    |
1_  |___| |___|
0_  
  | | | | | | |
  0 1 2 3 4 5 6

Das heißt natürlich nicht, dass keine Bugs vorhanden wären ;-)

Man kann das Verfahren auch auf Polygone mit nicht achsenparallelen
Seiten erweitern. Dann muss die Schleife zum einen über alle Eckpunkte
laufen, zum anderen ist dann die Abfrage poly[i].x>point.x durch einen
komplizierteren Ausdruck zu ersetzen, der prüft, ob der Punkt links oder
rechts von der (möglicherweisen schiefen) Polygonseite liegt. Dazu sind
einige zusätzliche Additionen und Multiplikationen erforderlich, die
teilweise in 16-Bit gerechnet werden müssen.

[yalu (Gast)]

Kleiner Nachtrag: Weitere 4 Bytes können eingespart werden, wenn man die
Zeilen

      cnt++;
    ...
  return cnt%2;

durch

      cnt=~cnt;
    ...
  return cnt;

ersetzt. Als Funktionswert erhält man jetzt statt der 1 eine 255, was
aber egal sein sollte, da er in einer logischen Abfrage beide als "wahr"
zählen.

[Johann L. (gjlayde)]

yalu schrob:
> Also geht es dir nicht primär darum, den Umlaufwinkel zu bestimmen,
> sondern tatsächlich herauszufinden, ob sich ein Punkt innmerhalb eines
> gegebenen Polygons befindet? Der Thread-Titel und dein Initialpost haben
> mich da etwas verunsichert :)

Für ne geschlossene Kurve ist das ja gleichbedeutend.

> Eine andere Möglichkeit des Punkt-im-Polygon-Tests besteht darin, dass
> man abzählt, wie oft eine in dem zu testenden Punkt beginnende
> Halbgerade die Polygonlinie durchquert. Ist die Anzahl ungerade, liegt
> der Punkt im Innern.

ARRRGH Tischkante-beiss, das ist genau das was ich gesucht hab!

Ich wusste doch daß ich wieder viel zu umständlich dachte...

> Ob es allerdings schneller ist als das von Karl
> Heinz genannte Verfahren, weiß ich nicht.

Ich würde sagen ja, weil man kein * oder / braucht.
Und es ist schneller als mein Tabellen-Verfahren.

Es braucht 60 Bytes freu.

void snake_umlauf (snake_t * s, uint8_t dir, uint8_t len)
{
    if (dir & 1)
        // NORD || SUED
        return;

    // Q: Testpunkt
    // P0/P1: Beginn/Ende des Mauerelements
    uint8_t yq = s->food.y;
    uint8_t y0 = s->draw.y;
    
    // Teststrahl soll nach oben gehen
    if (y0 < yq)
        return;
        
    // Laenge des Mauerelements
    uint8_t mlen = len << WALL_SKIP;
    
    // Trifft der Strahl die Mauer?    
    uint8_t xq = s->food.x;
    uint8_t x0 = s->draw.x;
    uint8_t x1 = x0;
    
    if (dir == SN_E)
    {
        // Ost
        x1 += mlen;
        if (xq < x0 || xq > x1)
            // Nein
            return;
    }
    else
    {
        // West
        x1 -= mlen;
        if (xq < x1 || xq > x0)
            // Nein
            return;
    }
    
    // Ja, der Strahl trifft
    // WEST: +1
    // Ost:  -1
    s->wall_umlauf -= dir -1;
}

...und da sie nur 1x verwendet wird braucht's nach demm Inlinen nur noch 
32 Bytes.

Und es ist wiklich die Umlaufzahl. Wenn man es über mehrere Kurven macht 
gehen auch Löcher und so :-))