wfm.c

图像置乱代码

#include"include.h"
#include"tch.h"

void wfm_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;
	FILE *fp=NULL;
	int g1, g2, gmin, min;
	static double Cdk, Cmk, Cbk;
	static int s1, s2, s3;
	double u1, u2, u3, sigma1, sigma2, sigma3;
	int **imageBuf1, **imageBuf2, **imageBuf3;
	double t1, t2;
	double **index;
	double sum, indexsum;
	
	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(i=0; i<85; i++)
	{
		Cdk+=(double)(i*rechist[i]);
		s1+=rechist[i];
	}	
	Cdk=Cdk/(double)s1;
printf("\nCdk=%f\n", Cdk);
	for(i=85; i<170; i++)
	{
		Cmk+=(double)(i*rechist[i]);
		s2+=rechist[i];
	}
	Cmk=Cmk/(double)s2;
printf("Cmk=%f\n", Cmk);
	for(i=170; i<255; i++)
	{
		Cbk+=(double)(i*rechist[i]);
		s3+=rechist[i];
	}
	Cbk=Cbk/(double)s3;
printf("Cbk=%f\n", Cbk);

	u1=Cdk;
	sigma1=(u1-25.0)/3.0;
printf("sigma1=%f\n", sigma1);	
	u3=Cbk;
	sigma3=(230.0-u3)/3.0;
printf("sigma3=%f\n", sigma3);	
	u2=Cmk;
	
for(i=(int)(u1+0.5), min=rechist[(int)(u1+0.5)]; i<=(int)(u2+0.5); i++)
	if(rechist[i]<min)
	{
		min=rechist[i];
		gmin=i;
	}	
printf("gmin=%d\n", gmin);
sigma2=(u2-(double)gmin)/3.0;
printf("sigma2=%f\n", sigma2);

for(t1=57.0; t1<130.0; t1+=0.1)
if(exp(-(t1-57.0)*(t1-57.0)/(2*10.66*10.66))<=exp(-(t1-130.0)*(t1-130.0)/(2*20.62*20.62)))
{printf("t1=%f\n", t1); break;}

for(t2=130.0; t2<193.0; t2+=0.1)
if(exp(-(t2-130.0)*(t2-130.0)/(2*20.62*20.62))<=exp(-(t2-193.0)*(t2-193.0)/(2*12.29*12.29)))
{printf("t2=%f\n", t2); 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]=exp(-(imageBuf2[i][j]-56.99)*(imageBuf2[i][j]-56.99)/(2*10.66*10.66));
		else if(imageBuf2[i][j]<=t2&&imageBuf2[i][j]>t1) index[i][j]=exp(-(imageBuf2[i][j]-129.86)*(imageBuf2[i][j]-129.86)/(2*20.62*20.62));
		else index[i][j]=exp(-(imageBuf2[i][j]-193.13)*(imageBuf2[i][j]-193.13)/(2*12.29*12.29));
	}
		
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)
		{
			sum=0.0;
			indexsum=0.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  
				{
					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];
	}	
/*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);
	
	sum=0.0;
	for(i=0; i<height; i++)
		for(j=0; j<width; j++)
		{
			if(imageBuf3[i][j]!=0)
			sum+=(double)(pow((imageBuf1[i][j]-imageBuf3[i][j]), 2));
			indexsum+=(double)(imageBuf1[i][j]*imageBuf1[i][j]);
		}	
printf("sum=%f\t, NMSE=%f\n", sum, sum/indexsum);			
indexsum=pow(256.0, 4);
printf("PSNR=%f\n", 10*log10(indexsum/sum));
	
	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);
}