/* AUTORIGHTS
This file is part of Ferret Toolkit.

This file is extracted from the Fourier image processing toolkit, and modified to work with the Ferret Toolkit.

Following is the copyright information of the original file:

 Copyright (C) 2007

 Developed by
 M. Emre Celebi
 Department of Computer Science
 Louisiana State University in Shreveport
 Shreveport, LA 71115, USA
 WWW:
 http://sourceforge.net/projects/fourier-ipal
 http://www.lsus.edu/faculty/~ecelebi/fourier.htm

 Project administrator: M. Emre Celebi


Ferret Toolkit is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software Foundation,
Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

*/

#include <cass.h>
/** 
 * @file srm.c
 * Routines for statistical region merging segmentation
 */

/* these macros are copied from src/image.h from original fourier package */

#define ERROR_RET( err_msg, err_no ) \
         do \
          { \
           fprintf ( stderr, "Error in %s: %s\n", func_name, ( err_msg ) ); \
           fflush ( stderr ); \
/*           if ( get_err_mode ( ) ) abort ( ); */\
           return ( err_no ) ; \
          } while ( 0 )

#define SET_FUNC_NAME( f_name ) static const char func_name[] = f_name

/**
 * @brief Determines whether or not a given pointer is NULL
 */

#define IS_NULL( x ) ( ( x ) == NULL )

/**
 * @brief Determines whether or not a given number is in the range of [0,255]
 * @warning This is not a safe macro
 */

#define IS_BYTE( x ) ( ( 0 <= ( x ) ) && ( ( x ) < 256 ) )

/**
 * @brief Determines whether or not a given floating-point number is close to 0.0
 */

#ifndef DBL_EPSILON
#define DBL_EPSILON 1e-9L
#endif

#define IS_ZERO( x ) ( fabs ( ( x ) ) < DBL_EPSILON )

/**
 * @brief Determines whether or not a given floating-point number is positive
 */

#define IS_POS( x ) ( ( x ) > DBL_EPSILON )

/**
 * @brief Determines whether or not a given floating-point number is negative
 */

#define IS_NEG( x ) ( ( x ) < -DBL_EPSILON )

#define IS_ODD( x ) ( ( x ) % 2 )
#define IS_EVEN( x ) ( !IS_ODD ( ( x ) ) )
#define IS_POS_ODD( x ) ( ( ( x ) > 0 ) && IS_ODD ( x ) )

/**
 * @brief Determines whether or not a given number is in the range of [0,1]
 * @warning This is not a safe macro
 */
#define IS_IN_0_1( x ) ( ! ( IS_NEG ( ( x )  ) || IS_POS ( ( x ) - 1.0 ) ) )


//#include <image.h>
//
#define NUM_GRAY	256

typedef struct
{
 int reg1, reg2, delta;
} RegionPair;

static int find_set ( const int *parent, int i );
static int union_set ( const int i, const int j, int *parent, int *rank );
static RegionPair *bucket_sort ( const RegionPair * pair, const int num_elems );

/** @cond INTERNAL_FUNCTION */

static int
find_set ( const int *parent, int i )
{
 while ( parent[i] != i )
  {
   i = parent[i];
  }

 return i;
}

/** @endcond INTERNAL_FUNCTION */

/** @cond INTERNAL_FUNCTION */

static int
union_set ( const int i, const int j, int *parent, int *rank )
{
 if ( rank[i] > rank[j] )
  {
   parent[j] = i;
   return i;
  }

 parent[i] = j;
 if ( rank[i] == rank[j] )
  {
   rank[j]++;
  }

 return j;
}

/** @endcond INTERNAL_FUNCTION */

/** @cond INTERNAL_FUNCTION */

