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

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

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