From 6e0b695e12a04a504f69da4eaba434e75ee95db8 Mon Sep 17 00:00:00 2001 From: giraffedata Date: Sun, 10 Jul 2011 03:06:42 +0000 Subject: cleanup git-svn-id: http://svn.code.sf.net/p/netpbm/code/trunk@1510 9d0c8265-081b-0410-96cb-a4ca84ce46f8 --- analyzer/pgmtexture.c | 1317 ++++++++++++++++++++++++++----------------------- 1 file changed, 713 insertions(+), 604 deletions(-) (limited to 'analyzer/pgmtexture.c') diff --git a/analyzer/pgmtexture.c b/analyzer/pgmtexture.c index 1ad89264..76f66c6c 100644 --- a/analyzer/pgmtexture.c +++ b/analyzer/pgmtexture.c @@ -1,4 +1,4 @@ -/* pgmtexture.c - calculate textural features on a portable graymap +/* pgmtexture.c - calculate textural features of a PGM image ** ** Author: James Darrell McCauley ** Texas Agricultural Experiment Station @@ -48,11 +48,12 @@ */ +#include #include #include "pm_c_util.h" -#include "pgm.h" #include "mallocvar.h" +#include "pgm.h" #define RADIX 2.0 #define EPSILON 0.000000001 @@ -72,32 +73,60 @@ #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;} +#define SWAP(a,b) do {float const y=(a);(a)=(b);(b)=y;} while (0) + + + +static float +sign(float const x, + float const y) { + + return y < 0 ? -fabs(x) : fabs(x); +} + + + +static bool const sortit = FALSE; -static bool sortit = FALSE; static float * -vector (int nl, int nh) -{ - float *v; +vector(unsigned int const nl, + unsigned int const nh) { + + float * v; + + assert(nh >= nl); 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) +matrix (unsigned int const nrl, + unsigned int const nrh, + unsigned int const ncl, + unsigned int const nch) { +/*---------------------------------------------------------------------------- + Allocate a float matrix with range [nrl..nrh][ncl..nch] -/* Allocates a float matrix with range [nrl..nrh][ncl..nch] */ -{ - int i; - float **m; + We do some seedy C here, subtracting an arbitrary integer from a pointer and + calling the result a pointer. It normally works because the only way we'll + use that pointer is by adding that same integer or something greater to it. + + The point of this is not to allocate memory for matrix elements that will + never be referenced (row < nrl or column < ncl). +-----------------------------------------------------------------------------*/ + unsigned int i; + float ** m; + + assert(nrh >= nrl); /* allocate pointers to rows */ MALLOCARRAY(m, (unsigned) (nrh - nrl + 1)); @@ -106,92 +135,102 @@ matrix (int nrl, int nrh, int ncl, int nch) m -= ncl; + assert (nch >= ncl); + /* allocate rows and set pointers to them */ - for (i = nrl; i <= nrh; i++) - { + 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; +results (const char * const name, + const float * const a) { + + unsigned int i; + + fprintf (stdout, "%s", name); - 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; +simplesrt (unsigned int const n, + float * const arr) { + + unsigned int j; float a; - for (j = 2; j <= n; j++) - { + for (j = 2; j <= n; ++j) { + unsigned int i; + a = arr[j]; - i = j - 1; - while (i > 0 && arr[i] > a) - { - arr[i + 1] = arr[i]; - i--; + i = j; + + while (i > 1 && arr[i-1] > a) { + arr[i] = arr[i-1]; + --i; } - arr[i + 1] = a; + arr[i] = a; } } + + static void -mkbalanced (float **a, int n) -{ - int last, j, i; - float s, r, g, f, c, sqrdx; +mkbalanced (float ** const a, + unsigned int const n) { + + float const sqrdx = SQR(RADIX); + + unsigned int last, i; + float s, r, g, f, c; - sqrdx = RADIX * RADIX; last = 0; - while (last == 0) - { + while (last == 0) { last = 1; - for (i = 1; i <= n; i++) - { + for (i = 1; i <= n; ++i) { + unsigned int j; r = c = 0.0; - for (j = 1; j <= n; j++) - if (j != i) - { + for (j = 1; j <= n; ++j) { + if (j != i) { c += fabs (a[j][i]); r += fabs (a[i][j]); } - if (c && r) - { + } + if (c && r) { g = r / RADIX; f = 1.0; s = c + r; - while (c < g) - { + while (c < g) { f *= RADIX; c *= sqrdx; } g = r * RADIX; - while (c > g) - { + while (c > g) { f /= RADIX; c /= sqrdx; } - if ((c + r) / f < 0.95 * s) - { + if ((c + r) / f < 0.95 * s) { + unsigned int j; last = 0; g = 1.0 / f; - for (j = 1; j <= n; j++) + for (j = 1; j <= n; ++j) a[i][j] *= g; - for (j = 1; j <= n; j++) + for (j = 1; j <= n; ++j) a[j][i] *= f; } } @@ -200,43 +239,43 @@ mkbalanced (float **a, int n) } + static void -reduction (float **a, int n) -{ - int m, j, i; - float y, x; +reduction (float ** const a, + unsigned int const n) { - for (m = 2; m < n; m++) - { + unsigned int m; + + for (m = 2; m < n; ++m) { + unsigned int j; + unsigned int i; + float x; x = 0.0; i = m; - for (j = m; j <= n; j++) - { - if (fabs (a[j][m - 1]) > fabs (x)) - { + 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 (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])) - { + if (x != 0.0) { + unsigned int i; + for (i = m + 1; i <= n; ++i) { + float y; + y = a[i][m - 1]; + if (y) { y /= x; a[i][m - 1] = y; - for (j = m; j <= n; j++) + for (j = m; j <= n; ++j) a[i][j] -= y * a[m][j]; - for (j = 1; j <= n; j++) + for (j = 1; j <= n; ++j) a[j][m] += y * a[j][i]; } } @@ -246,128 +285,140 @@ reduction (float **a, int n) +static float +norm(float ** const a, + unsigned int const n) { + + float anorm; + unsigned int i; + + for (i = 2, anorm = fabs(a[1][1]); i <= n; ++i) { + unsigned int j; + for (j = (i - 1); j <= n; ++j) + anorm += fabs(a[i][j]); + } + return anorm; +} + + + 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) - { +hessenberg(float ** const a, + unsigned int const n, + float * const wr, + float * const wi) { + + float const anorm = norm(a, n); + + int nn; + float t; + + assert(n >= 1); + + for (nn = n, t = 0.0; nn >= 1; ) { + unsigned int its; + int l; its = 0; - do - { - for (l = nn; l >= 2; l--) - { + do { + float x; + for (l = nn; l >= 2; --l) { + float s; 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; } + assert(nn >= 1); x = a[nn][nn]; - if (l == 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)); + } else { + float w, y; + y = a[nn - 1][nn - 1]; /* initial value */ + w = a[nn][nn - 1] * a[nn - 1][nn]; /* initial value */ + if (l == (nn - 1)) { + float const p = 0.5 * (y - x); + float const q = p * p + w; + float const 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; + if (q >= 0.0) { + float const z2 = p + sign(z, p); + wr[nn - 1] = wr[nn] = x + z2; + if (z2) + wr[nn] = x - w / z2; wi[nn - 1] = wi[nn] = 0.0; - } - else - { + } else { wr[nn - 1] = wr[nn] = x + p; wi[nn - 1] = -(wi[nn] = z); } nn -= 2; - } - else - { + } else { + int i, k, m; + float p, q, r; if (its == 30) pm_error("Too many iterations to required " "to find %s. Giving up", F14); - if (its == 10 || its == 20) - { + if (its == 10 || its == 20) { + int i; + float s; t += x; - for (i = 1; i <= nn; i++) + for (i = 1; i <= nn; ++i) a[i][i] -= x; - s = fabs (a[nn][nn - 1]) + fabs (a[nn - 1][nn - 2]); + 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]; + for (m = (nn - 2); m >= l; --m) { + float const z = a[m][m]; + float s, u, v; 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); + 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) + 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 (u + v == v) break; } - for (i = m + 2; i <= nn; i++) - { + 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) - { + for (k = m; k <= nn - 1; ++k) { + float s; + 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))) - { + 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) - { + s = sign(sqrt(SQR(p) + SQR(q) + SQR(r)), p); + if (s) { + int const mmin = nn < k + 3 ? nn : k + 3; + float z; + int j; + if (k == m) { if (l != m) a[k][k - 1] = -a[k][k - 1]; - } - else + } else a[k][k - 1] = -s * x; p += s; x = p / s; @@ -375,23 +426,18 @@ hessenberg (float **a, int n, float wr[], float wi[]) z = r / s; q /= p; r /= p; - for (j = k; j <= nn; j++) - { + for (j = k; j <= nn; ++j) { p = a[k][j] + q * a[k + 1][j]; - if (k != (nn - 1)) - { + 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++) - { + for (i = l; i <= mmin; ++i) { p = x * a[i][k] + y * a[i][k + 1]; - if (k != (nn - 1)) - { + if (k != (nn - 1)) { p += z * a[i][k + 2]; a[i][k + 2] -= p * r; } @@ -409,427 +455,485 @@ hessenberg (float **a, int n, float wr[], float wi[]) 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]; - +f1_a2m(float ** const p, + unsigned int const ng) { +/*---------------------------------------------------------------------------- + Angular Second Moment + + 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. +-----------------------------------------------------------------------------*/ + unsigned int i; + float sum; + + for (i = 0, sum = 0.0; i < ng; ++i) { + unsigned int j; + 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) +static float +f2_contrast(float ** const p, + unsigned int const ng) { +/*---------------------------------------------------------------------------- + Contrast + + 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. +-----------------------------------------------------------------------------*/ + unsigned int n; + float bigsum; + + for (n = 0, bigsum = 0.0; n < ng; ++n) { + unsigned int i; + float sum; + for (i = 0, sum = 0.0; i < ng; ++i) { + unsigned int j; + for (j = 0; j < ng; ++j) { if ((i - j) == n || (j - i) == n) - sum += P[i][j]; - bigsum += n * n * sum; - - sum = 0; + sum += p[i][j]; + } + } + bigsum += SQR(n) * sum; } 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) +static float +f3_corr(float ** const p, + unsigned int const ng) { +/*---------------------------------------------------------------------------- + Correlation + + This correlation feature is a measure of gray-tone linear-dependencies in + the image. +-----------------------------------------------------------------------------*/ + unsigned int i; + float sumSqrx, sumSqry, tmp; + float * px; + float meanx, meany, stddevx, stddevy; + + sumSqrx = 0.0; sumSqry = 0.0; + meanx = 0.0; meany = 0.0; + + 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]; - + /* 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) { + unsigned int j; + 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; + for (i = 0; i < ng; ++i) { + meanx += px[i] * i; + sumSqrx += px[i] * SQR(i); } + meany = meanx; - sum_sqry = sum_sqrx; - stddevx = sqrt (sum_sqrx - (meanx * meanx)); + sumSqry = sumSqrx; + stddevx = sqrt(sumSqrx - (SQR(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]; - + for (i = 0, tmp = 0; i < ng; ++i) { + unsigned int j; + 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]; +static float +f4_var (float ** const p, + unsigned int const ng) { +/*---------------------------------------------------------------------------- + Sum of Squares: Variance +-----------------------------------------------------------------------------*/ + unsigned int i; + float mean, var; + + for (i = 0, mean = 0.0; i < ng; ++i) { + unsigned int j; + for (j = 0; j < ng; ++j) + mean += i * p[i][j]; + } + for (i = 0, var = 0.0; i < ng; ++i) { + unsigned int j; + 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)); +static float +f5_idm (float ** const p, + unsigned int const ng) { +/*---------------------------------------------------------------------------- + Inverse Difference Moment +-----------------------------------------------------------------------------*/ + unsigned int i; + float idm; + + for (i = 0, idm = 0.0; i < ng; ++i) { + unsigned int j; + 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]; +static float +f6_savg (float ** const p, + unsigned int const ng) { +/*---------------------------------------------------------------------------- + Sum Average +-----------------------------------------------------------------------------*/ + float pxpy[2 * PGM_MAXMAXVAL]; + unsigned int i; + float savg; + + for (i = 0; i <= 2 * ng; ++i) + pxpy[i] = 0.0; + + for (i = 0; i < ng; ++i) { + unsigned int j; + for (j = 0; j < ng; ++j) + pxpy[i + j + 2] += p[i][j]; + } + for (i = 2, savg = 0.0; 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]; +static float +f7_svar (float ** const p, + unsigned int const ng, + float const s) { +/*---------------------------------------------------------------------------- + Sum Variance +-----------------------------------------------------------------------------*/ + float pxpy[2 * PGM_MAXMAXVAL]; + unsigned int i; + float var; + + for (i = 0; i <= 2 * ng; ++i) + pxpy[i] = 0; + + for (i = 0; i < ng; ++i) { + unsigned int j; + for (j = 0; j < ng; ++j) + pxpy[i + j + 2] += p[i][j]; + } + for (i = 2, var = 0.0; 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); +static float +f8_sentropy (float ** const p, + unsigned int const ng) { +/*---------------------------------------------------------------------------- + Sum Entropy +-----------------------------------------------------------------------------*/ + float pxpy[2 * PGM_MAXMAXVAL]; + unsigned int i; + float sentropy; + + for (i = 0; i <= 2 * ng; ++i) + pxpy[i] = 0; + + for (i = 0; i < ng; ++i) { + unsigned int j; + for (j = 0; j < ng; ++j) + pxpy[i + j + 2] += p[i][j]; + } + for (i = 2, sentropy = 0.0; 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); +static float +f9_entropy (float ** const p, + unsigned int const ng) { +/*---------------------------------------------------------------------------- + Entropy +-----------------------------------------------------------------------------*/ + unsigned int i; + float entropy; + + for (i = 0, entropy = 0.0; i < ng; ++i) { + unsigned int j; + 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]; +static float +f10_dvar (float ** const p, + unsigned int const ng) { +/*---------------------------------------------------------------------------- + Difference Variance +-----------------------------------------------------------------------------*/ + float pxpy[2 * PGM_MAXMAXVAL]; + unsigned int i; + unsigned int sqrNg; + float sum, sumSqr, var; + + for (i = 0; i <= 2 * ng; ++i) + pxpy[i] = 0; + + for (i = 0; i < ng; ++i) { + unsigned int j; + 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]; + for (i = 0, sum = 0.0, sumSqr = 0.0; i < ng; ++i) { + sum += pxpy[i]; + sumSqr += pxpy[i] * pxpy[i]; } - tmp = Ng * Ng; - var = ((tmp * sum_sqr) - (sum * sum)) / (tmp * tmp); + sqrNg = SQR(ng); + var = (sqrNg * sumSqr - SQR(sum)) / SQR(sqrNg); 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); +static float +f11_dentropy (float ** const p, + unsigned int const ng) { +/*---------------------------------------------------------------------------- + Difference Entropy +-----------------------------------------------------------------------------*/ + float pxpy[2 * PGM_MAXMAXVAL]; + unsigned int i; + float sum; + + for (i = 0; i <= 2 * ng; ++i) + pxpy[i] = 0; + + for (i = 0; i < ng; ++i) { + unsigned int j; + for (j = 0; j < ng; ++j) + pxpy[abs(i - j)] += p[i][j]; + } + for (i = 0, sum = 0.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); +static float +f12_icorr (float ** const p, + unsigned int const ng) { +/*---------------------------------------------------------------------------- + Information Measures of Correlation +-----------------------------------------------------------------------------*/ + unsigned int i; + float * px; + float * py; + float hx, hy, hxy, hxy1, hxy2; + + 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) { + unsigned int j; + 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); - } + hx = 0.0; hy = 0.0; hxy = 0.0; hxy1 = 0.0; hxy2 = 0.0; + for (i = 0; i < ng; ++i) { + unsigned int j; + 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); + 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)); + 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); +static float +f13_icorr (float ** const p, + unsigned int const ng) { +/*---------------------------------------------------------------------------- + Information Measures of Correlation +-----------------------------------------------------------------------------*/ + unsigned int i; + float * px; + float * py; + float hx, hy, hxy, hxy1, hxy2; + + 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) { + unsigned int j; + 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); - } + hx = 0.0; hy = 0.0; hxy = 0.0; hxy1 = 0.0; hxy2 = 0.0; + for (i = 0; i < ng; ++i) { + unsigned int j; + 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) - { + 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))))); + 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); +static float +f14_maxcorr (float ** const p, + unsigned int const ng) { +/*---------------------------------------------------------------------------- + The Maximal Correlation Coefficient +-----------------------------------------------------------------------------*/ + unsigned int i; + float *px, *py; + float ** q; + float * x; + float * iy; + float 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]; + for (i = 0; i < ng; ++i) { + unsigned int j; + 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]; + /* Compute the Q matrix */ + for (i = 0; i < ng; ++i) { + unsigned int j; + for (j = 0; j < ng; ++j) { + unsigned int k; + 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); + mkbalanced(q, ng); /* Reduction to Hessenberg Form */ - reduction (Q, Ng); + reduction(q, ng); /* Finding eigenvalue for nonsymetric matrix using QR algorithm */ - hessenberg (Q, Ng, x, iy); + 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) + simplesrt(ng, x); + + /* Return 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]); + + 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]; +main (int argc, const char ** argv) { + + FILE * ifP; + gray ** grays; + unsigned int tone[PGM_MAXMAXVAL]; + unsigned int r0, r45, r90; + unsigned int d; + unsigned int x, y; + unsigned int row; + int rows, cols; + int argn; + unsigned int itone; + unsigned int toneCt; + float ** p_matrix0, ** p_matrix45, ** p_matrix90, ** p_matrix135; + float a2m[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; + unsigned int i; const char * const usage = "[-d ] [pgmfile]"; - - pgm_init( &argc, argv ); + pm_proginit(&argc, argv); argn = 1; @@ -839,7 +943,7 @@ main (int argc, char *argv[]) { if ( argv[argn][1] == 'd' ) { ++argn; - if ( argn == argc || sscanf( argv[argn], "%d", &d ) != 1 ) + if ( argn == argc || sscanf( argv[argn], "%u", &d ) != 1 ) pm_usage( usage ); } else @@ -849,205 +953,210 @@ main (int argc, char *argv[]) { if ( argn < argc ) { - ifp = pm_openr( argv[argn] ); + ifP = pm_openr( argv[argn] ); ++argn; } else - ifp = stdin; + ifP = stdin; if ( argn != argc ) pm_usage( usage ); - grays = pgm_readpgm (ifp, &cols, &rows, &maxval); - pm_close (ifp); + d = 1; - 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); + grays = pgm_readpgm (ifP, &cols, &rows, &maxval); + pm_close (ifP); /* 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 (i = 0; i <= PGM_MAXMAXVAL; ++i) + tone[i] = -1; + for (row = 0; row < rows; ++row) { + unsigned int col; 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); + } + for (i = 0, toneCt = 0; i <= PGM_MAXMAXVAL; ++i) { + if (tone[i] != -1) + ++toneCt; + } + pm_message("(Image has %u graylevels.)", toneCt); /* Collapse array, taking out all zero values */ - for (row = 0, itone = 0; row <= PGM_MAXMAXVAL; row++) + 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; + p_matrix0 = matrix (0, toneCt, 0, toneCt); + p_matrix45 = matrix (0, toneCt, 0, toneCt); + p_matrix90 = matrix (0, toneCt, 0, toneCt); + p_matrix135 = matrix (0, toneCt, 0, toneCt); + + for (row = 0; row < toneCt; ++row) { + unsigned int col; + for (col = 0; col < toneCt; ++col) { + p_matrix0 [row][col] = p_matrix45 [row][col] = 0; + p_matrix90[row][col] = p_matrix135[row][col] = 0; } + } + if (d > cols) + pm_error("Image is narrower (%u columns) " + "than specified distance (%u)", cols, d); /* 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) - { + pm_message("Computing spatial dependence matrix..."); + for (row = 0; row < rows; ++row) { + unsigned int col; + for (col = 0; col < cols; ++col) { + unsigned int angle; + for (angle = 0, x = 0; angle <= 135; angle += 45) { while (tone[x] != grays[row][col]) - x++; - if (angle == 0 && col + d < cols) - { + ++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]++; + ++y; + ++p_matrix0[x][y]; + ++p_matrix0[y][x]; } - if (angle == 90 && row + d < rows) - { + if (angle == 90 && row + d < rows) { y = 0; while (tone[y] != grays[row + d][col]) - y++; - P_matrix90[x][y]++; - P_matrix90[y][x]++; + ++y; + ++p_matrix90[x][y]; + ++p_matrix90[y][x]; } - if (angle == 45 && row + d < rows && col - d >= 0) - { + if (angle == 45 && row + d < rows && col >= d) { y = 0; while (tone[y] != grays[row + d][col - d]) - y++; - P_matrix45[x][y]++; - P_matrix45[y][x]++; + ++y; + ++p_matrix45[x][y]; + ++p_matrix45[y][x]; } - if (angle == 135 && row + d < rows && col + d < cols) - { + 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]++; + ++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; + 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; + for (itone = 0; itone < toneCt; ++itone) { + unsigned int jtone; + for (jtone = 0; jtone < toneCt; ++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); - - svar[0] = f7_svar (P_matrix0, tones, savg[0]); - svar[1] = f7_svar (P_matrix45, tones, savg[1]); - svar[2] = f7_svar (P_matrix90, tones, savg[2]); - svar[3] = f7_svar (P_matrix135, tones, savg[3]); - results (F7, svar); - - 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); - 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); + } + pm_message(" ...done."); + + pm_message("Computing textural features ..."); + + fprintf(stdout, "\n"); + fprintf(stdout, + "%s 0 45 90 135 Avg\n", + BL); + + a2m[0] = f1_a2m(p_matrix0, toneCt); + a2m[1] = f1_a2m(p_matrix45, toneCt); + a2m[2] = f1_a2m(p_matrix90, toneCt); + a2m[3] = f1_a2m(p_matrix135, toneCt); + results(F1, a2m); + + contrast[0] = f2_contrast(p_matrix0, toneCt); + contrast[1] = f2_contrast(p_matrix45, toneCt); + contrast[2] = f2_contrast(p_matrix90, toneCt); + contrast[3] = f2_contrast(p_matrix135, toneCt); + results(F2, contrast); + + + corr[0] = f3_corr(p_matrix0, toneCt); + corr[1] = f3_corr(p_matrix45, toneCt); + corr[2] = f3_corr(p_matrix90, toneCt); + corr[3] = f3_corr(p_matrix135, toneCt); + results(F3, corr); + + var[0] = f4_var(p_matrix0, toneCt); + var[1] = f4_var(p_matrix45, toneCt); + var[2] = f4_var(p_matrix90, toneCt); + var[3] = f4_var(p_matrix135, toneCt); + results(F4, var); + + + idm[0] = f5_idm(p_matrix0, toneCt); + idm[1] = f5_idm(p_matrix45, toneCt); + idm[2] = f5_idm(p_matrix90, toneCt); + idm[3] = f5_idm(p_matrix135, toneCt); + results(F5, idm); + + savg[0] = f6_savg(p_matrix0, toneCt); + savg[1] = f6_savg(p_matrix45, toneCt); + savg[2] = f6_savg(p_matrix90, toneCt); + savg[3] = f6_savg(p_matrix135, toneCt); + results(F6, savg); + + svar[0] = f7_svar(p_matrix0, toneCt, savg[0]); + svar[1] = f7_svar(p_matrix45, toneCt, savg[1]); + svar[2] = f7_svar(p_matrix90, toneCt, savg[2]); + svar[3] = f7_svar(p_matrix135, toneCt, savg[3]); + results(F7, svar); + + sentropy[0] = f8_sentropy(p_matrix0, toneCt); + sentropy[1] = f8_sentropy(p_matrix45, toneCt); + sentropy[2] = f8_sentropy(p_matrix90, toneCt); + sentropy[3] = f8_sentropy(p_matrix135, toneCt); + results(F8, sentropy); + + entropy[0] = f9_entropy(p_matrix0, toneCt); + entropy[1] = f9_entropy(p_matrix45, toneCt); + entropy[2] = f9_entropy(p_matrix90, toneCt); + entropy[3] = f9_entropy(p_matrix135, toneCt); + results(F9, entropy); + + dvar[0] = f10_dvar(p_matrix0, toneCt); + dvar[1] = f10_dvar(p_matrix45, toneCt); + dvar[2] = f10_dvar(p_matrix90, toneCt); + dvar[3] = f10_dvar(p_matrix135, toneCt); + results(F10, dvar); + + dentropy[0] = f11_dentropy(p_matrix0, toneCt); + dentropy[1] = f11_dentropy(p_matrix45, toneCt); + dentropy[2] = f11_dentropy(p_matrix90, toneCt); + dentropy[3] = f11_dentropy(p_matrix135, toneCt); 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"); + icorr[0] = f12_icorr(p_matrix0, toneCt); + icorr[1] = f12_icorr(p_matrix45, toneCt); + icorr[2] = f12_icorr(p_matrix90, toneCt); + icorr[3] = f12_icorr(p_matrix135, toneCt); + results(F12, icorr); + + icorr[0] = f13_icorr(p_matrix0, toneCt); + icorr[1] = f13_icorr(p_matrix45, toneCt); + icorr[2] = f13_icorr(p_matrix90, toneCt); + icorr[3] = f13_icorr(p_matrix135, toneCt); + results(F13, icorr); + + maxcorr[0] = f14_maxcorr(p_matrix0, toneCt); + maxcorr[1] = f14_maxcorr(p_matrix45, toneCt); + maxcorr[2] = f14_maxcorr(p_matrix90, toneCt); + maxcorr[3] = f14_maxcorr(p_matrix135, toneCt); + results(F14, maxcorr); + + pm_message(" ...done."); return 0; } -- cgit 1.4.1