/* Hole Filler                   Version 2003.04.22              */
/* Grow hole edges one pixel at a time by interpolation.         */
/*                                                               */
/* Only fills void pixels that are next to non-void pixels.      */
/* Interpolates in eight directions using weighted average.      */
/* Iterates on the entire image until all voids are filled.      */
/* Thus grows edges by interpolation - best of both worlds.      */
/* Smoothing option is slow but effective (two-byte data only).  */

/* There is a lot of experimental code here.                     */
/* The program was intended to do two iterations in each pass    */
/* through the code and the screen output indicates that.        */
/* Input A --> Output B becomes Input B --> Output A in memory.  */
/* For experimental purposes this was bypassed so that the code  */
/* could be written just once.                                   */
/* Instead, for now, Output B is copied to Input A for the next  */
/* iteration.                                                    */


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

int number_input();
int yes_no_input();

main(void)
{
	FILE 			*fp1, *fp2;
	char			file1name[400], file2name[400], pause[400];
	unsigned long	filesize;
	long			lines, samps, infinit, temp;
	long			i, j, k, ii, jj, bytes;
	long			iii, jjj;
	long			count, test, prevlines, prevpix, test1, test2, thresh;
	long			bell, cycles, leftovers, smooth;
	float			howlong;
	float			nsteps, ssteps, wsteps, esteps;
	float			nwsteps, nesteps, swsteps, sesteps;
	float			north, south, east, west;
	float			nw, sw, ne, se, elev, range;
	float			nnesteps, enesteps, esesteps, ssesteps;
	float			sswsteps, wswsteps, wnwsteps, nnwsteps;
	float			nne, ene, ese, sse;
	float			ssw, wsw, wnw, nnw;	
	float			sum;
	
	time_t			t1, t2;
	
	unsigned char	**dns1 , **dns2 , **dns3 , *buffy ;
	unsigned short	**dns01, **dns02, **dns03, *buffy2;
	
	printf( "\n" ); 
	
	printf( ">>> Hole Filler                   Version 2003.04.22              <<<\n" ); 
	printf( ">>> Grow hole edges one pixel at a time by interpolation.         <<<\n" );
	printf( ">>>                                                               <<<\n" );
	printf( ">>> Only fills void pixels that are next to non-void pixels.      <<<\n" );
	printf( ">>> Interpolates in eight directions using weighted average.      <<<\n" );
	printf( ">>> Iterates on the entire image until all voids are filled.      <<<\n" );
	printf( ">>> Thus grows edges by interpolation - best of both worlds.      <<<\n" );
	printf( ">>> Smoothing option is slow but effective (two-byte data only).  <<<\n" );
	
	retry1:
	printf( "\nEnter INPUT IMAGE file NAME: " );
	gets(file1name);
	if ( (fp1 = fopen( file1name, "rb")) == NULL ) 
		{ printf( "*** No such file found! *** \a\n" ); 
		  goto retry1; }
	
	fseek( fp1, 0, SEEK_END );
	filesize = ftell( fp1 );
	fseek( fp1, 0, SEEK_SET );
	
	count = 0;
	
	printf( "\nData Type ( Default = ONE-BYTE ), type '2' for two-byte data: ");
	bytes = number_input();  if ( bytes > 1 ) bytes = 2; else bytes = 1;
	if ( bytes == 1 )   printf( "ONE-BYTE data specified.\n");
				 else { printf( "TWO-BYTE data specified.\n"); goto twobyte; }
	

	retry2:	
	do { printf( "\nNumber of SAMPS in image: ");
		  samps = number_input();
		  if (samps < 1) printf( "*** Must be positive integer. *** \a\n" );
	    } while ( samps < 1 );
	
	if ( filesize % samps ) 
	{ 
		printf( "*** Number of SAMPS is INCORRECT for file.  *** \n" );
		printf( "***  It has %ld pixels.  ", filesize );
		printf( "That is %6.5f lines.  ***\n\a", (float)filesize / samps );
		goto retry2;
	} 
	
	lines = filesize / samps;
	printf( "Number of LINES in image = %ld \n", lines );
	infinit = lines + samps;
	
	/* No smoothing yet for one-byte.
	
	retry4:
	printf( "\nSMOOTH the output with a boxfilter? ( Default = Yes, else type 'no'): ");
	smooth = yes_no_input();
	if ( smooth == 0 )  printf("Data patches will be SMOOTHED. \n");
	if ( smooth == 1 )  printf("Data patches will NOT be smoothed. \n");
	
	*/
	
	retry3:
	printf( "\nEnter OUTPUT IMAGE file NAME: " );
	gets(file2name);  
	if ( strcmp(file1name,file2name) == 0 )
	{
		printf("*** Cannot overwrite input file! ***\n\a");
		goto retry3;
	}
	if ( strlen( file2name ) == 0 ) { printf("\a"); goto retry3; }
	if ( (fp2 = fopen( file2name, "rb")) != NULL ) 
	{
		fclose( fp2 );
		printf( "*** File already EXISTS! *** \a\n\n" ); 
		printf( "Hit 'return' to REPLACE, type 'no' to rename: "); 
		if ( yes_no_input() > 0 ) goto retry3;
	}
	howlong = filesize / 44347.0 / 60 ;  /* 44347 = empirical speed factor */





	printf("\nPull down menu to 'Quit', else hit 'return':"); gets(pause);
	
	time(&t1);  printf( "\nStart                "); printf( ctime(&t1) );	

	fp2 = fopen( file2name, "wb" );
		
		
	dns1 = (unsigned char **) malloc( lines * sizeof( char * ) );
	if (dns1 == NULL) { printf("Not enough RAM.\a\n"); goto end; }
	
	for ( i = 0; i < lines; i ++ )				/* Get input RAM */
	{
		dns1[ i ] = ( unsigned char * ) malloc( samps );
		if (dns1[ i ] == NULL) { printf("Not enough RAM.\a\n"); goto end; }
	}	

	
	/*vvvvvvvvvvvv  Start Buffering of Input vvvvvvvvvvvv*/

	buffy = (unsigned char *) malloc( samps*100 );	
	if (buffy == NULL) { printf("Not enough RAM.\a\n"); goto end; }
		
	cycles = lines / 100;
	for ( k = 0; k < cycles; k++ )
	{
		fread( buffy, 1, 100*samps, fp1 );
		for ( i = 0; i < 100; i++ )
		{
			prevlines = ( k * 100 ) + i;
			prevpix   = i * samps;
			for ( j = 0; j < samps; j++ )
				dns1 [ prevlines ][ j ] = buffy [ prevpix + j ];
		}
	}

	leftovers = lines - ( cycles * 100 );
	fread( buffy, 1, samps * leftovers, fp1 );
	for ( i = 0; i < leftovers; i++ )
	{
		prevlines = ( cycles * 100 ) + i;
		prevpix   = i * samps;
		for ( j = 0; j < samps; j++ )
			dns1 [ prevlines ][ j ] = buffy [ prevpix + j ];
	}

	/*^^^^^^^^^^^^^  End Buffering of Input ^^^^^^^^^^^^^*/


	dns2 = (unsigned char **) malloc( lines * sizeof( char * ) );
	if (dns2 == NULL) { printf("Not enough RAM.\a\n"); goto end; }
	
	for ( i = 0; i < lines; i ++ )				/*  Get output RAM */
	{
		dns2[ i ] = ( unsigned char * ) malloc( samps );
		if (dns2[ i ] == NULL) { printf("Not enough RAM.\a\n"); goto end; }
	}	

	for ( i = 0; i < lines; i++ )
		for ( j = 0; j < samps; j++ )
			dns2[i][j] = 0;
	

/***************** Cycling of Run Pairs *********************/
	
	count += 2;
	printf( "\n" );
	printf( "Iterations %5ld   and %5ld.", count-1, count );
	fflush( stdout );
	goto afterfirst1;
	

	again1:	
	
	count += 2;
	
	printf( "\b\b\b\b\b\b\b\b\b\b\b\b\b\b" );
	printf( "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b" );
	printf( "Iterations %5ld   and %5ld.", count-1, count );
	fflush( stdout );
	
	afterfirst1:


/***************** FIRST RUN *********************/

/*
  fprint ( "First  Application:  " ); time(&t2); printf( ctime(&t2) );
*/
	
	for ( i = 0; i < lines; i++ )
	{
		for ( j = 0; j < samps; j++ )								
		{
			if ( dns1[i][j] > 0 ) { dns2[i][j] = dns1[i][j]; continue; } /* If not a void pixel */

			test = 0;		
			for ( ii = i-1; ii < i+2; ii++ ) /* Find it to be a void edge */
			{
				for ( jj = j-1; jj < j+2; jj++ )
				{
					if (( ii < 0 ) || ( ii >= lines )) continue;
					if (( jj < 0 ) || ( jj >= samps )) continue;
					if ( dns1[ii][jj] > 0 ) test = 1;			
				}
			}
			if ( test == 0 ) { dns2[i][j] = 0; continue; }  /* If not a void edge */
			
			
			/* For void edge pixels: */			
			
			nsteps = 0.0; ssteps = 0.0; esteps = 0.0; wsteps = 0.0;
			nwsteps = 0.0; nesteps = 0.0; swsteps = 0.0; sesteps = 0.0;
			north = 0.0; south = 0.0; east = 0.0; west = 0.0;
			nw = 0.0; ne = 0.0; sw = 0.0; se = 0.0;


			
			for ( ii = i - 1; ii >= 0; ii-- )    /* North */
			{
				if ( ii <= 0 ) { north = 0; nsteps = 0; break; }
				if ( dns1[ii][j] > 0 ) { north = dns1[ii][j]; nsteps = i - ii; break; }
			}
					
			for ( ii = i + 1; ii < lines; ii++ ) /* South */
			{
				if ( ii >= lines-1 ) { south = 0; ssteps = 0; break; }
				if ( dns1[ii][j] > 0 ) { south = dns1[ii][j]; ssteps = ii - i; break; }
			}
				
			for ( jj = j - 1; jj >= 0; jj-- )    /* West */
			{
				if ( jj <= 0 ) { west = 0; wsteps = 0; break; }
				if ( dns1[i][jj] > 0 ) { west = dns1[i][jj]; wsteps = j - jj; break; }
			}
				
			for ( jj = j + 1; jj < samps; jj++ ) /* East */
			{
				if ( jj >= samps-1 ) { east = 0; esteps = 0; break; }
				if ( dns1[i][jj] > 0 ) { east = dns1[i][jj]; esteps = jj - j; break; }
			}
				
						

			/* Southeast */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i + k; jj = j + k;
				if (( ii == lines ) || ( jj == samps )) { se = 0; sesteps = 0; break; }
				if ( dns1[ii][jj] > 0 ) { se = dns1[ii][jj]; sesteps = 1.4142 * k; break; }
			}
			
			
			/* Northeast */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i - k; jj = j + k;
				if (( ii == -1 ) || ( jj == samps )) { ne = 0; nesteps = 0; break; }
				if ( dns1[ii][jj] > 0 ) { ne = dns1[ii][jj]; nesteps = 1.4142 * k; break; }
			}
			
			/* Northwest */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i - k; jj = j - k;
				if (( ii == -1 ) || ( jj == -1 )) { nw = 0; nwsteps = 0; break; }
				if ( dns1[ii][jj] > 0 ) { nw = dns1[ii][jj]; nwsteps = 1.4142 * k; break; }
			}
			
			
			/* Southwest */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i + k; jj = j - k;
				if (( ii == lines ) || ( jj == -1 )) { sw = 0; swsteps = 0; break; }
				if ( dns1[ii][jj] > 0 ) { sw = dns1[ii][jj]; swsteps = 1.4142 * k; break; }
			}
			
		
			elev = 0; range = 0;
			if (  nsteps > 0.5 )  { elev += north /  nsteps; range +=  1.00 /  nsteps; }
			if (  ssteps > 0.5 )  { elev += south /  ssteps; range +=  1.00 /  ssteps; }
			if (  esteps > 0.5 )  { elev +=  east /  esteps; range +=  1.00 /  esteps; }
			if (  wsteps > 0.5 )  { elev +=  west /  wsteps; range +=  1.00 /  wsteps; }
			if ( nwsteps > 0.5 )  { elev +=    nw / nwsteps; range +=  1.00 / nwsteps; }
			if ( nesteps > 0.5 )  { elev +=    ne / nesteps; range +=  1.00 / nesteps; }
			if ( swsteps > 0.5 )  { elev +=    sw / swsteps; range +=  1.00 / swsteps; }
			if ( sesteps > 0.5 )  { elev +=    se / sesteps; range +=  1.00 / sesteps; }
			
										
			temp = ( elev / range ) + 0.5 ;
			if ( temp > 255 ) temp = 255;
			if ( temp <   0 ) temp =   0;

			dns2[i][j] = temp ;			

		}
	}


/***************** SECOND RUN *********************/

/*
fprintf ( stdout, "Second Application:  " ); time(&t2); printf( ctime(&t2) );
*/
	
	for ( i = 0; i < lines; i++ )
	{
		for ( j = 0; j < samps; j++ )								
		{
			if ( dns2[i][j] > 0 ) { dns1[i][j] = dns2[i][j]; continue; } /* If not a void pixel */

			test = 0;		
			for ( ii = i-1; ii < i+2; ii++ ) /* Find it to be a void edge */
			{
				for ( jj = j-1; jj < j+2; jj++ )
				{
					if (( ii < 0 ) || ( ii >= lines )) continue;
					if (( jj < 0 ) || ( jj >= samps )) continue;
					if ( dns2[ii][jj] > 0 ) test = 1;			
				}
			}
			if ( test == 0 ) { dns1[i][j] = 0; continue; }  /* If not a void edge */
			
			
			/* For void edge pixels: */			
			
			nsteps = 0.0; ssteps = 0.0; esteps = 0.0; wsteps = 0.0;
			nwsteps = 0.0; nesteps = 0.0; swsteps = 0.0; sesteps = 0.0;
			north = 0.0; south = 0.0; east = 0.0; west = 0.0;
			nw = 0.0; ne = 0.0; sw = 0.0; se = 0.0;


			
			for ( ii = i - 1; ii >= 0; ii-- )    /* North */
			{
				if ( ii <= 0 ) { north = 0; nsteps = 0; break; }
				if ( dns2[ii][j] > 0 ) { north = dns2[ii][j]; nsteps = i - ii; break; }
			}
					
			for ( ii = i + 1; ii < lines; ii++ ) /* South */
			{
				if ( ii >= lines-1 ) { south = 0; ssteps = 0; break; }
				if ( dns2[ii][j] > 0 ) { south = dns2[ii][j]; ssteps = ii - i; break; }
			}
				
			for ( jj = j - 1; jj >= 0; jj-- )    /* West */
			{
				if ( jj <= 0 ) { west = 0; wsteps = 0; break; }
				if ( dns2[i][jj] > 0 ) { west = dns2[i][jj]; wsteps = j - jj; break; }
			}
				
			for ( jj = j + 1; jj < samps; jj++ ) /* East */
			{
				if ( jj >= samps-1 ) { east = 0; esteps = 0; break; }
				if ( dns2[i][jj] > 0 ) { east = dns2[i][jj]; esteps = jj - j; break; }
			}
				
						

			/* Southeast */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i + k; jj = j + k;
				if (( ii == lines ) || ( jj == samps )) { se = 0; sesteps = 0; break; }
				if ( dns2[ii][jj] > 0 ) { se = dns2[ii][jj]; sesteps = 1.4142 * k; break; }
			}
			
			
			/* Northeast */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i - k; jj = j + k;
				if (( ii == -1 ) || ( jj == samps )) { ne = 0; nesteps = 0; break; }
				if ( dns2[ii][jj] > 0 ) { ne = dns2[ii][jj]; nesteps = 1.4142 * k; break; }
			}
			
			/* Northwest */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i - k; jj = j - k;
				if (( ii == -1 ) || ( jj == -1 )) { nw = 0; nwsteps = 0; break; }
				if ( dns2[ii][jj] > 0 ) { nw = dns2[ii][jj]; nwsteps = 1.4142 * k; break; }
			}
			
			
			/* Southwest */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i + k; jj = j - k;
				if (( ii == lines ) || ( jj == -1 )) { sw = 0; swsteps = 0; break; }
				if ( dns2[ii][jj] > 0 ) { sw = dns2[ii][jj]; swsteps = 1.4142 * k; break; }
			}
			
		
			elev = 0; range = 0;
			if (  nsteps > 0.5 )  { elev += north /  nsteps; range +=  1.00 /  nsteps; }
			if (  ssteps > 0.5 )  { elev += south /  ssteps; range +=  1.00 /  ssteps; }
			if (  esteps > 0.5 )  { elev +=  east /  esteps; range +=  1.00 /  esteps; }
			if (  wsteps > 0.5 )  { elev +=  west /  wsteps; range +=  1.00 /  wsteps; }
			if ( nwsteps > 0.5 )  { elev +=    nw / nwsteps; range +=  1.00 / nwsteps; }
			if ( nesteps > 0.5 )  { elev +=    ne / nesteps; range +=  1.00 / nesteps; }
			if ( swsteps > 0.5 )  { elev +=    sw / swsteps; range +=  1.00 / swsteps; }
			if ( sesteps > 0.5 )  { elev +=    se / sesteps; range +=  1.00 / sesteps; }
			
										
			temp = ( elev / range ) + 0.5 ;
			if ( temp > 255 ) temp = 255;
			if ( temp <   0 ) temp =   0;

			dns1[i][j] = temp ;			

		}
	}


