median.cpp

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>

double
sec(void)
{
    struct timespec ts;

    clock_gettime(CLOCK_MONOTONIC, &ts);

    return ts.tv_sec + ts.tv_nsec/1000000000.0;
}

void 
test(unsigned char *out,
     unsigned char *in,
     int pitch,
     int w, int h)
{

#pragma omp parallel for
    for (int yi=3; yi<h-3; yi++) {
        int hist[256];
        int hist16[16];

        memset(hist, 0, sizeof(hist));
        memset(hist16, 0, sizeof(hist16));

        unsigned char *line0 = in + (yi-3) * pitch;
        unsigned char *line1 = in + (yi-2) * pitch;
        unsigned char *line2 = in + (yi-1) * pitch;
        unsigned char *line3 = in + (yi+0) * pitch;
        unsigned char *line4 = in + (yi+1) * pitch;
        unsigned char *line5 = in + (yi+2) * pitch;
        unsigned char *line6 = in + (yi+3) * pitch;

        unsigned char *out_line = out + yi*pitch;

        for (int xi=0; xi<6; xi++) {
            hist[line0[xi]]++;
            hist[line1[xi]]++;
            hist[line2[xi]]++;
            hist[line3[xi]]++;
            hist[line4[xi]]++;
            hist[line5[xi]]++;
            hist[line6[xi]]++;

            hist16[line0[xi]>>4]++;
            hist16[line1[xi]>>4]++;
            hist16[line2[xi]>>4]++;
            hist16[line3[xi]>>4]++;
            hist16[line4[xi]>>4]++;
            hist16[line5[xi]>>4]++;
            hist16[line6[xi]>>4]++;
        }

        for (int xi=3; xi<w-3; xi++) {
            int med = 25;       // ceil(49/2)
            int sum16 = 0;
            int i= 0;

            unsigned int v0 = line0[xi+3];
            unsigned int v1 = line1[xi+3];
            unsigned int v2 = line2[xi+3];
            unsigned int v3 = line3[xi+3];
            unsigned int v4 = line4[xi+3];
            unsigned int v5 = line5[xi+3];
            unsigned int v6 = line6[xi+3];

            hist[v0]++;
            hist[v1]++;
            hist[v2]++;
            hist[v3]++;
            hist[v4]++;
            hist[v5]++;
            hist[v6]++;

            hist16[v0>>4]++;
            hist16[v1>>4]++;
            hist16[v2>>4]++;
            hist16[v3>>4]++;
            hist16[v4>>4]++;
            hist16[v5>>4]++;
            hist16[v6>>4]++;

            for (int i16=0; i16<16; i16++) {
                int cur = hist16[i16];
                if (sum16 + cur >= med) {
                    int sum = sum16;

                    for (i=i16*16; ;i++) {
                        sum += hist[i];
                        if (sum >= med) {
                            goto next;
                        }
                    }

                }

                sum16 += cur;
            }

        next:
            out_line[xi] = i;

            v0 = line0[xi-3];
            v1 = line1[xi-3];
            v2 = line2[xi-3];
            v3 = line3[xi-3];
            v4 = line4[xi-3];
            v5 = line5[xi-3];
            v6 = line6[xi-3];

            hist[v0]--;
            hist[v1]--;
            hist[v2]--;
            hist[v3]--;
            hist[v4]--;
            hist[v5]--;
            hist[v6]--;

            hist16[v0>>4]--;
            hist16[v1>>4]--;
            hist16[v2>>4]--;
            hist16[v3>>4]--;
            hist16[v4>>4]--;
            hist16[v5>>4]--;
            hist16[v6>>4]--;
        }
    }
}
     

int main(int argc, char **argv)
{
    int w=1920, h=1080;
    int pitch=1920;
    unsigned char *in = (unsigned char*)malloc(pitch * h);
    unsigned char *out = (unsigned char*)malloc(pitch * h);

    if (argc >= 3) {
        pitch = w = atoi(argv[1]);
        h = atoi(argv[2]);
    }

    memset(in, 255, pitch*h);
    memset(out, 0, pitch*h);

    test(out, in, pitch, w, h);

    double t0 = sec();
    test(out, in, pitch, w, h);
    double t1 = sec();

    printf("%f\n", t1-t0);
}