/*
 * Realiza a filtragem da media a partir da convolucao entre uma imagem ruidosa
 * com uma mascara de tamanho nxn (n, impar, eh um parametro de entrada)
 *
 * n = 3 (default)
 */
 
#define MAX
#include "lib.h"
#include <math.h>
#include <time.h>

IMAGE gauss(IMAGE x, int n);

IMAGE gauss(IMAGE x, int n)
{
	int i,j,mi,mj;
    double *umask, sigma, g;
    double **mask;
    IMAGE y;
    
    //mask creation
    umask = (double *)malloc(sizeof(double)*n);
    mask = (double **)malloc(sizeof(double *)*n);
    sigma = ((double)n-1)/4;
    for (i = 0; i < n; i++) {
        umask[i] = (1 / (sqrt(2*PI) * sigma) * exp(-((i - ((double)n-1)/2)*(i - ((double)n-1)/2))/(2*sigma*sigma)));
    }

    //multiplication
    for (i = 0; i < n; i++) {
        mask[i] = (double *)malloc(sizeof(double)*n);
        for (j = 0; j < n; j++) {
            mask[i][j] = umask[i] * umask[j];
        }
    }
            
    y = newimage(x->info->nr, x->info->nc);
	for (i = n/2; i < x->info->nr-n/2; i++) {
        for (j = n/2; j < x->info->nc-n/2; j++) {
            g = 0;
            //calculation
            for (mi = -n/2; mi <= n/2; mi++) {
                for (mj = -n/2; mj <= n/2; mj++) {
                    g += (double)(mask[mi + n/2][mj + n/2] * x->data[i + mi][j + mj]);
                }
            }
            y->data[i][j] = (unsigned char)g;
        }
    }
    
    return(y);
}   
 
main(int argc, char *argv[])
{
    int n = 3;
	IMAGE x=0;
    time_t t0, t1;

	if (argc < 2)
	{
	  printf("Usage: avg <image> [<mask_size>]\n");
	  exit(1);
	}

	x = Input_PBM(argv[1]);
	if (x == 0)
	{
	  printf("No input image ('%s')\n", argv[1]);
	  exit(2);
	}

    if (argc > 2) {
        sscanf(argv[2], "%d", &n);
        if (n % 2 <= 0) {
            printf("Mask size must be an odd natural number\n");
            exit(2);
        }
    }

    printf("Average filter with mask of size %d\n", n);

    t0 = time(NULL);
	x = gauss(x, n);
    t1 = time(NULL);

	Output_PBM(x, "gauss.pgm");
	printf("Output is in file 'gauss.pgm'\n");
	printf("Processing time: %lf seconds\n", difftime(t1,t0) / 1000000);
}