A sphere with three ocean waves in differing shades of blue. Next to the sphere, there is the bolded text "Ocean Color". Under that, there are two acronyms, separated by a vertical pipe: OB.DAAC (Ocean Biology Distributed Active Archive Center) and OBPG (Ocean Biology Processing Group). Lorem ipsum dolor sit amet
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) ||
500  (l1rec->l1file->sensorID == VIIRSJ1) ||
501  (l1rec->l1file->sensorID == VIIRSJ2)) &&
502  (l1rec->scn_fmt == 0)) {
503  ip_scan = l1rec->pixnum[ip]; /* scan pixel */
504  int spix = ip_scan - ip;
505  filt_dist = -nx / 2;
506  viirs_pxcvt_agdel(ip_scan, filt_dist, &ipmod);
507  ipmod -= spix;
508  ip1 = MIN(MAX(0, ipmod), npix - 1);
509 
510  filt_dist = nx / 2;
511  viirs_pxcvt_agdel(ip_scan, filt_dist, &ipmod);
512  ipmod -= spix;
513  ip2 = MAX(MIN(npix - 1, ipmod), 0);
514  } else {
515  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
516  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
517  }
518 
519  /* If the pixel is already flagged, skip it. */
520  if (*(pf + ip))
521  continue;
522 
523  /* Search the ROI for the desired flag */
524  for (j = is1; j <= is2; j++) for (i = ip1; i <= ip2; i++) {
525  k = (j - is1) * nx + (i - ip1);
526  if (kernel[k] && *(pflag[j] + i)) {
527  l1rec->stlight[ip] = ON; /* set straylight flag in l1rec */
528  goto next_pixel;
529  }
530  }
531  next_pixel:
532  ;
533  }
534 }
535 
536 void fclean(l1qstr *l1que, int32_t nx, int32_t ny, int flag, char kernel[], l1str *l1rec) {
537  int32_t nscan = l1que->nq;
538  int32_t npix = l1rec->npix;
539  int32_t ip, ip1, ip2;
540  int32_t is, is1, is2;
541  int32_t i, j, k;
542  char *pflag[NQMAX];
543  char *pf;
544  uint32_t flagID = pow(2L, flag);
545  int32_t cnt = 0;
546  int32_t tot = 0;
547 
548  /* Get pointer to desired l1rec flag */
549  switch (flagID) {
550  case HILT: pf = l1rec->hilt;
551  break;
552  case LAND: pf = l1rec->land;
553  break;
554  case CLOUD: pf = l1rec->cloud;
555  break;
556  case HIGLINT: pf = l1rec->glint;
557  break;
558  default:
559  fprintf(stderr,
560  "-E- %s line %d: The specified flag, %d (%d), is not supported\n",
561  __FILE__, __LINE__, flag + 1, flagID);
562  fprintf(stderr, "for dilation filtering.\n");
563  exit(FATAL_ERROR);
564  }
565 
566  /* Build array of pointers to equivalent queue flag */
567  for (is = 0; is < nscan; is++)
568  switch (flagID) {
569  case HILT: pflag[is] = l1que->r[is].hilt;
570  break;
571  case LAND: pflag[is] = l1que->r[is].land;
572  break;
573  case CLOUD: pflag[is] = l1que->r[is].cloud;
574  break;
575  case HIGLINT: pflag[is] = l1que->r[is].glint;
576  break;
577  }
578 
579  /* Compute queue scan limits for the ROI */
580  is = nscan / 2;
581  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
582  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
583 
584  /* Loop over each pixel of the central scan */
585 
586  for (ip = 0; ip < l1rec->npix; ip++) {
587 
588  /* compute pixel neighbor limits for the ROI */
589  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
590  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
591 
592  /* If the pixel is already flagged, just mask it and go on. */
593  if (*(pf + ip)) {
594  l1rec->mask[ip] = ON;
595  continue;
596  }
597 
598  /* Search the ROI for the desired flag */
599  for (j = is1; j <= is2; j++) for (i = ip1; i <= ip2; i++) {
600  k = (j - is1) * nx + (i - ip1);
601  if (kernel[k]) {
602  tot++;
603  if (*(pflag[j] + i))
604  cnt++;
605  }
606  }
607 
608  /* If every surrounding pixel is flagged, then flag the center pixel */
609  if (tot - cnt == 1) {
610  *(pf + ip) = ON; /* set flag in l1rec */
611  l1rec->mask[ip] = ON; /* set mask in l1rec */
612  }
613  }
614 }
615 
616 void fLTmean(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t minfill, char kernel[],
617  l1str *l1rec) {
618  int32_t nscan = l1que->nq;
619  int32_t npix = l1rec->npix;
620  int32_t ip, ip1, ip2;
621  int32_t is, is1, is2;
622  int32_t ii, iib;
623  int32_t i, j, k;
624  float x;
625  int32_t cnt;
626 
627  /* define range of scan numbers within queue */
628  is = nscan / 2;
629  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
630  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
631 
632  /* */
633  /* Loop through each pixel in the scan */
634  /* */
635  for (ip = 0; ip < l1rec->npix; ip++) {
636 
637  iib = ip * l1rec->l1file->nbands + ib;
638 
639  /* if the pixel is already masked, skip it */
640  if (l1rec->mask[ip] || l1rec->Lt[iib] == badval)
641  continue;
642 
643  /* define pixel range around the central pixel */
644  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
645  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
646 
647 
648  /* */
649  /* Compute the mean over the pixel and scan range */
650  /* and replace the central pixel value with it. */
651  /* */
652  x = 0.0;
653  cnt = 0;
654  for (i = ip1; i <= ip2; i++) {
655  for (j = is1; j <= is2; j++) {
656  ii = i * l1rec->l1file->nbands + ib;
657  k = (j - is1) * nx + (i - ip1);
658  if (kernel[k] && !l1que->r[j].mask[i] && l1que->r[j].Lt[ii] != badval) {
659  x += l1que->r[j].Lt[ii];
660  cnt++;
661  }
662  }
663  }
664  if (cnt >= minfill)
665  l1rec->Lt[iib] = x / cnt;
666  else {
667  l1rec->filter[ip] = ON;
668  l1rec->mask[ip] = ON;
669  }
670  }
671 }
672 
673 void fLTRmean(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t minfill, char kernel[],
674  l1str *l1rec) {
675  int32_t nscan = l1que->nq;
676  int32_t npix = l1rec->npix;
677  int32_t ip, ip1, ip2;
678  int32_t is, is1, is2;
679  int32_t ii, iib;
680  int32_t i, j, k;
681  float x;
682  int32_t cnt;
683 
684  /* define range of scan numbers within queue */
685  is = nscan / 2;
686  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
687  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
688 
689  /* */
690  /* Loop through each pixel in the scan */
691  /* */
692  for (ip = 0; ip < l1rec->npix; ip++) {
693 
694  iib = ip * l1rec->l1file->nbands + ib;
695 
696  /* if the pixel is already masked, skip it */
697  if (l1rec->mask[ip] || l1rec->Lt[iib] == badval)
698  continue;
699 
700  /* define pixel range around the central pixel */
701  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
702  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
703 
704 
705  /* */
706  /* Compute the mean over the pixel and scan range */
707  /* and replace the central pixel value with it. */
708  /* */
709  x = 0.0;
710  cnt = 0;
711  for (i = ip1; i <= ip2; i++) {
712  for (j = is1; j <= is2; j++) {
713  ii = i * l1rec->l1file->nbands + ib;
714  k = (j - is1) * nx + (i - ip1);
715  if (kernel[k] && !l1que->r[j].mask[i] && l1que->r[j].Lt[ii] != badval) {
716  x += (l1que->r[j].Lt[ii] - l1que->r[j].Lr[ii]);
717  cnt++;
718  }
719  }
720  }
721  if (cnt >= minfill)
722  l1rec->Lt[iib] = x / cnt + l1rec->Lr[iib];
723  else {
724  l1rec->filter[ip] = ON;
725  l1rec->mask[ip] = ON;
726  }
727  }
728 }
729 
730 int compfloat(float *x, float *y) {
731  if (*x < *y)
732  return (-1);
733  else
734  return ( 1);
735 }
736 
737 int compfnode(fnode *x, fnode *y) {
738  if (x->v < y->v)
739  return (-1);
740  else
741  return ( 1);
742 }
743 
744 void fLTRreject(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t minfill, char kernel[],
745  l1str *l1rec) {
746  static fnode *medx = NULL;
747  static fnode *mad = NULL;
748  static int32_t len = 0;
749 
750  int32_t nscan = l1que->nq;
751  int32_t npix = l1rec->npix;
752  int32_t ip, ip1, ip2;
753  int32_t is, is1, is2;
754  int32_t ii, iib;
755  int32_t i, j, k;
756  float x;
757  int32_t cnt;
758  float median;
759  float madev;
760 
761  /* define range of scan numbers within queue */
762  is = nscan / 2;
763  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
764  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
765 
766  /* allocate sufficient workspace for the filter size */
767  if (medx == NULL || len < nx * ny) {
768  len = nx*ny;
769  if (medx != NULL) free(medx);
770  if ((medx = (fnode *) malloc(len * sizeof (fnode))) == NULL) {
771  printf("-E- %s %d: Unable to allocate workspace for median\n",
772  __FILE__, __LINE__);
773  return;
774  }
775  }
776  if (mad == NULL || len < nx * ny) {
777  len = nx*ny;
778  if (mad != NULL) free(mad);
779  if ((mad = (fnode *) malloc(len * sizeof (fnode))) == NULL) {
780  printf("-E- %s %d: Unable to allocate workspace for median absolute deviation\n",
781  __FILE__, __LINE__);
782  return;
783  }
784  }
785 
786  /* */
787  /* Loop through each pixel in the scan */
788  /* */
789  for (ip = 0; ip < l1rec->npix; ip++) {
790 
791  iib = ip * l1rec->l1file->nbands + ib;
792 
793  /* if the pixel is already masked, skip it */
794  if (l1rec->mask[ip] || l1rec->Lt[iib] == badval)
795  continue;
796 
797  /* define pixel range around the central pixel */
798  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
799  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
800 
801 
802  /* */
803  /* Compute the mean over the pixel and scan range */
804  /* and replace the central pixel value with it. */
805  /* */
806  cnt = 0;
807  for (i = ip1; i <= ip2; i++) {
808  for (j = is1; j <= is2; j++) {
809  ii = i * l1rec->l1file->nbands + ib;
810  k = (j - is1) * nx + (i - ip1);
811  if (kernel[k] && !l1que->r[j].mask[i] && l1que->r[j].Lt[ii] != badval) {
812  medx[cnt].v = l1que->r[j].Lt[ii] - l1que->r[j].Lr[ii];
813  medx[cnt].i = i;
814  medx[cnt].j = j;
815  cnt++;
816  }
817  }
818  }
819  if (cnt >= minfill && cnt >= 2) {
820  qsort(medx, cnt, sizeof (fnode),
821  (int (*)(const void *, const void *)) compfnode);
822 
823  median = medx[cnt / 2].v;
824 
825  cnt = 0;
826  for (i = ip1; i <= ip2; i++) {
827  for (j = is1; j <= is2; j++) {
828  ii = i * l1rec->l1file->nbands + ib;
829  k = (j - is1) * nx + (i - ip1);
830  if (kernel[k] && !l1que->r[j].mask[i] && l1que->r[j].Lt[ii] != badval) {
831  x = (l1que->r[j].Lt[ii] - l1que->r[j].Lr[ii]) - median;
832  mad[cnt].v = fabs(x);
833  mad[cnt].i = i;
834  mad[cnt].j = j;
835  cnt++;
836  }
837  }
838  }
839  qsort(mad, cnt, sizeof (fnode),
840  (int (*)(const void *, const void *)) compfnode);
841 
842  madev = 1.4826 * mad[cnt / 2].v; // const. 1.4826 makes MAD a consistent
843  // esitmator - equiv. to stdev
844 
845  if (fabs((l1rec->Lt[iib] - l1rec->Lr[iib] - median) / madev) >= 3) {
846  l1rec->filter[ip] = ON;
847  //l1rec->mask[ip] = ON;
848  }
849  }
850  }
851 }
852 
853 void fEPSmean(l1qstr *l1que, int32_t nx, int32_t ny, int ib1, int32_t minfill, char kernel[],
854  l1str *l1rec) {
855  int32_t nscan = l1que->nq;
856  int32_t npix = l1rec->npix;
857  int32_t ip, ip1, ip2;
858  int32_t is, is1, is2;
859  int32_t ii, iib;
860  int32_t i, j, k;
861  float x1;
862  float x2;
863  int32_t cnt;
864  float r_bar;
865  float La1;
866  float La2;
867  float La;
868  int ib2 = l1rec->l1file->nbands - 1;
869 
870  /* define range of scan numbers within queue */
871  is = nscan / 2;
872  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
873  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
874 
875  /* */
876  /* Loop through each pixel in the scan */
877  /* */
878  for (ip = 0; ip < l1rec->npix; ip++) {
879 
880  iib = ip * l1rec->l1file->nbands;
881 
882  /* if the pixel is already masked, skip it */
883  if (l1rec->mask[ip] ||
884  l1rec->Lt[iib + ib1] == badval ||
885  l1rec->Lt[iib + ib2] == badval)
886  continue;
887 
888  /* define pixel range around the central pixel */
889  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
890  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
891 
892  /* */
893  /* Compute the mean over the pixel and scan range */
894  /* and replace the central pixel value with it. */
895  /* */
896  x1 = 0.0;
897  x2 = 0.0;
898  cnt = 0;
899  La = 0.0;
900  for (i = ip1; i <= ip2; i++) {
901  for (j = is1; j <= is2; j++) {
902  k = (j - is1) * nx + (i - ip1);
903  if (kernel[k] && !l1que->r[j].mask[i]) {
904  La1 = 0.0;
905  ii = i * l1rec->l1file->nbands + ib1;
906  if (l1que->r[j].Lt[ii] != badval) {
907  La1 = (l1que->r[j].Lt[ii]
908  / l1que->r[j].tg_sol[ii]
909  / l1que->r[j].tg_sen[ii]
910  - l1que->r[j].tLf[ii]
911  - l1que->r[j].Lr[ii])
912  / l1que->r[j].t_o2[ii];
913  }
914  La2 = 0.0;
915  ii = i * l1rec->l1file->nbands + ib2;
916  if (l1que->r[j].Lt[ii] != badval) {
917  La2 = (l1que->r[j].Lt[ii]
918  / l1que->r[j].tg_sol[ii]
919  / l1que->r[j].tg_sen[ii]
920  - l1que->r[j].tLf[ii]
921  - l1que->r[j].Lr[ii])
922  / l1que->r[j].t_o2[ii];
923  }
924 
925  /* Note: we may be averaging negative values (dark pixels) */
926  x1 += La1;
927  x2 += La2;
928  cnt++;
929 
930  /* Center-pixel long-wavelength value */
931  if (i == ip && j == is) {
932  La = La2;
933  }
934  }
935  }
936  }
937 
938  if (cnt >= minfill) {
939  ii = ip * l1rec->l1file->nbands + ib1;
940  /* if region is anomalously dark, we'll just leave it alone */
941  /* here and handle it elsewhere */
942  if (x1 > 0 && x2 > 0) {
943  r_bar = x1 / x2; /* local average band ratio */
944  La1 = La*r_bar; /* modify La1 at center_pixel */
945  l1rec->Lt[ii] = (La1 * l1rec->t_o2[ii] + l1rec->Lr[ii] + l1rec->tLf[ii])
946  * l1rec->tg_sol[ii]
947  * l1rec->tg_sen[ii];
948  }
949  } else {
950  l1rec->filter[ip] = ON;
951  l1rec->mask[ip] = ON;
952  }
953  }
954 }
955 
956 void fBTdetavg(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t id, char kernel[],
957  l1str *l1rec) {
958  int32_t ndet = 10; /* need to generalize this later */
959  int32_t nscan = l1que->nq;
960  int32_t npix = l1rec->npix;
961  int32_t detnum = l1rec->detnum;
962  int32_t scanum = l1rec->iscan;
963  int32_t ip, ip1, ip2;
964  int32_t is, is1, is2;
965  int32_t ii, ipb;
966  int32_t i, j, k;
967  float x;
968  int32_t cnt;
969 
970  /* only apply to a specific detector number */
971  if (detnum != id) return;
972 
973  /* define range of scan numbers within queue */
974  is = nscan / 2;
975  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
976  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
977 
978  /* limit to detector boundaries */
979  if (detnum == 0) is1 = is + 1;
980  if (detnum == ndet - 1) is2 = is - 1;
981 
982  printf("%d %d %d %d\n", ib, id, is1, is2);
983 
984  /* */
985  /* Loop through each pixel in the scan */
986  /* */
987  for (ip = 0; ip < l1rec->npix; ip++) {
988 
989  ipb = ip * NBANDSIR + ib;
990 
991  /* define pixel range around the central pixel */
992  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
993  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
994 
995  /* */
996  /* Compute the mean over the pixel and scan range */
997  /* and replace the central pixel value with it. */
998  /* */
999  x = 0.0;
1000  cnt = 0;
1001  for (i = ip1; i <= ip2; i++) {
1002  for (j = is1; j <= is2; j++) {
1003  if (j == scanum) continue;
1004  ii = i * NBANDSIR + ib;
1005  k = (j - is1) * nx + (i - ip1);
1006  if (kernel[k] && !l1que->r[j].mask[i] && l1que->r[j].Ltir[ii] != badval) {
1007  x += l1que->r[j].Bt[ii];
1008  cnt++;
1009  }
1010  }
1011  }
1012  if (cnt > 0) {
1013  l1rec->Bt[ipb] = x / cnt;
1014  l1rec->filter[ip] = ON;
1015  } else {
1016  l1rec->filter[ip] = ON;
1017  l1rec->mask[ip] = ON;
1018  }
1019  }
1020 }
1021 
1022 void fLTmed(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t minfill, char kernel[],
1023  l1str *l1rec) {
1024  static float *x = NULL;
1025  static int32_t len = 0;
1026 
1027  int32_t nscan = l1que->nq;
1028  int32_t npix = l1rec->npix;
1029  int32_t ip, ip1, ip2;
1030  int32_t is, is1, is2;
1031  int32_t ii, iib;
1032  int32_t i, j, k;
1033  int32_t cnt;
1034 
1035  /* allocate sufficient workspace for the filter size */
1036  if (x == NULL || len < nx * ny) {
1037  len = nx*ny;
1038  if (x != NULL) free(x);
1039  if ((x = (float *) malloc(len * sizeof (float))) == NULL) {
1040  printf("-E- %s %d: Unable to allocate workspace for median\n",
1041  __FILE__, __LINE__);
1042  return;
1043  }
1044  }
1045 
1046  /* define range of scan numbers within queue */
1047  is = nscan / 2;
1048  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
1049  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
1050 
1051  /* */
1052  /* Loop through each pixel in the scan */
1053  /* */
1054  for (ip = 0; ip < l1rec->npix; ip++) {
1055 
1056  iib = ip * l1rec->l1file->nbands + ib;
1057 
1058  /* if the pixel is already masked, skip it */
1059  if (l1rec->mask[ip] || l1rec->Lt[iib] == badval)
1060  continue;
1061 
1062  /* define pixel range around the central pixel */
1063  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
1064  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
1065 
1066 
1067  /* */
1068  /* Compute the mean over the pixel and scan range */
1069  /* and replace the central pixel value with it. */
1070  /* */
1071  cnt = 0;
1072  for (i = ip1; i <= ip2; i++) {
1073  for (j = is1; j <= is2; j++) {
1074  ii = i * l1rec->l1file->nbands + ib;
1075  k = (j - is1) * nx + (i - ip1);
1076  if (kernel[k] && !l1que->r[j].mask[i] && l1que->r[j].Lt[ii] != badval) {
1077  x[cnt] = l1que->r[j].Lt[ii];
1078  cnt++;
1079  }
1080  }
1081  }
1082  if (cnt >= minfill) {
1083  qsort(x, cnt, sizeof (float),
1084  (int (*)(const void *, const void *)) compfloat);
1085  l1rec->Lt[iib] = x[cnt / 2];
1086  } else {
1087  l1rec->filter[ip] = ON;
1088  l1rec->mask[ip] = ON;
1089  }
1090  }
1091 }
1092 
1093 void fLTRmed(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t minfill, char kernel[],
1094  l1str *l1rec) {
1095  static fnode *x = NULL;
1096  static int32_t len = 0;
1097 
1098  int32_t nscan = l1que->nq;
1099  int32_t npix = l1rec->npix;
1100  int32_t ip, ip1, ip2;
1101  int32_t is, is1, is2;
1102  int32_t ii, iib;
1103  int32_t i, j, k;
1104  int32_t cnt;
1105 
1106  /* allocate sufficient workspace for the filter size */
1107  if (x == NULL || len < nx * ny) {
1108  len = nx*ny;
1109  if (x != NULL) free(x);
1110  if ((x = (fnode *) malloc(len * sizeof (fnode))) == NULL) {
1111  printf("-E- %s %d: Unable to allocate workspace for median\n",
1112  __FILE__, __LINE__);
1113  return;
1114  }
1115  }
1116 
1117  /* define range of scan numbers within queue */
1118  is = nscan / 2;
1119  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
1120  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
1121 
1122  /* */
1123  /* Loop through each pixel in the scan */
1124  /* */
1125  for (ip = 0; ip < l1rec->npix; ip++) {
1126 
1127  iib = ip * l1rec->l1file->nbands + ib;
1128 
1129  /* if the pixel is already masked, skip it */
1130  if (l1rec->mask[ip] || l1rec->Lt[iib] == badval)
1131  continue;
1132 
1133  /* define pixel range around the central pixel */
1134  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
1135  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
1136 
1137  /* */
1138  /* Accumulate non-masked values over the pixel and scan */
1139  /* */
1140  cnt = 0;
1141  for (i = ip1; i <= ip2; i++) {
1142  for (j = is1; j <= is2; j++) {
1143  ii = i * l1rec->l1file->nbands + ib;
1144  k = (j - is1) * nx + (i - ip1);
1145  if (kernel[k] && !l1que->r[j].mask[i] && l1que->r[j].Lt[ii] != badval) {
1146  x[cnt].v = l1que->r[j].Lt[ii] - l1que->r[j].Lr[ii];
1147  x[cnt].i = i;
1148  x[cnt].j = j;
1149  cnt++;
1150  }
1151  }
1152  }
1153  if (cnt >= minfill) {
1154  qsort(x, cnt, sizeof (fnode),
1155  (int (*)(const void *, const void *)) compfnode);
1156 
1157  l1rec->Lt[iib] = x[cnt / 2].v + l1rec->Lr[iib];
1158 
1159  } else {
1160  l1rec->filter[ip] = ON;
1161  l1rec->mask[ip] = ON;
1162  }
1163  }
1164 }
1165 
1166 void fLTRiqmean(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t minfill, char kernel[],
1167  l1str *l1rec) {
1168  static fnode *x = NULL;
1169  static int32_t len = 0;
1170 
1171  int32_t nscan = l1que->nq;
1172  int32_t npix = l1rec->npix;
1173  int32_t ip, ip1, ip2;
1174  int32_t is, is1, is2;
1175  int32_t ii, iib;
1176  int32_t i, j, k;
1177  int32_t cnt, num;
1178  float val;
1179 
1180  /* allocate sufficient workspace for the filter size */
1181  if (x == NULL || len < nx * ny) {
1182  len = nx*ny;
1183  if (x != NULL) free(x);
1184  if ((x = (fnode *) malloc(len * sizeof (fnode))) == NULL) {
1185  printf("-E- %s %d: Unable to allocate workspace for median\n",
1186  __FILE__, __LINE__);
1187  return;
1188  }
1189  }
1190 
1191  /* define range of scan numbers within queue */
1192  is = nscan / 2;
1193  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
1194  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
1195 
1196  /* */
1197  /* Loop through each pixel in the scan */
1198  /* */
1199  for (ip = 0; ip < l1rec->npix; ip++) {
1200 
1201  iib = ip * l1rec->l1file->nbands + ib;
1202 
1203  /* if the pixel is already masked, skip it */
1204  if (l1rec->mask[ip] || l1rec->Lt[iib] == badval)
1205  continue;
1206 
1207  /* define pixel range around the central pixel */
1208  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
1209  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
1210 
1211  /* */
1212  /* Accumulate non-masked values over the pixel and scan */
1213  /* */
1214  cnt = 0;
1215  for (i = ip1; i <= ip2; i++) {
1216  for (j = is1; j <= is2; j++) {
1217  ii = i * l1rec->l1file->nbands + ib;
1218  k = (j - is1) * nx + (i - ip1);
1219  if (kernel[k] && !l1que->r[j].mask[i] && l1que->r[j].Lt[ii] != badval) {
1220  x[cnt].v = l1que->r[j].Lt[ii] - l1que->r[j].Lr[ii];
1221  x[cnt].i = i;
1222  x[cnt].j = j;
1223  cnt++;
1224  }
1225  }
1226  }
1227  if (cnt >= minfill) {
1228  qsort(x, cnt, sizeof (fnode),
1229  (int (*)(const void *, const void *)) compfnode);
1230 
1231  num = 0;
1232  val = 0.0;
1233 
1234  for (i = cnt / 4; i < 3 * cnt / 4; i++) {
1235  num++;
1236  val += x[i].v;
1237  }
1238  val /= num;
1239 
1240  l1rec->Lt[iib] = val + l1rec->Lr[iib];
1241 
1242  } else {
1243  l1rec->filter[ip] = ON;
1244  l1rec->mask[ip] = ON;
1245  }
1246  }
1247 }
1248 
1249 void fEPSiqmean(l1qstr *l1que, int32_t nx, int32_t ny, int ib1, int32_t minfill, char kernel[],
1250  l1str *l1rec) {
1251  static fnode *x1 = NULL;
1252  static fnode *x2 = NULL;
1253  static int32_t len = 0;
1254 
1255  int32_t nscan = l1que->nq;
1256  int32_t npix = l1rec->npix;
1257  int32_t ip, ip1, ip2;
1258  int32_t is, is1, is2;
1259  int32_t ii, iib;
1260  int32_t i, j, k;
1261  int32_t cnt, num;
1262  float r_bar;
1263  float La1;
1264  float La2;
1265  float La;
1266  int ib2 = l1rec->l1file->nbands - 1;
1267 
1268  /* allocate sufficient workspace for the filter size */
1269  if (x1 == NULL || len < nx * ny) {
1270  len = nx*ny;
1271  if (x1 != NULL) free(x1);
1272  if ((x1 = (fnode *) malloc(len * sizeof (fnode))) == NULL) {
1273  printf("-E- %s %d: Unable to allocate workspace for median\n",
1274  __FILE__, __LINE__);
1275  return;
1276  }
1277  }
1278  if (x2 == NULL || len < nx * ny) {
1279  len = nx*ny;
1280  if (x2 != NULL) free(x2);
1281  if ((x2 = (fnode *) malloc(len * sizeof (fnode))) == NULL) {
1282  printf("-E- %s %d: Unable to allocate workspace for median\n",
1283  __FILE__, __LINE__);
1284  return;
1285  }
1286  }
1287 
1288  /* define range of scan numbers within queue */
1289  is = nscan / 2;
1290  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
1291  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
1292 
1293  /* */
1294  /* Loop through each pixel in the scan */
1295  /* */
1296  for (ip = 0; ip < l1rec->npix; ip++) {
1297 
1298  iib = ip * l1rec->l1file->nbands;
1299 
1300  /* if the pixel is already masked, skip it */
1301  if (l1rec->mask[ip] ||
1302  l1rec->Lt[iib + ib1] <= 0.0 ||
1303  l1rec->Lt[iib + ib2] <= 0.0)
1304  continue;
1305 
1306  /* define pixel range around the central pixel */
1307  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
1308  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
1309 
1310 
1311  /* */
1312  /* Accumulate non-masked values over the pixel and scan */
1313  /* */
1314  cnt = 0;
1315  La = 0.0;
1316  for (i = ip1; i <= ip2; i++) {
1317  for (j = is1; j <= is2; j++) {
1318 
1319  k = (j - is1) * nx + (i - ip1);
1320  if (kernel[k] &&
1321  !l1que->r[j].mask[i] &&
1322  l1que->r[j].Lt[i * l1rec->l1file->nbands + ib1] > 0.0 &&
1323  l1que->r[j].Lt[i * l1rec->l1file->nbands + ib2] > 0.0) {
1324 
1325  La1 = 0.0;
1326  ii = i * l1rec->l1file->nbands + ib1;
1327  La1 = (l1que->r[j].Lt[ii]
1328  / l1que->r[j].tg_sol[ii]
1329  / l1que->r[j].tg_sen[ii]
1330  - l1que->r[j].tLf[ii]
1331  - l1que->r[j].Lr[ii])
1332  / l1que->r[j].t_o2[ii];
1333 
1334  La2 = 0.0;
1335  ii = i * l1rec->l1file->nbands + ib2;
1336  La2 = (l1que->r[j].Lt[ii]
1337  / l1que->r[j].tg_sol[ii]
1338  / l1que->r[j].tg_sen[ii]
1339  - l1que->r[j].tLf[ii]
1340  - l1que->r[j].Lr[ii])
1341  / l1que->r[j].t_o2[ii];
1342 
1343  /* Center-pixel long-wavelength value */
1344  if (i == ip && j == is)
1345  La = La2;
1346 
1347  x1[cnt].v = La1;
1348  x1[cnt].i = i;
1349  x1[cnt].j = j;
1350 
1351  x2[cnt].v = La2;
1352  x2[cnt].i = i;
1353  x2[cnt].j = j;
1354  cnt++;
1355 
1356  }
1357  }
1358  }
1359 
1360  if (cnt >= minfill) {
1361 
1362  qsort(x1, cnt, sizeof (fnode),
1363  (int (*)(const void *, const void *)) compfnode);
1364  qsort(x2, cnt, sizeof (fnode),
1365  (int (*)(const void *, const void *)) compfnode);
1366 
1367  num = 0;
1368  La1 = 0.0;
1369  La2 = 0.0;
1370 
1371  for (i = cnt / 4; i <= 3 * cnt / 4; i++) {
1372  num++;
1373  La1 += x1[i].v;
1374  La2 += x2[i].v;
1375  }
1376  La1 /= num;
1377  La2 /= num;
1378 
1379  /* if region is anomalously dark, we'll just leave it alone */
1380  /* here and handle it elsewhere */
1381 
1382  if (La1 > 0.0 && La2 > 0.0 && La1 / La2 < 2.0) {
1383  ii = ip * l1rec->l1file->nbands + ib1;
1384  r_bar = La1 / La2; /* local average band ratio */
1385  La1 = La*r_bar; /* modify La1 at center_pixel */
1386  l1rec->Lt[ii] = (La1 * l1rec->t_o2[ii] + l1rec->Lr[ii] + l1rec->tLf[ii])
1387  * l1rec->tg_sol[ii]
1388  * l1rec->tg_sen[ii];
1389  }
1390 
1391  } else {
1392  l1rec->filter[ip] = ON;
1393  l1rec->mask[ip] = ON;
1394  }
1395  }
1396 }
1397 
1398 void filter(fctlstr *fctl, l1qstr *l1que, l1str *l1rec, int32_t dscan) {
1399  int i;
1400  int32_t func, band, nx, ny, minfill;
1401  char *kernel;
1402 
1403  for (i = 0; i < fctl->nfilt; i++) {
1404 
1405  func = fctl->f[i].func;
1406  band = fctl->f[i].band;
1407  nx = fctl->f[i].nx;
1408  ny = fctl->f[i].ny;
1409  minfill = fctl->f[i].minfill;
1410  kernel = fctl->f[i].kernel;
1411 
1412  switch (func) {
1413  case FDILATE: fdilate(l1que, nx, ny, band, kernel, l1rec);
1414  break;
1415  case FSTLIGHT: fstlight(l1que, nx, ny, dscan, band, kernel, l1rec);
1416  break;
1417  case FCLEAN: fclean(l1que, nx, ny, band, kernel, l1rec);
1418  break;
1419  case FLTMEAN: fLTmean(l1que, nx, ny, band, minfill, kernel, l1rec);
1420  break;
1421  case FBTDETAVG: fBTdetavg(l1que, nx, ny, band, minfill, kernel, l1rec);
1422  break;
1423  case FLTRMEAN: fLTRmean(l1que, nx, ny, band, minfill, kernel, l1rec);
1424  break;
1425  case FLTRIQMEAN: fLTRiqmean(l1que, nx, ny, band, minfill, kernel, l1rec);
1426  break;
1427  case FLTMED: fLTmed(l1que, nx, ny, band, minfill, kernel, l1rec);
1428  break;
1429  case FLTRMED: fLTRmed(l1que, nx, ny, band, minfill, kernel, l1rec);
1430  break;
1431  case FTEST: fLTRiqmean(l1que, nx, ny, band, minfill, kernel, l1rec);
1432  break;
1433  case FEPSMEAN: fEPSiqmean(l1que, nx, ny, band, minfill, kernel, l1rec);
1434  break;
1435  case FLTRREJECT: fLTRreject(l1que, nx, ny, band, minfill, kernel, l1rec);
1436  break;
1437  default:
1438  printf("-E- %s %d: unknown filter function %d\n",
1439  __FILE__, __LINE__, func);
1440  return;
1441  }
1442 
1443  }
1444 
1445  return;
1446 }
1447 
1448 
FTEST
#define FTEST
Definition: filter.h:18
MAX
#define MAX(A, B)
Definition: swl0_utils.h:25
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:1249
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:1166
fnode_str::v
float v
Definition: filter.c:16
ON
#define ON
Definition: l1.h:44
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:853
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:730
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:737
fLTRmed
void fLTRmed(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t minfill, char kernel[], l1str *l1rec)
Definition: filter.c:1093
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:956
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:56
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:1398
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:1022
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:536
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:616
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
VIIRSJ2
#define VIIRSJ2
Definition: sensorDefs.h:44
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:744
fLTRmean
void fLTRmean(l1qstr *l1que, int32_t nx, int32_t ny, int ib, int32_t minfill, char kernel[], l1str *l1rec)
Definition: filter.c:673
abs
#define abs(a)
Definition: misc.h:90
i
int i
Definition: decode_rs.h:71
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:28
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

How to Use Earthdata Search

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

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster