Fastest way to calculate Euclidean distance in c

Question

I need to calculate euclidean distance between two points in the fastest way possible. In C.

My code is this and seems a little bit slow:

float distance(int py, int px, int jy, int jx){
     return sqrtf((float)((px)*(px)+(py)*(py)));
}

Thanks in advance.

EDIT:

I'm sorry I hadn't been clear. I'd better specify the context: I'm working with images and I need the euclidean distance from each pixel to all the other pixels. So I have to calculate it a lot of times. I can't use the square of the distance. I'll add a little bit more code to be more clear:

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)));
}

That's what I have t do. And I have to do it with all the pixels (px, py)

EDIT2: I'm sorry i wasn't clear but I tried to keep the question as general as I could. I'm writing a program to process images with an algorithm. The big problem is the time because I have to do it really really fast. Now what I need to optimize is this function: `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;

} `

This function is called in a loop for every pixel(px; py) so it's called a lot of times and it takes al lot of time to compute this. The rfun function is not optimizable, because is already really simple and fast. what I need to do is to make faster the distance function.

1)I tried hypotf but it was slower than the distance function

2)I increased the optimization settings of the compiler and that made the process 2x faster!

3) I tried with a macro #define DIST(x, y) sqrtf((float)((x)*(x)+(y)*(y))) but nothing changed (as I expected)

EDIT3:

In the end I found that the fastest way was to calculate all the possible distances and store them in an array before starting computing in a loop the normC function. To make it faster I calculated the inverses of the distances so that I could use the product and not the quotient:

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<sizeX; jx++) {
            if (jx==px && jy==py) {
                ;
            }
            else if (py>=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<sizeX; jx++) {
            if (jx==px && jy==py) {
                ;
            }
            else if (py>=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<sizeX; jx++) {
            if (jx==px && jy==py) {
                ;
            }
            else if (py>=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);
}

Answer 1

The square root is an expensive function. If you just care about how distances compare (is this distance less than this distance, equal, etc.), you don't need the square root.

It's basically the same reason why many Digital Signal Processing frameworks (X-Midas, Midas 2k, PicklingTools) have magnitude (which is basically the same distance computation for complex numbers) and magnitude2 (which elides the square root). Sometimes all you care about is how things compare, not necessarily the actual value.

Answer 2

try something like

double dx = (jx-qx);
double dy = (jy-qy);
double dist = sqrt(dx*dx + dy*dy);

if you can live with just the square (which is useful when you just want to do something like sort by distance, you can use the much more efficient

double dx = (jx-qx);
double dy = (jy-qy);
double dist = dx*dx + dy*dy;

Answer 3

It depends on what you mean fast.

First, as other people suggested that you might adjust your requirement and live with distance squared in some situations.

Next, you would not benefit much on squeezing and extra cycle on the operation, but instead consider to speed up when you need to process huge amount of such operations - vectorization, optimizing memory layout, cache utilization, parallel computing, etc. But those topic require more understanding on specific program needs.

Answer 4

First, in your loop you call "distance" twice in a row, with the same argument:

{
  num += foo / distance(...);
  den += bar / distance(...);
}

Declaring a new variable "tmp" and changing above code to

{
  tmp = 1 / distance(...);
  num += foo * distance;
  den += bar * distance;
}

should speed up your code considerably, although this is probably the factor 2 speed increase you are seeing with more compiler optimization. (Please find out!) Also note that I swapped the divisions for multiplications - that is cheaper in terms of computing time as well.

Second, some compilers have optimizations for 1/x^2 and maybe even 1/sqrt(x). In order to allow your compiler to make use of such optimizations, don't calculate the distance, but the inverse distance. So you return statement should be

return 1 / sqrtf(...);

For many optimizations, the CPU you are running the code on is important, too. Especially so if you wish to vectorize (do you have SSE? AVX?) or parallelize (how many cores do you have available? Are they busy doing other stuff in the meantime or can you use them?).

It's hard to be sure without trying it, but it looks like you'd be better off vectorizing rather than parallelizing (running on separate cores -- the overhead will probably eat up any performance gains you might get).

Answer 5

if you are working on images you may use the Fast exact Euclidean distance (FEED) transforms. https://ieeexplore.ieee.org/abstract/document/6717981/

Answer 6

optimal loop

np=length(pixelsCand);
dc=zeros(np);


for(ir=1:np)
    for(ic=ir+1:np)
        %dc(ir,ir)=0
        dc(ir,ic)=sqrt(sum((pixelsCand(ir,:)-pixelsCand(ic,:)).^2));
        dc(ic,ir)= dc(ir,ic);
    end
end