对角扁平矩阵的邻域索引计算

我有一个二维矩阵存储在沿对角线的平面缓冲区中。例如，一个4x4矩阵的索引应该像这样分散:

 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

我们假设这里永远不会发生溢出，也没有绕行。

我有一个涉及大量三角数和三角根计算的解决方案。它还包含许多分支，如果可能的话，我更愿意用代数替换(这将在分散控制流昂贵的gpu上运行)。我的解决方案是有效的，但很长。我觉得一定有一种更简单，计算量更少的方法。

如果有人能给这个特殊的问题/表述一个名字，也许会对我有所帮助。

我可以张贴我的完整解决方案，如果有人感兴趣，但正如我所说，这是非常长，相对复杂的这样一个简单的任务。简而言之，我的解决方案是:

将原始索引转换为更大的三角形矩阵，以避免处理2个三角形(例如13将变成17)

对于4x4矩阵，这将是:

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为"下三角形")，并在内部设置一个循环，通过在上三角形中添加1和在下三角形中减去1来累积偏移量(如果这有意义的话)。但对于我的解决方案，必须不惜一切代价避免循环，特别是具有不平衡行程计数的循环(再次，非常对gpu不利)。

所以我正在寻找更多的代数解决方案，而不是算法解决方案。

创建查找表:

同样，由于GPU的原因，最好避免构建查找表并在其中进行随机访问(非常昂贵)。代数解法更可取。

矩阵:

矩阵的大小是已知的。
现在我只考虑方阵，但矩形的解决方案也是不错的。
正如我的例子中函数的名称所暗示的那样，将解决方案扩展到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("%dn", getNGonalNeighbor(15,-1,-1));
    printf("%dn", getNGonalNeighbor(15, 0,-1));
    printf("%dn", getNGonalNeighbor(15,-1, 0));
    printf("%dn", getNGonalNeighbor(11,-2,-1));
    printf("%dn", getNGonalNeighbor(0, -1,-1));
    return 0;
}

不要使用table version.

#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("%dn", getNGonalNeighbor(15,-1,-1));
    printf("%dn", getNGonalNeighbor(15, 0,-1));
    printf("%dn", getNGonalNeighbor(15,-1, 0));
    printf("%dn", getNGonalNeighbor(11,-2,-1));
    printf("%dn", getNGonalNeighbor(0, -1,-1));
    return 0;
}

我实际上已经有元素来解决它在我的代码的其他地方。正如BLUEPIXY的解决方案所暗示的那样，我正在使用分散/聚集操作，我已经为布局转换实现了这些操作。

该解决方案基本上重建矩阵中给定元素的原始(x,y)索引，应用索引偏移量并将结果转换回转换后的布局。它将正方形分成两个三角形，并根据它属于哪个三角形来调整计算。

它几乎完全是一个代数变换:它不使用循环和表查找，内存占用小，分支少。代码可能还可以进一步优化。

下面是代码的草稿:

#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("%dn", getNGonalNeighbor(15,-1,-1)); // should return 11
    printf("%dn", getNGonalNeighbor(15, 0,-1)); // should return 14
    printf("%dn", getNGonalNeighbor(15,-1, 0)); // should return 13
    printf("%dn", getNGonalNeighbor(11,-2,-1)); // should return 1
}

输出确实是:

如果你认为这个解决方案看起来过于复杂和低效，我提醒你这里的目标是GPU，与内存访问相比，计算成本几乎为零，所有索引计算都是使用大规模并行架构同时计算的。

相关内容

最新更新

热门标签：