/***************** CHECK FOR REMAINING ZEROES *********************/


	for ( i = 0; i < lines; i++ )
		for ( j = 0; j < samps; j++ )
			if ( dns1 [i][j] == 0 ) goto again1;
		

	
	/*vvvvvvvvvvvv  Start Buffering of Output vvvvvvvvvvvv*/

	cycles = lines / 100;
	for ( k = 0; k < cycles; k++ )
	{
		for ( i = 0; i < 100; i++ )
		{
			prevlines = ( k * 100 ) + i;
			prevpix   = i * samps;
			for ( j = 0; j < samps; j++ )
				buffy [ prevpix + j ] = dns1 [ prevlines ][ j ];
		}
		fwrite( buffy, 1, 100*samps, fp2 );
	}

	leftovers = lines - ( cycles * 100 );
	for ( i = 0; i < leftovers; i++ )
	{
		prevlines = ( cycles * 100 ) + i;
		prevpix   = i * samps;
		for ( j = 0; j < samps; j++ )
			buffy [ prevpix + j ] = dns1 [ prevlines ][ j ];
	}
	fwrite( buffy, 1, samps * leftovers, fp2 );

	
	/*^^^^^^^^^^^^^  End Buffering of Output ^^^^^^^^^^^^^*/

	
	fclose( fp1 ); fclose( fp2 );
			
	time(&t2);	printf ( "\n\nEnd                  "); printf( ctime(&t2) );
	printf( "\n                           Total time:        ");
	printf( "%4.2f minutes. \n", difftime(t2,t1)/60 );	
	goto end;
	





	
