obliczenia
Rysowanie dużych wypełnianych elips
Artykuł archiwalny z magazynu dyskowego „Izviestia” (1998). Opisuje szybki algorytm rysowania dużych wypełnianych elips, przekraczających możliwości blittera. Przykładowy kod jest przeznaczony dla komputerów Amiga i napisany w asemblerze procesorów M68000.
Grzegorz Kraszewski

Systematyczni czytelnicy „Izviestii” pamiętają być może mój artykuł na temat rysowania wypełnianych figur w poprzednim numerze. Wspomniałem tam chyba o wszystkich funkcjach graphics.library dotyczących bezpośrednio rysowania wypełnianych figur. Cóż więc jeszcze można w tym temacie wymyślić? Okazuje się, że dużo. Jak wiadomo wypełniane elipsy rysujemy funkcją AreaEllipse(). Promienie elipsy podajemy w rejestrach D2 i D3 procesora jako liczby 16-bitowe. Można by więc przypuszczać, że maksymalny promień elipsy to 32 767 pikseli (promienie są w inkludach podane jako liczby WORD, a więc ze znakiem). Tak jednak nie jest z dwóch przyczyn. Obie wynikają z zasady działania funkcji AreaEllipse(). Przed jej wywołaniem musimy do RastPortu dołączyć TmpRas, czyli obszar pamięci, na którym procesor narysuje brzeg elipsy, a blitter ją zamaluje. Szkopuł jednak w tym, że na TmpRas musi zmieścić się cała elipsa, niezależnie od tego, jaki jej kawałek widać na ekranie. Po narysowaniu na TmpRasie elipsa jest kopiowana do RastPortu i jednocześnie przycinana do jego granic. Rozmiar TmpRasa można orientacyjnie określić jako promień_x × promień_y / 2. Jeżeli rysujemy np. elipsę na cały ekran w HiResLaced (640×512), to TmpRas zajmie drobne 40 kB. Ale jeżeli pokusimy się o narysowanie koła (które jest, jakby nie patrzeć, szczególnym przypadkiem elipsy) o promieniu 1000 pikseli, to nagle potrzeba 489 kB jakże cennego chip-ramu... Wreszcie elipsa o obu promieniach równych 2048 pikseli potrzebuje równiutkie 2 MB chipu, a nie widziałem póki co Amigi która miała by więcej. Druga sprawa to rozmiar okna blittera, który jest też skończony i z tego, co wiem ogranicza nas do elips o promieniach do 512 pikseli na blitterze OCS/ECS i zdaje się, że AGA nie wnosi nic nowego w tym temacie.

No dobrze, ale po co nam właściwie takie wielkie elipsy? Przecież nie wlezie to na ekran! Ale wyobraźmy sobie jakiś program do grafiki wektorowej. Narysowaliśmy małą elipsę i powiększamy obraz. W końcu wypada nam narysować na ekranie jakiś kawałek elipsy o promieniu 5000 pikseli. I co? Można oczywiście wyliczyć brzeg tego kawałka, narysować go i zamalować resztę. Wobec tego drugi przykład: trzeba wydrukować tą elipsę na drukarce w skali 1:1. Jeżeli na monitorze (którego „naturalna” rozdzielczość w trybie HighRes Laced i ekranie 14 cali wynosi około 75 dpi) elipsa ma skromny promień 200 pikseli, to na drukarce w 360 dpi mamy już 960. A przecież chip-ram nie jest cały nasz, może użytkownik chce sobie w czasie pracy posłuchać modułu, albo otworzyć kilka ekranów innych programów. Dlatego warto, by duże elipsy rysować nie zużywając za wiele pamięci chip. Okazuje się, że można je rysować przy zerowym zapotrzebowaniu na chip-ram!

Najpierw zajmijmy się problemem wyznaczenia brzegu elipsy. Pierwsze podejście to zastosowanie funkcji trygonometrycznych. Jak wiadomo dla elipsy mamy:

