denoise_3sigma.c

图像置乱代码

/********************************************************************
 *
 * image repaired tool
 *
 ********************************************************************/
#include <tina/all_tina.h>
#include "bmp.h"
static int iminarg1,iminarg2;
#define IMIN(a,b) (iminarg1 = (a),iminarg2 = (b),(iminarg1) < (iminarg2) ? (iminarg1) : (iminarg2))


#define nMAX 256*256
static double gdata[nMAX];

static Tv *tv = NULL;
static Tv *tv1=NULL;

#define MAXPATHLEN 64
static char directory_name[MAXPATHLEN];
static char file_name[MAXPATHLEN];
static char pathname[MAXPATHLEN];
static void *pdir=0, *pfile=0;

static 	Imrect *read_image(void)
{
	Imrect *srcIm;
/*----------------Read the image---------------------------*/
  	(void) strip_spaces(file_name);
  	(void) strip_spaces(directory_name);
  	(void) string_append(pathname, directory_name, "/", file_name, NULL);
	srcIm = BmpToImrect(pathname);
/*------------------------------------------------*/
	return(srcIm);
}
/*
static void denoise_3sigma_proc(void)
{
	
}
*/


double get_mean(data,n) 
double data[];	int n;
{
	//printf("get_mean   %d\n",n);
	int i;	double mean=0.0;
	for(i=0; i<n; i++)  mean += data[i];
	return (mean /= (double)n);
}

double get_var(data,n)
double data[];   int n;
{
	//printf("get_var   %d\n",n);
	int i;  
	double ave,var;
	ave = get_mean(data,n);
	var = 0.0;
	for(i=0; i<n; i++) 
		var += (data[i]-ave)*(data[i]-ave);
	return (var /= (double)(n-1));

/*	double ep,s;
	var = ep = 0.0;
	for (i = 0;i < n;i ++)
	{
		s = data[i] - ave;
		ep += s;
		var += s * s;
	}
	return (var - ep * ep/n)/(n-1);
*/
}

/*-------ascending numerical order-------*/
static void sorting(data,n)
double data[];	int n;
{
	int i, j, min;	
	double dmin;

	for(i=0; i<n-1; i++)
	{
		min = i;	
		dmin = data[i];
		for(j=i+1; j<n; j++)
			if(data[j]<dmin)	
			{
				min = j;	
				dmin = data[j];
			}
		data[min] = data[i];	
		data[i] = dmin;	

	}
}

double get_med(data,n)
double data[];   int n;
{
	int i; 
/*-----------------------------------*/
	double sorted[n];
	//printf("%d \n ",n);
	for(i = 0;i < n;i ++)
	{
		sorted[i] = data[i];
	}
	sorting(sorted,n);
	if(fmod(n,2) == 1)
		return sorted[((n-1)/2)];
	else 
		return ((sorted[n/2] + sorted[n/2-1])/2);
/*-----------------------------------
	sorting(data,n);
	printf("sorting done!\n");
	if(fmod(n,2) == 1)
		return data[((n-1)/2)];
	else 
		return ((data[n/2] + data[n/2-1])/2);
/*-----------------------------------*/
}

/*-------------------------------------------------------------------*/
static void Get_ave_var_proc()
{

	Imrect *srcIm, *destIm;
	int width, height;
	int type;
	BYTE pix;

	int lx, ux, ly, uy;
	Imregion *roi;
	if (stack_check_types(IMRECT, NULL) == false)
	{
		error("Warning : wrong types on stack", warning);
		return;
	}
	srcIm = (Imrect *) stack_pop(&type);
	width = srcIm->width;
    	height = srcIm->height;

	roi = srcIm->region;
	if (roi == NULL) return;
	lx = roi->lx;
	ux = roi->ux;
	ly = roi->ly;
	uy = roi->uy;
	int i,j;
	for (i = 0; i < height  ; i++)
	{
		for (j = 0; j< width ; j++)
		{
			IM_PIX_GET(srcIm, i, j, pix);
			gdata[256*i+j] = pix;
			//printf("gdata[%d]  =  %f\n",256*i+j,gdata[256*i+j]);
		}
	}
	double ave,var,med;
	ave = var = med =0;
	ave = get_mean(gdata,nMAX);
	var = get_var(gdata,nMAX);
	med = get_med(gdata,nMAX);
	printf("ave = %f    var = %f  med = %f \n",ave,var,med);
}
/*-------------------------------------------------------------------*/
#define NR  	20
#define R0  	2097151
#define LAMBDA  2045
#define MOD 	2097152

