From 1fd361a1ea06e44286c213ca1f814f49306fdc43 Mon Sep 17 00:00:00 2001 From: giraffedata Date: Sat, 19 Aug 2006 03:12:28 +0000 Subject: Create Subversion repository git-svn-id: http://svn.code.sf.net/p/netpbm/code/trunk@1 9d0c8265-081b-0410-96cb-a4ca84ce46f8 --- analyzer/pgmtexture.c | 1047 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1047 insertions(+) create mode 100644 analyzer/pgmtexture.c (limited to 'analyzer/pgmtexture.c') diff --git a/analyzer/pgmtexture.c b/analyzer/pgmtexture.c new file mode 100644 index 00000000..38eab114 --- /dev/null +++ b/analyzer/pgmtexture.c @@ -0,0 +1,1047 @@ +/* pgmtexture.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. +** +** 05 Oct 05 - Marc Breithecker +** Fix calculation or normalizing constants for d > 1. +*/ + +#include + +#include "pm_c_util.h" +#include "pgm.h" +#include "mallocvar.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;} + + +static bool sortit = FALSE; + +static float * +vector (int nl, int nh) +{ + float *v; + + MALLOCARRAY(v, (unsigned) (nh - nl + 1)); + if (v == NULL) + pm_error("Unable to allocate memory for a vector."); + return v - nl; +} + + +static float ** +matrix (int nrl, int nrh, int ncl, int nch) + +/* Allocates a float matrix with range [nrl..nrh][ncl..nch] */ +{ + int i; + float **m; + + /* allocate pointers to rows */ + MALLOCARRAY(m, (unsigned) (nrh - nrl + 1)); + if (m == NULL) + pm_error("Unable to allocate memory for a matrix."); + + m -= ncl; + + /* allocate rows and set pointers to them */ + for (i = nrl; i <= nrh; i++) + { + MALLOCARRAY(m[i], (unsigned) (nch - ncl + 1)); + if (m[i] == NULL) + pm_error("Unable to allocate memory for a matrix row."); + m[i] -= ncl; + } + /* return pointer to array of pointers to rows */ + return m; +} + +static void +results (const char * const c, const float * const a) +{ + int i; + + DOT; + fprintf (stdout, "%s", c); + for (i = 0; i < 4; ++i) + fprintf (stdout, "% 1.3e ", a[i]); + fprintf (stdout, "% 1.3e\n", (a[0] + a[1] + a[2] + a[3]) / 4); +} + +static void +simplesrt (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; + } +} + +static void +mkbalanced (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; + } + } + } + } +} + + +static void +reduction (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]; + } + } + } + } +} + + + +static void +hessenberg (float **a, int n, float wr[], float wi[]) + +{ + 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) + pm_error("Too many iterations to required " + "to find %s. Giving up", F14); + 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); + } +} + + + +static float +f1_asm (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. + */ +} + + +static float +f2_contrast (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. + */ +} + +static float +f3_corr (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 = 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; + } + 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]; + + return (tmp - meanx * meany) / (stddevx * stddevy); + /* + * This correlation feature is a measure of gray-tone linear-dependencies + * in the image. + */ +} + + +static float +f4_var (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) + var += (i + 1 - mean) * (i + 1 - mean) * P[i][j]; + + return var; +} + +static float +f5_idm (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; +} + +static float +Pxpy[2 * PGM_MAXMAXVAL]; + +static float +f6_savg (float **P, int Ng) + +/* Sum Average */ +{ + int i, j; + 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) + Pxpy[i + j + 2] += P[i][j]; + for (i = 2; i <= 2 * Ng; ++i) + savg += i * Pxpy[i]; + + return savg; +} + + +static float +f7_svar (float **P, int Ng, float S) { +/* Sum Variance */ + int i, j; + 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) + Pxpy[i + j + 2] += P[i][j]; + + for (i = 2; i <= 2 * Ng; ++i) + var += (i - S) * (i - S) * Pxpy[i]; + + return var; +} + +static float +f8_sentropy (float **P, int Ng) + +/* Sum Entropy */ +{ + int i, j; + 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) + sentropy -= Pxpy[i] * log10 (Pxpy[i] + EPSILON); + + return sentropy; +} + + +static float +f9_entropy (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); + + return -entropy; +} + + +static float +f10_dvar (float **P, int Ng) + +/* Difference Variance */ +{ + int i, j, tmp; + 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 += Pxpy[i]; + sum_sqr += Pxpy[i] * Pxpy[i]; + } + tmp = Ng * Ng; + var = ((tmp * sum_sqr) - (sum * sum)) / (tmp * tmp); + + return var; +} + +static float +f11_dentropy (float **P, int Ng) + +/* Difference Entropy */ +{ + int i, j; + float sum = 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); + + return -sum; +} + +static float +f12_icorr (float **P, int Ng) + +/* Information Measures of Correlation */ +{ + int i, j; + float *px, *py; + float hx = 0, hy = 0, hxy = 0, hxy1 = 0, hxy2 = 0; + + px = vector (0, Ng); + py = 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); + hxy2 -= px[i] * py[j] * log10 (px[i] * py[j] + EPSILON); + hxy -= P[i][j] * log10 (P[i][j] + EPSILON); + } + + /* Calculate entropies of px and py - is this right? */ + for (i = 0; i < Ng; ++i) + { + hx -= px[i] * log10 (px[i] + EPSILON); + hy -= py[i] * log10 (py[i] + EPSILON); + } +/* fprintf(stderr,"hxy1=%f\thxy=%f\thx=%f\thy=%f\n",hxy1,hxy,hx,hy); */ + return ((hxy - hxy1) / (hx > hy ? hx : hy)); +} + +static float +f13_icorr (float **P, int Ng) + +/* Information Measures of Correlation */ +{ + int i, j; + float *px, *py; + float hx = 0, hy = 0, hxy = 0, hxy1 = 0, hxy2 = 0; + + px = vector (0, Ng); + py = 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); + hxy2 -= px[i] * py[j] * log10 (px[i] * py[j] + EPSILON); + hxy -= P[i][j] * log10 (P[i][j] + EPSILON); + } + + /* Calculate entropies of px and py */ + for (i = 0; i < Ng; ++i) + { + hx -= px[i] * log10 (px[i] + EPSILON); + hy -= py[i] * log10 (py[i] + EPSILON); + } + /* fprintf(stderr,"hx=%f\thxy2=%f\n",hx,hxy2); */ + return (sqrt (fabs (1 - exp (-2.0 * (hxy2 - hxy))))); +} + +static float +f14_maxcorr (float **P, int Ng) + +/* Returns the Maximal Correlation Coefficient */ +{ + int i, j, k; + float *px, *py, **Q; + float *x, *iy, tmp; + + px = vector (0, Ng); + py = vector (0, Ng); + Q = matrix (1, Ng + 1, 1, Ng + 1); + x = vector (1, Ng); + iy = 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 */ + hessenberg (Q, Ng, x, iy); + if (sortit) + 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]; + return sqrt (x[Ng - 1]); +} + +int +main (int argc, char *argv[]) { + FILE *ifp; + register gray **grays; + int tone[PGM_MAXMAXVAL], R0, R45, R90, angle, d = 1, x, y; + int argn, rows, cols, row, col; + int itone, jtone, tones; + 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]; + gray maxval; + const char * const usage = "[-d ] [pgmfile]"; + + + pgm_init( &argc, argv ); + + argn = 1; + + /* Check for flags. */ + if ( argn < argc && argv[argn][0] == '-' ) + { + if ( argv[argn][1] == 'd' ) + { + ++argn; + if ( argn == argc || sscanf( argv[argn], "%d", &d ) != 1 ) + pm_usage( usage ); + } + else + pm_usage( usage ); + ++argn; + } + + if ( argn < argc ) + { + ifp = pm_openr( argv[argn] ); + ++argn; + } + else + ifp = stdin; + + if ( argn != argc ) + pm_usage( usage ); + + grays = pgm_readpgm (ifp, &cols, &rows, &maxval); + pm_close (ifp); + + if (maxval > PGM_MAXMAXVAL) + pm_error("The maxval of the image (%d) is too high. \n" + "This program's maximum is %d.", maxval, PGM_MAXMAXVAL); + + /* 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) + tone[grays[row][col]] = grays[row][col]; + for (row = PGM_MAXMAXVAL, tones = 0; row >= 0; --row) + if (tone[row] != -1) + tones++; + pm_message("(Image has %d graylevels.)", 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 = matrix (0, tones, 0, tones); + P_matrix45 = matrix (0, tones, 0, tones); + P_matrix90 = matrix (0, tones, 0, tones); + P_matrix135 = 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; + } + + /* Find gray-tone spatial dependence matrix */ + fprintf (stderr, "(Computing spatial dependence matrix..."); + for (row = 0; row < rows; ++row) + for (col = 0; col < cols; ++col) + for (x = 0, angle = 0; angle <= 135; angle += 45) + { + while (tone[x] != grays[row][col]) + x++; + if (angle == 0 && col + d < cols) + { + y = 0; + while (tone[y] != grays[row][col + d]) + y++; + P_matrix0[x][y]++; + P_matrix0[y][x]++; + } + if (angle == 90 && row + d < rows) + { + y = 0; + while (tone[y] != grays[row + d][col]) + y++; + P_matrix90[x][y]++; + P_matrix90[y][x]++; + } + if (angle == 45 && row + d < rows && col - d >= 0) + { + y = 0; + while (tone[y] != grays[row + d][col - d]) + y++; + P_matrix45[x][y]++; + P_matrix45[y][x]++; + } + if (angle == 135 && row + d < rows && col + d < cols) + { + y = 0; + while (tone[y] != grays[row + d][col + d]) + y++; + P_matrix135[x][y]++; + P_matrix135[y][x]++; + } + } + /* Gray-tone spatial dependence matrices are complete */ + + /* Find normalizing constants */ + R0 = 2 * rows * (cols - d); + R45 = 2 * (rows - d) * (cols - d); + R90 = 2 * (rows - d) * cols; + + /* 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] /= R45; + } + + fprintf (stderr, " done.)\n"); + fprintf (stderr, "(Computing textural features"); + fprintf (stdout, "\n"); + DOT; + fprintf (stdout, + "%s 0 45 90 135 Avg\n", + BL); + + 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); + results (F1, ASM); + + 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); + results (F2, contrast); + + + 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); + results (F3, corr); + + 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); + results (F4, var); + + + 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); + results (F5, idm); + + 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); + results (F6, savg); + + 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); + 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]); + results (F7, svar); + results (F8, sentropy); + + 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); + results (F9, entropy); + + 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); + results (F10, dvar); + + 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); + results (F11, dentropy); + + 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); + results (F12, icorr); + + 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); + results (F13, icorr); + + 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); + results (F14, maxcorr); + + + fprintf (stderr, " done.)\n"); + + return 0; +} -- cgit 1.4.1