x = a⋅cos(t)
y = b⋅sin(t)

gdzie kąt t zmienia się od 0 do 360 stopni. Dla danego x można znaleźć y w sposób następujący:

y(x) = ±b⋅sin(arccos(x/a))

Niestety funkcje trygonometryczne, to nie jest dobry sposób na przyspieszenie programu. Co prawda koprocesor rulez, ale nie każdy go ma. No to może inaczej: równanie kanoniczne elipsy to

(x/a)² + (y/b)² = 1, a stąd
y = ±b⋅sqrt(1 - (x/a)²)

To już lepiej, ale pierwiastek dalej oznacza koproca, albo szybkość Windowsów na 386 :-). Jeżeli ktoś wpadł na pomysł stablicowania sinusa, czy pierwiastka, to radzę o tym zapomnieć, chyba, że chce przeznaczyć kilkaset kB na tablicę. To niestety nie demo i dokładność musi być znacznie większa. Na szczęście jest jeszcze jeden algorytm, zwany decyzyjnym. Wymaga on wyłącznie dodawania, odejmowania i mnożenia. Jego wersja dla kół była swego czasu opisana w „Amigowcu”. Wykorzystujemy tu, że ekran komputera nie jest czymś ciągłym, ale składa się z pikseli. Jeżeli więc znajdujemy się w jakimś punkcie na brzegu elipsy, to mamy osiem pikseli na które możemy się przenieść, wystarczy podjąć decyzję (stąd nazwa algorytmu), który z nich jest najlepszy. A jak się podejmuje decyzję? Punkt, który idealnie leży na brzegu elipsy spełnia równanie kanoniczne napisane powyżej. Można je przekształcić do postaci:

x²b² + y²a² − a²b² = 0

Ponieważ elipsa jest krzywa, a piksele prostokątne, to nie będą one dokładnie spełniać równania i lewa strona będzie różna od zera. Wystarczy teraz policzyć odchyłki od zera dla wszystkich ośmiu pikseli, a następnie wybrać ten, dla którego wartość bezwzględna odchyłki jest najmniejsza. Można sobie jeszcze sprawę uprościć zauważając, że elipsa to figura symetryczna i wystarczy jedynie wykreślić jej ćwiartkę, reszta to odbicia symetryczne względem osi x i y. Będąc na danej ćwiartce do dyspozycji mamy już tylko trzy kierunki ruchu (nie będziemy przecież wracać po już narysowanych pikselach) i odchyłkę liczymy tylko trzy razy. Gdybyśmy rysowali tylko koła, możnaby wykorzystać symetrię koła względem prostych pod kątem 45 stopni i podzielić koło na 8 części. Wtedy w każdej części obowiązują już tylko 2 kierunki. Wróćmy jednak do elips. Pozostaje pytanie jak znaleźć ten pierwszy punkt? Nic prostszego, wystarczy wybrać jeden z czterech punktów leżących na osiach elipsy. Oto więc sposób, w jaki możemy narysować brzeg elipsy, przy założeniu, że wyliczamy prawą dolną ćwiartkę, a pozostałe rysujemy symetrycznie.

  1. Startujemy z punktu (0, b).
  2. Obliczamy odchyłki dla punktu (x + 1, y), (x, y − 1) i (x + 1, y − 1).
  3. Wybieramy punkt o najmniejszej odchyłce i stawiamy tam kropkę.
  4. Stawiamy kropki w punktach symetrycznych, tzn. (−x, y), (x, −y) i (−x, −y).
  5. Jeżeli y > 0, skocz do punktu 2.

Aby elipsa była zamalowana, a nie pusta w środku zamiast stawiać piksele rysujemy linie od punktu (−x, y) do (x, y) i od (−x, −y) do (x, −y). Warto zauważyć, że linie trzeba rysować tylko wtedy, gdy zmienia się współrzędna y. Oto procedura w języku C:

void AreaBigEllipse (struct RastPort *rp,
  WORD x, WORD y, WORD a, WORD b)
{
  WORD dx, dy;
  float ab2, err1, err2, err2;

  /* wyrażenie stałe - liczę tylko raz */

  ab2 = a * a * b * b;

  /* punkt początkowy */

  dx = 0;
  dy = b;

  while (dy > 0)
  {
    err1 = (dx + 1) * (dx + 1) * b * b;
    err2 = (dy - 1) * (dy - 1) * a * a;
    err3 = err1 + dy * dy * a * a - ab2;
    err1 += err2 - ab2;
    err2 += dx * dx * b * b - ab2;

    /* wartości bezwzględne */

    err1 = SPAbs(err1);
    err2 = SPAbs(err2);
    err3 = SPAbs(err3);

    if ((err3 <= err2) && (err3 <= err1)) dx++;
    else
    {
      dy--;
      if (err2 > err1) dx++;
      RectFill(rp, x - dx, y + dy, x + dx, y + dy);
      RectFill(rp, x - dx, y - dy, x + dx, y - dy);
    }
  }
}

Jak widać użyłem tu liczb zmiennoprzecinkowych, gdyż iloczyn czterech liczb szesnastobitowych jest poza zakresem liczb całkowitych. Spowalnia to nieco działanie procedury. Zamiast rysowania linii użyłem RectFill(), który jest szybszy z dwóch powodów. Po pierwsze amigowy blitter szybciej rysuje prostokąty niż linie, a po drugie w RectFill() są krótsze niż w Draw() obliczenia przycinające linię do granic RastPortu, no i odpada wywołanie Move(). Niestety to jeszcze nie koniec problemów, bowiem procedura ta ma ukrytą wadę. Przy odpowiednio dużych elipsach ich brzegi stają się nierówne i poszarpane, a potem wycinane w sinusoidalne fale. Wynika to z niedostatecznej precyzji liczb float. Można zastosować liczby double, ale to jeszcze bardziej spowolni procedurę. Można też użyć arytmetyki całkowitej na liczbach 64-bitowych. Zaletą jest 100% dokładność i uniezależnienie się od koprocesora. Wady? Trzeba się niestety przeprosić z asemblerem. Po kilkugodzinnej walce powstała ostateczna wersja procedury:

;Procedura AreaBigEllipse().
;WEJŚCIE:a1 - adres RastPortu
;        d0 - x
;        d1 - y
;        d2 - a
;        d3 - b   (tak samo jak w systemowej
;                 AreaEllipse)

    INCLUDE "utility/utility_lib.i"
    INCLUDE "graphics/graphics_lib.i"

    XREF    _UtilityBase
    XREF    _GfxBase

    XDEF    _AreaBigEllipse

_AreaBigEllipse:
    MOVEM.L d0-d7/a0/a2-a6,-(sp)
    LINK    a5,#-40
    MOVE.W  d0,d4             ;x
    MOVE.W  d1,d5             ;y
    CLR.W   d6                ;dx = 0
    MOVE.W  d3,d7             ;dy = b
    MULS    d2,d2             ;a²
    MULS    d3,d3             ;b²
    MOVE.L  d2,d0
    MOVE.L  d3,d1
    MOVEA.L _UtilityBase,a6
    JSR     _LVOSMult64(a6)   ;a²b²
    MOVE.L  d1,-40(a5)
    MOVE.L  d0,-36(a5)