double gauss_random(i, mean,sigma)
int i;		/* set to 0 for the first random number, non-zero for all others */
double mean, sigma;
{
    double uniform_random(), sum;   int j;
 
    if(i==0)  return(mean + sigma * uniform_random(-1.0, 1.0, R0, LAMBDA, MOD));
	for(sum=0.0, j=0; j<NR; j++)
		sum += uniform_random(-1.0, 1.0, 0, LAMBDA, MOD);
	return ((double)(mean + sigma * sum*sqrt(3.0/(double)NR)) );
}
 
/*-------------------------------------------------------------------*/
double uniform_random(a, b, ar0, lamda, mod)
double a, b;
int ar0, lamda, mod;
{
    static  double x;
    long    entier;

    if(ar0 != 0) x = (double)ar0/(double)mod;
    x *= (double)lamda;
    entier = (x > 0.0)?		floor(x)  :  ceil(x);
    x -= (double)entier;
    return(x*(b-a)+a);
}
 /*---------------------------------------------------------------------------
	show an example 
	use median to estimate the sigma when there are some outline points.
 ---------------------------------------------------------------------------*/
static void Gauss_proc()
{
	#define G_N 5000
	double gauss1[G_N+500];
	double gauss4[G_N+500];
	double gaussG_N_1[G_N];
	double gaussG_N_4[G_N];
	double fabs_gauss1[G_N+500];
	double fabs_gauss4[G_N+500];
	double fabs_gaussG_N_1[G_N];
	double fabs_gaussG_N_4[G_N];
	int i;
	int Totol = G_N+500;
	int TotolG_N = G_N;

	i = 0;
	gauss1[i] = fabs(gauss_random(0, 0.0, sqrt(1.0)));
	fabs_gauss1[i] = fabs(gauss1[i]);
	gaussG_N_1[i] = gauss1[i];
	fabs_gaussG_N_1[i] = fabs(gaussG_N_1[i]);
	for (i = 1 ;i < G_N ;i ++)
	{
		gauss1[i] = gauss_random(1, 0.0, sqrt(1.0));
		fabs_gauss1[i] = fabs(gauss1[i]);
		gaussG_N_1[i] = gauss1[i];
		fabs_gaussG_N_1[i] = fabs(gaussG_N_1[i]);
	}
	for (i = G_N ;i < G_N+500 ;i ++)
	{
		gauss1[i] = gauss_random(1, 3.1, sqrt(0.5));
		fabs_gauss1[i] = fabs(gauss1[i]);
	}

	double aveG_N_1,varG_N_1,medG_N_1,sdevG_N_1,fabs_medG_N_1;
	aveG_N_1 = varG_N_1 = medG_N_1 = sdevG_N_1 = 0.0;
	aveG_N_1 = get_mean(gaussG_N_1,TotolG_N);
	varG_N_1 = get_var(gaussG_N_1,TotolG_N);
	sdevG_N_1 = sqrt(varG_N_1);
	medG_N_1 = get_med(gaussG_N_1,TotolG_N);
	fabs_medG_N_1 = get_med(fabs_gaussG_N_1,TotolG_N);
	printf("(0,1)->5000 \n aveG_N_1 = %f  sdevG_N_1 = %f   medG_N_1 = %f   fabs_medG_N_1 = %f\n\n",aveG_N_1,sdevG_N_1,medG_N_1,fabs_medG_N_1);

	double ave1,var1,med1,sdev1,fabs_med1;
	ave1 = var1 = med1 = sdev1 = 0.0;
	ave1 = get_mean(gauss1,Totol);
	var1 = get_var(gauss1,Totol);
	sdev1 = sqrt(var1);
	med1 = get_med(gauss1,Totol);
	fabs_med1 = get_med(fabs_gauss1,Totol);
	printf("(0,1)->5000 + (3.1,0.5)->500\n ave1 = %f       sdev1 = %f       med1 = %f        fabs_med1 = %f \n\n",ave1,sdev1,med1,fabs_med1);



/////////////////////////////
	//i = 0;
	//gauss4[i] = fabs(gauss_random(0, 0.0, 4.0));
	//gaussG_N_4[i] = gauss4[i];
	for (i = 0;i < G_N;i ++)
	{
		gauss4[i] = gauss_random(1, 0.0, sqrt(4.0));
		fabs_gauss4[i] = fabs(gauss4[i]);
		gaussG_N_4[i] = gauss4[i];
		fabs_gaussG_N_4[i] = fabs(gaussG_N_4[i]);
	}
	for (i = G_N ;i < G_N+500 ;i ++)
	{
		gauss4[i] = gauss_random(1, 6.1, sqrt(0.5));
		fabs_gauss4[i] = fabs(gauss4[i]);
	}

	double aveG_N_4,varG_N_4,medG_N_4,sdevG_N_4,fabs_medG_N_4;
	aveG_N_4 = varG_N_4 = medG_N_4 = sdevG_N_4 = 0.0;
	aveG_N_4 = get_mean(gaussG_N_4,TotolG_N);
	varG_N_4 = get_var(gaussG_N_4,TotolG_N);
	sdevG_N_4 = sqrt(varG_N_4);
	medG_N_4 = get_med(gaussG_N_4,TotolG_N);
	fabs_medG_N_4 = get_med(fabs_gaussG_N_4,TotolG_N);
	printf("(0,4)->5000 \n aveG_N_4 = %f   sdevG_N_4 = %f   medG_N_4 = %f    fabs_medG_N_4 = %f\n\n",aveG_N_4,sdevG_N_4,medG_N_4,fabs_medG_N_4);

	double ave4,var4,med4,sdev4,fabs_med4;
	ave4 = var4 = med4 = sdev4 = 0.0;
	ave4 = get_mean(gauss4,Totol);
	var4 = get_var(gauss4,Totol);
	sdev4 = sqrt(var4);
	med4 = get_med(gauss4,Totol);
	fabs_med4 = get_med(fabs_gauss4,Totol);
	printf("(0,4)->5000 + (6.1,0.5)->500 \n ave4 = %f       sdev4 = %f       med4 = %f        fabs_med4 = %f\n\n",ave4,sdev4,med4,fabs_med4);
	

	#define Area  36
	int gauss_histg1[Area];
	int gauss_histg4[Area];
	int gauss_histgG_N_1[Area];
	int gauss_histgG_N_4[Area];
	int rn;
	for (i = 0; i < Area; i++)
		gauss_histg1[i] = 0;
	for (i = 0; i < Totol; i ++)
	{
		rn = floor(gauss1[i]);
		gauss_histg1[rn+15]++;
	}

	for (i = 0; i < Area; i++)
		gauss_histg4[i] = 0;
	for (i = 0; i < Totol; i ++)
	{
		rn = floor(gauss4[i]);
		gauss_histg4[rn+15]++;
	}

	for (i = 0; i < Area; i++)
		gauss_histgG_N_1[i] = 0;
	for (i = 0; i < TotolG_N; i ++)
	{
		rn = floor(gauss1[i]);
		gauss_histgG_N_1[rn+15]++;
	}

	for (i = 0; i < Area; i++)
		gauss_histgG_N_4[i] = 0;
	for (i = 0; i < TotolG_N; i ++)
	{
		rn = floor(gauss4[i]);
		gauss_histgG_N_4[rn+15]++;
	}
//////////
	FILE *fp;
     	if((fp=fopen("gauss/gaussG_N_1.txt","w"))==NULL)
    	{
         	format("Can't open the file!\n");
         	exit(0);
    	} 
	for (i = 0;i < TotolG_N;i ++)
		fprintf(fp,"%d  %f\n",i,gaussG_N_1[i]);
     	fclose(fp);
///
     	if((fp=fopen("gauss/gauss1.txt","w"))==NULL)
    	{
         	format("Can't open the file!\n");
         	exit(0);
    	} 
	for (i = 0;i < Totol;i ++)
		fprintf(fp,"%d  %f\n",i,gauss1[i]);
     	fclose(fp);
//////////
     	if((fp=fopen("gauss/gaussG_N_4.txt","w"))==NULL)
    	{
         	format("Can't open the file!\n");
         	exit(0);
    	} 
	for (i = 0;i < TotolG_N;i ++)
		fprintf(fp,"%d  %f\n",i,gaussG_N_4[i]);
     	fclose(fp);
////
     	if((fp=fopen("gauss/gauss4.txt","w"))==NULL)
    	{
         	format("Can't open the file!\n");
         	exit(0);
    	} 
	for (i = 0;i < Totol;i ++)
		fprintf(fp,"%d  %f\n",i,gauss4[i]);
     	fclose(fp);
//////////
     	if((fp=fopen("gauss/gauss_histg1.txt","w"))==NULL)
    	{
         	format("Can't open the file!\n");
         	exit(0);
    	} 
	for (i = 0;i < Area;i ++)
		fprintf(fp,"%d  %d\n",i-15,gauss_histg1[i]);
     	fclose(fp);
//////////
     	if((fp=fopen("gauss/gauss_histg4.txt","w"))==NULL)
    	{
         	format("Can't open the file!\n");
         	exit(0);
    	} 
	for (i = 0;i < Area;i ++)
		fprintf(fp,"%d  %d\n",i-15,gauss_histg4[i]);
     	fclose(fp);

     	if((fp=fopen("gauss/gauss_histgG_N_1.txt","w"))==NULL)
    	{
         	format("Can't open the file!\n");
         	exit(0);
    	} 
	for (i = 0;i < Area;i ++)
		fprintf(fp,"%d  %d\n",i-15,gauss_histgG_N_1[i]);
     	fclose(fp);

     	if((fp=fopen("gauss/gauss_histgG_N_4.txt","w"))==NULL)
    	{
         	format("Can't open the file!\n");
         	exit(0);
    	} 
	for (i = 0;i < Area;i ++)
		fprintf(fp,"%d  %d\n",i-15,gauss_histgG_N_4[i]);
     	fclose(fp);
	
	printf("-----ok!-----\n");
}

