快速查找图像中最近的非黑色像素

Question

我有一个二维图像，像素随机且稀疏地分散。
给定图像上的一个点，我需要找到到不属于背景颜色（黑色）的最近像素的距离。
最快的方法是什么？

我能想到的唯一方法是为像素构建一个 kd 树。 但我真的想避免如此昂贵的预处理。 另外，kd-tree 似乎给了我比我需要的更多的东西。 我只需要到某个东西的距离，而不关心这个东西是什么。

Answer 1

正如 Pyro 所说，搜索一个正方形的周长，每次从原始点移出一个像素（即一次将宽度和高度增加两个像素）。 当您击中非黑色像素时，您会计算距离（这是您的第一个昂贵的计算），然后继续向外搜索，直到框的宽度是到第一个找到点的距离的两倍（超出此范围的任何点都不可能更近）比您原来找到的像素）。 保存您在这部分中找到的所有非黑点，然后计算它们的每个距离，看看它们中是否有任何一个比您的原始点更近。

在理想的查找中，您只需进行一次昂贵的距离计算。

更新：因为您在这里计算像素到像素的距离（而不是任意精度的浮点位置），所以您可以通过使用预先计算的查找表（只是高度-by-width 数组），以 x 和 y 的函数形式给出距离。 一个 100x100 的数组实际上需要 40K 的内存，并且覆盖了原始点周围 200x200 的正方形，并且节省了您为找到的每个彩色像素进行昂贵的距离计算（无论是毕达哥拉斯还是矩阵代数）的成本。 该数组甚至可以预先计算并作为资源嵌入到您的应用程序中，以节省您的初始计算时间（这可能是严重的过度杀伤力）。

更新 2：此外，还有一些方法可以优化正方形周长的搜索。 您的搜索应该从与轴相交的四个点开始，并向角落一次移动一个像素（您有 8 个移动搜索点，这很容易使这变得比它值得的麻烦，具体取决于您的应用程序的要求）。 一旦找到彩色像素，就无需继续向角移​​动，因为其余点都离原点更远。

在第一个找到的像素之后，您可以使用查找表进一步将所需的额外搜索区域限制为最小，以确保每个搜索点比找到的点更近（再次从轴开始，并在达到距离限制时停止） ）。 如果您必须动态计算每个距离，那么第二次优化的成本可能会太高。

如果最近的像素位于 200x200 框内（或任何适合您的数据的尺寸），您将仅在由像素界定的圆圈内进行搜索，仅进行查找和比较。

Answer 2

就我个人而言，我会忽略 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
   }
}

Tada！ r^2 计算仅进行位移位、加法和减法:)

当然，在任何像样的现代 CPU 上计算 r^2 = dx*dx + dy*dy 可能和这个一样快......

Answer 3

好吧，这听起来很有趣。
我做了一个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;
}

Answer 4

是的，最近邻居搜索很好，但不能保证您会找到“最近的”。 每次移出一个像素将产生方形搜索 - 对角线将比水平/垂直更远。 如果这很重要，您需要验证 - 继续扩展，直到绝对水平距离大于“找到的”像素，然后计算找到的所有非黑色像素的距离。

Answer 5

您没有指定如何测量距离。 我假设 L1（直线），因为它更容易； 也许这些想法可以针对 L2（欧几里德）进行修改。

如果您只对相对较少的像素执行此操作，则只需从源像素开始螺旋式向外搜索，直到找到非黑色像素。

如果您要对其中的许多/全部执行此操作，那么这样如何：构建一个图像大小的二维数组，其中每个单元格存储到最近的非黑色像素的距离（如果需要，还存储该像素的坐标） ）。 进行四行扫描：从左到右、从右到左、从下到上、从上到下。 考虑从左到右的扫描； 扫描时，保留包含每行中最后一个非黑色像素的一维列，并使用到该像素的距离和/或坐标标记二维数组中的每个单元。 O(n^2)。

另外，kd 树也太过分了； 你可以使用四叉树。 只是比我的行扫描更难编码一点，内存多一点（但不到两倍），而且可能更快。

Answer 6

我研究过并且可能会坚持使用的另一种方法：利用 Bresenham 圆算法。
它的速度快得惊人，因为它可以帮助您节省任何距离比较！
您实际上只是在目标点周围绘制越来越大的圆圈，这样当您第一次遇到非黑色像素时，您会自动知道它是最近的，从而节省任何进一步的检查。
我尚未验证的是布雷森汉姆圆是否会捕获每个像素，但这对我的情况来说不是问题，因为我的像素将出现在某种类型的斑点中。

Answer 7

搜索“最近邻居搜索”，Google 中的前两个链接应该可以帮助您。

如果您只对每个图像 1 个像素执行此操作，我认为您最好的选择就是线性搜索，一次向外 1 个像素宽度的框。 如果您的搜索框是方形的，您将无法获取找到的第一个点。 你必须要小心

Answer 8

我确信这可以做得更好，但这里有一些代码可以搜索中心像素周围的正方形的周长，首先检查中心，然后向角落移动。 如果未找到像素，则周长（半径）将扩展，直到达到半径限制或找到像素。 第一个实现是一个围绕中心点进行简单螺旋的循环，但如前所述，它找不到绝对最近的像素。 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;

}

Answer 9

您可以结合多种方法来加快速度。

• 加速像素查找的一种方法是使用我所说的空间查找图。 它基本上是该块中像素的下采样图（例如 8x8 像素，但这是一种权衡）。 值可以是“无像素集”、“部分像素集”、“所有像素集”。 这样，一次读取就可以判断块/单元是否已满、部分满或空。
• 扫描中心周围的盒子/矩形可能并不理想，因为有许多距离很远的像素/单元。 我使用圆形绘制算法（Bresenham）来减少开销。
• 读取原始像素值可以水平批量进行，例如一个字节（对于 8x8 或其倍数的单元大小）、双字或长字节。 这应该会让你再次获得明显的加速。
• 您还可以使用多个级别的“空间查找地图”，这又是一个权衡。

对于距离计算，可以使用提到的查找表，但它是（缓存）带宽与计算速度的权衡（例如，我不知道它在 GPU 上的执行情况）。

Answer 10

我会做一个简单的查找表 - 对于每个像素，预先计算到最近的非黑色像素的距离，并将该值存储在与相应像素相同的偏移量中。 当然，这样你就需要更多的内存。