next up previous contents
Next: Feature Extraction - ad Up: HeLa Data Previous: mb_texture.c

cvip_pgmtexture.c

NOTE: The following code has been modified from its original form in order to correct errors in some calculations.

 /* pgmtxtur.c - calculate textural features on a portable graymap
**
** Author: James Darrell McCauley
**         Texas Agricultural Experiment Station
**         Department of Agricultural Engineering
**         Texas A&M University
**         College Station, Texas 77843-2117 USA
**
** Code written partially taken from pgmtofs.c in the PBMPLUS package
** by Jef Poskanzer.
**
** Algorithms for calculating features (and some explanatory comments) are
** taken from:
**
**   Haralick, R.M., K. Shanmugam, and I. Dinstein. 1973. Textural features
**   for image classification.  IEEE Transactions on Systems, Man, and
**   Cybertinetics, SMC-3(6):610-621.
**
** Copyright (C) 1991 Texas Agricultural Experiment Station, employer for
** hire of James Darrell McCauley
**
** 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.
**
** THE TEXAS AGRICULTURAL EXPERIMENT STATION (TAES) AND THE TEXAS A&M
** UNIVERSITY SYSTEM (TAMUS) MAKE NO EXPRESS OR IMPLIED WARRANTIES
** (INCLUDING BY WAY OF EXAMPLE, MERCHANTABILITY) WITH RESPECT TO ANY
** ITEM, AND SHALL NOT BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL
** OR CONSEQUENTAL DAMAGES ARISING OUT OF THE POSESSION OR USE OF
** ANY SUCH ITEM. LICENSEE AND/OR USER AGREES TO INDEMNIFY AND HOLD
** TAES AND TAMUS HARMLESS FROM ANY CLAIMS ARISING OUT OF THE USE OR
** POSSESSION OF SUCH ITEMS.
** 
** Modification History:
** 24 Jun 91 - J. Michael Carstensen  supplied fix for 
**             correlation function.
**
** Aug. 7 96 - Wenxing Li: huge memory leaks are fixed.

   23 Nov 98 - M. Boland : Compile with the following for use with Matlab
                 under Red Hat Linux 5.1 (i.e. use libc5 instead of glibc)

gcc -c -B /usr/libc5/usr/lib/gcc-lib/ -nostdinc -nostdinc++ \
-I/usr/libc5/usr/include \
-I/usr/libc5/usr/lib/gcc-lib/i386-linux/2.7.2.1/include \
-I/home/boland/Matlab/Mex/Include -ansi cvip_pgmtexture.c

ar r libmb_cvip.a cvip_pgmtexture.o

>> mex -f gccopts.sh -lmb_cvip -L/home/boland/Matlab/Mex  mb_texture.c


   29 Nov 98 - M. Boland : Modified calculations to produce the same values
                 as the kharalick routine in Khoros.  Some feature 
		 calculations below were wrong, others made different 
		 assumptions (the Haralick paper is not always explicit).

*/

#include 
#include "ppgm.h"
#include "CVIPtexture.h"

#define RADIX 2.0
#define EPSILON 0.000000001
#define BL  "Angle                 "
#define F1  "Angular Second Moment "
#define F2  "Contrast              "
#define F3  "Correlation           "
#define F4  "Variance              "
#define F5  "Inverse Diff Moment   "
#define F6  "Sum Average           "
#define F7  "Sum Variance          "
#define F8  "Sum Entropy           "
#define F9  "Entropy               "
#define F10 "Difference Variance   "
#define F11 "Difference Entropy    "
#define F12 "Meas of Correlation-1 "
#define F13 "Meas of Correlation-2 "
#define F14 "Max Correlation Coeff "

#define SIGN(x,y) ((y)<0 ? -fabs(x) : fabs(x))
#define DOT fprintf(stderr,".")
#define SWAP(a,b) {y=(a);(a)=(b);(b)=y;}


 

void results (),  mkbalanced (), reduction (), simplesrt ();
int hessenberg ();
float f1_asm (), f2_contrast (), f3_corr (), f4_var (), f5_idm (),
 f6_savg (), f7_svar (), f8_sentropy (), f9_entropy (), f10_dvar (),
 f11_dentropy (), f12_icorr (), f13_icorr (), f14_maxcorr (), 
  *pgm_vector (), **pgm_matrix ();



