Problems with Gaussian filter application

Question

The function should only apply the Gaussian filter to a PPM image, but the result is several overlaid filtered images. Any suggestions ?

[Edit] image before filter | image after filter | expected result

pixel** gauss_filter(pixel** picture, int* L, int* C)
{    
    float novopixel = 0.;
    float peso = 0.;        
    int i, j = 0;        
    int a, b = 0;


    float filter[5][5] =    {{(2.),(4.),(5.),(4.),(2.)},
                            {(4.),(9.),(12.),(9.),(4.)},
                            {(5.),(12.),(15.),(12.),(5.)},
                            {(4.),(9.),(12.),(9.),(4.)},
                            {(2.),(4.),(5.),(4.),(2.)}};

    pixel** new_picture = (pixel**)malloc(*L * sizeof(pixel*)); 

    for(i = 0; i < *L; i++)
    {       
      new_picture[i] = (pixel*) malloc(*C * sizeof(pixel));

      for(j = 0; j < *C; j++)
      {   
        novopixel = 0.;
        peso = 0.;

        for(a = -2; a <= 2; a++)
        {               
            for(b = -2; b <= 2; b++)
            {
                if((i+a >= 0 && i+a < *L) && (j+b >= 0 && j+b < *C))
                { 
                    //Estou usando uma imagem em tons de cinza
                    //por isso calculo um canal e repito o valor
                    //nos outros canais

                    novopixel += (picture[i+a][j+b].r * filter[a+2][b+2]);
                    peso += filter[a+2][b+2];
                }                       
            }
        }

        novopixel = novopixel/peso;

        new_picture[i][j].r = (int)novopixel;
        new_picture[i][j].g = (int)novopixel;
        new_picture[i][j].b = (int)novopixel;
      }
    }

    free(picture);

    picture = new_picture;

    return picture;
}

//======================================================================

pixel** create_picture(FILE *arq, int* L, int* C)
{
    int i, j = 0;
    //armazena tipo do arquivo
    char tipo[3];
    //buffer para leitura do arquivo
    char line[64];
    //controla a leitura do cabecalho
    int flag = 0;
    int componente = 0;

    while ( fgets(line, sizeof line, arq) )
    {   //primeira linha: tipo
        if(flag == 0)
        {
            sscanf(line, "%s", tipo);
            printf("\ntipo: %s",tipo);
            flag += 1;
        }
        //terceira linha: dimensoes 
        if(flag == 2)
        {
            sscanf(line, "%d %d", &*L, &*C);
            printf("\nlinhas: %d colunas: %d", *L, *C);
            flag +=1;
        }
        //valor de componente
        if(flag == 3)
        {
            fgets(line, sizeof line, arq);
            sscanf(line, "%d", &componente);
            printf("\ncomponente: %d", componente);
            if(componente != 255)
            {
                fprintf(stderr, "\nComponente de imagem RGB invalido (%d)\n", componente);
                exit(1);
            }
            break;
        }
        //segunda linha: comentarios
        if(flag == 1)
        {
            fgets(line, sizeof line, arq);
            flag += 1;
        }
    }

    if (tipo[0] != 'P' || tipo[1] != '3')
    {
        fprintf(stderr, "Formato de arquivo invalido ('%c%c')\n", tipo[0], tipo[1]);
        exit(1);
    } 

    pixel** picture = (pixel**)malloc(*L * sizeof(pixel*));

    for (i = 0; i < *L; i++)
    { 
         picture[i] = (pixel*) malloc(*C * sizeof(pixel));

         for (j = 0; j < *C; j++)
         {  
            fgets(line, sizeof line, arq);
            sscanf(line, "%d", &picture[i][j].r);
            fgets(line, sizeof line, arq);          
            sscanf(line, "%d", &picture[i][j].g);
            fgets(line, sizeof line, arq);          
            sscanf(line, "%d", &picture[i][j].b);
         }
    }

    return picture;
}
//=============================================

void output(pixel** picture, int* L, int* C, char* arquivo, char* etapa)
{
    int i = 0;
    int j = 0;

    char name[30] = "";
    strcat(name, etapa);
    strcat(name, arquivo);

    FILE *file = NULL;
    file = fopen(name, "w");

    //formato da imagem
    fprintf(file, "P3\n");

    //comentario
    fprintf(file, "# Arquivo PPM de saida\n");

    //tamanho da imagem
    fprintf(file, "%d %d\n", *L, *C);

    //componente RGB
    fprintf(file, "255\n");

    for(i = 0; i < *L; i++)
    {
        for(j = 0; j < *C; j++) 
        {
            fprintf(file, "%d\n", picture[i][j].r);
            fprintf(file, "%d\n", picture[i][j].g);
            fprintf(file, "%d\n", picture[i][j].b);         
        }

    }

    fclose(file);
}

