/* pgmminkowsky.c - read a portable graymap and calculate the Minkowski
** Integrals as a function of the threshold.
**
** Copyright (C) 2000 by Luuk van Dijk/Mind over Matter
**
** Based on pgmhist.c,
** Copyright (C) 1989 by Jef Poskanzer.
**
** Permission to use, copy, modify, and distribute this software and its
** documentation for any purpose and without fee is hereby granted, provided
** that the above copyright notice appear in all copies and that both that
** copyright notice and this permission notice appear in supporting
** documentation.  This software is provided "as is" without express or
** implied warranty.
*/

#include "pgm.h"
#include "mallocvar.h"


#define MAX2(a,b) ( ( (a)>(b) ) ? (a) : (b) )
#define MAX4(a,b,c,d) MAX2( MAX2((a),(b)), MAX2((c),(d)) )

int main( int argc, char** argv ){

  FILE *ifp;

  gray maxval;
  int cols, rows, format;

  gray* prevrow;
  gray* thisrow;
  gray* tmprow;

  int* countTile;
  int* countEdgeX;
  int* countEdgeY;
  int* countVertex;

  int i, col, row;

  int maxtiles, maxedgex, maxedgey, maxvertex;
  int area, perimeter, eulerchi;

  double l2inv, linv;

  /*
   * parse arg and initialize
   */

  pgm_init( &argc, argv );

  if ( argc > 2 ) pm_usage( "[pgmfile]" );

  if ( argc == 2 )
    ifp = pm_openr( argv[1] );
  else
    ifp = stdin;

  /*
   * initialize
   */

  pgm_readpgminit( ifp, &cols, &rows, &maxval, &format );

  prevrow = pgm_allocrow( cols );
  thisrow = pgm_allocrow( cols );

  MALLOCARRAY(countTile   , maxval + 1 );
  MALLOCARRAY(countEdgeX  , maxval + 1 );
  MALLOCARRAY(countEdgeY  , maxval + 1 );
  MALLOCARRAY(countVertex , maxval + 1 );

  if (countTile == NULL || countEdgeX == NULL || countEdgeY == NULL ||
      countVertex == NULL)
      pm_error( "out of memory" );

  for ( i = 0; i <= maxval; i++ ) countTile[i]   = 0;
  for ( i = 0; i <= maxval; i++ ) countEdgeX[i]  = 0;
  for ( i = 0; i <= maxval; i++ ) countEdgeY[i]  = 0;
  for ( i = 0; i <= maxval; i++ ) countVertex[i] = 0;




  /* first row */

  pgm_readpgmrow( ifp, thisrow, cols, maxval, format );

  /* tiles */

  for ( col = 0; col < cols; ++col ) ++countTile[thisrow[col]];

  /* y-edges */

  for ( col = 0; col < cols; ++col ) ++countEdgeY[thisrow[col]];

  /* x-edges */

  ++countEdgeX[thisrow[0]];

  for ( col = 0; col < cols-1; ++col )
    ++countEdgeX[ MAX2(thisrow[col], thisrow[col+1]) ];

  ++countEdgeX[thisrow[cols-1]];

  /* shortcut: for the first row, countVertex == countEdgeX */

  ++countVertex[thisrow[0]];

  for ( col = 0; col < cols-1; ++col )
    ++countVertex[ MAX2(thisrow[col], thisrow[col+1]) ];

  ++countVertex[thisrow[cols-1]];



  for ( row = 1; row < rows; ++row ){

    tmprow = prevrow;
    prevrow = thisrow;
    thisrow = tmprow;

    pgm_readpgmrow( ifp, thisrow, cols, maxval, format );

    /* tiles */

    for ( col = 0; col < cols; ++col ) ++countTile[thisrow[col]];

    /* y-edges */

    for ( col = 0; col < cols; ++col )
      ++countEdgeY[ MAX2(thisrow[col], prevrow[col]) ];
    /* x-edges */

    ++countEdgeX[thisrow[0]];

    for ( col = 0; col < cols-1; ++col )
      ++countEdgeX[ MAX2(thisrow[col], thisrow[col+1]) ];

    ++countEdgeX[thisrow[cols-1]];

    /* vertices */

    ++countVertex[ MAX2(thisrow[0],prevrow[0]) ];

    for ( col = 0; col < cols-1; ++col )
      ++countVertex[
        MAX4(thisrow[col], thisrow[col+1], prevrow[col], prevrow[col+1])
      ];

    ++countVertex[ MAX2(thisrow[cols-1],prevrow[cols-1]) ];

  } /* for row */

  /* now thisrow contains the top row*/

  /* tiles and x-edges have been counted, now upper
     y-edges and top vertices remain */

  /* y-edges */

  for ( col = 0; col < cols; ++col ) ++countEdgeY[ thisrow[col] ];

  /* vertices */

  ++countVertex[thisrow[0]];

  for ( col = 0; col < cols-1; ++col )
    ++countVertex[ MAX2(thisrow[col],thisrow[col+1]) ];

  ++countVertex[ thisrow[cols-1] ];


  /* cleanup */

  maxtiles =  rows    * cols;
  maxedgex =  rows    * (cols+1);
  maxedgey = (rows+1) *  cols;
  maxvertex= (rows+1) * (cols+1);

  l2inv = 1.0/maxtiles;
  linv  = 0.5/(rows+cols);

  /* And print it. */
  printf( "#threshold\t tiles\tx-edges\ty-edges\tvertices\n" );
  printf( "#---------\t -----\t-------\t-------\t--------\n" );
  for ( i = 0; i <= maxval; i++ ){

    if( !(countTile[i] || countEdgeX[i] || countEdgeY[i] || countVertex[i] ) )
      continue; /* skip empty slots */

    area      = maxtiles;
    perimeter = 2*maxedgex + 2*maxedgey - 4*maxtiles;
    eulerchi  = maxtiles - maxedgex - maxedgey + maxvertex;

    printf( "%f\t%6d\t%7d\t%7d\t%8d\t%g\t%g\t%6d\n", (float) i/(1.0*maxval),
        maxtiles, maxedgex, maxedgey, maxvertex,
        area*l2inv, perimeter*linv, eulerchi
        );


    maxtiles -= countTile[i];
    maxedgex -= countEdgeX[i];
    maxedgey -= countEdgeY[i];
    maxvertex-= countVertex[i];

    /*  i, countTile[i], countEdgeX[i], countEdgeY[i], countVertex[i] */

  }

  /* these should be zero: */
  printf( "#  check:\t%6d\t%7d\t%7d\t%8d\n",
          maxtiles, maxedgex, maxedgey, maxvertex );

  pm_close( ifp );

  exit( 0 );

} /*main*/