void prefilt()

in lear_gist-1.2/gist.c [286:415]
94 lines of code
13 McCabe index (conditional complexity)

/*static*/ void prefilt(image_t *src, int fc)
{
    fftw_lock();

    int i, j;

    /* Log */
    for(j = 0; j < src->height; j++)
    {
        for(i = 0; i < src->width; i++) {
            src->data[j*src->stride+i] = log(src->data[j*src->stride+i]+1.0f);
        }
    }

    image_t *img_pad = image_add_padding(src, 5);

    /* Get sizes */
    int width = img_pad->width;
    int height = img_pad->height;
    int stride = img_pad->stride;

    /* Alloc memory */
    float *fx  = (float *) fftwf_malloc(width*height*sizeof(float));
    float *fy  = (float *) fftwf_malloc(width*height*sizeof(float));
    float *gfc = (float *) fftwf_malloc(width*height*sizeof(float));
    fftwf_complex *in1 = (fftwf_complex *) fftwf_malloc(width*height*sizeof(fftwf_complex));
    fftwf_complex *in2 = (fftwf_complex *) fftwf_malloc(width*height*sizeof(fftwf_complex));
    fftwf_complex *out = (fftwf_complex *) fftwf_malloc(width*height*sizeof(fftwf_complex));

    /* Build whitening filter */
    float s1 = fc/sqrt(log(2));
    for(j = 0; j < height; j++)
    {
        for(i = 0; i < width; i++)
        {
            in1[j*width + i][0] = img_pad->data[j*stride+i];
            in1[j*width + i][1] = 0.0f;

            fx[j*width + i] = (float) i - width/2.0f;
            fy[j*width + i] = (float) j - height/2.0f;

            gfc[j*width + i] = exp(-(fx[j*width + i]*fx[j*width + i] + fy[j*width + i]*fy[j*width + i]) / (s1*s1));
        }
    }

    fftshift(gfc, width, height);

    /* FFT */
    fftwf_plan fft1 = fftwf_plan_dft_2d(width, height, in1, out, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_unlock();
    fftwf_execute(fft1);
    fftw_lock();

    /* Apply whitening filter */
    for(j = 0; j < height; j++)
    {
        for(i = 0; i < width; i++)
        {
            out[j*width+i][0] *= gfc[j*width + i];
            out[j*width+i][1] *= gfc[j*width + i];
        }
    }

    /* IFFT */
    fftwf_plan ifft1 = fftwf_plan_dft_2d(width, height, out, in2, FFTW_BACKWARD, FFTW_ESTIMATE);
    fftw_unlock();
    fftwf_execute(ifft1);
    fftw_lock();

    /* Local contrast normalisation */
    for(j = 0; j < height; j++)
    {
        for(i = 0; i < width; i++)
        {
            img_pad->data[j*stride + i] -= in2[j*width+i][0] / (width*height);

            in1[j*width + i][0] = img_pad->data[j*stride + i] * img_pad->data[j*stride + i];
            in1[j*width + i][1] = 0.0f;
        }
    }

    /* FFT */
    fftwf_plan fft2 = fftwf_plan_dft_2d(width, height, in1, out, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_unlock();
    fftwf_execute(fft2);
    fftw_lock();

    /* Apply contrast normalisation filter */
    for(j = 0; j < height; j++)
    {
        for(i = 0; i < width; i++)
        {
            out[j*width+i][0] *= gfc[j*width + i];
            out[j*width+i][1] *= gfc[j*width + i];
        }
    }

    /* IFFT */
    fftwf_plan ifft2 = fftwf_plan_dft_2d(width, height, out, in2, FFTW_BACKWARD, FFTW_ESTIMATE);
    fftw_unlock();
    fftwf_execute(ifft2);
    fftw_lock();

    /* Get result from contrast normalisation filter */
    for(j = 0; j < height; j++)
    {
        for(i = 0; i < width; i++) {
            img_pad->data[j*stride+i] = img_pad->data[j*stride + i] / (0.2f+sqrt(sqrt(in2[j*width+i][0]*in2[j*width+i][0]+in2[j*width+i][1]*in2[j*width+i][1]) / (width*height)));
        }
    }

    image_rem_padding(src, img_pad, 5);

    /* Free */
    fftwf_destroy_plan(fft1);
    fftwf_destroy_plan(fft2);
    fftwf_destroy_plan(ifft1);
    fftwf_destroy_plan(ifft2);

    image_delete(img_pad);

    fftwf_free(in1);
    fftwf_free(in2);
    fftwf_free(out);
    fftwf_free(fx);
    fftwf_free(fy);
    fftwf_free(gfc);

    fftw_unlock();
}