/*

	img8bit.c 11/22/91

	from -- reduce.c  3/3/89

	cl /AL img8bit.c

	reduce image1 /c 64 /m 10 -- sets max colors from 255 to 64
   	                          and min hits to use to 10


	last update 11/25/91

	by Russell A. Ambroziak, Ph. D.

	U. S. Geological Survey
	Office of Energy and Marine Geology
	National Center , STOP 915
	Reston, VA 22092

	Tel:  703-648-6168
	Fax:  703-648-5464


*/

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#define MAX_X 1024
#define MAX_Y 1024
/* Some constants used by RAM32 board */
/* Macro to turn on/off concurrent mode */

#define R      6367.65
#define PI		3.1415927
#define DEG_RAD 0.01745329251				/* converts degrees to radians */
#define KM_DEG 111.136458114
#define MID_LAT 26.5
#define MAXCOL 8191
#define RED 36
#define GRN 44
#define BLU 27
#define MAXLEVEL 256
#define TIMES 3

int VGA_Color[16][3]=
{
	{105, 90, 75},		/* 240 BUFF */
	{  0,  0,  0},		/* 241 BLACK */
	{ 70, 70, 70},		/* 242 D_GREY */
	{127,127,127},		/* 243 M_GREY */
	{191,191,191},		/* 244 L_GREY */
	{255,255,255},		/* 245 WHITE */
	{255,  0,  0},		/* 246 RED */
	{255,128,  0},		/* 247 ORANGE */
	{191,192,  0},		/* 248 YELLOW */
	{  0,191,  0},		/* 249 GREEN */
	{  0,191,191},		/* 250 CYAN */
	{  0,  0,255},		/* 251 BLUE */
	{191,  0,191},		/* 252 MAGENTA */
	{ 96, 28, 14},		/* 253 BROWN */
	{  0, 91, 91},		/* 254 D_CYAN */
	{  0,128,  0}		/* 255 D_GREEN */
};
int Num_Color,Min=1;
double Lat0,Lon0;
double X0,Y0,Xs,Ys,Scale;
float MaxLon,MinLon,MaxLat,MinLat;
unsigned char Inbuff[BUFSIZ],Inbuff2[BUFSIZ],Outbuff[BUFSIZ];
char InFileName[100],OutFileName[100];
FILE *FpIn;
int Row,Col;
unsigned char ImgBuff[4][1024];
int WordBuff[1024];
long int offset;
int Dval[8][2]=
{
	{ 0,-1},
	{ 1,-1},
	{ 1, 0},
	{ 1, 1},
	{ 0, 1},
	{-1, 1},
	{-1, 0},
	{-1,-1}
};


char Color_Name[3][6]={"red  ","green","blue "};
double AB[3][2]=
{
	{ 3.48007, -0.00897 },
	{ 3.31476, -0.01058 },
	{ 3.52924,  0.00050 }
};
unsigned char RGB_Div[3][MAXLEVEL],RGB_Mid[3][GRN];
int Rgb_Div[3][MAXLEVEL];
int Del[3];

int Lut[3][256];
float Flut[4][256];
int huge Hist[RED][GRN][BLU];
int huge Hist2[RED][GRN][BLU];
unsigned char huge Table[RED][GRN][BLU];
unsigned char Buffer[4][MAXCOL];
char Dither[MAXCOL*4];
double Beta[3][4]=
{
	{ 2.437, -0.00901,  0.00748,  0.00046},
	{ 2.531,  0.00437, -0.01062,  0.00173},
	{ 2.751,  0.00149,  0.00459, -0.00068}
};
int Max_Color=255;
char Name[100];

float rnum();
void limit_area();
int linear1();
int var_hist();
void make_lut();
void fill_table();
void make_img();
double check_sd();
void average_lut();
void make_div();
void measure_div();
void really_dumb();
int make_dither();
void output();

int boot_nine();
int clear_screen();
void UseCursor();


main(argc,argv)