TEXTURE * Extract_Texture_Features(int distance, 
				   register gray **grays, 
				   int rows, 
				   int cols, 
				   TEXTURE_FEATURE_MAP *feature_usage)  
{
  FILE *ifp;
  register gray  *gP;
  int tone[PGM_MAXMAXVAL+1], R0, R45, R90, R135, angle, d = 1, x, y;
  int argn, bps, padright, row, col, i;
  int itone, jtone, tones,g_val;
  float **P_matrix0, **P_matrix45, **P_matrix90, **P_matrix135;
  float ASM[4], contrast[4], corr[4], var[4], idm[4], savg[4];
  float sentropy[4], svar[4], entropy[4], dvar[4], dentropy[4];
  float icorr[4], maxcorr[4];
  float *Tp;
  gray nmaxval, maxval;
  char *usage = "[-d ] [pgmfile]";
  TEXTURE *Texture;

    Texture = (TEXTURE *) calloc(1,sizeof(TEXTURE));
    if(!Texture) {
        printf("\nERROR in TEXTURE structure allocate\n");
        exit(1);
    }

  d = distance; 

   /* Determine the number of different gray scales (not maxval) */
  for (row = PGM_MAXMAXVAL; row >= 0; --row)
    tone[row] = -1;
  for (row = rows - 1; row >= 0; --row)
    for (col = 0; col < cols; ++col)
      {
   /*   if (grays[row][col])   If gray value equal 0 don't include */		
        tone[grays[row][col]] = grays[row][col];
      }	
  
 for (row = PGM_MAXMAXVAL, tones = 0; row >= 0; --row)
    if (tone[row] != -1)
      tones++;
 /* fprintf (stderr, "(Image has %d graylevels.)\n", tones); */

  /* Collapse array, taking out all zero values */
  for (row = 0, itone = 0; row <= PGM_MAXMAXVAL; row++)
    if (tone[row] != -1)
      tone[itone++] = tone[row];
  /* Now array contains only the gray levels present (in ascending order) */

  /* Allocate memory for gray-tone spatial dependence matrix */
  P_matrix0 = pgm_matrix (0, tones, 0, tones);
  P_matrix45 = pgm_matrix (0, tones, 0, tones);
  P_matrix90 = pgm_matrix (0, tones, 0, tones);
  P_matrix135 = pgm_matrix (0, tones, 0, tones);
  for (row = 0; row < tones; ++row)
    for (col = 0; col < tones; ++col)
    {
      P_matrix0[row][col] = P_matrix45[row][col] = 0;
      P_matrix90[row][col] = P_matrix135[row][col] = 0;
    }

   R0 = 0;
   R45 = 0;
   R90 = 0;
   R135 = 0;	
  /* Find gray-tone spatial dependence matrix */
 /* fprintf (stderr, "(Computing spatial dependence matrix..."); */
 
  for (row = 0; row < rows; ++row)
    for (col = 0; col < cols; ++col)
      if (grays[row][col])  /* if value anything other than zero */
      for (x = 0, angle = 0; angle <= 135; angle += 45)
      {
	while (tone[x] != grays[row][col])
	  x++;
	/* M. Boland if (angle == 0 && col + d < cols)  */
	/* M. Boland - include neighbor only if != 0 */
	if (angle == 0 && col + d < cols && grays[row][col + d]) 
	{
	  y = 0;
	  while (tone[y] != grays[row][col + d])
	    y++;
  	  P_matrix0[x][y]++;
 	  P_matrix0[y][x]++;
  	  /* R0++;  M. Boland 25 Nov 98 */
	  R0+=2 ;
	}
	/* M. Boland if (angle == 90 && row + d < rows) */
	/* M. Boland - include neighbor only if != 0 */
	else if (angle == 90 && row + d < rows && grays[row + d][col])
	{
	  y = 0;
	  while (tone[y] != grays[row + d][col])
	    y++;
	    P_matrix90[x][y]++;
	    P_matrix90[y][x]++;
   	    /* R90++;  M. Boland 25 Nov 98 */
	    R90+=2 ;
	}
	/* M. Boland if (angle == 45 && row + d < rows && col - d >= 0) */
	/* M. Boland - include neighbor only if != 0 */
	else if (angle == 45 && row + d < rows && col - d >= 0 
		 && grays[row + d][col - d])
	{
	  y = 0;
	  while (tone[y] != grays[row + d][col - d])
	    y++;
  	  P_matrix45[x][y]++;
	  P_matrix45[y][x]++;
	  /* R45++;  M. Boland 25 Nov 98 */
	  R45+=2 ;
	}
	/* M. Boland if (angle == 135 && row + d < rows && col + d < cols) */
	else if (angle == 135 && row + d < rows && col + d < cols 
		 && grays[row + d][col + d])
	{
	  y = 0;
	  while (tone[y] != grays[row + d][col + d])
	    y++;
	  P_matrix135[x][y]++;
	  P_matrix135[y][x]++;
	  /* R135++;  M. Boland 25 Nov 98 */
	  R135+=2 ;
	}

     }	
  /* Gray-tone spatial dependence matrices are complete */

  /* Find normalizing constants */
/* R0 = 2 * rows * (cols - 1);
  R45 = 2 * (rows - 1) * (cols - 1);
  R90 = 2 * (rows - 1) * cols;
  R135 = R45;
*/

  /* Normalize gray-tone spatial dependence matrix */
  for (itone = 0; itone < tones; ++itone)
    for (jtone = 0; jtone < tones; ++jtone)
    {
      P_matrix0[itone][jtone] /= R0;
      P_matrix45[itone][jtone] /= R45;
      P_matrix90[itone][jtone] /= R90;
      P_matrix135[itone][jtone] /= R135;
    }

/*  fprintf (stderr, " done.)\n"); */
/*  fprintf (stderr, "(Computing textural features"); */
/*  fprintf (stdout, "\n"); */
/*  DOT; */
/*  fprintf (stdout,
"%s         0         45         90        135        Avg       Range\n",
	   BL);
*/
  if (feature_usage->ASM)
	{
  	ASM[0] = f1_asm (P_matrix0, tones);
  	ASM[1] = f1_asm (P_matrix45, tones);
  	ASM[2] = f1_asm (P_matrix90, tones);
  	ASM[3] = f1_asm (P_matrix135, tones);
	}
  Tp = &Texture->ASM[0];
  results (Tp, F1, ASM);

  if (feature_usage->contrast)
	{
  	contrast[0] = f2_contrast (P_matrix0, tones);
  	contrast[1] = f2_contrast (P_matrix45, tones);
  	contrast[2] = f2_contrast (P_matrix90, tones);
  	contrast[3] = f2_contrast (P_matrix135, tones);
	}
  Tp = &Texture->contrast[0];
  results (Tp, F2, contrast);

  if (feature_usage->correlation)
	{
 	corr[0] = f3_corr (P_matrix0, tones);
 	corr[1] = f3_corr (P_matrix45, tones);
  	corr[2] = f3_corr (P_matrix90, tones);
  	corr[3] = f3_corr (P_matrix135, tones);
	}
  Tp = &Texture->correlation[0];
  results (Tp, F3, corr);

  if (feature_usage->variance)
	{
  	var[0] = f4_var (P_matrix0, tones);
  	var[1] = f4_var (P_matrix45, tones);
  	var[2] = f4_var (P_matrix90, tones);
  	var[3] = f4_var (P_matrix135, tones);
	}
  Tp = &Texture->variance[0];
  results (Tp, F4, var); 

  if (feature_usage->IDM)
	{
  	idm[0] = f5_idm (P_matrix0, tones);
  	idm[1] = f5_idm (P_matrix45, tones);
  	idm[2] = f5_idm (P_matrix90, tones);
  	idm[3] = f5_idm (P_matrix135, tones);
	}
  Tp = &Texture->IDM[0];
  results (Tp, F5, idm); 

  if (feature_usage->sum_avg)
	{
  	savg[0] = f6_savg (P_matrix0, tones);
  	savg[1] = f6_savg (P_matrix45, tones);
  	savg[2] = f6_savg (P_matrix90, tones);
  	savg[3] = f6_savg (P_matrix135, tones);
	}
  Tp = &Texture->sum_avg[0];
  results (Tp, F6, savg); 

  if (feature_usage->sum_var)
	{
	sentropy[0] = f8_sentropy (P_matrix0, tones);
  	sentropy[1] = f8_sentropy (P_matrix45, tones);
  	sentropy[2] = f8_sentropy (P_matrix90, tones);
  	sentropy[3] = f8_sentropy (P_matrix135, tones);
	}
  if (feature_usage->sum_entropy)
	{  
	svar[0] = f7_svar (P_matrix0, tones, sentropy[0]);
  	svar[1] = f7_svar (P_matrix45, tones, sentropy[1]);
  	svar[2] = f7_svar (P_matrix90, tones, sentropy[2]);
  	svar[3] = f7_svar (P_matrix135, tones, sentropy[3]);
	}
  Tp = &Texture->sum_var[0];
  results (Tp, F7, svar); 
  Tp = &Texture->sum_entropy[0];
  results (Tp, F8, sentropy); 

  if (feature_usage->entropy)
	{
  	entropy[0] = f9_entropy (P_matrix0, tones);
  	entropy[1] = f9_entropy (P_matrix45, tones);
  	entropy[2] = f9_entropy (P_matrix90, tones);
  	entropy[3] = f9_entropy (P_matrix135, tones);
	}
  Tp = &Texture->entropy[0];
  results (Tp, F9, entropy); 

  if (feature_usage->diff_var)
	{  
	dvar[0] = f10_dvar (P_matrix0, tones);
 	dvar[1] = f10_dvar (P_matrix45, tones);
  	dvar[2] = f10_dvar (P_matrix90, tones);
  	dvar[3] = f10_dvar (P_matrix135, tones);
	}
  Tp = &Texture->diff_var[0];
  results (Tp, F10, dvar);

  if (feature_usage->diff_entropy)
	{  
	dentropy[0] = f11_dentropy (P_matrix0, tones);
  	dentropy[1] = f11_dentropy (P_matrix45, tones);
  	dentropy[2] = f11_dentropy (P_matrix90, tones);
  	dentropy[3] = f11_dentropy (P_matrix135, tones);
	}
  Tp = &Texture->diff_entropy[0];
  results (Tp, F11, dentropy);

   if (feature_usage->meas_corr1)
	{ 
	icorr[0] = f12_icorr (P_matrix0, tones);
  	icorr[1] = f12_icorr (P_matrix45, tones);
  	icorr[2] = f12_icorr (P_matrix90, tones);
  	icorr[3] = f12_icorr (P_matrix135, tones);
	}
  Tp = &Texture->meas_corr1[0];
  results (Tp, F12, icorr);

  if (feature_usage->meas_corr2)
	{  
	icorr[0] = f13_icorr (P_matrix0, tones);
  	icorr[1] = f13_icorr (P_matrix45, tones);
  	icorr[2] = f13_icorr (P_matrix90, tones);
  	icorr[3] = f13_icorr (P_matrix135, tones);
	}
  Tp = &Texture->meas_corr2[0];
  results (Tp, F13, icorr);

  if (feature_usage->max_corr_coef)
	{  
	maxcorr[0] = f14_maxcorr (P_matrix0, tones);
	maxcorr[1] = f14_maxcorr (P_matrix45, tones);
	maxcorr[2] = f14_maxcorr (P_matrix90, tones);
        maxcorr[3] = f14_maxcorr (P_matrix135, tones);
	}
  /* M. Boland - 24 Nov 98 */
  else {
    maxcorr[0] = 0 ;
    maxcorr[1] = 0 ;
    maxcorr[2] = 0 ;
    maxcorr[3] = 0 ;
  }

  Tp = &Texture->max_corr_coef[0];
  results (Tp, F14, maxcorr);

  for (i=0; i<=tones; i++) free(P_matrix0[i]);
  for (i=0; i<=tones; i++) free(P_matrix45[i]);
  for (i=0; i<=tones; i++) free(P_matrix90[i]);
  for (i=0; i<=tones; i++) free(P_matrix135[i]);

  free(P_matrix0);
  free(P_matrix45);
  free(P_matrix90);
  free(P_matrix135);

/*  fprintf (stderr, " done.)\n"); */
  return (Texture);
 /* exit (0);*/
}

