source: issm/branches/trunk-jpl-damage/src/c/modules/HoleFillerx/HoleFillerx.cpp@ 11427

Last change on this file since 11427 was 9313, checked in by Eric.Larour, 14 years ago

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