/***********************************************************************/
#define OP_NAME 	"ppmadapmedian"
#define VERSION 	"1.05"
#define DATE    	"31.01.98" 
#define AUTHOR  	"Stefan Diener"
/**********************************************************************/

#include <stdlib.h>
#include <stdarg.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <sys/types.h>

#include <STIMP/ppm.c>

struct PPM_Info source, desti;
static int m, n;
static int edge=1, breite, mitte;

double mittel(unsigned char *quelle, int z, int s)
{
  int k, l, summe=0;

  /* Summe des Fensters */
  for (k=-edge; k<=edge; k++)   /* Zeilen */
    for (l=-edge; l<=edge; l++)   /* Spalten */
      summe+=quelle[(z+k)*n+s+l];

  /* Summe/Anzahl berechnen */
  return (double)summe/(double)(breite*breite);
}

int median(unsigned char *quelle, int z, int s)
{
  static unsigned char bucket[256];
  int k, l, summe=0;

  /* Buckets leeren */
  for (k=0; k<256; k++) bucket[k]=0;

  /* Buckets fuellen */
  for (k=-edge; k<=edge; k++)   /* Zeilen */
    for (l=-edge; l<=edge; l++)   /* Spalten */
      bucket[quelle[(z+k)*n+s+l]]++;

  /* Mitte suchen */
  for (k=0; k<256; k++)
  {
    summe+=bucket[k];
    if (summe>mitte) return k;
  }

  /* durchgelaufen ??? */
  return 255;
}

void Do_It(void)
{
  char kandidat[27];
  double resultat[27];
  int i, j, k, l, a;
  double r, g, b, h, s, u, v, rr, gg, bb, mini;
  unsigned char matrix[3][3];
  unsigned char *srcR, *srcG, *srcB, *dstR, *dstG, *dstB;

  srcR = source.redData;
  srcG = source.greenData;
  srcB = source.blueData;
  dstR = desti.redData;
  dstG = desti.greenData;
  dstB = desti.blueData;
  desti.maxval=source.maxval;

  for (i=edge; i<m-edge; i++)   /* Zeilen */
  {
    for (j=edge; j<n-edge; j++)   /* Spalten */
    {
      /* Mediane in jeder Komponente bestimmen */
      matrix[0][0]=median(srcR, i, j);
      matrix[1][1]=median(srcG, i, j);
      matrix[2][2]=median(srcB, i, j);

      /* Median-Matrix fuellen */
      for (k=-edge; k<=edge; k++)   /* Zeilen */
        for (l=-edge; l<=edge; l++)   /* Spalten */
        {
          if (matrix[0][0]==srcR[(i+k)*n+j+l])
          {
            matrix[1][0]=srcG[(i+k)*n+j+l];
            matrix[2][0]=srcB[(i+k)*n+j+l];
          }
          if (matrix[1][1]==srcG[(i+k)*n+j+l])
          {
            matrix[0][1]=srcR[(i+k)*n+j+l];
            matrix[2][1]=srcB[(i+k)*n+j+l];
          }
          if (matrix[2][2]==srcB[(i+k)*n+j+l])
          {
            matrix[0][2]=srcR[(i+k)*n+j+l];
            matrix[1][2]=srcG[(i+k)*n+j+l];
          }
        }

      /* Mittelwerte in jedem Kanal bilden */
      r=mittel(srcR, i, j);
      g=mittel(srcG, i, j);
      b=mittel(srcB, i, j);

      /* RGB -> HSI: auf Singularitaeten pruefen */
      if ((2.0*b-r-g==0.0) || ((r-g)*(r-g)+r*g+2.0*b*(b-r-g)<0.0))
      {
        /* Bedingungen nicht erfuellt => Originalpixel verwenden */
        r=(double) srcR[i*n+j];
        g=(double) srcG[i*n+j];
        b=(double) srcB[i*n+j];

        /* Bedingungen nochmal pruefen */
        if ((2.0*b-r-g==0.0) || ((r-g)*(r-g)+r*g+2.0*b*(b-r-g)<0.0))
        {
          /* auch nicht erfuellt => abbrechen */
          r=0.0;  g=0.0;  b=0.0;
        }
      }

      if ((r!=0.0) || (g!=0.0) || (b!=0.0))
      {
        /* RGB -> HSI: Transformation */
        u=2.0*b-r-g;
        v=r-g;
        h=atan2(v,u);
        s=0.4082482905*sqrt(u*u+v*v);

        /* Listen leeren */
        for (k=0; k<27; k++)
        {
          kandidat[k]=0;
          resultat[k]=0.0;
        }

        /* H-Transformation fuer alle 27 Kandidaten berechnen */
        for (k=0; k<3; k++)   /* 1. Zeile */
          for (l=0; l<3; l++)   /* 2. Zeile */
            for (a=0; a<3; a++)   /* 3. Zeile */
            {
              rr=(double) matrix[0][k];
              gg=(double) matrix[1][l];
              bb=(double) matrix[2][a];

              /* Bedingungen ueberpruefen */
              if ((2.0*bb-rr-gg!=0.0) && ((rr-gg)*(rr-gg)+rr*gg+2.0*bb*(bb-rr-gg)>=0.0))
              {
                u=2.0*bb-rr-gg;
                v=rr-gg;
                resultat[a+l*3+k*9]=fabs(atan2(v,u)-h);
              }
              else resultat[a+l*3+k*9]=1000000.0;
            }

        /* Minimum suchen */
        mini=resultat[0];
        for (k=1; k<27; k++)
          if (resultat[k]<mini) mini=resultat[k];

        /* Kandidaten auswaehlen */
        l=0;
        for (k=0; k<27; k++)
          if (resultat[k]==mini)
          {
            kandidat[k]=1;
            l++;
          }

        if (l>1)
        {
          /* S-Transformation fuer restliche Kandidaten berechnen */
          for (k=0; k<3; k++)   /* 1. Zeile */
            for (l=0; l<3; l++)   /* 2. Zeile */
              for (a=0; a<3; a++)   /* 3. Zeile */
                if (kandidat[a+l*3+k*9]==1)
                {
                  rr=(double) matrix[0][k];
                  gg=(double) matrix[1][l];
                  bb=(double) matrix[2][a];

                  u=2.0*bb-rr-gg;
                  v=rr-gg;
                  resultat[a+l*3+k*9]=fabs(0.4082482905*sqrt(u*u+v*v)-s);
                }

          /* Minimum suchen */
          mini=1000000.0;
          for (k=0; k<27; k++)
            if ((kandidat[k]==1) && (resultat[k]<mini)) mini=resultat[k];

          /* Kandidaten auswaehlen */
          l=0;
          for (k=0; k<27; k++)
            if ((kandidat[k]==1) && (resultat[k]==mini)) l++;
           else kandidat[k]=0;

           if (l>1)
           {
            /* I-Transformation fuer restliche Kandidaten berechnen */
            for (k=0; k<3; k++)   /* 1. Zeile */
              for (l=0; l<3; l++)   /* 2. Zeile */
                for (a=0; a<3; a++)   /* 3. Zeile */
                  if (kandidat[a+l*3+k*9]==1)
                  {
                    rr=(double) matrix[0][k];
                    gg=(double) matrix[1][l];
                    bb=(double) matrix[2][a];

                    resultat[a+l*3+k*9]=(rr+gg+bb)/3.0;
                  }

            /* Maximum suchen */
            mini=0.0;   /* Maximum heisst ausnahmsweise mini */
            for (k=0; k<27; k++)
              if ((kandidat[k]==1) && (resultat[k]>mini)) mini=resultat[k];

            /* Kandidaten auswaehlen */
            for (k=0; k<27; k++)
              if ((kandidat[k]==0) || (resultat[k]!=mini)) kandidat[k]=0;
           }
        }

        /* es bleibt nur ein Kandidat uebrig */
        /* den Kandidaten suchen */
        for (k=0; k<27; k++)
          if (kandidat[k]==1)
          {
            l=k;
            k=27;
          }

        /* Daten aus Median-Matrix ins Zielbild kopieren */
        k=l/9;  l=l%9;
        dstR[(i-edge)*(n-2*edge)+j-edge]=matrix[0][k];
        k=l/3;  l=l%3;
        dstG[(i-edge)*(n-2*edge)+j-edge]=matrix[1][k];
        dstB[(i-edge)*(n-2*edge)+j-edge]=matrix[2][l];
      }
      else   /* Singularitaet */
      {
        /* Quellpixel einfach in Zielbild kopieren */
        dstR[(i-edge)*(n-2*edge)+j-edge]=srcR[i*n+j];
        dstG[(i-edge)*(n-2*edge)+j-edge]=srcG[i*n+j];
        dstB[(i-edge)*(n-2*edge)+j-edge]=srcB[i*n+j];
      }
    }
  }
}