int argc;
char *argv[];
{
	int i,j,k,a,b,c;
	int xc1=0,yc1=0,xs=512,ys=512,speed=32,maxbox=1024;
	int xc2=0,yc2=0;
	int xct=0,yct=0,xst=64,yst=64;
	char ans,filename[40];
	int min=1,numcol,num,seed=35,iop;
	FILE *fpbout;
	float sd=5.0,fx;

	srand(seed);
	if(argc<=1)
		print_help();
	if(argc>1)
		strcpy(InFileName,argv[1]);
	else
	{
		printf("Give input file name of '.cmg' file.\n\n");
		scanf("%s",InFileName);
	}
	for(k=0;k<strlen(InFileName);k++)
		if(InFileName[k]=='.')
			InFileName[k]='\0';
	strcat(InFileName,".cmg");
	FpIn=fopen(InFileName,"rb");
	if(!FpIn)
	{
		printf("Could not open file '%s'\n\n",InFileName);
		exit(0);
	}
/* Get out file name */
	if(argc>2)
		strcpy(OutFileName,argv[2]);
	else
	{
		printf("Give output file name.\n\n");
		scanf("%s",OutFileName);
		for(i=0;i<strlen(OutFileName);i++)
			if(OutFileName[i]=='.')
				OutFileName[i]='\0';
		for(i=0;i<strlen(OutFileName);i++)
			if(OutFileName[i]=='*'||OutFileName[i]=='?')
			{
				printf("Output file name is not legal -- '%s'\n",OutFileName);
				exit(0);
			}
		if((OutFileName[0]>='a'&&OutFileName[0]<='z')||
			(OutFileName[0]>='A'&&OutFileName[0]<='Z')||
			(OutFileName[0]>='0'&&OutFileName[0]<='9'))
		{
			printf("Output file name = '%s'\n",OutFileName);
		}
		else
		{
			printf("Output file name is not legal -- '%s'\n",OutFileName);
			exit(0);
		}
	}
	if(argc>3)
	{
		if(sscanf(argv[3],"%f",&fx)<1)
		{
			printf("'%s' is an illegal value for dither constant.\n",argv[3]);
			printf("No change made -- dither constant = %5.2f\n",sd);
		}
		else if(fx<0.0||fx>16.0)
		{
			printf("'%5.2f' is an illegal value for dither constant.\n",fx);
			printf("No change made -- dither constant = %5.2f\n",sd);
		}
		else
		{
			printf("dither constant changed from %5.2f to %5.2f\n",sd,fx);
			sd=fx;
		}
	}
	else
		printf("dither constant = %5.2f\n",sd);
	for(i=2;i<argc;i++)
	{
		if(argv[i][0]=='/'&&argv[i][1]=='c')
		{
			sscanf(argv[i+1],"%d",&Max_Color);
			printf("Max Colors reset from 255 to %d.\n\n",Max_Color);
		}
		if(argv[i][0]=='/'&&argv[i][1]=='m')
		{
			sscanf(argv[i+1],"%d",&Min);
			if(Min<1)
				Min=1;
			printf("Min points to keep set from 1 to %d.\n\n",Min);
			min=Min;
		}
	}
	printf("Processing 24 bit image '%s' into 8 bit image '%s'.\n",
		InFileName,OutFileName);
/* Add Dither */
	make_dither(FpIn,sd);
	fcloseall();
	FpIn=fopen("dither.cmg","rb");
	fread((char *)&Row,sizeof(int),1,FpIn);
	fread((char *)&Col,sizeof(int),1,FpIn);
	printf("There are %d rows and %d columns.\n",Row,Col);
	printf("Scaling Red, Green and Blue to fit image colors.\n");
	if(Row<1||Col<1)
		exit(0);
/*

	else
	{
		printf("You must give a file name and may give options:\n\n");
		printf("reduce image1 /c 64 /m 10 -- sets max colors from 255 to 64\n");
		printf("                             and min hits to use to 10\n\n");
		exit(0);	
	}

*/
	numcol=linear1(Buffer,Hist,Lut,Table,min);
	Num_Color=numcol;
	var_hist(Buffer,Hist2,1);
	printf("Now filling color assignment Table.\n");
	fill_table(Table,Hist,Hist2,min);
	printf("\nAveraging lookup table colors:\n");
	average_lut(Lut,Flut,Table,Buffer);
	if(Max_Color<=240)
	{
		for(k=0;k<16;k++)
		{
			Lut[0][k+240]=VGA_Color[k][0];
			Lut[1][k+240]=VGA_Color[k][1];
			Lut[2][k+240]=VGA_Color[k][2];
		}
	}
	printf("\nCreating image file\n");
	make_img(Table,Lut,Buffer,1);
	unlink("dither.cmg");
}

/***********************************************************************
**
**
**
**
************************************************************************/

int print_help()

