| 1 |
/* code pulled from ana for a standalone despike for radiation hits */ |
| 2 |
#include <stdio.h> |
| 3 |
#include <stdlib.h> |
| 4 |
#include <math.h> |
| 5 |
#include <string.h> |
| 6 |
#include <sys/types.h> |
| 7 |
#include <sys/time.h> |
| 8 |
#define ALN2I 1.442695022 |
| 9 |
#define TINY 1.0e-5 |
| 10 |
#define ABS(x) ((x)>=0?(x):-(x)) |
| 11 |
// double systime(); /* note that systime is defined elsewhere in the ana environment, remove |
| 12 |
// all refs for standalone (or add func systime) */ |
| 13 |
/* 3/17/2010 - interface change to I*4, also work on locations to remove duplicates (some |
| 14 |
pixels are despiked twice when iterated) and do not use output array, modify the input array */ |
| 15 |
/* 9/29/2010 - some mods to remove negative values from the data */ |
| 16 |
/* 11/20/10 - now want to keep neg values, so remove those mods and make sure |
| 17 |
we don't get too many extra spikes by adjusting absmin (from 1 to 4), we |
| 18 |
already had a test to avoid despiking neg pixels from 9/29/2010 mods */ |
| 19 |
/*---------------------------------------------------------------------------*/ |
| 20 |
int aia_despike( |
| 21 |
int *array, unsigned char *mask, int nx, int ny, float frac, int level, int niter, /* void **outarray, */ |
| 22 |
int *badblobs, int sizeofbads, |
| 23 |
int *nspikes, int **oldvalues, int **spikelocs, int **newvalues); |
| 24 |
int mask_erode(unsigned char *pbase, unsigned char *qbase, int n, int m); /* erosion function */ |
| 25 |
/*------------------------------------------------------------------------- */ |
| 26 |
int aia_despike( /* despike function, modified heavily for AIA */ |
| 27 |
int *array, unsigned char *mask, int nx, int ny, float frac, int level, int niter, /* void **outarray, */ |
| 28 |
int *badblobs, int sizeofbads, |
| 29 |
int *nspikes, int **oldvalues, int **spikelocs, int **newvalues) |
| 30 |
/* the call is |
| 31 |
int *ptr, *oldvalues, *newvalues; |
| 32 |
int nx, ny, level=7, niter=3, *spikelocs; |
| 33 |
float frac = 0.5; |
| 34 |
unsigned char *mask; |
| 35 |
|
| 36 |
stat = despike(ptr, mask, nx, ny, frac, level, niter, |
| 37 |
&nspikes, &oldvalues, &spikelocs, &newvalues); |
| 38 |
|
| 39 |
note that a NULL passed for mask implies full image (a fake mask is generated internally) |
| 40 |
a 0 for sizeofbads implies no bad pixel correction |
| 41 |
|
| 42 |
note that the nx by ny I*4 array in ptr will be despiked thereby modifying the original, |
| 43 |
the original can be recovered using the offsets in spikelocs and the original values |
| 44 |
in oldvalues, the new values are also available in the newvalues array, these 3 arrays |
| 45 |
are malloc'ed here and the caller is responsible for free'ing them |
| 46 |
|
| 47 |
The bad pixels are represented as bad "blobs" to allow for possible clusters. The badblobs |
| 48 |
array is a null terminated packed list of these, the size of the list must also be |
| 49 |
passed in sizeofbads to avoid any overflows (and allow for a quick check for none), |
| 50 |
the format of the list is [nperim, perim_offsets, npix, pix_offsets] in series for |
| 51 |
each "blob", nperim is the # of perimeter points (offsets in ptr) to use for averaging |
| 52 |
to obtain a value to substitute for the bad pixels. perim_offsets is the list |
| 53 |
of nperim offsets. npix is the count of bad pixels to be replaced by this average and |
| 54 |
pix_offsets is the npix offsets. Note that all bad pixels are replaced by the same value |
| 55 |
(no gradients or such) computed from the average of the perimeter. The length of each |
| 56 |
in the the list is nperim + npix + 2. The code just steps through. |
| 57 |
*/ |
| 58 |
{ |
| 59 |
int iq, result_sym, n, m, sum, nc; |
| 60 |
int lognb2, jj, jc; |
| 61 |
float absmin = 4.0; |
| 62 |
int nxc, nyc, offset, *lastiterendaddr; |
| 63 |
int save_niter, ntotal, npix, nslocal, ndups, itercount; |
| 64 |
float fsum, cfrac, tq, fdif; |
| 65 |
int *p, *ptr, *p1, *p2, *p3, *spikeold, *oldptr; |
| 66 |
unsigned char *mptr, *mp, mskv, *mskbase, *eroded, *eroded2; |
| 67 |
int *spikebufadd, *sptradd; |
| 68 |
int *spikenew, *newptr, *base; |
| 69 |
int arr[16], *pps, *ss, *spikestarts[20], *spikeends[20]; |
| 70 |
//double t1, t2, t3, t4; |
| 71 |
//t1 = systime(); |
| 72 |
lognb2 = (log((double) 16)*ALN2I+TINY); /* needed for median sort */ |
| 73 |
|
| 74 |
/* always an int array for this routine */ |
| 75 |
if (nx < 5 || ny < 5 ) { |
| 76 |
printf("dimensions must be 5 or greater, nx, ny = %d %d\n", nx,ny); |
| 77 |
return -1; } |
| 78 |
/* bad pixel correction (if any), uses blob list as noted above */ |
| 79 |
if (sizeofbads > 0 && badblobs != 0 ) { |
| 80 |
/* correction done in place, if badblobs not defined correctly, a memory access violation |
| 81 |
could occur */ |
| 82 |
int countdown = sizeofbads; |
| 83 |
int *pblob, nperim, npix, xq, offset; |
| 84 |
float acc; |
| 85 |
pblob = badblobs; |
| 86 |
while (1) { |
| 87 |
acc = 0.0; |
| 88 |
nperim = *pblob++; |
| 89 |
if (nperim <= 0 || nperim > countdown) break; |
| 90 |
n = nperim; |
| 91 |
/* this takes the average of the perimeter values, a median might be better */ |
| 92 |
while (n--) { offset = *pblob++; acc += *(array + offset); } |
| 93 |
xq = rintf(acc/(float) nperim); |
| 94 |
npix = *pblob++; |
| 95 |
if (npix <= 0 || npix > countdown) break; |
| 96 |
n = npix; |
| 97 |
while (n--) { offset = *pblob++; *(array + offset) = xq; } |
| 98 |
countdown = countdown - nperim - npix - 2; |
| 99 |
if (countdown <= 0) break; |
| 100 |
} |
| 101 |
} |
| 102 |
|
| 103 |
/* flip image */ |
| 104 |
{ |
| 105 |
int i, *toprow, *botrow, *tmprow; |
| 106 |
botrow = array; |
| 107 |
toprow = botrow + 4095*4096; |
| 108 |
tmprow = (int *) malloc(sizeof(int)*4096); |
| 109 |
for (i = 0; i < 2048; i++) { |
| 110 |
memcpy(tmprow, toprow, 16384); |
| 111 |
memcpy(toprow, botrow, 16384); |
| 112 |
memcpy(botrow, tmprow, 16384); |
| 113 |
botrow += 4096; |
| 114 |
toprow -= 4096; |
| 115 |
} |
| 116 |
free(tmprow); |
| 117 |
} |
| 118 |
/* if no mask, we actually make one full of 1's, partly for testing */ |
| 119 |
eroded = malloc(nx*ny*sizeof(unsigned char)); |
| 120 |
if (!eroded) { printf("malloc error in local mask copy\n"); return 1; } |
| 121 |
mp = eroded; |
| 122 |
/* unfortunately, we need to check array for 0x80000000 values and avoid using |
| 123 |
these for despike decisions and we must not modify them. This is done by zeroing |
| 124 |
the mask for the corresponding points before erosion */ |
| 125 |
p = array; |
| 126 |
/* 9/29/2010 - also catch any other negative values while we are looping and zero them */ |
| 127 |
if (mask == 0) { |
| 128 |
int nq = nx * ny; |
| 129 |
while (nq--) { |
| 130 |
/* 9/29/2010 - change here to remove negative values (convert to 0) */ |
| 131 |
/* 11/20/10 - try leaving the negative values in */ |
| 132 |
//if (*p == 0x80000000) { *mp++ = 0; } else { *mp++ = 1; if (*p < 0) *p = 0; } |
| 133 |
if (*p == 0x80000000) { *mp++ = 0; } else { *mp++ = 1; } |
| 134 |
p++; |
| 135 |
} |
| 136 |
/* mask gets set to "eroded" further down */ |
| 137 |
} else { |
| 138 |
/* have a real mask, need to erode it after checking for 0x80000000's */ |
| 139 |
int nq = nx * ny; |
| 140 |
mp = mask; |
| 141 |
while (nq--) { |
| 142 |
/* 9/29/2010 - change here to remove negative values (convert to 0) */ |
| 143 |
/* 11/20/10 - try leaving the negative values in */ |
| 144 |
//if (*p == 0x80000000) { *mp++ = 0; } else { if (*p < 0) *p = 0; } |
| 145 |
if (*p == 0x80000000) { *mp++ = 0; } |
| 146 |
p++; |
| 147 |
} |
| 148 |
eroded2 = malloc(nx*ny*sizeof(unsigned char)); |
| 149 |
if (!eroded2) { printf("malloc error in local mask copy\n"); return 1; } |
| 150 |
mask_erode(mask, eroded2, nx, ny); |
| 151 |
mask_erode(eroded2, eroded, nx, ny); |
| 152 |
/* leaves final result in eroded, so free eroded2 now and redefine mask */ |
| 153 |
free(eroded2); |
| 154 |
} |
| 155 |
mask = mp = eroded; |
| 156 |
ptr = base = array; |
| 157 |
/* and the mask */ |
| 158 |
mptr = mskbase = mask; |
| 159 |
cfrac = 1.0/(1.0 + frac); /* 3/2/2010 changed from 1.0 - frac which had less range */ |
| 160 |
nc = ntotal = nslocal = 0; |
| 161 |
niter = ABS(niter); |
| 162 |
/* 4/7/2010 - add a check for niter < 0 to avoid runaways */ |
| 163 |
if (niter > 20 || niter < 0) { printf("DESPIKE - error, excessive # of iterations = %d\n", |
| 164 |
niter); return 1; } |
| 165 |
|
| 166 |
/* add internal iteration 10/8/98 */ |
| 167 |
save_niter = niter; |
| 168 |
spikeold = oldptr = malloc(nx*ny*sizeof(int)); /* that allows every pixel */ |
| 169 |
spikenew = newptr = malloc(nx*ny*sizeof(int)); /* that allows every pixel */ |
| 170 |
spikebufadd = sptradd = malloc(nx*ny*sizeof(int)); /* that allows every pixel */ |
| 171 |
lastiterendaddr = spikebufadd; |
| 172 |
ndups = 0; |
| 173 |
//t2 = systime(); |
| 174 |
npix = 0; |
| 175 |
itercount = 0; |
| 176 |
spikestarts[itercount] = spikebufadd; |
| 177 |
/* note that niter = 0 should not do any despike */ |
| 178 |
while (niter--) { |
| 179 |
ptr = array + 2*nx; |
| 180 |
/* and the mask */ |
| 181 |
mptr = mask + 2*nx; |
| 182 |
/* skip the outer edges */ |
| 183 |
m = ny-4; |
| 184 |
while (m--) { |
| 185 |
/* skip the outer edges */ |
| 186 |
p = ptr; |
| 187 |
p += 2; |
| 188 |
/* and the mask */ |
| 189 |
mp = mptr + 2; |
| 190 |
n = nx-4; |
| 191 |
p2 = p - 1; |
| 192 |
p1 = p2 - nx; |
| 193 |
p3 = p2 + nx; |
| 194 |
while (n--) { |
| 195 |
npix++; |
| 196 |
/* watch out for 0's, lots in test images */ |
| 197 |
/* 3/22/2010 - mask check finally added for each position */ |
| 198 |
/* 9/29/2010 - avoid processing negative values as well as zeroes */ |
| 199 |
if ( (*p > 0) && *mp ) { |
| 200 |
/* add the 8 around this point */ |
| 201 |
tq = (cfrac * (float) *p); |
| 202 |
sum = *p1 + *(p1+1) + *(p1+2) + *p2 + *(p2+2) + *p3 + *(p3+1) + *(p3+2); |
| 203 |
/* note the absmin term, this is to avoid problems with data that has large |
| 204 |
zones of zeroes, we then get too many hits and erosion of edges, mostly a problem |
| 205 |
with artifical data */ |
| 206 |
fsum = (float) sum * 0.125; |
| 207 |
/* now the test */ |
| 208 |
if ( (fsum < tq) && ((tq-fsum) > absmin) ) { /* we have a bad one, zap it */ |
| 209 |
nc++; |
| 210 |
/* load up sort array and find the desired one */ |
| 211 |
ss = arr; pps = p - 2*nx -2; |
| 212 |
*ss++ = *pps++; *ss++ = *pps++; *ss++ = *pps++; *ss++ = *pps++; *ss++ = *pps++; |
| 213 |
*ss++ = *(p - nx -2); *ss++ = *(p - nx +2); |
| 214 |
*ss++ = *(p -2); *ss++ = *(p +2); |
| 215 |
*ss++ = *(p + nx -2); *ss++ = *(p + nx +2); |
| 216 |
pps = p + 2 *nx - 2; |
| 217 |
*ss++ = *pps++; *ss++ = *pps++; *ss++ = *pps++; *ss++ = *pps++; *ss++ = *pps++; |
| 218 |
/* this is basically a median filter (adjustable using level value) to find |
| 219 |
a "good" value from the surroundings, we do a median rather than an average to |
| 220 |
try to avoid some effects of spike overflow on the CCD */ |
| 221 |
/* a built in sorter here */ |
| 222 |
{ int nn,m,j,i,n=16; |
| 223 |
int t; |
| 224 |
m=n; |
| 225 |
for (nn=1;nn<=lognb2;nn++) { |
| 226 |
m >>= 1; |
| 227 |
for (j=m;j<n;j++) { |
| 228 |
i=j-m; |
| 229 |
t=arr[j]; |
| 230 |
while (i >= 0 && arr[i] > t) { |
| 231 |
arr[i+m]=arr[i]; |
| 232 |
i -= m; |
| 233 |
} |
| 234 |
arr[i+m]=t; |
| 235 |
} |
| 236 |
} |
| 237 |
} |
| 238 |
/* 3/11/2010 - log into various arrays */ |
| 239 |
*oldptr++ = *p; offset = p - base; |
| 240 |
//if (nslocal == 0) printf("sptradd, offset = %d, %d\n", sptradd, offset); |
| 241 |
*sptradd++ = offset; |
| 242 |
nslocal++; /* yet another spike counter */ |
| 243 |
/* now get the indicated one using level, also save into spikenew */ |
| 244 |
*newptr++ = arr[level]; |
| 245 |
} |
| 246 |
} |
| 247 |
p++; mp++; p1++; p2++; p3++; |
| 248 |
|
| 249 |
} |
| 250 |
/* skip the outer edges */ |
| 251 |
//p = p + 2; |
| 252 |
ptr = ptr + nx; |
| 253 |
/* and the mask */ |
| 254 |
//mp = mp + 2; |
| 255 |
mptr = mptr + nx; |
| 256 |
} |
| 257 |
/* skip the last 2 rows */ |
| 258 |
/* end of this iteration, reconfigure for next one, note that for |
| 259 |
a single iteration we don't need out2 */ |
| 260 |
// printf("despike got %d spikes in iteration %d out of %d pixels\n", nc, itercount, npix); |
| 261 |
// printf("nslocal = %d\n", nslocal); |
| 262 |
ntotal += nc; nc = 0; |
| 263 |
npix = 0; |
| 264 |
if (niter < (save_niter - 1)) { |
| 265 |
int *ip, *jp, *iplast = spikebufadd; |
| 266 |
int ioff, joff, delta1, delta2, previter; |
| 267 |
/* check the range spikebufadd to lastiterendaddr against lastiterendaddr+1 to sptradd-1, |
| 268 |
the idea here is to check if any of our new spikes have the same location from a previous |
| 269 |
iteration, this isn't a n*m operation because each set is monotonic, note however that |
| 270 |
we have to scan each iteration set separately (which is annoying) */ |
| 271 |
int *jprevstart[20]; |
| 272 |
spikestarts[itercount] = spikeends[itercount-1] + 1; |
| 273 |
for (previter=0; previter<itercount; previter++) { jprevstart[previter] = spikestarts[previter]; } |
| 274 |
for (ip=spikestarts[itercount]; ip<sptradd; ip++) { |
| 275 |
ioff = *ip; |
| 276 |
/* now compare with all the previous ones */ |
| 277 |
for (previter=0; previter<itercount; previter++) { |
| 278 |
|
| 279 |
for (jp=jprevstart[previter]; jp<=spikeends[previter]; jp++) { |
| 280 |
joff = *jp; |
| 281 |
if (joff == ioff) { |
| 282 |
/* we want the new instance of this location to be the one used, it needs |
| 283 |
the earlier old value, it already has the latest new value */ |
| 284 |
delta2 = jp - spikebufadd; |
| 285 |
delta1 = ip - spikebufadd; |
| 286 |
*(spikeold + delta1) = *(spikeold + delta2); |
| 287 |
/* and zap the original location */ |
| 288 |
*jp = -1; |
| 289 |
ndups++; |
| 290 |
break; |
| 291 |
} |
| 292 |
if (joff > ioff) { |
| 293 |
/* because both are monotonic, can't be any more now */ |
| 294 |
|
| 295 |
break; |
| 296 |
} |
| 297 |
} |
| 298 |
/* reset the start to avoid wasting time on check for next jp since both are monotonic */ |
| 299 |
jprevstart[previter] = jp; |
| 300 |
} |
| 301 |
} |
| 302 |
//printf("# of dups = %d\n", ndups); |
| 303 |
} |
| 304 |
spikeends[itercount] = sptradd - 1; |
| 305 |
// printf("setting spikestarts[itercount], spikeends[itercount] = %d, %d\n", spikestarts[itercount], |
| 306 |
// spikeends[itercount]); |
| 307 |
/* and we now set the new values from this iteration in the input array */ |
| 308 |
{ |
| 309 |
int *ip, ioff, *snew, nq; |
| 310 |
ip=spikestarts[itercount]; |
| 311 |
snew = spikenew + (ip - spikebufadd); |
| 312 |
while ( ip<sptradd ) { |
| 313 |
ioff = *ip++; |
| 314 |
array[ioff] = *snew++; |
| 315 |
} |
| 316 |
} |
| 317 |
itercount++; |
| 318 |
} |
| 319 |
/* wrap up the spike log */ |
| 320 |
*nspikes = nslocal; |
| 321 |
// printf("# of spikes logged = %d, dups = %d\n", nslocal, ndups); |
| 322 |
*nspikes = nslocal - ndups; |
| 323 |
// printf("# of unique spikes logged = %d\n", *nspikes); |
| 324 |
n = (nslocal - ndups)*sizeof(int); |
| 325 |
/* the duplicates are tagged, copy only the unique ones for export */ |
| 326 |
if (n > 0) { |
| 327 |
int i; |
| 328 |
int *p1, *p2, *p3, *q1, *q2, *q3; |
| 329 |
*oldvalues = p1 = malloc(n); |
| 330 |
q1 = spikeold; |
| 331 |
*newvalues = p2 = malloc(n); |
| 332 |
q2 = spikenew; |
| 333 |
*spikelocs = p3 = malloc(n); |
| 334 |
q3 = spikebufadd; |
| 335 |
for (i=0;i<nslocal;i++) { |
| 336 |
if (*q3 > 0) { *p1++ = *q1; *p2++ = *q2; *p3++ = *q3; } |
| 337 |
q1++; q2++; q3++; |
| 338 |
} |
| 339 |
} else { |
| 340 |
*oldvalues = 0; *spikelocs = 0; *newvalues = 0; |
| 341 |
} |
| 342 |
/* free up (reverse order) */ |
| 343 |
free(spikebufadd); free(spikenew); free(spikeold); free(eroded); |
| 344 |
|
| 345 |
// printf("despike got %d spikes in %d iterations\n", ntotal, save_niter); |
| 346 |
// t3 = systime(); |
| 347 |
|
| 348 |
// printf("despike times, setup %10.6fs, iters %10.6fs, total %10.6fs\n", t2-t1,t3-t2,t3-t1); |
| 349 |
return 0; |
| 350 |
} |
| 351 |
/*------------------------------------------------------------------------- */ |
| 352 |
/* a standalone erode function specialized for eroding the AIA masks prior to |
| 353 |
despiking */ |
| 354 |
/*------------------------------------------------------------------------- */ |
| 355 |
int mask_erode(unsigned char *pbase, unsigned char *qbase, int n, int m) /* erosion function */ |
| 356 |
{ |
| 357 |
unsigned char *p, *q, *p_endline, *qabove, *qbelow; |
| 358 |
int iq, mq, type; |
| 359 |
double t1, t2, t3, t4; |
| 360 |
// t1 = systime(); |
| 361 |
if ( n < 3 || m < 3 ) { |
| 362 |
printf("ERODE: array too small, nx, ny = %d %d\n", n, m); |
| 363 |
return -1; } |
| 364 |
|
| 365 |
bcopy(pbase, qbase, n*m); |
| 366 |
p = pbase; |
| 367 |
q = qbase; |
| 368 |
|
| 369 |
/* the edges and corners are done as special cases */ |
| 370 |
/* first corner */ |
| 371 |
if (*p == 0) { /* zap the 3 neighbors */ |
| 372 |
*q++; *q = 0; q = qbase + n; *q++ = 0; *q = 0;} |
| 373 |
/* starting row */ |
| 374 |
p_endline = pbase + n - 2; |
| 375 |
p++; |
| 376 |
// t2 = systime(); |
| 377 |
while (1) { |
| 378 |
if (*p == 0) { |
| 379 |
/* got a hit, this means we need to clear 6 pixels in the output image */ |
| 380 |
q = (qbase - pbase - 1) + p; /* this is the q just before the corresponding p */ |
| 381 |
qabove = q + n; |
| 382 |
*q++ = 0; q++; *q++ = 0; |
| 383 |
*qabove++ = 0; *qabove++ = 0; *qabove++ = 0; |
| 384 |
if (p >= p_endline) break; |
| 385 |
p++; |
| 386 |
/* if we continue to get consecutive hits, just need to set next 2 */ |
| 387 |
while (*p == 0) { |
| 388 |
if (p >= p_endline) break; |
| 389 |
*q++ = 0; *qabove++ = 0; |
| 390 |
p++; } |
| 391 |
} |
| 392 |
if (p >= p_endline) break; |
| 393 |
p++; |
| 394 |
} |
| 395 |
p++; |
| 396 |
/* last point in starting row, set independent of previous so we may |
| 397 |
set some pixels twice */ |
| 398 |
if (*p == 0) { |
| 399 |
q = (qbase - pbase - 1) + p; /* this is the q just before the corresponding p */ |
| 400 |
qabove = q + n; |
| 401 |
*q = 0; |
| 402 |
*qabove++ = 0; *qabove = 0; |
| 403 |
} |
| 404 |
p++; |
| 405 |
|
| 406 |
|
| 407 |
/* central area */ |
| 408 |
/* now the interior rows */ |
| 409 |
mq = m - 2; |
| 410 |
while (mq-- > 0) { |
| 411 |
/* start row, not top or bottom */ |
| 412 |
/* left hand edge */ |
| 413 |
if (*p == 0) { /* set 6 points */ |
| 414 |
q = (qbase - pbase - 1) + 1 + p; /* q at p */ |
| 415 |
qabove = q + n; qbelow = q - n; |
| 416 |
*q++; *q = 0; |
| 417 |
*qabove++ = 0; *qabove = 0; |
| 418 |
*qbelow++ = 0; *qbelow = 0; |
| 419 |
} |
| 420 |
p_endline = p + n - 2; |
| 421 |
p++; |
| 422 |
|
| 423 |
/* done with left edge, now the middle */ |
| 424 |
|
| 425 |
while (1) { |
| 426 |
if (*p == 0) { |
| 427 |
/* got a hit, this means we need to clear 8 pixels in the output image */ |
| 428 |
q = (qbase - pbase - 1) + p; |
| 429 |
qabove = q + n; qbelow = q - n; |
| 430 |
*q++ = 0; q++; *q++ = 0; |
| 431 |
*qabove++ = 0; *qabove++ = 0; *qabove++ = 0; |
| 432 |
*qbelow++ = 0; *qbelow++ = 0; *qbelow++ = 0; |
| 433 |
if (p >= p_endline) break; |
| 434 |
p++; |
| 435 |
/* if we continue to get consecutive hits, just need to clear next 3 */ |
| 436 |
while (*p == 0) { |
| 437 |
if (p >= p_endline) break; |
| 438 |
*q++ = 0; *qabove++ = 0; *qbelow++ = 0; |
| 439 |
p++; } |
| 440 |
} |
| 441 |
if (p >= p_endline) break; |
| 442 |
p++; |
| 443 |
} |
| 444 |
p++; |
| 445 |
|
| 446 |
|
| 447 |
/* the last point in row */ |
| 448 |
if (*p == 0) { |
| 449 |
q = (qbase - pbase - 1) + p; /* this is the q just before the corresponding p */ |
| 450 |
qabove = q + n; qbelow = q - n; |
| 451 |
*q = 0; |
| 452 |
*qabove++ = 0; *qabove++ = 0; |
| 453 |
*qbelow++ = 0; *qbelow++ = 0; |
| 454 |
} |
| 455 |
p++; |
| 456 |
|
| 457 |
} |
| 458 |
|
| 459 |
/* at last, the last row */ |
| 460 |
/* left hand edge */ |
| 461 |
if (*p == 0) { /* set 4 points */ |
| 462 |
q = (qbase - pbase - 1) + 1 + p; /* q at p */ |
| 463 |
qbelow = q - n; |
| 464 |
q++; *q = 0; |
| 465 |
*qbelow++ = 0; *qbelow = 0; |
| 466 |
} |
| 467 |
p_endline = p + n - 2; |
| 468 |
p++; |
| 469 |
|
| 470 |
/* now the middle of the last row */ |
| 471 |
|
| 472 |
while (1) { |
| 473 |
if (*p == 0) { |
| 474 |
/* got a hit, this means we need to clear 9 pixels in the output image */ |
| 475 |
q = (qbase - pbase - 1) + p; |
| 476 |
qbelow = q - n; |
| 477 |
*q++ = 0; q++; *q++ = 0; |
| 478 |
*qbelow++ = 0; *qbelow++ = 0; *qbelow++ = 0; |
| 479 |
if (p >= p_endline) break; |
| 480 |
p++; |
| 481 |
/* if we continue to get consecutive hits, just need to set next 3 */ |
| 482 |
while (*p == 0) { |
| 483 |
if (p >= p_endline) break; |
| 484 |
*q++ = 0; *qbelow++ = 0; p++; } |
| 485 |
} |
| 486 |
if (p >= p_endline) break; |
| 487 |
p++; |
| 488 |
} |
| 489 |
p++; |
| 490 |
/* printf("pbase, p_endline, p = %#x, %#x, %#x\n", pbase, p_endline, p); */ |
| 491 |
/* the last point in last row */ |
| 492 |
/* printf("pbase, p = %#x, %#x\n", pbase, p); */ |
| 493 |
if (*p == 0) { |
| 494 |
q = (qbase - pbase - 1) + p; /* this is the q just before the corresponding p */ |
| 495 |
qbelow = q - n; |
| 496 |
*q = 0; |
| 497 |
*qbelow++ = 0; *qbelow = 0; |
| 498 |
} |
| 499 |
// t3 = systime(); |
| 500 |
// printf("AIA erode total time = %10.2f ms, for setup = %10.2f ms, part 2 = %10.2f ms\n", |
| 501 |
// 1.E3*(t3-t1), 1.E3*(t2-t1), 1.E3*(t3-t2)); |
| 502 |
return 0; |
| 503 |
} |