void four1(data, nn, isign)
float data[];
int nn, isign;
{
    int n, mmax, m, j, istep, i;
    double wtemp, wr, wpr, wpi, wi, theta;
    float tempr, tempi;
    
    n = nn << 1;
    j = 1;
    for (i = 1; i < n; i += 2) {
	if (j > i) {
	    tempr = data[j];     data[j] = data[i];     data[i] = tempr;
	    tempr = data[j+1]; data[j+1] = data[i+1]; data[i+1] = tempr;
	}
	m = n >> 1;
	while (m >= 2 && j > m) {
	    j -= m;
	    m >>= 1;
	}
	j += m;
    }
    mmax = 2;
    while (n > mmax) {
	istep = 2*mmax;
	theta = TWOPI/(isign*mmax);
	wtemp = sin(0.5*theta);
	wpr = -2.0*wtemp*wtemp;
	wpi = sin(theta);
	wr = 1.0;
	wi = 0.0;
	for (m = 1; m < mmax; m += 2) {
	    for (i = m; i <= n; i += istep) {
		j =i + mmax;
		tempr = wr*data[j]   - wi*data[j+1];
		tempi = wr*data[j+1] + wi*data[j];
		data[j]   = data[i]   - tempr;
		data[j+1] = data[i+1] - tempi;
		data[i] += tempr;
		data[i+1] += tempi;
	    }
	    wr = (wtemp = wr)*wpr - wi*wpi + wr;
	    wi = wi*wpr + wtemp*wpi + wi;
	}
	mmax = istep;
    }
}