float f1_asm (P, Ng)
  float **P;
  int Ng;

/* Angular Second Moment */
{
  int i, j;
  float sum = 0;

  for (i = 0; i < Ng; ++i)
    for (j = 0; j < Ng; ++j)
      sum += P[i][j] * P[i][j];

  return sum;

  /*
   * The angular second-moment feature (ASM) f1 is a measure of homogeneity
   * of the image. In a homogeneous image, there are very few dominant
   * gray-tone transitions. Hence the P matrix for such an image will have
   * fewer entries of large magnitude.
   */
}


float f2_contrast (P, Ng)
  float **P;
  int Ng;

/* Contrast */
{
  int i, j, n;
  float sum = 0, bigsum = 0;

  for (n = 0; n < Ng; ++n)
  {
    for (i = 0; i < Ng; ++i)
      for (j = 0; j < Ng; ++j)
	if ((i - j) == n || (j - i) == n)
	  sum += P[i][j];
    bigsum += n * n * sum;

    sum = 0;
  }
  return bigsum;

  /*
   * The contrast feature is a difference moment of the P matrix and is a
   * measure of the contrast or the amount of local variations present in an
   * image.
   */
}

float f3_corr (P, Ng)
  float **P;
  int Ng;

/* Correlation */
{
  int i, j;
  float sum_sqrx = 0, sum_sqry = 0, tmp, *px;
  float meanx =0 , meany = 0 , stddevx, stddevy;

  px = pgm_vector (0, Ng);
  for (i = 0; i < Ng; ++i)
    px[i] = 0;

  /*
   * px[i] is the (i-1)th entry in the marginal probability matrix obtained
   * by summing the rows of p[i][j]
   */
  for (i = 0; i < Ng; ++i)
    for (j = 0; j < Ng; ++j)
      px[i] += P[i][j];


  /* Now calculate the means and standard deviations of px and py */
  /*- fix supplied by J. Michael Christensen, 21 Jun 1991 */
  /*- further modified by James Darrell McCauley, 16 Aug 1991 
   *     after realizing that meanx=meany and stddevx=stddevy
   */
  for (i = 0; i < Ng; ++i)
  {
    meanx += px[i]*i;
    sum_sqrx += px[i]*i*i;
  }
  /* M. Boland meanx = meanx/(sqrt(Ng)); */
  meany = meanx;
  sum_sqry = sum_sqrx;
  stddevx = sqrt (sum_sqrx - (meanx * meanx));
  stddevy = stddevx;

  /* Finally, the correlation ... */
  for (tmp = 0, i = 0; i < Ng; ++i)
    for (j = 0; j < Ng; ++j)
      tmp += i*j*P[i][j];

 free(px); 
 return (tmp - meanx * meany) / (stddevx * stddevy);
  /*
   * This correlation feature is a measure of gray-tone linear-dependencies
   * in the image.
   */
}


