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

Last change on this file since 1904 was 1904, checked in by Eric.Larour, 16 years ago

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