int main(int argc, char **argv)
/* Hauptprogramm */
{
  int i;

  /* offizielle Begruessung */
  PrintOpening(argc,argv);

  /* Uebergebene Parameter auswerten */
  for (i=1; i<argc; i++)
  {
    if ((argv[i][0]=='-') && argv[i][1])
    {
      switch (argv[i][1])
      {
        case '3': edge=1;
                      break;

        case '5': edge=2;
                      break;

        case '7': edge=3;
                      break;

        case '9': edge=4;
                      break;

        case 'v': beVerbose=FALSE;
                      break;

        default: PrintMessage("Unknown parameter: %s", argv[i]);
                     Hilfe();
                     exit(-1);
                     break;
      }
    }

    if (argv[i][0]=='+')
    {
      switch (argv[i][1])
      {
        case 'v': beVerbose=TRUE;
                      break;

        default: PrintMessage("Unknown parameter: %s", argv[i]);
                     Hilfe();
                     exit(-1);
                     break;
      }
    }
  }

  /* Mindestzahl der Argumente pruefen */
  if (argc<3)
  {
    PrintMessage("Wrong number of arguments !");
    Hilfe();
    exit(-1);
  }

  /* Anzahl der Dateinamen überprüfen */
  if (FilenameCount(argc, argv)!=2)
  {
    PrintMessage("Wrong number of file names !");
    Hilfe();
    exit(-1);
  }

  if (ReadPPMFile(GetFilename(1,argc,argv),&source)==0)
  {
    m=source.height;
    n=source.width;

    if ((m<=2*edge) || (n<=2*edge)) PrintMessage("The source picture is too small !!!");
    else
    {
      if (CreatePPMArray(m-2*edge,n-2*edge,&desti)==0)
      {
        PrintMessage("Working ...");

        breite=2*edge+1;
        mitte=2*edge*(edge+1);
        Do_It();

        WritePPMFile(GetFilename(2,argc,argv),&desti);
        FreePPMArray(&desti);
      }
    }

    FreePPMArray(&source);
  }

  PrintClosing();
  exit(0);
}