/*=============================  TWO-BYTE VERSION ================================*/
	
	twobyte:
	
	retry02:	
	do { printf( "\nNumber of SAMPS in image: ");
		  samps = number_input();
		  if (samps < 1) printf( "*** Must be positive integer. *** \a\n" );
	    } while ( samps < 1 );
	
	if ( filesize % samps ) 
	{ 
		printf( "*** Number of SAMPS is INCORRECT for file.  *** \n" );
		printf( "***  It has %ld pixels.  ", filesize );
		printf( "That is %6.5f lines.  ***\n\a", (float)filesize / samps );
		goto retry02;
	} 
	
	lines = filesize / ( samps * sizeof(short) );
	printf( "Number of LINES in image = %ld \n", lines );
	infinit = lines + samps;	


	
	retry04:
	printf( "\nSMOOTH the output with a boxfilter? ( Default = Yes, else type 'no'): ");
	smooth = yes_no_input();
	if ( smooth == 0 )  printf("Data patches will be SMOOTHED. \n");
	if ( smooth == 1 )  printf("Data patches will NOT be smoothed. \n");
	


		
	retry03:
	printf( "\nEnter OUTPUT IMAGE file NAME: " );
	gets(file2name);  
	if ( strcmp(file1name,file2name) == 0 )
	{
		printf("*** Cannot overwrite input file! ***\n\a");
		goto retry03;
	}
	if ( strlen( file2name ) == 0 ) { printf("\a"); goto retry03; }
	if ( (fp2 = fopen( file2name, "rb")) != NULL ) 
	{
		fclose( fp2 );
		printf( "*** File already EXISTS! *** \a\n\n" ); 
		printf( "Hit 'return' to REPLACE, type 'no' to rename: "); 
		if ( yes_no_input() > 0 ) goto retry03;
	}
	howlong = filesize / 44347.0 / 60 ;  /* 44347 = empirical speed factor */




	printf("\nPull down menu to 'Quit', else hit 'return':"); gets(pause);
	
	time(&t1);  printf( "\nStart                "); printf( ctime(&t1) );	

	fp2 = fopen( file2name, "wb" );
		
		
	dns01 = (unsigned short **) malloc( lines * sizeof( short * ) );
	if (dns01 == NULL) { printf("Not enough RAM.\a\n"); goto end; }
	
	for ( i = 0; i < lines; i ++ )				/* Get input file */
	{
		dns01[ i ] = ( unsigned short * ) malloc( samps * sizeof(short) );
		if (dns01[ i ] == NULL) { printf("Not enough RAM.\a\n"); goto end; }
	}	

	
	/*vvvvvvvvvvvv  Start Buffering of Input vvvvvvvvvvvv*/

	buffy2 = (unsigned short *) malloc( samps*100 * sizeof(short) );	
	if (buffy2 == NULL) { printf("Not enough RAM.\a\n"); goto end; }
		
	cycles = lines / 100;
	for ( k = 0; k < cycles; k++ )
	{
		fread( buffy2, sizeof(short), 100*samps, fp1 );
		for ( i = 0; i < 100; i++ )
		{
			prevlines = ( k * 100 ) + i;
			prevpix   = i * samps;
			for ( j = 0; j < samps; j++ )
				dns01 [ prevlines ][ j ] = buffy2 [ prevpix + j ];
		}
	}

	leftovers = lines - ( cycles * 100 );
	fread( buffy2, sizeof(short), samps * leftovers, fp1 );
	for ( i = 0; i < leftovers; i++ )
	{
		prevlines = ( cycles * 100 ) + i;
		prevpix   = i * samps;
		for ( j = 0; j < samps; j++ )
			dns01 [ prevlines ][ j ] = buffy2 [ prevpix + j ];
	}

	/*^^^^^^^^^^^^^  End Buffering of Input ^^^^^^^^^^^^^*/
	

	dns02 = (unsigned short **) malloc( lines * sizeof( short * ) );
	if (dns02 == NULL) { printf("Not enough RAM.\a\n"); goto end; }
	
	for ( i = 0; i < lines; i ++ )				/*  Get output RAM */
	{
		dns02[ i ] = ( unsigned short * ) malloc( samps * sizeof(short) );
		if (dns02[ i ] == NULL) { printf("Not enough RAM.\a\n"); goto end; }
	}	

	
/********** NEW: Eliminate loner pixels ***************/	

