home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 4 / DATAFILE_PDCD4.iso / utilities / utilsm / netpbmsca / pgm / c / pgmtexture < prev    next >
Encoding:
Text File  |  1993-10-04  |  23.8 KB  |  1,078 lines
  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
  91. main (argc, argv)
  92.   int argc;
  93.   char *argv[];
  94. {
  95.   FILE *ifp;
  96.   register gray **grays, *gP;
  97.   int tone[PGM_MAXMAXVAL], R0, R45, R90, angle, d = 1, x, y;
  98.   int argn, rows, cols, bps, padright, 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, nmaxval;
  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.   /* Determine the number of different gray scales (not maxval) */
  141.   for (row = PGM_MAXMAXVAL; row >= 0; --row)
  142.     tone[row] = -1;
  143.   for (row = rows - 1; row >= 0; --row)
  144.     for (col = 0; col < cols; ++col)
  145.       tone[grays[row][col]] = grays[row][col];
  146.   for (row = PGM_MAXMAXVAL, tones = 0; row >= 0; --row)
  147.     if (tone[row] != -1)
  148.       tones++;
  149.   fprintf (stderr, "(Image has %d graylevels.)\n", tones);
  150.  
  151.   /* Collapse array, taking out all zero values */
  152.   for (row = 0, itone = 0; row <= PGM_MAXMAXVAL; row++)
  153.     if (tone[row] != -1)
  154.       tone[itone++] = tone[row];
  155.   /* Now array contains only the gray levels present (in ascending order) */
  156.  
  157.   /* Allocate memory for gray-tone spatial dependence matrix */
  158.   P_matrix0 = matrix (0, tones, 0, tones);
  159.   P_matrix45 = matrix (0, tones, 0, tones);
  160.   P_matrix90 = matrix (0, tones, 0, tones);
  161.   P_matrix135 = matrix (0, tones, 0, tones);
  162.   for (row = 0; row < tones; ++row)
  163.     for (col = 0; col < tones; ++col)
  164.     {
  165.       P_matrix0[row][col] = P_matrix45[row][col] = 0;
  166.       P_matrix90[row][col] = P_matrix135[row][col] = 0;
  167.     }
  168.  
  169.   /* Find gray-tone spatial dependence matrix */
  170.   fprintf (stderr, "(Computing spatial dependence matrix...");
  171.   for (row = 0; row < rows; ++row)
  172.     for (col = 0; col < cols; ++col)
  173.       for (x = 0, angle = 0; angle <= 135; angle += 45)
  174.       {
  175.     while (tone[x] != grays[row][col])
  176.       x++;
  177.     if (angle == 0 && col + d < cols)
  178.     {
  179.       y = 0;
  180.       while (tone[y] != grays[row][col + d])
  181.         y++;
  182.       P_matrix0[x][y]++;
  183.       P_matrix0[y][x]++;
  184.     }
  185.     if (angle == 90 && row + d < rows)
  186.     {
  187.       y = 0;
  188.       while (tone[y] != grays[row + d][col])
  189.         y++;
  190.       P_matrix90[x][y]++;
  191.       P_matrix90[y][x]++;
  192.     }
  193.     if (angle == 45 && row + d < rows && col - d >= 0)
  194.     {
  195.       y = 0;
  196.       while (tone[y] != grays[row + d][col - d])
  197.         y++;
  198.       P_matrix45[x][y]++;
  199.       P_matrix45[y][x]++;
  200.     }
  201.     if (angle == 135 && row + d < rows && col + d < cols)
  202.     {
  203.       y = 0;
  204.       while (tone[y] != grays[row + d][col + d])
  205.         y++;
  206.       P_matrix135[x][y]++;
  207.       P_matrix135[y][x]++;
  208.     }
  209.       }
  210.   /* Gray-tone spatial dependence matrices are complete */
  211.  
  212.   /* Find normalizing constants */
  213.   R0 = 2 * rows * (cols - 1);
  214.   R45 = 2 * (rows - 1) * (cols - 1);
  215.   R90 = 2 * (rows - 1) * cols;
  216.  
  217.   /* Normalize gray-tone spatial dependence matrix */
  218.   for (itone = 0; itone < tones; ++itone)
  219.     for (jtone = 0; jtone < tones; ++jtone)
  220.     {
  221.       P_matrix0[itone][jtone] /= R0;
  222.       P_matrix45[itone][jtone] /= R45;
  223.       P_matrix90[itone][jtone] /= R90;
  224.       P_matrix135[itone][jtone] /= R45;
  225.     }
  226.  
  227.   fprintf (stderr, " done.)\n");
  228.   fprintf (stderr, "(Computing textural features");
  229.   fprintf (stdout, "\n");
  230.   DOT;
  231.   fprintf (stdout,
  232.        "%s         0         45         90        135        Avg\n",
  233.        BL);
  234.  
  235.   ASM[0] = f1_asm (P_matrix0, tones);
  236.   ASM[1] = f1_asm (P_matrix45, tones);
  237.   ASM[2] = f1_asm (P_matrix90, tones);
  238.   ASM[3] = f1_asm (P_matrix135, tones);
  239.   results (F1, ASM);
  240.  
  241.   contrast[0] = f2_contrast (P_matrix0, tones);
  242.   contrast[1] = f2_contrast (P_matrix45, tones);
  243.   contrast[2] = f2_contrast (P_matrix90, tones);
  244.   contrast[3] = f2_contrast (P_matrix135, tones);
  245.   results (F2, contrast);
  246.  
  247.  
  248.   corr[0] = f3_corr (P_matrix0, tones);
  249.   corr[1] = f3_corr (P_matrix45, tones);
  250.   corr[2] = f3_corr (P_matrix90, tones);
  251.   corr[3] = f3_corr (P_matrix135, tones);
  252.   results (F3, corr);
  253.  
  254.   var[0] = f4_var (P_matrix0, tones);
  255.   var[1] = f4_var (P_matrix45, tones);
  256.   var[2] = f4_var (P_matrix90, tones);
  257.   var[3] = f4_var (P_matrix135, tones);
  258.   results (F4, var);
  259.  
  260.  
  261.   idm[0] = f5_idm (P_matrix0, tones);
  262.   idm[1] = f5_idm (P_matrix45, tones);
  263.   idm[2] = f5_idm (P_matrix90, tones);
  264.   idm[3] = f5_idm (P_matrix135, tones);
  265.   results (F5, idm);
  266.  
  267.   savg[0] = f6_savg (P_matrix0, tones);
  268.   savg[1] = f6_savg (P_matrix45, tones);
  269.   savg[2] = f6_savg (P_matrix90, tones);
  270.   savg[3] = f6_savg (P_matrix135, tones);
  271.   results (F6, savg);
  272.  
  273.   sentropy[0] = f8_sentropy (P_matrix0, tones);
  274.   sentropy[1] = f8_sentropy (P_matrix45, tones);
  275.   sentropy[2] = f8_sentropy (P_matrix90, tones);
  276.   sentropy[3] = f8_sentropy (P_matrix135, tones);
  277.   svar[0] = f7_svar (P_matrix0, tones, sentropy[0]);
  278.   svar[1] = f7_svar (P_matrix45, tones, sentropy[1]);
  279.   svar[2] = f7_svar (P_matrix90, tones, sentropy[2]);
  280.   svar[3] = f7_svar (P_matrix135, tones, sentropy[3]);
  281.   results (F7, svar);
  282.   results (F8, sentropy);
  283.  
  284.   entropy[0] = f9_entropy (P_matrix0, tones);
  285.   entropy[1] = f9_entropy (P_matrix45, tones);
  286.   entropy[2] = f9_entropy (P_matrix90, tones);
  287.   entropy[3] = f9_entropy (P_matrix135, tones);
  288.   results (F9, entropy);
  289.  
  290.   dvar[0] = f10_dvar (P_matrix0, tones);
  291.   dvar[1] = f10_dvar (P_matrix45, tones);
  292.   dvar[2] = f10_dvar (P_matrix90, tones);
  293.   dvar[3] = f10_dvar (P_matrix135, tones);
  294.   results (F10, dvar);
  295.  
  296.   dentropy[0] = f11_dentropy (P_matrix0, tones);
  297.   dentropy[1] = f11_dentropy (P_matrix45, tones);
  298.   dentropy[2] = f11_dentropy (P_matrix90, tones);
  299.   dentropy[3] = f11_dentropy (P_matrix135, tones);
  300.   results (F11, dentropy);
  301.  
  302.   icorr[0] = f12_icorr (P_matrix0, tones);
  303.   icorr[1] = f12_icorr (P_matrix45, tones);
  304.   icorr[2] = f12_icorr (P_matrix90, tones);
  305.   icorr[3] = f12_icorr (P_matrix135, tones);
  306.   results (F12, icorr);
  307.  
  308.   icorr[0] = f13_icorr (P_matrix0, tones);
  309.   icorr[1] = f13_icorr (P_matrix45, tones);
  310.   icorr[2] = f13_icorr (P_matrix90, tones);
  311.   icorr[3] = f13_icorr (P_matrix135, tones);
  312.   results (F13, icorr);
  313.  
  314.   maxcorr[0] = f14_maxcorr (P_matrix0, tones);
  315.   maxcorr[1] = f14_maxcorr (P_matrix45, tones);
  316.   maxcorr[2] = f14_maxcorr (P_matrix90, tones);
  317.   maxcorr[3] = f14_maxcorr (P_matrix135, tones);
  318.   results (F14, maxcorr);
  319.  
  320.  
  321.   fprintf (stderr, " done.)\n");
  322.   exit (0);
  323. }
  324.  
  325. float f1_asm (P, Ng)
  326.   float **P;
  327.   int Ng;
  328.  
  329. /* Angular Second Moment */
  330. {
  331.   int i, j;
  332.   float sum = 0;
  333.  
  334.   for (i = 0; i < Ng; ++i)
  335.     for (j = 0; j < Ng; ++j)
  336.       sum += P[i][j] * P[i][j];
  337.  
  338.   return sum;
  339.  
  340.   /*
  341.    * The angular second-moment feature (ASM) f1 is a measure of homogeneity
  342.    * of the image. In a homogeneous image, there are very few dominant
  343.    * gray-tone transitions. Hence the P matrix for such an image will have
  344.    * fewer entries of large magnitude.
  345.    */
  346. }
  347.  
  348.  
  349. float f2_contrast (P, Ng)
  350.   float **P;
  351.   int Ng;
  352.  
  353. /* Contrast */
  354. {
  355.   int i, j, n;
  356.   float sum = 0, bigsum = 0;
  357.  
  358.   for (n = 0; n < Ng; ++n)
  359.   {
  360.     for (i = 0; i < Ng; ++i)
  361.       for (j = 0; j < Ng; ++j)
  362.     if ((i - j) == n || (j - i) == n)
  363.       sum += P[i][j];
  364.     bigsum += n * n * sum;
  365.  
  366.     sum = 0;
  367.   }
  368.   return bigsum;
  369.  
  370.   /*
  371.    * The contrast feature is a difference moment of the P matrix and is a
  372.    * measure of the contrast or the amount of local variations present in an
  373.    * image.
  374.    */
  375. }
  376.  
  377. float f3_corr (P, Ng)
  378.   float **P;
  379.   int Ng;
  380.  
  381. /* Correlation */
  382. {
  383.   int i, j;
  384.   float sum_sqrx = 0, sum_sqry = 0, tmp, *px;
  385.   float meanx =0 , meany = 0 , stddevx, stddevy;
  386.  
  387.   px = vector (0, Ng);
  388.   for (i = 0; i < Ng; ++i)
  389.     px[i] = 0;
  390.  
  391.   /*
  392.    * px[i] is the (i-1)th entry in the marginal probability matrix obtained
  393.    * by summing the rows of p[i][j]
  394.    */
  395.   for (i = 0; i < Ng; ++i)
  396.     for (j = 0; j < Ng; ++j)
  397.       px[i] += P[i][j];
  398.  
  399.  
  400.   /* Now calculate the means and standard deviations of px and py */
  401.   /*- fix supplied by J. Michael Christensen, 21 Jun 1991 */
  402.   /*- further modified by James Darrell McCauley, 16 Aug 1991 
  403.    *     after realizing that meanx=meany and stddevx=stddevy
  404.    */
  405.   for (i = 0; i < Ng; ++i)
  406.   {
  407.     meanx += px[i]*i;
  408.     sum_sqrx += px[i]*i*i;
  409.   }
  410.   meany = meanx;
  411.   sum_sqry = sum_sqrx;
  412.   stddevx = sqrt (sum_sqrx - (meanx * meanx));
  413.   stddevy = stddevx;
  414.  
  415.   /* Finally, the correlation ... */
  416.   for (tmp = 0, i = 0; i < Ng; ++i)
  417.     for (j = 0; j < Ng; ++j)
  418.       tmp += i*j*P[i][j];
  419.  
  420.   return (tmp - meanx * meany) / (stddevx * stddevy);
  421.   /*
  422.    * This correlation feature is a measure of gray-tone linear-dependencies
  423.    * in the image.
  424.    */
  425. }
  426.  
  427.  
  428. float f4_var (P, Ng)
  429.   float **P;
  430.   int Ng;
  431.  
  432. /* Sum of Squares: Variance */
  433. {
  434.   int i, j;
  435.   float mean = 0, var = 0;
  436.  
  437.   /*- Corrected by James Darrell McCauley, 16 Aug 1991
  438.    *  calculates the mean intensity level instead of the mean of
  439.    *  cooccurrence matrix elements 
  440.    */
  441.   for (i = 0; i < Ng; ++i)
  442.     for (j = 0; j < Ng; ++j)
  443.       mean += i * P[i][j];
  444.  
  445.   for (i = 0; i < Ng; ++i)
  446.     for (j = 0; j < Ng; ++j)
  447.       var += (i + 1 - mean) * (i + 1 - mean) * P[i][j];
  448.  
  449.   return var;
  450. }
  451.  
  452. float f5_idm (P, Ng)
  453.   float **P;
  454.   int Ng;
  455.  
  456. /* Inverse Difference Moment */
  457. {
  458.   int i, j;
  459.   float idm = 0;
  460.  
  461.   for (i = 0; i < Ng; ++i)
  462.     for (j = 0; j < Ng; ++j)
  463.       idm += P[i][j] / (1 + (i - j) * (i - j));
  464.  
  465.   return idm;
  466. }
  467.  
  468. float Pxpy[2 * PGM_MAXMAXVAL];
  469.  
  470. float f6_savg (P, Ng)
  471.   float **P;
  472.   int Ng;
  473.  
  474. /* Sum Average */
  475. {
  476.   int i, j;
  477.   extern float Pxpy[2 * PGM_MAXMAXVAL];
  478.   float savg = 0;
  479.  
  480.   for (i = 0; i <= 2 * Ng; ++i)
  481.     Pxpy[i] = 0;
  482.  
  483.   for (i = 0; i < Ng; ++i)
  484.     for (j = 0; j < Ng; ++j)
  485.       Pxpy[i + j + 2] += P[i][j];
  486.   for (i = 2; i <= 2 * Ng; ++i)
  487.     savg += i * Pxpy[i];
  488.  
  489.   return savg;
  490. }
  491.  
  492.  
  493. float f7_svar (P, Ng, S)
  494.   float **P, S;
  495.   int Ng;
  496.  
  497. /* Sum Variance */
  498. {
  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, tmp;
  593.   extern float Pxpy[2 * PGM_MAXMAXVAL];
  594.   float sum = 0, sum_sqr = 0, var = 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, tmp;
  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. }
  1078.