OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
filter.c
Go to the documentation of this file.
1 /* ---------------------------------------------------------- */
2 /* filter.c - Filtering control module for MSl12 */
3 /* */
4 /* Written By: B. Franz, SAIC GSC, NASA/SIMBIOS, October 1998 */
5 /* ---------------------------------------------------------- */
6 
7 #include "l12_proto.h"
8 #include <ctype.h>
9 #include <stdint.h>
10 
11 static float badval = BAD_FLT;
12 
13 typedef struct fnode_str {
14  int32_t i;
15  int32_t j;
16  float v;
17 } fnode;
18 
19 void fctl_init(fctlstr *fctl) {
20  fctl->nfilt = 0;
21  fctl->nscan = 1;
22  fctl->npix = 1;
23 }
24 
25 void get_kernel(filstr *f) {
26  int type = 0;
27  int32_t i;
28  int32_t j;
29  int32_t nnx;
30  int32_t nx = f->nx;
31  int32_t ny = f->ny;
32 
33  if (nx == -1) {
34 
35  type = 1;
36  f->nx = nx = ny;
37  }
38 
39  if (nx == -2) {
40  type = 2;
41  f->nx = nx = 1;
42  }
43 
44  if ((f->kernel = (char *) malloc(nx * ny * sizeof (char))) == NULL) {
45  printf("-E- %s %d: Unable to allocate workspace of %d bytes for filter kernel\n",
46  __FILE__, __LINE__, nx * ny);
47  exit(1);
48  }
49 
50  for (j = 0; j < ny; j++)
51  for (i = 0; i < nx; i++)
52  f->kernel[j * nx + i] = 0;
53 
54  if (type == 0) {
55 
56  /* square */
57 
58  for (j = 0; j < ny; j++)
59  for (i = 0; i < nx; i++)
60  f->kernel[j * nx + i] = 1;
61  if (f->minfill <= 0)
62  f->minfill = nx * ny / 2 + 1;
63  else
64  f->minfill = MIN(MAX(f->minfill, 1), nx * ny / 2 + 1);
65 
66  } else if (type == 2) {
67 
68  /* along-track */
69 
70  for (j = 0; j < ny; j++)
71  for (i = 0; i < nx; i++)
72  f->kernel[j * nx + i] = 1;
73  if (f->minfill <= 0)
74  f->minfill = nx * ny / 2 + 1;
75 
76  } else {
77 
78  /* diamond */
79 
80  for (j = 0; j < ny; j++) {
81  nnx = nx - 2 * abs(ny / 2 - j);
82  for (i = nx / 2 - nnx / 2; i <= nx / 2 + nnx / 2; i++)
83  f->kernel[j * nx + i] = 1;
84  }
85  if (f->minfill <= 0)
86  f->minfill = nx * ny / 4 + 1;
87  else
88  f->minfill = MIN(MAX(f->minfill, 1), nx * ny / 4 + 1);
89  }
90 }
91 
92 int fctl_set(fctlstr *fctl, int32_t npix, char *fname,
93  int32_t band, int32_t nx, int32_t ny, int32_t minfill, int32_t nbands) {
94  int i = fctl->nfilt;
95  int32_t j;
96  int32_t func = -1;
97 
98  /* filter width can't exceed number of pixels per scan */
99  nx = MIN(nx, npix);
100 
101  /* force filter window to odd number */
102  if (nx >= 0)
103  nx = MAX(1, (nx / 2)*2 + 1);
104  ny = MAX(1, (ny / 2)*2 + 1);
105 
106 
107  if (strstr(fname, "dilate"))
108  func = FDILATE;
109  else if (strstr(fname, "stlight"))
110  func = FSTLIGHT;
111  else if (strstr(fname, "clean"))
112  func = FCLEAN;
113  else if (strstr(fname, "ltrmed"))
114  func = FLTRMED;
115  else if (strstr(fname, "ltrmean"))
116  func = FLTRMEAN;
117  else if (strstr(fname, "ltriqmean"))
118  func = FLTRIQMEAN;
119  else if (strstr(fname, "ltmed"))
120  func = FLTMED;
121  else if (strstr(fname, "ltmean"))
122  func = FLTMEAN;
123  else if (strstr(fname, "btdetavg"))
124  func = FBTDETAVG;
125  else if (strstr(fname, "test"))
126  func = FTEST;
127  else if (strstr(fname, "epsmean"))
128  func = FEPSMEAN;
129  else if (strstr(fname, "ltrreject"))
130  func = FLTRREJECT;
131 
132 
133 
134  switch (func) {
135  case FDILATE:
136  if (band < 1 || band > L1_NFLAGS) {
137  if (want_verbose)
138  printf("-E- %s %d: bogus band number %d for filter %s\n",
139  __FILE__, __LINE__, band, fname);
140  return (0);
141  }
142  if (want_verbose)
143  printf("Setting %d x %d dilation filter on %s mask\n",
144  nx, ny, l2_flag_lname[band - 1]);
145  band--;
146  break;
147  case FSTLIGHT:
148  if (band < 1 || band > L1_NFLAGS) {
149  if (want_verbose)
150  printf("-E- %s %d: bogus band number %d for filter %s\n",
151  __FILE__, __LINE__, band, fname);
152  return (0);
153  }
154  if (want_verbose)
155  printf("Setting %d x %d straylight filter on %s mask\n",
156  nx, ny, l2_flag_lname[band - 1]);
157  band--;
158  break;
159  case FCLEAN:
160  if (band < 1 || band > L1_NFLAGS) {
161  if (want_verbose)
162  printf("-E- %s %d: bogus band number %d for filter %s\n",
163  __FILE__, __LINE__, band, fname);
164  return (0);
165  }
166  if (want_verbose)
167  printf("Setting %d x %d cleaning filter on %s mask\n",
168  nx, ny, l2_flag_lname[band - 1]);
169  band--;
170  break;
171  case FLTMEAN:
172  if (want_verbose)
173  printf("Setting %d x %d averaging filter on Lt(%d)\n",
174  nx, ny, band);
175  if (band < 1 || band > nbands) {
176  if (want_verbose)
177  printf("-E- %s %d: bogus band number %d for filter %s\n",
178  __FILE__, __LINE__, band, fname);
179  return (0);
180  }
181  band--;
182  break;
183  case FLTRMEAN:
184  if (want_verbose)
185  printf("Setting %d x %d averaging filter on Lt(%d)-Lr(%d)\n",
186  nx, ny, band, band);
187  if (band < 1 || band > nbands) {
188  if (want_verbose)
189  printf("-E- %s %d: bogus band number %d for filter %s\n",
190  __FILE__, __LINE__, band, fname);
191  return (0);
192  }
193  band--;
194  break;
195  case FLTRREJECT:
196  if (want_verbose)
197  printf("Setting %d x %d rejection filter on Lt(%d)-Lr(%d)\n",
198  nx, ny, band, band);
199  if (band < 1 || band > nbands) {
200  if (want_verbose)
201  printf("-E- %s %d: bogus band number %d for filter %s\n",
202  __FILE__, __LINE__, band, fname);
203  return (0);
204  }
205  band--;
206  break;
207  case FLTRIQMEAN:
208  if (want_verbose)
209  printf("Setting %d x %d interquartile averaging filter on Lt(%d)-Lr(%d)\n",
210  nx, ny, band, band);
211  if (band < 1 || band > nbands) {
212  if (want_verbose)
213  printf("-E- %s %d: bogus band number %d for filter %s\n",
214  __FILE__, __LINE__, band, fname);
215  return (0);
216  }
217  band--;
218  break;
219  case FLTMED:
220  if (want_verbose)
221  printf("Setting %d x %d median filter on Lt(%d)\n",
222  nx, ny, band);
223  if (band < 1 || band > nbands) {
224  if (want_verbose)
225  printf("-E- %s %d: bogus band number %d for filter %s\n",
226  __FILE__, __LINE__, band, fname);
227  return (0);
228  }
229  band--;
230  break;
231  case FLTRMED:
232  if (want_verbose)
233  printf("Setting %d x %d median filter on Lt(%d)-Lr(%d)\n",
234  nx, ny, band, band);
235  if (band < 1 || band > nbands) {
236  if (want_verbose)
237  printf("-E- %s %d: bogus band number %d for filter %s\n",
238  __FILE__, __LINE__, band, fname);
239  return (0);
240  }
241  band--;
242  break;
243  case FEPSMEAN:
244  if (want_verbose)
245  printf("Setting %d x %d smoothing filter on epsilon(%d)\n",
246  nx, ny, band);
247  if (band < 1 || band + 1 > nbands) {
248  if (want_verbose)
249  printf("-E- %s %d: bogus band number %d for filter %s\n",
250  __FILE__, __LINE__, band, fname);
251  return (0);
252  }
253  band--;
254  break;
255  case FBTDETAVG:
256  if (want_verbose)
257  printf("Setting %d x %d filter on IR band %d for detector %d\n",
258  nx, ny, band, minfill);
259  band--; /* band number (0-based) */
260  minfill--; /* detector number (0-based)*/
261  break;
262  case FTEST:
263  if (want_verbose)
264  printf("Setting %d x %d test filter on %d\n",
265  nx, ny, band);
266  band--;
267  break;
268  default:
269  printf("-E- %s %d: unknown filter function %s\n",
270  __FILE__, __LINE__, fname);
271  return (0);
272  }
273 
274  fctl->nfilt++;
275  ;
276  fctl->npix = MAX(nx, fctl->npix);
277  fctl->nscan = MAX(ny, fctl->nscan);
278 
279  fctl->f[i].func = func;
280  fctl->f[i].band = band;
281  fctl->f[i].nx = nx;
282  fctl->f[i].ny = ny;
283  fctl->f[i].minfill = minfill;
284 
285  get_kernel(&fctl->f[i]);
286 
287  if (want_verbose) {
288  printf("\nFilter Kernel");
289  for (j = 0; j < fctl->f[i].nx * fctl->f[i].ny; j++) {
290  if ((j % fctl->f[i].nx) == 0) printf("\n");
291  printf("%d ", fctl->f[i].kernel[j]);
292  }
293  printf("\n\nMinimum fill set to %d pixels\n", fctl->f[i].minfill);
294  printf("\n\n");
295  }
296 
297  return (1);
298 }
299 
300 int rdfilter(char *file, fctlstr *fctl, int32_t nbands) {
301  FILE *fp;
302  char line[80];
303  char *p1, *p2;
304  char fname[80];
305  int band;
306  int nscan;
307  int32_t npix;
308  int32_t minfill;
309  int i;
310 
311  if (want_verbose)
312  printf("Opening filter file %s\n", file);
313 
314  if ((fp = fopen(file, "r")) == NULL) {
315  printf("The specified filter file (%s) was not found.\n", file);
316  exit(1);
317  }
318 
319  while (fgets(line, 80, fp)) {
320  if (line[0] == '#' || line[0] == '\n')
321  continue;
322 
323  p1 = line;
324  if (!(p2 = strchr(p1, ','))) {
325  printf("-E- %s %d: filter parsing error on %s. "
326  "Expecting comma-separated list.\n",
327  __FILE__, __LINE__, file);
328  exit(1);
329  }
330 
331  memset(fname, '\0', sizeof (fname));
332  strncpy(fname, p1, p2 - p1);
333  for (i = 0; i < (p2 - p1); i++)
334  fname[i] = tolower(fname[i]);
335 
336  p1 = p2 + 1;
337  if (!(p2 = strchr(p1, ',')))
338  continue;
339  band = atoi(p1);
340 
341  p1 = p2 + 1;
342  if (!(p2 = strchr(p1, ',')))
343  continue;
344  npix = atoi(p1);
345 
346  p1 = p2 + 1;
347  if (!(p2 = strchr(p1, ',')))
348  continue;
349  nscan = atoi(p1);
350 
351  p1 = p2 + 1;
352  minfill = atoi(p1);
353 
354  fctl_set(fctl, npix, fname, band, npix, nscan, minfill, nbands);
355 
356  }
357 
358  return (0);
359 }
360 
361 void fdilate(l1qstr *l1que, int32_t nx, int32_t ny, int flag, char kernel[], l1str *l1rec) {
362  int32_t nscan = l1que->nq;
363  int32_t npix = l1rec->npix;
364  int32_t ip, ip1, ip2;
365  int32_t is, is1, is2;
366  int32_t i, j, k;
367  char *pflag[NQMAX];
368  char *pf;
369  uint32_t flagID = pow(2L, flag);
370 
371  /* Get pointer to desired l1rec flag */
372  switch (flagID) {
373  case HILT: pf = l1rec->hilt;
374  break;
375  case LAND: pf = l1rec->land;
376  break;
377  case CLOUD: pf = l1rec->cloud;
378  break;
379  case HIGLINT: pf = l1rec->glint;
380  break;
381  default:
382  fprintf(stderr,
383  "-E- %s line %d: The specified flag, %d (%d), is not supported\n",
384  __FILE__, __LINE__, flag + 1, flagID);
385  fprintf(stderr, "for dilation filtering.\n");
386  exit(FATAL_ERROR);
387  }
388 
389  /* Build array of pointers to equivalent queue flag */
390  for (is = 0; is < nscan; is++)
391  switch (flagID) {
392  case HILT: pflag[is] = l1que->r[is].hilt;
393  break;
394  case LAND: pflag[is] = l1que->r[is].land;
395  break;
396  case CLOUD: pflag[is] = l1que->r[is].cloud;
397  break;
398  case HIGLINT: pflag[is] = l1que->r[is].glint;
399  break;
400  }
401 
402  /* Compute queue scan limits for the ROI */
403  is = nscan / 2;
404  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
405  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
406 
407  /* Loop over each pixel of the central scan */
408 
409  for (ip = 0; ip < l1rec->npix; ip++) {
410 
411  /* compute pixel neighbor limits for the ROI */
412  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
413  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
414 
415  /* If the pixel is already flagged, just mask it and go on. */
416  /* Also mask the surrounding pixels in the queue records. */
417  /* This ensures that the dilation masking will be visible */
418  /* to later smoothing processes. */
419 
420  if (*(pf + ip)) {
421  l1rec->mask[ip] = ON;
422  for (j = is1; j <= is2; j++) for (i = ip1; i <= ip2; i++) {
423  k = (j - is1) * nx + (i - ip1);
424  if (kernel[k])
425  l1que->r[j].mask[i] = ON;
426  }
427  continue;
428  }
429 
430  /* Search the ROI for the desired flag */
431  for (j = is1; j <= is2; j++) for (i = ip1; i <= ip2; i++) {
432  k = (j - is1) * nx + (i - ip1);
433  if (kernel[k] && *(pflag[j] + i)) {
434  *(pf + ip) = ON; /* set flag in l1rec */
435  l1rec->mask[ip] = ON; /* set mask in l1rec */
436  goto next_pixel;
437  }
438  }
439  next_pixel:
440  ;
441  }
442 }
443 
444 void fstlight(l1qstr *l1que, int32_t nx, int32_t ny, int32_t dscan, int flag, char kernel[], l1str *l1rec)
445 /*
446  * W. Robinson, SAIC, 16 Aug 2011 - modify for VIIRS aggregation zones
447  */ {
448  int32_t nscan = l1que->nq;
449  int32_t npix = l1rec->npix;
450  int32_t ip, ip1, ip2;
451  int32_t is, is1, is2;
452  int32_t i, j, k;
453  int32_t ip_scan, filt_dist, ipmod;
454  char *pflag[NQMAX];
455  char *pf;
456  uint32_t flagID = pow(2L, flag);
457 
458  /* Get pointer to desired l1rec flag */
459  switch (flagID) {
460  case HILT: pf = l1rec->hilt;
461  break;
462  case LAND: pf = l1rec->land;
463  break;
464  case CLOUD: pf = l1rec->cloud;
465  break;
466  case HIGLINT: pf = l1rec->glint;
467  break;
468  default:
469  fprintf(stderr,
470  "-E- %s line %d: The specified flag, %d (%d), is not supported\n",
471  __FILE__, __LINE__, flag + 1, flagID);
472  fprintf(stderr, "for dilation filtering.\n");
473  exit(FATAL_ERROR);
474  }
475 
476  /* Build array of pointers to equivalent queue flag */
477  for (is = 0; is < nscan / dscan; is++)
478  switch (flagID) {
479  case HILT: pflag[is] = l1que->r[is].hilt;
480  break;
481  case LAND: pflag[is] = l1que->r[is].land;
482  break;
483  case CLOUD: pflag[is] = l1que->r[is].cloud;
484  break;
485  case HIGLINT: pflag[is] = l1que->r[is].glint;
486  break;
487  }
488 
489  /* Compute queue scan limits for the ROI */
490  is = nscan / 2 / dscan;
491  is1 = MIN(MAX(0, is - ny / dscan / 2), nscan / dscan - 1);
492  is2 = MAX(MIN(nscan / dscan - 1, is + ny / dscan / 2), 0);
493  /* Loop over each pixel of the central scan */
494 
495  for (ip = 0; ip < l1rec->npix; ip++) {
496 
497  /* compute pixel neighbor limits for the ROI */
498  /* for aggregated VIIRS, limits are aggregation zone dependent */
499  if (((l1rec->l1file->sensorID == VIIRSN) || (l1rec->l1file->sensorID == VIIRSJ1)) &&
500  (l1rec->scn_fmt == 0)) {
501  ip_scan = l1rec->pixnum[ip]; /* scan pixel */
502  int spix = ip_scan - ip;
503  filt_dist = -nx / 2;
504  viirs_pxcvt_agdel(ip_scan, filt_dist, &ipmod);
505  ipmod -= spix;
506  ip1 = MIN(MAX(0, ipmod), npix - 1);
507 
508  filt_dist = nx / 2;
509  viirs_pxcvt_agdel(ip_scan, filt_dist, &ipmod);
510  ipmod -= spix;
511  ip2 = MAX(MIN(npix - 1, ipmod), 0);
512  } else {
513  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
514  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
515  }
516 
517  /* If the pixel is already flagged, skip it. */
518  if (*(pf + ip))
519  continue;
520 
521  /* Search the ROI for the desired flag */
522  for (j = is1; j <= is2; j++) for (i = ip1; i <= ip2; i++) {
523  k = (j - is1) * nx + (i - ip1);
524  if (kernel[k] && *(pflag[j] + i)) {
525  l1rec->stlight[ip] = ON; /* set straylight flag in l1rec */
526  goto next_pixel;
527  }
528  }
529  next_pixel:
530  ;
531  }
532 }
533 
534 void fclean(l1qstr *l1que, int32_t nx, int32_t ny, int flag, char kernel[], l1str *l1rec) {
535  int32_t nscan = l1que->nq;
536  int32_t npix = l1rec->npix;
537  int32_t ip, ip1, ip2;
538  int32_t is, is1, is2;
539  int32_t i, j, k;
540  char *pflag[NQMAX];
541  char *pf;
542  uint32_t flagID = pow(2L, flag);
543  int32_t cnt = 0;
544  int32_t tot = 0;
545 
546  /* Get pointer to desired l1rec flag */
547  switch (flagID) {
548  case HILT: pf = l1rec->hilt;
549  break;
550  case LAND: pf = l1rec->land;
551  break;
552  case CLOUD: pf = l1rec->cloud;
553  break;
554  case HIGLINT: pf = l1rec->glint;
555  break;
556  default:
557  fprintf(stderr,
558  "-E- %s line %d: The specified flag, %d (%d), is not supported\n",
559  __FILE__, __LINE__, flag + 1, flagID);
560  fprintf(stderr, "for dilation filtering.\n");
561  exit(FATAL_ERROR);
562  }
563 
564  /* Build array of pointers to equivalent queue flag */
565  for (is = 0; is < nscan; is++)
566  switch (flagID) {
567  case HILT: pflag[is] = l1que->r[is].hilt;
568  break;
569  case LAND: pflag[is] = l1que->r[is].land;
570  break;
571  case CLOUD: pflag[is] = l1que->r[is].cloud;
572  break;
573  case HIGLINT: pflag[is] = l1que->r[is].glint;
574  break;
575  }
576 
577  /* Compute queue scan limits for the ROI */
578  is = nscan / 2;
579  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
580  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
581 
582  /* Loop over each pixel of the central scan */
583 
584  for (ip = 0; ip < l1rec->npix; ip++) {
585 
586  /* compute pixel neighbor limits for the ROI */
587  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
588  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
589 
590  /* If the pixel is already flagged, just mask it and go on. */
591  if (*(pf + ip)) {
592  l1rec->mask[ip] = ON;
593  continue;
594  }
595 
596  /* Search the ROI for the desired flag */
597  for (j = is1; j <= is2; j++) for (i = ip1; i <= ip2; i++) {
598  k = (j - is1) * nx + (i - ip1);
599  if (kernel[k]) {
600  tot++;
601  if (*(pflag[j] + i))
602  cnt++;
603  }
604  }
605 
606  /* If every surrounding pixel is flagged, then flag the center pixel */
607  if (tot - cnt == 1) {
608  *(pf + ip) = ON; /* set flag in l1rec */
609  l1rec->mask[ip] = ON; /* set mask in l1rec */
610  }
611  }
612 }
613 
614 void fLTmean(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t minfill, char kernel[],
615  l1str *l1rec) {
616  int32_t nscan = l1que->nq;
617  int32_t npix = l1rec->npix;
618  int32_t ip, ip1, ip2;
619  int32_t is, is1, is2;
620  int32_t ii, iib;
621  int32_t i, j, k;
622  float x;
623  int32_t cnt;
624 
625  /* define range of scan numbers within queue */
626  is = nscan / 2;
627  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
628  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
629 
630  /* */
631  /* Loop through each pixel in the scan */
632  /* */
633  for (ip = 0; ip < l1rec->npix; ip++) {
634 
635  iib = ip * l1rec->l1file->nbands + ib;
636 
637  /* if the pixel is already masked, skip it */
638  if (l1rec->mask[ip] || l1rec->Lt[iib] == badval)
639  continue;
640 
641  /* define pixel range around the central pixel */
642  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
643  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
644 
645 
646  /* */
647  /* Compute the mean over the pixel and scan range */
648  /* and replace the central pixel value with it. */
649  /* */
650  x = 0.0;
651  cnt = 0;
652  for (i = ip1; i <= ip2; i++) {
653  for (j = is1; j <= is2; j++) {
654  ii = i * l1rec->l1file->nbands + ib;
655  k = (j - is1) * nx + (i - ip1);
656  if (kernel[k] && !l1que->r[j].mask[i] && l1que->r[j].Lt[ii] != badval) {
657  x += l1que->r[j].Lt[ii];
658  cnt++;
659  }
660  }
661  }
662  if (cnt >= minfill)
663  l1rec->Lt[iib] = x / cnt;
664  else {
665  l1rec->filter[ip] = ON;
666  l1rec->mask[ip] = ON;
667  }
668  }
669 }
670 
671 void fLTRmean(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t minfill, char kernel[],
672  l1str *l1rec) {
673  int32_t nscan = l1que->nq;
674  int32_t npix = l1rec->npix;
675  int32_t ip, ip1, ip2;
676  int32_t is, is1, is2;
677  int32_t ii, iib;
678  int32_t i, j, k;
679  float x;
680  int32_t cnt;
681 
682  /* define range of scan numbers within queue */
683  is = nscan / 2;
684  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
685  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
686 
687  /* */
688  /* Loop through each pixel in the scan */
689  /* */
690  for (ip = 0; ip < l1rec->npix; ip++) {
691 
692  iib = ip * l1rec->l1file->nbands + ib;
693 
694  /* if the pixel is already masked, skip it */
695  if (l1rec->mask[ip] || l1rec->Lt[iib] == badval)
696  continue;
697 
698  /* define pixel range around the central pixel */
699  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
700  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
701 
702 
703  /* */
704  /* Compute the mean over the pixel and scan range */
705  /* and replace the central pixel value with it. */
706  /* */
707  x = 0.0;
708  cnt = 0;
709  for (i = ip1; i <= ip2; i++) {
710  for (j = is1; j <= is2; j++) {
711  ii = i * l1rec->l1file->nbands + ib;
712  k = (j - is1) * nx + (i - ip1);
713  if (kernel[k] && !l1que->r[j].mask[i] && l1que->r[j].Lt[ii] != badval) {
714  x += (l1que->r[j].Lt[ii] - l1que->r[j].Lr[ii]);
715  cnt++;
716  }
717  }
718  }
719  if (cnt >= minfill)
720  l1rec->Lt[iib] = x / cnt + l1rec->Lr[iib];
721  else {
722  l1rec->filter[ip] = ON;
723  l1rec->mask[ip] = ON;
724  }
725  }
726 }
727 
728 int compfloat(float *x, float *y) {
729  if (*x < *y)
730  return (-1);
731  else
732  return ( 1);
733 }
734 
735 int compfnode(fnode *x, fnode *y) {
736  if (x->v < y->v)
737  return (-1);
738  else
739  return ( 1);
740 }
741 
742 void fLTRreject(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t minfill, char kernel[],
743  l1str *l1rec) {
744  static fnode *medx = NULL;
745  static fnode *mad = NULL;
746  static int32_t len = 0;
747 
748  int32_t nscan = l1que->nq;
749  int32_t npix = l1rec->npix;
750  int32_t ip, ip1, ip2;
751  int32_t is, is1, is2;
752  int32_t ii, iib;
753  int32_t i, j, k;
754  float x;
755  int32_t cnt;
756  float median;
757  float madev;
758 
759  /* define range of scan numbers within queue */
760  is = nscan / 2;
761  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
762  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
763 
764  /* allocate sufficient workspace for the filter size */
765  if (medx == NULL || len < nx * ny) {
766  len = nx*ny;
767  if (medx != NULL) free(medx);
768  if ((medx = (fnode *) malloc(len * sizeof (fnode))) == NULL) {
769  printf("-E- %s %d: Unable to allocate workspace for median\n",
770  __FILE__, __LINE__);
771  return;
772  }
773  }
774  if (mad == NULL || len < nx * ny) {
775  len = nx*ny;
776  if (mad != NULL) free(mad);
777  if ((mad = (fnode *) malloc(len * sizeof (fnode))) == NULL) {
778  printf("-E- %s %d: Unable to allocate workspace for median absolute deviation\n",
779  __FILE__, __LINE__);
780  return;
781  }
782  }
783 
784  /* */
785  /* Loop through each pixel in the scan */
786  /* */
787  for (ip = 0; ip < l1rec->npix; ip++) {
788 
789  iib = ip * l1rec->l1file->nbands + ib;
790 
791  /* if the pixel is already masked, skip it */
792  if (l1rec->mask[ip] || l1rec->Lt[iib] == badval)
793  continue;
794 
795  /* define pixel range around the central pixel */
796  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
797  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
798 
799 
800  /* */
801  /* Compute the mean over the pixel and scan range */
802  /* and replace the central pixel value with it. */
803  /* */
804  cnt = 0;
805  for (i = ip1; i <= ip2; i++) {
806  for (j = is1; j <= is2; j++) {
807  ii = i * l1rec->l1file->nbands + ib;
808  k = (j - is1) * nx + (i - ip1);
809  if (kernel[k] && !l1que->r[j].mask[i] && l1que->r[j].Lt[ii] != badval) {
810  medx[cnt].v = l1que->r[j].Lt[ii] - l1que->r[j].Lr[ii];
811  medx[cnt].i = i;
812  medx[cnt].j = j;
813  cnt++;
814  }
815  }
816  }
817  if (cnt >= minfill && cnt >= 2) {
818  qsort(medx, cnt, sizeof (fnode),
819  (int (*)(const void *, const void *)) compfnode);
820 
821  median = medx[cnt / 2].v;
822 
823  cnt = 0;
824  for (i = ip1; i <= ip2; i++) {
825  for (j = is1; j <= is2; j++) {
826  ii = i * l1rec->l1file->nbands + ib;
827  k = (j - is1) * nx + (i - ip1);
828  if (kernel[k] && !l1que->r[j].mask[i] && l1que->r[j].Lt[ii] != badval) {
829  x = (l1que->r[j].Lt[ii] - l1que->r[j].Lr[ii]) - median;
830  mad[cnt].v = fabs(x);
831  mad[cnt].i = i;
832  mad[cnt].j = j;
833  cnt++;
834  }
835  }
836  }
837  qsort(mad, cnt, sizeof (fnode),
838  (int (*)(const void *, const void *)) compfnode);
839 
840  madev = 1.4826 * mad[cnt / 2].v; // const. 1.4826 makes MAD a consistent
841  // esitmator - equiv. to stdev
842 
843  if (fabs((l1rec->Lt[iib] - l1rec->Lr[iib] - median) / madev) >= 3) {
844  l1rec->filter[ip] = ON;
845  //l1rec->mask[ip] = ON;
846  }
847  }
848  }
849 }
850 
851 void fEPSmean(l1qstr *l1que, int32_t nx, int32_t ny, int ib1, int32_t minfill, char kernel[],
852  l1str *l1rec) {
853  int32_t nscan = l1que->nq;
854  int32_t npix = l1rec->npix;
855  int32_t ip, ip1, ip2;
856  int32_t is, is1, is2;
857  int32_t ii, iib;
858  int32_t i, j, k;
859  float x1;
860  float x2;
861  int32_t cnt;
862  float r_bar;
863  float La1;
864  float La2;
865  float La;
866  int ib2 = l1rec->l1file->nbands - 1;
867 
868  /* define range of scan numbers within queue */
869  is = nscan / 2;
870  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
871  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
872 
873  /* */
874  /* Loop through each pixel in the scan */
875  /* */
876  for (ip = 0; ip < l1rec->npix; ip++) {
877 
878  iib = ip * l1rec->l1file->nbands;
879 
880  /* if the pixel is already masked, skip it */
881  if (l1rec->mask[ip] ||
882  l1rec->Lt[iib + ib1] == badval ||
883  l1rec->Lt[iib + ib2] == badval)
884  continue;
885 
886  /* define pixel range around the central pixel */
887  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
888  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
889 
890  /* */
891  /* Compute the mean over the pixel and scan range */
892  /* and replace the central pixel value with it. */
893  /* */
894  x1 = 0.0;
895  x2 = 0.0;
896  cnt = 0;
897  La = 0.0;
898  for (i = ip1; i <= ip2; i++) {
899  for (j = is1; j <= is2; j++) {
900  k = (j - is1) * nx + (i - ip1);
901  if (kernel[k] && !l1que->r[j].mask[i]) {
902  La1 = 0.0;
903  ii = i * l1rec->l1file->nbands + ib1;
904  if (l1que->r[j].Lt[ii] != badval) {
905  La1 = (l1que->r[j].Lt[ii]
906  / l1que->r[j].tg_sol[ii]
907  / l1que->r[j].tg_sen[ii]
908  - l1que->r[j].tLf[ii]
909  - l1que->r[j].Lr[ii])
910  / l1que->r[j].t_o2[ii];
911  }
912  La2 = 0.0;
913  ii = i * l1rec->l1file->nbands + ib2;
914  if (l1que->r[j].Lt[ii] != badval) {
915  La2 = (l1que->r[j].Lt[ii]
916  / l1que->r[j].tg_sol[ii]
917  / l1que->r[j].tg_sen[ii]
918  - l1que->r[j].tLf[ii]
919  - l1que->r[j].Lr[ii])
920  / l1que->r[j].t_o2[ii];
921  }
922 
923  /* Note: we may be averaging negative values (dark pixels) */
924  x1 += La1;
925  x2 += La2;
926  cnt++;
927 
928  /* Center-pixel long-wavelength value */
929  if (i == ip && j == is) {
930  La = La2;
931  }
932  }
933  }
934  }
935 
936  if (cnt >= minfill) {
937  ii = ip * l1rec->l1file->nbands + ib1;
938  /* if region is anomalously dark, we'll just leave it alone */
939  /* here and handle it elsewhere */
940  if (x1 > 0 && x2 > 0) {
941  r_bar = x1 / x2; /* local average band ratio */
942  La1 = La*r_bar; /* modify La1 at center_pixel */
943  l1rec->Lt[ii] = (La1 * l1rec->t_o2[ii] + l1rec->Lr[ii] + l1rec->tLf[ii])
944  * l1rec->tg_sol[ii]
945  * l1rec->tg_sen[ii];
946  }
947  } else {
948  l1rec->filter[ip] = ON;
949  l1rec->mask[ip] = ON;
950  }
951  }
952 }
953 
954 void fBTdetavg(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t id, char kernel[],
955  l1str *l1rec) {
956  int32_t ndet = 10; /* need to generalize this later */
957  int32_t nscan = l1que->nq;
958  int32_t npix = l1rec->npix;
959  int32_t detnum = l1rec->detnum;
960  int32_t scanum = l1rec->iscan;
961  int32_t ip, ip1, ip2;
962  int32_t is, is1, is2;
963  int32_t ii, ipb;
964  int32_t i, j, k;
965  float x;
966  int32_t cnt;
967 
968  /* only apply to a specific detector number */
969  if (detnum != id) return;
970 
971  /* define range of scan numbers within queue */
972  is = nscan / 2;
973  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
974  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
975 
976  /* limit to detector boundaries */
977  if (detnum == 0) is1 = is + 1;
978  if (detnum == ndet - 1) is2 = is - 1;
979 
980  printf("%d %d %d %d\n", ib, id, is1, is2);
981 
982  /* */
983  /* Loop through each pixel in the scan */
984  /* */
985  for (ip = 0; ip < l1rec->npix; ip++) {
986 
987  ipb = ip * NBANDSIR + ib;
988 
989  /* define pixel range around the central pixel */
990  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
991  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
992 
993  /* */
994  /* Compute the mean over the pixel and scan range */
995  /* and replace the central pixel value with it. */
996  /* */
997  x = 0.0;
998  cnt = 0;
999  for (i = ip1; i <= ip2; i++) {
1000  for (j = is1; j <= is2; j++) {
1001  if (j == scanum) continue;
1002  ii = i * NBANDSIR + ib;
1003  k = (j - is1) * nx + (i - ip1);
1004  if (kernel[k] && !l1que->r[j].mask[i] && l1que->r[j].Ltir[ii] != badval) {
1005  x += l1que->r[j].Bt[ii];
1006  cnt++;
1007  }
1008  }
1009  }
1010  if (cnt > 0) {
1011  l1rec->Bt[ipb] = x / cnt;
1012  l1rec->filter[ip] = ON;
1013  } else {
1014  l1rec->filter[ip] = ON;
1015  l1rec->mask[ip] = ON;
1016  }
1017  }
1018 }
1019 
1020 void fLTmed(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t minfill, char kernel[],
1021  l1str *l1rec) {
1022  static float *x = NULL;
1023  static int32_t len = 0;
1024 
1025  int32_t nscan = l1que->nq;
1026  int32_t npix = l1rec->npix;
1027  int32_t ip, ip1, ip2;
1028  int32_t is, is1, is2;
1029  int32_t ii, iib;
1030  int32_t i, j, k;
1031  int32_t cnt;
1032 
1033  /* allocate sufficient workspace for the filter size */
1034  if (x == NULL || len < nx * ny) {
1035  len = nx*ny;
1036  if (x != NULL) free(x);
1037  if ((x = (float *) malloc(len * sizeof (float))) == NULL) {
1038  printf("-E- %s %d: Unable to allocate workspace for median\n",
1039  __FILE__, __LINE__);
1040  return;
1041  }
1042  }
1043 
1044  /* define range of scan numbers within queue */
1045  is = nscan / 2;
1046  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
1047  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
1048 
1049  /* */
1050  /* Loop through each pixel in the scan */
1051  /* */
1052  for (ip = 0; ip < l1rec->npix; ip++) {
1053 
1054  iib = ip * l1rec->l1file->nbands + ib;
1055 
1056  /* if the pixel is already masked, skip it */
1057  if (l1rec->mask[ip] || l1rec->Lt[iib] == badval)
1058  continue;
1059 
1060  /* define pixel range around the central pixel */
1061  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
1062  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
1063 
1064 
1065  /* */
1066  /* Compute the mean over the pixel and scan range */
1067  /* and replace the central pixel value with it. */
1068  /* */
1069  cnt = 0;
1070  for (i = ip1; i <= ip2; i++) {
1071  for (j = is1; j <= is2; j++) {
1072  ii = i * l1rec->l1file->nbands + ib;
1073  k = (j - is1) * nx + (i - ip1);
1074  if (kernel[k] && !l1que->r[j].mask[i] && l1que->r[j].Lt[ii] != badval) {
1075  x[cnt] = l1que->r[j].Lt[ii];
1076  cnt++;
1077  }
1078  }
1079  }
1080  if (cnt >= minfill) {
1081  qsort(x, cnt, sizeof (float),
1082  (int (*)(const void *, const void *)) compfloat);
1083  l1rec->Lt[iib] = x[cnt / 2];
1084  } else {
1085  l1rec->filter[ip] = ON;
1086  l1rec->mask[ip] = ON;
1087  }
1088  }
1089 }
1090 
1091 void fLTRmed(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t minfill, char kernel[],
1092  l1str *l1rec) {
1093  static fnode *x = NULL;
1094  static int32_t len = 0;
1095 
1096  int32_t nscan = l1que->nq;
1097  int32_t npix = l1rec->npix;
1098  int32_t ip, ip1, ip2;
1099  int32_t is, is1, is2;
1100  int32_t ii, iib;
1101  int32_t i, j, k;
1102  int32_t cnt;
1103 
1104  /* allocate sufficient workspace for the filter size */
1105  if (x == NULL || len < nx * ny) {
1106  len = nx*ny;
1107  if (x != NULL) free(x);
1108  if ((x = (fnode *) malloc(len * sizeof (fnode))) == NULL) {
1109  printf("-E- %s %d: Unable to allocate workspace for median\n",
1110  __FILE__, __LINE__);
1111  return;
1112  }
1113  }
1114 
1115  /* define range of scan numbers within queue */
1116  is = nscan / 2;
1117  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
1118  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
1119 
1120  /* */
1121  /* Loop through each pixel in the scan */
1122  /* */
1123  for (ip = 0; ip < l1rec->npix; ip++) {
1124 
1125  iib = ip * l1rec->l1file->nbands + ib;
1126 
1127  /* if the pixel is already masked, skip it */
1128  if (l1rec->mask[ip] || l1rec->Lt[iib] == badval)
1129  continue;
1130 
1131  /* define pixel range around the central pixel */
1132  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
1133  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
1134 
1135  /* */
1136  /* Accumulate non-masked values over the pixel and scan */
1137  /* */
1138  cnt = 0;
1139  for (i = ip1; i <= ip2; i++) {
1140  for (j = is1; j <= is2; j++) {
1141  ii = i * l1rec->l1file->nbands + ib;
1142  k = (j - is1) * nx + (i - ip1);
1143  if (kernel[k] && !l1que->r[j].mask[i] && l1que->r[j].Lt[ii] != badval) {
1144  x[cnt].v = l1que->r[j].Lt[ii] - l1que->r[j].Lr[ii];
1145  x[cnt].i = i;
1146  x[cnt].j = j;
1147  cnt++;
1148  }
1149  }
1150  }
1151  if (cnt >= minfill) {
1152  qsort(x, cnt, sizeof (fnode),
1153  (int (*)(const void *, const void *)) compfnode);
1154 
1155  l1rec->Lt[iib] = x[cnt / 2].v + l1rec->Lr[iib];
1156 
1157  } else {
1158  l1rec->filter[ip] = ON;
1159  l1rec->mask[ip] = ON;
1160  }
1161  }
1162 }
1163 
1164 void fLTRiqmean(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t minfill, char kernel[],
1165  l1str *l1rec) {
1166  static fnode *x = NULL;
1167  static int32_t len = 0;
1168 
1169  int32_t nscan = l1que->nq;
1170  int32_t npix = l1rec->npix;
1171  int32_t ip, ip1, ip2;
1172  int32_t is, is1, is2;
1173  int32_t ii, iib;
1174  int32_t i, j, k;
1175  int32_t cnt, num;
1176  float val;
1177 
1178  /* allocate sufficient workspace for the filter size */
1179  if (x == NULL || len < nx * ny) {
1180  len = nx*ny;
1181  if (x != NULL) free(x);
1182  if ((x = (fnode *) malloc(len * sizeof (fnode))) == NULL) {
1183  printf("-E- %s %d: Unable to allocate workspace for median\n",
1184  __FILE__, __LINE__);
1185  return;
1186  }
1187  }
1188 
1189  /* define range of scan numbers within queue */
1190  is = nscan / 2;
1191  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
1192  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
1193 
1194  /* */
1195  /* Loop through each pixel in the scan */
1196  /* */
1197  for (ip = 0; ip < l1rec->npix; ip++) {
1198 
1199  iib = ip * l1rec->l1file->nbands + ib;
1200 
1201  /* if the pixel is already masked, skip it */
1202  if (l1rec->mask[ip] || l1rec->Lt[iib] == badval)
1203  continue;
1204 
1205  /* define pixel range around the central pixel */
1206  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
1207  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
1208 
1209  /* */
1210  /* Accumulate non-masked values over the pixel and scan */
1211  /* */
1212  cnt = 0;
1213  for (i = ip1; i <= ip2; i++) {
1214  for (j = is1; j <= is2; j++) {
1215  ii = i * l1rec->l1file->nbands + ib;
1216  k = (j - is1) * nx + (i - ip1);
1217  if (kernel[k] && !l1que->r[j].mask[i] && l1que->r[j].Lt[ii] != badval) {
1218  x[cnt].v = l1que->r[j].Lt[ii] - l1que->r[j].Lr[ii];
1219  x[cnt].i = i;
1220  x[cnt].j = j;
1221  cnt++;
1222  }
1223  }
1224  }
1225  if (cnt >= minfill) {
1226  qsort(x, cnt, sizeof (fnode),
1227  (int (*)(const void *, const void *)) compfnode);
1228 
1229  num = 0;
1230  val = 0.0;
1231 
1232  for (i = cnt / 4; i < 3 * cnt / 4; i++) {
1233  num++;
1234  val += x[i].v;
1235  }
1236  val /= num;
1237 
1238  l1rec->Lt[iib] = val + l1rec->Lr[iib];
1239 
1240  } else {
1241  l1rec->filter[ip] = ON;
1242  l1rec->mask[ip] = ON;
1243  }
1244  }
1245 }
1246 
1247 void fEPSiqmean(l1qstr *l1que, int32_t nx, int32_t ny, int ib1, int32_t minfill, char kernel[],
1248  l1str *l1rec) {
1249  static fnode *x1 = NULL;
1250  static fnode *x2 = NULL;
1251  static int32_t len = 0;
1252 
1253  int32_t nscan = l1que->nq;
1254  int32_t npix = l1rec->npix;
1255  int32_t ip, ip1, ip2;
1256  int32_t is, is1, is2;
1257  int32_t ii, iib;
1258  int32_t i, j, k;
1259  int32_t cnt, num;
1260  float r_bar;
1261  float La1;
1262  float La2;
1263  float La;
1264  int ib2 = l1rec->l1file->nbands - 1;
1265 
1266  /* allocate sufficient workspace for the filter size */
1267  if (x1 == NULL || len < nx * ny) {
1268  len = nx*ny;
1269  if (x1 != NULL) free(x1);
1270  if ((x1 = (fnode *) malloc(len * sizeof (fnode))) == NULL) {
1271  printf("-E- %s %d: Unable to allocate workspace for median\n",
1272  __FILE__, __LINE__);
1273  return;
1274  }
1275  }
1276  if (x2 == NULL || len < nx * ny) {
1277  len = nx*ny;
1278  if (x2 != NULL) free(x2);
1279  if ((x2 = (fnode *) malloc(len * sizeof (fnode))) == NULL) {
1280  printf("-E- %s %d: Unable to allocate workspace for median\n",
1281  __FILE__, __LINE__);
1282  return;
1283  }
1284  }
1285 
1286  /* define range of scan numbers within queue */
1287  is = nscan / 2;
1288  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
1289  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
1290 
1291  /* */
1292  /* Loop through each pixel in the scan */
1293  /* */
1294  for (ip = 0; ip < l1rec->npix; ip++) {
1295 
1296  iib = ip * l1rec->l1file->nbands;
1297 
1298  /* if the pixel is already masked, skip it */
1299  if (l1rec->mask[ip] ||
1300  l1rec->Lt[iib + ib1] <= 0.0 ||
1301  l1rec->Lt[iib + ib2] <= 0.0)
1302  continue;
1303 
1304  /* define pixel range around the central pixel */
1305  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
1306  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
1307 
1308 
1309  /* */
1310  /* Accumulate non-masked values over the pixel and scan */
1311  /* */
1312  cnt = 0;
1313  La = 0.0;
1314  for (i = ip1; i <= ip2; i++) {
1315  for (j = is1; j <= is2; j++) {
1316 
1317  k = (j - is1) * nx + (i - ip1);
1318  if (kernel[k] &&
1319  !l1que->r[j].mask[i] &&
1320  l1que->r[j].Lt[i * l1rec->l1file->nbands + ib1] > 0.0 &&
1321  l1que->r[j].Lt[i * l1rec->l1file->nbands + ib2] > 0.0) {
1322 
1323  La1 = 0.0;
1324  ii = i * l1rec->l1file->nbands + ib1;
1325  La1 = (l1que->r[j].Lt[ii]
1326  / l1que->r[j].tg_sol[ii]
1327  / l1que->r[j].tg_sen[ii]
1328  - l1que->r[j].tLf[ii]
1329  - l1que->r[j].Lr[ii])
1330  / l1que->r[j].t_o2[ii];
1331 
1332  La2 = 0.0;
1333  ii = i * l1rec->l1file->nbands + ib2;
1334  La2 = (l1que->r[j].Lt[ii]
1335  / l1que->r[j].tg_sol[ii]
1336  / l1que->r[j].tg_sen[ii]
1337  - l1que->r[j].tLf[ii]
1338  - l1que->r[j].Lr[ii])
1339  / l1que->r[j].t_o2[ii];
1340 
1341  /* Center-pixel long-wavelength value */
1342  if (i == ip && j == is)
1343  La = La2;
1344 
1345  x1[cnt].v = La1;
1346  x1[cnt].i = i;
1347  x1[cnt].j = j;
1348 
1349  x2[cnt].v = La2;
1350  x2[cnt].i = i;
1351  x2[cnt].j = j;
1352  cnt++;
1353 
1354  }
1355  }
1356  }
1357 
1358  if (cnt >= minfill) {
1359 
1360  qsort(x1, cnt, sizeof (fnode),
1361  (int (*)(const void *, const void *)) compfnode);
1362  qsort(x2, cnt, sizeof (fnode),
1363  (int (*)(const void *, const void *)) compfnode);
1364 
1365  num = 0;
1366  La1 = 0.0;
1367  La2 = 0.0;
1368 
1369  for (i = cnt / 4; i <= 3 * cnt / 4; i++) {
1370  num++;
1371  La1 += x1[i].v;
1372  La2 += x2[i].v;
1373  }
1374  La1 /= num;
1375  La2 /= num;
1376 
1377  /* if region is anomalously dark, we'll just leave it alone */
1378  /* here and handle it elsewhere */
1379 
1380  if (La1 > 0.0 && La2 > 0.0 && La1 / La2 < 2.0) {
1381  ii = ip * l1rec->l1file->nbands + ib1;
1382  r_bar = La1 / La2; /* local average band ratio */
1383  La1 = La*r_bar; /* modify La1 at center_pixel */
1384  l1rec->Lt[ii] = (La1 * l1rec->t_o2[ii] + l1rec->Lr[ii] + l1rec->tLf[ii])
1385  * l1rec->tg_sol[ii]
1386  * l1rec->tg_sen[ii];
1387  }
1388 
1389  } else {
1390  l1rec->filter[ip] = ON;
1391  l1rec->mask[ip] = ON;
1392  }
1393  }
1394 }
1395 
1396 void filter(fctlstr *fctl, l1qstr *l1que, l1str *l1rec, int32_t dscan) {
1397  int i;
1398  int32_t func, band, nx, ny, minfill;
1399  char *kernel;
1400 
1401  for (i = 0; i < fctl->nfilt; i++) {
1402 
1403  func = fctl->f[i].func;
1404  band = fctl->f[i].band;
1405  nx = fctl->f[i].nx;
1406  ny = fctl->f[i].ny;
1407  minfill = fctl->f[i].minfill;
1408  kernel = fctl->f[i].kernel;
1409 
1410  switch (func) {
1411  case FDILATE: fdilate(l1que, nx, ny, band, kernel, l1rec);
1412  break;
1413  case FSTLIGHT: fstlight(l1que, nx, ny, dscan, band, kernel, l1rec);
1414  break;
1415  case FCLEAN: fclean(l1que, nx, ny, band, kernel, l1rec);
1416  break;
1417  case FLTMEAN: fLTmean(l1que, nx, ny, band, minfill, kernel, l1rec);
1418  break;
1419  case FBTDETAVG: fBTdetavg(l1que, nx, ny, band, minfill, kernel, l1rec);
1420  break;
1421  case FLTRMEAN: fLTRmean(l1que, nx, ny, band, minfill, kernel, l1rec);
1422  break;
1423  case FLTRIQMEAN: fLTRiqmean(l1que, nx, ny, band, minfill, kernel, l1rec);
1424  break;
1425  case FLTMED: fLTmed(l1que, nx, ny, band, minfill, kernel, l1rec);
1426  break;
1427  case FLTRMED: fLTRmed(l1que, nx, ny, band, minfill, kernel, l1rec);
1428  break;
1429  case FTEST: fLTRiqmean(l1que, nx, ny, band, minfill, kernel, l1rec);
1430  break;
1431  case FEPSMEAN: fEPSiqmean(l1que, nx, ny, band, minfill, kernel, l1rec);
1432  break;
1433  case FLTRREJECT: fLTRreject(l1que, nx, ny, band, minfill, kernel, l1rec);
1434  break;
1435  default:
1436  printf("-E- %s %d: unknown filter function %d\n",
1437  __FILE__, __LINE__, func);
1438  return;
1439  }
1440 
1441  }
1442 
1443  return;
1444 }
1445 
1446 
FTEST
#define FTEST
Definition: filter.h:18
MAX
#define MAX(A, B)
Definition: swl0_utils.h:26
MIN
#define MIN(x, y)
Definition: rice.h:169
kernel
subroutine kernel(is, mu, rm, xpl, psl, bp)
Definition: 6sm1.f:4700
j
int j
Definition: decode_rs.h:73
L
#define L(lambda, T)
Definition: PreprocessP.h:185
fEPSiqmean
void fEPSiqmean(l1qstr *l1que, int32_t nx, int32_t ny, int ib1, int32_t minfill, char kernel[], l1str *l1rec)
Definition: filter.c:1247
FDILATE
#define FDILATE
Definition: filter.h:7
NBANDSIR
#define NBANDSIR
Definition: filehandle.h:23
FLTMED
#define FLTMED
Definition: filter.h:9
rdfilter
int rdfilter(char *file, fctlstr *fctl, int32_t nbands)
Definition: filter.c:300
NULL
#define NULL
Definition: decode_rs.h:63
FBTDETAVG
#define FBTDETAVG
Definition: filter.h:16
l1rec
read l1rec
Definition: HOWTO_Add_a_sensor.txt:72
VIIRSN
#define VIIRSN
Definition: sensorDefs.h:23
fnode_str::j
int32_t j
Definition: filter.c:15
band
PARAM_TYPE_NONE Default value No parameter is buried in the product name name_prefix is case insensitive string compared to the product name PARAM_TYPE_VIS_WAVE The visible wavelength bands from the sensor are buried in the product name The product name is compared by appending and name_suffix ie aph_412_giop where prod_ix will be set to PARAM_TYPE_IR_WAVE same search method as PARAM_TYPE_VIS_WAVE except only wavelength above are looped through but prod_ix is still based ie aph_2_giop for the second band
Definition: HOWTO_Add_a_product.txt:42
fdilate
void fdilate(l1qstr *l1que, int32_t nx, int32_t ny, int flag, char kernel[], l1str *l1rec)
Definition: filter.c:361
viirs_pxcvt_agdel
void viirs_pxcvt_agdel(int, int, int *)
Definition: viirs_pxcvt.c:106
fLTRiqmean
void fLTRiqmean(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t minfill, char kernel[], l1str *l1rec)
Definition: filter.c:1164
fnode_str::v
float v
Definition: filter.c:16
ON
#define ON
Definition: l1.h:43
FLTRIQMEAN
#define FLTRIQMEAN
Definition: filter.h:14
is
this program makes no use of any feature of the SDP Toolkit that could generate such a then geolocation is calculated at that and then aggregated up to Resolved feature request Bug by adding three new int8 SDSs for each high resolution offsets between the high resolution geolocation and a bi linear interpolation extrapolation of the positions This can be used to reconstruct the high resolution geolocation Resolved Bug by delaying cumulation of gflags until after validation of derived products Resolved Bug by setting Latitude and Longitude to the correct fill resolving to support Near Real Time because they may be unnecessary if use of entrained ephemeris and attitude data is turned resolving bug report Corrected to filter out Aqua attitude records with missing status helping resolve bug MOD_PR03 will still correctly write scan and pixel data that does not depend upon the start thereby resolving MODur00108 Changed header guard macro names to avoid reserved name resolving MODur00104 Maneuver list file for Terra satellite was updated to include data from Jecue DuChateu Maneuver list files for both Terra and Aqua were updated to include two maneuvers from recent IOT weekly reports The limits for Z component of angular momentum was and to set the fourth GEO scan quality flag for a scan depending upon whether or not it occurred during one of them Added _FillValue attributes to many and changed the fill value for the sector start times from resolving MODur00072 Writes boundingcoordinate metadata to L1A archived metadata For PERCENT *ECS change to treat rather resolving GSFcd01518 Added a LogReport Changed to create the Average Temperatures vdata even if the number of scans is
Definition: HISTORY.txt:301
NQMAX
#define NQMAX
Definition: l12_parms.h:12
nscan
int32 nscan
Definition: l1_czcs_hdf.c:19
FCLEAN
#define FCLEAN
Definition: filter.h:13
fEPSmean
void fEPSmean(l1qstr *l1que, int32_t nx, int32_t ny, int ib1, int32_t minfill, char kernel[], l1str *l1rec)
Definition: filter.c:851
file
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude resolving resolving GSFcd00179 Corrected handling of fill values for[Sensor|Solar][Zenith|Azimuth] resolving MODxl01751 Changed to validate LUT version against a value retrieved from the resolving MODxl02056 Changed to calculate Solar Diffuser angles without adjustment for estimated post launch changes in the MODIS orientation relative to incidentally resolving defects MODxl01766 Also resolves MODxl01947 Changed to ignore fill values in SCI_ABNORM and SCI_STATE rather than treating them as resolving MODxl01780 Changed to use spacecraft ancillary data to recognise when the mirror encoder data is being set by side A or side B and to change calculations accordingly This removes the need for seperate LUTs for Side A and Side B data it makes the new LUTs incompatible with older versions of the and vice versa Also resolves MODxl01685 A more robust GRing algorithm is being which will create a non default GRing anytime there s even a single geolocated pixel in a granule Removed obsolete messages from seed file
Definition: HISTORY.txt:413
names::if
character(len=1000) if
Definition: names.f90:13
l12_proto.h
FLTMEAN
#define FLTMEAN
Definition: filter.h:8
l1que
l1qstr l1que
Definition: getl1rec.c:7
compfloat
int compfloat(float *x, float *y)
Definition: filter.c:728
FLTRMEAN
#define FLTRMEAN
Definition: filter.h:10
f
double precision function f(R1)
Definition: tmd.lp.f:1454
compfnode
int compfnode(fnode *x, fnode *y)
Definition: filter.c:735
fLTRmed
void fLTRmed(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t minfill, char kernel[], l1str *l1rec)
Definition: filter.c:1091
hico.value_locate.x
list x
Definition: value_locate.py:64
FLTRREJECT
#define FLTRREJECT
Definition: filter.h:17
color_dtdb.y
y
Definition: color_dtdb.py:430
fctl_init
void fctl_init(fctlstr *fctl)
Definition: filter.c:19
fstlight
void fstlight(l1qstr *l1que, int32_t nx, int32_t ny, int32_t dscan, int flag, char kernel[], l1str *l1rec)
Definition: filter.c:444
func
subroutine func(x, conec, n, bconecno, bn, units, u, inno, i, outno, o, Input, Targ, p, sqerr)
Definition: ffnet.f:287
FEPSMEAN
#define FEPSMEAN
Definition: filter.h:12
fBTdetavg
void fBTdetavg(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t id, char kernel[], l1str *l1rec)
Definition: filter.c:954
FATAL_ERROR
#define FATAL_ERROR
Definition: swl0_parms.h:5
FSTLIGHT
#define FSTLIGHT
Definition: filter.h:15
want_verbose
int want_verbose
Definition: genutils_globals.c:3
FLTRMED
#define FLTRMED
Definition: filter.h:11
fix_mac_rpath.line
line
Definition: fix_mac_rpath.py:55
fnode_str::i
int32_t i
Definition: filter.c:14
CLOUD
#define CLOUD
Definition: l2_flags.h:20
filter
void filter(fctlstr *fctl, l1qstr *l1que, l1str *l1rec, int32_t dscan)
Definition: filter.c:1396
LAND
#define LAND
Definition: l2_flags.h:12
HILT
#define HILT
Definition: l2_flags.h:15
fLTmed
void fLTmed(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t minfill, char kernel[], l1str *l1rec)
Definition: filter.c:1020
BAD_FLT
#define BAD_FLT
Definition: jplaeriallib.h:19
nbands
int32_t nbands
Definition: atrem_corl1v3.h:190
spix
int32 spix
Definition: l1_czcs_hdf.c:21
fabs
#define fabs(a)
Definition: misc.h:93
fclean
void fclean(l1qstr *l1que, int32_t nx, int32_t ny, int flag, char kernel[], l1str *l1rec)
Definition: filter.c:534
HIGLINT
#define HIGLINT
Definition: l2_flags.h:14
fLTmean
void fLTmean(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t minfill, char kernel[], l1str *l1rec)
Definition: filter.c:614
strncpy
PARAM_TYPE_NONE Default value No parameter is buried in the product name name_prefix is case insensitive string compared to the product name PARAM_TYPE_VIS_WAVE The visible wavelength bands from the sensor are buried in the product name The product name is compared by appending and name_suffix ie aph_412_giop where prod_ix will be set to PARAM_TYPE_IR_WAVE same search method as PARAM_TYPE_VIS_WAVE except only wavelength above are looped through but prod_ix is still based ie aph_2_giop for the second and prod_ix set to PARAM_TYPE_INT name_prefix is compared with the beginning of the product name If name_suffix is not empty the it must match the end of the product name The characters right after the prefix are read as an integer and prod_ix is set to that number strncpy(l2prod->name_prefix, "myprod", UNITLEN)
fctl_set
int fctl_set(fctlstr *fctl, int32_t npix, char *fname, int32_t band, int32_t nx, int32_t ny, int32_t minfill, int32_t nbands)
Definition: filter.c:92
fnode_str
Definition: filter.c:13
fLTRreject
void fLTRreject(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t minfill, char kernel[], l1str *l1rec)
Definition: filter.c:742
fLTRmean
void fLTRmean(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t minfill, char kernel[], l1str *l1rec)
Definition: filter.c:671
abs
#define abs(a)
Definition: misc.h:90
i
int i
Definition: decode_rs.h:71
val
msiBandIdx val
Definition: l1c_msi.cpp:34
VIIRSJ1
#define VIIRSJ1
Definition: sensorDefs.h:37
viirs_ler_luts::p1
real p1
Definition: viirs_ler_luts.f95:37
L1_NFLAGS
#define L1_NFLAGS
Definition: filehandle.h:21
k
int k
Definition: decode_rs.h:73
npix
int npix
Definition: get_cmp.c:27
get_kernel
void get_kernel(filstr *f)
Definition: filter.c:25


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Mission Quality Monitor

HELP

Level 1 & 2 Browser Help

Level 3 & 4 Browser Help

File Search Help

Search & Download Help

Earthdata Forum

Join Mailing List

How to Cite Data

CONNECT WITH US
COMPLIANCE

Web Privacy Policy

Data & Information Policy

Communications Policy

Accessibility Policy

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster