【发布时间】：2015-05-08 17:55:56
【问题描述】：
我有一个 2D 图像随机且稀疏地散布着像素。
给定图像上的一个点，我需要找到到不在背景颜色（黑色）中的最近像素的距离。
最快的方法是什么？
我能想到的唯一方法是为像素构建一个 kd 树。但我真的很想避免如此昂贵的预处理。此外，似乎 kd-tree 给了我比我需要的更多。我只需要与某物的距离，我不在乎这东西是什么。
【问题讨论】：
标签： performance graphics 2d pixel
【发布时间】：2015-05-08 17:55:56
【问题描述】：
就个人而言，我会忽略 MusiGenesis 关于查找表的建议。
计算像素之间的距离并不很昂贵，特别是对于这个初始测试，您不需要实际距离，因此无需取平方根。您可以使用距离^2，即：
r^2 = dx^2 + dy^2
此外，如果您一次向外移动一个像素，请记住：
(n + 1)^2 = n^2 + 2n + 1
或者如果 nx 是当前值，而 ox 是之前的值：
nx^2 = ox^2 + 2ox + 1
= ox^2 + 2(nx - 1) + 1
= ox^2 + 2nx - 1
=> nx^2 += 2nx - 1
很容易看出这是如何工作的：
1^2 = 0 + 2*1 - 1 = 1
2^2 = 1 + 2*2 - 1 = 4
3^2 = 4 + 2*3 - 1 = 9
4^2 = 9 + 2*4 - 1 = 16
5^2 = 16 + 2*5 - 1 = 25
etc...
因此，在每次迭代中，您因此只需要保留一些中间变量：
int dx2 = 0, dy2, r2;
for (dx = 1; dx < w; ++dx) { // ignoring bounds checks
dx2 += (dx << 1) - 1;
dy2 = 0;
for (dy = 1; dy < h; ++dy) {
dy2 += (dy << 1) - 1;
r2 = dx2 + dy2;
// do tests here
}
}
多田！ r^2 计算只有位移、加法和减法:)
当然，在任何体面的现代 CPU 上计算 r^2 = dx*dx + dy*dy 可能和这个一样快...
【讨论】：
正如 Pyro 所说，搜索一个正方形的周长，从原始点一次移动一个像素（即，一次将宽度和高度增加两个像素）。当你击中一个非黑色像素时，你计算距离（这是你的第一个昂贵的计算），然后继续向外搜索，直到你的框的宽度是到第一个找到的点的距离的两倍（超出这个的任何点都不可能更近比你原来找到的像素）。保存您在此部分中找到的所有非黑点，然后计算它们的距离，看看它们是否比您的原始点更近。
在理想的查找中，您只需进行一次昂贵的距离计算。
更新：因为您在此处计算像素到像素的距离（而不是任意精度的浮点位置），您可以通过使用预先计算的查找表（只是一个高宽数组），以将距离作为 x 和 y 的函数。一个 100x100 的数组基本上要花费 40K 的内存并覆盖原始点周围的 200x200 方格，并且为您节省了为找到的每个彩色像素进行昂贵的距离计算（无论是毕达哥拉斯还是矩阵代数）的成本。这个数组甚至可以预先计算并作为资源嵌入到您的应用中，以节省您的初始计算时间（这可能是严重的矫枉过正）。
更新 2：此外，还有一些方法可以优化搜索正方形周长。您的搜索应该从与轴相交的四个点开始，并且一次向角落移动一个像素（您有 8 个移动的搜索点，这很容易使这比它的价值更麻烦，这取决于您的应用程序的要求）。一旦您找到一个彩色像素，就无需继续向角落移动，因为其余点都离原点更远。
在找到第一个像素之后，您可以通过使用查找表进一步将所需的额外搜索区域限制为最小，以确保每个搜索点都比找到的点更近（再次从轴开始，当距离已达到上限）。如果您必须即时计算每个距离，那么第二次优化可能会过于昂贵。
如果最近的像素在 200x200 框内（或任何适合您的数据的大小），您将只在该像素限定的圆圈内搜索，只进行查找和 比较。
【讨论】：
-
我们真的需要一直到两倍的盒子大小来确保没有更近的点吗？我认为达到 ceil(d * sqrt(2)) 的盒子大小就足够了。因为对于大小为 d 的盒子，就距离而言，最远点是 d * sqrt(2)。因此，任何比这更大的框都不需要检查。
您没有指定测量距离的方式。我将假设 L1（直线），因为它更容易；可能这些想法可以针对 L2（欧几里得）进行修改。
如果您只对相对较少的像素进行此操作，则只需从源像素以螺旋状向外搜索，直到找到非黑色像素。
如果您要为许多/全部执行此操作，那么如何：构建一个与图像大小相同的二维数组，其中每个单元存储到最近的非黑色像素的距离（以及必要时的坐标该像素的）。进行四行扫描：从左到右、从右到左、从下到上和从上到下。考虑从左到右扫过；扫描时，保留一维列，其中包含每行中看到的最后一个非黑色像素，并用到该像素的距离和/或坐标标记二维数组中的每个单元格。 O(n^2)。
另外，k-d 树是多余的；你可以使用四叉树。只比我的线扫描更难编码，更多的内存（但不到两倍），并且可能更快。
【讨论】：
-
我看不出扫描算法是如何正确的。如果最近的像素位于某个对角线方向，则此方法将找不到它。它只会找到从该点开始的两个轴上的像素。
搜索“最近邻搜索”，Google 中的前两个链接应该对您有所帮助。
如果您只为每张图像 1 个像素执行此操作，我认为您最好的选择只是线性搜索，每次向外 1 个像素宽度的框。如果您的搜索框是方形的，则您无法获取找到的第一个点。一定要小心
【讨论】：
是的，最近邻搜索很好，但不能保证您会找到“最近”。每次移出一个像素将产生方形搜索 - 对角线将比水平/垂直更远。如果这很重要，您需要验证 - 继续扩展，直到绝对水平距离大于“找到”像素，然后计算所有已定位的非黑色像素的距离。
【讨论】：
好的，这听起来很有趣。 我做了一个 c++ 版本的解决方案，我不知道这是否对你有帮助。我认为它工作得足够快，因为它在 800*600 矩阵上几乎是即时的。如果您有任何问题，请随时提出。
抱歉我犯了任何错误，这是一个 10 分钟的代码... 这是一个迭代版本（我也打算制作一个递归版本，但我改变了主意）。 该算法可以通过不向点数组添加任何点来改进，该点距起点的距离大于 min_dist，但这涉及计算每个像素（尽管它是颜色）到起点的距离。
希望有帮助
//(c++ version)
#include<iostream>
#include<cmath>
#include<ctime>
using namespace std;
//ITERATIVE VERSION
//picture witdh&height
#define width 800
#define height 600
//indexex
int i,j;
//initial point coordinates
int x,y;
//variables to work with the array
int p,u;
//minimum dist
double min_dist=2000000000;
//array for memorising the points added
struct point{
int x;
int y;
} points[width*height];
double dist;
bool viz[width][height];
// direction vectors, used for adding adjacent points in the "points" array.
int dx[8]={1,1,0,-1,-1,-1,0,1};
int dy[8]={0,1,1,1,0,-1,-1,-1};
int k,nX,nY;
//we will generate an image with white&black pixels (0&1)
bool image[width-1][height-1];
int main(){
srand(time(0));
//generate the random pic
for(i=1;i<=width-1;i++)
for(j=1;j<=height-1;j++)
if(rand()%10001<=9999) //9999/10000 chances of generating a black pixel
image[i][j]=0;
else image[i][j]=1;
//random coordinates for starting x&y
x=rand()%width;
y=rand()%height;
p=1;u=1;
points[1].x=x;
points[1].y=y;
while(p<=u){
for(k=0;k<=7;k++){
nX=points[p].x+dx[k];
nY=points[p].y+dy[k];
//nX&nY are the coordinates for the next point
//if we haven't added the point yet
//also check if the point is valid
if(nX>0&&nY>0&&nX<width&&nY<height)
if(viz[nX][nY] == 0 ){
//mark it as added
viz[nX][nY]=1;
//add it in the array
u++;
points[u].x=nX;
points[u].y=nY;
//if it's not black
if(image[nX][nY]!=0){
//calculate the distance
dist=(x-nX)*(x-nX) + (y-nY)*(y-nY);
dist=sqrt(dist);
//if the dist is shorter than the minimum, we save it
if(dist<min_dist)
min_dist=dist;
//you could save the coordinates of the point that has
//the minimum distance too, like sX=nX;, sY=nY;
}
}
}
p++;
}
cout<<"Minimum dist:"<<min_dist<<"\n";
return 0;
}
【讨论】：
我确信这可以做得更好，但这里有一些代码搜索围绕中心像素的正方形的周长，首先检查中心并移向角落。如果未找到像素，则扩展周长（半径），直到达到半径限制或找到像素。第一个实现是一个围绕中心点进行简单螺旋的循环，但如前所述，它找不到绝对最近的像素。 SomeBigObjCStruct 在循环内的创建非常缓慢 - 将其从循环中移除使其足够好，并且使用了螺旋方法。但无论如何，这是这个实现 - 请注意，几乎没有进行任何测试。
这一切都是通过整数加减法完成的。
- (SomeBigObjCStruct *)nearestWalkablePoint:(SomeBigObjCStruct)point {
typedef struct _testPoint { // using the IYMapPoint object here is very slow
int x;
int y;
} testPoint;
// see if the point supplied is walkable
testPoint centre;
centre.x = point.x;
centre.y = point.y;
NSMutableData *map = [self getWalkingMapDataForLevelId:point.levelId];
// check point for walkable (case radius = 0)
if(testThePoint(centre.x, centre.y, map) != 0) // bullseye
return point;
// radius is the distance from the location of point. A square is checked on each iteration, radius units from point.
// The point with y=0 or x=0 distance is checked first, i.e. the centre of the side of the square. A cursor variable
// is used to move along the side of the square looking for a walkable point. This proceeds until a walkable point
// is found or the side is exhausted. Sides are checked until radius is exhausted at which point the search fails.
int radius = 1;
BOOL leftWithinMap = YES, rightWithinMap = YES, upWithinMap = YES, downWithinMap = YES;
testPoint leftCentre, upCentre, rightCentre, downCentre;
testPoint leftUp, leftDown, rightUp, rightDown;
testPoint upLeft, upRight, downLeft, downRight;
leftCentre = rightCentre = upCentre = downCentre = centre;
int foundX = -1;
int foundY = -1;
while(radius < 1000) {
// radius increases. move centres outward
if(leftWithinMap == YES) {
leftCentre.x -= 1; // move left
if(leftCentre.x < 0) {
leftWithinMap = NO;
}
}
if(rightWithinMap == YES) {
rightCentre.x += 1; // move right
if(!(rightCentre.x < kIYMapWidth)) {
rightWithinMap = NO;
}
}
if(upWithinMap == YES) {
upCentre.y -= 1; // move up
if(upCentre.y < 0) {
upWithinMap = NO;
}
}
if(downWithinMap == YES) {
downCentre.y += 1; // move down
if(!(downCentre.y < kIYMapHeight)) {
downWithinMap = NO;
}
}
// set up cursor values for checking along the sides of the square
leftUp = leftDown = leftCentre;
leftUp.y -= 1;
leftDown.y += 1;
rightUp = rightDown = rightCentre;
rightUp.y -= 1;
rightDown.y += 1;
upRight = upLeft = upCentre;
upRight.x += 1;
upLeft.x -= 1;
downRight = downLeft = downCentre;
downRight.x += 1;
downLeft.x -= 1;
// check centres
if(testThePoint(leftCentre.x, leftCentre.y, map) != 0) {
foundX = leftCentre.x;
foundY = leftCentre.y;
break;
}
if(testThePoint(rightCentre.x, rightCentre.y, map) != 0) {
foundX = rightCentre.x;
foundY = rightCentre.y;
break;
}
if(testThePoint(upCentre.x, upCentre.y, map) != 0) {
foundX = upCentre.x;
foundY = upCentre.y;
break;
}
if(testThePoint(downCentre.x, downCentre.y, map) != 0) {
foundX = downCentre.x;
foundY = downCentre.y;
break;
}
int i;
for(i = 0; i < radius; i++) {
if(leftWithinMap == YES) {
// LEFT Side - stop short of top/bottom rows because up/down horizontal cursors check that line
// if cursor position is within map
if(i < radius - 1) {
if(leftUp.y > 0) {
// check it
if(testThePoint(leftUp.x, leftUp.y, map) != 0) {
foundX = leftUp.x;
foundY = leftUp.y;
break;
}
leftUp.y -= 1; // moving up
}
if(leftDown.y < kIYMapHeight) {
// check it
if(testThePoint(leftDown.x, leftDown.y, map) != 0) {
foundX = leftDown.x;
foundY = leftDown.y;
break;
}
leftDown.y += 1; // moving down
}
}
}
if(rightWithinMap == YES) {
// RIGHT Side
if(i < radius - 1) {
if(rightUp.y > 0) {
if(testThePoint(rightUp.x, rightUp.y, map) != 0) {
foundX = rightUp.x;
foundY = rightUp.y;
break;
}
rightUp.y -= 1; // moving up
}
if(rightDown.y < kIYMapHeight) {
if(testThePoint(rightDown.x, rightDown.y, map) != 0) {
foundX = rightDown.x;
foundY = rightDown.y;
break;
}
rightDown.y += 1; // moving down
}
}
}
if(upWithinMap == YES) {
// UP Side
if(upRight.x < kIYMapWidth) {
if(testThePoint(upRight.x, upRight.y, map) != 0) {
foundX = upRight.x;
foundY = upRight.y;
break;
}
upRight.x += 1; // moving right
}
if(upLeft.x > 0) {
if(testThePoint(upLeft.x, upLeft.y, map) != 0) {
foundX = upLeft.x;
foundY = upLeft.y;
break;
}
upLeft.y -= 1; // moving left
}
}
if(downWithinMap == YES) {
// DOWN Side
if(downRight.x < kIYMapWidth) {
if(testThePoint(downRight.x, downRight.y, map) != 0) {
foundX = downRight.x;
foundY = downRight.y;
break;
}
downRight.x += 1; // moving right
}
if(downLeft.x > 0) {
if(testThePoint(upLeft.x, upLeft.y, map) != 0) {
foundX = downLeft.x;
foundY = downLeft.y;
break;
}
downLeft.y -= 1; // moving left
}
}
}
if(foundX != -1 && foundY != -1) {
break;
}
radius++;
}
// build the return object
if(foundX != -1 && foundY != -1) {
SomeBigObjCStruct *foundPoint = [SomeBigObjCStruct mapPointWithX:foundX Y:foundY levelId:point.levelId];
foundPoint.z = [self zWithLevelId:point.levelId];
return foundPoint;
}
return nil;
}
【讨论】：
您可以结合多种方式来加快速度。
- 一种加速像素查找的方法是使用我所说的空间查找图。它基本上是该块中像素的下采样图（例如 8x8 像素，但它是一种折衷）。值可以是“无像素设置”“部分像素设置”“所有像素设置”。这样一读就可以判断一个块/单元格是满的、部分满的还是空的。
- 围绕中心扫描一个框/矩形可能并不理想，因为有许多像素/单元距离很远。我使用圆形绘制算法 (Bresenham) 来减少开销。
- 读取原始像素值可以水平批量进行，例如一个字节（对于 8x8 或其倍数的单元格大小）、dword 或 long。这应该会再次给您带来显着的加速。
- 您还可以使用多个级别的“空间查找图”，这又是一种权衡。
对于距离计算，可以使用提到的查找表，但它是（缓存）带宽与计算速度的权衡（例如，我不知道它在 GPU 上的表现如何）。
【讨论】：
我会做一个简单的查找表 - 对于每个像素，预先计算到最近的非黑色像素的距离，并将值存储在与相应像素相同的偏移量中。当然，这样你会需要更多的内存。
【讨论】：
-
重点是避免冗长的处理操作，这些操作是 O(num-of-points)
-
好的，我认为与查找距离相比，预处理会少得多。