#include <cstdio>

/* Perform extended weather filtering (per OMB Tech Note 120) on
     SSMI data for sea ice use */
/* Robert Grumbine 4 June 1997 */
 
float gr37(const ssmis *map, const int i, const int j, 
           const int nx, const int ny, const int range);
float gr22(const ssmis *map, const int i, const int j, 
           const int nx, const int ny, const int range);

/* Variables for tracking which filters get used */
/* 16 March 2004 */
extern int efilt_37_n, efilt_37_s;
extern int efilt_22_n, efilt_22_s;

int ssmis_newfilt(ssmis *nmap, ssmis *smap) {
  int i, j;
  float g37[NY_NORTH][NX_NORTH];
  float g22[NY_NORTH][NX_NORTH];
  bool debug;
  unsigned char nconc[NY_NORTH][NX_NORTH], sconc[NY_SOUTH][NX_SOUTH];
  int range;
 
  #ifdef VERBOSE
    debug = true;
  #else
    debug = false;
  #endif
  #ifdef HIRES
    range = 1;
  #else
    range = 0;
  #endif

  ssmis_getfld(nmap, NX_NORTH*NY_NORTH, &nconc[0][0], &g37[0][0], SSMIS_BAR_CONC);
  ssmis_getfld(smap, NX_SOUTH*NY_SOUTH, &sconc[0][0], &g37[0][0], SSMIS_BAR_CONC);

/* Find the northern hemisphere gradient ratio */
  for (j = 0; j < NY_NORTH  ; j++) {
    for (i = 0; i < NX_NORTH  ; i++) {
      g37[j][i] = gr37(nmap, i, j, NX_NORTH, NY_NORTH, range);
      g22[j][i] = gr22(nmap, i, j, NX_NORTH, NY_NORTH, range);
    }
  }
/* Loop over all points.  If, in any case, the 2 pt averaged gradient ratio 
     for 37-19 v is greater than the cut off, then filter out the ice 
     concentration */

  for (j = 1; j < NY_NORTH - 1 ; j++) {
    for (i = 1; i < NX_NORTH - 1 ; i++) {
      if (nconc[j][i] != 0 && nconc[j][i] != BAD_DATA) {
        /* Include this as either second check (low resolution) or
           an appropriately-ranged check (high resolution 
           Robert Grumbine 16 March 2004 */
        if (g37[j][i] > SSMIS_GR37LIM) {
           efilt_37_n += 1;
           nconc[j][i] = WEATHER;
        }
        if (g22[j][i] > SSMIS_GR22LIM) {
           efilt_22_n += 1;
           nconc[j][i] = WEATHER;
        }

        if (nconc[j+1][i] != BAD_DATA) {
          if (g37[j+1][i] + g37[j][i] > 2*SSMIS_GR37LIM ) {
            efilt_37_n += 1;
            nconc[j][i] = WEATHER;
            if (debug) printf("1 resetting %3d %3d \n",i,j);
            continue;
          }
        }
        if (nconc[j-1][i] != BAD_DATA) {
          if (g37[j-1][i] + g37[j][i] > 2*SSMIS_GR37LIM ) {
            efilt_37_n += 1;
            nconc[j][i] = WEATHER;
            if (debug) printf("2 resetting %3d %3d \n",i,j);
            continue;
          }
        }
        if (nconc[j][i-1] != BAD_DATA) {
          if (g37[j][i-1] + g37[j][i] > 2*SSMIS_GR37LIM ) {
            efilt_37_n += 1;
            nconc[j][i] = WEATHER;
            if (debug) printf("3 resetting %3d %3d \n",i,j);
            continue;
          }
        }
        if (nconc[j][i+1] != BAD_DATA) {
          if (g37[j][i+1] + g37[j][i] > 2*SSMIS_GR37LIM ) {
            efilt_37_n += 1;
            nconc[j][i] = WEATHER;
            if (debug) printf("4 resetting %3d %3d \n",i,j);
            continue;
          }
        }
      
      } /* End of filtration testing */

    }
  }

  for (j = 0; j < NY_NORTH; j++) {
    for (i = 0; i < NX_NORTH; i++) {
       nmap[j*NX_NORTH +i].bar_conc = nconc[j][i];
    }
  }

/* Need to put southern filtering in here */
  for (j = 0; j < NY_SOUTH  ; j++) {
    for (i = 0; i < NX_SOUTH  ; i++) {
      g37[j][i] = gr37(smap, i, j, NX_SOUTH, NY_SOUTH, range);
      g22[j][i] = gr22(smap, i, j, NX_SOUTH, NY_SOUTH, range);
    }
  }
/* Note above that we've used the same array for both north and south
   gradients */
  for (j = 1; j < NY_SOUTH - 1; j++) {
    for (i = 1; i < NX_SOUTH - 1; i++) {
      if (sconc[j][i] != 0 && sconc[j][i] != BAD_DATA) {
        /* Include this as either second check (low resolution) or
           an appropriately-ranged check (high resolution 
           Robert Grumbine 16 March 2004 */
        if (g37[j][i] > SSMIS_GR37LIM) {
           efilt_37_s += 1;
           sconc[j][i] = WEATHER;
        }
        if (g22[j][i] > SSMIS_GR22LIM) {
           efilt_22_s += 1;
           sconc[j][i] = WEATHER;
        }

        if (sconc[j+1][i] != BAD_DATA) {
          if (g37[j+1][i] + g37[j][i] > 2*SSMIS_GR37LIM ) {
            efilt_37_s += 1;
            sconc[j][i] = WEATHER;
            if (debug) printf("1 resetting %3d %3d \n",i,j);
            continue;
          }
        }
        if (sconc[j-1][i] != BAD_DATA) {
          if (g37[j-1][i] + g37[j][i] > 2*SSMIS_GR37LIM ) {
            efilt_37_s += 1;
            sconc[j][i] = WEATHER;
            if (debug) printf("2 resetting %3d %3d \n",i,j);
            continue;
          }
        }
        if (sconc[j][i-1] != BAD_DATA) {
          if (g37[j][i-1] + g37[j][i] > 2*SSMIS_GR37LIM ) {
            efilt_37_s += 1;
            sconc[j][i] = WEATHER;
            if (debug) printf("3 resetting %3d %3d \n",i,j);
            continue;
          }
        }
        if (sconc[j][i+1] != BAD_DATA) {
          if (g37[j][i+1] + g37[j][i] > 2*SSMIS_GR37LIM ) {
            efilt_37_s += 1;
            sconc[j][i] = WEATHER;
            if (debug) printf("4 resetting %3d %3d \n",i,j);
            continue;
          }
        }
      
      } /* End of filtration testing */

    }
  }

  for (j = 0; j < NY_SOUTH; j++) {
    for (i = 0; i < NX_SOUTH; i++) {
       smap[j*NX_SOUTH + i].bar_conc = sconc[j][i];
    }
  }

  return 0;
}


