/* 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*/