nhwfm.c

图像置乱代码

#include <tina/all_tina.h>
#include"tch.h"
#define pi 3.141592654

static void sort(int array[], int n)
{
	int i, j;
	int t;

	for(i=0; i<n-1; i++)
		for(j=i+1; j<n; j++)
			if(array[i]>array[j])
			{
				t=array[i];
				array[i]=array[j];
				array[j]=t;
			}
/*	for(i=0; i<n; i++)		
	printf("%d\t", array[i]);		
	printf("\n");*/
}


void nhwfm_proc(void)
{
	Imrect *im;
	Imregion *roi;
	int lx, ux, ly, uy;
	int height, width;
	int type;
	int i, j, k, t, n;
	int *hRow, *arRow;
	int hist[256]={0};
	int rechist[256];
	double temphist[256];
	int r;
	static int s1, s2, s3;
	double u1, u2, u3, sigma1, sigma2, sigma3;
	int **imageBuf1, **imageBuf2, **imageBuf3;
	double t1, t2, t3;
	double **index;
	double sum, indexsum;
	int cnt;
	int array[9];
	
	if(stack_check_types(IMRECT, NULL)==false)		
	{
		error("imcalc_sum: wrong types on stack", warning);
		return;
	}

	im=(Imrect *)stack_pop(&type);

	roi=im->region;
	if(roi==NULL) return;
	lx=roi->lx;
	ux=roi->ux;
	ly=roi->ly;
	uy=roi->uy;

	width=im->width;
	height=im->height;
	hRow=(int *)malloc(sizeof(int)*(ux-lx));
	imageBuf1=int_matrix_alloc(height, width);
	imageBuf2=int_matrix_alloc(height, width);
	imageBuf3=int_matrix_alloc(height, width);
	index=db_matrix(height, width);

	for(i=ly; i<uy; ++i)
	{
		im_get_row(hRow, im, i, lx, ux);
			for(j=lx; j<ux; ++j)
				imageBuf1[i][j]=hRow[j];
	}

	for(i=0; i<height; i++)
		for(j=0; j<width; j++)
		{
			n=imageBuf1[i][j];
			hist[n]++;
		}

	r=(int)((double)hist[255]/(double)(height*width)*100+0.5);
	printf("noise rate is %d%\n", r);
	printf("%f\n", r/100.0);

	for(i=0; i<255; i++)
		temphist[i]=(double)hist[i]/(1-r/100.0);
	temphist[255]=0.0;	

/*	fixed_spline_smoothing(temphist, 256, 30);*/

/*	printf("\n");
	for(i=0; i<256; i++)
	{
		rechist[i]=(int)(temphist[i]+0.5);
		if(i%8==0) printf("\n");
		printf("%d,%d\t", i,rechist[i]);
	}	
*/
for(t1=49.0; t1<130.0; t1+=0.1)
if(1401.682928*pow(pi*10.3*10.3, -0.25)*exp(-(t1-49.0)*(t1-49.0)/(2*39.9*39.9))<=4774.235939*pow(pi*39.9*39.9, -0.25)*exp(-(t1-130.0)*(t1-130.0)/(2*39.9*39.9)))
{printf("t1=%f\n", t1); break;}

for(t2=130.0; t2<207.0; t2+=0.1)
if(4774.235939*pow(pi*39.9*39.9, -0.25)*exp(-(t2-130.0)*(t2-130.0)/(2*39.9*39.9))<=469.964662*pow(pi*8.7*8.7, -0.25)*exp(-(t2-207.0)*(t2-207.0)/(2*8.7*8.7)))
{printf("t2=%f\n", t2); break;}

for(t3=207.0; t3<255; t3++)
if(469.964662*pow(pi*8.7*8.7, -0.25)*exp(-(t3-207.0)*(t3-207.0)/(2*8.7*8.7))<=4774.235939*pow(pi*39.9*39.9, -0.25)*exp(-(t3-130.0)*(t3-130.0)/(2*39.9*39.9)))
{printf("t3=%f\n", t3); break;}

Send(imageBuf1, imageBuf2, height, width);		
for(i=0; i<height; i++)
	for(j=0; j<width; j++)
		if(imageBuf2[i][j]==255) imageBuf2[i][j]=0;

for(i=0; i<height; i++)
	for(j=0; j<width; j++)
	{
		if(imageBuf2[i][j]<=t1) index[i][j]=1401.682928*pow(pi*10.3*10.3, -0.25)*exp(-(imageBuf2[i][j]-49.0)*(imageBuf2[i][j]-49.0)/(2*10.3*10.3));
		else if(imageBuf2[i][j]<=t2&&imageBuf2[i][j]>t1) index[i][j]=4774.235939*pow(pi*39.9*39.9, -0.25)*exp(-(imageBuf2[i][j]-130.0)*(imageBuf2[i][j]-130.0)/(2*39.9*39.9));
		else if (imageBuf2[i][j]>=t2&&imageBuf2[i][j]<=t3) index[i][j]=469.964662*pow(pi*8.7*8.7, -0.25)*exp(-(imageBuf2[i][j]-207.0)*(imageBuf2[i][j]-207.0)/(2*8.7*8.7));
		else  index[i][j]=4774.235939*pow(pi*39.9*39.9, -0.25)*exp(-(imageBuf2[i][j]-130.0)*(imageBuf2[i][j]-130.0)/(2*39.9*39.9));
	}
		
for(i=0; i<height; i++)
	for(j=0; j<width; j++)
		if(imageBuf2[i][j]==0) index[i][j]=0.0;
		
for(i=0; i<height; i++)
	for(j=0; j<width; j++)
	{
		if(imageBuf2[i][j]==0)
		{
			n=0;
			for(k=-1; k<=1; k++)
			for(t=-1; t<=1; t++)
			{
				if((i+k)<0||(j+t)<0||(i+k)>=height||(j+t)>=width) continue;
				else  if(imageBuf2[i+k][j+t]!=0)
				{
					array[n]=imageBuf2[i+k][j+t];
					n++;
				}	
			}
			if(n!=0)
			{
				sort(array, n);
				if(n%2==0) imageBuf3[i][j]=(int)((array[(n-1)/2]+array[n/2])/2);	
				else imageBuf3[i][j]=array[(n-1)/2];
			}	
			else
			{
				sum=0.0;
				indexsum=0.0;
				for(k=-2; k<=2; k++)
				for(t=-2; t<=2; t++)
				{
					if((i+k)<0||(j+t)<0||(i+k)>=height||(j+t)>=width) continue;
					else
					{
						sum+=(double)imageBuf2[i+k][j+t]*index[i+k][j+t];
						indexsum+=index[i+k][j+t];
					}
				}		
				if(indexsum!=0.0)
					imageBuf3[i][j]=(int)(sum/indexsum+0.5);
				else imageBuf3[i][j]=imageBuf2[i][j];	
			}		
		}    
		else imageBuf3[i][j]=imageBuf2[i][j];
	}	
/*for(i=0; i<height; i++)
	for(j=0; j<width; j++)
		printf("%d\t", imageBuf3[i][j]);
*/	
	for(i=0; i<height; i++)
		for(j=0; j<width; j++)
			im_put_pix(imageBuf3[i][j], im, i, j);
/*-----------------------------add-start--------------------------*/
	Imrect *comp_orig = NULL;
	int pix;
	char filename[64];
	char string[25];
	gcvt(r,5,string);
	strcpy(filename,"result/test/nhwfm_");
	strcat(filename, string);
	strcat(filename,".bmp");
	comp_orig = ReadGIF("/home/wzz/newimscramdegree/result/Lena.gif",1);
	if (comp_orig == NULL)
	{
		printf("read image error.");
		return;
	}
	for(i = 0;i < height; i ++)
		for(j = 0;j < width; j++)
		{
			IM_PIX_GET(comp_orig, i, j, pix);
			imageBuf1[i][j] = pix;
		}
/*---------------------------------end---------------------------*/
	double mae = 0.0;
	sum=0.0;
	for(i=0; i<height; i++)
		for(j=0; j<width; j++)
		{
			sum+=(double)(pow((imageBuf1[i][j]-imageBuf3[i][j]), 2));
			indexsum+=(double)(imageBuf1[i][j]*imageBuf1[i][j]);	
			mae += fabs(imageBuf1[i][j]-imageBuf3[i][j]);
		}	
	mae = mae/(height*width);
printf("sum=%f\t, indexsum=%f\n", sum, indexsum);			
printf("NMSE=%f\n", sum/indexsum);			
indexsum=pow(256.0, 4);
printf("PSNR=%f MAE = %f\n", 10*log10(indexsum/sum),mae);
			
	ImrectToBMP8(im,filename);//add
	//mono_image_set(im);
	//stack_push(im, IMRECT, im_free);
	
	free_int_matrix(imageBuf1, height);
	free_int_matrix(imageBuf2, height);
	free_int_matrix(imageBuf3, height);
	free_db_matrix(index, height);
	free(hRow);
}