Кривые Гильберта и Серпинского.

Кривые Гильберта и Серпинского, или снова рекурсия

Д.М.Златопольский

В статье [1] был приведен ряд рекурсивных алгоритмов для построения изображений, составленных из повторяющихся фигур. Ряд аналогичных программ представлен в [2]. C использованием рекурсии можно получить также фигуры, показанные на рис. 1 и 2. Первая из них называется кривая Гильберта, вторая — кривая Серпинского.

Рис. 1

Рис. 2

Эти кривые связаны с любопытным понятием теории функций, а именно — всюду плотными кривыми [7]. Кривая на плоскости называется всюду плотной в некоторой области, если она проходит через любую сколь угодно малую окрестность каждой точки этой области. Несколько упрощенно можно считать, что всюду плотные кривые целиком заполняют указанную область. Известные математики Гильберт и Серпинский построили примеры всюду плотных кривых. Хотя эти примеры различны, схема получения соответствующих кривых одинакова. По определенному правилу строятся кривые (соответственно Гильберта и Серпинского) первого, второго, ..., n-го порядка, вписанные в заданный квадрат. При неограниченном возрастании n они стремятся к некоторой предельной кривой, которая является всюду плотной в заданном квадрате.

В этом номере мы рассмотрим алгоритм построения кривой Гильберта, а в следующем — кривой Серпинского.

Кривая Гильберта первого порядка, обозначаемая H1, похожа на изображение буквы П, вычерченной в виде трех сторон квадрата, как показано на рис. 3а. На рис. 3б изображена кривая Гильберта второго порядка H2. Видно, что кривая H2 состоит из кривых H1, ориентированных в разные стороны (вправо, вверх и влево). Кривые H1, составляющие кривую H2, соединены тремя отрезками прямых, называемых связками (на рис. 3б они вычерчены утолщенными линиями). В действительности эти отрезки должны иметь одинаковую толщину с другими линиями, утолщенными они показаны единственно с целью демонстрации способа получения H2 из Н1.

Рис. 3. Кривые Гильберта
первого, второго и третьего порядков

Аналогично фигуру H3 (рис. 3в) можно рассматривать как состоящую из четырех кривых H2 (ориентированных в разные стороны) и трех связок.

Заметим, что отрезки, образующие линию H1, можно рассматривать как связки, соединяющие 4 точки — кривые Гильберта нулевых порядков.

Таким образом, кривую Гильберта i-го порядка Hi можно получить из четырех кривых Hi–1, ориентированных в разные стороны, и трех связок. Если процедуры рисования кривых Hi, ориентированных вверх, вниз, влево и вправо, обозначить соответственно GU(i), GD(i), GL(i) и GR(i), то можно составить следующие рекурсивные схемы построения этих кривых:

GU(i): GR(i—1) GU(i—1) GU(i—1) GL(i—1)

GR(i): GU(i—1) GR(i—1) GR(i—1) GD(i—1)

GD(i): GL(i—1) GD(i—1)¬ GD(i—1) GR(i—1)

GL(i): GD(i—1) GL(i—1) GL(i—1) GU(i—1)

Обозначим через h длину элементарного отрезка прямых в кривых H_i. Тогда, например, процедура GU в программе на школьном алгоритмическом языке
[3, 4] может быть оформлена следующим образом:

алг GU(арг цел i)
нач
   если i>0
   то
    GR(i-1)
    LineUp
   GU(i-1)
    LineRight
    GU(i-1)
    LineDown
    GL(i-1)
все
кон

Здесь LineUp, LineRight, LineLeft, LineDown — процедуры рисования связок, направленных соответственно вверх, вправо, влево и вниз (напомним, что ось Y на экране направлена сверху вниз).

алг LineUp
нач
| вектор(0, -d)
кон

алг LineRight
нач
| вектор(d, 0)
кон

алг LineLeft
нач
вектор(-h, 0)
кон

алг LineDown
нач
| вектор(0, d)
кон

Аналогично можно оформить процедуры рисования кривых Гильберта, ориентированных вниз, влево и вправо:

алг GD(арг цел i)
нач
    если i>0
     то
       GL(i-1)
      LineDown
      GD(i-1)
      LineLeft
     GD(i-1)
      LineUp
      GR(i-1)
   все
кон

алг GL(арг цел i)
нач
   если i>0
    то
      GD(i-1)
     LineLeft
     GL(i-1)
     LineDown
    GL(i-1)
     LineRight
    GU(i-1)
