w32tex
About: TeX Live provides a comprehensive TeX system including all the major TeX-related programs, macro packages, and fonts that are free software. Windows sources.
  Fossies Dox: w32tex-src.tar.xz  ("unofficial" and yet experimental doxygen-generated source code documentation)  

pgmtexture.c
Go to the documentation of this file.
1 /* pgmtxtur.c - calculate textural features on a portable graymap
2 **
3 ** Author: James Darrell McCauley
4 ** Texas Agricultural Experiment Station
5 ** Department of Agricultural Engineering
6 ** Texas A&M University
7 ** College Station, Texas 77843-2117 USA
8 **
9 ** Code written partially taken from pgmtofs.c in the PBMPLUS package
10 ** by Jef Poskanzer.
11 **
12 ** Algorithms for calculating features (and some explanatory comments) are
13 ** taken from:
14 **
15 ** Haralick, R.M., K. Shanmugam, and I. Dinstein. 1973. Textural features
16 ** for image classification. IEEE Transactions on Systems, Man, and
17 ** Cybertinetics, SMC-3(6):610-621.
18 **
19 ** Copyright (C) 1991 Texas Agricultural Experiment Station, employer for
20 ** hire of James Darrell McCauley
21 **
22 ** Permission to use, copy, modify, and distribute this software and its
23 ** documentation for any purpose and without fee is hereby granted, provided
24 ** that the above copyright notice appear in all copies and that both that
25 ** copyright notice and this permission notice appear in supporting
26 ** documentation. This software is provided "as is" without express or
27 ** implied warranty.
28 **
29 ** THE TEXAS AGRICULTURAL EXPERIMENT STATION (TAES) AND THE TEXAS A&M
30 ** UNIVERSITY SYSTEM (TAMUS) MAKE NO EXPRESS OR IMPLIED WARRANTIES
31 ** (INCLUDING BY WAY OF EXAMPLE, MERCHANTABILITY) WITH RESPECT TO ANY
32 ** ITEM, AND SHALL NOT BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL
33 ** OR CONSEQUENTAL DAMAGES ARISING OUT OF THE POSESSION OR USE OF
34 ** ANY SUCH ITEM. LICENSEE AND/OR USER AGREES TO INDEMNIFY AND HOLD
35 ** TAES AND TAMUS HARMLESS FROM ANY CLAIMS ARISING OUT OF THE USE OR
36 ** POSSESSION OF SUCH ITEMS.
37 **
38 ** Modification History:
39 ** 24 Jun 91 - J. Michael Carstensen <jmc@imsor.dth.dk> supplied fix for
40 ** correlation function.
41 */
42 
43 #include <math.h>
44 #include "pgm.h"
45 
46 #define RADIX 2.0
47 #define EPSILON 0.000000001
48 #define BL "Angle "
49 #define F1 "Angular Second Moment "
50 #define F2 "Contrast "
51 #define F3 "Correlation "
52 #define F4 "Variance "
53 #define F5 "Inverse Diff Moment "
54 #define F6 "Sum Average "
55 #define F7 "Sum Variance "
56 #define F8 "Sum Entropy "
57 #define F9 "Entropy "
58 #define F10 "Difference Variance "
59 #define F11 "Difference Entropy "
60 #define F12 "Meas of Correlation-1 "
61 #define F13 "Meas of Correlation-2 "
62 #define F14 "Max Correlation Coeff "
63 
64 #define SIGN(x,y) ((y)<0 ? -fabs(x) : fabs(x))
65 #define DOT fprintf(stderr,".")
66 #define SWAP(a,b) {y=(a);(a)=(b);(b)=y;}
67 
68 float f1_asm ARGS((float **P, int Ng));
69 float f2_contrast ARGS((float **P, int Ng));
70 float f3_corr ARGS((float **P, int Ng));
71 float f4_var ARGS((float **P, int Ng));
72 float f5_idm ARGS((float **P, int Ng));
73 float f6_savg ARGS((float **P, int Ng));
74 float f7_svar ARGS((float **P, int Ng, float S));
75 float f8_sentropy ARGS((float **P, int Ng));
76 float f9_entropy ARGS((float **P, int Ng));
77 float f10_dvar ARGS((float **P, int Ng));
78 float f11_dentropy ARGS((float **P, int Ng));
79 float f12_icorr ARGS((float **P, int Ng));
80 float f13_icorr ARGS((float **P, int Ng));
81 float f14_maxcorr ARGS((float **P, int Ng));
82 float *vector ARGS((int nl, int nh));
83 float **matrix ARGS((int nrl, int nrh, int ncl, int nch));
84 void results ARGS((char *c, float *a));
85 void simplesrt ARGS((int n, float arr[]));
86 void mkbalanced ARGS((float **a, int n));
87 void reduction ARGS((float **a, int n));
88 void hessenberg ARGS((float **a, int n, float wr[], float wi[]));
89 
90 int
92  int argc;
93  char *argv[];
94 {
95  FILE *ifp;
96  register gray **grays;
97  int tone[PGM_MAXMAXVAL], R0, R45, R90, angle, d = 1, x, y;
98  int argn, rows, cols, row, col;
99  int itone, jtone, tones;
100  float **P_matrix0, **P_matrix45, **P_matrix90, **P_matrix135;
101  float ASM[4], contrast[4], corr[4], var[4], idm[4], savg[4];
102  float sentropy[4], svar[4], entropy[4], dvar[4], dentropy[4];
103  float icorr[4], maxcorr[4];
104  gray maxval;
105  char *usage = "[-d <d>] [pgmfile]";
106 
107 
108  pgm_init( &argc, argv );
109 
110  argn = 1;
111 
112  /* Check for flags. */
113  if ( argn < argc && argv[argn][0] == '-' )
114  {
115  if ( argv[argn][1] == 'd' )
116  {
117  ++argn;
118  if ( argn == argc || sscanf( argv[argn], "%d", &d ) != 1 )
119  pm_usage( usage );
120  }
121  else
122  pm_usage( usage );
123  ++argn;
124  }
125 
126  if ( argn < argc )
127  {
128  ifp = pm_openr( argv[argn] );
129  ++argn;
130  }
131  else
132  ifp = stdin;
133 
134  if ( argn != argc )
135  pm_usage( usage );
136 
137  grays = pgm_readpgm (ifp, &cols, &rows, &maxval);
138  pm_close (ifp);
139 
140  if (maxval > PGM_MAXMAXVAL)
141  pm_error("The maxval of the image (%d) is too high. \n"
142  "This program's maximum is %d.", maxval, PGM_MAXMAXVAL);
143 
144  /* Determine the number of different gray scales (not maxval) */
145  for (row = PGM_MAXMAXVAL; row >= 0; --row)
146  tone[row] = -1;
147  for (row = rows - 1; row >= 0; --row)
148  for (col = 0; col < cols; ++col)
149  tone[grays[row][col]] = grays[row][col];
150  for (row = PGM_MAXMAXVAL, tones = 0; row >= 0; --row)
151  if (tone[row] != -1)
152  tones++;
153  fprintf (stderr, "(Image has %d graylevels.)\n", tones);
154 
155  /* Collapse array, taking out all zero values */
156  for (row = 0, itone = 0; row <= PGM_MAXMAXVAL; row++)
157  if (tone[row] != -1)
158  tone[itone++] = tone[row];
159  /* Now array contains only the gray levels present (in ascending order) */
160 
161  /* Allocate memory for gray-tone spatial dependence matrix */
162  P_matrix0 = matrix (0, tones, 0, tones);
163  P_matrix45 = matrix (0, tones, 0, tones);
164  P_matrix90 = matrix (0, tones, 0, tones);
165  P_matrix135 = matrix (0, tones, 0, tones);
166  for (row = 0; row < tones; ++row)
167  for (col = 0; col < tones; ++col)
168  {
169  P_matrix0[row][col] = P_matrix45[row][col] = 0;
170  P_matrix90[row][col] = P_matrix135[row][col] = 0;
171  }
172 
173  /* Find gray-tone spatial dependence matrix */
174  fprintf (stderr, "(Computing spatial dependence matrix...");
175  for (row = 0; row < rows; ++row)
176  for (col = 0; col < cols; ++col)
177  for (x = 0, angle = 0; angle <= 135; angle += 45)
178  {
179  while (tone[x] != grays[row][col])
180  x++;
181  if (angle == 0 && col + d < cols)
182  {
183  y = 0;
184  while (tone[y] != grays[row][col + d])
185  y++;
186  P_matrix0[x][y]++;
187  P_matrix0[y][x]++;
188  }
189  if (angle == 90 && row + d < rows)
190  {
191  y = 0;
192  while (tone[y] != grays[row + d][col])
193  y++;
194  P_matrix90[x][y]++;
195  P_matrix90[y][x]++;
196  }
197  if (angle == 45 && row + d < rows && col - d >= 0)
198  {
199  y = 0;
200  while (tone[y] != grays[row + d][col - d])
201  y++;
202  P_matrix45[x][y]++;
203  P_matrix45[y][x]++;
204  }
205  if (angle == 135 && row + d < rows && col + d < cols)
206  {
207  y = 0;
208  while (tone[y] != grays[row + d][col + d])
209  y++;
210  P_matrix135[x][y]++;
211  P_matrix135[y][x]++;
212  }
213  }
214  /* Gray-tone spatial dependence matrices are complete */
215 
216  /* Find normalizing constants */
217  R0 = 2 * rows * (cols - 1);
218  R45 = 2 * (rows - 1) * (cols - 1);
219  R90 = 2 * (rows - 1) * cols;
220 
221  /* Normalize gray-tone spatial dependence matrix */
222  for (itone = 0; itone < tones; ++itone)
223  for (jtone = 0; jtone < tones; ++jtone)
224  {
225  P_matrix0[itone][jtone] /= R0;
226  P_matrix45[itone][jtone] /= R45;
227  P_matrix90[itone][jtone] /= R90;
228  P_matrix135[itone][jtone] /= R45;
229  }
230 
231  fprintf (stderr, " done.)\n");
232  fprintf (stderr, "(Computing textural features");
233  fprintf (stdout, "\n");
234  DOT;
235  fprintf (stdout,
236  "%s 0 45 90 135 Avg\n",
237  BL);
238 
239  ASM[0] = f1_asm (P_matrix0, tones);
240  ASM[1] = f1_asm (P_matrix45, tones);
241  ASM[2] = f1_asm (P_matrix90, tones);
242  ASM[3] = f1_asm (P_matrix135, tones);
243  results (F1, ASM);
244 
245  contrast[0] = f2_contrast (P_matrix0, tones);
246  contrast[1] = f2_contrast (P_matrix45, tones);
247  contrast[2] = f2_contrast (P_matrix90, tones);
248  contrast[3] = f2_contrast (P_matrix135, tones);
249  results (F2, contrast);
250 
251 
252  corr[0] = f3_corr (P_matrix0, tones);
253  corr[1] = f3_corr (P_matrix45, tones);
254  corr[2] = f3_corr (P_matrix90, tones);
255  corr[3] = f3_corr (P_matrix135, tones);
256  results (F3, corr);
257 
258  var[0] = f4_var (P_matrix0, tones);
259  var[1] = f4_var (P_matrix45, tones);
260  var[2] = f4_var (P_matrix90, tones);
261  var[3] = f4_var (P_matrix135, tones);
262  results (F4, var);
263 
264 
265  idm[0] = f5_idm (P_matrix0, tones);
266  idm[1] = f5_idm (P_matrix45, tones);
267  idm[2] = f5_idm (P_matrix90, tones);
268  idm[3] = f5_idm (P_matrix135, tones);
269  results (F5, idm);
270 
271  savg[0] = f6_savg (P_matrix0, tones);
272  savg[1] = f6_savg (P_matrix45, tones);
273  savg[2] = f6_savg (P_matrix90, tones);
274  savg[3] = f6_savg (P_matrix135, tones);
275  results (F6, savg);
276 
277  sentropy[0] = f8_sentropy (P_matrix0, tones);
278  sentropy[1] = f8_sentropy (P_matrix45, tones);
279  sentropy[2] = f8_sentropy (P_matrix90, tones);
280  sentropy[3] = f8_sentropy (P_matrix135, tones);
281  svar[0] = f7_svar (P_matrix0, tones, sentropy[0]);
282  svar[1] = f7_svar (P_matrix45, tones, sentropy[1]);
283  svar[2] = f7_svar (P_matrix90, tones, sentropy[2]);
284  svar[3] = f7_svar (P_matrix135, tones, sentropy[3]);
285  results (F7, svar);
286  results (F8, sentropy);
287 
288  entropy[0] = f9_entropy (P_matrix0, tones);
289  entropy[1] = f9_entropy (P_matrix45, tones);
290  entropy[2] = f9_entropy (P_matrix90, tones);
291  entropy[3] = f9_entropy (P_matrix135, tones);
292  results (F9, entropy);
293 
294  dvar[0] = f10_dvar (P_matrix0, tones);
295  dvar[1] = f10_dvar (P_matrix45, tones);
296  dvar[2] = f10_dvar (P_matrix90, tones);
297  dvar[3] = f10_dvar (P_matrix135, tones);
298  results (F10, dvar);
299 
300  dentropy[0] = f11_dentropy (P_matrix0, tones);
301  dentropy[1] = f11_dentropy (P_matrix45, tones);
302  dentropy[2] = f11_dentropy (P_matrix90, tones);
303  dentropy[3] = f11_dentropy (P_matrix135, tones);
304  results (F11, dentropy);
305 
306  icorr[0] = f12_icorr (P_matrix0, tones);
307  icorr[1] = f12_icorr (P_matrix45, tones);
308  icorr[2] = f12_icorr (P_matrix90, tones);
309  icorr[3] = f12_icorr (P_matrix135, tones);
310  results (F12, icorr);
311 
312  icorr[0] = f13_icorr (P_matrix0, tones);
313  icorr[1] = f13_icorr (P_matrix45, tones);
314  icorr[2] = f13_icorr (P_matrix90, tones);
315  icorr[3] = f13_icorr (P_matrix135, tones);
316  results (F13, icorr);
317 
318  maxcorr[0] = f14_maxcorr (P_matrix0, tones);
319  maxcorr[1] = f14_maxcorr (P_matrix45, tones);
320  maxcorr[2] = f14_maxcorr (P_matrix90, tones);
321  maxcorr[3] = f14_maxcorr (P_matrix135, tones);
322  results (F14, maxcorr);
323 
324 
325  fprintf (stderr, " done.)\n");
326  exit (0);
327 }
328 
329 float f1_asm (P, Ng)
330  float **P;
331  int Ng;
332 
333 /* Angular Second Moment */
334 {
335  int i, j;
336  float sum = 0;
337 
338  for (i = 0; i < Ng; ++i)
339  for (j = 0; j < Ng; ++j)
340  sum += P[i][j] * P[i][j];
341 
342  return sum;
343 
344  /*
345  * The angular second-moment feature (ASM) f1 is a measure of homogeneity
346  * of the image. In a homogeneous image, there are very few dominant
347  * gray-tone transitions. Hence the P matrix for such an image will have
348  * fewer entries of large magnitude.
349  */
350 }
351 
352 
353 float f2_contrast (P, Ng)
354  float **P;
355  int Ng;
356 
357 /* Contrast */
358 {
359  int i, j, n;
360  float sum = 0, bigsum = 0;
361 
362  for (n = 0; n < Ng; ++n)
363  {
364  for (i = 0; i < Ng; ++i)
365  for (j = 0; j < Ng; ++j)
366  if ((i - j) == n || (j - i) == n)
367  sum += P[i][j];
368  bigsum += n * n * sum;
369 
370  sum = 0;
371  }
372  return bigsum;
373 
374  /*
375  * The contrast feature is a difference moment of the P matrix and is a
376  * measure of the contrast or the amount of local variations present in an
377  * image.
378  */
379 }
380 
381 float f3_corr (P, Ng)
382  float **P;
383  int Ng;
384 
385 /* Correlation */
386 {
387  int i, j;
388  float sum_sqrx = 0, sum_sqry = 0, tmp, *px;
389  float meanx =0 , meany = 0 , stddevx, stddevy;
390 
391  px = vector (0, Ng);
392  for (i = 0; i < Ng; ++i)
393  px[i] = 0;
394 
395  /*
396  * px[i] is the (i-1)th entry in the marginal probability matrix obtained
397  * by summing the rows of p[i][j]
398  */
399  for (i = 0; i < Ng; ++i)
400  for (j = 0; j < Ng; ++j)
401  px[i] += P[i][j];
402 
403 
404  /* Now calculate the means and standard deviations of px and py */
405  /*- fix supplied by J. Michael Christensen, 21 Jun 1991 */
406  /*- further modified by James Darrell McCauley, 16 Aug 1991
407  * after realizing that meanx=meany and stddevx=stddevy
408  */
409  for (i = 0; i < Ng; ++i)
410  {
411  meanx += px[i]*i;
412  sum_sqrx += px[i]*i*i;
413  }
414  meany = meanx;
415  sum_sqry = sum_sqrx;
416  stddevx = sqrt (sum_sqrx - (meanx * meanx));
417  stddevy = stddevx;
418 
419  /* Finally, the correlation ... */
420  for (tmp = 0, i = 0; i < Ng; ++i)
421  for (j = 0; j < Ng; ++j)
422  tmp += i*j*P[i][j];
423 
424  return (tmp - meanx * meany) / (stddevx * stddevy);
425  /*
426  * This correlation feature is a measure of gray-tone linear-dependencies
427  * in the image.
428  */
429 }
430 
431 
432 float f4_var (P, Ng)
433  float **P;
434  int Ng;
435 
436 /* Sum of Squares: Variance */
437 {
438  int i, j;
439  float mean = 0, var = 0;
440 
441  /*- Corrected by James Darrell McCauley, 16 Aug 1991
442  * calculates the mean intensity level instead of the mean of
443  * cooccurrence matrix elements
444  */
445  for (i = 0; i < Ng; ++i)
446  for (j = 0; j < Ng; ++j)
447  mean += i * P[i][j];
448 
449  for (i = 0; i < Ng; ++i)
450  for (j = 0; j < Ng; ++j)
451  var += (i + 1 - mean) * (i + 1 - mean) * P[i][j];
452 
453  return var;
454 }
455 
456 float f5_idm (P, Ng)
457  float **P;
458  int Ng;
459 
460 /* Inverse Difference Moment */
461 {
462  int i, j;
463  float idm = 0;
464 
465  for (i = 0; i < Ng; ++i)
466  for (j = 0; j < Ng; ++j)
467  idm += P[i][j] / (1 + (i - j) * (i - j));
468 
469  return idm;
470 }
471 
472 float Pxpy[2 * PGM_MAXMAXVAL];
473 
474 float f6_savg (P, Ng)
475  float **P;
476  int Ng;
477 
478 /* Sum Average */
479 {
480  int i, j;
481  extern float Pxpy[2 * PGM_MAXMAXVAL];
482  float savg = 0;
483 
484  for (i = 0; i <= 2 * Ng; ++i)
485  Pxpy[i] = 0;
486 
487  for (i = 0; i < Ng; ++i)
488  for (j = 0; j < Ng; ++j)
489  Pxpy[i + j + 2] += P[i][j];
490  for (i = 2; i <= 2 * Ng; ++i)
491  savg += i * Pxpy[i];
492 
493  return savg;
494 }
495 
496 
497 float f7_svar (float **P, int Ng, float S) {
498 /* Sum Variance */
499  int i, j;
500  extern float Pxpy[2 * PGM_MAXMAXVAL];
501  float var = 0;
502 
503  for (i = 0; i <= 2 * Ng; ++i)
504  Pxpy[i] = 0;
505 
506  for (i = 0; i < Ng; ++i)
507  for (j = 0; j < Ng; ++j)
508  Pxpy[i + j + 2] += P[i][j];
509 
510  for (i = 2; i <= 2 * Ng; ++i)
511  var += (i - S) * (i - S) * Pxpy[i];
512 
513  return var;
514 }
515 
516 float f8_sentropy (P, Ng)
517  float **P;
518  int Ng;
519 
520 /* Sum Entropy */
521 {
522  int i, j;
523  extern float Pxpy[2 * PGM_MAXMAXVAL];
524  float sentropy = 0;
525 
526  for (i = 0; i <= 2 * Ng; ++i)
527  Pxpy[i] = 0;
528 
529  for (i = 0; i < Ng; ++i)
530  for (j = 0; j < Ng; ++j)
531  Pxpy[i + j + 2] += P[i][j];
532 
533  for (i = 2; i <= 2 * Ng; ++i)
534  sentropy -= Pxpy[i] * log10 (Pxpy[i] + EPSILON);
535 
536  return sentropy;
537 }
538 
539 
540 float f9_entropy (P, Ng)
541  float **P;
542  int Ng;
543 
544 /* Entropy */
545 {
546  int i, j;
547  float entropy = 0;
548 
549  for (i = 0; i < Ng; ++i)
550  for (j = 0; j < Ng; ++j)
551  entropy += P[i][j] * log10 (P[i][j] + EPSILON);
552 
553  return -entropy;
554 }
555 
556 
557 float f10_dvar (P, Ng)
558  float **P;
559  int Ng;
560 
561 /* Difference Variance */
562 {
563  int i, j, tmp;
564  extern float Pxpy[2 * PGM_MAXMAXVAL];
565  float sum = 0, sum_sqr = 0, var = 0;
566 
567  for (i = 0; i <= 2 * Ng; ++i)
568  Pxpy[i] = 0;
569 
570  for (i = 0; i < Ng; ++i)
571  for (j = 0; j < Ng; ++j)
572  Pxpy[abs (i - j)] += P[i][j];
573 
574  /* Now calculate the variance of Pxpy (Px-y) */
575  for (i = 0; i < Ng; ++i)
576  {
577  sum += Pxpy[i];
578  sum_sqr += Pxpy[i] * Pxpy[i];
579  }
580  tmp = Ng * Ng;
581  var = ((tmp * sum_sqr) - (sum * sum)) / (tmp * tmp);
582 
583  return var;
584 }
585 
586 float f11_dentropy (P, Ng)
587  float **P;
588  int Ng;
589 
590 /* Difference Entropy */
591 {
592  int i, j;
593  extern float Pxpy[2 * PGM_MAXMAXVAL];
594  float sum = 0;
595 
596  for (i = 0; i <= 2 * Ng; ++i)
597  Pxpy[i] = 0;
598 
599  for (i = 0; i < Ng; ++i)
600  for (j = 0; j < Ng; ++j)
601  Pxpy[abs (i - j)] += P[i][j];
602 
603  for (i = 0; i < Ng; ++i)
604  sum += Pxpy[i] * log10 (Pxpy[i] + EPSILON);
605 
606  return -sum;
607 }
608 
609 float f12_icorr (P, Ng)
610  float **P;
611  int Ng;
612 
613 /* Information Measures of Correlation */
614 {
615  int i, j;
616  float *px, *py;
617  float hx = 0, hy = 0, hxy = 0, hxy1 = 0, hxy2 = 0;
618 
619  px = vector (0, Ng);
620  py = vector (0, Ng);
621 
622  /*
623  * px[i] is the (i-1)th entry in the marginal probability matrix obtained
624  * by summing the rows of p[i][j]
625  */
626  for (i = 0; i < Ng; ++i)
627  {
628  for (j = 0; j < Ng; ++j)
629  {
630  px[i] += P[i][j];
631  py[j] += P[i][j];
632  }
633  }
634 
635  for (i = 0; i < Ng; ++i)
636  for (j = 0; j < Ng; ++j)
637  {
638  hxy1 -= P[i][j] * log10 (px[i] * py[j] + EPSILON);
639  hxy2 -= px[i] * py[j] * log10 (px[i] * py[j] + EPSILON);
640  hxy -= P[i][j] * log10 (P[i][j] + EPSILON);
641  }
642 
643  /* Calculate entropies of px and py - is this right? */
644  for (i = 0; i < Ng; ++i)
645  {
646  hx -= px[i] * log10 (px[i] + EPSILON);
647  hy -= py[i] * log10 (py[i] + EPSILON);
648  }
649 /* fprintf(stderr,"hxy1=%f\thxy=%f\thx=%f\thy=%f\n",hxy1,hxy,hx,hy); */
650  return ((hxy - hxy1) / (hx > hy ? hx : hy));
651 }
652 
653 float f13_icorr (P, Ng)
654  float **P;
655  int Ng;
656 
657 /* Information Measures of Correlation */
658 {
659  int i, j;
660  float *px, *py;
661  float hx = 0, hy = 0, hxy = 0, hxy1 = 0, hxy2 = 0;
662 
663  px = vector (0, Ng);
664  py = vector (0, Ng);
665 
666  /*
667  * px[i] is the (i-1)th entry in the marginal probability matrix obtained
668  * by summing the rows of p[i][j]
669  */
670  for (i = 0; i < Ng; ++i)
671  {
672  for (j = 0; j < Ng; ++j)
673  {
674  px[i] += P[i][j];
675  py[j] += P[i][j];
676  }
677  }
678 
679  for (i = 0; i < Ng; ++i)
680  for (j = 0; j < Ng; ++j)
681  {
682  hxy1 -= P[i][j] * log10 (px[i] * py[j] + EPSILON);
683  hxy2 -= px[i] * py[j] * log10 (px[i] * py[j] + EPSILON);
684  hxy -= P[i][j] * log10 (P[i][j] + EPSILON);
685  }
686 
687  /* Calculate entropies of px and py */
688  for (i = 0; i < Ng; ++i)
689  {
690  hx -= px[i] * log10 (px[i] + EPSILON);
691  hy -= py[i] * log10 (py[i] + EPSILON);
692  }
693 /* fprintf(stderr,"hx=%f\thxy2=%f\n",hx,hxy2); */
694  return (sqrt (abs (1 - exp (-2.0 * (hxy2 - hxy)))));
695 }
696 
697 float f14_maxcorr (P, Ng)
698  float **P;
699  int Ng;
700 
701 /* Returns the Maximal Correlation Coefficient */
702 {
703  int i, j, k;
704  float *px, *py, **Q;
705  float *x, *iy, tmp;
706 
707  px = vector (0, Ng);
708  py = vector (0, Ng);
709  Q = matrix (1, Ng + 1, 1, Ng + 1);
710  x = vector (1, Ng);
711  iy = vector (1, Ng);
712 
713  /*
714  * px[i] is the (i-1)th entry in the marginal probability matrix obtained
715  * by summing the rows of p[i][j]
716  */
717  for (i = 0; i < Ng; ++i)
718  {
719  for (j = 0; j < Ng; ++j)
720  {
721  px[i] += P[i][j];
722  py[j] += P[i][j];
723  }
724  }
725 
726  /* Find the Q matrix */
727  for (i = 0; i < Ng; ++i)
728  {
729  for (j = 0; j < Ng; ++j)
730  {
731  Q[i + 1][j + 1] = 0;
732  for (k = 0; k < Ng; ++k)
733  Q[i + 1][j + 1] += P[i][k] * P[j][k] / px[i] / py[k];
734  }
735  }
736 
737  /* Balance the matrix */
738  mkbalanced (Q, Ng);
739  /* Reduction to Hessenberg Form */
740  reduction (Q, Ng);
741  /* Finding eigenvalue for nonsymetric matrix using QR algorithm */
742  hessenberg (Q, Ng, x, iy);
743  /* simplesrt(Ng,x); */
744  /* Returns the sqrt of the second largest eigenvalue of Q */
745  for (i = 2, tmp = x[1]; i <= Ng; ++i)
746  tmp = (tmp > x[i]) ? tmp : x[i];
747  return sqrt (x[Ng - 1]);
748 }
749 
750 float *vector (nl, nh)
751  int nl, nh;
752 {
753  float *v;
754 
755  v = (float *) malloc ((unsigned) (nh - nl + 1) * sizeof (float));
756  if (!v)
757  fprintf (stderr, "memory allocation failure"), exit (1);
758  return v - nl;
759 }
760 
761 
762 float **matrix (nrl, nrh, ncl, nch)
763  int nrl, nrh, ncl, nch;
764 
765 /* Allocates a float matrix with range [nrl..nrh][ncl..nch] */
766 {
767  int i;
768  float **m;
769 
770  /* allocate pointers to rows */
771  m = (float **) malloc ((unsigned) (nrh - nrl + 1) * sizeof (float *));
772  if (!m)
773  fprintf (stderr, "memory allocation failure"), exit (1);
774  m -= ncl;
775 
776  /* allocate rows and set pointers to them */
777  for (i = nrl; i <= nrh; i++)
778  {
779  m[i] = (float *) malloc ((unsigned) (nch - ncl + 1) * sizeof (float));
780  if (!m[i])
781  fprintf (stderr, "memory allocation failure"), exit (2);
782  m[i] -= ncl;
783  }
784  /* return pointer to array of pointers to rows */
785  return m;
786 }
787 
788 void results (c, a)
789  char *c;
790  float *a;
791 {
792  int i;
793 
794  DOT;
795  fprintf (stdout, "%s", c);
796  for (i = 0; i < 4; ++i)
797  fprintf (stdout, "% 1.3e ", a[i]);
798  fprintf (stdout, "% 1.3e\n", (a[0] + a[1] + a[2] + a[3]) / 4);
799 }
800 
801 void simplesrt (n, arr)
802  int n;
803  float arr[];
804 {
805  int i, j;
806  float a;
807 
808  for (j = 2; j <= n; j++)
809  {
810  a = arr[j];
811  i = j - 1;
812  while (i > 0 && arr[i] > a)
813  {
814  arr[i + 1] = arr[i];
815  i--;
816  }
817  arr[i + 1] = a;
818  }
819 }
820 
821 void mkbalanced (a, n)
822  float **a;
823  int n;
824 {
825  int last, j, i;
826  float s, r, g, f, c, sqrdx;
827 
828  sqrdx = RADIX * RADIX;
829  last = 0;
830  while (last == 0)
831  {
832  last = 1;
833  for (i = 1; i <= n; i++)
834  {
835  r = c = 0.0;
836  for (j = 1; j <= n; j++)
837  if (j != i)
838  {
839  c += fabs (a[j][i]);
840  r += fabs (a[i][j]);
841  }
842  if (c && r)
843  {
844  g = r / RADIX;
845  f = 1.0;
846  s = c + r;
847  while (c < g)
848  {
849  f *= RADIX;
850  c *= sqrdx;
851  }
852  g = r * RADIX;
853  while (c > g)
854  {
855  f /= RADIX;
856  c /= sqrdx;
857  }
858  if ((c + r) / f < 0.95 * s)
859  {
860  last = 0;
861  g = 1.0 / f;
862  for (j = 1; j <= n; j++)
863  a[i][j] *= g;
864  for (j = 1; j <= n; j++)
865  a[j][i] *= f;
866  }
867  }
868  }
869  }
870 }
871 
872 
873 void reduction (a, n)
874  float **a;
875  int n;
876 {
877  int m, j, i;
878  float y, x;
879 
880  for (m = 2; m < n; m++)
881  {
882  x = 0.0;
883  i = m;
884  for (j = m; j <= n; j++)
885  {
886  if (fabs (a[j][m - 1]) > fabs (x))
887  {
888  x = a[j][m - 1];
889  i = j;
890  }
891  }
892  if (i != m)
893  {
894  for (j = m - 1; j <= n; j++)
895  SWAP (a[i][j], a[m][j])
896  for (j = 1; j <= n; j++)
897  SWAP (a[j][i], a[j][m])
898  a[j][i] = a[j][i];
899  }
900  if (x)
901  {
902  for (i = m + 1; i <= n; i++)
903  {
904  if ((y = a[i][m - 1]))
905  {
906  y /= x;
907  a[i][m - 1] = y;
908  for (j = m; j <= n; j++)
909  a[i][j] -= y * a[m][j];
910  for (j = 1; j <= n; j++)
911  a[j][m] += y * a[j][i];
912  }
913  }
914  }
915  }
916 }
917 
918 void hessenberg (a, n, wr, wi)
919  float **a, wr[], wi[];
920  int n;
921 
922 {
923  int nn, m, l, k, j, its, i, mmin;
924  float z, y, x, w, v, u, t, s, r, q, p, anorm;
925 
926  anorm = fabs (a[1][1]);
927  for (i = 2; i <= n; i++)
928  for (j = (i - 1); j <= n; j++)
929  anorm += fabs (a[i][j]);
930  nn = n;
931  t = 0.0;
932  while (nn >= 1)
933  {
934  its = 0;
935  do
936  {
937  for (l = nn; l >= 2; l--)
938  {
939  s = fabs (a[l - 1][l - 1]) + fabs (a[l][l]);
940  if (s == 0.0)
941  s = anorm;
942  if ((float) (fabs (a[l][l - 1]) + s) == s)
943  break;
944  }
945  x = a[nn][nn];
946  if (l == nn)
947  {
948  wr[nn] = x + t;
949  wi[nn--] = 0.0;
950  }
951  else
952  {
953  y = a[nn - 1][nn - 1];
954  w = a[nn][nn - 1] * a[nn - 1][nn];
955  if (l == (nn - 1))
956  {
957  p = 0.5 * (y - x);
958  q = p * p + w;
959  z = sqrt (fabs (q));
960  x += t;
961  if (q >= 0.0)
962  {
963  z = p + SIGN (z, p);
964  wr[nn - 1] = wr[nn] = x + z;
965  if (z)
966  wr[nn] = x - w / z;
967  wi[nn - 1] = wi[nn] = 0.0;
968  }
969  else
970  {
971  wr[nn - 1] = wr[nn] = x + p;
972  wi[nn - 1] = -(wi[nn] = z);
973  }
974  nn -= 2;
975  }
976  else
977  {
978  if (its == 30)
979  fprintf (stderr,
980  "Too many iterations to required to find %s\ngiving up\n",
981  F14), exit (1);
982  if (its == 10 || its == 20)
983  {
984  t += x;
985  for (i = 1; i <= nn; i++)
986  a[i][i] -= x;
987  s = fabs (a[nn][nn - 1]) + fabs (a[nn - 1][nn - 2]);
988  y = x = 0.75 * s;
989  w = -0.4375 * s * s;
990  }
991  ++its;
992  for (m = (nn - 2); m >= l; m--)
993  {
994  z = a[m][m];
995  r = x - z;
996  s = y - z;
997  p = (r * s - w) / a[m + 1][m] + a[m][m + 1];
998  q = a[m + 1][m + 1] - z - r - s;
999  r = a[m + 2][m + 1];
1000  s = fabs (p) + fabs (q) + fabs (r);
1001  p /= s;
1002  q /= s;
1003  r /= s;
1004  if (m == l)
1005  break;
1006  u = fabs (a[m][m - 1]) * (fabs (q) + fabs (r));
1007  v = fabs (p) * (fabs (a[m - 1][m - 1]) + fabs (z) + fabs (a[m + 1][m + 1]));
1008  if ((float) (u + v) == v)
1009  break;
1010  }
1011  for (i = m + 2; i <= nn; i++)
1012  {
1013  a[i][i - 2] = 0.0;
1014  if (i != (m + 2))
1015  a[i][i - 3] = 0.0;
1016  }
1017  for (k = m; k <= nn - 1; k++)
1018  {
1019  if (k != m)
1020  {
1021  p = a[k][k - 1];
1022  q = a[k + 1][k - 1];
1023  r = 0.0;
1024  if (k != (nn - 1))
1025  r = a[k + 2][k - 1];
1026  if ((x = fabs (p) + fabs (q) + fabs (r)))
1027  {
1028  p /= x;
1029  q /= x;
1030  r /= x;
1031  }
1032  }
1033  if ((s = SIGN (sqrt (p * p + q * q + r * r), p)))
1034  {
1035  if (k == m)
1036  {
1037  if (l != m)
1038  a[k][k - 1] = -a[k][k - 1];
1039  }
1040  else
1041  a[k][k - 1] = -s * x;
1042  p += s;
1043  x = p / s;
1044  y = q / s;
1045  z = r / s;
1046  q /= p;
1047  r /= p;
1048  for (j = k; j <= nn; j++)
1049  {
1050  p = a[k][j] + q * a[k + 1][j];
1051  if (k != (nn - 1))
1052  {
1053  p += r * a[k + 2][j];
1054  a[k + 2][j] -= p * z;
1055  }
1056  a[k + 1][j] -= p * y;
1057  a[k][j] -= p * x;
1058  }
1059  mmin = nn < k + 3 ? nn : k + 3;
1060  for (i = l; i <= mmin; i++)
1061  {
1062  p = x * a[i][k] + y * a[i][k + 1];
1063  if (k != (nn - 1))
1064  {
1065  p += z * a[i][k + 2];
1066  a[i][k + 2] -= p * r;
1067  }
1068  a[i][k + 1] -= p * q;
1069  a[i][k] -= p;
1070  }
1071  }
1072  }
1073  }
1074  }
1075  } while (l < nn - 1);
1076  }
1077 }
double __cdecl exp(double _X)
double __cdecl log10(double _X)
q
Definition: afm2pl.c:2287
int nl
Definition: afm2tfm.c:885
int nh
Definition: afm2tfm.c:885
static gray maxval
Definition: asciitopgm.c:38
#define n
Definition: t4ht.c:1290
int z
Definition: dviconv.c:26
int w
Definition: dviconv.c:26
int v
Definition: dviconv.c:10
double sqrt()
int sscanf()
mpz_t * f
Definition: gen-fib.c:34
#define s
Definition: afcover.h:80
#define t
Definition: afcover.h:96
static char usage[]
Definition: giftopnm.c:59
paragraph P
#define c(n)
Definition: gpos-common.c:150
#define a(n)
Definition: gpos-common.c:148
#define d(n)
Definition: gpos-common.c:151
int col
Definition: gsftopk.c:443
#define mmin(a, b)
Definition: blocksort.c:594
small capitals from c petite p scientific f u
Definition: afcover.h:88
small capitals from c petite p
Definition: afcover.h:72
small capitals from c petite p scientific i
Definition: afcover.h:80
void exit()
kerning y
Definition: ttdriver.c:212
#define fprintf
Definition: mendex.h:64
#define fabs(x)
Definition: cpascal.h:211
#define malloc
Definition: alloca.c:91
angle
Definition: cordic.py:17
float x
Definition: cordic.py:15
int k
Definition: otp-parser.c:70
unsigned ncl
Definition: parse_ofm.c:64
void pm_usage(char *usage)
Definition: libpbm1.c:343
FILE * pm_openr(char *name)
Definition: libpbm1.c:600
static int rows
Definition: pbmclean.c:15
static int cols
Definition: pbmmask.c:21
#define abs(a)
Definition: pbmplus.h:225
#define ARGS(alist)
Definition: pbmplus.h:235
gray ** pgm_readpgm(FILE *file, int *colsP, int *rowsP, gray *maxvalP)
Definition: libpgm1.c:155
void pgm_init(int *argcP, argv)
Definition: libpgm1.c:19
#define PGM_MAXMAXVAL
Definition: pgm.h:35
unsigned int gray
Definition: pgm.h:10
#define EPSILON
Definition: pgmtexture.c:47
#define SIGN(x, y)
Definition: pgmtexture.c:64
void reduction()
float ** matrix()
float f1_asm()
float f13_icorr()
#define F5
Definition: pgmtexture.c:53
#define F9
Definition: pgmtexture.c:57
#define F14
Definition: pgmtexture.c:62
void simplesrt()
float f3_corr()
#define RADIX
Definition: pgmtexture.c:46
#define F2
Definition: pgmtexture.c:50
#define F1
Definition: pgmtexture.c:49
float f14_maxcorr()
float * vector()
#define F3
Definition: pgmtexture.c:51
float f4_var()
#define F7
Definition: pgmtexture.c:55
#define F4
Definition: pgmtexture.c:52
void results()
float Pxpy[2 *255]
Definition: pgmtexture.c:472
float f8_sentropy()
#define DOT
Definition: pgmtexture.c:65
float f10_dvar()
float f9_entropy()
#define F11
Definition: pgmtexture.c:59
#define F10
Definition: pgmtexture.c:58
void mkbalanced()
float f6_savg()
float f7_svar()
#define SWAP(a, b)
Definition: pgmtexture.c:66
float f11_dentropy()
float f5_idm()
#define F12
Definition: pgmtexture.c:60
#define BL
Definition: pgmtexture.c:48
#define F6
Definition: pgmtexture.c:54
float f2_contrast()
#define F8
Definition: pgmtexture.c:56
void hessenberg()
float f12_icorr()
#define F13
Definition: pgmtexture.c:61
int main(int argc, argv)
Definition: pgmtexture.c:91
integer nn[24]
Definition: pmxab.c:90
#define pm_error
Definition: png22pnm.c:118
#define pm_close(file)
Definition: png22pnm.c:120
static int32_t last
Definition: ppagelist.c:29
int g
Definition: ppmqvga.c:68
int r
Definition: ppmqvga.c:68
static int row
Definition: ps2pk.c:587
Definition: dvips.h:235
Definition: sh.h:1782
Definition: sed.h:50
#define FILE
Definition: t1stdio.h:34
int j
Definition: t4ht.c:1589
m
Definition: tex4ht.c:3990
FILE * ifp
Definition: t1asm.c:88
@ S
Definition: ubidiimp.h:53
#define argv
Definition: xmain.c:270
#define argc
Definition: xmain.c:269
static zic_t corr[50]
Definition: zic.c:413
#define argn