static RegionPair *
bucket_sort ( const RegionPair * pair, const int num_elems )
{
 int ih;
 int *histo;
 int *cum_histo;
 RegionPair *sorted;

 sorted = ( RegionPair * ) calloc ( num_elems, sizeof ( RegionPair ) );
 if ( IS_NULL ( sorted ) )
  {
   return NULL;
  }

 histo = ( int * ) calloc ( NUM_GRAY, sizeof ( int ) );
 cum_histo = ( int * ) malloc ( NUM_GRAY * sizeof ( int ) );

 /* Calculate histogram */
 for ( ih = 0; ih < num_elems; ih++ )
  {
   histo[pair[ih].delta]++;
  }

 /* Calculate cumulative histogram */
 cum_histo[0] = 0;
 for ( ih = 1; ih < NUM_GRAY; ih++ )
  {
   cum_histo[ih] = cum_histo[ih - 1] + histo[ih - 1];
  }

 /* Perform bucket sort */
 for ( ih = 0; ih < num_elems; ih++ )
  {
   sorted[cum_histo[pair[ih].delta]++] = pair[ih];
  }

 free ( histo );
 free ( cum_histo );

 return sorted;
}

/** @endcond INTERNAL_FUNCTION */

/** @cond INTERNAL_MACRO */
#define MIN_2( x, y ) ( ( ( x ) > ( y ) ) ? ( y ) : ( x ) )

/** @endcond INTERNAL_MACRO */

/** @cond INTERNAL_MACRO */
#define MAX_2( x, y ) ( ( ( x ) > ( y ) ) ? ( x ) : ( y ) )

/** @endcond INTERNAL_MACRO */

/** @cond INTERNAL_MACRO */
#define MAX_3( x, y, z ) MAX_2 ( ( x ), MAX_2( ( y ), ( z ) ) )

/** @endcond INTERNAL_MACRO */

/** 
 * @brief Statistical Region Merging segmentation algorithm 
 * 
 * @param[in] in_img Image pointer { rgb }
 * @param[in] Q_value Q parameter of the algorithm (higher values result in 
 *                    more regions) { positive }
 * @param[in] size_factor Regions smaller than SIZE_FACTOR * (# pixels in the
 *                        image) pixels will be merged with others { [0,1] }
 *
 * @return Segmented image or NULL
 * @reco For a typical image try Q_VALUE = 32 and SIZE_FACTOR = 0.001
 *
 * @ref Nock R. and Nielsen F. (2004) "Statistical Region Merging" IEEE 
 *      Trans. on Pattern Analysis and Machine Intelligence, 26(11): 1452-1458
 *      http://www.univ-ag.fr/~rnock/Articles/Drafts/tpami04-nn.pdf
 *
 * @author M. Emre Celebi
 * @date 06.22.2007
 */


#ifndef byte
typedef uchar byte;
#endif

typedef uchar pixel_t[3];

extern double Q_value;
extern double size_factor;