/*	
	for ( i = 0; i < lines; i++ )
	{
		for ( j = 0; j < samps; j++ )								
		{
			dns02[i][j] = dns01[i][j];
			test = 0;
			if ((( i == 0 ) || ( dns01[i-1][j] < 1 ) ) && ( ( i == lines-1 ) || ( dns01[i+1][j] < 1 ))) test += 1;
			if ((( j == 0 ) || ( dns01[i][j-1] < 1 ) ) && ( ( j == samps-1 ) || ( dns01[i][j+1] < 1 ))) test += 1;
			if ( test > 0 ) dns02[i][j] = 0;
		}
	}	
	
	for ( i = 0; i < lines; i++ )
	{
		for ( j = 0; j < samps; j++ )								
		{
			dns01[i][j] = dns02[i][j];
		}
	}	
*/	
	
	
/***************** Cycling of Run Pairs *********************/
	


	count += 2;
	printf( "\n" );
	printf( "Iterations %5ld   and %5ld.", count-1, count );
	fflush( stdout );
	goto afterfirst2;
	

	again2:	
	
	count += 2;
	
	printf( "\b\b\b\b\b\b\b\b\b\b\b\b\b\b" );
	printf( "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b" );
	printf( "Iterations %5ld   and %5ld.", count-1, count );
	fflush( stdout );
	
	afterfirst2:




/***************** FIRST RUN *********************/

/*
fprintf ( stdout, "First  Application:  " ); time(&t2); printf( ctime(&t2) );
*/
	
	for ( i = 0; i < lines; i++ )
	{
		for ( j = 0; j < samps; j++ )								
		{
			if ( dns01[i][j] > 0 ) { dns02[i][j] = dns01[i][j]; continue; } /* If not a void pixel */


			/* Find it to be a void edge */

/* Sides Only Test1 */
/*
			test1 = 0;
			if (( i >       0 ) && ( dns01[i-1][j] > 0 )) test1 += 1;
			if (( i < lines-1 ) && ( dns01[i+1][j] > 0 )) test1 += 1;
			if (( j >       0 ) && ( dns01[i][j-1] > 0 )) test1 += 1;
			if (( j < samps-1 ) && ( dns01[i][j+1] > 0 )) test1 += 1;

			if ( test1 < 1 ) { dns02[i][j] = 0; continue; } 
*/
/* Diagnals Only Test2 */

/*			test2 = 0;
			if (( i > 0 ) && ( j > 0 ) && ( dns01[i-1][j-1] > 0 )) test2 += 1;
			if (( i < lines-1 ) && ( j < samps-1 ) && ( dns01[i+1][j+1] > 0 )) test2 += 1;
			if (( i > 0 ) && ( j < samps-1 ) && ( dns01[i-1][j+1] > 0 )) test2 += 1;
			if (( i < lines-1 ) && ( j > 0 ) && ( dns01[i+1][j-1] > 0 )) test2 += 1;

			if ( test2 < 2 ) { dns02[i][j] = 0; continue; }

			if ((  test1 < 1 ) && ( test2 < 2 )) { dns02[i][j] = 0; continue; }
*/






/* Count All Neighbors Test */

			test = 0;		
			for ( ii = i-1; ii < i+2; ii++ ) 
			{
				for ( jj = j-1; jj < j+2; jj++ )
				{
					if (( ii < 0 ) || ( ii >= lines )) continue;
					if (( jj < 0 ) || ( jj >= samps )) continue;
					if ( dns01[ii][jj] > 0 ) test += 1;			
				}
			}
			if ( test < 3 ) { dns02[i][j] = 0; continue; } 







			
			
			/* For void edge pixels: */			
			
			nsteps = 0.0; ssteps = 0.0; esteps = 0.0; wsteps = 0.0;
			nwsteps = 0.0; nesteps = 0.0; swsteps = 0.0; sesteps = 0.0;
			north = 0.0; south = 0.0; east = 0.0; west = 0.0;
			nw = 0.0; ne = 0.0; sw = 0.0; se = 0.0;

			nnesteps = 0.0; enesteps = 0.0; esesteps = 0.0; ssesteps = 0.0;
			sswsteps = 0.0; wswsteps = 0.0; wnwsteps = 0.0; nnwsteps = 0.0;
			nne = 0.0; ene = 0.0; ese = 0.0; sse = 0.0;
			ssw = 0.0; wsw = 0.0; wnw = 0.0; nnw = 0.0;


/** NSEW **/
			
			for ( ii = i - 1; ii >= 0; ii-- )    /* North */
			{
				if ( ii <= 0 ) { north = 0; nsteps = 0; break; }
				if ( dns01[ii][j] > 0 ) { north = dns01[ii][j]; nsteps = i - ii; break; }
			}
					
			for ( ii = i + 1; ii < lines; ii++ ) /* South */
			{
				if ( ii >= lines-1 ) { south = 0; ssteps = 0; break; }
				if ( dns01[ii][j] > 0 ) { south = dns01[ii][j]; ssteps = ii - i; break; }
			}
				
			for ( jj = j - 1; jj >= 0; jj-- )    /* West */
			{
				if ( jj <= 0 ) { west = 0; wsteps = 0; break; }
				if ( dns01[i][jj] > 0 ) { west = dns01[i][jj]; wsteps = j - jj; break; }
			}
				
			for ( jj = j + 1; jj < samps; jj++ ) /* East */
			{
				if ( jj >= samps-1 ) { east = 0; esteps = 0; break; }
				if ( dns01[i][jj] > 0 ) { east = dns01[i][jj]; esteps = jj - j; break; }
			}

				
/** Diagonals  **/						

			/* Southeast */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i + k; jj = j + k;
				if (( ii == lines ) || ( jj == samps )) { se = 0; sesteps = 0; break; }
				if ( dns01[ii][jj] > 0 ) { se = dns01[ii][jj]; sesteps = 1.4142 * k; break; }
			}
			
			
			/* Northeast */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i - k; jj = j + k;
				if (( ii == -1 ) || ( jj == samps )) { ne = 0; nesteps = 0; break; }
				if ( dns01[ii][jj] > 0 ) { ne = dns01[ii][jj]; nesteps = 1.4142 * k; break; }
			}
			
			/* Northwest */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i - k; jj = j - k;
				if (( ii == -1 ) || ( jj == -1 )) { nw = 0; nwsteps = 0; break; }
				if ( dns01[ii][jj] > 0 ) { nw = dns01[ii][jj]; nwsteps = 1.4142 * k; break; }
			}
			
			
			/* Southwest */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i + k; jj = j - k;
				if (( ii == lines ) || ( jj == -1 )) { sw = 0; swsteps = 0; break; }
				if ( dns01[ii][jj] > 0 ) { sw = dns01[ii][jj]; swsteps = 1.4142 * k; break; }
			}