float f4_var (P, Ng)
  float **P;
  int Ng;

/* Sum of Squares: Variance */
{
  int i, j;
  float mean = 0, var = 0;

  /*- Corrected by James Darrell McCauley, 16 Aug 1991
   *  calculates the mean intensity level instead of the mean of
   *  cooccurrence matrix elements 
   */
  for (i = 0; i < Ng; ++i)
    for (j = 0; j < Ng; ++j)
      mean += i * P[i][j];

  for (i = 0; i < Ng; ++i)
    for (j = 0; j < Ng; ++j)
      /*  M. Boland - var += (i + 1 - mean) * (i + 1 - mean) * P[i][j]; */
      var += (i - mean) * (i - mean) * P[i][j];

  return var;
}

float f5_idm (P, Ng)
  float **P;
  int Ng;

/* Inverse Difference Moment */
{
  int i, j;
  float idm = 0;

  for (i = 0; i < Ng; ++i)
    for (j = 0; j < Ng; ++j)
      idm += P[i][j] / (1 + (i - j) * (i - j));

  return idm;
}

float Pxpy[2 * PGM_MAXMAXVAL];

float f6_savg (P, Ng)
  float **P;
  int Ng;

/* Sum Average */
{
  int i, j;
  extern float Pxpy[2 * PGM_MAXMAXVAL];
  float savg = 0;

  for (i = 0; i <= 2 * Ng; ++i)
    Pxpy[i] = 0;

  for (i = 0; i < Ng; ++i)
    for (j = 0; j < Ng; ++j)
      /* M. Boland Pxpy[i + j + 2] += P[i][j]; */
      /* Indexing from 2 instead of 0 is inconsistent with rest of code*/
      Pxpy[i + j] += P[i][j];
  /* M. Boland for (i = 2; i <= 2 * Ng; ++i) */
  /* Indexing from 2 instead of 0 is inconsistent with rest of code*/
  for (i = 0; i <= (2 * Ng - 2); ++i)
    savg += i * Pxpy[i];

  return savg;
}


