/************************************************************************/
#define OP_NAME      "ppmcolmedian"
#define VERSION      "1.02"
#define DATE         "31.01.98"
#define AUTHOR       "Stefan Diener"
/************************************************************************/

#include <stdio.h>
#include <stdarg.h>
#include <stdlib.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;

#define quad(x) (x)*(x)

void Do_It(void)
{
  int i, j, k, l, posx, poss, index, abstand1, abstand2, x, y;
  double r, g, b, h, s, u, v, su, mini, Re[81], Im[81];
  unsigned char *srcR,*srcG, *srcB, *dstR, *dstG, *dstB;

  /* Zeiger auf Bilddaten holen */
  srcR=source.redData;
  srcG=source.greenData;
  srcB=source.blueData;
  dstR=desti.redData;
  dstG=desti.greenData;
  dstB=desti.blueData;

  /* gleiches maxval */
  desti.maxval=source.maxval;
  
  /* Schleife ueber alle Punkte des Zielbildes */
  /* => Randproblem unterdruecken */
  for (i=edge; i<(n-edge); i++)
    for (l=edge; l<(m-edge); l++)
    {
      index=l*n+i;   /* Index im Quellbild, Filtermittelpunkt */

      /* Filtermaske abarbeiten */
      posx=0;
      for (j=-edge; j<=edge; j++)
        for (k=-edge; k<=edge; k++)
        {
          /* Punkt aus Originalbild kopieren */
          poss=index+j*n+k;
          r=(double) srcR[poss];
          g=(double) srcG[poss];
          b=(double) srcB[poss];

          /* 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))
          {
            /* RGB -> HSI: Transformation */
            u=2.0*b-r-g;
            v=r-g;
            h=atan2(v,u);
            s=0.4082482905*sqrt(u*u+v*v);

            Re[posx]=s*cos(h);
            Im[posx]=s*sin(h);
          }
          else
          {
            /* grosse Werte => Aussortieren bei Minimumsuche */
            Re[posx]=1000.1;
            Im[posx]=1000.1;
          }
          posx++;
        }

      /* minimale Summe der Quadrate finden */
      mini=6000000.0;
      posx=0;
      for (k=0; k<=4*edge*(edge+1); k++)
      {
        /* keine Singularitaetspunkte betrachten */
        if (Re[k]<1000.0)
        {
          su=0.0;
          for (j=0; j<=4*edge*(edge+1); j++)
          {
            /* keine Singularitaetspunkte betrachten */
            if (Re[j]<1000.0) su+=quad(Re[k]-Re[j])+quad(Im[k]-Im[j]);
          }
          if (su<=mini)
          {
            if (su==mini)
            {
              x=posx/(2*edge+1)-edge;
              y=posx%(2*edge+1)-edge;
              poss=index+x*n+y;
              abstand2=quad(srcR[poss]-srcR[index])+quad(srcG[poss]-srcG[index])+quad(srcB[poss]-srcB[index]);
              if (abstand2<abstand1)
              {
                abstand1=abstand2;
                posx=k;
              }
            }
            else
            {
              mini=su;
              posx=k;
              x=posx/(2*edge+1)-edge;
              y=posx%(2*edge+1)-edge;
              poss=index+x*n+y;
              abstand1=quad(srcR[poss]-srcR[index])+quad(srcG[poss]-srcG[index])+quad(srcB[poss]-srcB[index]);
            }
          }
        }
      }

      if (posx==0)
      {
        /* nur Singularitaetspunkte vorhanden => Originalpixel verwenden */
        j=0;
        k=0;
      }
      else
      {
        /* ausgewaehlten Pixel verwenden */
        j=posx/(2*edge+1)-edge;
        k=posx%(2*edge+1)-edge;
      }

      /* Pixel ins Zielbild kopieren */
      posx=(l-edge)*(n-2*edge)+i-edge;
      poss=index+j*n+k;
      dstR[posx]=srcR[poss];
      dstG[posx]=srcG[poss];
      dstB[posx]=srcB[poss];
    }
}

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 ...");
        Do_It();

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

    FreePPMArray(&source);
  }

  PrintClosing();
  exit(0);
}