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