#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <algorithm>
#define size 256
#define max(a,b) (((a) > (b)) ? (a) : (b))
#define min(a,b) (((a) < (b)) ? (a) : (b))
using namespace std;
void nlm(unsigned char (*I)[size][3],unsigned char (*J)[size][3],int sear,int sim,int h);
void gau_kernel(float *,int simw);
int main(int argc, char *argv[])
{
unsigned char image[size][size][3];
unsigned char result[size][size][3];
int sw=21,simw=7,h=2;
int pp;
// Write you code HERE
char path[1025];
cout<<"path of lena_noisy:";
cin>>path;
cout<<"if it is color please input 3, else input 1:";
cin>>pp;
int sz = size*size*pp;
FILE *fp;
fp = fopen(path,"rb");
unsigned char *pData = new unsigned char[sz];
int num = fread(pData,1,sz,fp);;
fclose(fp);
int i,j;
//unsigned char * pTemp;
for (i=0; i<size; ++i)
{
for (j=0; j<size; ++j)
{
image[i][j][0] = pData[i*size*pp+j*pp];
image[i][j][1] = pData[i*size*pp+j*pp+1];
image[i][j][2] = pData[i*size*pp+j*pp+2];
}
}
nlm(image,result,21,7,2);
for (i=0; i<size; ++i)
{
for (j=0; j<size; ++j)
{
pData[i*size*pp+pp*3] = result[i][j][0];
pData[i*size*pp+pp*3+1] = result[i][j][1];
pData[i*size*pp+pp*3+2] = result[i][j][2];
}
}
cout<<"input output path of lena_noisy:";
cin>>path;
fp = fopen(path,"wb");
num = fwrite(pData,1,sz,fp);
fclose(fp);
delete []pData;
//nlm(image,result,sw,simw,h);
}
void nlm(unsigned char (*I)[size][3],unsigned char (*J)[size][3],int sw,int simw,int h)
{
////////////////////////////////////////////////////////////////////////////////////////
// This code is based on follwoing paper:
// I. "A non-local algorithm for image denoising", Buades, A. and Coll, B. and Morel,
// J.M.,Proceedings of the 2005 IEEE Computer Society Conference on Computer Vision
// and Pattern Recognition, vol. 2, pp. 60-65, 2005
////////////////////////////////////////////////////////////////////////////////////////
//
// "sw", "simw" can be odd numbers ONLY
// "sw" is the search window
// "simw" is the similarity window
// "h" is the filter paprameter as defined in [I]
//
////////////////////////////////////////////////////////////////////////////////////////
float *kernel;
int i,j,k,x,y,a,b,c;
int temp_i,temp_j,temp_x,temp_y;
int rmin,rmax,cmin,cmax;
int win=(sw-1)/2;
int s=(simw-1)/2;
float *weight[3];
float sum_weight[3]={0,0,0};
float temp_pixel[3];
float temp_diff;
int index_sw,index_simw;
FILE *file;
weight[0]=(float *)calloc(sw*sw,sizeof(float));
weight[1]=(float *)calloc(sw*sw,sizeof(float));
weight[2]=(float *)calloc(sw*sw,sizeof(float));
kernel=(float *)calloc(simw*simw,sizeof(float));
gau_kernel(kernel,simw);
h=h*h; // This is done because in [I], h^2 is used for calculation
for(i=0;i<size;i++) {
for(j=0;j<size;j++) {
rmin=max(i-win,0);
rmax=min(i+win,size-1);
cmin=max(j-win,0);
cmax=min(j+win,size-1);
index_sw=0;
temp_pixel[0]=0;
temp_pixel[1]=0;
temp_pixel[2]=0;
sum_weight[0]=0;
sum_weight[1]=0;
sum_weight[2]=0;
//Start searching within the search window
for(x=rmin;x<rmax;x++) {
for(y=cmin;y<cmax;y++) {
index_simw=0;
weight[0][index_sw]=0;
weight[1][index_sw]=0;
weight[2][index_sw]=0;
//In the similiarity window:
for(a=-s;a<s+1;a++) {
temp_i=max(0,i+a);
temp_i=min(temp_i,size-1);
temp_x=max(0,x+a);
temp_x=min(temp_x,size-1);
for(b=-s;b<s+1;b++) {
temp_j=max(0,j+b);
temp_j=min(temp_j,size-1);
temp_y=max(0,y+b);
temp_y=min(temp_y,size-1);
for(k=0;k<3;k++) {
temp_diff=(float)I[temp_x][temp_y][k]-(float)I[temp_i][temp_j][k];
temp_diff=temp_diff*temp_diff;
temp_diff*=kernel[index_simw];
weight[k][index_sw]+=temp_diff;
}
index_simw++;
}
}
for (k=0;k<3;k++) {
weight[k][index_sw]=exp(-weight[k][index_sw]/h);
sum_weight[k]+=weight[k][index_sw];
temp_pixel[k]+=weight[k][index_sw]*(float)I[x][y][k];//temporary pixel value before normalizing
}
index_sw++;
}
}
//Finishing calcuating all weights within a window
for(k=0;k<3;k++) {
for(c=0;c<index_sw;c++)
J[i][j][k]=(unsigned char)(temp_pixel[k]/sum_weight[k]);
}
}
}
}
void gau_kernel(float *k,int w)
{
int i,j,d;
float v;
int f;
f=(w-1)/2;
for (d=1;d<f+1;d++) {
v=(2*d+1)*(2*d+1);
v=1/v;
for (i=-d;i<d+1;i++) {
for(j=-d;j<d+1;j++)
k[(f-i)*w+f-j]=k[(f-i)*w+f-j]+v;
}
}
}