newscramdegree.c

图像置乱代码

		if(pp[i] > 0.000001)
		{
			printf("pp[%d] = %f \n",i,pp[i]);
			//fprintf(fp1,"pp[%d] = %f \n",i,pp[i]);
			fprintf(fp1," %d   %f \n",i,pp[i]);
		}
		pmean = pmean + i*pp[i];
	}
	printf("pmean = %f\n",pmean);
	fprintf(fp1,"pmean = %f\n",pmean);

	pdegree = 0;
	ped1 = 0;
	ped2 = 0;
	for(i = 0;i < pMASK_size;i ++)
	{
		pdegree += (double)fabs(i - pMEAN_n)*pp[i];
		if(i - pMEAN_n < 0.000001)
			ped1 +=  (double)(i - pMEAN_n)*pp[i];
		if(i - pMEAN_n > 0.000001)
			ped2 +=  (double)(i - pMEAN_n)*pp[i];
	}
	printf("pdegree = %f  ped1 = %f  ped2 = %f ped1+ped2 = %f fabs(ped1+ped2)/(pMASK_size-1- pMEAN_n) = %f\n",pdegree,ped1,ped2,ped1+ped2,fabs(ped1+ped2)/(pMASK_size-1- pMEAN_n));
	fprintf(fp1,"pdegree = %f  ped1 = %f  ped2 = %f ped1+ped2 = %f fabs(ped1+ped2)/(pMASK_size-1- pMEAN_n) = %f\n",pdegree,ped1,ped2,ped1+ped2,fabs(ped1+ped2)/(pMASK_size-1- pMEAN_n));

	var =0.0;
	for(i=0; i< pMASK_size; i++) 
		var += (i - pMEAN_n)*(i - pMEAN_n)*pp[i];
	//var /= (double)(pMASK_size-1);
	var = sqrt(var);
	Var = (double)(pMASK_size-1)*pp_m*(1-pp_m);
	//Var = sqrt(Var);
	Varp = pp_m*(1-pp_m)/(double)(pMASK_size-1);
	Varp = sqrt(Varp);
	printf("var = %f  Var*Var = %f Var = %f Varp = %f n*pp_m = %f n*(1-pp_m) = %f\n",var,Var,sqrt(Var),Varp,(double)(pMASK_size-1)*pp_m,(double)(pMASK_size-1)*(1-pp_m));
	fprintf(fp1,"var = %f  Var*Var = %f Var = %f Varp = %f n*pp_m = %f n*(1-pp_m) = %f\n",var,Var,sqrt(Var),Varp,(double)(pMASK_size-1)*pp_m,(double)(pMASK_size-1)*(1-pp_m));
	printf("\n");
	fprintf(fp1,"\n");
}
fclose(fp1);
/*----------------------------------------------------------------*/


	double var,Var,Varp;
	var =0.0;
	for(i=0; i< MASK_size; i++) 
		var += (i - MEAN_n)*(i - MEAN_n)*p[i];
	var = sqrt(var);
	Var = (double)(MASK_size-1)*p_m*(1-p_m);
	//Var = sqrt(Var);
	Varp = p_m*(1-p_m)/(double)(MASK_size-1);
	Varp = sqrt(Varp);
	printf("var = %f  Var*Var = %f Var = %f Varp = %f n*p_m = %f n*(1-p_m) = %f\n",var,Var,sqrt(Var),Varp,(double)(MASK_size-1)*p_m,(double)(MASK_size-1)*(1-p_m));

/*----------------------------------------------------------------*/

/*----------------------------------------------------------------
	double meanpc;
	meanpc = 0.0;
	for(i = (int)(ed1+MEAN_n)+1;i < ed2+MEAN_n;i ++)
		meanpc += p[i];
	printf("meanpc = %f\n",meanpc);
/*----------------------------------------------------------------*/
//+===================================================

	int k; double DSF, GSF, DDR, Ds;
	/*FILE * fp1 = fopen("zhiluandu.txt", "wb");*/
/*	if (fp1==NULL)
	{
		error("file open error ", warning );
		return;
	}
*/
//+
	gcvt(n_scram,5,string);
	strcat(filename,string);
	strcpy(filename2, filename);
	strcat(filename2,"_times012.txt");
	//fp3 = fopen(filename2,"wb");
	//strcat(filename3,"_detail.txt");
	strcpy(filename3,filename);
	strcat(filename3,"_mean_n.txt");
	fp4 = fopen(filename3,"wb");
	strcpy(filename4,filename);
	strcat(filename4,"_mean_nn.txt");
	//fp2 = fopen(filename4,"wb");
