void color_prefilt()

in lear_gist-1.2/gist.c [419:592]
135 lines of code
13 McCabe index (conditional complexity)

void color_prefilt(color_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->c1[j*src->width+i] = log(src->c1[j*src->width+i]+1.0f);
            src->c2[j*src->width+i] = log(src->c2[j*src->width+i]+1.0f);
            src->c3[j*src->width+i] = log(src->c3[j*src->width+i]+1.0f);
        }
    }

    color_image_t *img_pad = color_image_add_padding(src, 5);

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

    /* 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 *ina1 = (fftwf_complex *) fftwf_malloc(width*height*sizeof(fftwf_complex));
    fftwf_complex *ina2 = (fftwf_complex *) fftwf_malloc(width*height*sizeof(fftwf_complex));
    fftwf_complex *ina3 = (fftwf_complex *) fftwf_malloc(width*height*sizeof(fftwf_complex));
    fftwf_complex *inb1 = (fftwf_complex *) fftwf_malloc(width*height*sizeof(fftwf_complex));
    fftwf_complex *inb2 = (fftwf_complex *) fftwf_malloc(width*height*sizeof(fftwf_complex));
    fftwf_complex *inb3 = (fftwf_complex *) fftwf_malloc(width*height*sizeof(fftwf_complex));
    fftwf_complex *out1 = (fftwf_complex *) fftwf_malloc(width*height*sizeof(fftwf_complex));
    fftwf_complex *out2 = (fftwf_complex *) fftwf_malloc(width*height*sizeof(fftwf_complex));
    fftwf_complex *out3 = (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++)
        {
            ina1[j*width + i][0] = img_pad->c1[j*width+i];
            ina2[j*width + i][0] = img_pad->c2[j*width+i];
            ina3[j*width + i][0] = img_pad->c3[j*width+i];
            ina1[j*width + i][1] = 0.0f;
            ina2[j*width + i][1] = 0.0f;
            ina3[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 fft11 = fftwf_plan_dft_2d(width, height, ina1, out1, FFTW_FORWARD, FFTW_ESTIMATE);
    fftwf_plan fft12 = fftwf_plan_dft_2d(width, height, ina2, out2, FFTW_FORWARD, FFTW_ESTIMATE);
    fftwf_plan fft13 = fftwf_plan_dft_2d(width, height, ina3, out3, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_unlock();
    fftwf_execute(fft11);
    fftwf_execute(fft12);
    fftwf_execute(fft13);
    fftw_lock();

    /* Apply whitening filter */
    for(j = 0; j < height; j++)
    {
        for(i = 0; i < width; i++)
        {
            out1[j*width+i][0] *= gfc[j*width + i];
            out2[j*width+i][0] *= gfc[j*width + i];
            out3[j*width+i][0] *= gfc[j*width + i];

            out1[j*width+i][1] *= gfc[j*width + i];
            out2[j*width+i][1] *= gfc[j*width + i];
            out3[j*width+i][1] *= gfc[j*width + i];
        }
    }

    /* IFFT */
    fftwf_plan ifft11 = fftwf_plan_dft_2d(width, height, out1, inb1, FFTW_BACKWARD, FFTW_ESTIMATE);
    fftwf_plan ifft12 = fftwf_plan_dft_2d(width, height, out2, inb2, FFTW_BACKWARD, FFTW_ESTIMATE);
    fftwf_plan ifft13 = fftwf_plan_dft_2d(width, height, out3, inb3, FFTW_BACKWARD, FFTW_ESTIMATE);
    fftw_unlock();
    fftwf_execute(ifft11);
    fftwf_execute(ifft12);
    fftwf_execute(ifft13);
    fftw_lock();

    /* Local contrast normalisation */
    for(j = 0; j < height; j++)
    {
        for(i = 0; i < width; i++)
        {
            img_pad->c1[j*width+i] -= inb1[j*width+i][0] / (width*height);
            img_pad->c2[j*width+i] -= inb2[j*width+i][0] / (width*height);
            img_pad->c3[j*width+i] -= inb3[j*width+i][0] / (width*height);

            float mean = (img_pad->c1[j*width+i] + img_pad->c2[j*width+i] + img_pad->c3[j*width+i])/3.0f;

            ina1[j*width+i][0] = mean*mean;
            ina1[j*width+i][1] = 0.0f;
        }
    }

    /* FFT */
    fftwf_plan fft21 = fftwf_plan_dft_2d(width, height, ina1, out1, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_unlock();
    fftwf_execute(fft21);
    fftw_lock();

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

    /* IFFT */
    fftwf_plan ifft2 = fftwf_plan_dft_2d(width, height, out1, inb1, 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++)
        {
            float val = sqrt(sqrt(inb1[j*width+i][0]*inb1[j*width+i][0]+inb1[j*width+i][1]*inb1[j*width+i][1]) / (width*height));

            img_pad->c1[j*width+i] /= (0.2f+val);
            img_pad->c2[j*width+i] /= (0.2f+val);
            img_pad->c3[j*width+i] /= (0.2f+val);
        }
    }

    color_image_rem_padding(src, img_pad, 5);

    /* Free */
    fftwf_destroy_plan(fft11);
    fftwf_destroy_plan(fft12);
    fftwf_destroy_plan(fft13);
    fftwf_destroy_plan(ifft11);
    fftwf_destroy_plan(ifft12);
    fftwf_destroy_plan(ifft13);
    fftwf_destroy_plan(fft21);
    fftwf_destroy_plan(ifft2);

    color_image_delete(img_pad);

    fftwf_free(ina1);
    fftwf_free(ina2);
    fftwf_free(ina3);
    fftwf_free(inb1);
    fftwf_free(inb2);
    fftwf_free(inb3);
    fftwf_free(out1);
    fftwf_free(out2);
    fftwf_free(out3);
    fftwf_free(fx);
    fftwf_free(fy);
    fftwf_free(gfc);

    fftw_unlock();
}