MODULE GDGeometry;

PROCEDURE InterpolateCoords(u, v: INT; r, s: INT): INT =
  BEGIN
    WITH t = FLOAT(r, LONG)/FLOAT(s, LONG) DO
      RETURN u + ROUND(t * FLOAT(v - u, LONG))
    END
  END InterpolateCoords;
  
PROCEDURE InterpolatePoints(u, v: Point; r, s: INT): Point =
  BEGIN
    WITH t = FLOAT(r, LONG)/FLOAT(s, LONG) DO
      RETURN Point{
        u[0] + ROUND(t * FLOAT(v[0] - u[0], LONG)),
        u[1] + ROUND(t * FLOAT(v[1] - u[1], LONG))
      }
    END
  END InterpolatePoints;
  
PROCEDURE SGN(x: INT): SIGN =
  BEGIN
    IF x < 0 THEN
      RETURN -1
    ELSIF x > 0 THEN
      RETURN +1
    ELSE
      RETURN 0
    END
  END SGN;

PROCEDURE LexCompare(READONLY p, q: Point): SIGN = 
  BEGIN 
    IF p[1] < q[1] THEN
      RETURN -1 
    ELSIF p[1] > q[1] THEN
      RETURN +1
    ELSIF p[0] < q[0] THEN
      RETURN -1
    ELSIF p[0] > q[0] THEN
      RETURN +1
    ELSE
      RETURN 0
    END
  END LexCompare;

PROCEDURE GetBounds(READONLY pt: Points): Box = 
  VAR xMin, xMax, yMin, yMax: INT;
  BEGIN
    IF NUMBER(pt) = 0 THEN RETURN Box{Interval{lo := +1, hi := -1}, ..} END;
    xMin := pt[0][0]; xMax := xMin;
    yMin := pt[0][1]; yMax := yMin;
    FOR i := 1 TO LAST(pt) DO
      xMin := MIN(xMin, pt[i][0]);
      xMax := MAX(xMax, pt[i][0]);
      yMin := MIN(yMin, pt[i][1]);
      yMax := MAX(yMax, pt[i][1]);
    END;
    RETURN Box{ 
      Interval{lo := xMin, hi := xMax},
      Interval{lo := yMin, hi := yMax}
    }
  END GetBounds;
  
PROCEDURE GetRange(READONLY pt: Points; axis: Axis): Interval = 
  VAR cMin, cMax: INT;
  BEGIN
    IF NUMBER(pt) = 0 THEN RETURN Interval{lo := +1, hi := -1} END;
    cMin := pt[0][axis]; cMax := cMin;
    FOR i := 1 TO LAST(pt) DO
      cMin := MIN(cMin, pt[i][axis]);
      cMax := MAX(cMax, pt[i][axis]);
    END;
    RETURN Interval{lo := cMin, hi := cMax}
  END GetRange;
  
PROCEDURE Medians(READONLY pt: Points; axis: Axis): Interval =
  VAR yLo, yHi, yM, yBgLo, ySmHi: INT;
      cM: INT;
  BEGIN
    <* ASSERT NUMBER(pt) > 0 *>
    WITH r = GetRange(pt, axis) DO
      yLo := r.lo-1;
      yHi := r.hi+1
    END;
    LOOP
      <* ASSERT yLo < yHi *>
      (* The median ordinate(s) must be contained in "(yLo..yHi)". *)
      IF yLo + 1 = yHi THEN 
        (* Nothing more to do: *)
        RETURN Interval{lo := yLo, hi := yHi}
      END;
      (* Split the interval: *)
      yM := (yLo + yHi) DIV 2;
      <* ASSERT yLo < yM AND yM < yHi *>
      (* Compute the unbalance between counts "<= yM" and ">= yM". *)
      (* Also computes the extremal ordinates in those two sets. *)
      cM := 0;
      yBgLo := yHi; ySmHi := yLo;
      FOR i := 0 TO LAST(pt) DO
        WITH yi = pt[i][axis] DO
          IF yi < yM THEN 
            cM := cM - 1; ySmHi := MAX(ySmHi, yi); 
          ELSIF yi > yM THEN
            cM := cM + 1; yBgLo := MIN(yBgLo, yi); 
          ELSE
            ySmHi := yi; yBgLo := yi;
          END
        END
      END;
      (* Decide what to do next: *)
      IF cM > 0 THEN
        (* The median(s) must be strictly "> yM" *)
        yLo := yM
      ELSIF cM < 0 THEN
        (* The median(s) must be strictly "< yM" *)
        yHi := yM
      ELSIF ySmHi < yBgLo THEN
        <* ASSERT ySmHi < yM AND yM < yBgLo *>
        (* No points on layer "yM", may pick any oridinate in "(ySmHi..yBgLo)": *)
        RETURN Interval{lo := ySmHi, hi := yBgLo}
      ELSE
        <* ASSERT ySmHi = yM AND yM = yBgLo *>
        (* The unique median is "yM" *)
        RETURN Interval{lo := yM, hi := yM}
      END
    END;
  END Medians;

PROCEDURE Median(READONLY pt: Points; axis: Axis): LONG =
  BEGIN
    WITH r = Medians(pt, axis) DO
      RETURN FLOAT(r.lo + r.hi, LONG) / 2.0d0
    END
  END Median;

PROCEDURE PickOutside(r: Interval; xBar: LONG; dir: SIGN): INT =
  VAR xLo, xHi: INT;
  BEGIN
    xLo := FLOOR(xBar); xHi := CEILING(xBar);
    IF r.lo <= r.hi THEN
      IF xLo <= r.hi THEN xLo := MIN(xLo, r.lo-1) END;
      IF xHi >= r.lo THEN xHi := MAX(xHi, r.hi+1) END;
    END;
    WITH 
      dLo = ABS(FLOAT(xLo,LONG)-xBar),
      dHi = ABS(FLOAT(xHi,LONG)-xBar)
    DO
      IF dLo < dHi THEN
        RETURN xLo
      ELSIF dLo > dHi THEN
        RETURN xHi
      ELSIF dir < 0 THEN
        RETURN xLo
      ELSIF dir > 0 THEN
        RETURN xHi
      ELSE
        <* ASSERT FALSE *>
      END
    END
  END PickOutside;

BEGIN
END GDGeometry.
(* Last edited on 2000-01-10 13:16:37 by stolfi *)