/*------------------------------------------------------------------------------
	printf("原始图矩阵:\n");
		for(i = 0; i < height; i++)
		{
			for(j = 0; j < width; j++)
			{
				printf("(%d %d) ",dest_r[i][j],dest_c[i][j]);
			}
			printf("\n");
		}
	printf("用于每次置乱的置乱矩阵:\n");
		for(i = 0; i < height; i++)
		{
			for(j = 0; j < width; j++)
			{
				r = trans_r[i][j];
				c = trans_c[i][j];
				src_r[r][c] = i;
				src_c[r][c] = j;
			}
		}
		for(i = 0; i < height; i++)
		{
			for(j = 0; j < width; j++)
			{
				printf("(%d %d) ",src_r[i][j],src_c[i][j]);
			}
			printf("\n");
		}
/*-----------------------------------------------------------------------------*/
//+
	strcpy(fnds,filename);
	strcat(fnds,"_DS");
	strcpy(fndstxt,fnds);
	strcat(fndstxt,".txt");
	//fp = fopen(fndstxt,"wb");
	char fnhistg[64];

	for(k = 0; k < n_scram; k++)//n_scram为置乱次数
	{
//+------------------------------------------------------
/*
	gcvt(k+1,5,string);
	strcpy(filename1, filename);
	strcat(filename1,string);
	strcat(filename1, ".txt");
	fp2 = fopen(filename1,"wb");
*/

		//printf("\n------------------%d次置乱------------------\n",k+1);
		//fprintf(fp4,"\n------------------%d次置乱------------------\n",k+1);
                 //计数器初始化
		for(i = 0; i < height; i ++)
			for(j = 0; j < width ; j ++)
			{
				count[i][j] = 0;
				countn[i][j] = 0;
			}
		for(i = 0; i < MASK_size; i ++)
		{
			count_histg[i] = 0;
			count_histgn[i] = 0;
		}
//+------------------------------------------------------
		for(i = 0; i < height; i++)
		{
			for(j = 0; j < width; j++)
			{
				tmp_r[i][j] = dest_r[ trans_r[i][j] ][ trans_c[i][j] ];
				tmp_c[i][j] = dest_c[ trans_r[i][j] ][ trans_c[i][j] ];
			}
		}

		for(i = 0; i < height; i++)
		{
			for(j = 0; j < width; j++)
			{
				dest_r[i][j] = tmp_r[i][j];
				dest_c[i][j] = tmp_c[i][j];
				IM_PIX_GET(srcIm, i, j, pix);
				IM_PIX_SET(destIm, dest_r[i][j], dest_c[i][j], pix);
			}
		}
//+=========================================================================

		for(i = 0; i < height; i++)
		{
			for(j = 0; j < width; j++)
			{
				r = dest_r[i][j];
				c = dest_c[i][j];
				src_r[r][c] = i;
				src_c[r][c] = j;
			}
		}
/*
		for(i = 0; i < height; i++)
		{
			for(j = 0; j < width; j++)
			{
				printf("(%d %d) ",src_r[i][j],src_c[i][j]);
			}
			printf("\n");
		}
*/
//fp1 = fopen("points.txt","wb");
/*-------------------------------------------------------------------------*/
	SUM = (width - 2*mask_half)*(height - 2*mask_half);//参与计算的pixel总数
		for(i =  mask_half; i < height - mask_half; i++)
			for(j =  mask_half; j < width - mask_half; j++) 
			{		
		        	X0 = src_r[i][j];
		        	Y0 = src_c[i][j];
				for(mask_row = i - mask_half; mask_row <= i + mask_half ; mask_row++)
				{	
				
					for(mask_col = j - mask_half; mask_col <= j + mask_half ; mask_col++)
					{
						if(mask_row == i && mask_col == j)
							continue;
						Xp =src_r[mask_row][mask_col];
						Yp =src_c[mask_row][mask_col];
						if((abs(X0 - Xp) <= mask_half) && (abs (Y0 - Yp) <= mask_half)) 
						{
							count[i][j]++;//统计(i,j)这个点的老邻居个数
						}
					}
				}
			}
		for(i = mask_half ; i < height - mask_half ; i ++)
			for(j = mask_half ; j < width - mask_half ; j ++)
			{  
				s = count[i][j];
				count_histg[s] ++; //统计有s个老邻居的pexel个数
				//fprintf(fp1,"%d    %d \n",i*height+j,s);
			}
/*----------------------------------------------------------------*/
//fclose(fp1);
/*----------------------------------------------------------------*/
	mean_n = 0.0;
	//mean_nn = 0.0;

	
		for(i = 0; i < MASK_size; i ++)
		{
			rp[i] = (double)count_histg[i] / SUM ;
			if(count_histg[i] != 0)
			{
				//printf("有 %d 个老邻居个数为 %d    有 %d 个老邻居概率为 %f \n",i,count_histg[i],i,(float)count_histg[i] / SUM);
				//fprintf(fp4,"有 %d 个老邻居个数为 %d    有 %d 个老邻居概率为 %f \n",i,count_histg[i],i,(float)count_histg[i] / SUM);
				//fprintf(fp2,"%d  %f\n",i,(float)count_histg[i] / SUM);
				mean_n += i * rp[i];
				//mean_nn += i * (float)count_histgn[i] / SUMn ;
			}
		}
//+========================================================================
/*----------------------------------------------------------------
		printf("%d  %f  %f  %f  %f  %f\n",k+1,(float)count_histgn[0] / SUMn, (float)count_histgn[1] / SUMn, (float)count_histgn[2] / SUMn,(float)count_histgn[3] / SUMn,(float)count_histgn[4] / SUMn);
		fprintf(fp2,"%d  %f  %f  %f  %f  %f\n",k+1,(float)count_histgn[0] / SUMn, (float)count_histgn[1] / SUMn, (float)count_histgn[2] / SUMn,(float)count_histgn[3] / SUMn,(float)count_histgn[4] / SUMn);
/*----------------------------------------------------------------*/
		printf("%d  %f \n",k+1,mean_n);
		//fprintf(fp4,"%d  %f ",k+1,mean_n);
/*----------------------------------------------------------------*/
	double rdegree = 0, red1 = 0, red2 = 0;
	double rvar = 0;
	double ED ;
	for(i = 0;i < MASK_size;i ++)
	{
		rvar += (i-mean_n)*(i-mean_n)*rp[i];
		rdegree += fabs(i - MEAN_n)*rp[i];
		if(i - MEAN_n < 0.000001)
			red1 +=  (i - MEAN_n)*rp[i];
		if(i - MEAN_n > 0.000001)
			red2 +=  (i - MEAN_n)*rp[i];
	}
	ED =  fabs(red1+red2)/(MASK_size-1 - MEAN_n);
	printf("ED = %f\nrvar = %f\n",ED,rvar);
/*---------------------------------------------------------------*/
	//fprintf(fp4," %f  %f \n",fabs(red1+red2),fabs(red1+red2)/(MASK_size-1 - MEAN_n));
	//printf("rdegree = %f  red1 = %f  red2 = %f red1+red2 = %f fabs(red1+red2)/ed2 = %f\n",rdegree,red1,red2,red1+red2,fabs(red1+red2)/ed2);
	//fprintf(fp4," %f %f  \n",fabs(red1+red2),fabs(red1+red2)/ed2);
	//fprintf(fp4," %f \n",mean_nn);	
/*-------------------------------------------------------------------*/
	double diff = 0.0;
	double similarity;
	double dot_product=0.0;
	double model_p=0.0,model_rp=0.0;
	for(i = 0;i < MASK_size;i ++)
	{
		diff += fabs(rp[i] - p[i]);
		model_p += p[i]*p[i];
		model_rp += rp[i]*rp[i];
		dot_product += p[i]*rp[i];
	}
	similarity = dot_product/sqrt(model_p*model_rp);
	printf("diff = %f    similarity = %f\n",diff/2,similarity);
/*-----------------------------------------------------------------*/

	fprintf(fp4,"%d  %f ",k+1,mean_n);
	fprintf(fp4," %f  %f  %f  %f  %f  %f\n",fabs(red1+red2),ED,diff/2,similarity,(ED+diff/2)/2,1-(ED+diff/2)/2);
/*-----------------------------------------------------------------*/
	double entropy = 0.0;
	double r_entropy = 0.0;
	for(i = 0;i < MASK_size;i ++)
	{
		if(p[i]>0)
			entropy += p[i]*log2(p[i]);
		if(rp[i]>0)
			r_entropy += rp[i]*log2(rp[i]);
	}
	entropy = - entropy;
	r_entropy = - r_entropy;
	printf("entropy = %f   r_entropy = %f \n",entropy,r_entropy);
/*----------------------------------------------------------------
	double rmeanpc;
	rmeanpc = 0.0;
	for(i = (int)(red1+MEAN_n)+1;i < red2+MEAN_n;i ++)
		rmeanpc += (double)count_histg[i]/SUM;
	printf("rmeanpc = %f\n",rmeanpc);
/*----------------------------------------------------------------*/
/*
	DSF = getDSF(dest_r, dest_c, height, width);
	GSF = getGSF(srcIm, destIm, height, width);
	DDR = getDDR(srcIm, destIm, height, width);
	Ds = DSF * GSF * DDR;
*/