qsort in C is an unexpected bottle neck in my program

Question

I am working on some code to align astronomical images. I am writing it in C99. For some reason, my algorithm for detecting the stars in an image was running much slower than expected. I basically need to ignore all pixels which are below the 99th percentile (because stars are just small bright points). To calculate the 99th percentile I applied a qsort to a copy of the pixels in the image. When I profiled the code, it said that it was spending 75% of the time executing the compare_quantum function used by qsort. The entire detection of stars takes about 3 seconds.

I had originally written the code in C++ and the same algorithm took about 0.2 seconds. I am guessing that the reason this is happening is that unlike C++, C can't just inline the call to the compare function like C++ can with std::sort.

I could write my own sort function, but I was just wondering if anyone had any other ideas to make this go faster. I have calls to qsort elsewhere in the code and I am thinking maybe I need to get rid of all of them.

I am using gcc 5.2. The first qsort in stars_map is the bottle neck. Also quantum_t is just a typedef for uint16_t.

int compare_quantum(const void* a, const void* b)
{
    return (*(quantum_t*)a > *(quantum_t*)b) - (*(quantum_t*)a < *(quantum_t*)b);
}

void star_register(quantum_t* grey, double* rowint, double* colint, double* lum, size_t row, size_t col, size_t w, size_t h)
{
    size_t gi = row * w + col;
    if(!(row >= 0 && col >= 0 && row < h && col < w && grey[gi]))
        return;

    *rowint += grey[gi] * row;
    *colint += grey[gi] * col;
    *lum += grey[gi];
    grey[gi] = 0;

    for(int dr = -1; dr <= 1; dr++)
    {
        for(int dc = -1; dc <= 1; dc++)
        {
            if(dc == 0 && dr == 0)
                continue;
            star_register(grey, rowint, colint, lum, row + dr, col + dc, w, h);
        }
    }
}

stars_t* stars_map(image_t* img, float detection_percentile)
{
    assert(img);

    quantum_t* grey = NULL;
    quantum_t* sorted = NULL;
    star_t* stars = NULL;
    size_t nstars = 0;
    size_t stars_alloc = 0;

    grey = malloc(sizeof(quantum_t) * img->w * img->h);
    if(grey == NULL) goto fail;

    sorted = malloc(sizeof(quantum_t) * img->w * img->h);
    if(sorted == NULL) goto fail;

    for(size_t i = 0; i < img->w * img->h; i++)
        sorted[i] = grey[i] = ((uint32_t)img->px[i].red + (uint32_t)img->px[i].green + (uint32_t)img->px[i].blue) / 3;
    //this qsort is the issue
    qsort(sorted, img->w * img->h, sizeof(quantum_t), compare_quantum);
    quantum_t cut = sorted[(size_t)(img->w * img->h * detection_percentile)];
    free(sorted);
    sorted = NULL;

    for(size_t i = 0; i < img->w * img->h; i++)
        grey[i] = clampq((int32_t)grey[i] - cut);

    for(size_t i = 0; i < img->h; i++)
    {
        for(size_t j = 0; j < img->w; j++)
        {
            if(grey[i * img->w + j])
            {
                if(nstars == stars_alloc)
                {
                    stars = realloc(stars, (stars_alloc += 500) * sizeof(star_t));
                    if(!stars) goto fail;
                }
                double rowint = 0.0;
                double colint = 0.0;
                double lum = 0.0;
                star_register(grey, &rowint, &colint, &lum, i, j, img->w, img->h);
                stars[nstars++] = (star_t){.x = colint / lum, .y = rowint / lum, .lum = lum};
            }
        }
    }

    free(grey);

    qsort(stars, nstars, sizeof(star_t), star_compare);

    stars_t* result = malloc(sizeof(stars_t) + nstars * sizeof(star_t));
    if(result == NULL) goto fail;
    result->npairs = nstars;
    memcpy(result->stars, stars, sizeof(star_t) * nstars);
    free(stars);


    return result;

    fail:
    if(grey) free(grey);
    if(sorted) free(sorted);
    if(stars) free(stars);

    return NULL;
}

I originally thought that the recursive call to star_register would be the performance hit, but it barely matters in the profile.

Answer 1

The issue was that I had forgotten that I was using std::nth_element not std::sort in the c++ version. That is why the code was slow. I wrote a qselect and now the entire program is about the same speed.

quantum_t quantum_qselect(quantum_t *v, size_t len, size_t k)
{
    size_t i, st;

    for(st = i = 0; i < len - 1; i++)
    {
        if(v[i] > v[len - 1])
            continue;
        swap(quantum_t, v[i], v[st]);
        st++;
    }

    swap(quantum_t, v[len - 1], v[st]);

    return k == st ? v[st] : st > k ? quantum_qselect(v, st, k) : quantum_qselect(v + st, len - st, k - st);
}