/*
 * 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
#define SQ(x) x*x
#include "lib.h"
#include <math.h>

//void errorstats(IMAGE x, IMAGE y, double *erms, double *snr)
//{
//	int i,j;
//    double e=0, s=0;
//    
//    if (!erms) erms = (double *)malloc(sizeof(double));
//    if (!snr) snr = (double *)malloc(sizeof(double));
//    
//    for (i=0; i<x->info->nr; i++) {
//        for (j = 0; j < x->info->nc; j++) {
//            e += (double)SQ((double)x->data[i][j]-(double)y->data[i][j]);
//            s += (double)SQ((double)x->data[i][j]);
//        }
//    }
//    s /= e;
//    e = sqrt(e / (x->info->nr * x->info->nc));
//    
//    erms[0] = e;
//    snr[0] = s;
//    
//    return;
//}   
// 
double erms(IMAGE x, IMAGE y)
{
	int i,j;
    double e=0;
    
    for (i=0; i<x->info->nr; i++) {
        for (j = 0; j < x->info->nc; j++) {
            e += (double)SQ((double)x->data[i][j]-(double)y->data[i][j]);
        }
    }
    e = sqrt(e / (double)(x->info->nr * x->info->nc));
    
    return(e);
}   
 
double snr(IMAGE x, IMAGE y)
{
	int i,j;
    double e=0,s=0;
    
    for (i=0; i<x->info->nr; i++) {
        for (j = 0; j < x->info->nc; j++) {
            e += (double)SQ((double)x->data[i][j]-(double)y->data[i][j]);
            s += (double)SQ((double)x->data[i][j]);
        }
    }
    s /= e;
    
    return(s);
}   
 
main(int argc, char *argv[])
{
    int n = 3;
	IMAGE x=0, y=0;
    double e, s;

	if (argc < 3)
	{
	  printf("Usage: error <imagem_filtrada> <imagem_ruido>\n");
	  exit(1);
	}

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

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

    if (x->info->nr != y->info->nr || x->info->nc != y->info->nc) {
        printf("Images must have the same dimensions\n");
        exit(2);
    }

	printf("Error stats...");
	//errorstats(x, y, erms, snr);
    e = erms(x,y);
    s = snr(x,y);

	printf("ERMS = %lf; SNR = %lf'\n", e,s);
}