home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / CH_5.5 / FITLINE / eigenline.c next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  2.1 KB  |  77 lines
  1. /* EIGENLINE:   function performs eigenvector fit to data points and
  2.  *            yields slope and intercept of line
  3.  *                      line = eigenline (pt, n, &m, &b)
  4.  *              If line=0, then fit has not been performed due to 1-pixel line.
  5.  */
  6.  
  7. #include <stdio.h>
  8. #include <stdlib.h>
  9. #include <math.h>
  10. #include <images.h>             /* for sun images */
  11.  
  12. #define INFINITY 1000000.0      /* infinity slope value */
  13.  
  14. short
  15. eigenline (pt, n, m, b)
  16.      struct point *pt;          /* input set of coordinate points */
  17.      long n;                    /* number of input points */
  18.      double *m, *b;             /* slope and intercept of line */
  19.  
  20. {
  21.   struct dpoint *data;          /* floating point data vector */
  22.   struct dpoint ptAvg;          /* average of data points */
  23.   double a11, a12, a21, a22;    /* matrix elements */
  24.   double temp;                  /* temporary variable */
  25.   double lambda1;               /* minimum eigenvalue */
  26.   long i;
  27.  
  28. /* check 1-point "line" */
  29.   if (n == 0)
  30.     return (0);
  31.  
  32. /* allocate floating point vector for data */
  33.   if ((data = (struct dpoint *) calloc (n, sizeof (struct dpoint))) == NULL) {
  34.     printf ("EIGENLINE: not enough memory -- sorry");
  35.     return (-1);
  36.   }
  37.   for (i = 0; i < n; i++) {
  38.     data[i].x = (double) pt[i].x;
  39.     data[i].y = (double) pt[i].y;
  40.   }
  41.  
  42. /* calculate average and translate each point by the average to around (0,0) */
  43.   ptAvg.x = ptAvg.y = 0.0;
  44.   for (i = 0; i < n; i++) {
  45.     ptAvg.x += data[i].x;
  46.     ptAvg.y += data[i].y;
  47.   }
  48.   ptAvg.x /= (double) n;
  49.   ptAvg.y /= (double) n;
  50.   for (i = 0; i < n; i++) {
  51.     data[i].x -= ptAvg.x;
  52.     data[i].y -= ptAvg.y;
  53.   }
  54.  
  55. /* determine matrix */
  56.   a11 = a12 = a21 = a22 = 0;
  57.   for (i = 0; i < n; i++) {
  58.     a11 += data[i].x * data[i].x;
  59.     a12 = a21 += data[i].x * data[i].y;
  60.     a22 += data[i].y * data[i].y;
  61.   }
  62.  
  63. /* determine minimum eigenvalue */
  64.   temp = (a11 + a22);
  65.   lambda1 = (temp - sqrt (temp * temp - 4.0 * a11 * a22 + 4.0 * a12 * a12))
  66.     / 2.0;
  67.  
  68. /* calculate slope and intercept */
  69.   if (a11 == lambda1)
  70.     *m = INFINITY;
  71.   else
  72.     *m = -a12 / (lambda1 - a11);
  73.   *b = ptAvg.y - *m * ptAvg.x;
  74.  
  75.   return (1);
  76. }
  77.