float f7_svar (P, Ng, S)
  float **P, S;
  int Ng;

/* Sum Variance */
{
  int i, j;
  extern float Pxpy[2 * PGM_MAXMAXVAL];
  float var = 0;

  for (i = 0; i <= 2 * Ng; ++i)
    Pxpy[i] = 0;

  for (i = 0; i < Ng; ++i)
    for (j = 0; j < Ng; ++j)
      /* M. Boland Pxpy[i + j + 2] += P[i][j]; */
      /* Indexing from 2 instead of 0 is inconsistent with rest of code*/
      Pxpy[i + j] += P[i][j];

  /*  M. Boland for (i = 2; i <= 2 * Ng; ++i) */
  /* Indexing from 2 instead of 0 is inconsistent with rest of code*/
  for (i = 0; i <= (2 * Ng - 2); ++i)
    var += (i - S) * (i - S) * Pxpy[i];

  return var;
}

float f8_sentropy (P, Ng)
  float **P;
  int Ng;

/* Sum Entropy */
{
  int i, j;
  extern float Pxpy[2 * PGM_MAXMAXVAL];
  float sentropy = 0;

  for (i = 0; i <= 2 * Ng; ++i)
    Pxpy[i] = 0;

  for (i = 0; i < Ng; ++i)
    for (j = 0; j < Ng; ++j)
      Pxpy[i + j + 2] += P[i][j];

  for (i = 2; i <= 2 * Ng; ++i)
    /*  M. Boland  sentropy -= Pxpy[i] * log10 (Pxpy[i] + EPSILON); */
    sentropy -= Pxpy[i] * log10 (Pxpy[i] + EPSILON)/log10(2.0) ;

  return sentropy;
}


float f9_entropy (P, Ng)
  float **P;
  int Ng;

/* Entropy */
{
  int i, j;
  float entropy = 0;

  for (i = 0; i < Ng; ++i)
    for (j = 0; j < Ng; ++j)
      /*      entropy += P[i][j] * log10 (P[i][j] + EPSILON); */
      entropy += P[i][j] * log10 (P[i][j] + EPSILON)/log10(2.0) ;

  return -entropy; 
}


float f10_dvar (P, Ng)
  float **P;
  int Ng;

/* Difference Variance */
{
  int i, j, tmp;
  extern float Pxpy[2 * PGM_MAXMAXVAL];
  float sum = 0, sum_sqr = 0, var = 0;

  for (i = 0; i <= 2 * Ng; ++i)
    Pxpy[i] = 0;

  for (i = 0; i < Ng; ++i)
    for (j = 0; j < Ng; ++j)
      Pxpy[abs (i - j)] += P[i][j];

  /* Now calculate the variance of Pxpy (Px-y) */
  for (i = 0; i < Ng; ++i)
  {
    sum += i * Pxpy[i] ;
    sum_sqr += i * i * Pxpy[i] ;
    /* M. Boland sum += Pxpy[i];
    sum_sqr += Pxpy[i] * Pxpy[i];*/
  }
  /*tmp = Ng * Ng ;  M. Boland - wrong anyway, should be Ng */
  /*var = ((tmp * sum_sqr) - (sum * sum)) / (tmp * tmp); */
  
  var = sum_sqr - sum*sum ;

  return var;
}

float f11_dentropy (P, Ng)
  float **P;
  int Ng;

/* Difference Entropy */
{
  int i, j, tmp;
  extern float Pxpy[2 * PGM_MAXMAXVAL];
  float sum = 0, sum_sqr = 0, var = 0;

  for (i = 0; i <= 2 * Ng; ++i)
    Pxpy[i] = 0;

  for (i = 0; i < Ng; ++i)
    for (j = 0; j < Ng; ++j)
      Pxpy[abs (i - j)] += P[i][j];

  for (i = 0; i < Ng; ++i)
    /*    sum += Pxpy[i] * log10 (Pxpy[i] + EPSILON); */
    sum += Pxpy[i] * log10 (Pxpy[i] + EPSILON)/log10(2.0) ;

  return -sum;
}

