计算欧氏距离的最快方法
我需要以最快的方式计算两点之间的欧氏距离。 在C.
我的代码是这样的,看起来有点慢:
float distance(int py, int px, int jy, int jx){ return sqrtf((float)((px)*(px)+(py)*(py))); }
提前致谢。
编辑:
对不起,我不清楚。 我最好指定上下文:我正在使用图像,我需要从每个像素到所有其他像素的欧几里德距离。 所以我必须计算很多次。 我不能使用距离的平方。 我将添加更多代码以便更清楚:
for (jy=0; jy<sizeY; jy++) { for (jx=0; jx<sizeX; jx++) { if (jx==px && jy==py) { ; } else{ num+=rfun(imgI[py][px].red-imgI[jy][jx].red)/distance(py, px, jy, jx); den+=RMAX/distance(py, px, jy, jx); } } } float distance(int py, int px, int jy, int jx){ return sqrtf((float)((px-jx)*(px-jx)+(py-jy)*(py-jy))); }
这就是我要做的。 我必须用所有像素(px,py)来做
编辑2:对不起,我不清楚,但我尽力保持这个问题的一般性。 我正在编写一个用算法处理图像的程序。 最大的问题是时间因为我必须真的非常快。 现在我需要优化的是这个函数:`float normC(int py,int px,int color,pixel ** imgI,int sizeY,int sizeX){
int jx, jy; float num=0, den=0; if (color==R) { for (jy=0; jy<sizeY; jy++) { for (jx=0; jx<sizeX; jx++) { if (jx==px && jy==py) { ; } else{ num+=rfun(imgI[py][px].red-imgI[jy][jx].red)/distance(py, px, jy, jx); den+=RMAX/distance(py, px, jy, jx); } } } } else if (color==B){ for (jy=0; jy<sizeY; jy++) { for (jx=0; jx<sizeX; jx++) { if (jx==px && jy==py) { ; } else{ num+=rfun(imgI[py][px].blue-imgI[jy][jx].blue)/distance(py, px, jy, jx); den+=RMAX/distance(py, px, jy, jx); } } } } else if (color==G){ for (jy=0; jy<sizeY; jy++) { for (jx=0; jx<sizeX; jx++) { if (jx==px && jy==py) { ; } else{ num+=rfun(imgI[py][px].green-imgI[jy][jx].green)/distance(py, px, jy, jx); den+=RMAX/distance(py, px, jy, jx); } } } } return num/den;
}
在每个像素(px; py)的循环中调用此函数,因此它被调用了很多次,并且需要花费大量时间来计算它。 rfun
函数不可优化,因为它已经非常简单快速。 我需要做的是加快距离function。
1)我试过hypotf
但它比distance
函数慢
2)我增加了编译器的优化设置,使进程快2倍!
3)我尝试使用宏#define DIST(x, y) sqrtf((float)((x)*(x)+(y)*(y)))
但没有任何改变(如我所料)
EDIT3:
最后,我发现最快的方法是计算所有可能的距离并将它们存储在一个数组中,然后再开始计算normC函数的循环。 为了使它更快,我计算了距离的倒数,以便我可以使用产品而不是商:
float** DIST DIST=malloc(500*sizeof(float*)); //allocating memory for the 2d array for (i=0; i<500; i++) { DIST[i]=malloc(500*sizeof(float)); } for(i=0; i<500; i++){ //storing the inverses of the distances for (p=0; p<500; p++) { DIST[i][p]= 1.0/sqrtf((float)((i)*(i)+(p)*(p))); } } float normC(int py, int px, int color, pixel** imgI, int sizeY, int sizeX){ int jx=0, jy=0; float num=0, den=0; if (color==R) { for (jy=0; jy<sizeY; jy++) { for (jx=0; jx=jy && px>=jx){ num+=rfun(imgI[py][px].red-imgI[jy][jx].red)*DIST[py-jy][px-jx]; den+=DIST[py-jy][px-jx]; } else if (jy>py && px>=jx){ num+=rfun(imgI[py][px].red-imgI[jy][jx].red)*DIST[jy-py][px-jx]; den+=DIST[jy-py][px-jx]; } else if (jy>py && jx>px){ num+=rfun(imgI[py][px].red-imgI[jy][jx].red)*DIST[jy-py][jx-px]; den+=DIST[jy-py][jx-px]; } } } } else if (color==B){ for (jy=0; jy<sizeY; jy++) { for (jx=0; jx=jy && px>=jx){ num+=rfun(imgI[py][px].blue-imgI[jy][jx].blue)*DIST[py-jy][px-jx]; den+=DIST[py-jy][px-jx]; } else if (jy>py && px>=jx){ num+=rfun(imgI[py][px].blue-imgI[jy][jx].blue)*DIST[jy-py][px-jx]; den+=DIST[jy-py][px-jx]; } else if (jy>py && jx>px){ num+=rfun(imgI[py][px].blue-imgI[jy][jx].blue)*DIST[jy-py][jx-px]; den+=DIST[jy-py][jx-px]; } } } } else if (color==G){ for (jy=0; jy<sizeY; jy++) { for (jx=0; jx=jy && px>=jx){ num+=rfun(imgI[py][px].green-imgI[jy][jx].green)*DIST[py-jy][px-jx]; den+=DIST[py-jy][px-jx]; } else if (jy>py && px>=jx){ num+=rfun(imgI[py][px].green-imgI[jy][jx].green)*DIST[jy-py][px-jx]; den+=DIST[jy-py][px-jx]; } else if (jy>py && jx>px){ num+=rfun(imgI[py][px].green-imgI[jy][jx].green)*DIST[jy-py][jx-px]; den+=DIST[jy-py][jx-px]; } } } } return num/(den*RMAX); }
平方根是一种昂贵的function。 如果你只关心距离的比较(这个距离是否小于这个距离,相等等),你就不需要平方根。
这与许多数字信号处理框架(X-Midas,Midas 2k,PicklingTools)具有幅度(对于复数而言基本上是相同的距离计算)和幅度2(其消除平方根)的原因基本相同。 有时你关心的是事物的比较,不一定是实际价值。
尝试类似的东西
double dx = (jx-qx); double dy = (jy-qy); double dist = sqrt(dx*dx + dy*dy);
如果你只能住在广场上(当你想做一些像距离排序的东西时很有用,你可以使用更高效的
double dx = (jx-qx); double dy = (jy-qy); double dist = dx*dx + dy*dy;
这取决于你的意思。
首先,正如其他人建议您可以调整您的要求并在某些情况下使用距离平方。
接下来,你不会在操作的压缩和额外周期上受益太多,而是考虑在需要处理大量此类操作时加速 – 矢量化,优化内存布局,缓存利用率,并行计算等。但是那些主题需要更多地了解具体的计划需求。
首先,在循环中,您使用相同的参数连续两次调用“距离”:
{ num += foo / distance(...); den += bar / distance(...); }
声明一个新变量“tmp”并将上面的代码更改为
{ tmp = 1 / distance(...); num += foo * distance; den += bar * distance; }
应该大大加快您的代码速度,尽管这可能是您通过更多编译器优化看到的因素2速度增加。 (请找出来!)另请注意,我将这些分区换成了乘法 – 这在计算时间方面也更便宜。
其次,一些编译器具有1 / x ^ 2甚至1 / sqrt(x)的优化。 为了让您的编译器能够使用这种优化,不要计算距离,而是计算反距离。 所以你应该返回声明
return 1 / sqrtf(...);
对于许多优化,您运行代码的CPU也很重要。 特别是如果你想要矢量化(你有SSE吗?AVX?)或并行化(你有多少核心可用?在此期间他们还忙着做其他的东西还是你可以使用它们?)。
如果不尝试它就很难确定,但看起来你最好不要进行矢量化而不是并行化(在不同的内核上运行 – 开销可能会占用你可能获得的任何性能提升)。
如果您正在处理图像,则可以使用快速精确欧氏距离(FEED)变换。 https://ieeexplore.ieee.org/abstract/document/6717981/