float gr37(const ssmis *map, const int i, const int j,
                            const int nx, const int ny, const int range)
{
   int index, ti, tj, count = 0;
   float t19v = 0.0, t37v = 0.0;

   if (range != 0) {
     for (tj = j-range ; tj < j+range ; tj++) {
       for (ti = i - range ; ti < i+range; ti++) {
         index = ti + tj*nx;
         if (index < 0 || index >= nx*ny) continue;
         if (map[index].t19v != 0 && map[index].t37v != 0) {
           count += 1;
           t19v += map[index].t19v;
           t37v += map[index].t37v;
         }
       }
     }

     t37v = t37v / count;
     t19v = t19v / count;
   }
   else {
     index = i + j*nx;
     t19v = map[index].t19v;
     t37v = map[index].t37v;
   }

   if (t19v != 0.0 && t37v != 0.0 ) {
     return ((t37v - t19v) / (t37v + t19v));
   }
   else {return 0.0;}

   return 0.0;
}

float gr22(const ssmis *map, const int i, const int j, 
                            const int nx, const int ny, const int range)
{
   int index, ti, tj, count = 0;
   float t19v = 0.0, t22v = 0.0;

   if (range != 0) {
     for (tj = j-range ; tj < j+range ; tj++) {
       for (ti = i - range ; ti < i+range; ti++) {
         index = ti + tj*nx;
         if (index < 0 || index >= nx*ny) continue; 
         if (map[index].t19v != 0 && map[index].t22v != 0) {
           count += 1;
           t19v += map[index].t19v;
           t22v += map[index].t22v;
         }
       }
     }
  
     t22v = t22v / count;
     t19v = t19v / count;
   }
   else {
     index = i + j*nx;
     t19v = map[index].t19v;
     t22v = map[index].t22v;
   }

   if (t19v != 0.0 && t22v != 0.0 ) {
     return ((t22v - t19v) / (t22v + t19v));
   }
   else {return 0.0;}
   
   return 0.0;
}