/** Other 8 **/

			/* NNE */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i - (k*2); jj = j + k;
				if (( ii <= -1 ) || ( jj >= samps )) { nne = 0; nnesteps = 0; break; }
				if ( dns01[ii][jj] > 0 ) { nne = dns01[ii][jj]; nnesteps = 2.2361 * k; break; }
			}
			
			
			/* ENE */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i - k; jj = j + (k*2);
				if (( ii <= -1 ) || ( jj >= samps )) { ene = 0; enesteps = 0; break; }
				if ( dns01[ii][jj] > 0 ) { ene = dns01[ii][jj]; enesteps = 2.2361 * k; break; }
			}
			
			/* ESE */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i + k; jj = j + (k*2);
				if (( ii >= lines ) || ( jj >= samps )) { ese = 0; esesteps = 0; break; }
				if ( dns01[ii][jj] > 0 ) { ese = dns01[ii][jj]; esesteps = 2.2361 * k; break; }
			}
			
			
			/* SSE */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i + (k*2); jj = j + k;
				if (( ii >= lines ) || ( jj >= samps )) { sse = 0; ssesteps = 0; break; }
				if ( dns01[ii][jj] > 0 ) { sse = dns01[ii][jj]; ssesteps = 2.2361 * k; break; }
			}



			/* SSW */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i + (k*2); jj = j - k;
				if (( ii >= lines ) || ( jj <= -1 )) { ssw = 0; sswsteps = 0; break; }
				if ( dns01[ii][jj] > 0 ) { ssw = dns01[ii][jj]; sswsteps = 2.2361 * k; break; }
			}
			
			
			/* WSW */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i + k; jj = j - (k*2);
				if (( ii >= lines ) || ( jj <= -1 )) { wsw = 0; wswsteps = 0; break; }
				if ( dns01[ii][jj] > 0 ) { wsw = dns01[ii][jj]; wswsteps = 2.2361 * k; break; }
			}
			
			/* WNW */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i - k; jj = j - (k*2);
				if (( ii <= -1 ) || ( jj <= -1 )) { wnw = 0; wnwsteps = 0; break; }
				if ( dns01[ii][jj] > 0 ) { wnw = dns01[ii][jj]; wnwsteps = 2.2361 * k; break; }
			}
			
			
			/* NNW */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i - (k*2); jj = j - k;
				if (( ii <= -1 ) || ( jj <= -1 )) { nnw = 0; nnwsteps = 0; break; }
				if ( dns01[ii][jj] > 0 ) { nnw = dns01[ii][jj]; nnwsteps = 2.2361 * k; break; }
			}
			


/* Avoid fills that use too many distant points */
/*	
thresh = 0;
			if (  nsteps > 5 ) thresh += 1;
			if (  ssteps > 5 ) thresh += 1;
			if (  esteps > 5 ) thresh += 1;
			if (  wsteps > 5 ) thresh += 1;

			if ( nwsteps > 5 ) thresh += 1;
			if ( nesteps > 5 ) thresh += 1;
			if ( swsteps > 5 ) thresh += 1;
			if ( sesteps > 5 ) thresh += 1;

			if ( nnesteps > 5 ) thresh += 1;
			if ( enesteps > 5 ) thresh += 1;
			if ( esesteps > 5 ) thresh += 1;
			if ( ssesteps > 5 ) thresh += 1;
			if ( sswsteps > 5 ) thresh += 1;
			if ( wswsteps > 5 ) thresh += 1;
			if ( wnwsteps > 5 ) thresh += 1;
			if ( nnwsteps > 5 ) thresh += 1;		
if ( thresh > 14 ) { dns02[i][j] = 0; continue; }			
*/			


			
/* TRY SQUARES -- BAD! */
/*			
			nsteps = pow(nsteps, 0.5);
			ssteps = pow(ssteps, 0.5); 
			esteps = pow(esteps, 0.5);
			wsteps = pow(wsteps, 0.5);
			nwsteps = pow(nwsteps, 0.5);
			nesteps = pow(nesteps, 0.5);
			swsteps = pow(swsteps, 0.5);
			sesteps = pow(sesteps, 0.5);
			nnesteps = pow(nnesteps, 0.5);
			enesteps = pow(enesteps, 0.5);
			esesteps = pow(esesteps, 0.5);
			ssesteps = pow(ssesteps, 0.5);
			sswsteps = pow(sswsteps, 0.5);
			wswsteps = pow(wswsteps, 0.5);
			wnwsteps = pow(wnwsteps, 0.5);
			nnwsteps = pow(nnwsteps, 0.5);
*/			

			
		
			elev = 0; range = 0;
/*NSEW*/			
			if (  nsteps > 0.5 )  { elev += north /  nsteps; range +=  1.00 /  nsteps; }
			if (  ssteps > 0.5 )  { elev += south /  ssteps; range +=  1.00 /  ssteps; }
			if (  esteps > 0.5 )  { elev +=  east /  esteps; range +=  1.00 /  esteps; }
			if (  wsteps > 0.5 )  { elev +=  west /  wsteps; range +=  1.00 /  wsteps; }
/*Diag*/
			if ( nwsteps > 0.5 )  { elev +=    nw / nwsteps; range +=  1.00 / nwsteps; }
			if ( nesteps > 0.5 )  { elev +=    ne / nesteps; range +=  1.00 / nesteps; }
			if ( swsteps > 0.5 )  { elev +=    sw / swsteps; range +=  1.00 / swsteps; }
			if ( sesteps > 0.5 )  { elev +=    se / sesteps; range +=  1.00 / sesteps; }
			