int
image_segment ( void **output, int *num_ccs, uchar *in_data_1d, int num_cols, int num_rows)
{
 SET_FUNC_NAME ( "srm" );
 byte *out_data;
 pixel_t **in_data_3d;
 int red, green, blue;
 int ir, ic, ik;
 int num_pixels;
 int num_pixels_t3;
 int num_rows_m1, num_cols_m1;
 int cnt, idx;
 int num_edges;
 int reg1, reg2;
 int min_reg_size;		/* minimum region size */
 int root, total_size;
 int *size;			/* holds the region sizes */
 int *parent;			/* holds the parents */
 int *rank;			/* holds the ranks of the trees */
 double log_delta;
 double threshold;		/* threshold for merge operation */
 double thresh_factor;		/* constant used in the computation of the threshold */
 double *red_mean;		/* holds the mean red values of the regions */
 double *green_mean;		/* holds the mean green values of the regions */
 double *blue_mean;		/* holds the mean blue values of the regions */
 RegionPair *pair;		/* holds the region pairs */
 RegionPair *sorted;		/* holds the sorted region pairs */

 int num_region;

 if ( IS_NEG ( Q_value ) )
  {
   warn( "Q value ( %f ) must be positive !", Q_value );
   return -1;
  }

 if ( !IS_IN_0_1 ( size_factor ) )
  {
   warn( "Size factor ( %f ) must be in [0,1] range !", size_factor );
   return -1;
  }

 num_pixels = num_rows * num_cols;
 num_pixels_t3 = 3 * num_pixels;

 num_rows_m1 = num_rows - 1;
 num_cols_m1 = num_cols - 1;

// !!!!!!!!!!!!!!!!!!!
 in_data_3d = type_array2matrix(pixel_t, num_rows, num_cols, (pixel_t *)in_data_1d);

 log_delta = 2.0 * log ( 6.0 * num_pixels );
 thresh_factor = ( NUM_GRAY * NUM_GRAY ) / ( 2.0 * Q_value );
 min_reg_size = size_factor * num_pixels;

 red_mean = ( double * ) malloc ( num_pixels * sizeof ( double ) );
 green_mean = ( double * ) malloc ( num_pixels * sizeof ( double ) );
 blue_mean = ( double * ) malloc ( num_pixels * sizeof ( double ) );
 size = ( int * ) malloc ( num_pixels * sizeof ( int ) );
 parent = ( int * ) malloc ( num_pixels * sizeof ( int ) );
 rank = ( int * ) calloc ( num_pixels, sizeof ( int ) );

 if ( IS_NULL ( red_mean ) || IS_NULL ( green_mean ) || IS_NULL ( blue_mean )
      || IS_NULL ( size ) || IS_NULL ( parent ) || IS_NULL ( rank ) )
  {
   ERROR_RET ( "Insufficient memory !", -1 );
  }

 /* In the beginning, each pixel forms one region */
 cnt = 0;
 for ( ik = 0; ik < num_pixels_t3; ik += 3 )
  {
   red_mean[cnt] = in_data_1d[ik];
   green_mean[cnt] = in_data_1d[ik + 1];
   blue_mean[cnt] = in_data_1d[ik + 2];
   size[cnt] = 1;
   parent[cnt] = cnt;
   cnt++;
  }

 idx = 0;
 num_edges = 2 * num_cols_m1 * num_rows_m1 + num_rows_m1 + num_cols_m1;
 pair = ( RegionPair * ) calloc ( num_edges, sizeof ( RegionPair ) );
 if ( IS_NULL ( pair ) )
  {
   ERROR_RET ( "Insufficient memory !", -1 );
  }

 cnt = 0;
 for ( ir = 0; ir < num_rows_m1; ir++ )
  {
   for ( ic = 0; ic < num_cols_m1; ic++ )
    {
     red = in_data_3d[ir][ic][0];
     green = in_data_3d[ir][ic][1];
     blue = in_data_3d[ir][ic][2];

     /* East neighbor */
     pair[idx].reg1 = cnt;
     pair[idx].reg2 = cnt + 1;
     pair[idx].delta = MAX_3 ( abs ( in_data_3d[ir][ic + 1][0] - red ),
			       abs ( in_data_3d[ir][ic + 1][1] - green ),
			       abs ( in_data_3d[ir][ic + 1][2] - blue ) );

     /* South neighbor */
     idx++;
     pair[idx].reg1 = cnt;
     pair[idx].reg2 = cnt + num_cols;
     pair[idx].delta = MAX_3 ( abs ( in_data_3d[ir + 1][ic][0] - red ),
			       abs ( in_data_3d[ir + 1][ic][1] - green ),
			       abs ( in_data_3d[ir + 1][ic][2] - blue ) );

     idx++;
     cnt++;
    }
   cnt++;
  }

 /* Last column of each row */
 cnt = num_cols_m1;
 for ( ir = 0; ir < num_rows_m1; ir++ )
  {
   pair[idx].reg1 = cnt;
   cnt += num_cols;
   pair[idx].reg2 = cnt;
   pair[idx].delta =
    MAX_3 ( abs
	    ( in_data_3d[ir + 1][num_cols_m1][0] -
	      in_data_3d[ir][num_cols_m1][0] ),
	    abs ( in_data_3d[ir + 1][num_cols_m1][1] -
		  in_data_3d[ir][num_cols_m1][1] ),
	    abs ( in_data_3d[ir + 1][num_cols_m1][2] -
		  in_data_3d[ir][num_cols_m1][2] ) );
   idx++;
  }

 /* Last row of each column */
 cnt = num_rows_m1 * num_cols;
 for ( ic = 0; ic < num_cols_m1; ic++ )
  {
   pair[idx].reg1 = cnt;
   cnt++;
   pair[idx].reg2 = cnt;
   pair[idx].delta =
    MAX_3 ( abs
	    ( in_data_3d[num_rows_m1][ic + 1][0] -
	      in_data_3d[num_rows_m1][ic][0] ),
	    abs ( in_data_3d[num_rows_m1][ic + 1][1] -
		  in_data_3d[num_rows_m1][ic][1] ),
	    abs ( in_data_3d[num_rows_m1][ic + 1][2] -
		  in_data_3d[num_rows_m1][ic][2] ) );
   idx++;
  }

 sorted = bucket_sort ( pair, num_edges );

 /* Merge similar regions */
 for ( ik = 0; ik < num_edges; ik++ )
  {
   reg1 = find_set ( parent, sorted[ik].reg1 );
   reg2 = find_set ( parent, sorted[ik].reg2 );

   if ( reg1 != reg2 )
    {
     threshold = sqrt ( thresh_factor
			*
			( ( ( MIN_2 ( NUM_GRAY, size[reg1] ) *
			      log ( 1.0 + size[reg1] ) +
			      log_delta ) / size[reg1] ) + ( ( MIN_2 ( NUM_GRAY,
								       size
								       [reg2] )
							       * log ( 1.0 +
								       size
								       [reg2] )
							       +
							       log_delta ) /
							     size[reg2] ) ) );

     if ( ( fabs ( red_mean[reg1] - red_mean[reg2] ) < threshold )
	  && ( fabs ( green_mean[reg1] - green_mean[reg2] ) < threshold )
	  && ( fabs ( blue_mean[reg1] - blue_mean[reg2] ) < threshold ) )
      {
       root = union_set ( reg1, reg2, parent, rank );
       total_size = size[reg1] + size[reg2];

       red_mean[root] =
	( size[reg1] * red_mean[reg1] +
	  size[reg2] * red_mean[reg2] ) / total_size;
       green_mean[root] =
	( size[reg1] * green_mean[reg1] +
	  size[reg2] * green_mean[reg2] ) / total_size;
       blue_mean[root] =
	( size[reg1] * blue_mean[reg1] +
	  size[reg2] * blue_mean[reg2] ) / total_size;
       size[root] = total_size;
      }
    }
  }

 /* Merge small regions */
 cnt = 0;
 for ( ir = 0; ir < num_rows; ir++ )
  {
   cnt++;
   for ( ic = 1; ic < num_cols; ic++ )
    {
     reg1 = find_set ( parent, cnt );
     reg2 = find_set ( parent, cnt - 1 );

     if ( ( reg1 != reg2 ) &&
	  ( ( size[reg2] < min_reg_size ) || ( size[reg1] < min_reg_size ) ) )
      {
       root = union_set ( reg1, reg2, parent, rank );
       total_size = size[reg1] + size[reg2];

       red_mean[root] =
	( size[reg1] * red_mean[reg1] +
	  size[reg2] * red_mean[reg2] ) / total_size;
       green_mean[root] =
	( size[reg1] * green_mean[reg1] +
	  size[reg2] * green_mean[reg2] ) / total_size;
       blue_mean[root] =
	( size[reg1] * blue_mean[reg1] +
	  size[reg2] * blue_mean[reg2] ) / total_size;
       size[root] = total_size;
      }
     cnt++;
    }
  }

 /* Allocate output image */
 out_data = type_calloc(uchar, num_pixels);


 /* reuse the rank array to map the regions to numbers 0~ num_region-1 */
 num_region = 0;
 for (ik = 0; ik < num_pixels; ik++) rank[ik] = -1;

 /* Assign each output pixel the mean color of its region */
 cnt = 0;
 for ( ik = 0; ik < num_pixels_t3; ik += 3 )
  {
   idx = find_set ( parent, cnt );
   if (rank[idx] < 0)
   {
	   rank[idx] = num_region;
	   num_region++;
   }
   out_data[cnt] = rank[idx];
   cnt++;
  }

 *num_ccs = num_region;

 free ( red_mean );
 free ( green_mean );
 free ( blue_mean );
 free ( size );
 free ( parent );
 free ( rank );
 free ( pair );
 free ( sorted );

 matrix_free_index(in_data_3d);

 *output = out_data;

 return 0;
}

#undef MIN_2
#undef MAX_2
#undef MAX_3

double Q_value = 128;
double size_factor = 0.005;