float f12_icorr (P, Ng)
  float **P;
  int Ng;

/* Information Measures of Correlation */
/* All /log10(2.0) added by M. Boland */
{
  int i, j, tmp;
  float *px, *py;
  float hx = 0, hy = 0, hxy = 0, hxy1 = 0, hxy2 = 0;

  px = pgm_vector (0, Ng);
  py = pgm_vector (0, Ng);

  /*
   * px[i] is the (i-1)th entry in the marginal probability matrix obtained
   * by summing the rows of p[i][j]
   */
  for (i = 0; i < Ng; ++i)
  {
    for (j = 0; j < Ng; ++j)
    {
      px[i] += P[i][j];
      py[j] += P[i][j];
    }
  }

  for (i = 0; i < Ng; ++i)
    for (j = 0; j < Ng; ++j)
    {
      hxy1 -= P[i][j] * log10 (px[i] * py[j] + EPSILON)/log10(2.0);
      hxy2 -= px[i] * py[j] * log10 (px[i] * py[j] + EPSILON)/log10(2.0);
      hxy -= P[i][j] * log10 (P[i][j] + EPSILON)/log10(2.0);
    }

  /* Calculate entropies of px and py - is this right? */
  for (i = 0; i < Ng; ++i)
  {
    hx -= px[i] * log10 (px[i] + EPSILON)/log10(2.0);
    hy -= py[i] * log10 (py[i] + EPSILON)/log10(2.0);
  }
/*  fprintf(stderr,"hxy1=%f\thxy=%f\thx=%f\thy=%f\n",hxy1,hxy,hx,hy); */
  free(px);
  free(py);
  return ((hxy - hxy1) / (hx > hy ? hx : hy));
}

float f13_icorr (P, Ng)
  float **P;
  int Ng;

/* Information Measures of Correlation */
/* All /log10(2.0) added by M. Boland */

{
  int i, j;
  float *px, *py;
  float hx = 0, hy = 0, hxy = 0, hxy1 = 0, hxy2 = 0;

  px = pgm_vector (0, Ng);
  py = pgm_vector (0, Ng);

  /*
   * px[i] is the (i-1)th entry in the marginal probability matrix obtained
   * by summing the rows of p[i][j]
   */
  for (i = 0; i < Ng; ++i)
  {
    for (j = 0; j < Ng; ++j)
    {
      px[i] += P[i][j];
      py[j] += P[i][j];
    }
  }

  for (i = 0; i < Ng; ++i)
    for (j = 0; j < Ng; ++j)
    {
      hxy1 -= P[i][j] * log10 (px[i] * py[j] + EPSILON)/log10(2.0);
      hxy2 -= px[i] * py[j] * log10 (px[i] * py[j] + EPSILON)/log10(2.0);
      hxy -= P[i][j] * log10 (P[i][j] + EPSILON)/log10(2.0);
    }

  /* Calculate entropies of px and py */
  for (i = 0; i < Ng; ++i)
  {
    hx -= px[i] * log10 (px[i] + EPSILON)/log10(2.0);
    hy -= py[i] * log10 (py[i] + EPSILON)/log10(2.0);
  }
/*  fprintf(stderr,"hx=%f\thxy2=%f\n",hx,hxy2); */
  free(px);
  free(py);
  return (sqrt (abs (1 - exp (-2.0 * (hxy2 - hxy)))));
}

float f14_maxcorr (P, Ng)
  float **P;
  int Ng;

/* Returns the Maximal Correlation Coefficient */
{
  int i, j, k;
  float *px, *py, **Q;
  float *x, *iy, tmp;
  float f;

  px = pgm_vector (0, Ng);
  py = pgm_vector (0, Ng);
  Q = pgm_matrix (1, Ng + 1, 1, Ng + 1);
  x = pgm_vector (1, Ng);
  iy = pgm_vector (1, Ng);

  /*
   * px[i] is the (i-1)th entry in the marginal probability matrix obtained
   * by summing the rows of p[i][j]
   */
  for (i = 0; i < Ng; ++i)
  {
    for (j = 0; j < Ng; ++j)
    {
      px[i] += P[i][j];
      py[j] += P[i][j];
    }
  }

  /* Find the Q matrix */
  for (i = 0; i < Ng; ++i)
  {
    for (j = 0; j < Ng; ++j)
    {
      Q[i + 1][j + 1] = 0;
      for (k = 0; k < Ng; ++k)
	Q[i + 1][j + 1] += P[i][k] * P[j][k] / px[i] / py[k];
    }
  }

  /* Balance the matrix */
  mkbalanced (Q, Ng);
  /* Reduction to Hessenberg Form */
  reduction (Q, Ng);
  /* Finding eigenvalue for nonsymetric matrix using QR algorithm */
  if (!hessenberg (Q, Ng, x, iy))
	{ for (i=1; i<=Ng+1; i++) free(Q[i]+1);
	  free(Q+1);
	  free((char *)px);
	  free((char *)py);
	  free((x+1));
	  free((iy+1));
	  return 0.0;
	  /* fixed for Linux porting,
	   * I don't know what should be returned
	   */
	}
  /* simplesrt(Ng,x); */
  /* Returns the sqrt of the second largest eigenvalue of Q */
  for (i = 2, tmp = x[1]; i <= Ng; ++i)
    tmp = (tmp > x[i]) ? tmp : x[i];

  f = sqrt(x[Ng - 1]);

 for (i=1; i<=Ng+1; i++) free(Q[i]+1);
 free(Q+1);
 free((char *)px); 
 free((char *)py); 
 free((x+1)); 
 free((iy+1)); 

 return f;
}

