source: issm/trunk/src/c/HoleFillerx/HoleFillerx.cpp@ 3110

Last change on this file since 3110 was 3110, checked in by Mathieu Morlighem, 15 years ago

no FUNC definition in the header file, always in the c/cpp file
File size: 12.2 KB
Line 
1/*!\file: HoleFillerx.cpp
2 * \brief plug missing values of a double* array, using interpolation techniques. This
3 * is an adaptation of Dr. Crippen's routine (Jet Propulsion Lab, see HoleFillerCrippen.c
4 * for more details.
5 */
6
7#undef __FUNCT__
8#define __FUNCT__ "HoleFillerx"
9
10#include <stdio.h>
11#include <time.h>
12#include "./HoleFillerx.h"
13#include "../shared/shared.h"
14
15int HoleFillerx(double** pimageout,double* image, int lines,int samps,int smooth){
16
17 FILE *fp1;
18 unsigned long filesize;
19 long infinit;
20 double temp;
21 long i, j, k, l, ii, jj, bytes;
22 long iii, jjj;
23 long count, test, prevlines, prevpix, test1, test2, thresh,counter;
24 long bell, cycles, leftovers;
25 float howlong;
26 float nsteps, ssteps, wsteps, esteps;
27 float nwsteps, nesteps, swsteps, sesteps;
28 float north, south, east, west;
29 float nw, sw, ne, se, elev, range;
30 float nnesteps, enesteps, esesteps, ssesteps;
31 float sswsteps, wswsteps, wnwsteps, nnwsteps;
32 float nne, ene, ese, sse;
33 float ssw, wsw, wnw, nnw;
34 float sum;
35
36 time_t t1, t2;
37
38 /*output: */
39 double* imageout=NULL;
40
41 double* image2=NULL;
42 double* image3=NULL;
43 double* image4=NULL;
44
45 /*Added to original routine: */
46 int imageoutsize;
47
48 #ifdef _ISSM_DEBUG_
49 if ( smooth == 1 ){
50 printf("Data patches will be SMOOTHED. \n");
51 }
52 if ( smooth == 0 ){
53 printf("Data patches will NOT be smoothed. \n");
54 }
55 #endif
56
57 /*^^^^^^^^^^^^^ Remove pixels close to the holes ^^^^^^^^^^^^^*/
58 image2 = (double*) xmalloc( lines*samps*sizeof(double));
59 memcpy(image2,image,lines*samps*sizeof(double));
60
61 for ( i = 0; i < lines; i++ ){
62 for ( j = 0; j < samps; j++ ){
63 if ( *(image+i*samps+j) == 0 ){
64 if ( (j > 3) && (j < samps-4) && (i > 3) && (i < lines-4)){
65 for ( k = 0; k < 4; k++ ){
66 for ( l = 0; l < 4; l++ ){
67 *(image2+samps*(i+k)+j+l)=0;
68 *(image2+samps*(i-k)+j+l)=0;
69 *(image2+samps*(i+k)+j-l)=0;
70 *(image2+samps*(i-k)+j-l)=0;
71
72 }
73 }
74 }
75 }
76 }
77 }
78
79 image3 = (double*) xmalloc( lines*samps*sizeof(double));
80 memcpy(image3,image2,lines*samps*sizeof(double));
81
82 count += 2;
83 #ifdef _DEBUG2_
84 counter=0;
85 for ( i = 0; i < lines; i++ ) {
86 for ( j = 0; j < samps; j++ ){
87 if ( *(image2+i*samps+j) != 0 ) {
88 counter++;
89 }
90 }
91 }
92 printf( "\n" );
93 printf("Number of zeroes remaining: %10ld",lines*samps-counter);
94 fflush( stdout );
95 #endif
96
97 #ifdef _ISSM_DEBUG_
98 printf( "\n" );
99 printf( "Iterations %5ld and %5ld.", count-1, count );
100 fflush( stdout );
101 #endif
102 goto afterfirst2;
103
104 again2:
105
106 count += 2;
107
108 #ifdef _DEBUG2_
109 counter=0;
110 for ( i = 0; i < lines; i++ ) {
111 for ( j = 0; j < samps; j++ ){
112 if ( *(image2+i*samps+j) != 0 ) {
113 counter++;
114 }
115 }
116 }
117 // 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
118 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" );
119 printf("Number of zeroes remaining: %10ld",lines*samps-counter);
120 fflush( stdout );
121 #endif
122
123 #ifdef _ISSM_DEBUG_
124 printf( "\b\b\b\b\b\b\b\b\b\b\b\b\b\b" );
125 printf( "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b" );
126 printf( "Iterations %5ld and %5ld.", count-1, count );
127 fflush( stdout );
128 #endif
129
130 afterfirst2:
131
132/***************** FIRST RUN *********************/
133/*
134fprintf ( stdout, "First Application: " ); time(&t2); printf( ctime(&t2) );
135*/
136 for ( i = 0; i < lines; i++ ){
137 for ( j = 0; j < samps; j++ ){
138 if ( *(image2+i*samps+j) != 0 ){
139 *(image3+samps*i+j) = *(image2+i*samps+j);
140 continue;
141 }
142
143 /* Count All Neighbors Test */
144 test = 0;
145 for ( ii = i-1; ii < i+2; ii++ ){
146 for ( jj = j-1; jj < j+2; jj++ ){
147 if (( ii < 0 ) || ( ii >= lines )) continue;
148 if (( jj < 0 ) || ( jj >= samps )) continue;
149 if ( *(image2+ii*samps+jj) != 0 ) test += 1;
150 }
151 }
152 if ( test < 3 ){
153 *(image3+samps*i+j) = 0;
154 continue;
155 }
156
157 /* For void edge pixels: */
158
159 nsteps = 0.0; ssteps = 0.0; esteps = 0.0; wsteps = 0.0;
160 nwsteps = 0.0; nesteps = 0.0; swsteps = 0.0; sesteps = 0.0;
161 north = 0.0; south = 0.0; east = 0.0; west = 0.0;
162 nw = 0.0; ne = 0.0; sw = 0.0; se = 0.0;
163
164 nnesteps = 0.0; enesteps = 0.0; esesteps = 0.0; ssesteps = 0.0;
165 sswsteps = 0.0; wswsteps = 0.0; wnwsteps = 0.0; nnwsteps = 0.0;
166 nne = 0.0; ene = 0.0; ese = 0.0; sse = 0.0;
167 ssw = 0.0; wsw = 0.0; wnw = 0.0; nnw = 0.0;
168
169
170/** NSEW **/
171 for ( ii = i - 1; ii >= 0; ii-- ){ /* North */
172 if ( ii <= 0 ) { north = 0; nsteps = 0; break;}
173 if ( *(image2+ii*samps+j) != 0 ){ north = *(image2+ii*samps+j); nsteps = i - ii; break;}
174 }
175
176 for ( ii = i + 1; ii < lines; ii++ ){ /* South */
177 if ( ii >= lines-1 ){ south = 0; ssteps = 0; break;}
178 if ( *(image2+ii*samps+j) != 0 ){ south = *(image2+ii*samps+j); ssteps = ii - i; break;}
179 }
180
181 for ( jj = j - 1; jj >= 0; jj-- ){ /* West */
182 if ( jj <= 0 ) { west = 0; wsteps = 0; break;}
183 if ( *(image2+i*samps+jj) != 0 ){ west = *(image2+i*samps+jj); wsteps = j - jj; break;}
184 }
185
186 for ( jj = j + 1; jj < samps; jj++ ){ /* East */
187 if ( jj >= samps-1 ){ east = 0; esteps = 0; break;}
188 if ( *(image2+i*samps+jj) != 0 ){ east = *(image2+i*samps+jj); esteps = jj - j; break;}
189 }
190
191
192/** Diagonals **/
193 /* Southeast */
194 for ( k = 1; k < infinit; k++ ){
195 ii = i + k; jj = j + k;
196 if (( ii == lines ) || ( jj == samps )){ se = 0; sesteps = 0; break;}
197 if ( *(image2+ii*samps+jj) != 0 ){ se = *(image2+ii*samps+jj); sesteps = 1.4142 * k; break;}
198 }
199
200 /* Northeast */
201 for ( k = 1; k < infinit; k++ ){
202 ii = i - k; jj = j + k;
203 if (( ii == -1 ) || ( jj == samps )){ ne = 0; nesteps = 0; break;}
204 if ( *(image2+ii*samps+jj) != 0 ){ ne = *(image2+ii*samps+jj); nesteps = 1.4142 * k; break;}
205 }
206
207 /* Northwest */
208 for ( k = 1; k < infinit; k++ ){
209 ii = i - k; jj = j - k;
210 if (( ii == -1 ) || ( jj == -1 )){ nw = 0; nwsteps = 0; break;}
211 if ( *(image2+ii*samps+jj) != 0 ){ nw = *(image2+ii*samps+jj); nwsteps = 1.4142 * k; break;}
212 }
213
214 /* Southwest */
215 for ( k = 1; k < infinit; k++ ){
216 ii = i + k; jj = j - k;
217 if (( ii == lines ) || ( jj == -1 )){ sw = 0; swsteps = 0; break;}
218 if ( *(image2+ii*samps+jj) != 0 ){ sw = *(image2+ii*samps+jj); swsteps = 1.4142 * k; break;}
219 }
220
221/** Other 8 **/
222 /* NNE */
223 for ( k = 1; k < infinit; k++ ){
224 ii = i - (k*2); jj = j + k;
225 if (( ii <= -1 ) || ( jj >= samps )){ nne = 0; nnesteps = 0; break;}
226 if ( *(image2+ii*samps+jj) != 0 ){ nne = *(image2+ii*samps+jj); nnesteps = 2.2361 * k; break;}
227 }
228
229 /* ENE */
230 for ( k = 1; k < infinit; k++ ){
231 ii = i - k; jj = j + (k*2);
232 if (( ii <= -1 ) || ( jj >= samps )){ ene = 0; enesteps = 0; break;}
233 if ( *(image2+ii*samps+jj) != 0 ){ ene = *(image2+ii*samps+jj); enesteps = 2.2361 * k; break;}
234 }
235
236 /* ESE */
237 for ( k = 1; k < infinit; k++ ){
238 ii = i + k; jj = j + (k*2);
239 if (( ii >= lines ) || ( jj >= samps )){ ese = 0; esesteps = 0; break;}
240 if ( *(image2+ii*samps+jj) != 0 ){ ese = *(image2+ii*samps+jj); esesteps = 2.2361 * k; break;}
241 }
242
243 /* SSE */
244 for ( k = 1; k < infinit; k++ ){
245 ii = i + (k*2); jj = j + k;
246 if (( ii >= lines ) || ( jj >= samps )){ sse = 0; ssesteps = 0; break;}
247 if ( *(image2+ii*samps+jj) != 0 ){ sse = *(image2+ii*samps+jj); ssesteps = 2.2361 * k; break;}
248 }
249
250 /* SSW */
251 for ( k = 1; k < infinit; k++ ){
252 ii = i + (k*2); jj = j - k;
253 if (( ii >= lines ) || ( jj <= -1 )){ ssw = 0; sswsteps = 0; break;}
254 if ( *(image2+ii*samps+jj) != 0 ){ ssw = *(image2+ii*samps+jj); sswsteps = 2.2361 * k; break;}
255 }
256
257 /* WSW */
258 for ( k = 1; k < infinit; k++ ){
259 ii = i + k; jj = j - (k*2);
260 if (( ii >= lines ) || ( jj <= -1 )){ wsw = 0; wswsteps = 0; break;}
261 if ( *(image2+ii*samps+jj) != 0 ){ wsw = *(image2+ii*samps+jj); wswsteps = 2.2361 * k; break;}
262 }
263
264 /* WNW */
265 for ( k = 1; k < infinit; k++ ){
266 ii = i - k; jj = j - (k*2);
267 if (( ii <= -1 ) || ( jj <= -1 )){ wnw = 0; wnwsteps = 0; break;}
268 if ( *(image2+ii*samps+jj) != 0 ){ wnw = *(image2+ii*samps+jj); wnwsteps = 2.2361 * k; break;}
269 }
270
271 /* NNW */
272 for ( k = 1; k < infinit; k++ ){
273 ii = i - (k*2); jj = j - k;
274 if (( ii <= -1 ) || ( jj <= -1 )){ nnw = 0; nnwsteps = 0; break;}
275 if ( *(image2+ii*samps+jj) != 0 ){ nnw = *(image2+ii*samps+jj); nnwsteps = 2.2361 * k; break;}
276 }
277
278 elev = 0; range = 0;
279/*NSEW*/
280 if ( nsteps > 0.5 ){ elev += north / nsteps; range += 1.00 / nsteps;}
281 if ( ssteps > 0.5 ){ elev += south / ssteps; range += 1.00 / ssteps;}
282 if ( esteps > 0.5 ){ elev += east / esteps; range += 1.00 / esteps;}
283 if ( wsteps > 0.5 ){ elev += west / wsteps; range += 1.00 / wsteps;}
284/*Diag*/
285 if ( nwsteps > 0.5 ){ elev += nw / nwsteps; range += 1.00 / nwsteps;}
286 if ( nesteps > 0.5 ){ elev += ne / nesteps; range += 1.00 / nesteps;}
287 if ( swsteps > 0.5 ){ elev += sw / swsteps; range += 1.00 / swsteps;}
288 if ( sesteps > 0.5 ){ elev += se / sesteps; range += 1.00 / sesteps;}
289
290/*Other 8*/
291 if ( nnesteps > 0.5 ){ elev += nne / nnesteps; range += 1.00 / nnesteps;}
292 if ( enesteps > 0.5 ){ elev += ene / enesteps; range += 1.00 / enesteps;}
293 if ( esesteps > 0.5 ){ elev += ese / esesteps; range += 1.00 / esesteps;}
294 if ( ssesteps > 0.5 ){ elev += sse / ssesteps; range += 1.00 / ssesteps;}
295 if ( sswsteps > 0.5 ){ elev += ssw / sswsteps; range += 1.00 / sswsteps;}
296 if ( wswsteps > 0.5 ){ elev += wsw / wswsteps; range += 1.00 / wswsteps;}
297 if ( wnwsteps > 0.5 ){ elev += wnw / wnwsteps; range += 1.00 / wnwsteps;}
298 if ( nnwsteps > 0.5 ){ elev += nnw / nnwsteps; range += 1.00 / nnwsteps;}
299
300 //temp = ( elev / range ) + 0.5 ;
301 temp = ( elev / range );
302 //+ 0.5 ;
303 //if ( temp > 10000 ) temp = 10000;
304 //if ( temp < 0 ) temp = 0;
305
306 #ifdef _DEBUG2_
307 //printf("%g %g %g \n",temp,elev,range);
308 #endif
309
310 *(image3+samps*i+j) = temp ;
311 }
312 }
313
314
315
316 for ( i = 0; i < lines; i++ ){
317 for ( j = 0; j < samps; j++ ){
318 *(image2+samps*i+j)= *(image3+samps*i+j);
319 }
320 }
321
322
323 for ( i = 0; i < lines; i++ ){
324 for ( j = 0; j < samps; j++ ){
325 if ( *(image2+samps*i+j) == 0 ) goto again2;
326 }
327 }
328
329 if ( smooth == 0 ) goto there2;
330
331
332/************************ SMOOTH THE RESULT ***********************/
333
334 image4 = (double*) xmalloc( lines*samps*sizeof(double));
335 memcpy(image4,image3,lines*samps*sizeof(double));
336
337
338 for ( i = 0; i < lines; i++ ) {
339 for ( j = 0; j < samps; j++ ) {
340 if ( *(image4+i*samps+j) != 0 ) { *(image3+i*samps+j) = *(image2+i*samps+j) ; continue; }
341
342
343 for ( k = 1; k < infinit; k++ ) { /* Find the smallest box size with data */
344 for ( ii = i-k; ii <= i+k; ii++ ){
345 if ( ii < 0 ) continue;
346 if ( ii >= lines ) break;
347 for ( jj = j-k; jj <= j+k; jj++ ) {
348 if ( jj < 0 ) continue;
349 if ( jj >= samps ) continue;
350 if ( *(image4+ii*samps+jj) != 0 ) goto k_nowset;
351 }
352 }
353 }
354
355 k_nowset:
356 k = k / 4; if ( k < 1 ) k = 1; /* Errrrr. Make it fourth size */
357 sum = 0;
358 for ( ii = i-k; ii <= i+k; ii++ ) {
359 if ( ii < 0 ) iii = -1 - ii; else iii = ii;
360 if ( ii >= lines ) iii = lines - (ii - lines + 1);
361 for ( jj = j-k; jj <= j+k; jj++ ) {
362 if ( jj < 0 ) jjj = -1 - jj; else jjj = jj;
363 if ( jj >= samps ) jjj = samps - (jj - samps + 1);
364 sum += *(image2+samps*iii+jjj) ;
365 }
366 }
367 k = (k*2)+1;
368 *(image3+samps*i+j) = ( sum / k / k ) ;//+ 0.5;
369 }
370
371 #ifdef _ISSM_DEBUG_
372 if ( !(i%100) ){
373 fprintf(stdout, "\b\b\b\b\b\b%5ld ", i ); /* Count lines on output screen */
374 fflush(stdout);
375 }
376 #endif
377 }
378
379 #ifdef _ISSM_DEBUG_
380 fprintf(stdout, "\b\b\b\b\b\b%5ld ", lines );
381 fflush(stdout);
382 #endif
383
384
385there2:
386
387
388 /*Allocate output image: */
389 imageout=(double*)xmalloc(samps*lines*sizeof(double));
390 memcpy(imageout,image3,lines*samps*sizeof(double));
391
392 time(&t2);
393 #ifdef _ISSM_DEBUG_
394 printf ( "\n\nEnd "); printf( ctime(&t2) );
395 printf( "\n Total time: ");
396 printf( "%4.2f minutes. \n", difftime(t2,t1)/60 );
397 #endif
398
399 #ifdef _DEBUG2_
400 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" );
401 printf("Number of zeroes remaining: 0\n\n");
402 printf ( "\n");
403 #endif
404
405 end:
406
407 #ifdef _ISSM_DEBUG_
408 printf( "Done.\a\n" );
409 #endif
410
411 /*Assign output pointers: */
412 *pimageout=imageout;
413
414 return 1;
415
416}
Note: See TracBrowser for help on using the repository browser.