все
кон

алг GR(арг цел i)
нач
    если i>0
     то
       GU(i-1)
      LineRight
     GR(i-1)
      LineUp
      GR(i-1)
      LineLeft
      GD(i-1)
     все
кон

Обратим внимание на то, что в приведенных процедурах рисования кривых Гильберта используется так называемая косвенная рекурсия — ситуация, когда процедура вызывает себе как вспомогательную не только непосредственно, а также и через другую процедуру [1].

Квадрат, в который вписывается кривая Гильберта, будем называть опорным, длину его стороны (в пикселях) обозначим через S. Обсудим теперь вопрос определения значения величины h в зависимости от порядка кривой n. Из рис. 3 видно, что при n=2 длина элементарного отрезка линии в три раза меньше стороны опорного квадрата, при n=3 — в семь раз. Отсюда получаем, что коэффициенты уменьшения для этих элементарных отрезков в фигурах H₁, H₂ , H₃, ... образуют ряд чисел 1, 3, 7, ..., то есть в общем случае коэффициент уменьшения для фигуры H_n может быть вычислен по формуле 2n—1.

Естественно расположить изображаемую кривую Гильберта по центру экрана. Для этого надо найти координаты x₀, y₀ начальной точки кривой. Проанализировав приведенные выше процедуры, можно убедиться, что при ориентации кривой вверх и влево она начинается с левой нижней точки опорного квадрата, т.е.

x₀=X_c — S/2; y₀ =Y_c+S/2,

а в остальных случаях — с правой верхней точки опорного квадрата, т.е. в этих случаях

x₀=X_c + S/2; y₀ =Y_c -- S/2,

Здесь Xc, Yc — координаты центра экрана.

Кроме того, удобно задавать размер опорного квадрата в процентах от высоты экрана, поскольку она всегда меньше ширины. Эту величину обозначим PrS. В программе построения кривой Гильберта используем, помимо указанных обозначений, еще переменную orient — число, определяющее ориентацию кривой (вверх — 1, вниз — 2, вправо — 3, влево — 4).

Школьный алгоритмический язык

цел h |Величину h описываем как глобальную
алг Кривая Гильберта
нач цел n, x0, y0, s, orient, Hscr, Wscr
     вещ PrS
     Hscr:=480 |Высота экрана
     Wscr:=640 |Ширина экрана
     |Вводим исходные данные для построения кривой Гильберта
    нц
       вывод нс, "Введите длину стороны опорного квадрата"
       вывод "в % от высоты экрана"
       ввод PrS
    кц при PrS<100
    вывод нс, "Введите порядок кривой"
    ввод n
    нц
      вывод нс, "Введите ориентацию кривой"
      вывод "(вверх - 1, вниз - 2,
      вправо - 3, влево - 4)"
      ввод orient
    кц при (orient>=1) и (orient<=4)
    S:=Int(PrS/100*Hscr)
      |Сторона опорного квадрата
    h:=div(S, 2**n-1) | Длина связок
     |Находим координаты начальной точки кривой
    если (orient=1) или (orient=3)
      то
        x0:=Div(Wscr,2)- Div(S,2)
        y0:=Div(Hscr,2)+ Div(S,2)
     иначе
        x0:=Div(Wscr,2)+ Div(S,2)
        y0:=Div(Hscr,2)- Div(S,2)
     все
   |Устанавливаем графический режим работы экрана
   видео(17) | VGA экран 640*480
   поз(x0, y0) |Начальная точка кривой
    |Рисуем соответствующий вариант кривой Гильберта
   выбор
     при orient=1: GU(n)
     при orient=2: GD(n)
     при orient=3: GR(n)
     иначе GL(n)
   все
   |Переходим в текстовый режим
   видео(0)
кон

Чего не хватает в этой программе? Хотелось бы наблюдать последовательность построения кривой, а для этого в конец процедур GU, GD, GR, GL необходимо включить процедуру задержки. Для школьного алгоритмического языка мы оставим указанные процедуры без изменений, а реализуем эту идею в программах на других языках программирования. Отметим, что параметр процедуры задержки надо подбирать экспериментально, поскольку его величина зависит от быстродействия компьютера и порядка кривой Гильберта.

Язык Паскаль

При реализации алгоритма построения кривой Гильберта на Паскале возникает проблема, связанная с тем, что все четыре процедуры рисования кривых Гильберта (направленных в разные стороны) используют друг друга в качестве вспомогательных. Это не дает возможности соблюсти правило, согласно которому каждый идентификатор перед употреблением должен быть описан. Выходом является так называемое опережающее описание процедур, обращение к которым фигурирует раньше их описания [4]. Такими процедурами являются GD и GU.

