home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / CH_5.4 / FITCRIT / thetafeats.c < prev    next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  57.9 KB  |  1,548 lines
  1. /* THETAFEAT:   function determines features and feature descriptions
  2.  *            from theta plot
  3.  *                      usage: nThetaFeats = thetafeat (data, edge, theta,
  4.  *                                              nTheta, wwPar)
  5.  *           -  Note that the <zCross> array contains a flag which
  6.  *              is UP, DOWN, ZERO for the *first sample* as read 
  7.  *              sequentially in the respective regions.
  8.  *           -  Note that if there is a cycle, the <data> array is
  9.  *              augmented by the wrap-around in function CYCLEAUGMENT.
  10.  *              However this array is NOT sent back to the calling program.
  11.  */
  12.  
  13. #include <stdio.h>
  14. #include <stdlib.h>
  15. #include <math.h>
  16. #include "pcc2.h"               /* wing-walk include file */
  17.  
  18. #define ROOT2 1.41421356237     /* square root of 2 */
  19. #if defined(WIN32)
  20. #define PI 3.1415927
  21. #endif
  22.  
  23. #define UP 3                    /* codes for zero ranges, above, and below */
  24. #define ZERO 2
  25. #define DOWN 1
  26. #define NOCROSSING 0
  27.  
  28. extern long maxCorner;          /* maximum arc length of corner curvature */
  29.     /* needed for zzWidth calc in FEATDETERM */
  30.  
  31. extern int linkcorner (int, struct point);
  32. extern int linkcurve (struct dpoint, struct point, struct point,
  33.                       struct dpoint, double, short);
  34.  
  35. long thetazz (struct theta *, long, short *, struct wwPar, long *, long *);
  36. long boundzz (struct theta *, long, short **, long *, struct point **,
  37.               struct edge *, long *, long *);
  38. long featdeterm (short *, long, long, struct theta *, long,
  39.                  struct point *, struct edge, struct wwPar,
  40.                  long *, long *, short);
  41. long featcorner (struct point *, long, long, long, long,
  42.                  struct edge, struct point *);
  43. long featcurve (struct point *, long, long, long, long, struct edge,
  44.                 double, double, struct dpoint *, struct point *,
  45.                 struct point *, struct dpoint *, double *, short *);
  46. long eventranspts (struct dpoint *, struct dpoint *,
  47.                    struct point, struct point, struct dpoint, double *);
  48. short cycleaugment (long, short **, long *, struct point **,
  49.                     struct edge *, long, long);
  50. long bounddeterm2 (short *, struct theta *, long, long,
  51.                    long, struct point *, struct edge, long *, long *);
  52. long bounddeterm (short *, struct theta *, long, long, long,
  53.                   struct point *, long *, long *);
  54. long featbound (struct point, struct point, struct point,
  55.                 double *, struct dpoint *, struct dpoint *, short *);
  56. long curvecorner (struct point, struct point, struct point,
  57.                   double, struct dpoint, struct dpoint *, short *);
  58. long feat2bound (struct point, struct point, double, struct dpoint,
  59.                  double *, struct dpoint *, short *);
  60. double radsearch (double, double, double);
  61. short curvedirn (struct point, struct point, struct dpoint, struct dpoint);
  62. long linesintrsct (struct point, struct dpoint, struct point,
  63.                    struct point, struct dpoint *);
  64.  
  65. long
  66. thetafeat (data, edge, theta, nTheta, wwPar, curveFlag)
  67.      struct point *data;        /* sequence of coordinates of lines */
  68.      struct edge edge;          /* lower and upper <data> indices */
  69.      struct theta *theta;       /* theta plot */
  70.      long nTheta;               /* no. values in theta plot */
  71.      struct wwPar wwPar;        /* wing-walk parameters */
  72.      short curveFlag;           /* =0 for curve fit by 2 lines; =1 for curve */
  73. {
  74.   short *zCross;                /* array of flags at locations of zero 
  75.                                  * * crossings on theta plot */
  76.   long zStart,                  /* starting zero crossing of peak */
  77.     zEnd,                       /* ending zero crossing of peak */
  78.     nThetaFeats,                /* no. of features found on theta plot */
  79.     i;
  80.   long nZCross,                 /* no. of values in zero-crossing plot */
  81.     zLow,                       /* line indices marking extent of feature to    */
  82.     zHigh;                      /*  halfway (+1) within line below feat 
  83.                                  * * and halfway (-1) within line after feat;
  84.                                  * * unless at either end where it goes to end - 1 */
  85.   long zInitial,                /* initial/final valid <zCross> indices */
  86.     zFinal;                     /*   inside of border values on theta plot */
  87.  
  88. /* find significant peak widths on theta plot */
  89.   nZCross = nTheta;
  90.   if ((zCross = (short *) calloc (nZCross, sizeof (*zCross))) == NULL) {
  91.     printf ("CALLOC: not enough memory -- sorry");
  92.     return (-1);
  93.   }
  94.  
  95.   thetazz (theta, nZCross, zCross, wwPar, &zInitial, &zFinal);
  96.  
  97.   if ((boundzz (theta, wwPar.l, &zCross, &nZCross,
  98.                 &data, &edge, &zInitial, &zFinal))
  99.       == 0)
  100.     return (0);
  101.  
  102. /* perform feature determination from peaks on theta plot */
  103.   for (zStart = 1; (zStart < nZCross)
  104.        && (zCross[zStart] == ZERO || zCross[zStart] == NOCROSSING);
  105.        zStart++);
  106.   for (i = zStart + 1, nThetaFeats = 0; i < nZCross; i++) {
  107.     if (zCross[i] != NOCROSSING) {
  108.       zEnd = i - 1;
  109.       if (zStart == zInitial || i == zFinal)
  110.         bounddeterm2 (zCross, theta, zStart, zEnd, nZCross, data,
  111.                       edge, &zLow, &zHigh);
  112.       else
  113.         featdeterm (zCross, zStart, zEnd, theta, nZCross, data,
  114.                     edge, wwPar, &zLow, &zHigh, curveFlag);
  115.       nThetaFeats++;
  116.       /* <= 1 added (formerly == 0) 6 Oct 88 */
  117.       zLow = (zLow <= 1) ? 1 : zLow + (zStart - zLow) / 2;
  118.       zHigh = (zHigh == nZCross - 1)
  119.         ? zHigh - 1 : zEnd + (zHigh - zEnd) / 2;
  120.       for (zStart = i; (zStart < nZCross)
  121.            && (zCross[zStart] == ZERO || zCross[zStart] == NOCROSSING);
  122.            zStart++);
  123.       i = zStart + 1;
  124.     }
  125.   }
  126.  
  127.   free (zCross);
  128.   return (nThetaFeats);
  129. }
  130.  
  131.  
  132. /* THETAZZ:     function finds significantly wide peak widths on theta plot
  133.  *                 usage: thetazz (theta, nTheta, zCross, wwPar,
  134.  *                                              &zInitial, &zFinal)
  135.  *              Note that although the peaks are represented by the
  136.  *              locations of the first and last non-zero region theta values
  137.  *              (iLow and i), the peak width is measured between zero
  138.  *              crossings, thus <iLow - 1> and <i + 1> below.
  139.  */
  140.  
  141. long
  142. thetazz (theta, nTheta, zCross, wwPar, zInitial, zFinal)
  143.      struct theta *theta;       /* theta plot */
  144.      long nTheta;               /* no. values in theta plot */
  145.      short *zCross;             /* array of zero crossings of theta plot */
  146.      struct wwPar wwPar;        /* wing-walk parameters */
  147.      long *zInitial,            /* initial/final valid <zCross> indices */
  148.       *zFinal;                  /*   inside of border values on theta plot */
  149. {
  150.   double zeroRange,             /* threshold for zero (strt line) region */
  151.     minZZ,                      /* minimum zero crossing length for peak */
  152.     lengthZZ,                   /* zero-crossing length */
  153.     atan2 ();
  154.   short range, rangeB4;         /* zero, up, or down regions */
  155.   int zStart,                   /* first zero crossing on theta plot */
  156.     iLow,                       /* lower index of zero crossing pair */
  157.     i;
  158.   struct narrowest {            /* narrowest zero crossing pair */
  159.     double length;              /* peak arc length */
  160.     long iLow, iHigh;           /* indices of start/end of peak on plot */
  161.   } narrowest;
  162.  
  163.   zeroRange = atan2 (2.0, wwPar.l - wwPar.m - 2.0 * ROOT2 + 2.0);
  164.  
  165. /* find first and last valid theta points, and set <zCross> to NOCROSSING
  166.  * in these border regions */
  167.   for (*zInitial = 0; *zInitial < nTheta; (*zInitial)++) {
  168.     if (theta[*zInitial].theta != BORDERTHETA)
  169.       break;
  170.     zCross[*zInitial] = NOCROSSING;
  171.   }
  172.   for (*zFinal = nTheta - 1; *zFinal >= 0; --(*zFinal)) {
  173.     if (theta[*zFinal].theta != BORDERTHETA)
  174.       break;
  175.     zCross[*zFinal] = NOCROSSING;
  176.   }
  177.  
  178. /* find zero crossings on theta plot */
  179.   if (theta[*zInitial].theta > zeroRange)
  180.     rangeB4 = UP;
  181.   else if (theta[*zInitial].theta < -zeroRange)
  182.     rangeB4 = DOWN;
  183.   else
  184.     rangeB4 = ZERO;
  185.   zCross[*zInitial] = rangeB4;
  186.   for (i = *zInitial + 1; i <= *zFinal; i++) {
  187.     if (theta[i].theta > zeroRange)
  188.       range = UP;
  189.     else if (theta[i].theta < -zeroRange)
  190.       range = DOWN;
  191.     else
  192.       range = ZERO;
  193.     zCross[i] = (range != rangeB4) ? range : NOCROSSING;
  194.     rangeB4 = range;
  195.   }
  196.  
  197. /* if any zero crossings are skipped (UP to DOWN, or vice versa without
  198.  * a ZERO in between) then put ZERO in between, or if the crossings are
  199.  * 1 or 2 samples apart then set lower crossing to NOCROSSING */
  200.   zStart = *zInitial;
  201.   for (i = *zInitial + 1; i <= *zFinal; i++) {
  202.     if (zCross[i] != NOCROSSING) {
  203.       if ((zCross[zStart] == UP && zCross[i] == DOWN)
  204.           || (zCross[zStart] == DOWN && zCross[i] == UP)) {
  205.         zCross[i - 1] = (zCross[i - 1] == NOCROSSING)
  206.           ? ZERO : NOCROSSING;
  207.       }
  208.       zStart = i;
  209.     }
  210.   }
  211.  
  212. /* eliminate zero crossings for peak widths which are too narrow */
  213.   minZZ = (double) (wwPar.l - 1) + EPSILON1;
  214.   while (1) {
  215.     narrowest.length = theta[nTheta - 1].length;
  216.     rangeB4 = zCross[*zInitial];
  217.     for (zStart = *zInitial + 1;
  218.          (zCross[zStart] == NOCROSSING) && zStart < *zFinal; zStart++);
  219.     if (zStart == *zFinal)
  220.       return (0);
  221.     for (iLow = zStart, i = zStart + 1; i <= *zFinal; i++) {
  222.       if (zCross[i] != NOCROSSING) {  /* see note above */
  223.         lengthZZ = theta[i + 1].length - theta[iLow - 1].length;
  224.         if (lengthZZ < narrowest.length && rangeB4 == zCross[i]) {
  225.           narrowest.length = theta[i + 1].length
  226.             - theta[iLow - 1].length;
  227.           narrowest.iLow = iLow;
  228.           narrowest.iHigh = i - 1;
  229.         }
  230.         rangeB4 = zCross[iLow];
  231.         iLow = i;
  232.       }
  233.     }
  234.     if (narrowest.length < minZZ)
  235.       zCross[narrowest.iLow] = zCross[narrowest.iHigh + 1]
  236.         = NOCROSSING;
  237.     else
  238.       break;
  239.   }
  240.  
  241.   return (0);
  242. }
  243.  
  244.  
  245. /* BOUNDZZ:     function adjusts zero crossings at bounds of line
  246.  *            or at seam of cycle
  247.  *                  usage: featFlag = boundzz (theta, wwParL
  248.  *                                          &zCross, &nZCross, &data, &edge,
  249.  *                                                       &zInitial, &zFinal)
  250.  *              Function returns 0 if cycle has no zero-crossings; 1 
  251.  *              otherwise.
  252.  */
  253.  
  254. long
  255. boundzz (theta, wwParL, zCross, nZCross, data, edge,
  256.          zInitial, zFinal)
  257.      struct theta *theta;       /* theta plot */
  258.      long wwParL;               /* ww length parameter */
  259.      short **zCross;            /* array of flags at zero-crossing locns */
  260.      long *nZCross;             /* no. values in zero-crossing plot */
  261.      struct point **data;       /* sequence of coordinates of lines */
  262.      struct edge *edge;         /* lower and upper <data> indices */
  263.      long *zInitial, *zFinal;   /* initial and final valid zCross indices */
  264. {
  265.   double minZZ,                 /* minimum zero-crossing length for peak */
  266.     lengthZZ;                   /* zero-crossing length */
  267.   long firstCross,              /* first, second, and last zero crossing */
  268.     secondCross, lastCross, featFlag,  /* =0 if no features in cycle; 1 otherwise */
  269.     i;
  270.  
  271.   minZZ = (double) (wwParL - 1) + EPSILON1;
  272.  
  273.   for (firstCross = *zInitial + 1; firstCross <= *zFinal; firstCross++)
  274.     if ((*zCross)[firstCross] != NOCROSSING)
  275.       break;
  276.   for (lastCross = *zFinal; lastCross >= *zInitial; --lastCross)
  277.     if ((*zCross)[lastCross] != NOCROSSING)
  278.       break;
  279.  
  280. /* if it is not a cycle, clean up zero-crossings at bounds */
  281. /*    if (cycleFlag == 0) { ??? */
  282.   if (1) {
  283.     featFlag = 1;
  284.  
  285.     /* eliminate peaks at beginning and end if too narrow */
  286.     lengthZZ = theta[firstCross + 1].length;
  287.     if (lengthZZ < minZZ) {
  288.       (*zCross)[*zInitial] = NOCROSSING;
  289.       *zInitial = firstCross;
  290.     }
  291.     lengthZZ = theta[*nZCross - 1].length - theta[lastCross - 1].length;
  292.     if (lengthZZ < minZZ)
  293.       (*zCross)[lastCross] = NOCROSSING;
  294.  
  295.     /* move the initial zero crossing to the first element and put a
  296.      * ZERO crossing at the final element if it ended with
  297.      * an UP/DOWN; this is for boundary feature determination */
  298.     (*zCross)[1] = (*zCross)[*zInitial];
  299.     (*zCross)[*zInitial] = NOCROSSING;
  300.     *zInitial = 1;
  301.     for (i = lastCross; ((*zCross)[i] == NOCROSSING && i >= *zInitial);
  302.          --i);
  303.     *zFinal = *nZCross - 1;
  304.     if ((*zCross)[i] != ZERO)
  305.       (*zCross)[*zFinal] = ZERO;
  306.   }
  307. /* if it is a cycle, clean up zero-crossings at bounds (= seam) */
  308.   else {
  309.     if ((*zCross)[*zInitial] != NOCROSSING)
  310.       firstCross = *zInitial;
  311.  
  312.     /* if the first and last zero-crossings are the same, remove first */
  313.     if (firstCross <= *zFinal
  314.         && (*zCross)[firstCross] == (*zCross)[lastCross]) {
  315.       (*zCross)[firstCross] = NOCROSSING;
  316.       for (firstCross = firstCross + 1; firstCross <= *zFinal;
  317.            firstCross++)
  318.         if ((*zCross)[firstCross] != NOCROSSING)
  319.           break;
  320.     }
  321.     /* if first and last zero crossing don't have zero between, insert it */
  322.     else if (firstCross <= *zFinal && (*zCross)[firstCross] != ZERO
  323.              && (*zCross)[lastCross] != ZERO) {
  324.       (*zCross)[firstCross + 1] = (*zCross)[firstCross];
  325.       (*zCross)[firstCross] = ZERO;
  326.     }
  327.     /* if a peak at the seam is too narrow, eliminate it */
  328.     for (secondCross = firstCross + 1; secondCross <= *zFinal;
  329.          secondCross++)
  330.       if ((*zCross)[secondCross] != NOCROSSING)
  331.         break;
  332.     if (secondCross < lastCross
  333.         && (*zCross)[secondCross] == (*zCross)[lastCross]) {
  334.       lengthZZ = theta[*zFinal].length - theta[lastCross - 1].length
  335.         + theta[firstCross + 1].length
  336.         - theta[*zInitial].length;
  337.       if (lengthZZ < minZZ) {
  338.         (*zCross)[firstCross] = NOCROSSING;
  339.         (*zCross)[secondCross] = NOCROSSING;
  340.         for (firstCross = secondCross + 1; firstCross <= *zFinal;
  341.              firstCross++)
  342.           if ((*zCross)[firstCross] != NOCROSSING)
  343.             break;
  344.         for (secondCross = firstCross + 1; secondCross <= *zFinal;
  345.              secondCross++)
  346.           if ((*zCross)[secondCross] != NOCROSSING)
  347.             break;
  348.       }
  349.     }
  350.   }
  351.   return (featFlag);
  352. }
  353.  
  354.  
  355.  
  356. /* FEATDETERM:  function determines the feature type and parameters
  357.  *            from the theta plot and other information
  358.  *                 usage: featdeterm (zCross, zStart, zEnd, theta, nTheta,
  359.  *                              data, edge, wwPar, &zLow, &zHigh, curveFlag)
  360.  */
  361.  
  362. long
  363. featdeterm (zCross, zStart, zEnd, theta, nTheta, data, edge, wwPar,
  364.             zLow, zHigh, curveFlag)
  365.      short *zCross;             /* array of flags at crossings on theta */
  366.      long zStart,               /* start and end zero crossing indices */
  367.        zEnd;                    /* (these are *within* peak) */
  368.      struct theta *theta;       /* theta plot */
  369.      long nTheta;               /* no. values in theta plot */
  370.      struct point *data;        /* data array */
  371.      struct edge edge;          /* lower and upper <data> indices */
  372.      struct wwPar wwPar;        /* wing-walk parameters */
  373.      long *zLow,                /* line indices of end of peak before,  */
  374.       *zHigh;                   /*   and start of peak after, feature peak */
  375.      short curveFlag;           /* =0 for curve fit by 2 lines; =1 for curve */
  376. {
  377.   long lD2,                     /* WMW-length divided by 2 */
  378.     i;
  379.  
  380.   long zLowB,                   /* line indices of peak before and after   */
  381.     zHighB;                     /*   including buffer of half ww length */
  382.  
  383.   double widthZZ,               /* peak width between zero crossings */
  384.     threshZZ;                   /* zero crossing width threshold between
  385.                                  * * curve and corner */
  386.   struct dpoint intrsct;        /* intersection of tangents to curve */
  387.   struct point corner,          /* corner coordinates */
  388.     trans1,                     /* transition coords into and out of curve */
  389.     trans2;
  390.   struct dpoint center;         /* center of curvature */
  391.   double radiusC;               /* radius of curvature */
  392.   short curveSweep;             /* indicates sweep of curve */
  393.   struct dpoint midArcPt;       /* middle coord on arc */
  394.   short curveDirn;
  395.  
  396.   widthZZ = theta[zEnd + 1].length - theta[zStart - 1].length;
  397.   threshZZ = (double) (wwPar.l + maxCorner - 1) + EPSILON2;
  398.  
  399. /* find zero crossings before and after peak */
  400.   for (i = zStart - 1; i >= 0; --i) {
  401.     if (zCross[i] == ZERO)
  402.       break;
  403.     else if (zCross[i] == UP || zCross[i] == DOWN)
  404.       printf ("FEATDETERM error: UP/DOWN crossing instead of ZERO\n");
  405.   }
  406.   *zLow = i;
  407.   for (i = zEnd + 2; i < nTheta; i++) {
  408.     if (zCross[i] == UP || zCross[i] == DOWN)
  409.       break;
  410.     else if (zCross[i] == ZERO)
  411.       printf ("FEATDETERM error: ZERO crossing instead of UP/DOWN\n");
  412.   }
  413.   *zHigh = i - 1;
  414.  
  415.   lD2 = wwPar.l / 2;            /* corrected from (L - M) / 2  6Jan88 */
  416.   zLowB = ((*zLow - lD2) < 0) ? 0 : (*zLow - lD2);
  417.   zStart = zStart - 1 + lD2;
  418.   zEnd = zEnd + 1 - lD2;
  419.   zHighB = ((*zHigh + lD2) >= nTheta) ? nTheta - 1 : (*zHigh + lD2);
  420.   if (zEnd < zStart)
  421.     zStart = zEnd = (zStart + zEnd) / 2;
  422.  
  423.   if (widthZZ <= threshZZ || zEnd <= zStart) {  /* 1Nov88, 2nd cond. added */
  424.     featcorner (data, zLowB, zStart, zEnd, zHighB, edge, &corner);
  425.     linkcorner (WWCORNER, corner);
  426.   }
  427.   else {
  428.     if ((featcurve (data, zLowB, zStart, zEnd, zHighB, edge, widthZZ,
  429.                     threshZZ, &intrsct, &trans1, &trans2, ¢er, &radiusC,
  430.                     &curveSweep)) == 0) {
  431.  
  432. /* the mid arc point are not reliable after <eventranspts> function
  433.  * if the sweep is <=180; therefore if this is so, use intersection pt;
  434.  * if the sweep is >180 then <eventranspts> shortens the longer length
  435.  * and using mid arc point is okay (13 Sept 88) */
  436.       if (curveSweep == CURVELT180)
  437.         midArcPt = intrsct;
  438.       else if (curveSweep == CURVEEQ180) {
  439.         midArcPt.x = (trans1.x + trans2.x) / 2;
  440.         midArcPt.y = (trans1.y + trans2.y) / 2;
  441.       }
  442.       else {
  443.         midArcPt.x = data[(zEnd + zStart) / 2].x;
  444.         midArcPt.y = data[(zEnd + zStart) / 2].y;
  445.       }
  446.       curveDirn = curvedirn (trans1, trans2, midArcPt, center);
  447.       if (curveFlag == 0) {
  448.         center.x = midArcPt.x;
  449.         center.y = midArcPt.y;
  450.       }
  451.       linkcurve (intrsct, trans1, trans2, center, radiusC, curveDirn);
  452.     }
  453.     else {                      /* asymmetric curve = shorter curve converted to corner */
  454.       corner.x = (long) (center.x + 0.5);
  455.       corner.y = (long) (center.y + 0.5);
  456.       linkcorner (WWCORNER, corner);
  457.     }
  458.   }
  459.   return (0);
  460. }
  461.  
  462.  
  463.  
  464. /* FEATCORNER:  function calculates corner parameters
  465.  *                    usage: featcorner (data, zLow, zStart, zEnd, zHigh,
  466.  *                                               edge, &coord)
  467.  *              Note that in this function the <zLow, zStart,...> 
  468.  *              variables have already been adjusted +- (L-1)/2 
  469.  *              and that they are indices of the <zCross> or <theta>
  470.  *              array, not of the data array.
  471.  */
  472.  
  473.  
  474. long
  475. featcorner (data, zLow, zStart, zEnd, zHigh, edge, coord)
  476.      struct point *data;        /* data array of coordinates */
  477.      long zLow,                 /* indices of line coords before peak */
  478.        zStart, zEnd,            /* indices of line coords after peak */
  479.        zHigh;
  480.      struct edge edge;          /* lower and upper <data> indices */
  481.      struct point *coord;       /* output location of corner */
  482. {
  483.   double d, e;                  /* temporary variables */
  484.   double deltaX1, deltaX2;      /* x difference along fore and aft lines */
  485.   double coordX, coordY;        /* floating point corner coordinates */
  486.   short parallel;               /* flag = 1 if lines parallel; else 0 */
  487.  
  488.   zLow += edge.iLow;
  489.   zStart += edge.iLow;
  490.   zEnd += edge.iLow;
  491.   zHigh += edge.iLow;
  492.  
  493.   deltaX1 = (double) (data[zStart].x - data[zLow].x);
  494.   deltaX2 = (double) (data[zHigh].x - data[zEnd].x);
  495.   parallel = 0;
  496.   if (deltaX1 == 0.0 && deltaX2 == 0.0)
  497.     parallel = 1;
  498.   else if (deltaX1 == 0.0) {
  499.     e = (double) (data[zHigh].y - data[zEnd].y) / deltaX2;
  500.     coord->x = data[zStart].x;
  501.     coord->y = (long) (e * (coord->x - (double) data[zEnd].x) + (double) data[zEnd].y);
  502.   }
  503.   else if (deltaX2 == 0.0) {
  504.     d = (double) (data[zStart].y - data[zLow].y) / deltaX1;
  505.     coord->x = data[zHigh].x;
  506.     coord->y = (long) (d * (coord->x - (double) data[zLow].x) + (double) data[zLow].y);
  507.   }
  508.   else {
  509.     d = (double) (data[zStart].y - data[zLow].y) / deltaX1;
  510.     e = (double) (data[zHigh].y - data[zEnd].y) / deltaX2;
  511.     if (d == e)
  512.       parallel = 1;
  513.     else {
  514.       coordX = (d * data[zLow].x - e * data[zEnd].x + data[zEnd].y
  515.                 - data[zLow].y) / (d - e);
  516.       coordY = d * (coordX - data[zLow].x) + data[zLow].y;
  517.       coord->x = (int) (coordX + 0.5);
  518.       coord->y = (int) (coordY + 0.5);
  519.     }
  520.   }
  521.   if (parallel == 1) {
  522.     coord->x = (data[zStart].x + data[zEnd].x) / 2;
  523.     coord->y = (data[zStart].y + data[zEnd].y) / 2;
  524.   }
  525.   return (0);
  526. }
  527.  
  528.  
  529.  
  530. /* FEATCURVE:   function calculates curve parameters
  531.  *                usage: cornerOrCurve = 
  532.  *                      featcurve (data, zLow, zStart, zEnd, zHigh, edge,
  533.  *                                       widthZZ, threshZZ, 
  534.  *                                      &intrsct, &trans1, &trans2,
  535.  *                                           ¢er, &radiusC, &curveSweep)
  536.  *            - Function returns a flag which is 0 if the feature is a curve
  537.  *              but -1 if the curve has been converted to a corner (see next).
  538.  *            - Note that if the curve is asymmetric, the further adjacent
  539.  *              line is extended to the same distance from the corner as
  540.  *              the other line. If the new zero crossing width after this
  541.  *              is less than the corner/curve threshold, then the feature
  542.  *              is a corner.
  543.  *            - Note that in this function the <zLow, zStart,...> 
  544.  *              variables have already been adjusted +- (L-1)/2 
  545.  *              and that they are indices of the <zCross> or <theta>
  546.  *              array, not of the data array.
  547.  *            - Note that the case where the transition slopes are equal
  548.  *              or they do not intersect are not taken care of here.
  549.  */
  550.  
  551. #include <math.h>
  552.  
  553. /*#define PI 3.14159265358979 */
  554.  
  555. long
  556. featcurve (data, zLow, zStart, zEnd, zHigh, edge, widthZZ, threshZZ,
  557.            intrsct, trans1, trans2, center, radiusC, curveSweep)
  558.      struct point *data;        /* data array of coordinates */
  559.      long zLow,                 /* indices of line coords before peak */
  560.        zStart, zEnd,            /* indices of line coords after peak */
  561.        zHigh;
  562.      struct edge edge;          /* lower and upper <data> indices */
  563.      double widthZZ,            /* peak width between zero crossings */
  564.        threshZZ;                /* corner/curve thresh zero crossing width */
  565.      struct dpoint *intrsct;    /* intersection of lines tangent to curve */
  566.      struct point *trans1,      /* transition coords into and out of curve */
  567.       *trans2;
  568.      struct dpoint *center;     /* center of curvature */
  569.      double *radiusC;           /* radius of curvature */
  570.      short *curveSweep;         /* indicates sweep of curve */
  571. {
  572.   double deltaX1, deltaX2,      /* x, y diff.s along fore and aft lines */
  573.     deltaY1, deltaY2, d, e,     /* used in equation for center */
  574.     distTI, distEI;             /* distances: trans and end to intersectn */
  575.   struct dpoint transD1,        /* (double) of transition point locations */
  576.     transD2;
  577.  
  578. /* add offset to curve locations to correspond to <data> array */
  579.   zLow += edge.iLow;
  580.   zStart += edge.iLow;
  581.   zEnd += edge.iLow;
  582.   zHigh += edge.iLow;
  583.  
  584. /* find corner of intersection of tangent lines to curve */
  585.   *curveSweep = CURVELT180;
  586.   deltaX1 = (double) (data[zStart].x - data[zLow].x);
  587.   deltaX2 = (double) (data[zHigh].x - data[zEnd].x);
  588.   deltaY1 = (double) (data[zStart].y - data[zLow].y);
  589.   deltaY2 = (double) (data[zHigh].y - data[zEnd].y);
  590.   if (deltaX1 == 0.0 && deltaX2 == 0.0)
  591.     *curveSweep = CURVEEQ180;
  592.   else if (deltaX1 == 0.0) {
  593.     e = deltaY2 / deltaX2;
  594.     intrsct->x = (double) data[zStart].x;
  595.     intrsct->y = e * (intrsct->x - data[zEnd].x) + data[zEnd].y;
  596.   }
  597.   else if (deltaX2 == 0.0) {
  598.     d = deltaY1 / deltaX1;
  599.     intrsct->x = (double) data[zHigh].x;
  600.     intrsct->y = d * (intrsct->x - data[zLow].x) + data[zLow].y;
  601.   }
  602.   else {
  603.     d = deltaY1 / deltaX1;
  604.     e = deltaY2 / deltaX2;
  605.     if (d == e)
  606.       *curveSweep = CURVEEQ180;
  607.     else {
  608.       intrsct->x = (d * data[zLow].x - e * data[zEnd].x + data[zEnd].y
  609.                     - data[zLow].y) / (d - e);
  610.       intrsct->y = d * (intrsct->x - data[zLow].x) + data[zLow].y;
  611.     }
  612.   }
  613.  
  614. /* find curve sweep, for curves which sweep greater than or less than 180 */
  615.   if (*curveSweep != CURVEEQ180) {
  616.     distTI = (intrsct->x - (double) data[zStart].x)
  617.       * (intrsct->x - (double) data[zStart].x)
  618.       + (intrsct->y - (double) data[zStart].y)
  619.       * (intrsct->y - (double) data[zStart].y);
  620.     distEI = (intrsct->x - (double) data[zLow].x)
  621.       * (intrsct->x - (double) data[zLow].x)
  622.       + (intrsct->y - (double) data[zLow].y)
  623.       * (intrsct->y - (double) data[zLow].y);
  624.     *curveSweep = (distEI > distTI) ? CURVELT180 : CURVEGT180;
  625.   }
  626.  
  627. /* adjust transition points so they are symmetric around middle of curve */
  628.   transD1.x = (double) data[zStart].x;
  629.   transD1.y = (double) data[zStart].y;
  630.   transD2.x = (double) data[zEnd].x;
  631.   transD2.y = (double) data[zEnd].y;
  632.   eventranspts (&transD1, &transD2, data[zLow], data[zHigh], *intrsct,
  633.                 &widthZZ);
  634.  
  635.   trans1->x = (transD1.x >= 0.0) ? (int) (transD1.x + 0.5) :
  636.     (int) (transD1.x - 0.5);
  637.   trans1->y = (transD1.y >= 0.0) ? (int) (transD1.y + 0.5) :
  638.     (int) (transD1.y - 0.5);
  639.   trans2->x = (transD2.x >= 0.0) ? (int) (transD2.x + 0.5) :
  640.     (int) (transD2.x - 0.5);
  641.   trans2->y = (transD2.y >= 0.0) ? (int) (transD2.y + 0.5) :
  642.     (int) (transD2.y - 0.5);
  643.  
  644. /* if the new zero-crossing  width is less than the corner/curve threshold,
  645.  * then return a corner */
  646.   if (widthZZ <= threshZZ) {
  647.     center->x = intrsct->x;
  648.     center->y = intrsct->y;
  649.     return (-1);
  650.   }
  651.  
  652. /* center of curvature */
  653.   if (deltaX1 == deltaX2);
  654. /*      printf ("FEATCURVE: ERROR:problems -curve >= 180 degrees in sweep\n"); */
  655.   else if (deltaX1 == 0.0) {
  656.     e = (double) (data[zHigh].y - data[zEnd].y) / deltaX2;
  657.     center->y = transD1.y;
  658.     center->x = transD2.x - e * (center->y - transD2.y);
  659.   }
  660.   else if (deltaX2 == 0.0) {
  661.     d = (double) (data[zStart].y - data[zLow].y) / deltaX1;
  662.     center->y = transD2.y;
  663.     center->x = transD1.x - d * (center->y - transD1.y);
  664.   }
  665.   else {
  666.     d = (double) (data[zStart].y - data[zLow].y) / deltaX1;
  667.     e = (double) (data[zHigh].y - data[zEnd].y) / deltaX2;
  668.     if (d == e);
  669. /*          printf ("FEATCURVE: transition line slopes equal\n"); */
  670.     else if (d == 0) {
  671.       center->x = transD1.x;
  672.       center->y = transD2.y - (center->x - transD2.x) / e;
  673.     }
  674.     else {
  675.       center->x = (d * e * transD2.y + d * transD2.x
  676.                    - d * e * transD1.y - e * transD1.x) / (d - e);
  677.       center->y = (transD1.x - center->x) / d + transD1.y;
  678.     }
  679.   }
  680.  
  681. /* radius of curvature */
  682.   *radiusC = sqrt ((center->x - transD1.x) * (center->x - transD1.x)
  683.                    + (center->y - transD1.y) * (center->y - transD1.y));
  684.  
  685.   return (0);
  686. }
  687.  
  688. /* EVENTRANSPTS:        function evens the transition points so they are
  689.  *                    symmetric around the middle of curve
  690.  *                              usage: eventranspts (&trans1, &trans2,
  691.  *                                              lineLow, lineHigh, intrsct,
  692.  *                                              &widthZZ)
  693.  *                      Note that the points are evened here so that
  694.  *                      the shorter is elongated so that its end is
  695.  *                      perpendicular to the longer. Other valid adjustments
  696.  *                      would be to shorten the longer to the shorter,
  697.  *                      or to adjust them both to the middle length.    
  698.  */
  699.  
  700. #define NOTEQUALSLOPES 0        /* 2 lines of different slopes */
  701. #define EQUALSLOPES 1           /* 2 lines of same slopes (not 0, infin.) */
  702. #define INFINITESLOPES 2        /* 2 lines of same infinite slopes */
  703. #define ZEROSLOPES 3            /* 2 lines of same 0 slopes */
  704.  
  705. long
  706. eventranspts (trans1, trans2, lineLow, lineHigh, intrsct, widthZZ)
  707.      struct dpoint *trans1, *trans2;  /* transition points to be evened */
  708.      struct point lineLow, lineHigh;  /* low/high points of transition lines */
  709.      struct dpoint intrsct;     /* intersection point of transition lines */
  710.      double *widthZZ;           /* zero crossing width */
  711. {
  712.   double deltaX1, deltaY1,      /* x,y increments of lines */
  713.     deltaX2, deltaY2, mPerpen, bPerpen;  /* slope, intercept for perpendicular line */
  714.   short equalSlopes;            /* flag showing relative slopes for lines */
  715.   long furtherPt;               /* extended point for parallel lines */
  716.   struct point pt;              /* other pt on perpendicular line */
  717.   struct dpoint intrsct1, intrsct2;  /* intersects with perpendicular */
  718.   double dist1, dist2, fraction;
  719.   struct point lineLow0;
  720.  
  721.   deltaX1 = trans1->x - (double) lineLow.x;
  722.   deltaX2 = (double) lineHigh.x - trans2->x;
  723.   deltaY1 = trans1->y - (double) lineLow.y;
  724.   deltaY2 = (double) lineHigh.y - trans2->y;
  725.  
  726.   if (deltaX1 == 0.0 && deltaX2 == 0.0)
  727.     equalSlopes = INFINITESLOPES;
  728.   else if (deltaY1 == 0.0 && deltaY2 == 0.0)
  729.     equalSlopes = ZEROSLOPES;
  730.   else if (deltaY1 / deltaX1 == deltaY2 / deltaX2)
  731.     equalSlopes = EQUALSLOPES;
  732.   else
  733.     equalSlopes = NOTEQUALSLOPES;
  734.  
  735.   switch (equalSlopes) {
  736.     /* adjust transition point for lines of different slopes */
  737.     /* this adjusts the further transition pt to be the same distance
  738.      * to the intersect for arcs of BOTH CURVELT180 and CURVEGT180 */
  739.   case NOTEQUALSLOPES:
  740.     dist1 = (trans1->x - intrsct.x) * (trans1->x - intrsct.x)
  741.       + (trans1->y - intrsct.y) * (trans1->y - intrsct.y);
  742.     dist2 = (trans2->x - intrsct.x) * (trans2->x - intrsct.x)
  743.       + (trans2->y - intrsct.y) * (trans2->y - intrsct.y);
  744.     if (dist1 > dist2) {
  745.       dist1 = sqrt (dist1);
  746.       dist2 = sqrt (dist2);
  747.       fraction = (dist1 - dist2) / dist1;
  748.       trans1->x = trans1->x + fraction * (intrsct.x - trans1->x);
  749.       trans1->y = trans1->y + fraction * (intrsct.y - trans1->y);
  750.       *widthZZ -= dist1 - dist2;
  751.     }
  752.     else if (dist2 > dist1) {
  753.       dist1 = sqrt (dist1);
  754.       dist2 = sqrt (dist2);
  755.       fraction = (dist2 - dist1) / dist2;
  756.       trans2->x = trans2->x + fraction * (intrsct.x - trans2->x);
  757.       trans2->y = trans2->y + fraction * (intrsct.y - trans2->y);
  758.       *widthZZ -= dist2 - dist1;
  759.     }
  760.     break;
  761.     /* adjust transition point for both lines of vertical slope */
  762.   case INFINITESLOPES:
  763.     if (trans1->y >= lineLow.y)
  764.       furtherPt = (trans1->y >= trans2->y) ? (long) trans1->y : (long) trans2->y;
  765.     else
  766.       furtherPt = (trans1->y <= trans2->y) ? (long) trans1->y : (long) trans2->y;
  767.     trans1->y = trans2->y = (double) furtherPt;
  768.     break;
  769.     /* adjust transition point for both lines of zero slope */
  770.   case ZEROSLOPES:
  771.     if (trans1->x >= lineLow.x)
  772.       furtherPt = (trans1->x >= trans2->x) ? (long) trans1->x : (long) trans2->x;
  773.     else
  774.       furtherPt = (trans1->x <= trans2->x) ? (long) trans1->x : (long) trans2->x;
  775.     trans1->x = trans2->x = (double) furtherPt;
  776.     break;
  777.     /* adjust transition point for both lines parallel */
  778.   case EQUALSLOPES:
  779.     mPerpen = -deltaX1 / deltaY1;
  780.     bPerpen = trans1->y - mPerpen * trans1->x;
  781.     pt.x = (long) trans1->x + 10000;  /* 10000 is arbitrary, but need big
  782.                                        * * number so int roundoff small */
  783.     pt.y = (long) (mPerpen * pt.x + bPerpen);
  784.     lineLow0.x = (long) trans2->x;
  785.     lineLow0.y = (long) trans2->y;
  786.     linesintrsct (pt, *trans1, lineLow0, lineHigh, &intrsct1);
  787.  
  788.     bPerpen = trans2->y - mPerpen * trans2->x;
  789.     pt.x = (long) trans2->x + 10000;  /* 10000 is arbitrary, but need big
  790.                                        * * number so int roundoff small */
  791.     pt.y = (long) (mPerpen * pt.x + bPerpen);
  792.     lineHigh.x = (long) trans1->x;
  793.     lineHigh.y = (long) trans1->y;
  794.     linesintrsct (pt, *trans2, lineLow, lineHigh, &intrsct2);
  795.     dist1 = (trans1->x - lineLow.x) * (trans1->x - lineLow.x)
  796.       + (trans1->y - lineLow.y) * (trans1->y - lineLow.y);
  797.     dist2 = (intrsct1.x - lineLow.x) * (intrsct1.x - lineLow.x)
  798.       + (intrsct1.y - lineLow.y) * (intrsct1.y - lineLow.y);
  799.     if (dist2 > dist1) {
  800.       trans2->x = intrsct1.x;
  801.       trans2->y = intrsct1.y;
  802.     }
  803.     else {
  804.       trans1->x = intrsct2.x;
  805.       trans1->y = intrsct2.y;
  806.     }
  807.     break;
  808.   }
  809.   return (0);
  810. }
  811.  
  812.  
  813. /* CYCLEAUGMENT: function rewrites the zero-crossing array (and corresponding
  814.  *             data array) so that the array begins at the first zero
  815.  *               crossing for the entire array, then the array is wrapped
  816.  *               around the beginning until the first UP/DOWN crossing
  817.  *               so that features can be determined in the FEATDETERM 
  818.  *               function
  819.  *                  usage: nFeatsFlag = cycleaugment (wwPar.l, 
  820.  *                                         &zCross, &nZCross, &data, &edge,
  821.  *                                                      zInitial, zFinal)
  822.  *           -  The flag <nFeatsFlag> is 0 if no crossings are found,
  823.  *              and 1 if there is at least one feature.
  824.  *           -  Note that the <zCross> plot has already been overlapped
  825.  *              by <zInitial> at the front and <zFinal> at the end.
  826.  *              Therefore the additional wrap-around here is done from
  827.  *              within this overlap.
  828.  */
  829.  
  830. short
  831. cycleaugment (wwParL, zCross, nZCross, data, edge, zInitial, zFinal)
  832.      long wwParL;               /* ww length parameter */
  833.      short **zCross;            /* array of flags at zero-crossing locns */
  834.      long *nZCross;             /* no. values in zero-crossing plot */
  835.      struct point **data;       /* sequence of coordinates of lines */
  836.      struct edge *edge;         /* lower and upper <data> indices */
  837.      long zInitial, zFinal;     /* initial and final valid zCross indices */
  838. {
  839.   long nNew;                    /* no. new points */
  840.   int firstZero, firstUpDown,   /* indices of first ZERO and first UP/DOWN */
  841.   /*    crossings on the zero-crossing plot */
  842.     i, j;
  843.   short *zCrossNew;             /* new zero-crossing plot */
  844.   struct point *dataNew;        /* new data array containing 
  845.                                  * * wrap-around cycle */
  846.  
  847. /* look for first ZERO crossing to start from */
  848.   for (firstZero = zInitial; firstZero <= zFinal; firstZero++)
  849.     if ((*zCross)[firstZero] == ZERO)
  850.       break;
  851.  
  852.   if (firstZero >= zFinal)
  853.     return (0);                 /* no ZEROs => circle */
  854.  
  855. /* look for first UP/DOWN crossing after initial point to end with */
  856.   for (firstUpDown = firstZero + 1; firstUpDown <= zFinal; firstUpDown++)
  857.     if ((*zCross)[firstUpDown] == UP
  858.         || (*zCross)[firstUpDown] == DOWN)
  859.       break;
  860.   if (firstUpDown > zFinal) {
  861.     for (firstUpDown = zInitial; firstUpDown < firstZero; firstUpDown++)
  862.       if ((*zCross)[firstUpDown] == UP
  863.           || (*zCross)[firstUpDown] == DOWN)
  864.         break;
  865.     if (firstUpDown >= firstZero)
  866.       return (0);               /* 1 crossing => circle */
  867.   }
  868.  
  869. /* allocate space and write zero-crossing plot augmented with wrap-around */
  870.   nNew = (zFinal - firstZero + 1) + (firstUpDown - zInitial + 1);
  871.   if (firstUpDown < firstZero)
  872.     nNew += zFinal - zInitial + 1;
  873.  
  874.   if ((zCrossNew = (short *) calloc (nNew, sizeof (*zCrossNew))) == NULL) {
  875.     printf ("CYCLEAUGMENT: CALLOC: not enough memory -- sorry");
  876.     return (-2);
  877.   }
  878.   if ((dataNew = (struct point *) calloc (nNew, sizeof (*dataNew))) == NULL) {
  879.     printf ("CYCLEAUGMENT: CALLOC not enough memory");
  880.     return (-3);
  881.   }
  882.  
  883.   for (i = firstZero, j = 0; i <= zFinal; i++) {
  884.     zCrossNew[j] = (*zCross)[i];
  885.     dataNew[j++] = (*data)[i];
  886.   }
  887.   for (i = zInitial; i < firstZero; i++) {
  888.     zCrossNew[j] = (*zCross)[i];
  889.     dataNew[j++] = (*data)[i];
  890.   }
  891.   if (firstUpDown > firstZero)
  892.     for (i = firstZero; i <= firstUpDown; i++) {
  893.       zCrossNew[j] = (*zCross)[i];
  894.       dataNew[j++] = (*data)[i];
  895.     }
  896.   else {
  897.     for (i = firstZero; i <= zFinal; i++) {
  898.       zCrossNew[j] = (*zCross)[i];
  899.       dataNew[j++] = (*data)[i];
  900.     }
  901.     for (i = zInitial; i <= firstUpDown; i++) {
  902.       zCrossNew[j] = (*zCross)[i];
  903.       dataNew[j++] = (*data)[i];
  904.     }
  905.   }
  906.   if (j != nNew)
  907.     printf ("CYCLEAUGMENT: oh nuts: j = %3d, nNew = %3d\n",
  908.             j, nNew);
  909.   free (*zCross);
  910.   (*zCross) = zCrossNew;
  911.   (*data) = dataNew;
  912.   (*edge).iHigh = nNew - 1;
  913.   *nZCross = nNew;
  914.  
  915.   return (1);
  916. }
  917.  
  918.  
  919.  
  920. /* BOUNDDETERM2: function determines boundary feature
  921.  *             former function, BOUNDDETERM, fit a curve as the boundary
  922.  *               feature; here, we fit two corners
  923.  *                 usage: bounddeterm2 (zCross, theta, zStart, zEnd, nTheta,
  924.  *                                       edge, data, &zLow, &zHigh)
  925.  */
  926.  
  927. #define NOBOUNDCURVE -1         /* flag indicates small curvature */
  928.  
  929. long
  930. bounddeterm2 (zCross, theta, zStart, zEnd, nTheta, data, edge, zLow, zHigh)
  931.      short *zCross;             /* array of flags at crossings on theta */
  932.      struct theta *theta;       /* theta plot */
  933.      long zStart,               /* start and end zero crossing indices */
  934.        zEnd;                    /* (these are *within* peak) */
  935.      long nTheta;               /* no. values in theta plot */
  936.      struct point *data;        /* data array */
  937.      struct edge edge;          /* lower and upper <data> indices */
  938.      long *zLow,                /* line indices of end of peak before,  */
  939.       *zHigh;                   /*   and start of peak after, feature peak */
  940. {
  941.   struct point linePt,          /* pt on line, not transition pt */
  942.     transPt,                    /* transition pt between curve/line */
  943.     endPt;                      /* final pt on curve */
  944.   struct dpoint midArcPt;       /* middle coord on arc */
  945.   struct point corner,          /* transition coords into and out of curve */
  946.     trans1, trans2;
  947.   int i;
  948.  
  949. /* for feature at beginning of line (or beginning and end), determine coords */
  950.  
  951.   *zLow += edge.iLow;
  952.   zStart += edge.iLow;
  953.   zEnd += edge.iLow;
  954.   *zHigh += edge.iLow;
  955.  
  956.   if (zStart == 1) {
  957.     *zLow = zStart;
  958.     endPt = data[zStart];
  959.     transPt = data[zEnd];
  960.     /* find zero crossing after peak */
  961.     for (i = zEnd + 2 - edge.iLow; i < nTheta; i++) {
  962.       if (zCross[i] == UP || zCross[i] == DOWN)
  963.         break;
  964.       else if (zCross[i] == ZERO)
  965.         printf
  966.           ("BOUNDDETERM error: ZERO crossing instead of UP/DOWN\n");
  967.     }
  968.     *zHigh = i - 1;
  969.     linePt = data[*zHigh];
  970.     trans1 = endPt;
  971.     trans2 = transPt;
  972.     corner.x = trans1.x;
  973.     corner.y = trans1.y;
  974.     linkcorner (WWCORNER, corner);
  975.  
  976.     midArcPt.x = data[(zEnd + zStart) / 2].x;
  977.     midArcPt.y = data[(zEnd + zStart) / 2].y;
  978.     corner.x = (long) (midArcPt.x + 0.5);
  979.     corner.y = (long) (midArcPt.y + 0.5);
  980.     corner.x = trans2.x;
  981.     corner.y = trans2.y;
  982.     linkcorner (WWCORNER, corner);
  983.   }
  984.  
  985. /* for a feature at end of line, determine coordinates */
  986.   else {
  987.     *zHigh = zEnd;
  988.     for (i = zStart - 1 - edge.iLow; i >= 0; --i) {
  989.       if (zCross[i] == ZERO)
  990.         break;
  991.       else if (zCross[i] == UP || zCross[i] == DOWN)
  992.         printf
  993.           ("BOUNDDETERM error: UP/DOWN crossing instead of ZERO\n");
  994.     }
  995.     *zLow = i;
  996.     linePt = data[*zLow];
  997.     transPt = data[zStart];
  998.     endPt = data[zEnd + 1];
  999.     trans1 = transPt;
  1000.     trans2 = endPt;
  1001.  
  1002.   }
  1003.  
  1004. /* determine feature parameters, and store them */
  1005.  
  1006.   return (0);
  1007. }
  1008.  
  1009.  
  1010.  
  1011.  
  1012. /* BOUNDDETERM: function determines boundary feature
  1013.  *               usage: bounddeterm (zCross, theta, zStart, zEnd, nTheta,
  1014.  *                                       data, &zLow, &zHigh)
  1015.  */
  1016.  
  1017. #define NOBOUNDCURVE -1         /* flag indicates small curvature */
  1018.  
  1019. long
  1020. bounddeterm (zCross, theta, zStart, zEnd, nTheta, data, zLow, zHigh)
  1021.      short *zCross;             /* array of flags at crossings on theta */
  1022.      struct theta *theta;       /* theta plot */
  1023.      long zStart,               /* start and end zero crossing indices */
  1024.        zEnd;                    /* (these are *within* peak) */
  1025.      long nTheta;               /* no. values in theta plot */
  1026.      struct point *data;        /* data array */
  1027.      long *zLow,                /* line indices of end of peak before,  */
  1028.       *zHigh;                   /*   and start of peak after, feature peak */
  1029. {
  1030.   struct point linePt,          /* pt on line, not transition pt */
  1031.     transPt,                    /* transition pt between curve/line */
  1032.     endPt;                      /* final pt on curve */
  1033.   struct dpoint intrsct;        /* intersection of tangents to curve */
  1034.   double arcLength;             /* length of arc of curve */
  1035.   short curveSweep;             /* indicates sweep of curve */
  1036.   struct dpoint midArcPt;       /* middle coord on arc */
  1037.   struct point trans1,          /* transition coords into and out of curve */
  1038.     trans2;
  1039.   struct dpoint center;         /* center of curvature */
  1040.   double radiusC;               /* radius of curvature */
  1041.   short curveDirn;
  1042.   int i;
  1043.  
  1044. /* for curve at beginning of line (or beginning and end), determine coords */
  1045.   if (zStart == 1) {
  1046.     *zLow = zStart;
  1047.     endPt = data[zStart];       /* 8Jan88 -- extra pt; 7Feb97 reverted back */
  1048.     transPt = data[zEnd];
  1049.     /* find zero crossing after peak */
  1050.     for (i = zEnd + 2; i < nTheta; i++) {
  1051.       if (zCross[i] == UP || zCross[i] == DOWN)
  1052.         break;
  1053.       else if (zCross[i] == ZERO)
  1054.         printf
  1055.           ("BOUNDDETERM error: ZERO crossing instead of UP/DOWN\n");
  1056.     }
  1057.     *zHigh = i - 1;
  1058.     linePt = data[*zHigh];
  1059.     trans1 = endPt;
  1060.     trans2 = transPt;
  1061.   }
  1062.  
  1063. /* for a curve at end of line, determine coordinates */
  1064.   else {
  1065.     *zHigh = zEnd;
  1066.     for (i = zStart - 1; i >= 0; --i) {
  1067.       if (zCross[i] == ZERO)
  1068.         break;
  1069.       else if (zCross[i] == UP || zCross[i] == DOWN)
  1070.         printf
  1071.           ("BOUNDDETERM error: UP/DOWN crossing instead of ZERO\n");
  1072.     }
  1073.     *zLow = i;
  1074.     linePt = data[*zLow];
  1075.     transPt = data[zStart];
  1076.     endPt = data[zEnd + 1];     /* 8Jan88 -- extra pt = better accuracy */
  1077.     trans1 = transPt;
  1078.     trans2 = endPt;
  1079.   }
  1080.  
  1081. /* determine curve parameters, and store them */
  1082.   midArcPt.x = data[(zEnd + zStart) / 2].x;
  1083.   midArcPt.y = data[(zEnd + zStart) / 2].y;
  1084.   /* if isolated arc with no boundary lines */
  1085.   if (zStart == 1 && zEnd == nTheta - 2) {
  1086.     arcLength = theta[zEnd + 1].length - theta[zStart - 1].length;
  1087.     if (feat2bound (transPt, endPt, arcLength, midArcPt, &radiusC, ¢er,
  1088.                     &curveSweep) == NOBOUNDCURVE)
  1089.       return (0);
  1090.   }
  1091.   /* if boundary curve with 1 boundary line */
  1092.   else {
  1093.     if (featbound (linePt, transPt, endPt, &radiusC, ¢er,
  1094.                    &intrsct, &curveSweep) == NOBOUNDCURVE)
  1095.       return (0);
  1096.   }
  1097.   if (((trans1.x == midArcPt.x) && (trans1.y == midArcPt.y))
  1098.       || ((trans2.x == midArcPt.x) && (trans2.y == midArcPt.y)));
  1099.   else
  1100.     curveDirn = curvedirn (trans1, trans2, midArcPt, center);
  1101.   linkcurve (intrsct, trans1, trans2, center, radiusC, curveDirn);
  1102.  
  1103.   return (0);
  1104. }
  1105.  
  1106.  
  1107. /* FEATBOUND:   function calculates curve parameters for boundary feature
  1108.  *                    usage: flag = featbound (linePt, transPt, endPt,
  1109.  *                                  &radius, ¢erPt, &jctPt, &curveSweep)
  1110.  *              Function returns -1 if there is no curve, rather a straight
  1111.  *              line.
  1112.  *
  1113.  */
  1114.  
  1115. #include <images.h>             /* for image format */
  1116.  
  1117. long
  1118. featbound (linePt, transPt, endPt, radius, centerPt, jctPt, curveSweep)
  1119.      struct point linePt,       /* pt on line, not transition pt */
  1120.        transPt,                 /* transition pt between curve/line */
  1121.        endPt;                   /* final pt on curve */
  1122.  
  1123.      double *radius;            /* output curve radius */
  1124.      struct dpoint *centerPt;   /* output curve center */
  1125.      struct dpoint *jctPt;      /* junction of curve tangents */
  1126.      short *curveSweep;         /* indicates sweep of curve arc */
  1127. {
  1128.   double denom,                 /* denominator */
  1129.     A, B,                       /* constants */
  1130.     sqrt ();
  1131.  
  1132. /* calculate curve parameters */
  1133.   denom = (double) (linePt.y - transPt.y);
  1134.   if (denom != 0.0) {
  1135.     A = (double) (transPt.x - linePt.x) / denom;
  1136.     B = (double) (transPt.x * linePt.x - transPt.x * transPt.x
  1137.                   + transPt.y * linePt.y - transPt.y * transPt.y) / denom;
  1138.     denom = 2.0 * (endPt.x + A * endPt.y - A * transPt.y - transPt.x);
  1139.     if (denom == 0.0)
  1140.       return (NOBOUNDCURVE);
  1141.     centerPt->x = (endPt.x * endPt.x + endPt.y * endPt.y
  1142.                    - 2.0 * B * endPt.y - transPt.x * transPt.x
  1143.                    - transPt.y * transPt.y + 2.0 * B * transPt.y)
  1144.       / denom;
  1145.     centerPt->y = A * centerPt->x + B;
  1146.   }
  1147.   else {
  1148.     centerPt->x = transPt.x;
  1149.     denom = 2.0 * (endPt.y - transPt.y);
  1150.     if (denom == 0.0)
  1151.       return (NOBOUNDCURVE);
  1152.     centerPt->y = (endPt.x * endPt.x - 2.0 * endPt.x * centerPt->x
  1153.                    + endPt.y * endPt.y - transPt.y * transPt.y
  1154.                    + 2.0 * transPt.x * centerPt->x
  1155.                    - transPt.x * transPt.x) / denom;
  1156.   }
  1157.  
  1158.   *radius = sqrt ((transPt.x - centerPt->x) * (transPt.x - centerPt->x)
  1159.                   + (transPt.y - centerPt->y) * (transPt.y - centerPt->y));
  1160.  
  1161.  
  1162. /* calculate intersection of tangents to curve */
  1163.   curvecorner (linePt, transPt, endPt, *radius, *centerPt,
  1164.                jctPt, curveSweep);
  1165.  
  1166.   return (0);
  1167. }
  1168.  
  1169.  
  1170. /* CURVECORNER: intersection point of lines tangent to the transition
  1171.  *            point and at the end point of an end curve feature
  1172.  *                      curvecorner (linePt, transPt, endPt, radius,
  1173.  *                                      centerPt, &jctPt, &curveSweep)
  1174.  */
  1175.  
  1176. #define PID2 1.570796327
  1177.  
  1178. long
  1179. curvecorner (linePt, transPt, endPt, radius, centerPt, jctPt, curveSweep)
  1180.      struct point linePt,       /* pt on line, not transition pt */
  1181.        transPt,                 /* transition pt between curve/line */
  1182.        endPt;                   /* final pt on curve */
  1183.      double radius;             /* output curve radius */
  1184.      struct dpoint centerPt;    /* center of curve */
  1185.      struct dpoint *jctPt;      /* jct pt of tangents to curve */
  1186.      short *curveSweep;         /* indicates sweep of curve arc */
  1187. {
  1188.   double denomE, denomT,        /* end, trans slope denominators */
  1189.     mTrans, mEnd,               /* end, transition line slopes */
  1190.     bTrans, bEnd,               /* end, transition line slopes */
  1191.     distTJSqr,                  /* sqr dist. jct to transition pt */
  1192.     distLJSqr;                  /* sqr dist. jct to line pt */
  1193.  
  1194.   denomE = endPt.y - centerPt.y;
  1195.   denomT = transPt.x - linePt.x;
  1196.  
  1197. /* parallel tangents */
  1198.   if (denomE == 0.0 && denomT == 0.0) {
  1199.     *curveSweep = CURVEEQ180;
  1200.     jctPt->x = (double) linePt.x;
  1201.     jctPt->y = (double) linePt.y;
  1202.     return (0);
  1203.   }
  1204.  
  1205. /* not parallel tangents */
  1206.   else if (denomE == 0.0) {
  1207.     jctPt->x = endPt.x;
  1208.     mTrans = (transPt.y - linePt.y) / denomT;
  1209.     bTrans = -mTrans * linePt.x + linePt.y;
  1210.     jctPt->y = mTrans * jctPt->x + bTrans;
  1211.   }
  1212.   else if (denomT == 0.0) {
  1213.     jctPt->x = transPt.x;
  1214.     mEnd = -(endPt.x - centerPt.x) / denomE;
  1215.     bEnd = endPt.y - mEnd * endPt.x;
  1216.     jctPt->y = mEnd * jctPt->x + bEnd;
  1217.   }
  1218.   else {
  1219.     mEnd = -(endPt.x - centerPt.x) / denomE;
  1220.     bEnd = -mEnd * endPt.x + endPt.y;
  1221.     mTrans = (transPt.y - linePt.y) / denomT;
  1222.     bTrans = -mTrans * linePt.x + linePt.y;
  1223.     if (mEnd == mTrans) {
  1224.       *curveSweep = CURVEEQ180;
  1225.       jctPt->x = (double) linePt.x;
  1226.       jctPt->y = (double) linePt.y;
  1227.       return (0);
  1228.     }
  1229.     jctPt->x = (bEnd - bTrans) / (mTrans - mEnd);
  1230.     jctPt->y = mEnd * jctPt->x + bEnd;
  1231.   }
  1232.  
  1233. /* check if tangents diverge or not */
  1234.   distTJSqr = (jctPt->x - transPt.x) * (jctPt->x - transPt.x)
  1235.     + (jctPt->y - transPt.y) * (jctPt->y - transPt.y);
  1236.   distLJSqr = (jctPt->x - linePt.x) * (jctPt->x - linePt.x)
  1237.     + (jctPt->y - linePt.y) * (jctPt->y - linePt.y);
  1238.   *curveSweep = (distLJSqr > distTJSqr) ? CURVELT180 : CURVEGT180;
  1239.  
  1240.   return (0);
  1241. }
  1242.  
  1243. /* FEAT2BOUND:  function calculates curve parameters for an isolated arc
  1244.  *            (i.e. a curve which is boundary on 2 sides)
  1245.  *                      usage: flag = feat2bound (strtPt, endPt, arcLength,
  1246.  *                                      midArcPt, &radius, ¢erPt,
  1247.  *                                       &curveSweep)
  1248.  *              Function returns a -1 if there is no curve.
  1249.  *
  1250.  *              When there are no bounding lines on the arc, we calculate
  1251.  *              the curve from the arc length. (Since arc length is not
  1252.  *              terribly reliable, we only use this method when lacking
  1253.  *              bounding lines for other curve features.)
  1254.  *
  1255.  *              Note: If the arc length is too large with respect to the
  1256.  *                    distance between the points, the radius is set equal
  1257.  *                    to half the circumference of the circle with the
  1258.  *                    transition-to-curve point diameter.
  1259.  */
  1260.  
  1261. #define MAXRAD 100.0            /* this should be put in ww.h, made a fct
  1262.                                  * * of something and generally cleaned up ?? */
  1263.  
  1264. long
  1265. feat2bound (strtPt, endPt, arcLength, midArcPt, radius, centerPt, curveSweep)
  1266.      struct point strtPt,       /* start pt of curve */
  1267.        endPt;                   /* final pt on curve */
  1268.      double arcLength;          /* curve arc length */
  1269.      struct dpoint midArcPt;    /* middle point on arc */
  1270.      double *radius;            /* output curve radius */
  1271.      struct dpoint *centerPt;   /* output curve center */
  1272.      short *curveSweep;         /* indicates sweep of curve arc */
  1273. {
  1274.   double dist,                  /* Euclid. dist start to end curve */
  1275.     beta,                       /* angle center to line and curve */
  1276.     perpend,                    /* dist center to <dist> line */
  1277.     xMidPt, yMidPt,             /* midpt: line end to curve end */
  1278.     A, B, C, D,                 /* constants */
  1279.     xCenter1, yCenter1,         /* potential center points */
  1280.     xCenter2, yCenter2, distSqr1,  /* sqr dist from line to potential */
  1281.     distSqr2,                   /*  center pts */
  1282.     temp;
  1283.  
  1284. /* initial calculations */
  1285.   xMidPt = (strtPt.x + endPt.x) * 0.5;
  1286.   yMidPt = (strtPt.y + endPt.y) * 0.5;
  1287.   dist = sqrt ((double) ((endPt.x - strtPt.x) * (endPt.x - strtPt.x)
  1288.                          + (endPt.y - strtPt.y) * (endPt.y - strtPt.y)));
  1289.  
  1290. /* calculate radius of curvature by iterative half interval search */
  1291.   *radius = radsearch (arcLength, dist, MAXRAD);
  1292.   if (*radius == NOBOUNDCURVE)
  1293.     return (0);
  1294.  
  1295. /* calculate angle of curve sweep */
  1296.   beta = arcLength / *radius;
  1297.   if (beta < PI)
  1298.     *curveSweep = CURVELT180;
  1299.   else if (beta > PI)
  1300.     *curveSweep = CURVEGT180;
  1301.   else
  1302.     *curveSweep = CURVEEQ180;
  1303.  
  1304. /* calculate center of curvature */
  1305.   if (strtPt.x == endPt.x) {
  1306.     yCenter1 = yCenter2 = yMidPt;
  1307.     A = (*radius * *radius
  1308.          - (strtPt.y - yCenter1) * (strtPt.y - yCenter1));
  1309.     if (A < 0.0)
  1310.       A = 0.0;
  1311.     xCenter1 = strtPt.x + A;
  1312.     xCenter2 = strtPt.x - A;
  1313.   }
  1314.   else if (strtPt.y == endPt.y) {
  1315.     xCenter1 = xCenter2 = xMidPt;
  1316.     A = (*radius * *radius
  1317.          - (strtPt.x - xCenter1) * (strtPt.x - xCenter1));
  1318.     if (A < 0.0)
  1319.       A = 0.0;
  1320.     yCenter1 = strtPt.y + A;
  1321.     yCenter2 = strtPt.y - A;
  1322.   }
  1323.   else {
  1324.     perpend = *radius * cos (beta * 0.5);
  1325.     A = (endPt.y - strtPt.y) * (endPt.y - strtPt.y);
  1326.     B = A + strtPt.x * strtPt.x + endPt.x * endPt.x
  1327.       - 2.0 * strtPt.x * endPt.x;
  1328.     C = -2.0 * A * xMidPt - 2.0 * strtPt.x * strtPt.x * xMidPt
  1329.       + 4.0 * strtPt.x * endPt.x * xMidPt
  1330.       - 2.0 * endPt.x * endPt.x * xMidPt;
  1331.     D = A * xMidPt * xMidPt + strtPt.x * strtPt.x * xMidPt * xMidPt
  1332.       + endPt.x * endPt.x * xMidPt * xMidPt
  1333.       - 2.0 * strtPt.x * endPt.x * xMidPt * xMidPt
  1334.       - A * perpend * perpend;
  1335.  
  1336.     temp = C * C - 4.0 * B * D;
  1337.     if (temp < 0.0)
  1338.       temp = 0.0;
  1339.     xCenter1 = (-C + sqrt (temp)) / (2.0 * B);
  1340.     yCenter1 = ((strtPt.x - endPt.x) * (xCenter1 - xMidPt))
  1341.       / (endPt.y - strtPt.y) + yMidPt;
  1342.     temp = C * C - 4.0 * B * D;
  1343.     if (temp < 0.0)
  1344.       temp = 0.0;
  1345.     xCenter2 = (-C - sqrt (temp)) / (2.0 * B);
  1346.     yCenter2 = ((strtPt.x - endPt.x) * (xCenter2 - xMidPt))
  1347.       / (endPt.y - strtPt.y) + yMidPt;
  1348.   }
  1349.  
  1350.   distSqr1 = (xCenter1 - midArcPt.x) * (xCenter1 - midArcPt.x)
  1351.     + (yCenter1 - midArcPt.y) * (yCenter1 - midArcPt.y);
  1352.   distSqr2 = (xCenter2 - midArcPt.x) * (xCenter2 - midArcPt.x)
  1353.     + (yCenter2 - midArcPt.y) * (yCenter2 - midArcPt.y);
  1354.  
  1355.   if (distSqr1 > distSqr2) {
  1356.     centerPt->x = xCenter1;
  1357.     centerPt->y = yCenter1;
  1358.   }
  1359.   else {
  1360.     centerPt->x = xCenter2;
  1361.     centerPt->y = yCenter2;
  1362.   }
  1363.  
  1364.   return (0);
  1365. }
  1366.  
  1367.  
  1368. /* RADSEARCH:   function calculates radius function by iterative
  1369.  *            half interval search
  1370.  *                      usage: radius = radsearch (arcLength, dist, maxRad)
  1371.  *                              double radius, radsearch ();
  1372.  *              Note: I don't know if the comparison of the <error> with
  1373.  *                    MAXERROR is the best test (should really be testing
  1374.  *                    the error in radius, but this has problems also). I
  1375.  *                    think this is okay, because of the search for
  1376.  *                    zero-crossing, preceding, which results in the
  1377.  *                    half-interval search having far fewer iterations.
  1378.  *                    In any event, by making the MAXERROR small, the
  1379.  *                    error in <radius> is decreased, and it seems to
  1380.  *                    converge pretty quickly (3-5 iterations).
  1381.  */
  1382.  
  1383. #define NINTERVALX 10           /* no. intervals to search 0-cross */
  1384. #define MAXERROR 0.001          /* error as fraction of radius */
  1385.  
  1386. double
  1387. radsearch (arcLength, dist, maxRad)
  1388.      double arcLength,          /* curve arc length */
  1389.        dist,                    /* dist beginning to end of curve */
  1390.        maxRad;                  /* maximum allowable radius */
  1391. {
  1392.   double minX, maxX,            /* x range */
  1393.     error,                      /* fct value - should tend to 0 */
  1394.     deltaX,                     /* interval in x */
  1395.     x,                          /* function variable */
  1396.     higherX,                    /* region contains x zero crossing */
  1397.     lowerX, sin ();
  1398.  
  1399.   long nIter;                   /* iteration number */
  1400.  
  1401.   minX = arcLength * 0.5 / maxRad;
  1402.   maxX = PI;
  1403.  
  1404. /* if arclength is too large, return largest possible radius */
  1405.   if (arcLength >= (PI * dist))
  1406.     return (PI * dist);
  1407. /* if curvature is too small, return flag indicating this */
  1408.   error = sin (minX) - dist * 0.5 / maxRad;
  1409.   if (error <= 0.0)
  1410.     return (NOBOUNDCURVE);
  1411.  
  1412. /* find zero crossing region */
  1413.   deltaX = (maxX - minX) / NINTERVALX;
  1414.   for (x = minX + deltaX; x <= maxX; x += deltaX) {
  1415.     error = sin (x) - (dist / arcLength) * x;
  1416.     if (error <= 0) {
  1417.       higherX = x;
  1418.       lowerX = higherX - deltaX;
  1419.       break;
  1420.     }
  1421.   }
  1422.  
  1423. /* perform iterative half interval search on region with zero crossing */
  1424.   x = (lowerX + higherX) * 0.5;
  1425.   deltaX = (higherX - x) * 0.5;
  1426.   error = sin (x) - (dist / arcLength) * x;
  1427.   nIter = 0;
  1428.   while (ABS (error) > MAXERROR) {  /* see note above */
  1429.     nIter++;
  1430.     if (error > 0) {
  1431.       x += deltaX;
  1432.       deltaX *= 0.5;
  1433.     }
  1434.     else if (error < 0) {
  1435.       x -= deltaX;
  1436.       deltaX *= 0.5;
  1437.     }
  1438.     error = sin (x) - (dist / arcLength) * x;
  1439.   }
  1440.  
  1441.   return ((arcLength * 0.5) / x);
  1442. }
  1443.  
  1444.  
  1445.  
  1446. /* CURVEDIRN:   function finds direction of curve and returns a 1 if the
  1447.  *            arc points are already ordered in counterclockwise (ccw)
  1448.  *              order or a -1 if they are ordered clockwise (cw)
  1449.  *                      usage: dirn = curvedirn (arc1, arc2, midArc, center)
  1450.  *
  1451.  */
  1452.  
  1453. /*#define PI 3.14159265358979 */
  1454. #define PIT2 6.283185307
  1455. #define PID2 1.570796327
  1456. #define PI3D2 4.71238898
  1457.  
  1458. short
  1459. curvedirn (arc1, arc2, midArc, center)
  1460.      struct point arc1, arc2;   /* boundary points of arc */
  1461.      struct dpoint midArc;      /* point on arc about midway from start */
  1462.      struct dpoint center;      /* center of curve */
  1463. {
  1464.   double angle1, angle2,        /* angle about center of arc1, arc2 */
  1465.     angleM,                     /* angle about center of mid arc point */
  1466.     theta12,                    /* angle from arc1 to arc2 ccw wrt center */
  1467.     theta1M;                    /* angle from arc1 to mid arc */
  1468.  
  1469. /* calculate the ccw angles from arc1 to midArc and arc2 around center */
  1470.   angle1 = atan2 ((double) arc1.y - center.y, (double) arc1.x - center.x);
  1471.   if (angle1 < 0.0)
  1472.     angle1 += PIT2;
  1473.   angleM = atan2 ((double) midArc.y - center.y, (double) midArc.x - center.x);
  1474.   if (angleM < 0.0)
  1475.     angleM += PIT2;
  1476.   angle2 = atan2 ((double) arc2.y - center.y, (double) arc2.x - center.x);
  1477.   if (angle2 < 0.0)
  1478.     angle2 += PIT2;
  1479.  
  1480.   theta1M = angleM - angle1;
  1481.   if (theta1M < 0.0)
  1482.     theta1M += PIT2;
  1483.   theta12 = angle2 - angle1;
  1484.   if (theta12 < 0.0)
  1485.     theta12 += PIT2;
  1486.  
  1487.   if (theta12 > theta1M)
  1488.     return (1);
  1489.   else
  1490.     return (-1);
  1491. }
  1492.  
  1493.  
  1494.  
  1495. /* LINESINTRSCT:        function calculates intersection point between
  1496.  *                    two lines
  1497.  *                          usage: flag = linesintrsct (line1Strt, line1End, 
  1498.  *                                            line2Strt, line2End, &intrsct)
  1499.  *                      Returns 0 if intersection, -1 if no intersection
  1500.  *                      and infinite slopes, -1 if no intersection and equal
  1501.  *                      (but not infinite) slopes.
  1502.  */
  1503.  
  1504. #define NOTEQUALSLOPES 0        /* so there is an intersection pt */
  1505. #define INFINITESLOPES2 -1      /* vertical slopes = no intersectn */
  1506. #define EQUALSLOPES2 -2         /* equal slopes = no intersection */
  1507.  
  1508. long
  1509. linesintrsct (line1Strt, line1End, line2Strt, line2End, intrsct)
  1510.      struct point line1Strt;
  1511.      struct dpoint line1End;    /* start/end coords of lines */
  1512.      struct point line2Strt, line2End;
  1513.      struct dpoint *intrsct;    /* intersection point */
  1514. {
  1515.   double deltaX1, deltaX2,      /* x increments of lines */
  1516.     m1, m2,                     /* slopes of lines */
  1517.     b1, b2;                     /* y-intersects of lines */
  1518.  
  1519. /* calculate intersection of lines */
  1520.   deltaX1 = (double) (line1Strt.x - line1End.x);
  1521.   deltaX2 = (double) (line2Strt.x - line2End.x);
  1522.   if (deltaX1 == 0.0 && deltaX2 == 0.0)
  1523.     return (INFINITESLOPES2);
  1524.   else if (deltaX1 == 0.0) {
  1525.     intrsct->x = (double) line1Strt.x;
  1526.     m2 = (double) (line2Strt.y - line2End.y) / deltaX2;
  1527.     b2 = -m2 * line2End.x + line2End.y;
  1528.     intrsct->y = m2 * intrsct->x + b2;
  1529.   }
  1530.   else if (deltaX2 == 0.0) {
  1531.     intrsct->x = (double) line2Strt.x;
  1532.     m1 = (line1Strt.y - line1End.y) / deltaX1;
  1533.     b1 = -m1 * line1End.x + line1End.y;
  1534.     intrsct->y = m1 * intrsct->x + b1;
  1535.   }
  1536.   else {
  1537.     m1 = (double) (line1Strt.y - line1End.y) / deltaX1;
  1538.     b1 = -m1 * line1Strt.x + line1Strt.y;
  1539.     m2 = (line2Strt.y - line2End.y) / deltaX2;
  1540.     b2 = -m2 * line2End.x + line2End.y;
  1541.     if (m1 == m2)
  1542.       return (EQUALSLOPES2);
  1543.     intrsct->x = (b1 - b2) / (m2 - m1);
  1544.     intrsct->y = m1 * intrsct->x + b1;
  1545.   }
  1546.   return (NOTEQUALSLOPES);
  1547. }
  1548.