/*Other 8*/
			if ( nnesteps > 0.5 )  { elev += nne / nnesteps; range +=  1.00 / nnesteps; }
			if ( enesteps > 0.5 )  { elev += ene / enesteps; range +=  1.00 / enesteps; }
			if ( esesteps > 0.5 )  { elev += ese / esesteps; range +=  1.00 / esesteps; }
			if ( ssesteps > 0.5 )  { elev += sse / ssesteps; range +=  1.00 / ssesteps; }
			if ( sswsteps > 0.5 )  { elev += ssw / sswsteps; range +=  1.00 / sswsteps; }
			if ( wswsteps > 0.5 )  { elev += wsw / wswsteps; range +=  1.00 / wswsteps; }
			if ( wnwsteps > 0.5 )  { elev += wnw / wnwsteps; range +=  1.00 / wnwsteps; }
			if ( nnwsteps > 0.5 )  { elev += nnw / nnwsteps; range +=  1.00 / nnwsteps; }
			
										
			temp = ( elev / range ) + 0.5 ;
			if ( temp > 10000 ) temp = 10000;
			if ( temp <     0 ) temp =     0;

			dns02[i][j] = temp ;			

		}
	}




goto there;





/***************** SECOND RUN *********************/

fprintf ( stdout, "Second Application:  " ); time(&t2); printf( ctime(&t2) );

	
	for ( i = 0; i < lines; i++ )
	{
		for ( j = 0; j < samps; j++ )								
		{
			if ( dns02[i][j] > 0 ) { dns01[i][j] = dns02[i][j]; continue; } /* If not a void pixel */

			test = 0;		
			for ( ii = i-1; ii < i+2; ii++ ) /* Find it to be a void edge */
			{
				for ( jj = j-1; jj < j+2; jj++ )
				{
					if (( ii < 0 ) || ( ii >= lines )) continue;
					if (( jj < 0 ) || ( jj >= samps )) continue;
					if ( dns02[ii][jj] > 0 ) test = 1;			
				}
			}
			if ( test == 0 ) { dns01[i][j] = 0; continue; }  /* If not a void edge */
			
			
			/* For void edge pixels: */			
			
			nsteps = 0.0; ssteps = 0.0; esteps = 0.0; wsteps = 0.0;
			nwsteps = 0.0; nesteps = 0.0; swsteps = 0.0; sesteps = 0.0;
			north = 0.0; south = 0.0; east = 0.0; west = 0.0;
			nw = 0.0; ne = 0.0; sw = 0.0; se = 0.0;


			
			for ( ii = i - 1; ii >= 0; ii-- )    /* North */
			{
				if ( ii <= 0 ) { north = 0; nsteps = 0; break; }
				if ( dns02[ii][j] > 0 ) { north = dns02[ii][j]; nsteps = i - ii; break; }
			}
					
			for ( ii = i + 1; ii < lines; ii++ ) /* South */
			{
				if ( ii >= lines-1 ) { south = 0; ssteps = 0; break; }
				if ( dns02[ii][j] > 0 ) { south = dns02[ii][j]; ssteps = ii - i; break; }
			}
				
			for ( jj = j - 1; jj >= 0; jj-- )    /* West */
			{
				if ( jj <= 0 ) { west = 0; wsteps = 0; break; }
				if ( dns02[i][jj] > 0 ) { west = dns02[i][jj]; wsteps = j - jj; break; }
			}
				
			for ( jj = j + 1; jj < samps; jj++ ) /* East */
			{
				if ( jj >= samps-1 ) { east = 0; esteps = 0; break; }
				if ( dns02[i][jj] > 0 ) { east = dns02[i][jj]; esteps = jj - j; break; }
			}
				
						

			/* Southeast */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i + k; jj = j + k;
				if (( ii == lines ) || ( jj == samps )) { se = 0; sesteps = 0; break; }
				if ( dns02[ii][jj] > 0 ) { se = dns02[ii][jj]; sesteps = 1.4142 * k; break; }
			}
			
			
			/* Northeast */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i - k; jj = j + k;
				if (( ii == -1 ) || ( jj == samps )) { ne = 0; nesteps = 0; break; }
				if ( dns02[ii][jj] > 0 ) { ne = dns02[ii][jj]; nesteps = 1.4142 * k; break; }
			}
			
			/* Northwest */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i - k; jj = j - k;
				if (( ii == -1 ) || ( jj == -1 )) { nw = 0; nwsteps = 0; break; }
				if ( dns02[ii][jj] > 0 ) { nw = dns02[ii][jj]; nwsteps = 1.4142 * k; break; }
			}
			
			
			/* Southwest */
			for ( k = 1; k < infinit; k++ )
			{
				ii = i + k; jj = j - k;
				if (( ii == lines ) || ( jj == -1 )) { sw = 0; swsteps = 0; break; }
				if ( dns02[ii][jj] > 0 ) { sw = dns02[ii][jj]; swsteps = 1.4142 * k; break; }
			}
			
		
			elev = 0; range = 0;
			if (  nsteps > 0.5 )  { elev += north /  nsteps; range +=  1.00 /  nsteps; }
			if (  ssteps > 0.5 )  { elev += south /  ssteps; range +=  1.00 /  ssteps; }
			if (  esteps > 0.5 )  { elev +=  east /  esteps; range +=  1.00 /  esteps; }
			if (  wsteps > 0.5 )  { elev +=  west /  wsteps; range +=  1.00 /  wsteps; }
			if ( nwsteps > 0.5 )  { elev +=    nw / nwsteps; range +=  1.00 / nwsteps; }
			if ( nesteps > 0.5 )  { elev +=    ne / nesteps; range +=  1.00 / nesteps; }
			if ( swsteps > 0.5 )  { elev +=    sw / swsteps; range +=  1.00 / swsteps; }
			if ( sesteps > 0.5 )  { elev +=    se / sesteps; range +=  1.00 / sesteps; }
			
										
			temp = ( elev / range ) + 0.5 ;
			if ( temp > 10000 ) temp = 10000;
			if ( temp <     0 ) temp =     0;

			dns01[i][j] = temp ;			

		}
	}


/***************** CHECK FOR REMAINING ZEROES *********************/


there:

for ( i = 0; i < lines; i++ )
for ( j = 0; j < samps; j++ )
dns01 [i][j] = dns02 [i][j];


if ( count > 3000 ) goto there2; /* Optional within-program early termination */


	for ( i = 0; i < lines; i++ )
		for ( j = 0; j < samps; j++ )
		if ( dns01 [i][j] == 0 ) goto again2;
		