Uses crt, graph;
const
del=5000;{Время задержки}
PATH=' ';
{Файлы *.BGI находятся в рабочем каталоге}
Var
d,r :integer;
n, orient :byte;
x0, y0, S, h,Hscr,Wscr : word;
PrS : real;
{Процедуры рисования связок. От последней точки (на нее указывает графический курсор) проводится вниз, вверх, влево, вправо отрезок длиной h пикселей. Напомним, что ось Y графического экрана направлена сверху вниз}
Procedure LineDown; begin Linerel(0, h) end;
Procedure LineUp; begin Linerel(0, -h) end;
Procedure LineLeft; begin Linerel(-h, 0) end;
Procedure LineRight; begin Linerel(h, 0) end;
{Опережающее описание процедур GD и GU, вызываемых до своего определения}
Procedure GD(i: byte); forward;
Procedure GU(i: byte); forward;
{Процедуры рисования четырех разновидностей кривых Гильберта}
Procedure GL (i: byte);
begin
if i > 0 then begin
     GD(i-1); LineLeft;
     GL(i-1); LineDown;
     GL(i-1); LineRight;
     GU(i-1); Delay(del);
   end
end;

Procedure GR(i: byte);
begin
    if i>0 then begin
    GU(i-1); LineRight;
    GR(i-1); LineUp;
    GR(i-1); LineLeft;
    GD(i-1); Delay(del);
   end
end;

Procedure GU;
{Параметр i процедуры GU указан при опережающем описании}
begin
if i>0 then begin
     GR(i-1); LineUp;
     GU(i-1); LineRight;
     GU(i-1); LineDown;
     GL(i-1); Delay(del);
   end
end;

Procedure GD;
{Параметр i процедуры GD указан при опережающем описании}
begin
    if i>0 then begin
       GL(i-1); LineDown;
       GD(i-1); LineLeft;
       GD(i-1); LineUp;
       GR(i-1); Delay(del);
    end
end;

Function Power2(n: byte): word;
{Возведение 2 в степень n}
var p,i: word;
begin
    p:=2; for i:=1 to n-1 do p:=p*2;
    Power2:=p
end;
BEGIN
   clrscr;{Чистка экрана}
   {Вводим исходные данные для построения кривой Гильберта}
   repeat
       write('Введите длину стороны опорного квадрата');
      write(' в % от высоты экрана ');
      readln(PrS);
   until PrS<100;
   write('Введите порядок кривой ');
   readln(n);
   repeat
     write('Введите ориентацию кривой ');
     write('вверх - 1, вниз - 2, вправо - 3, влево - 4 ');
     readln(orient);
   until (orient>=1) and (orient<=4);
   d:=detect;
   initgraph(d, r, PATH);
     {Пеpеход в гpафический pежим}
   Hscr:=GetMaxY+1;{Высота экрана}
   Wscr:=GetMaxX+1;{Ширина экрана}
   S:=round(PrS/100*Hscr); {Сторона квадрата}
   h:=round(S/(Power2(n)-1)); {длина связок}
    {Находим координаты начальной точки кривой.
      Для ориентации: вверх и вправо начальная точка - левая нижняя точка квадрата;
       для ориентации: вниз и влево - правая верхняя точка квадрата}
   Case orient of
      1,3:{ориентация: вверх или вправо}
        begin
            x0:=Wscr div 2 - S div 2;
            y0:=Hscr div 2 + S div 2;
       end;
      2,4:{ориентация: вниз или влево}
       begin
           x0:=Wscr div 2 + S div 2;
           y0:=Hscr div 2 - S div 2;
       end;
   end;{Case orient}
     {Графический курсор устанавливаем в начальную точку}
   moveto(x0, y0);
    {Рисуем соответствующий вариант кривой Гильберта}
    case orient of
       1: GU(n); 2: GD(n); 3: GR(n); 4: GL(n);
    end {case orient};
       readln; {Выход - нажатием клавиши Enter}
    closegraph {Переходим в текстовый режим}
END.

Язык Бейсик (вариант QuickBasic [5])

