/*********************************************************/
/*   program 16: statistical filters                     */
/*********************************************************/

#define  XSIZE  640
#define  YSIZE  480
#define  XMASK  5
#define  YMASK  5

#define  RED    0
#define  GREEN  1
#define  BLUE   2

#define  MIN    1
#define  MED    2
#define  MAX    3


main()
  {
  unsigned char *ir,*ig,*ib;
  unsigned char *or,*og,*ob;
  unsigned char *i,*i1,*i2,*i3;
  unsigned char *o,*o1;
  unsigned char *pg,*p;
  unsigned char min,med,max;
  long  int  x,y,c;
  long  int  xm,ym;
  long  int  filter;

  /* indicate type of filter to perform */
  filter = MAX;

  /* allocate input and output memory buffers */
  ir = (unsigned char *) malloc (XSIZE*YSIZE);
  ig = (unsigned char *) malloc (XSIZE*YSIZE);
  ib = (unsigned char *) malloc (XSIZE*YSIZE);
  or = (unsigned char *) malloc (XSIZE*YSIZE);
  og = (unsigned char *) malloc (XSIZE*YSIZE);
  ob = (unsigned char *) malloc (XSIZE*YSIZE);
  pg = (unsigned char *) malloc (XMASK*YMASK);

  /* read the input image */
  read_image("input",ir,ig,ib);

  /* process the red, green, and blue planes one at a time */
  for (c = RED; c <= BLUE; c++)
    {
    /* initialize pointers for current color plane */
    if (c == RED  ) i = ir;
    if (c == GREEN) i = ig;
    if (c == BLUE ) i = ib;
    if (c == RED  ) o = or;
    if (c == GREEN) o = og;
    if (c == BLUE ) o = ob;

    /* process all lines in the image */
    for (y = (YMASK/2); y < (YSIZE-(YMASK/2)); y++)
      {
      i1 = i + (y*XSIZE) + (XMASK/2);
      o1 = o + (y*XSIZE) + (XMASK/2);

      /* process all pixels in the line */
      for (x = (XMASK/2); x < (XSIZE-(XMASK/2)); x++)
        {
        i2 = i1 - ((YMASK/2)*XSIZE) - (XMASK/2);
        p  = pg;

        /* get the input pixels for this group */
        for (ym = 0; ym < YMASK; ym++)
          {
          i3 = i2 + (ym*XSIZE);
          for (xm = 0; xm < XMASK; xm++)
            {
            *p++ = *i3++;
            }
          }

        /* sort the pixels by intensity value */
        sort(pg,(XMASK*YMASK),&min,&med,&max);

        /* save correct value to perform requested filter */
        switch (filter)
          {
          case MIN:
            *o1++ = min;
            break;
          case MED:
            *o1++ = min;
            break;
          case MAX:
            *o1++ = min;
            break;
          default:
            break;
          }
        i1++;
        }
      }
    }

  /* write the output image (one plane) */
  write_image("stat_max.bw",or,og,ob);

  /* free memory buffers */
  free(ir);
  free(ig);
  free(ib);
  free(or);
  free(og);
  free(ob);
  free(pg);
  }


/**** sort function ****/
sort(pixel_group,elements,min,med,max)
unsigned char *pixel_group;
long int elements;
unsigned char *min,*med,*max;
  {
  unsigned char *p;
  unsigned char t;
  long int e,n;

  /* perform bubble sort until all elements are in order */
  n = 1;
  while (n > 0)
    {
    p = pixel_group;
    n = 0;
    for (e = 1; e < elements; e++)
      {
      if (*p > *(p+1))
        {
        t = *p;
        *p = *(p+1);
        *(p+1) = t;
        n++;
        }
      p++;
      }
    }

  /* save minimum, median, and maximum values */
  *min = pixel_group[0];
  *med = pixel_group[(elements+1)/2-1];
  *max = pixel_group[elements-1];

  return;
}