MainLoop:
    TST.W d7                  ;dy = 0?
    BEQ.W   JumpOut    
    MOVE.W  d7,d0             ;dy
    MULS    d0,d0             ;dy²
    MOVE.L  d2,d1             ;a²
    MOVEA.L _UtilityBase,a6
    JSR     _LVOSMult64(a6)   ;dy²a²
    MOVE.L  d1,-32(a5)
    MOVE.L  d0,-28(a5)

    MOVE.W  d6,d0             ;dx
    MULS    d0,d0             ;dx²
    MOVE.L  d3,d1             ;b²
    JSR     _LVOSMult64(a6)   ;dx²b²
    MOVE.L  d1,-24(a5)
    MOVE.L  d0,-20(a5)

    MOVE.W  d7,d0             ;dy
    SUBQ.W  #1,d0             ;dy - 1
    MULS    d0,d0             ;(dy - 1)²
    MOVE.L  d2,d1             ;a²
    JSR     _LVOSMult64(a6)   ;(dy - 1)²a²
    MOVE.L  d1,-16(a5)
    MOVE.L  d0,-12(a5)

    MOVE.W  d6,d0             ;dx
    ADDQ.W  #1,d0             ;dx + 1
    MULS    d0,d0             ;(dx + 1)²
    MOVE.L  d3,d1             ;b²
    JSR     _LVOSMult64(a6)   ;(dx + 1)²b²
    MOVE.L  d1,-8(a5)
    MOVE.L  d0,-4(a5)

    LEA     -28(a5),a0        ;adres dy²a²
    LEA     -4(a5),a2         ;adres (dx + 1)²b²
    MOVE.L  (a2),d0
    ADD.L   d0,(a0)
    ADDX.L  -(a2),-(a0)       ;dy²a² + (dx + 1)²b²
    LEA     -28(a5),a0
    LEA     -36(a5),a2        ;adres a²b²
    MOVE.L  (a2),d0
    SUB.L   d0,(a0)           ;dy²a² + (dx + 1)²b²
    SUBX.L  -(a2),-(a0)       ;- a²b² (błąd dla
    BPL.S   Err3Positive      ;kierunku prawo)
    NEG.L   4(a0)
    NEGX.L  (a0)         

Err3Positive:
    LEA     -20(a5),a0        ;adres dx²b²
    LEA     -12(a5),a2        ;adres (dy - 1)²a²
    MOVE.L  (a2),d0
    ADD.L   d0,(a0)
    ADDX.L  -(a2),-(a0)       ;(dy - 1)²a² + dx²b²
    LEA     -20(a5),a0
    LEA     -36(a5),a2        ;adres a²b²
    MOVE.L  (a2),d0
    SUB.L   d0,(a0)           ;(dy - 1)²a² + dx²b²
    SUBX.L  -(a2),-(a0)       ;- a²b² (błąd dla
    BPL.S   Err2Positive      ;kierunku góra)
    NEG.L   4(a0)
    NEGX.L  (a0)

Err2Positive:
    LEA     -12(a5),a0        ;adres (dy - 1)²a²
    LEA     -4(a5),a2         ;adres (dx + 1)²b²
    MOVE.L  (a2),d0
    ADD.L   d0,(a0)
    ADDX.L  -(a2),-(a0)       ;(dy - 1)²a² + (dx + 1)²b²
    LEA     -12(a5),a0
    LEA     -36(a5),a2        ;adres a²b²
    MOVE.L  (a2),d0
    SUB.L   d0,(a0)           ;(dy - 1)²a² + (dx + 1)²b²
    SUBX.L  -(a2),-(a0)       ;- a²b² (błąd dla
    BPL.S   Err1Positive      ;kierunku prawo-góra)
    NEG.L   4(a0)
    NEGX.L  (a0)

Err1Positive:
    MOVE.L  -16(a5),d0        ;blok porównań (który
    CMP.L   -24(a5),d0        ;błąd jest mniejszy)
    BGT.S   Err1GrtErr2
    BLT.S   Err2GrtErr1
    MOVE.L  -12(a5),d0
    CMP.L   -20(a5),d0
    BLE.S   Err2GrtErr1

Err1GrtErr2:
    MOVE.L  -24(a5),d0
    CMP.L   -32(a5),d0
    BGT.S   Err3Best
    BLT.S   Err2Best
    MOVE.L  -20(a5),d0
    CMP.L   -28(a5),d0
    BGE.S   Err3Best
    BRA.S   Err2Best