DECLARE SUB Delay (del&) 'Задеpжка
'Процедуры рисования четырех разновидностей 'кривых Гильберта
DECLARE SUB GL (i AS INTEGER)
DECLARE SUB GR (i AS INTEGER)
DECLARE SUB GU (i AS INTEGER)
DECLARE SUB GD (i AS INTEGER)
'Процедуры рисования связок
DECLARE SUB LineDown ()
DECLARE SUB LineUp ()
DECLARE SUB LineLeft ()
DECLARE SUB LineRight ()
'Описания переменных
DIM PrS AS SINGLE, S AS SINGLE, x0 AS SINGLE, y0 AS SINGLE
DIM n AS INTEGER, orient AS INTEGER
DIM SHARED h AS SINGLE
'Константы
DIM SHARED del&: del& = 200000
'паpаметp задеpжки
Hscr! = 480: Wscr! = 640
'Высота и шиpина экpана
'Основной алгоритм
CLS 'Чистка экрана
'Вводим исходные данные для построения 'кривой Гильберта
DO
   PRINT "Введите длину стороны опорного
   квадрата";
   INPUT " в % от высоты экрана ", PrS
LOOP UNTIL PrS < 100
INPUT "Введите порядок кривой ", n
DO
   PRINT "Введите ориентацию кривой: ";
   INPUT "вверх - 1, вниз - 2, вправо - 3,
    влево - 4 ", orient
LOOP UNTIL orient >= 1 AND orient <= 4
S = PrS / 100! * Hscr! 'Сторона квадрата
h = S / (2 ^ n - 1) 'длина связок
'Находим координаты начальной точки кривой.
'Для ориентации: вверх и вправо начальная 'точка - левая нижняя точка квадрата;
'для ориентации: вниз и влево -
'правая верхняя точка квадрата
IF orient = 1 OR orient = 3 THEN
         x0 = Wscr! / 2 - S / 2
         y0 = Hscr! / 2 + S / 2
     ELSE
         x0 = Wscr! / 2 + S / 2
         y0 = Hscr! / 2 - S / 2
END IF
'Пеpеход в гpафический pежим для монитора 'VGA. Экpан 640*480
SCREEN (12)
'Графический курсор устанавливаем в 'начальную точку
PSET (x0, y0)
'Рисуем соответствующий вариант кривой 'Гильберта
SELECT CASE orient
       CASE 1: CALL GU(n)
       CASE 2: CALL GD(n)
       CASE 3: CALL GR(n)
       CASE 4: CALL GL(n)
END SELECT
DO: LOOP UNTIL INKEY$ = CHR$(13)
'Выход - нажатием клавиши Enter
SCREEN 0 'Переходим в текстовый режим
END
SUB Delay (del&)'Процедура задержки
      DIM i&
      FOR i& = 1 TO del&: NEXT i&
END SUB
SUB GD (i AS INTEGER)
     IF i > 0 THEN
         CALL GL(i - 1): CALL LineDown
         CALL GD(i - 1): CALL LineLeft
         CALL GD(i - 1): CALL LineUp
         CALL GR(i - 1): CALL Delay(del&)
     END IF
END SUB
SUB GL (i AS INTEGER)
   IF i > 0 THEN
        CALL GD(i - 1): CALL LineLeft
        CALL GL(i - 1): CALL LineDown
        CALL GL(i - 1): CALL LineRight
        CALL GU(i - 1): CALL Delay(del&)
   END IF
END SUB
SUB GR (i AS INTEGER)
    IF i > 0 THEN
        CALL GU(i - 1): CALL LineRight
        CALL GR(i - 1): CALL LineUp
        CALL GR(i - 1): CALL LineLeft
        CALL GD(i - 1): CALL Delay(del&)
    END IF
END SUB
SUB GU (i AS INTEGER)
     IF i > 0 THEN
         CALL GR(i - 1): CALL LineUp
         CALL GU(i - 1): CALL LineRight
         CALL GU(i - 1): CALL LineDown
         CALL GL(i - 1): CALL Delay(del&)
     END IF
END SUB
SUB LineDown
      LINE -STEP(0, h)
END SUB
SUB LineLeft
      LINE -STEP(-h, 0)
END SUB
SUB LineRight
      LINE -STEP(h, 0)
END SUB
SUB LineUp
     LINE -STEP(0, -h)
END SUB

Язык Си