{
printf("Format:\n\n");
printf(
"img8bit [infile] [outfile] <dither> </m min points> </c max colors>\n\n");
printf("Where:   infile -- name of input file of type '.cmg'\n");
printf("        outfile -- name of output file set\n");
printf("         dither -- amount of dithering:  default = 5 \n");
printf(
"     min points -- minimum number of pixels of a color to use\n");
printf(
"     max colors -- maximum number of colors in final image:  default = 256\n");
printf("\n");
printf("Example:\n\n");
printf("img8bit infile outfile 7 /m 10 /c 240\n");
	exit(0);
}





/***********************************************************************
**
**
**
**
************************************************************************/

int make_dither(fpin,sd)

FILE *fpin;
float sd;

{
	unsigned int i,j,k;
	long lk;
	int row,col,n,m,numwrote,val,numread;
	FILE *fpout;
	char cr=0xd,*d;

	if(sd<=0.0||sd>20.0)
	{
		printf("Bad standard deviation -- %f\n\n",sd);
		exit(0);
	}
	fread((char *)&row,sizeof(int),1,fpin);
	fread((char *)&col,sizeof(int),1,fpin);
	if(row<=0||col<=0)
	{
		printf("illegal row or column:  row = %d col = %d\n\n",row,col);
		exit(0);
	}
	printf("row = %d  col = %d\n",row,col);
	fpout=fopen("dither.cmg","wb");
	if(!fpout)
	{
		printf("Could not open 'dither.cmg'\n\n");
		exit(0);
	}
	fwrite((char *)&row,sizeof(int),1,fpout);
	fwrite((char *)&col,sizeof(int),1,fpout);
	printf("Loading Dither constants.\n");
	for(k=0;k<MAXCOL*4-1;k++)
		Dither[k]=rnum(0)*sd;
	printf("Creating Dithered image.\n");
	for(i=0;i<row;i++)
	{
		printf("%5d of %d%c",i+1,row,cr);
		for(k=0;k<3;k++)
			numread=fread((char *)Buffer[k],sizeof(char),col,fpin);
		m=rand()/10;
		if(m<0)
		{
			printf("rand() < 0\n\n");
			exit(0);
		}
		d=Dither+m;
		for(j=0;j<col;j++)
		{
			for(k=0;k<3;k++)
			{
				val=Buffer[k][j]+*d++;
				if(val>255)val=255;
				if(val<0)val=0;
				Buffer[k][j]=val;
			}
		}
		for(k=0;k<3;k++)
		{
			numwrote=fwrite((char *)Buffer[k],sizeof(char),col,fpout);
			if(numwrote<col)
			{
				printf("Disk full -- program ending.\n\n");
				unlink("dither.cmg");
				exit(0);
			}
		}
	}
	rewind(fpin);
}





/***********************************************************************
**
**
**
***********************************************************************/

void make_div(max)

int max;

{
	int i,j,k,j1,j2,n;
	double cmin[3],cmax[3],delc[3],maxc=max,d1,d2,dm;


	for(k=0;k<3;k++)
		for(j=0;j<256;j++)
			Rgb_Div[k][j]=-1;

	for(k=0;k<3;k++)
	{
		cmin[k]=log(AB[k][0])/AB[k][1];
		cmax[k]=log(AB[k][0]+maxc*AB[k][1])/AB[k][1];
		delc[k]=Del[k];
		delc[k]=(cmax[k]-cmin[k])/delc[k];
		for(i=0;i<Del[k];i++)
		{
			d1=i;
			d2=i+1;
			d1=cmin[k]+(d1*delc[k]);
			d2=cmin[k]+(d2*delc[k]);
			d1=(exp(AB[k][1]*d1)-AB[k][0])/AB[k][1];
			d2=(exp(AB[k][1]*d2)-AB[k][0])/AB[k][1];
			j1=d1+0.5;j2=d2+0.5;
			if(j1<0)j1=0;
			if(j2>max)j2=max;
			dm=(d1+d2)/2.0+0.5;
			RGB_Mid[k][i]=dm;
			for(j=j1;j<=j2;j++)
				Rgb_Div[k][j]=i;
		}
	}
	for(k=0;k<3;k++)
		for(j=0;j<256;j++)
		{
			if(Rgb_Div[k][j]>=0)
				RGB_Div[k][j]=Rgb_Div[k][j];
			else
			{
				printf("\n\nk = %d  j = %d has negative value\n\n",k,j);
				getch();
			}
		}

}

/***********************************************************************
**
**
**
***********************************************************************/

void measure_div(max)

int max;

