Index: /issm/trunk/src/c/modules/Gradjx/Gradjx.cpp =================================================================== --- /issm/trunk/src/c/modules/Gradjx/Gradjx.cpp (revision 3901) +++ /issm/trunk/src/c/modules/Gradjx/Gradjx.cpp (revision 3901) @@ -0,0 +1,35 @@ +/*!\file Gradjx + * \brief: compute inverse method gradient + */ + +#include "./Gradjx.h" + +#include "../shared/shared.h" +#include "../include/include.h" +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +void Gradjx( Vec* pgrad_g, int numberofnodes, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters, + int analysis_type,int sub_analysis_type,int control_type){ + + /*output: */ + Vec grad_g=NULL; + + /*First, get elements and loads configured: */ + elements->Configure(elements,loads, nodes,vertices, materials,parameters); + nodes->Configure(elements,loads, nodes,vertices, materials,parameters); + parameters->Configure(elements,loads, nodes,vertices, materials,parameters); + + /*Allocate grad_g: */ + grad_g=NewVec(numberofnodes); + + /*Compute gradients: */ + elements->Gradj(grad_g,analysis_type,sub_analysis_type,control_type); + + /*Assemble vector: */ + VecAssemblyBegin(grad_g); + VecAssemblyEnd(grad_g); + + /*Assign output pointers: */ + *pgrad_g=grad_g; +} Index: /issm/trunk/src/c/modules/Gradjx/Gradjx.h =================================================================== --- /issm/trunk/src/c/modules/Gradjx/Gradjx.h (revision 3901) +++ /issm/trunk/src/c/modules/Gradjx/Gradjx.h (revision 3901) @@ -0,0 +1,16 @@ +/*!\file: Gradjx.h + * \brief header file for inverse methods gradient computation + */ + +#ifndef _GRADJX_H +#define _GRADJX_H + +#include "../DataSet/DataSet.h" +#include "../objects/objects.h" + +/* local prototypes: */ +void Gradjx( Vec* pgrad_g, int numberofnodes, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters, + int analysis_type,int sub_analysis_type,int control_type); + +#endif /* _GRADJX_H */ + Index: /issm/trunk/src/c/modules/HoleFillerx/HoleFillerCrippen.c =================================================================== --- /issm/trunk/src/c/modules/HoleFillerx/HoleFillerCrippen.c (revision 3901) +++ /issm/trunk/src/c/modules/HoleFillerx/HoleFillerCrippen.c (revision 3901) @@ -0,0 +1,1296 @@ +/* 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 +#include +#include + + +#include +#include + +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. */ + Index: /issm/trunk/src/c/modules/HoleFillerx/HoleFillerx.cpp =================================================================== --- /issm/trunk/src/c/modules/HoleFillerx/HoleFillerx.cpp (revision 3901) +++ /issm/trunk/src/c/modules/HoleFillerx/HoleFillerx.cpp (revision 3901) @@ -0,0 +1,372 @@ +/*!\file: HoleFillerx.cpp + * \brief plug missing values of a double* array, using interpolation techniques. This + * is an adaptation of Dr. Crippen's routine (Jet Propulsion Lab, see HoleFillerCrippen.c + * for more details. + */ + +#include +#include +#include "./HoleFillerx.h" +#include "../shared/shared.h" + +int HoleFillerx(double** pimageout,double* image, int lines,int samps,int smooth){ + + FILE *fp1; + unsigned long filesize; + long infinit; + double temp; + long i, j, k, l, ii, jj, bytes; + long iii, jjj; + long count, test, prevlines, prevpix, test1, test2, thresh,counter; + long bell, cycles, leftovers; + 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; + + /*output: */ + double* imageout=NULL; + + double* image2=NULL; + double* image3=NULL; + double* image4=NULL; + + /*Added to original routine: */ + int imageoutsize; + + /*^^^^^^^^^^^^^ Remove pixels close to the holes ^^^^^^^^^^^^^*/ + image2 = (double*) xmalloc( lines*samps*sizeof(double)); + memcpy(image2,image,lines*samps*sizeof(double)); + + for ( i = 0; i < lines; i++ ){ + for ( j = 0; j < samps; j++ ){ + if ( *(image+i*samps+j) == 0 ){ + if ( (j > 3) && (j < samps-4) && (i > 3) && (i < lines-4)){ + for ( k = 0; k < 4; k++ ){ + for ( l = 0; l < 4; l++ ){ + *(image2+samps*(i+k)+j+l)=0; + *(image2+samps*(i-k)+j+l)=0; + *(image2+samps*(i+k)+j-l)=0; + *(image2+samps*(i-k)+j-l)=0; + + } + } + } + } + } + } + + image3 = (double*) xmalloc( lines*samps*sizeof(double)); + memcpy(image3,image2,lines*samps*sizeof(double)); + + count += 2; + #ifdef _DEBUG2_ + counter=0; + for ( i = 0; i < lines; i++ ) { + for ( j = 0; j < samps; j++ ){ + if ( *(image2+i*samps+j) != 0 ) { + counter++; + } + } + } + printf( "\n" ); + printf("Number of zeroes remaining: %10ld",lines*samps-counter); + fflush( stdout ); + #endif + + goto afterfirst2; + + again2: + + count += 2; + + #ifdef _DEBUG2_ + counter=0; + for ( i = 0; i < lines; i++ ) { + for ( j = 0; j < samps; j++ ){ + if ( *(image2+i*samps+j) != 0 ) { + counter++; + } + } + } + // n u m b e r o f z e r o e s r e m a i n i n g : 1 2 3 4 5 6 7 8 9 10 + printf( "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b" ); + printf("Number of zeroes remaining: %10ld",lines*samps-counter); + fflush( stdout ); + #endif + + 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 ( *(image2+i*samps+j) != 0 ){ + *(image3+samps*i+j) = *(image2+i*samps+j); + 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 ( *(image2+ii*samps+jj) != 0 ) test += 1; + } + } + if ( test < 3 ){ + *(image3+samps*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 ( *(image2+ii*samps+j) != 0 ){ north = *(image2+ii*samps+j); nsteps = i - ii; break;} + } + + for ( ii = i + 1; ii < lines; ii++ ){ /* South */ + if ( ii >= lines-1 ){ south = 0; ssteps = 0; break;} + if ( *(image2+ii*samps+j) != 0 ){ south = *(image2+ii*samps+j); ssteps = ii - i; break;} + } + + for ( jj = j - 1; jj >= 0; jj-- ){ /* West */ + if ( jj <= 0 ) { west = 0; wsteps = 0; break;} + if ( *(image2+i*samps+jj) != 0 ){ west = *(image2+i*samps+jj); wsteps = j - jj; break;} + } + + for ( jj = j + 1; jj < samps; jj++ ){ /* East */ + if ( jj >= samps-1 ){ east = 0; esteps = 0; break;} + if ( *(image2+i*samps+jj) != 0 ){ east = *(image2+i*samps+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 ( *(image2+ii*samps+jj) != 0 ){ se = *(image2+ii*samps+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 ( *(image2+ii*samps+jj) != 0 ){ ne = *(image2+ii*samps+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 ( *(image2+ii*samps+jj) != 0 ){ nw = *(image2+ii*samps+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 ( *(image2+ii*samps+jj) != 0 ){ sw = *(image2+ii*samps+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 ( *(image2+ii*samps+jj) != 0 ){ nne = *(image2+ii*samps+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 ( *(image2+ii*samps+jj) != 0 ){ ene = *(image2+ii*samps+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 ( *(image2+ii*samps+jj) != 0 ){ ese = *(image2+ii*samps+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 ( *(image2+ii*samps+jj) != 0 ){ sse = *(image2+ii*samps+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 ( *(image2+ii*samps+jj) != 0 ){ ssw = *(image2+ii*samps+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 ( *(image2+ii*samps+jj) != 0 ){ wsw = *(image2+ii*samps+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 ( *(image2+ii*samps+jj) != 0 ){ wnw = *(image2+ii*samps+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 ( *(image2+ii*samps+jj) != 0 ){ nnw = *(image2+ii*samps+jj); nnwsteps = 2.2361 * k; break;} + } + + 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 ; + temp = ( elev / range ); + //+ 0.5 ; + //if ( temp > 10000 ) temp = 10000; + //if ( temp < 0 ) temp = 0; + + #ifdef _DEBUG2_ + //printf("%g %g %g \n",temp,elev,range); + #endif + + *(image3+samps*i+j) = temp ; + } + } + + + + for ( i = 0; i < lines; i++ ){ + for ( j = 0; j < samps; j++ ){ + *(image2+samps*i+j)= *(image3+samps*i+j); + } + } + + + for ( i = 0; i < lines; i++ ){ + for ( j = 0; j < samps; j++ ){ + if ( *(image2+samps*i+j) == 0 ) goto again2; + } + } + + if ( smooth == 0 ) goto there2; + + +/************************ SMOOTH THE RESULT ***********************/ + + image4 = (double*) xmalloc( lines*samps*sizeof(double)); + memcpy(image4,image3,lines*samps*sizeof(double)); + + + for ( i = 0; i < lines; i++ ) { + for ( j = 0; j < samps; j++ ) { + if ( *(image4+i*samps+j) != 0 ) { *(image3+i*samps+j) = *(image2+i*samps+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 ( *(image4+ii*samps+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 += *(image2+samps*iii+jjj) ; + } + } + k = (k*2)+1; + *(image3+samps*i+j) = ( sum / k / k ) ;//+ 0.5; + } + + } + + +there2: + + + /*Allocate output image: */ + imageout=(double*)xmalloc(samps*lines*sizeof(double)); + memcpy(imageout,image3,lines*samps*sizeof(double)); + + time(&t2); + + #ifdef _DEBUG2_ + printf( "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b" ); + printf("Number of zeroes remaining: 0\n\n"); + printf ( "\n"); + #endif + + end: + + /*Assign output pointers: */ + *pimageout=imageout; + + return 1; + +} Index: /issm/trunk/src/c/modules/HoleFillerx/HoleFillerx.h =================================================================== --- /issm/trunk/src/c/modules/HoleFillerx/HoleFillerx.h (revision 3901) +++ /issm/trunk/src/c/modules/HoleFillerx/HoleFillerx.h (revision 3901) @@ -0,0 +1,16 @@ +/* + HoleFillerx.h +*/ + + +#ifndef HOLEFILLERX_H +#define HOLEFILLERX_H + +//#define _DEBUG_ +//#define _DEBUG2_ + +/* local prototypes: */ +int HoleFillerx(double** pimageout,double* imagein, int lines,int samp,int smooth); + +#endif /* HOLEFILLERX_H */ + Index: /issm/trunk/src/c/modules/InputToResultx/InputToResultx.cpp =================================================================== --- /issm/trunk/src/c/modules/InputToResultx/InputToResultx.cpp (revision 3901) +++ /issm/trunk/src/c/modules/InputToResultx/InputToResultx.cpp (revision 3901) @@ -0,0 +1,33 @@ +/*!\file InputToResultx + * \brief: average field throfieldh thickness + */ + +#include "./InputToResultx.h" + +#include "../shared/shared.h" +#include "../objects/objects.h" +#include "../include/include.h" +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +void InputToResultx(Result** presult,DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials,Parameters* parameters,int enum_type,int id, double time, int step){ + + /*Intermediary output vector*/ + Mat solution=NULL; + Result* result=NULL; + + /*First, get elements*/ + elements->Configure(elements,loads, nodes,vertices, materials,parameters); + + /*Then extrude vertically the new inputs*/ + elements->InputToResult(&solution,enum_type); + + /*Create Result and clean up*/ + /*Mat support in result to be added!!!!!!*/ + //result=new Result(id,time,step,enum_type,solution); + MatFree(&solution); + + /*Assign output pointer*/ + *presult=result; + +} Index: /issm/trunk/src/c/modules/InputToResultx/InputToResultx.h =================================================================== --- /issm/trunk/src/c/modules/InputToResultx/InputToResultx.h (revision 3901) +++ /issm/trunk/src/c/modules/InputToResultx/InputToResultx.h (revision 3901) @@ -0,0 +1,14 @@ +/*!\file: InputToResultx.h + * \brief header file for averaging a field throfieldh thickness + */ + +#ifndef _INPUTTORESULTX_H +#define _INPUTTORESULTX_H + +#include "../DataSet/DataSet.h" + +/* local prototypes: */ +void InputToResultx(Result** presult,DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials,Parameters* parameters,int enum_type,int id, double time, int step); + +#endif /* _INPUTTORESULTX_H */ + Index: /issm/trunk/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp =================================================================== --- /issm/trunk/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp (revision 3901) +++ /issm/trunk/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp (revision 3901) @@ -0,0 +1,90 @@ +/*!\file: InterpFromGridToMeshx.cpp + * \brief "c" core code for interpolating values from a structured grid. + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "./InterpFromGridToMeshx.h" +#include "../shared/shared.h" +#include "../include/include.h" + +int InterpFromGridToMeshx( Vec* pdata_mesh,double* x_in, int x_rows, double* y_in, int y_rows, double* data, int M, int N, double* x_mesh, double* y_mesh, int nods,double default_value) { + + + /*output: */ + Vec data_mesh=NULL; + + /*Intermediary*/ + double* x=NULL; + double* y=NULL; + double x_grid,y_grid; + int i; + + /*threading: */ + InterpFromGridToMeshxThreadStruct gate; + int num=1; + + #ifdef _MULTITHREADING_ + num=_NUMTHREADS_; + #endif + + /*Some checks on arguments: */ + if ((M<=2) || (N<=2) || (nods<=0)){ + ISSMERROR("nothing to be done according to the dimensions of input matrices and vectors."); + } + if (x_in[1]-x_in[0]<0){ + ISSMERROR("x coordinate vectore should be increasing.\n use Matlab's command x=flipud(x), also flip the data matrix data=fliplr(data)"); + } + if (y_in[1]-y_in[0]<0){ + ISSMERROR("y coordinate vectore should be increasing.\n use Matlab's command y=flipud(y), also flip the data matrix data=flipud(data)"); + } + + /*Allocate output vector: */ + data_mesh=NewVec(nods); + + /*Find out what kind of coordinates (x_in,y_in) have been given is input*/ + if(N==(x_rows-1) && M==(y_rows-1)){ + + /*The coordinates given in input describe the contour of each pixel. Take the center of each pixel*/ + x=(double*)xmalloc(N*sizeof(double)); + y=(double*)xmalloc(M*sizeof(double)); + for (i=0;igate; + my_thread=handle->id; + num_threads=handle->num; + + /*recover parameters :*/ + x_mesh=gate->x_mesh; + y_mesh=gate->y_mesh; + x_rows=gate->x_rows; + y_rows=gate->y_rows; + x=gate->x; + y=gate->y; + nods=gate->nods; + data_mesh=gate->data_mesh; + data=gate->data; + default_value=gate->default_value; + M=gate->M; + N=gate->N; + + /*partition loop across threads: */ + PartitionRange(&i0,&i1,nods,num_threads,my_thread); + + /*Linear (triangle) interpolation: */ + for ( i=i0; i