В языке Си для функций, обращение к которым предшествует их определению, надо указать прототип [7]. Такими функциями являются GD и GU.
#include <stdio.h>
#include <graphics.h>
#include <math.h>
#include <conio.h>
#include <dos.h>
#include<stdlib.h>
#define PATH ""
/* файлы *.BGI находятся в pабочем каталоге */
#define del 5000 /* Время задержки */
int h;
/* Прототипы функций GD и GU, вызываемых
до своего определения */
void GD(int); void GU(int);
/* Округление вещественного числа до ближайшего целого */
int round(float a)
{return (int)floor(a+0.5);}
/*Функции рисования связок. От последней точки (на нее указывает графический курсор)
проводится вниз, вверх, влево, вправо отрезок длиной h пикселей. Напомним, что ось Y графического экрана направлена сверху вниз */
void LineDown() {linerel(0, h);}
void LineUp() {linerel(0, -h);}
void LineLeft() {linerel(-h,0);}
void LineRight() {linerel(h,0);}
/*Функции рисования четырех разновидностей кривых Гильберта*/
void GL(int i)
{if (i>0)
          {GD(i-1); LineLeft();
            GL(i-1); LineDown();
            GL(i-1); LineRight();
            GU(i-1); delay(del);
           }
}
void GR(int i)
{if (i>0)
           {GU(i-1); LineRight();
             GR(i-1); LineUp();
             GR(i-1); LineLeft();
             GD(i-1); delay(del);
            }
}
void GU(int i)
{
if (i>0)
         {GR(i-1); LineUp();
           GU(i-1); LineRight();
           GU(i-1); LineDown();
           GL(i-1); delay(del);
          }
}
void GD(int i)
{if (i>0)
         {GL(i-1); LineDown();
           GD(i-1); LineLeft();
           GD(i-1); LineUp();
           GR(i-1); delay(del);
          }
}
void main()
{int d=DETECT,r,n,orient,x0, y0, S,Hscr,Wscr;
float PrS;
clrscr();/*Чистка экрана*/
/*Вводим исходные данные для построения кривой Гильберта*/
do
          {printf("\nВведите длину стороны опорного квадрата");
            printf(" в % от высоты экрана ");
            scanf("%f",&PrS);
          }
while (PrS>=100);
printf("\nВведите порядок кривой "); scanf("%d",&n);
do
           {printf("\nВведите ориентацию кривой ");
             printf("вверх - 1, вниз - 2, вправо - 3,
             влево - 4 ");
             scanf("%d",&orient);
           }
while (orient<1 || orient>4);
initgraph(&d, &r, "");
/* Пеpеход в гpафический pежим */
Hscr=getmaxy()+1;/*Высота экрана*/
Wscr=getmaxx()+1;/*Ширина экрана*/
S=round(PrS/100*Hscr); /*Сторона квадрата*/
h=round(S/(pow(2,n)-1));/*длина связок*/
/*Находим координаты начальной точки кривой. Для ориентации: вверх и вправо начальная точка - левая нижняя точка квадрата; для ориентации: вниз и влево -
правая верхняя точка квадрата */
if (orient == 1 || orient == 3)
           {x0=Wscr/2 - S/2; y0=Hscr/2 + S/2;}
       else {x0=Wscr/2 + S/2; y0=Hscr/2 - S/2;}
/*Графический курсор устанавливаем в начальную точку*/
moveto(x0, y0);
/*Рисуем соответствующий вариант кривой Гильберта*/
switch (orient)
           {case 1: GU(n); break;
             case 2: GD(n); break;
             case 3: GR(n); break;
             case 4: GL(n); break;
           }
getch(); /*Выход - нажатием любой клавиши*/
closegraph();
/*Переходим в текстовый режим*/
}

ЛИТЕРАТУРА

1. Златопольский Д.М. Рекурсия. Информатика, 1996, № 7, 8.
2. Островский С.Л., Гольдшлаг О.Я. Фрактальные кривые. Информатика, 1995, № 23.
3. Кушниренко А.Г., Лебедев Г.В., Сворень Р.А. Основы информатики и вычислительной техники. М.: Просвещение, 1990.
4. Эпиктетов М.Г. Почему школьный алгоритмический? Информатика, 1995, № 24.
5. Фаронов В.В. Программирование на персональных ЭВМ в среде ТУРБО-ПАСКАЛЬ. М.: Изд-во МГТУ, 1992.
6. Зельднер Г.А. QuickBasic 4.5. М.: ABF, 1994.
7. Александров П.С. Введение в общую теорию множеств и функций. М.: Гостехиздат, 1948.
8. Романовская Л.М. и др. Программирование в среде Си. М.: Финансы и статистика, 1992.

Окончание в следующем номере