{
	int i,j,k,j1,j2,n;
	double cmin[3],cmax[3],delc[3],maxc=max,d1,d2,dm;

	printf("Maximum discernable colors:\n");
	for(k=0;k<3;k++)
	{
		cmin[k]=log(AB[k][0])/AB[k][1];
		cmax[k]=log(AB[k][0]+maxc*AB[k][1])/AB[k][1];
		delc[k]=Del[k];
		delc[k]=(cmax[k]-cmin[k])/delc[k];
		printf("%lf -- %s\n",cmax[k]-cmin[k],Color_Name[k]);
	}
}


/********************************************************************
**
**
**
********************************************************************* */

void average_lut(lut,flut,table,buffer)

int lut[3][256];
float flut[4][256];
unsigned char huge table[RED][GRN][BLU],buffer[4][MAXCOL];
{
	int i,j,k,m,n;
	int sa=255/Del[0]+1,sb=255/Del[1]+1,sc=255/Del[2]+1;
	char cr=0xd;

	for(n=0;n<256;n++)
		for(m=0;m<4;m++)
		{
			flut[m][n]=0.0;
			lut[m][n]=255;
		}

	make_div(255);

	rewind(FpIn);
	fread((char *)&Row,sizeof(int),1,FpIn);
	fread((char *)&Col,sizeof(int),1,FpIn);
	for(i=0;i<Row;i++)
	{
		printf("%d of %d lines%c",i+1,Row,cr);
		for(k=0;k<3;k++)
			fread((char *)buffer[k],sizeof(char),Col,FpIn);
		for(j=0;j<Col;j++)
		{
			n=table[RGB_Div[0][buffer[0][j]]]
				 	 [RGB_Div[1][buffer[1][j]]]
				 	 [RGB_Div[2][buffer[2][j]]];
			flut[3][n]+=1.0;
			for(m=0;m<3;m++)
				flut[m][n]+=buffer[m][j];
		}
	}

	for(n=0;n<Num_Color;n++)
	{
		if(flut[3][n]>0.0)
			for(m=0;m<3;m++)
				lut[m][n]=flut[m][n]/flut[3][n]+.5;
	}
}



/***********************************************************************
**
**	iop = 1 to output 8-bit file and lut
**
***********************************************************************/

void 	make_img(table,lut,buffer,iop)

int lut[3][256],iop;
unsigned char huge table[RED][GRN][BLU],buffer[4][MAXCOL];

{
	int i,j,k,m;
	int a,b,c;
	int na,nb,nc,cval;
	int sa=255/Del[0]+1,sb=255/Del[1]+1,sc=255/Del[2]+1;
	FILE *fpimg,*fppal,*fplbl,*fpchn;
	char filename[40],palname[40],cr=0xd,string[100];

	rewind(FpIn);
	fread((char *)&Row,sizeof(int),1,FpIn);
	fread((char *)&Col,sizeof(int),1,FpIn);
	if(iop==1)
	{
		strcpy(Name,OutFileName);
		for(k=0;k<strlen(Name);k++)
			if(Name[k]=='.')
				Name[k]='\0';
		strcpy(palname,Name);
		strcat(palname,".pal");
		fppal=fopen(palname,"wt");
		strcpy(filename,Name);
		strcat(filename,".lbl");
		fplbl=fopen(filename,"wt");
		strcpy(filename,Name);
		strcat(filename,".dat");
		fpimg=fopen(filename,"wb");
		if(!fpimg)
		{
			printf("Can't open '%s' for output.\n",filename);
			exit(0);
		}
		fprintf(fplbl,"FILE_TYPE          = IMAGE\n");
		fprintf(fplbl,"IMAGE_LINES        = %4d\n",Row);
		fprintf(fplbl,"LINE_SAMPLES       = %4d\n",Col);
		fprintf(fplbl,"IMAGE_POINTER      = '%s'\n",filename);
		fprintf(fplbl,"PAL_POINTER        = '%s'\n",palname);
		fprintf(fplbl,"PALETTE_POINTER    = '%s'\n",palname);
/* include channel info if available */
		strcpy(Name,InFileName);
		for(k=0;k<strlen(Name);k++)
			if(Name[k]=='.')
				Name[k]='\0';
		sprintf(string,"%s.chn",Name);
		fpchn=fopen(string,"rt");
		if(fpchn)
		{
			while(fgets(string,99,fpchn))
			{
				fprintf(fplbl,"%s",string);
				printf("%s",string);
			}
			fclose(fpchn);
		}
		else
		{
			printf("No channel information ('%s') found.\n\n",string);
		}
		fclose(fplbl);
		for(i=0;i<256;i++)
			fprintf(fppal,"%3d %3d %3d %3d\n",i,lut[0][i],lut[1][i],lut[2][i]);
		fclose(fppal);
	}
	na=Del[0];nb=Del[1];nc=Del[2];
	make_div(255);
	for(i=0;i<Row;i++)
	{
		for(k=0;k<3;k++)
			fread((char *)buffer[k],sizeof(char),Col,FpIn);
		for(j=0;j<Col;j++)
		{
			a=RGB_Div[0][buffer[0][j]];
			b=RGB_Div[1][buffer[1][j]];
			c=RGB_Div[2][buffer[2][j]];
			m=table[a][b][c];
			buffer[0][j]=lut[0][m];
			buffer[1][j]=lut[1][m];
			buffer[2][j]=lut[2][m];
			buffer[4][j]=m;
		}
		printf("%d of %d lines%c",i+1,Row,cr);
		if(iop==1)
		{
			if(fwrite((char *)buffer[4],sizeof(char),Col,fpimg)<Col)
			{
				printf("Disk is full!\n\n");
				fclose(fpimg);
				fclose(fplbl);
				fclose(fppal);
				exit(0);
			}
		}
	}
	fclose(fpimg);
	fclose(fplbl);
	fclose(fppal);
}