//============================================

int main()
{
    FILE *arq = NULL;
    pixel** imagem = NULL;
    int linhas = 0;
    int colunas = 0;

    printf("\nEspecifique o arquivo de imagem a ser analisado: ");
    //armazena o nome do arquivo
    char *arquivo = (char*) malloc(25 * sizeof(char));
    scanf("%s", arquivo);

    //Abre o arquivo, e, se nao for possivel, envia mensagem de erro
    if ((arq = fopen(arquivo, "r")) == NULL)
    {
        fprintf(stderr, "\nErro.\nImpossivel abrir arquivo especificado.\n");
        exit(1);
    }   

    //imagem original
    imagem = create_picture(arq, &linhas, &colunas);
    output(imagem, &linhas, &colunas, arquivo, "original_");
    //imagem filtrada
    imagem = gauss_filter(imagem, &linhas, &colunas);
    output(imagem, &linhas, &colunas, arquivo, "gauss_");

    free(arquivo);  
    free(arq);  
    free(imagem);

    return 0;
}

Answer 1

The code of your question is an image processing algorithm called Filtro Espacial.

Filtragem Espacial refers to the image plane, involves direct manipulation of the image pixels using a Máscara Espacial, also known as kernel.

To Filtragem Espacial consists in the displacement of the mask H on the original image I where the sum of products in each region is calculated (v,u), composing a new filtered image I':

The values of the masks are coefficients and the filtering process is similar to a mathematical operation called Convolução.

When it is said that a Filtro Espacial is Gaussian means the mask used in the Convolução is a discrete representation of a Função de Gauss:

In essence, the implementation of its Filtro Espacial is correct! Your code error is at the moment when you read and write the file header .ppm, Notice that you are reading the image dimensions in reverse. You are swapping the number of rows with the number of columns in the image.

Follows a code tested and commented able to apply a Filtro Espacial Gaussiano de Smooth in an image in format .ppm:

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

#define LINE_MAX_LEN    (64)

/* Representa um pixel da imagem */
typedef struct pixel_s
{
    unsigned char red;
    unsigned char green;
    unsigned char blue;
} pixel_t;


/* Representa uma imagem */
typedef struct image_s
{
    pixel_t ** buf;
    int ncols;
    int nrows;
} image_t;


image_t * image_create( int ncols, int nrows )
{
    int i = 0;

    image_t * img = (image_t*) calloc( 1, sizeof(image_t) );

    img->buf = (pixel_t**) calloc( nrows, sizeof(pixel_t*) );

    for( i = 0; i < nrows; i++ )
        img->buf[i] = (pixel_t*) calloc( ncols, sizeof(pixel_t) );

    img->ncols = ncols;
    img->nrows = nrows;

    return img;
}

void image_destroy( image_t * img )
{
    int i = 0;

    for( i = 0; i < img->nrows; i++ )
        free(img->buf[i]);

    free(img->buf);
    free(img);
}

pixel_t * image_read_pixel( image_t * img, int col, int row )
{
    /* Nao permite leitura fora das fronteiras da imagem original */
    if( col >= img->ncols ) col = img->ncols - 1;
    if( row >= img->nrows ) row = img->nrows - 1;
    if( col < 0 ) col = 0;
    if( row < 0 ) row = 0;

    return &img->buf[row][col];
}

image_t * image_gaussian_smooth_filter( image_t * img )
{
    int x, y, col, row, newpx;
    pixel_t * px;
    int sum, div;

    int kernel[5][5] = {{ 2,  4,  5,  4, 2 },
                        { 4,  9, 12,  9, 4 },
                        { 5, 12, 15, 12, 5 },
                        { 4,  9, 12,  9, 4 },
                        { 2,  4,  5,  4, 2 }};

    /* Debug */
    printf( "Filtrando Imagem: tipo=P3; nrows=%d; ncols=%d\n", img->nrows, img->ncols );

    image_t * newimg = image_create( img->ncols, img->nrows );

    for( row = 0; row < img->nrows; row++ )
    {
        for( col = 0; col < img->ncols; col++ )
        {
            sum = 0;
            div = 0;

            for( y = 0; y < 5; y++ )
            {
                for( x = 0; x < 5; x++ )
                {
                    px = image_read_pixel( img,  col + (x - 2), row + (y - 2) );
                    sum += ( px->red *  kernel[y][x] );
                    div += kernel[y][x];
                }
            }

            newpx = sum / div;

            newimg->buf[row][col].red = newpx;
            newimg->buf[row][col].green = newpx;
            newimg->buf[row][col].blue = newpx;
        }
    }

    return newimg;
}