float *pgm_vector (nl, nh)
  int nl, nh;
{
  float *v;
  int    i;

  v = (float *) malloc ((unsigned) (nh - nl + 1) * sizeof (float));
  if (!v)
    fprintf (stderr, "memory allocation failure (pgm_vector) "), exit (1);

  for (i=0; i<=(nh-nl); i++) v[i]=0;
  return v - nl;
}


float **pgm_matrix (nrl, nrh, ncl, nch)
  int nrl, nrh, ncl, nch;

/* Allocates a float matrix with range [nrl..nrh][ncl..nch] */
{
  int i;
  float **m;

  /* allocate pointers to rows */
  m = (float **) malloc ((unsigned) (nrh - nrl + 1) * sizeof (float *));
  if (!m)
    fprintf (stderr, "memory allocation failure (pgm_matrix 1) "), exit (1);
  m -= ncl;

  /* allocate rows and set pointers to them */
  for (i = nrl; i <= nrh; i++)
  {
    m[i] = (float *) malloc ((unsigned) (nch - ncl + 1) * sizeof (float));
    if (!m[i])
      fprintf (stderr, "memory allocation failure (pgm_matrix 2) "), 
	exit (2);
    m[i] -= ncl;
  }
  /* return pointer to array of pointers to rows */
  return m;
}

void results (Tp, c, a)
  float *Tp;
  char *c;
  float *a;
{
  int i;
  float max, min;
  max = a[0];
  min = a[0];
/*  DOT;
  fprintf (stdout, "%s", c);
*/  for (i = 0; i < 4; ++i, *Tp++)
    {	
    if (a[i] <= min)
	min = a[i];
    if (a[i] > max)
	max = a[i];
  /*  fprintf (stdout, "% 1.3e ", a[i]); */
    *Tp = a[i];
    }	
/*  fprintf (stdout, "% 1.3e  % 1.3e\n", 
    (a[0] + a[1] + a[2] + a[3]) / 4,max-min); */
  *Tp = (a[0] + a[1] + a[2] + a[3]) / 4;
  *Tp++;
  *Tp = max - min;
 
  	
}

void simplesrt (n, arr)
  int n;
  float arr[];
{
  int i, j;
  float a;

  for (j = 2; j <= n; j++)
  {
    a = arr[j];
    i = j - 1;
    while (i > 0 && arr[i] > a)
    {
      arr[i + 1] = arr[i];
      i--;
    }
    arr[i + 1] = a;
  }
}

void mkbalanced (a, n)
  float **a;
  int n;
{
  int last, j, i;
  float s, r, g, f, c, sqrdx;

  sqrdx = RADIX * RADIX;
  last = 0;
  while (last == 0)
  {
    last = 1;
    for (i = 1; i <= n; i++)
    {
      r = c = 0.0;
      for (j = 1; j <= n; j++)
	if (j != i)
	{
	  c += fabs (a[j][i]);
	  r += fabs (a[i][j]);
	}
      if (c && r)
      {
	g = r / RADIX;
	f = 1.0;
	s = c + r;
	while (c < g)
	{
	  f *= RADIX;
	  c *= sqrdx;
	}
	g = r * RADIX;
	while (c > g)
	{
	  f /= RADIX;
	  c /= sqrdx;
	}
	if ((c + r) / f < 0.95 * s)
	{
	  last = 0;
	  g = 1.0 / f;
	  for (j = 1; j <= n; j++)
	    a[i][j] *= g;
	  for (j = 1; j <= n; j++)
	    a[j][i] *= f;
	}
      }
    }
  }
}


void reduction (a, n)
  float **a;
  int n;
{
  int m, j, i;
  float y, x;

  for (m = 2; m < n; m++)
  {
    x = 0.0;
    i = m;
    for (j = m; j <= n; j++)
    {
      if (fabs (a[j][m - 1]) > fabs (x))
      {
	x = a[j][m - 1];
	i = j;
      }
    }
    if (i != m)
    {
      for (j = m - 1; j <= n; j++)
	SWAP (a[i][j], a[m][j])  
	for (j = 1; j <= n; j++)
	  SWAP (a[j][i], a[j][m]) 
	  a[j][i] = a[j][i];
    }
    if (x)
    {
      for (i = m + 1; i <= n; i++)
      {
	if (y = a[i][m - 1])
	{
	  y /= x;
	  a[i][m - 1] = y;
	  for (j = m; j <= n; j++)
	    a[i][j] -= y * a[m][j];
	  for (j = 1; j <= n; j++)
	    a[j][m] += y * a[j][i];
	}
      }
    }
  }
}

