Кривая Гильберта vs Z-order

[Stochkas]

очень интересно. спасибо.

[DimPal]

Если интересно могу предложить альтернативу кривой Гильберта — собственный велосипед. Алгоритм оптимизирован для случайной выборки.

[zzeng]

Конечно, интересно. Оптимизирован по какому критерию?

[DimPal]

Для 2D случая кривая выглядит как то так:

private const	d	: int = 4;
private const	dm	: int = (1 << d) - 1;
private var		xx_ : int;
private var		yy_ : int;
public function get_GCurve_XY(i: int) : void
{
	var j	: int;
	var v	: int;
	var b	: int;
	var f	: int;
	xx_ = yy_ = b = 0;
	f = i;
	for (j = d - 1; j >= 0; --j)
	{
		f += (1 << (j << 1)) >> 1;
		v = ((f >> (j << 1)) + b) & 3;
		b = v ^ 2;
		v  ^= b >> 1;
		xx_ |= (v &  1) << j;
		yy_ |= (v >> 1) << j;
	}
}
//.............................................................................
public function DP_Curve() : void
{
	var px	: int;
	var py	: int;
	var x	: int;
	var y	: int;
	var i	: int;
	var nn	: int = (1 << (d << 1));
	for (i = 0; i < nn; ++i)
	{
		get_GCurve_XY(i);
		if (0 != i)
			draw_Line(px, py, xx_, yy_);
		px = xx_;
		py = yy_;
	}
}

[DimPal]

Этот алгоритм можно перенести на N-мерный случай.
Кстати, для Z-order функция выглядит так:

public function get_ZCurve_XY(i: int) : void
{
	var v : int;
	v = (i&0x99999999) | ((i&0x44444444) >> 1) | ((i << 1)&0x44444444);
	v = (v&0xc3c3c3c3) | ((v&0x30303030) >> 2) | ((v << 2)&0x30303030);
	v = (v&0xF00FF00F) | ((v&0x0F000F00) >> 4) | ((v << 4)&0x0F000F00);
	v = (v&0xFF0000FF) | ((v&0x00FF0000) >> 4) | ((v << 4)&0x00FF0000);
	xx_ = v & 0xFFFF;
	yy_ = v >> 16;
}

[DimPal]

Извините, ошибочка… Вот так правильно:

public function get_ZCurve_XY(i: int) : void
{
	var v : int;
	v = (i&0x99999999) | ((i&0x44444444) >> 1) | ((i << 1)&0x44444444);
	v = (v&0xc3c3c3c3) | ((v&0x30303030) >> 2) | ((v << 2)&0x30303030);
	v = (v&0xF00FF00F) | ((v&0x0F000F00) >> 4) | ((v << 4)&0x0F000F00);
	v = (v&0xFF0000FF) | ((v&0x00FF0000) >> 8) | ((v << 8)&0x00FF0000);
	xx_ = v & 0xFFFF;
	yy_ = v >> 16;
}

[zzeng]

Симпатично. Проблема с этой кривой в том, что ее невозможно одним разрезом разделить на два непрерывных интервала значений. А значит работать с ней крайне неудобно алгоритмически. Т.е. можно разделить по диагонали, но на самой диагонали окажутся точки с обеих половин.

[DimPal]

Можно конечно и по диагонали. А можно и по окто-дереву. Кривая замыкается в кольцо и начало можно назначать произвольно. Если взять начало в середине получатся те же свойства что и у кривой Гилберта. Кстати, функцию get_GCurve_XY можно еще примерно в два раза ускорить.

[zzeng]

Первый раз можно по середине. А потом всё равно по диагонали придётся.
Либо резать на 4 части. В двумерном случае, на 256 частей в 8-мерном.

[DimPal]

Да, понял. На больших размерах четверть квадов побито диагоналями. Ну тогда вот вам утешительный приз:

private function get_Gilbert_XY(i: int) : void
{
	var j	: int;
	var x	: int;
	var y	: int;
	var d	: int;
	var f	: int;
	xx_ = yy_ = f = 0;
	for (j = curve_d_ - 1; j >= 0; --j)
	{
		d = ((i >> (j << 1)) & 3) | f;
		f = f ^ ((0xa0004 >> ((d & 3) * 5)) & 0x1c);
		xx_ |= ((0x3c93c66c & (1 << d)) >> d) << j;
		yy_ |= ((0x6c3993c6 & (1 << d)) >> d) << j;
	}
}

[DimPal]

Даже еще проще:

private function get_Gilbert_XY(i: int) : void
{
	var j	: int;
	var d	: int;
	var f	: int;
	i <<= 2;
	xx_ = yy_ = f = 0;
	for (j = d - 1; j >= 0; --j)
	{
		d = ((i >> (j << 1)) & 12) | f;
		f ^= (12289 >> (d & 12)) & 3;
		xx_ |= ((37740 >> d) & 1) << j;
		yy_ |= ((25500 >> d) & 1) << j;
	}
}

Причем если использовать массив заранее рассчитанных значений (до любой выбранной глубины) то количество итераций в цикле можно серьезно уменьшить.

[zzeng]

Хорошо, когда всё вмещается в 32/64 разряда.

[DimPal]

А сколько разрядов используется в ваших задачах?

[zzeng]

32 * число размерностей

[DimPal]

Ну возьмем для примера 3D.
Значит в наивной реализации 32 прохода по циклу. А в оптимизированном варианте можно сделать, скажем, за 4 прохода, правда нужен индекс порядка 64Mb размером.

[zzeng]

Я использовал алгоритм A. R. Butz, «Alternative Algorithm for Hilbert's Space-Filling Curve», IEEE Trans. Comp., April, 1971, pp 424-426.
Насколько понимаю, это компромисс между числом проходов и объемом предвычисленных данных.

[DimPal]

До меня как до жирафа доходило что у вас многомерный случай, и используемая библиотека оптимизирована именно под это. Уверен что для 2D/3D/4D моя функция «уделает» hilbert_i2c() в тестах (не хотите проверить?), но:
— под бОльшие мерности пространства оптимизации меняются
— мой трехмерный мозг туго работает с N-мерностями выше трех (это усложняет визуальную отладку)

[DimPal]

И на счет быстродействия пару замечаний:
— на фоне скорости чтения/записи памяти математические операции для процессора это раз плюнуть. Даже sin/cos вычисляются за почти за пару тактов.
— время доступа для кэшей процессора register/L1/L2/… для каждого следующего уровня отличается на порядок: 1 такт/10 тактов/100 тактов/1000 тактов. Поэтому не факт что использование предвычисленных данных будет давать значительный прирост, возможно даже и наоборот (но это не точно, практика — критерий истины).

[DimPal]

Пояснения к примененным методам для функции построения кривой Гилберта:

1. Минимальный конструктивный элемент я обозначил для себя термином «симплекс».
   Размер симплекса 2 в степени количества измерений. Например для 3D это кубик размером 2x2x2. Эмпирически мной установлено, что элементы симплекса обходятся по коду Грея, например для 3D каждому из трех бит соответствует координатная ось.
2. Симплексы подвергаются рекурсивным трансформациям, чтобы конец предыдущего и начало следующего составили один шаг. Эмпирически мной установлено, что из всех возможных трансформаций используется только 2 в степени измерений (для 3D случая 8 из 16 возможных). Причем если визуально упростить симплексы соединив отрезком только начало и конец, то множество возможных трансформаций заполняют кривую построенную по коду грея. Четные/нечетные отрезки будут направлены в разные стороны.
3. Рекурсивно трансформируемый элемент обходится по коду Грея. Эмпирически мной установлено, что индекс трансформации тоже можно получить по коду Грея, но тогда кривая замкнется. Что бы этого не происходило первый и последний элемент нужно развернуть «рогами» наружу
4. Ну а дальше банальная бинарная арифметика в помощь. После оптимизаций получилась функция get_Gilbert_XY(i: int)

На практике может понадобиться обойти элементы кривой Гилберта по очереди. Можно решить задачу «в лоб» как в приведенном мной примере или рекурсивными вызовами, но существуют и более оптимальные способы это сделать, правда это совсем уже другая история…

[zzeng]

Спасибо