/********************************************************************
**
**
**
********************************************************************* */

void fill_table(table,hist,hist2,min)

unsigned char huge table[RED][GRN][BLU];
int huge hist[RED][GRN][BLU],huge hist2[RED][GRN][BLU],min;

{
	int i,j,k,a,b,c,hit;
	int na,nb,nc,a1,a2,b1,b2,c1,c2;

	na=Del[0];nb=Del[1];nc=Del[2];
	for(a=0;a<na;a++)
	{
		for(b=0;b<nb;b++)
		{
			for(c=0;c<nc;c++)
			{
				if(hist[a][b][c]<min&&hist2[a][b][c]>0.0)
				{									/* find closest point */
					hit=0;
					if(a>0)a1=a-1;
					else a1=0;
					if(a<na-1)a2=a+1;
					else a2=na-1;
					if(b>0)b1=b-1;
					else b1=0;
					if(b<nb-1)b2=b+1;
					else b1=b2-1;
					if(c>0)c1=c-1;
					else c1=0;
					if(c<nc-1)c2=c+1;
					else c2=nc-1;

					do
					{
						if(a1>0)a1-=1;
						if(a2<na-1)a2+=1;
						if(b1>0)b1-=1;
						if(b2<na-1)b2+=1;
						if(c1>0)c1-=1;
						if(c2<nc-1)c2+=1;

						for(i=a1;i<=a2;i++)
							for(j=b1;j<=b2;j++)
								for(k=c1;k<=c2;k++)
									if(hist[i][j][k]>=min)
									{
										hit=1;
										table[a][b][c]=table[i][j][k];
									}
					}while(hit==0);
				}
			}
		}
	}
}




/********************************************************************
**
**
**
********************************************************************* */

int linear1(buffer,hist,lut,table,min)

int huge hist[RED][GRN][BLU],*lut,min;
unsigned char buffer[4][MAXCOL];
unsigned char huge *table;

{
	int i,j,k;
	int num;
	int *phist;
	float fdel[3],ddel[3],max=GRN;
	double mc3=Max_Color,power=0.33333333333;
	int enough=Max_Color;

	enough*=19;
	enough/=20;
	printf("enough colors = %d\n",enough);
	mc3=pow(mc3,power);
	Del[0]=ddel[0]=fdel[0]=RED;
	Del[1]=ddel[1]=fdel[1]=GRN;
	Del[2]=ddel[2]=fdel[2]=BLU;
	printf("%2d red x %2d grn x %2d blu -- ",Del[0],Del[1],Del[2]);
	num=var_hist(buffer,hist,min);
	printf("%d colors\n",num);

	for(k=0;k<3;k++)
	{	
		fdel[k]=fdel[k]/26.0*mc3;
		ddel[k]/=max;
		Del[k]=fdel[k];
/*		printf("Del = %d  fdel = %f  ddel= %f\n",Del[k],fdel[k],ddel[k]);*/
	}
	measure_div(255);
	do
	{
		printf("%2d red x %2d grn x %2d blu -- ",Del[0],Del[1],Del[2]);
		num=var_hist(buffer,hist,min);
		printf("%d colors\n",num);
		if(num<=enough/4)
			for(k=0;k<3;k++)
			{
				fdel[k]+=ddel[k]*6;
				Del[k]=fdel[k];
			}
		else if(num<=enough/3)
			for(k=0;k<3;k++)
			{
				fdel[k]+=ddel[k]*4;
				Del[k]=fdel[k];
			}
		else if(num<=enough/2)
			for(k=0;k<3;k++)
			{
				fdel[k]+=ddel[k]*3;
				Del[k]=fdel[k];
			}
		else if(num<=enough*3/4)
			for(k=0;k<3;k++)
			{
				fdel[k]+=ddel[k]*2;
				Del[k]=fdel[k];
			}
		else if(num<enough)
			for(k=0;k<3;k++)
			{
				fdel[k]+=ddel[k];
				Del[k]=fdel[k];
			}
	}while(num<enough&&Del[1]<GRN-1);
	k=2;
	if(num>Max_Color)
		do
		{
			Del[k]-=1;
			if(k>0)k-=1;
			else k=2;
			num=var_hist(buffer,hist,min);
			printf("Image has %d colors at  %d x %d x %d\n",
							num,Del[0],Del[1],Del[2]);
		}while(num>Max_Color);
	make_lut(hist,lut,table,min);

	return(num);
}