Err2GrtErr1:
    MOVE.L  -16(a5),d0
    CMP.L   -32(a5),d0
    BGT.S   Err3Best
    BLT.S   Err1Best
    MOVE.L  -12(a5),d0
    CMP.L   -28(a5),d0
    BLE.S   Err1Best
     
Err3Best:
    ADDQ.W  #1,d6             ;dx++
    BRA.W   MainLoop

Err1Best:
    ADDQ.W #1,d6              ;dx++

Err2Best:
    SUBQ.W #1,d7              ;dy--
    MOVEA.L d2,a3             ;przechowaj a² i b²
    MOVEA.L d3,a4
    MOVE.L  a1,-(sp)          ;przechowaj adres RPortu
    MOVEA.L _GfxBase,a6       ;w a1 jest adres RPortu
    MOVE.W  d4,d0
    SUB.W   d6,d0             ;x - dx
    MOVE.W  d4,d2
    ADD.W   d6,d2             ;x + dx
    MOVE.W  d5,d1
    ADD.W   d7,d1             ;y + dy
    MOVE.W  d1,d3
    JSR     _LVORectFill(a6)
    MOVE.W  d4,d0
    SUB.W   d6,d0             ;x - dx
    MOVE.W  d4,d2
    ADD.W   d6,d2             ;x + dx
    MOVE.W  d5,d1
    SUB.W   d7,d1             ;y - dy
    MOVE.W  d1,d3
    MOVEA.L (sp),a1
    JSR     _LVORectFill(a6)
    MOVE.L  a3,d2             ;odtwórz a² i b²
    MOVE.L  a4,d3
    MOVEA.L (sp)+,a1          ;i adres RastPortu
    BRA.W   MainLoop

JumpOut:
    UNLK    a5
    MOVEM.L (sp)+,d0-d7/a0/a2-a6
    RTS

    END

Procedura została przeze mnie praktycznie wypróbowana, spisuje się bardzo dobrze, jest szybka i krótka (równiutkie 400 bajtów). Teraz słów kilka o arytmetyce 64-bitowej. Dodawanie i odejmowanie wykorzystuje rozkazy rozszerzonej precyzji ADDX i SUBX. Podobnie negowanie liczby potrzebne przy obliczaniu wartości bezwzględnej, wykonane jest przez NEGX. Do mnożenia wykorzystałem procedurkę SMult64() z utility.library. Uniezależnia ona procedurę od procesora, dla 020 i 030 jest to po prostu zwykłe MULS.L d0,d0:d1, dla innych (zdaje się że w 040 i 060 nie ma tych instrukcji, nie mówiąc oczywiście o 68000/010) mały podprogramik. Niestety pociąga to za sobą wymaganie Kickstartu 3.0+. Jeżeli ktoś chce, może do kodu wstawić odpowiednią instrukcję lub procedurę na stałe. Zmienne 64-bitowe są przechowywane w ramce stosu stworzonej rozkazem LINK, ułatwia to użycie na nich rozkazów ADDX i SUBX. Ponieważ w RectFill() blitter kreśli linie bezpośrednio na bitplanach, procedura nie wymaga ani jednego bajtu dodatkowego chip-ramu. Dodatkowo przy elipsach większych niż mniej więcej pół ekranu (DblPAL HighRes NoFilcker 640×512) procedura jest szybsza od systemowej (!!!) i przewaga wzrasta w miarę zwiększania elipsy (sprawdzone na A1200 z 030@50 MHz i 8 MB fast). Procedura oczywiście pozwala na narysowanie elipsy o promieniach do 32767 pikseli, co nawet na naświetlarce poligraficznej o rozdzielczości 2400 dpi da figurę o średnicy 70 cm.

Przy korzystaniu z AreaBigEllipse() warto pamiętać o kilku rzeczach:

Oczywiście kod procedury jest Public Domain.