void realft(data,n,isign)
float data[];
int n,isign;
{
    int i, i1, i2, i3, i4, n2p3;
    float c1 = 0.5, c2, h1r, h1i, h2r, h2i;
    double wr, wi, wpr, wpi, wtemp, theta;
    void four1();

    theta = 3.141592653589793/(double) (n>>1);
    if (isign == 1) {
	c2 = -0.5;
	four1(data, n>>1, 1);
    } 
    else {
	c2 = 0.5;
	theta = -theta;
    }
    wtemp = sin(0.5*theta);
    wpr = -2.0*wtemp*wtemp;
    wpi = sin(theta);
    wr = 1.0+wpr;
    wi = wpi;
    n2p3 = 2*n+3;
    for (i = 2; i <= (n>>2); i++) {
	i4 = 1 + (i3 = n2p3 - (i2 = 1 + ( i1 = i + i - 1)));
	h1r =  c1*(data[i1] + data[i3]);
	h1i =  c1*(data[i2] - data[i4]);
	h2r = -c2*(data[i2] + data[i4]);
	h2i =  c2*(data[i1] - data[i3]);
	data[i1] =  h1r + wr*h2r - wi*h2i;
	data[i2] =  h1i + wr*h2i + wi*h2r;
	data[i3] =  h1r - wr*h2r + wi*h2i;
	data[i4] = -h1i + wr*h2i + wi*h2r;
	wr = (wtemp = wr)*wpr - wi*wpi+wr;
	wi = wi*wpr + wtemp*wpi + wi;
    }
    if (isign == 1) {
	data[1] = (h1r = data[1]) + data[2];
	data[2] = h1r - data[2];
    } else {
	data[1] = c1*((h1r = data[1]) + data[2]);
	data[2] = c1*(h1r - data[2]);
	four1(data, n>>1, -1);
    }
}