if ( smooth == 1 ) goto there2;		
		
		
/************************ SMOOTH THE RESULT ***********************/		
		
		
		
	dns03 = (unsigned short **) malloc( lines * sizeof( short * ) );
	if (dns03 == NULL) { printf("Not enough RAM.\a\n"); goto end; }
	
	for ( i = 0; i < lines; i ++ )				/*  Get output RAM */
	{
		dns03[ i ] = ( unsigned short * ) malloc( samps * sizeof(short) );
		if (dns03[ i ] == NULL) { printf("Not enough RAM.\a\n"); goto end; }
	}	
	
		
	fseek( fp1, 0, SEEK_SET ); /* Reset pointer to Input Data */
	
	/* Load original data into RAM Array 3 (Need data / Non-Data Info) */
	
	/*vvvvvvvvvvvv  Start Buffering of Input vvvvvvvvvvvv*/

	buffy2 = (unsigned short *) malloc( samps*100 * sizeof(short) );	
	if (buffy2 == NULL) { printf("Not enough RAM.\a\n"); goto end; }
		
	cycles = lines / 100;
	for ( k = 0; k < cycles; k++ )
	{
		fread( buffy2, sizeof(short), 100*samps, fp1 );
		for ( i = 0; i < 100; i++ )
		{
			prevlines = ( k * 100 ) + i;
			prevpix   = i * samps;
			for ( j = 0; j < samps; j++ )
				dns03 [ prevlines ][ j ] = buffy2 [ prevpix + j ];
		}
	}

	leftovers = lines - ( cycles * 100 );
	fread( buffy2, sizeof(short), samps * leftovers, fp1 );
	for ( i = 0; i < leftovers; i++ )
	{
		prevlines = ( cycles * 100 ) + i;
		prevpix   = i * samps;
		for ( j = 0; j < samps; j++ )
			dns03 [ prevlines ][ j ] = buffy2 [ prevpix + j ];
	}

	/*^^^^^^^^^^^^^  End Buffering of Input ^^^^^^^^^^^^^*/
		
	fprintf(stdout, "\n\nSmoothing Line Counter ( %5ld total ):       ", lines );


	for ( i = 0; i < lines; i++ )
	{
		for ( j = 0; j < samps; j++ )
		{
			if ( dns03 [i][j] > 0 ) { dns02 [i][j] = dns01 [i][j]; continue; }



			for ( k = 1; k < infinit; k++ )  /* Find the smallest box size with data */
			{
				for ( ii = i-k; ii <= i+k; ii++ )
				{
					if ( ii <      0 ) continue;
					if ( ii >= lines ) break;
					for ( jj = j-k; jj <= j+k; jj++ )
					{
						if ( jj <      0 ) continue;
						if ( jj >= samps ) continue;
						if ( dns03[ii][jj] > 0 ) goto k_nowset;
					}
				}
			}
			
			k_nowset:
			k = k / 4; if ( k < 1 ) k = 1; /* Errrrr.  Make it fourth size */
			sum = 0;
			for ( ii = i-k; ii <= i+k; ii++ )
			{
				if ( ii <      0 ) iii = -1 - ii; else iii = ii;
				if ( ii >= lines ) iii = lines - (ii - lines + 1);
				for ( jj = j-k; jj <= j+k; jj++ )
				{
					if ( jj <      0 ) jjj = -1 - jj; else jjj = jj;
					if ( jj >= samps ) jjj = samps - (jj - samps + 1);
					sum += dns01[iii][jjj];
				}
			}
			k = (k*2)+1;
			dns02[i][j] = ( sum / k / k ) + 0.5;


			
		}
		if ( !(i%100) )
		fprintf(stdout, "\b\b\b\b\b\b%5ld ", i );   /* Count lines on output screen */
		fflush(stdout);	
	}
	fprintf(stdout, "\b\b\b\b\b\b%5ld ", lines );
	fflush(stdout);	
		
		
		
		
		
		
		
there2:		
		
	
	/*vvvvvvvvvvvv  Start Buffering of Output vvvvvvvvvvvv*/

	cycles = lines / 100;
	for ( k = 0; k < cycles; k++ )
	{
		for ( i = 0; i < 100; i++ )
		{
			prevlines = ( k * 100 ) + i;
			prevpix   = i * samps;
			for ( j = 0; j < samps; j++ )
				buffy2 [ prevpix + j ] = dns02 [ prevlines ][ j ];
		}
		fwrite( buffy2, sizeof(short), 100*samps, fp2 );
	}

	leftovers = lines - ( cycles * 100 );
	for ( i = 0; i < leftovers; i++ )
	{
		prevlines = ( cycles * 100 ) + i;
		prevpix   = i * samps;
		for ( j = 0; j < samps; j++ )
			buffy2 [ prevpix + j ] = dns02 [ prevlines ][ j ];
	}
	fwrite( buffy2, sizeof(short), samps * leftovers, fp2 );

	
	/*^^^^^^^^^^^^^  End Buffering of Output ^^^^^^^^^^^^^*/

	
	fclose( fp1 ); fclose( fp2 );
			
	time(&t2);	printf ( "\n\nEnd                  "); printf( ctime(&t2) );
	printf( "\n                           Total time:        ");
	printf( "%4.2f minutes. \n", difftime(t2,t1)/60 );		
	

	end:
	printf( "Done.\a\n" );
	
	return(0);
	
}

/*==================================================================================*/

int number_input()

{
	char buffer[16];
	char	*p;
	
	p=buffer;
	while ( (*p = getc(stdin) ) != '\n') p++;
	*p = '\0';
	return( atoi (buffer) );
}



/*==================================================================================*/

int yes_no_input()

{
	char buffer[16];
	char	*p;
	
	p=buffer;
	while ( (*p = getc(stdin) ) != '\n') p++;
	*p = '\0';
	if ( strlen(buffer) == 0 ) return ( 0 );
	return( 1 );
}




/*  EDITS:                                                              */
/*                                                                      */
/*  03.04.02	Modified from Zero Fill Cycle.                          */