int hessenberg (a, n, wr, wi)
  float **a, wr[], wi[];
  int n;

{
  int nn, m, l, k, j, its, i, mmin;
  float z, y, x, w, v, u, t, s, r, q, p, anorm;

  anorm = fabs (a[1][1]);
  for (i = 2; i <= n; i++)
    for (j = (i - 1); j <= n; j++)
      anorm += fabs (a[i][j]);
  nn = n;
  t = 0.0;
  while (nn >= 1)
  {
    its = 0;
    do
    {
      for (l = nn; l >= 2; l--)
      {
	s = fabs (a[l - 1][l - 1]) + fabs (a[l][l]);
	if (s == 0.0)
	  s = anorm;
	if ((float) (fabs (a[l][l - 1]) + s) == s)
	  break;
      }
      x = a[nn][nn];
      if (l == nn)
      {
	wr[nn] = x + t;
	wi[nn--] = 0.0;
      }
      else
      {
	y = a[nn - 1][nn - 1];
	w = a[nn][nn - 1] * a[nn - 1][nn];
	if (l == (nn - 1))
	{
	  p = 0.5 * (y - x);
	  q = p * p + w;
	  z = sqrt (fabs (q));
	  x += t;
	  if (q >= 0.0)
	  {
	    z = p + SIGN (z, p); 
	    wr[nn - 1] = wr[nn] = x + z;
	    if (z)
	      wr[nn] = x - w / z;
	    wi[nn - 1] = wi[nn] = 0.0;
	  }
	  else
	  {
	    wr[nn - 1] = wr[nn] = x + p;
	    wi[nn - 1] = -(wi[nn] = z);
	  }
	  nn -= 2;
	}
	else
	{
	  if (its == 30)
	    {
/*	    fprintf (stderr, 
  "Too many iterations to required to find %s\ngiving up\n", F14);  */
	     return 0; /*exit (1);*/
	     }			
	  if (its == 10 || its == 20)
	  {
	    t += x;
	    for (i = 1; i <= nn; i++)
	      a[i][i] -= x;
	    s = fabs (a[nn][nn - 1]) + fabs (a[nn - 1][nn - 2]);
	    y = x = 0.75 * s;
	    w = -0.4375 * s * s;
	  }
	  ++its;
	  for (m = (nn - 2); m >= l; m--)
	  {
	    z = a[m][m];
	    r = x - z;
	    s = y - z;
	    p = (r * s - w) / a[m + 1][m] + a[m][m + 1];
	    q = a[m + 1][m + 1] - z - r - s;
	    r = a[m + 2][m + 1];
	    s = fabs (p) + fabs (q) + fabs (r);
	    p /= s;
	    q /= s;
	    r /= s;
	    if (m == l)
	      break;
	    u = fabs (a[m][m - 1]) * (fabs (q) + fabs (r));
	    v = fabs (p) * (fabs (a[m - 1][m - 1]) + 
			    fabs (z) + fabs (a[m + 1][m + 1]));
	    if ((float) (u + v) == v)
	      break;
	  }
	  for (i = m + 2; i <= nn; i++)
	  {
	    a[i][i - 2] = 0.0;
	    if (i != (m + 2))
	      a[i][i - 3] = 0.0;
	  }
	  for (k = m; k <= nn - 1; k++)
	  {
	    if (k != m)
	    {
	      p = a[k][k - 1];
	      q = a[k + 1][k - 1];
	      r = 0.0;
	      if (k != (nn - 1))
		r = a[k + 2][k - 1];
	      if (x = fabs (p) + fabs (q) + fabs (r))
	      {
		p /= x;
		q /= x;
		r /= x;
	      }
	    }
	    if (s = SIGN (sqrt (p * p + q * q + r * r), p)) 
	    {
	      if (k == m)
	      {
		if (l != m)
		  a[k][k - 1] = -a[k][k - 1];
	      }
	      else
		a[k][k - 1] = -s * x;
	      p += s;
	      x = p / s;
	      y = q / s;
	      z = r / s;
	      q /= p;
	      r /= p;
	      for (j = k; j <= nn; j++)
	      {
		p = a[k][j] + q * a[k + 1][j];
		if (k != (nn - 1))
		{
		  p += r * a[k + 2][j];
		  a[k + 2][j] -= p * z;
		}
		a[k + 1][j] -= p * y;
		a[k][j] -= p * x;
	      }
	      mmin = nn < k + 3 ? nn : k + 3;
	      for (i = l; i <= mmin; i++)
	      {
		p = x * a[i][k] + y * a[i][k + 1];
		if (k != (nn - 1))
		{
		  p += z * a[i][k + 2];
		  a[i][k + 2] -= p * r;
		}
		a[i][k + 1] -= p * q;
		a[i][k] -= p;
	      }
	    }
	  }
	}
      }
    } while (l < nn - 1);
  }
return 1;
}
 



Copyright ©1999 Michael V. Boland
1999-09-18