/********************************************************************
**
**
**
********************************************************************* */

void make_lut(hist,lut,table,min)

int huge hist[RED][GRN][BLU],lut[3][256],min;
unsigned char huge table[RED][GRN][BLU];

{
	int a,b,c,n=0,k;
	int na,nb,nc;
	int sa=255/Del[0]+1,sb=255/Del[1]+1,sc=255/Del[2]+1;

	na=Del[0];nb=Del[1];nc=Del[2];

	for(a=0;a<RED;a++)
		for(b=0;b<GRN;b++)
			for(c=0;c<BLU;c++)
				table[a][b][c]=0;

	for(a=0;a<na;a++)
		for(b=0;b<nb;b++)
			for(c=0;c<nc;c++)
				if(hist[a][b][c]>=min)
				{
					table[a][b][c]=n;
					lut[0][n]=RGB_Mid[0][a];
					lut[1][n]=RGB_Mid[1][b];
					lut[2][n]=RGB_Mid[2][c];
					for(k=0;k<3;k++)
					{
						if(lut[k][n]>255)lut[k][n]=255;
						if(lut[k][n]<0)  lut[k][n]=0;
					}
					n+=1;
					if(n>255)return;
				}
	return;
}


/********************************************************************
**
**
**
********************************************************************* */

int var_hist(buffer,hist,min)

int huge hist[RED][GRN][BLU];
unsigned char buffer[4][MAXCOL];

{
	int i,j,k,a,b,c;
	int sa=255/Del[0]+1,sb=255/Del[1]+1,sc=255/Del[2]+1;
	int num=0;
	int na,nb,nc;

	na=Del[0];nb=Del[1];nc=Del[2];

	for(a=0;a<na;a++)
		for(b=0;b<nb;b++)
			for(c=0;c<nc;c++)
				hist[a][b][c]=0;

	rewind(FpIn);
	fread((char *)&Row,sizeof(int),1,FpIn);
	fread((char *)&Col,sizeof(int),1,FpIn);

	make_div(255);
	for(i=0;i<Row;i++)
	{
		for(k=0;k<3;k++)
			fread((char *)buffer[k],sizeof(char),Col,FpIn);
		for(j=0;j<Col;j++)
			if(hist[RGB_Div[0][buffer[0][j]]]
				 [RGB_Div[1][buffer[1][j]]]
				 [RGB_Div[2][buffer[2][j]]]<32000)
				hist[RGB_Div[0][buffer[0][j]]]
					 [RGB_Div[1][buffer[1][j]]]
					 [RGB_Div[2][buffer[2][j]]]+=1;
	}

	for(a=0;a<na;a++)
		for(b=0;b<nb;b++)
			for(c=0;c<nc;c++)
				if(hist[a][b][c]>=min)
					num+=1;

	return(num);
}

/***********************************************************************
**
**
**
***********************************************************************/

float rnum(iop)

int iop;

{
	int i;
	int seed,numsum=TIMES;
	long lnum;
	float rn=0.0,num,fnum=numsum,scale;

	scale=numsum;
	scale=sqrt(scale*3.0);
	

	if(iop==1)
	{
		printf("give int number\n");
		scanf("%d",&seed);
		srand(seed);
	}
	for(i=0;i<numsum;i++)
	{
		lnum=rand();
		num=lnum;
		num-=16383.5;
		num/=16383.5;
		rn+=num;
	}
	return(rn/fnum*scale);

}