Krzywa HILBERTa

Kod źródłowy: mixyh.c

Zapełnianie kwadratu krzywą łamaną

Najlepiej zobaczyć pierwsze kroki rekurencyjnej konstrukcji. Celem jest konstrukcja nierekurencyjnego algorytmu pozwalającego przekształcać numer załamania na krzywej na parę liczb określających jego położenie w dwóch wymiarach.

0 == lbn  
    \x 0
    y+--+
    0| 0|
     +--+

Przypadek 0 == lbn jest trywialny i nie wnosi wiele.

1 == lbn
    \x 0 1
    y+-----+
    0| 0-1 |
     |   | |
    1| 3-2 |
     +-----+

Przypadek 1 == lbn prezentuje podstawowy schemat pokrywania powierzchni krzywą. Możemy z niego między innymi zobaczyć, że fragment numer 3 krzywej znajmuje miejsce o  współrzędnych (x,y)[3]==(0,1).

2 == lbn
    \x  0  1  2  3
    y+-------------+
    0|  0  3--4--5 |
     |  |  |     | |
    1|  1--2  7--6 |
     |        |    |
    2| 14-13  8--9 |
     |  |  |     | |
    3| 15 12-11-10 |
     +-------------+

Przypadek 2 == lbn jest już bardzo pouczający. Widać z niego w jaki sposób następuje wypełnianie przestrzeni dwuwymiarowej przy zwiększeniu rozdzielczości. Dwukrotne zwiększenie rozdzielczości wymaga czterokrotnego zwiększenia długości krzywej pokrywającej powierzchnię kwadratu.

Weźmy jakikolwiek fragment krzywej, na przykład numer 12. Jego współrzędne wynoszą: (x,y)[12]=(1,3). Zapiszmy to używając dwójkowej reprezentacji liczb: (x,y)[1100]=(01,11).

Przygądnijmy się jeszcze jednemu, bardziej złożonemu przypadkowi, aby pojąć występujące to regularności:

3 == lbn
    \x  0  1  2  3  4  5  6  7
    y+-------------------------+
    0|  0--1 14-15-16 19-20-21 |
     |     |  |     |  |     | |
    1|  3--2 13-12 17-18 23-22 |
     |  |        |       |     |
    2|  4  7--8 11 30-29 24-25 |
     |  |  |  |  |  |  |     | |
    3|  5--6  9-10 31 28-27-26 |
     |              |          |
    4| 58-57 54-53 32 35-36-37 |
     |  |  |  |  |  |  |     | |
    5| 59 56-55 52 33-34 39-38 |
     |  |        |        |    |
    6| 60-61 50-51 46-45 40-41 |
     |     |  |     |  |     | |
    7| 63-62 49-48-47 44-43-42 |
     +-------------------------+

Tym razem weźmy fragment 49. Jego wsółrzędne wynoszą: (x,y)[49]==(2,7), w reprezentacji dwójkowej jest to: (x,y)[110001]==(010,111).

Fragment w większej rozdzielczościach był wybierany w tej samej ćwiartce (następnie szestnastce). Numer takiego fragmenty dziedziczy bardziej znaczące bity dwójkowej reprezentacji. Tylko dwa najmniej znaczące bity określają położenie fragmenty w  obszarze najdrobniejszej rozdzielczosci. To samo odnosi się do współrzednych fragmentu, tu tylko jeden najmniej znaczący bit każdej ze wsółrzędnych nie jest dziedziczony z mniejszej rozdzielczości.

Zatem, aby obliczyć wsółrzędne zadanego fragmentu krzywej o długości 2^(2*lbn), trzeba jego numer w reprezentacji dwójkowej podzielić na lbn dwubitowych części. Każda dwubitowa część opisuje położenie w ćwiartce na odpowiednim poziomie rozdzielczości.

Zadanie to realizuje ten fragment kodu:

	while(i < lbn) bt[i++] = (int)count & 0x03, count >>=2;

Dochodzimy tu do najtrudniejszej i najistotniejszej części algorytmu. Znaczwenie dwubitowej części (czyli położenie w ćwiartce) zalezy od położenia w mniejszej skali (czyli wartości bardziej znaczących części dwubitowych). Ta zależność jest zakodowana w dwóch tablicach. Tablica hilb opisuje cztery orientację schematu pokrywania.

	static short hilb[4][4] = {
		{ 0, 1, 3, 2 },
		{ 0, 2, 3, 1 },
		{ 3, 1, 0, 2 },
		{ 3, 2, 0, 1 }
	};

To znaczy:

	hilb[0]     hilb[1]     hilb[2]     hilb[3]

	0[0]--1[1]  0[0]  3[1]  2[0]--1[1]  2[0]--3[1] 
	      |     |     |     |     |     |
	3[2]--2[3]  1[2]--2[3]  3[2]  0[3]  1[2]--0[3]

Tablica turn opisuje zmianę orientacji ćwiartki względem orientacji zawierającej ją ćwiartki:

	static short turn[4] = { 1, 0, 0, 2 };

Liczby tablicy turn są indeksami wierszy tablicy hilb. Mówią, że przy 2==lbn orentacje komórki podstawowej są takie, jak opisują { turn[1], turn[0], turn[0], turn[2] }. Odtwarzanie orientacji realizuje ta linia kodu:

		short bb = bt[i];
		bb = bb[hilb[kk ^= turn[bb]]];

Wreszcie rozseparowanie bitów do współrzędnych:

		x = (x << 1) | (bb & 1);
		y = (y << 1) | (bb >> 1);

---

(C) Aleksander Nabagło