Расчет индекса соседей для диагонально сплюснутой матрицы
У меня 2D-матрица хранится в плоском буфере по диагонали. Например, у матрицы 4х4 индексы будут разбросаны примерно так:
0 2 5 9
1 4 8 12
3 7 11 14
6 10 13 15
С этим представлением, каков наиболее эффективный способ вычислить индекс соседнего элемента с учетом исходного индекса и смещения X / Y? Например:
// return the index of a neighbor given an offset
int getNGonalNeighbor(const size_t index,
const int x_offset,
const int y_offset){
//...
}
// for the array above:
getNGonalNeighbor(15,-1,-1); // should return 11
getNGonalNeighbor(15, 0,-1); // should return 14
getNGonalNeighbor(15,-1, 0); // should return 13
getNGonalNeighbor(11,-2,-1); // should return 1
Здесь мы предполагаем, что переполнение никогда не происходит и нет обхода.
У меня есть решение, включающее много треугольное число и вычисления треугольного корня. Он также содержит много ветвей, которые я бы предпочел заменить на алгебру, если это возможно (это будет работать на графических процессорах, где расходящийся поток управления стоит дорого). Мое решение работает, но очень долго. Я чувствую, что должен быть намного более простой и менее сложный способ сделать это.
Может быть, это помогло бы мне, если бы кто-то мог поставить название для этой конкретной проблемы / представления.
Я могу опубликовать свое полное решение, если кому-то интересно, но, как я уже сказал, оно очень длинное и относительно сложное для такой простой задачи. Короче говоря, мое решение делает:
- переведите исходный индекс в большую треугольную матрицу, чтобы не иметь дело с 2 треугольниками (например, 13 станет 17)
Для матрицы 4х4 это будет:
0 2 5 9 14 20 27
1 4 8 13 19 26
3 7 12 18 25
6 11 17 24
10 16 23
15 22
21
- Рассчитайте индекс диагонали соседа в этом представлении, используя манхэттенское расстояние смещения и треугольный корень индекса.
- рассчитать положение соседа по этой диагонали, используя смещение
- перевести обратно в исходное представление, удалив отступы.
По некоторым причинам это самое простое решение, которое я мог придумать.
Редактировать:
имея цикл для накопления смещения:
Я понимаю, что, учитывая свойства номеров треугольников, было бы легче разделить матрицу на два треугольника (назовем 0–9 «верхним треугольником» и 10–15 «нижним треугольником») и получить цикл с тестом внутри накапливать смещение, добавляя единицу в верхнем треугольнике и вычитая единицу в нижнем (если это имеет смысл). Но для моего решения следует избегать циклов любой ценой, особенно циклов с несбалансированным количеством отключений (опять же, очень плохо для графических процессоров).
Так Я больше ищу алгебраическое решение, чем алгоритмическое.
Создание таблицы поиска:
Опять же, из-за графического процессора предпочтительно избегать построения таблицы поиска и иметь случайный доступ к ней (очень дорого). Алгебраическое решение является предпочтительным.
Свойства матрицы:
- Размер матрицы известен.
- Пока я рассматриваю только квадратную матрицу, но было бы неплохо и решение для прямоугольных.
- Как следует из названия функции в моем примере, расширение решения на N-мерные объемы (следовательно, N-гональное уплощение) также будет большим плюсом.
Решение
Таблица поиска
#include <stdio.h>
#define SIZE 16
#define SIDE 4 //sqrt(SIZE)
int table[SIZE];
int rtable[100];// {x,y| x<=99, y<=99 }
void setup(){
int i, x, y, xy, index;//xy = x + y
x=y=xy=0;
for(i=0;i<SIZE;++i){
table[i]= index= x*10 + y;
rtable[x*10+y]=i;
x = x + 1; y = y - 1;//right up
if(y < 0 || x >= SIDE){
++xy;
x = 0;
y = xy;;
while(y>=SIDE){
++x;
--y;
}
}
}
}
int getNGonalNeighbor(int index, int offsetX, int offsetY){
int x,y;
x=table[index] / 10 + offsetX;
y=table[index] % 10 + offsetY;
if(x < 0 || x >= SIDE || y < 0 || y >= SIDE) return -1; //ERROR
return rtable[ x*10+y ];
}
int main() {
int i;
setup();
printf("%d\n", getNGonalNeighbor(15,-1,-1));
printf("%d\n", getNGonalNeighbor(15, 0,-1));
printf("%d\n", getNGonalNeighbor(15,-1, 0));
printf("%d\n", getNGonalNeighbor(11,-2,-1));
printf("%d\n", getNGonalNeighbor(0, -1,-1));
return 0;
}
не используйте настольную версию.
#include <stdio.h>
#define SIZE 16
#define SIDE 4
void num2xy(int index, int *offsetX, int *offsetY){
int i, x, y, xy;//xy = x + y
x=y=xy=0;
for(i=0;i<SIZE;++i){
if(i == index){
*offsetX = x;
*offsetY = y;
return;
}
x = x + 1; y = y - 1;//right up
if(y < 0 || x >= SIDE){
++xy;
x = 0;
y = xy;;
while(y>=SIDE){
++x;
--y;
}
}
}
}
int xy2num(int offsetX, int offsetY){
int i, x, y, xy, index;//xy = x + y
x=y=xy=0;
for(i=0;i<SIZE;++i){
if(offsetX == x && offsetY == y) return i;
x = x + 1; y = y - 1;//right up
if(y < 0 || x >= SIDE){
++xy;
x = 0;
y = xy;;
while(y>=SIDE){
++x;
--y;
}
}
}
return -1;
}
int getNGonalNeighbor(int index, int offsetX, int offsetY){
int x,y;
num2xy(index, &x, &y);
return xy2num(x + offsetX, y + offsetY);
}
int main() {
printf("%d\n", getNGonalNeighbor(15,-1,-1));
printf("%d\n", getNGonalNeighbor(15, 0,-1));
printf("%d\n", getNGonalNeighbor(15,-1, 0));
printf("%d\n", getNGonalNeighbor(11,-2,-1));
printf("%d\n", getNGonalNeighbor(0, -1,-1));
return 0;
}
Другие решения
У меня на самом деле уже были элементы, чтобы решить это где-то еще в моем коде. Как подсказало решение BLUEPIXY, я использую операции разброса / сбора, которые я уже реализовал для преобразования макета.
Это решение в основном восстанавливает оригинал (x,y) Индекс данного элемента в матрице, применяет смещение индекса и переводит результат обратно в преобразованный макет. Он разбивает квадрат на 2 треугольника и корректирует вычисления в зависимости от того, к какому треугольнику он принадлежит.
Это почти полностью алгебраическое преобразование: оно не использует цикл и поиск по таблице, имеет небольшой объем памяти и небольшое разветвление. Код может быть оптимизирован в дальнейшем.
Вот черновик кода:
#include <stdio.h>
#include <math.h>
// size of the matrix
#define SIZE 4
// triangle number of X
#define TRIG(X) (((X) * ((X) + 1)) >> 1)
// triangle root of X
#define TRIROOT(X) ((int)(sqrt(8*(X)+1)-1)>>1);
// return the index of a neighbor given an offset
int getNGonalNeighbor(const size_t index,
const int x_offset,
const int y_offset){
// compute largest upper triangle index
const size_t upper_triangle = TRIG(SIZE);
// position of the actual element of index
unsigned int x = 0,y = 0;
// adjust the index depending of upper/lower triangle.
const size_t adjusted_index = index < upper_triangle ?
index :
SIZE * SIZE - index - 1;
// compute triangular root
const size_t triroot = TRIROOT(adjusted_index);
const size_t trig = TRIG(triroot);
const size_t offset = adjusted_index - trig;
// upper triangle
if(index < upper_triangle){
x = offset;
y = triroot-offset;
}
// lower triangle
else {
x = SIZE - offset - 1;
y = SIZE - (trig + triroot + 1 - adjusted_index);
}
// adjust the offset
x += x_offset;
y += y_offset;
// manhattan distance
const size_t man_dist = x+y;
// calculate index using triangular number
return TRIG(man_dist) +
(man_dist >= SIZE ? x - (man_dist - SIZE + 1) : x) -
(man_dist > SIZE ? 2* TRIG(man_dist - SIZE) : 0);
}
int main(){
printf("%d\n", getNGonalNeighbor(15,-1,-1)); // should return 11
printf("%d\n", getNGonalNeighbor(15, 0,-1)); // should return 14
printf("%d\n", getNGonalNeighbor(15,-1, 0)); // should return 13
printf("%d\n", getNGonalNeighbor(11,-2,-1)); // should return 1
}
И вывод действительно:
11
14
13
1
Если вы думаете, что это решение выглядит сложным и неэффективным, я напомню вам, что целью здесь является GPU, где вычисления практически ничего не стоят по сравнению с обращениями к памяти, и все вычисления индекса вычисляются одновременно с использованием массивно параллельных архитектур.