jdm5548@tamsun.tamu.edu (James Darrell McCauley) (04/30/91)
This is pgmtxtur - an addition to Jef Poskanser's PBMPLUS. It calculates
14 textural features of a grayscale image using the (statistical) approach
found in:
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.
You need the PBMPLUS package to compile.
James Darrell McCauley, Grad Res Asst, Spatial Analysis Lab
Dept of Ag Engr, Texas A&M Univ, College Station, TX 77843-2117, USA
(jdm5548@diamond.tamu.edu, jdm5548@tamagen.bitnet)
--cut here--
#!/bin/sh
# This is a shell archive (produced by shar 3.49)
# To extract the files from this archive, save it to a file, remove
# everything above the "!/bin/sh" line above, and type "sh file_name".
#
# made 04/30/1991 00:29 UTC by jdm5548@amber
# Source directory /u/class/jdm5548/tmp
#
# existing files will NOT be overwritten unless -c is specified
#
# This shar contains:
# length mode name
# ------ ---------- ------------------------------------------
# 503 -rw-r--r-- README
# 2849 -rw-r--r-- pgmtxtur.1
# 24184 -rw-r--r-- pgmtxtur.c
#
# ============= README ==============
if test -f 'README' -a X"$1" != X"-c"; then
echo 'x - skipping README (File already exists)'
else
echo 'x - extracting README (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'README' &&
This is pgmtxtur - an addition to Jef Poskanser's PBMPLUS. It calculates
14 textural features of a grayscale image using the (statistical) approach
found in:
X
X Haralick, R.M., K. Shanmugam, and I. Dinstein. 1973. Textural features
X for image classification. IEEE Transactions on Systems, Man, and
X Cybertinetics, SMC-3(6):610-621.
X
You need the PBMPLUS package to compile:
X 1. Add pgmtxtur as a math target in pbmplus/pgm/Makefile
X and put the sources in that directory.
X 2. 'make pgmtxtur'
X
SHAR_EOF
chmod 0644 README ||
echo 'restore of README failed'
Wc_c="`wc -c < 'README'`"
test 503 -eq "$Wc_c" ||
echo 'README: original size 503, current size' "$Wc_c"
fi
# ============= pgmtxtur.1 ==============
if test -f 'pgmtxtur.1' -a X"$1" != X"-c"; then
echo 'x - skipping pgmtxtur.1 (File already exists)'
else
echo 'x - extracting pgmtxtur.1 (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'pgmtxtur.1' &&
.TH pgmtxtur 1 "27 Apr 1991"
.SH NAME
pgmtxtur - calculate textural features on a portable graymap
.SH SYNOPSIS
.B pgmtxtur
.RB [ -d
.IR d ]
.RI [ pgmfile ]
.SH DESCRIPTION
Reads a portable graymap as input. Calculates textural features
based on spatial dependence matrices at 0, 45, 90, and 135 degrees for
a distance
.IR d
(default = 1). Textural features include:
.IP
(1) Angular Second Moment,
.br
(2) Contrast,
.br
(3) Correlation,
.br
(4) Variance,
.br
(5) Inverse Difference Moment,
.br
(6) Sum Average,
.br
(7) Sum Variance,
.br
(8) Sum Entropy,
.br
(9) Entropy,
.br
(10) Difference Variance,
.br
(11) Difference Entropy,
.br
(12, 13) Information Measures of Correlation, and
.br
(14) Maximal Correlation Coefficient.
.PP
Algorithm taken from:
.br
Haralick, R.M., K. Shanmugam, and I. Dinstein. 1973. Textural features
for image classification.
.I IEEE Transactions on Systems, Man,
.I and Cybertinetics,
SMC-3(6):610-621.
.SH BUGS
The program can run incredibly slow for large images (larger than 64 x 64)
and command line options are limited.
The method for finding (14) the maximal correlation coefficient, which
requires finding the second largest eigenvalue of a matrix Q, does not
always converge.
.SH "SEE ALSO"
pgm(5), pnmcut(1),
.I IEEE Transactions on Systems, Man,
.I and Cybertinetics,
SMC-3(6):610-621.
.SH AUTHOR
Copyright (C) 1991 by 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.
.\"
.\" Original copyright notice follows:
.\" Copyright (C) 1991 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.
X
SHAR_EOF
chmod 0644 pgmtxtur.1 ||
echo 'restore of pgmtxtur.1 failed'
Wc_c="`wc -c < 'pgmtxtur.1'`"
test 2849 -eq "$Wc_c" ||
echo 'pgmtxtur.1: original size 2849, current size' "$Wc_c"
fi
# ============= pgmtxtur.c ==============
if test -f 'pgmtxtur.c' -a X"$1" != X"-c"; then
echo 'x - skipping pgmtxtur.c (File already exists)'
else
echo 'x - extracting pgmtxtur.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'pgmtxtur.c' &&
/*-pgmtxtur.c - calculate textural features on a portable graymap. March 1991.
**
** 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 Poskanser (Copyright (C) 1991 by Jef Poskanser), whose original
** copyright notice follows:
**
** 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.
**
** 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.
*/
X
#include "pgm.h"
#include<math.h>
X
#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 "
X
#define SIGN(x,y) ((y)<0 ? -fabs(x) : fabs(x))
#define DOT fprintf(stderr,".")
#define SWAP(a,b) {y=(a);(a)=(b);(b)=y;}
X
void results (), hessenberg (), mkbalanced (), reduction (), simplesrt ();
float f1_asm (), f2_contrast (), f3_corr (), f4_var (), f5_idm (),
X f6_savg (), f7_svar (), f8_sentropy (), f9_entropy (), f10_dvar (),
X f11_dentropy (), f12_icorr (), f13_icorr (), f14_maxcorr (), *vector (),
X **matrix ();
X
X
void main (argc, argv)
X int argc;
X char *argv[];
{
X FILE *ifp;
X register gray **grays, *gP;
X int tone[PGM_MAXMAXVAL], R0, R45, R90, angle, d = 1, x, y;
X int argn, rows, cols, bps, padright, row, col;
X int itone, jtone, tones;
X float **P_matrix0, **P_matrix45, **P_matrix90, **P_matrix135;
X float ASM[4], contrast[4], corr[4], var[4], idm[4], savg[4];
X float sentropy[4], svar[4], entropy[4], dvar[4], dentropy[4];
X float icorr[4], maxcorr[4];
X gray maxval, nmaxval;
X char *usage = "[-d <d>] [pgmfile]";
X
X pgm_init( &argc, argv );
X
X argn = 1;
X
X /* Check for flags. */
X if ( argn < argc && argv[argn][0] == '-' )
X {
X if ( argv[argn][1] == 'd' )
X {
X ++argn;
X if ( argn == argc || sscanf( argv[argn], "%d", &d ) != 1 )
X pm_usage( usage );
X }
X else
X pm_usage( usage );
X ++argn;
X }
X
X if ( argn < argc )
X {
X ifp = pm_openr( argv[argn] );
X ++argn;
X }
X else
X ifp = stdin;
X
X if ( argn != argc )
X pm_usage( usage );
X
X
X grays = pgm_readpgm (ifp, &cols, &rows, &maxval);
X pm_close (ifp);
X
X /* Determine the number of different gray scales (not maxval) */
X for (row = PGM_MAXMAXVAL; row >= 0; --row)
X tone[row] = -1;
X for (row = rows - 1; row >= 0; --row)
X for (col = 0; col < cols; ++col)
X tone[grays[row][col]] = grays[row][col];
X for (row = PGM_MAXMAXVAL, tones = 0; row >= 0; --row)
X if (tone[row] != -1)
X tones++;
X fprintf (stderr, "(Image has %d graylevels.)\n", tones);
X
X /* Collapse array, taking out all zero values */
X for (row = 0, itone = 0; row <= PGM_MAXMAXVAL; row++)
X if (tone[row] != -1)
X tone[itone++] = tone[row];
X /* Now array contains only the gray levels present (in ascending order) */
X
X /* Allocate memory for gray-tone spatial dependence matrix */
X P_matrix0 = matrix (0, tones, 0, tones);
X P_matrix45 = matrix (0, tones, 0, tones);
X P_matrix90 = matrix (0, tones, 0, tones);
X P_matrix135 = matrix (0, tones, 0, tones);
X for (row = 0; row < tones; ++row)
X for (col = 0; col < tones; ++col)
X {
X P_matrix0[row][col] = P_matrix45[row][col] = 0;
X P_matrix90[row][col] = P_matrix135[row][col] = 0;
X }
X
X /* Find gray-tone spatial dependence matrix */
X fprintf (stderr, "(Computing spatial dependence matrix...");
X for (row = 0; row < rows; ++row)
X for (col = 0; col < cols; ++col)
X for (x = 0, angle = 0; angle <= 135; angle += 45)
X {
X while (tone[x] != grays[row][col])
X x++;
X if (angle == 0 && col + d < cols)
X {
X y = 0;
X while (tone[y] != grays[row][col + d])
X y++;
X P_matrix0[x][y]++;
X P_matrix0[y][x]++;
X }
X if (angle == 90 && row + d < rows)
X {
X y = 0;
X while (tone[y] != grays[row + d][col])
X y++;
X P_matrix90[x][y]++;
X P_matrix90[y][x]++;
X }
X if (angle == 45 && row + d < rows && col - d >= 0)
X {
X y = 0;
X while (tone[y] != grays[row + d][col - d])
X y++;
X P_matrix45[x][y]++;
X P_matrix45[y][x]++;
X }
X if (angle == 135 && row + d < rows && col + d < cols)
X {
X y = 0;
X while (tone[y] != grays[row + d][col + d])
X y++;
X P_matrix135[x][y]++;
X P_matrix135[y][x]++;
X }
X }
X /* Gray-tone spatial dependence matrices are complete */
X
X /* Find normalizing constants */
X R0 = 2 * rows * (cols - 1);
X R45 = 2 * (rows - 1) * (cols - 1);
X R90 = 2 * (rows - 1) * cols;
X
X /* Normalize gray-tone spatial dependence matrix */
X for (itone = 0; itone < tones; ++itone)
X for (jtone = 0; jtone < tones; ++jtone)
X {
X P_matrix0[itone][jtone] /= R0;
X P_matrix45[itone][jtone] /= R45;
X P_matrix90[itone][jtone] /= R90;
X P_matrix135[itone][jtone] /= R45;
X }
X
X fprintf (stderr, " done.)\n");
X fprintf (stderr, "(Computing textural features");
X fprintf (stdout, "\n");
X DOT;
X fprintf (stdout,
X "%s 0 45 90 135 Avg\n",
X BL);
X ASM[0] = f1_asm (P_matrix0, tones);
X ASM[1] = f1_asm (P_matrix45, tones);
X ASM[2] = f1_asm (P_matrix90, tones);
X ASM[3] = f1_asm (P_matrix135, tones);
X results (F1, ASM);
X
X contrast[0] = f2_contrast (P_matrix0, tones);
X contrast[1] = f2_contrast (P_matrix45, tones);
X contrast[2] = f2_contrast (P_matrix90, tones);
X contrast[3] = f2_contrast (P_matrix135, tones);
X results (F2, contrast);
X
X corr[0] = f3_corr (P_matrix0, tones);
X corr[1] = f3_corr (P_matrix45, tones);
X corr[2] = f3_corr (P_matrix90, tones);
X corr[3] = f3_corr (P_matrix135, tones);
X results (F3, corr);
X
X var[0] = f4_var (P_matrix0, tones);
X var[1] = f4_var (P_matrix45, tones);
X var[2] = f4_var (P_matrix90, tones);
X var[3] = f4_var (P_matrix135, tones);
X results (F4, var);
X
X idm[0] = f5_idm (P_matrix0, tones);
X idm[1] = f5_idm (P_matrix45, tones);
X idm[2] = f5_idm (P_matrix90, tones);
X idm[3] = f5_idm (P_matrix135, tones);
X results (F5, idm);
X
X savg[0] = f6_savg (P_matrix0, tones);
X savg[1] = f6_savg (P_matrix45, tones);
X savg[2] = f6_savg (P_matrix90, tones);
X savg[3] = f6_savg (P_matrix135, tones);
X results (F6, savg);
X
X sentropy[0] = f8_sentropy (P_matrix0, tones);
X sentropy[1] = f8_sentropy (P_matrix45, tones);
X sentropy[2] = f8_sentropy (P_matrix90, tones);
X sentropy[3] = f8_sentropy (P_matrix135, tones);
X svar[0] = f7_svar (P_matrix0, tones, sentropy[0]);
X svar[1] = f7_svar (P_matrix45, tones, sentropy[1]);
X svar[2] = f7_svar (P_matrix90, tones, sentropy[2]);
X svar[3] = f7_svar (P_matrix135, tones, sentropy[3]);
X results (F7, svar);
X results (F8, sentropy);
X
X entropy[0] = f9_entropy (P_matrix0, tones);
X entropy[1] = f9_entropy (P_matrix45, tones);
X entropy[2] = f9_entropy (P_matrix90, tones);
X entropy[3] = f9_entropy (P_matrix135, tones);
X results (F9, entropy);
X
X dvar[0] = f10_dvar (P_matrix0, tones);
X dvar[1] = f10_dvar (P_matrix45, tones);
X dvar[2] = f10_dvar (P_matrix90, tones);
X dvar[3] = f10_dvar (P_matrix135, tones);
X results (F10, dvar);
X
X dentropy[0] = f11_dentropy (P_matrix0, tones);
X dentropy[1] = f11_dentropy (P_matrix45, tones);
X dentropy[2] = f11_dentropy (P_matrix90, tones);
X dentropy[3] = f11_dentropy (P_matrix135, tones);
X results (F11, dentropy);
X
X icorr[0] = f12_icorr (P_matrix0, tones);
X icorr[1] = f12_icorr (P_matrix45, tones);
X icorr[2] = f12_icorr (P_matrix90, tones);
X icorr[3] = f12_icorr (P_matrix135, tones);
X results (F12, icorr);
X
X icorr[0] = f13_icorr (P_matrix0, tones);
X icorr[1] = f13_icorr (P_matrix45, tones);
X icorr[2] = f13_icorr (P_matrix90, tones);
X icorr[3] = f13_icorr (P_matrix135, tones);
X results (F13, icorr);
X
X maxcorr[0] = f14_maxcorr (P_matrix0, tones);
X maxcorr[1] = f14_maxcorr (P_matrix45, tones);
X maxcorr[2] = f14_maxcorr (P_matrix90, tones);
X maxcorr[3] = f14_maxcorr (P_matrix135, tones);
X results (F14, maxcorr);
X
X
X fprintf (stderr, " done.)\n");
X exit (0);
}
X
float f1_asm (P, Ng)
X float **P;
X int Ng;
X
/* Angular Second Moment */
{
X int i, j;
X float sum = 0;
X
X for (i = 0; i < Ng; ++i)
X for (j = 0; j < Ng; ++j)
X sum += P[i][j] * P[i][j];
X
X return sum;
X
X /*
X * The angular second-moment feature (ASM) f1 is a measure of homogeneity
X * of the image. In a homogeneous image, there are very few dominant
X * gray-tone transitions. Hence the P matrix for such an image will have
X * fewer entries of large magnitude.
X */
}
X
X
float f2_contrast (P, Ng)
X float **P;
X int Ng;
X
/* Contrast */
{
X int i, j, n;
X float sum = 0, bigsum = 0;
X
X for (n = 0; n < Ng; ++n)
X {
X for (i = 0; i < Ng; ++i)
X for (j = 0; j < Ng; ++j)
X if ((i - j) == n || (j - i) == n)
X sum += P[i][j];
X bigsum += n * n * sum;
X
X sum = 0;
X }
X return bigsum;
X
X /*
X * The contrast feature is a difference moment of the P matrix and is a
X * measure of the contrast or the amount of local variations present in an
X * image.
X */
}
X
float f3_corr (P, Ng)
X float **P;
X int Ng;
X
/* Correlation */
{
X int i, j;
X float sumx = 0, sum_sqrx = 0, sumy = 0, sum_sqry = 0, tmp, *px, *py;
X float meanx, meany, stddevx, stddevy;
X float margpx, margpy;
X
X px = vector (0, Ng);
X py = vector (0, Ng);
X for (i = 0; i < Ng; ++i)
X px[i] = py[i] = 0;
X
X /*
X * px[i] is the (i-1)th entry in the marginal probability matrix obtained
X * by summing the rows of p[i][j]
X */
X for (i = 0; i < Ng; ++i)
X {
X for (j = 0; j < Ng; ++j)
X {
X px[i] += P[i][j];
X py[j] += P[i][j];
X }
X }
X
X /* Now calculate the means and standard deviations of px and py */
X for (i = 0; i < Ng; ++i)
X {
X sumx += px[i];
X sum_sqrx += (px[i] * px[i]);
X sumy += py[i];
X sum_sqry += (py[i] * py[i]);
X }
X
X tmp = Ng * Ng;
X meanx = sumx / tmp;
X stddevx = sqrt (((tmp * sum_sqrx) - (sumx * sumx)) / (tmp * tmp - 1));
X meany = sumy / tmp;
X stddevy = sqrt (((tmp * sum_sqry) - (sumy * sumy)) / (tmp * tmp - 1));
X
X /* Finally, the correlation ... */
X for (tmp = 0, i = 0; i < Ng; ++i)
X for (j = 0; j < Ng; ++j)
X tmp += (i + 1) * (j + 1) * P[i][j];
X
X if (tmp - meanx * meany)
X return ((tmp - meanx * meany) / (stddevx * stddevy));
X else
X return 0;
X /*
X * This correlation feature is a measure of gray-tone linear-dependencies
X * in the image.
X */
}
X
float f4_var (P, Ng)
X float **P;
X int Ng;
X
/* Sum of Squares: Variance */
{
X int i, j;
X float mean = 0, var = 0;
X
X /* Find population mean */
X for (i = 0; i < Ng; ++i)
X for (j = 0; j < Ng; ++j)
X mean += P[i][j];
X mean /= (Ng * Ng);
X
X /* fprintf(stderr,"f4:mean=%f\n",mean); */
X for (i = 0; i < Ng; ++i)
X for (j = 0; j < Ng; ++j)
X var += (i + 1 - mean) * (i + 1 - mean) * P[i][j];
X
X return var;
}
X
float f5_idm (P, Ng)
X float **P;
X int Ng;
X
/* Inverse Difference Moment */
{
X int i, j;
X float idm = 0;
X
X for (i = 0; i < Ng; ++i)
X for (j = 0; j < Ng; ++j)
X idm += P[i][j] / (1 + (i - j) * (i - j));
X
X return idm;
}
X
float Pxpy[2 * PGM_MAXMAXVAL];
X
float f6_savg (P, Ng)
X float **P;
X int Ng;
X
/* Sum Average */
{
X int i, j;
X extern float Pxpy[2 * PGM_MAXMAXVAL];
X float savg = 0;
X
X for (i = 0; i <= 2 * Ng; ++i)
X Pxpy[i] = 0;
X
X for (i = 0; i < Ng; ++i)
X for (j = 0; j < Ng; ++j)
X Pxpy[i + j + 2] += P[i][j];
X for (i = 2; i <= 2 * Ng; ++i)
X savg += i * Pxpy[i];
X
X return savg;
}
X
X
float f7_svar (P, Ng, S)
X float **P, S;
X int Ng;
X
/* Sum Variance */
{
X int i, j;
X extern float Pxpy[2 * PGM_MAXMAXVAL];
X float var = 0;
X
X for (i = 0; i <= 2 * Ng; ++i)
X Pxpy[i] = 0;
X
X for (i = 0; i < Ng; ++i)
X for (j = 0; j < Ng; ++j)
X Pxpy[i + j + 2] += P[i][j];
X
X for (i = 2; i <= 2 * Ng; ++i)
X var += (i - S) * (i - S) * Pxpy[i];
X
X return var;
}
X
float f8_sentropy (P, Ng)
X float **P;
X int Ng;
X
/* Sum Entropy */
{
X int i, j;
X extern float Pxpy[2 * PGM_MAXMAXVAL];
X float sentropy = 0;
X
X for (i = 0; i <= 2 * Ng; ++i)
X Pxpy[i] = 0;
X
X for (i = 0; i < Ng; ++i)
X for (j = 0; j < Ng; ++j)
X Pxpy[i + j + 2] += P[i][j];
X
X for (i = 2; i <= 2 * Ng; ++i)
X sentropy -= Pxpy[i] * log10 (Pxpy[i] + EPSILON);
X
X return sentropy;
}
X
X
float f9_entropy (P, Ng)
X float **P;
X int Ng;
X
/* Entropy */
{
X int i, j;
X float entropy = 0;
X
X for (i = 0; i < Ng; ++i)
X for (j = 0; j < Ng; ++j)
X entropy += P[i][j] * log10 (P[i][j] + EPSILON);
X
X return -entropy;
}
X
X
float f10_dvar (P, Ng)
X float **P;
X int Ng;
X
/* Difference Variance */
{
X int i, j, tmp;
X extern float Pxpy[2 * PGM_MAXMAXVAL];
X float sum = 0, sum_sqr = 0, var = 0;
X
X for (i = 0; i <= 2 * Ng; ++i)
X Pxpy[i] = 0;
X
X for (i = 0; i < Ng; ++i)
X for (j = 0; j < Ng; ++j)
X Pxpy[abs (i - j)] += P[i][j];
X
X /* Now calculate the variance of Pxpy (Px-y) */
X for (i = 0; i < Ng; ++i)
X {
X sum += Pxpy[i];
X sum_sqr += Pxpy[i] * Pxpy[i];
X }
X tmp = Ng * Ng;
X var = ((tmp * sum_sqr) - (sum * sum)) / (tmp * tmp);
X
X return var;
}
X
float f11_dentropy (P, Ng)
X float **P;
X int Ng;
X
/* Difference Entropy */
{
X int i, j, tmp;
X extern float Pxpy[2 * PGM_MAXMAXVAL];
X float sum = 0, sum_sqr = 0, var = 0;
X
X for (i = 0; i <= 2 * Ng; ++i)
X Pxpy[i] = 0;
X
X for (i = 0; i < Ng; ++i)
X for (j = 0; j < Ng; ++j)
X Pxpy[abs (i - j)] += P[i][j];
X
X for (i = 0; i < Ng; ++i)
X sum += Pxpy[i] * log10 (Pxpy[i] + EPSILON);
X
X return -sum;
}
X
float f12_icorr (P, Ng)
X float **P;
X int Ng;
X
/* Information Measures of Correlation */
{
X int i, j, tmp;
X float *px, *py;
X float hx = 0, hy = 0, hxy = 0, hxy1 = 0, hxy2 = 0;
X
X px = vector (0, Ng);
X py = vector (0, Ng);
X
X /*
X * px[i] is the (i-1)th entry in the marginal probability matrix obtained
X * by summing the rows of p[i][j]
X */
X for (i = 0; i < Ng; ++i)
X {
X for (j = 0; j < Ng; ++j)
X {
X px[i] += P[i][j];
X py[j] += P[i][j];
X }
X }
X
X for (i = 0; i < Ng; ++i)
X for (j = 0; j < Ng; ++j)
X {
X hxy1 -= P[i][j] * log10 (px[i] * py[j] + EPSILON);
X hxy2 -= px[i] * py[j] * log10 (px[i] * py[j] + EPSILON);
X hxy -= P[i][j] * log10 (P[i][j] + EPSILON);
X }
X
X /* Calculate entropies of px and py - is this right? */
X for (i = 0; i < Ng; ++i)
X {
X hx -= px[i] * log10 (px[i] + EPSILON);
X hy -= py[i] * log10 (py[i] + EPSILON);
X }
X fprintf(stderr,"hxy1=%f\thxy=%f\thx=%f\thy=%f\n",hxy1,hxy,hx,hy);
X return ((hxy - hxy1) / (hx > hy ? hx : hy));
}
X
float f13_icorr (P, Ng)
X float **P;
X int Ng;
X
/* Information Measures of Correlation */
{
X int i, j;
X float *px, *py;
X float hx = 0, hy = 0, hxy = 0, hxy1 = 0, hxy2 = 0;
X
X px = vector (0, Ng);
X py = vector (0, Ng);
X
X /*
X * px[i] is the (i-1)th entry in the marginal probability matrix obtained
X * by summing the rows of p[i][j]
X */
X for (i = 0; i < Ng; ++i)
X {
X for (j = 0; j < Ng; ++j)
X {
X px[i] += P[i][j];
X py[j] += P[i][j];
X }
X }
X
X for (i = 0; i < Ng; ++i)
X for (j = 0; j < Ng; ++j)
X {
X hxy1 -= P[i][j] * log10 (px[i] * py[j] + EPSILON);
X hxy2 -= px[i] * py[j] * log10 (px[i] * py[j] + EPSILON);
X hxy -= P[i][j] * log10 (P[i][j] + EPSILON);
X }
X
X /* Calculate entropies of px and py */
X for (i = 0; i < Ng; ++i)
X {
X hx -= px[i] * log10 (px[i] + EPSILON);
X hy -= py[i] * log10 (py[i] + EPSILON);
X }
X fprintf(stderr,"hx=%f\thxy2=%f\n",hx,hxy2);
X return (sqrt (abs (1 - exp (-2.0 * (hxy2 - hxy)))));
}
X
float f14_maxcorr (P, Ng)
X float **P;
X int Ng;
X
/* Returns the Maximal Correlation Coefficient */
{
X int i, j, k;
X float *px, *py, **Q;
X float *x, *iy, tmp;
X
X px = vector (0, Ng);
X py = vector (0, Ng);
X Q = matrix (1, Ng + 1, 1, Ng + 1);
X x = vector (1, Ng);
X iy = vector (1, Ng);
X
X /*
X * px[i] is the (i-1)th entry in the marginal probability matrix obtained
X * by summing the rows of p[i][j]
X */
X for (i = 0; i < Ng; ++i)
X {
X for (j = 0; j < Ng; ++j)
X {
X px[i] += P[i][j];
X py[j] += P[i][j];
X }
X }
X
X /* Find the Q matrix */
X for (i = 0; i < Ng; ++i)
X {
X for (j = 0; j < Ng; ++j)
X {
X Q[i + 1][j + 1] = 0;
X for (k = 0; k < Ng; ++k)
X Q[i + 1][j + 1] += P[i][k] * P[j][k] / px[i] / py[k];
X }
X }
X
X /* Balance the matrix */
X mkbalanced (Q, Ng);
X /* Reduction to Hessenberg Form */
X reduction (Q, Ng);
X /* Finding eigenvalue for nonsymetric matrix using QR algorithm */
X hessenberg (Q, Ng, x, iy);
X /* simplesrt(Ng,x); */
X /* Returns the sqrt of the second largest eigenvalue of Q */
X for (i = 2, tmp = x[1]; i <= Ng; ++i)
X tmp = (tmp > x[i]) ? tmp : x[i];
X return sqrt (x[Ng - 1]);
}
X
float *vector (nl, nh)
X int nl, nh;
{
X float *v;
X
X v = (float *) malloc ((unsigned) (nh - nl + 1) * sizeof (float));
X if (!v)
X fprintf (stderr, "memory allocation failure"), exit (1);
X return v - nl;
}
X
X
float **matrix (nrl, nrh, ncl, nch)
X int nrl, nrh, ncl, nch;
X
/* Allocates a float matrix with range [nrl..nrh][ncl..nch] */
{
X int i;
X float **m;
X
X /* allocate pointers to rows */
X m = (float **) malloc ((unsigned) (nrh - nrl + 1) * sizeof (float *));
X if (!m)
X fprintf (stderr, "memory allocation failure"), exit (1);
X m -= ncl;
X
X /* allocate rows and set pointers to them */
X for (i = nrl; i <= nrh; i++)
X {
X m[i] = (float *) malloc ((unsigned) (nch - ncl + 1) * sizeof (float));
X if (!m[i])
X fprintf (stderr, "memory allocation failure"), exit (2);
X m[i] -= ncl;
X }
X /* return pointer to array of pointers to rows */
X return m;
}
X
void results (c, a)
X char *c;
X float *a;
{
X int i;
X
X DOT;
X fprintf (stdout, "%s", c);
X for (i = 0; i < 4; ++i)
X fprintf (stdout, "% 1.3e ", a[i]);
X fprintf (stdout, "% 1.3e\n", (a[0] + a[1] + a[2] + a[3]) / 4);
}
X
void simplesrt (n, arr)
X int n;
X float arr[];
{
X int i, j;
X float a;
X
X for (j = 2; j <= n; j++)
X {
X a = arr[j];
X i = j - 1;
X while (i > 0 && arr[i] > a)
X {
X arr[i + 1] = arr[i];
X i--;
X }
X arr[i + 1] = a;
X }
}
X
void mkbalanced (a, n)
X float **a;
X int n;
{
X int last, j, i;
X float s, r, g, f, c, sqrdx;
X
X sqrdx = RADIX * RADIX;
X last = 0;
X while (last == 0)
X {
X last = 1;
X for (i = 1; i <= n; i++)
X {
X r = c = 0.0;
X for (j = 1; j <= n; j++)
X if (j != i)
X {
X c += fabs (a[j][i]);
X r += fabs (a[i][j]);
X }
X if (c && r)
X {
X g = r / RADIX;
X f = 1.0;
X s = c + r;
X while (c < g)
X {
X f *= RADIX;
X c *= sqrdx;
X }
X g = r * RADIX;
X while (c > g)
X {
X f /= RADIX;
X c /= sqrdx;
X }
X if ((c + r) / f < 0.95 * s)
X {
X last = 0;
X g = 1.0 / f;
X for (j = 1; j <= n; j++)
X a[i][j] *= g;
X for (j = 1; j <= n; j++)
X a[j][i] *= f;
X }
X }
X }
X }
}
X
X
void reduction (a, n)
X float **a;
X int n;
{
X int m, j, i;
X float y, x;
X
X for (m = 2; m < n; m++)
X {
X x = 0.0;
X i = m;
X for (j = m; j <= n; j++)
X {
X if (fabs (a[j][m - 1]) > fabs (x))
X {
X x = a[j][m - 1];
X i = j;
X }
X }
X if (i != m)
X {
X for (j = m - 1; j <= n; j++)
X SWAP (a[i][j], a[m][j])
X for (j = 1; j <= n; j++)
X SWAP (a[j][i], a[j][m])
X a[j][i] = a[j][i];
X }
X if (x)
X {
X for (i = m + 1; i <= n; i++)
X {
X if (y = a[i][m - 1])
X {
X y /= x;
X a[i][m - 1] = y;
X for (j = m; j <= n; j++)
X a[i][j] -= y * a[m][j];
X for (j = 1; j <= n; j++)
X a[j][m] += y * a[j][i];
X }
X }
X }
X }
}
X
void hessenberg (a, n, wr, wi)
X float **a, wr[], wi[];
X int n;
X
{
X int nn, m, l, k, j, its, i, mmin;
X float z, y, x, w, v, u, t, s, r, q, p, anorm;
X
X anorm = fabs (a[1][1]);
X for (i = 2; i <= n; i++)
X for (j = (i - 1); j <= n; j++)
X anorm += fabs (a[i][j]);
X nn = n;
X t = 0.0;
X while (nn >= 1)
X {
X its = 0;
X do
X {
X for (l = nn; l >= 2; l--)
X {
X s = fabs (a[l - 1][l - 1]) + fabs (a[l][l]);
X if (s == 0.0)
X s = anorm;
X if ((float) (fabs (a[l][l - 1]) + s) == s)
X break;
X }
X x = a[nn][nn];
X if (l == nn)
X {
X wr[nn] = x + t;
X wi[nn--] = 0.0;
X }
X else
X {
X y = a[nn - 1][nn - 1];
X w = a[nn][nn - 1] * a[nn - 1][nn];
X if (l == (nn - 1))
X {
X p = 0.5 * (y - x);
X q = p * p + w;
X z = sqrt (fabs (q));
X x += t;
X if (q >= 0.0)
X {
X z = p + SIGN (z, p);
X wr[nn - 1] = wr[nn] = x + z;
X if (z)
X wr[nn] = x - w / z;
X wi[nn - 1] = wi[nn] = 0.0;
X }
X else
X {
X wr[nn - 1] = wr[nn] = x + p;
X wi[nn - 1] = -(wi[nn] = z);
X }
X nn -= 2;
X }
X else
X {
X if (its == 30)
X fprintf (stderr,
X "Too many iterations to required to find %s\ngiving up\n",
X F14), exit (1);
X if (its == 10 || its == 20)
X {
X t += x;
X for (i = 1; i <= nn; i++)
X a[i][i] -= x;
X s = fabs (a[nn][nn - 1]) + fabs (a[nn - 1][nn - 2]);
X y = x = 0.75 * s;
X w = -0.4375 * s * s;
X }
X ++its;
X for (m = (nn - 2); m >= l; m--)
X {
X z = a[m][m];
X r = x - z;
X s = y - z;
X p = (r * s - w) / a[m + 1][m] + a[m][m + 1];
X q = a[m + 1][m + 1] - z - r - s;
X r = a[m + 2][m + 1];
X s = fabs (p) + fabs (q) + fabs (r);
X p /= s;
X q /= s;
X r /= s;
X if (m == l)
X break;
X u = fabs (a[m][m - 1]) * (fabs (q) + fabs (r));
X v = fabs (p) * (fabs (a[m - 1][m - 1]) + fabs (z) + fabs (a[m + 1][m + 1]));
X if ((float) (u + v) == v)
X break;
X }
X for (i = m + 2; i <= nn; i++)
X {
X a[i][i - 2] = 0.0;
X if (i != (m + 2))
X a[i][i - 3] = 0.0;
X }
X for (k = m; k <= nn - 1; k++)
X {
X if (k != m)
X {
X p = a[k][k - 1];
X q = a[k + 1][k - 1];
X r = 0.0;
X if (k != (nn - 1))
X r = a[k + 2][k - 1];
X if (x = fabs (p) + fabs (q) + fabs (r))
X {
X p /= x;
X q /= x;
X r /= x;
X }
X }
X if (s = SIGN (sqrt (p * p + q * q + r * r), p))
X {
X if (k == m)
X {
X if (l != m)
X a[k][k - 1] = -a[k][k - 1];
X }
X else
X a[k][k - 1] = -s * x;
X p += s;
X x = p / s;
X y = q / s;
X z = r / s;
X q /= p;
X r /= p;
X for (j = k; j <= nn; j++)
X {
X p = a[k][j] + q * a[k + 1][j];
X if (k != (nn - 1))
X {
X p += r * a[k + 2][j];
X a[k + 2][j] -= p * z;
X }
X a[k + 1][j] -= p * y;
X a[k][j] -= p * x;
X }
X mmin = nn < k + 3 ? nn : k + 3;
X for (i = l; i <= mmin; i++)
X {
X p = x * a[i][k] + y * a[i][k + 1];
X if (k != (nn - 1))
X {
X p += z * a[i][k + 2];
X a[i][k + 2] -= p * r;
X }
X a[i][k + 1] -= p * q;
X a[i][k] -= p;
X }
X }
X }
X }
X }
X } while (l < nn - 1);
X }
}
SHAR_EOF
chmod 0644 pgmtxtur.c ||
echo 'restore of pgmtxtur.c failed'
Wc_c="`wc -c < 'pgmtxtur.c'`"
test 24184 -eq "$Wc_c" ||
echo 'pgmtxtur.c: original size 24184, current size' "$Wc_c"
fi
exit 0
---end of shar archive