int image_save( const char * file, image_t * img )
{
    int col = 0;
    int row = 0;

    FILE * pf = fopen( file , "w" );

    if(!pf)
        return -1;

    fprintf( pf, "P3\n");
    fprintf( pf, "# Arquivo PPM de saida\n");
    fprintf( pf, "%d %d\n", img->ncols, img->nrows );
    fprintf( pf, "255\n" );

    /* Debug */
    printf( "Gravando Imagem: arquivo=%s; tipo=P3; nrows=%d; ncols=%d; componente=255\n", file, img->nrows, img->ncols );

    for( row = 0; row < img->nrows; row++)
    {
        for( col = 0; col < img->ncols; col++ )
        {
            fprintf( pf, "%u\n", img->buf[row][col].red );
            fprintf( pf, "%u\n", img->buf[row][col].green );
            fprintf( pf, "%u\n", img->buf[row][col].blue );
        }
    }

    fclose(pf);

    return 0;
}

image_t * image_load( const char * file )
{
    int col = 0;
    int row = 0;
    char tipo[3];
    char line[LINE_MAX_LEN];
    int header_index = 0;
    int componente = 0;
    int nrows = 0;
    int ncols = 0;
    image_t * img = NULL;

    FILE * pf = fopen( file, "r" );

    if(!pf)
        return NULL;

    /* Carega o Header do Arquivo PPM */
    while( fgets( line, LINE_MAX_LEN, pf ) )
    {
        if( header_index == 0 )
        {
            sscanf( line, "%s", tipo );
            header_index++;
        }
        else if(header_index == 1 )
        {
            header_index++;
        }
        else if(header_index == 2 )
        {
            sscanf( line, "%d %d", &ncols, &nrows );
            header_index++;
        }
        else if( header_index == 3 )
        {
            sscanf(line, "%d", &componente);
            break;
        }
    }

    /* Debug */
    printf( "Carregando Imagem: arquivo=%s; tipo=%s; nrows=%d; ncols=%d; componente=%d\n", file, tipo, nrows, ncols, componente );

    /* Verifica se o tipo do arquivo eh P3 */
    if( tipo[0] != 'P' || tipo[1] != '3')
    {
        fclose(pf);
        return NULL;
    }

    /* Cria uma imagem em branco */
    img = image_create( ncols, nrows );

    /* Carrega imagem */
    for( row = 0; row < nrows; row++)
    {
        for( col = 0; col < ncols; col++)
        {
            /* Componente Vermelho */
            fgets(line, LINE_MAX_LEN, pf );
            sscanf( line, "%hhu", &img->buf[row][col].red );

            /* Componente Verde */
            fgets( line, LINE_MAX_LEN, pf );
            sscanf( line, "%hhu", &img->buf[row][col].green );

            /* Componente Azul */
            fgets( line, LINE_MAX_LEN, pf );
            sscanf( line, "%hhu", &img->buf[row][col].blue );
        }
    }

    fclose(pf);

    return img;
}

int main( void )
{
    /* Carrega imagem a partir do arquivo */
    image_t * original = image_load( "imagem_cinza.ppm" );

    if(!original)
    {
        printf("Erro ao carregar imagem PPM.\n");
        return 1;
    }

    /* Aplica filtro */
    image_t * filtered = image_gaussian_smooth_filter( original );

    /* Grava imagem filtrada em arquivo */
    int ret = image_save("imagem_filtrada.ppm", filtered );

    if(ret < 0)
    {
        printf("Erro ao salvar imagem PPM.\n");

        /* Libera memoria ocupada pelas imagens */
        image_destroy(filtered);
        image_destroy(original);

        return 1;
    }

    /* Libera memoria ocupada pelas imagens */
    image_destroy(filtered);
    image_destroy(original);

    return 0;
}

Exit:

Carregando Imagem: arquivo=imagem_cinza.ppm; tipo=P3; nrows=480; ncols=610; componente=255
Filtrando Imagem: tipo=P3; nrows=480; ncols=610
Gravando Imagem: arquivo=imagem_filtrada.ppm; tipo=P3; nrows=480; ncols=610; componente=255

Original Image: (Download . ppm)

Filtered Image: (Download . ppm)

References:

1. http://dcm.ffclrp.usp.br/~murta/PIM/PIM_2.pdf

2. http://www.facom.ufu.br/~backes/gsi058/Aula06-Filtragemespacial.pdf

3. http://www.coe.utah.edu/~cs4640/slides/Lecture5.pdf