-- CGStrokeImpl.mesa
-- Last changed by Doug Wyatt, September 18, 1982 6:58 pm

DIRECTORY
  CGArea USING [Empty, Ref],
  CGClipper USING [GenerateLine, Load, Ref],
  CGContext USING [Ref],
  CGCubic USING [Bezier, Flat, Split],
  CGDevice USING [Ref],
  CGMatrix USING [Map, Ref],
  CGPath USING [Empty, MapAndFilter],
  CGPrivate USING [Context],
  CGReducer USING [Ref, Vertex, Generate],
  CGSource USING [Ref],
  CGVector USING [Add, In, Max, Min, Sub, Dot],
  GraphicsBasic USING [Path, StrokeEnds, Vec],
  Real USING [SqRt];

CGStrokeImpl: CEDAR PROGRAM
IMPORTS CGArea, CGClipper, CGCubic, CGMatrix, CGPath, CGReducer, V: CGVector, Real
EXPORTS CGPrivate = {

Context: TYPE = CGPrivate.Context;
ContextData: TYPE = CGContext.Ref;
Path: TYPE = GraphicsBasic.Path;
StrokeEnds: TYPE = GraphicsBasic.StrokeEnds;
Vec: TYPE = GraphicsBasic.Vec;
Bezier: TYPE = CGCubic.Bezier;

DrawStroke: PUBLIC PROC[self: Context, path: Path, width: REAL,
  closed: BOOLEAN, ends: StrokeEnds] = {
  ctx: ContextData _ NARROW[self.data];
  IF CGPath.Empty[path] THEN RETURN;
  IF width=0 THEN DrawThinStroke[ctx, path, closed]
  ELSE DrawThickStroke[ctx, path, width, closed, ends];
  };

maxdepth: NAT = 10;

Subdivide: PROC[b: Bezier, vertex: PROC[Vec], tolerance: REAL _ 1.5, depth: NAT _ 0] = {
  IF depth>=maxdepth OR CGCubic.Flat[b, tolerance] THEN vertex[b.b3]
  ELSE { b1, b2: Bezier;
    [b1, b2] _ CGCubic.Split[b];
    Subdivide[b1, vertex, tolerance, depth+1]; -- first half
    Subdivide[b2, vertex, tolerance, depth+1]; -- second half
    };
  };

DrawThinStroke: PROC[ctx: ContextData, path: Path, closed: BOOLEAN] = {
  matrix: CGMatrix.Ref _ ctx.matrix;
  clipper: CGClipper.Ref _ ctx.clipper;
  area: CGArea.Ref _ ctx.area;
  first, last: Vec;
  Map: PROC[v: Vec] RETURNS[Vec] = { RETURN[CGMatrix.Map[matrix, v]] };
  Move: PROC[v: Vec] = { first _ last _ v };
  Line: PROC[v: Vec] = { CGClipper.GenerateLine[clipper, last, v, area]; last _ v };
  Curve: PROC[v1, v2, v3: Vec] = { Subdivide[[last,v1,v2,v3], Line] };
  Close: PROC = { IF closed THEN Line[first] };
  CGPath.MapAndFilter[path, Map, Move, Line, Curve, Close];
  IF NOT CGArea.Empty[area] THEN {
    device: CGDevice.Ref _ ctx.device;
    src: CGSource.Ref _ ctx.src;
    fat: BOOLEAN _ src.fat;
    src.fat _ TRUE;
    device.Show[device,area,src,ctx.matrix];
    src.fat _ fat;
    };
  };

DrawThickStroke: PROC[ctx: ContextData, path: Path, width: REAL,
  closed: BOOLEAN, ends: StrokeEnds] = {
  matrix: CGMatrix.Ref _ ctx.matrix;
  clipper: CGClipper.Ref _ ctx.clipper;
  reducer: CGReducer.Ref _ ctx.reducer;
  firstVertex: BOOLEAN _ TRUE;
  Map: PROC[v: Vec] RETURNS[Vec] = { RETURN[CGMatrix.Map[matrix, v]] };
  Vertex: PROC[v: Vec] = {
    IF firstVertex THEN { CGClipper.Load[clipper, reducer]; firstVertex _ FALSE };
    CGReducer.Vertex[reducer, v];
    };
  Close: PROC = {
    area: CGArea.Ref _ ctx.area;
    IF firstVertex THEN RETURN;
    CGReducer.Generate[reducer, area];
    IF NOT CGArea.Empty[area] THEN {
      src: CGSource.Ref _ ctx.src;
      device: CGDevice.Ref _ ctx.device;
      device.Show[device, area, src, matrix];
      };
    firstVertex _ TRUE;
    };
  GenerateStroke[path, width, closed, ends, Map, Vertex, Close];
  };

Mag: PROC[v: Vec] RETURNS [REAL] = INLINE {
  RETURN [Real.SqRt[v.x * v.x + v.y * v.y]]
  };

Normalize: PROC[v: Vec] RETURNS [Vec] = INLINE {
  m: REAL _ Mag[v];
  IF m = 0 THEN RETURN[v] ELSE RETURN[[v.x/m, v.y/m]]
  };

Midpoint: PROC[a, b: Vec] RETURNS[Vec] = INLINE {
  RETURN[[(a.x+b.x)/2, (a.y+b.y)/2]] };

Miter: PROC[v1, v2: Vec] RETURNS[REAL] = {
  s, a: REAL; v3: Vec;
  v3 _ Normalize[V.Add[v1, v2]];
  s _ V.Dot[[v1.y, -v1.x], v3]; -- compute sin of angle
  a _ s / Real.SqRt[1 - s*s]; -- compute horz vec for angle between v1 & v3
  RETURN[a];
  };

-- Three point flatness test:
Flat: PROC[v1, v2, v3: Vec, epsilon: REAL] RETURNS[b: BOOL] = {
  dx, dy: REAL;
  d1, d, bl, bh: Vec;
  oh: Vec = [0.5, 0.5];
  bh _ V.Add[V.Max[v1, v3], oh];
  bl _ V.Sub[V.Min[v1, v3], oh];
  IF NOT V.In[v2, bl, bh] THEN RETURN[FALSE];
  d1 _ V.Sub[v2, v1];
  d _ V.Sub[v3, v1];
  dx _ ABS[d.x]; dy _ ABS[d.y];
  IF dx + dy < 1 THEN RETURN[TRUE];
  IF dy < dx THEN { dydx: REAL _ d.y / d.x; RETURN[ABS[d1.y - d1.x * dydx] < epsilon] }
  ELSE { dxdy: REAL _ d.x / d.y; RETURN[ABS[d1.x - d1.y * dxdy] < epsilon] };
  };

GenerateStroke: PUBLIC PROC[path: Path, width: REAL, closed: BOOLEAN, ends: StrokeEnds,
  map: PROC[Vec] RETURNS[Vec], vertex: PROC[Vec], close: PROC] = {
  h: REAL _ ABS[width]/2; -- half width
  hsq: REAL _ h*h;
  smallSin: REAL _ 0.001*hsq;
  line: BOOLEAN; -- is current segment a line?
  b0, b1, b2, b3: Vec; -- control points for current segment
  u0, u3: Vec; -- half width normal vectors at b0, b3
  p0, q0, p3, q3: Vec; -- corners of stroke boundary for current segment
  MiterFlags: TYPE = RECORD[p, q: BOOLEAN];
  m0, m3: MiterFlags;
  fb0, fb1, fb2, fb3, fu0, fu3, fp3, fq3: Vec; -- saved values for deferred segment
  fline: BOOLEAN; fm3: MiterFlags; -- ditto
  begin, defer: BOOLEAN;
  Cap: PROC[b, d, p, q: Vec] = {
    SELECT ends FROM
      butt => { };
      square => {
        v0: Vec _ [b.x-d.y, b.y+d.x];
        v3: Vec _ [b.x+d.y, b.y-d.x];
        v1: Vec _ [v0.x+d.x, v0.y+d.y];
        v2: Vec _ [v3.x+d.x, v3.y+d.y];
        vertex[p]; vertex[map[v1]]; vertex[map[v2]]; vertex[q]; close[];
        };
      round => {
        tolerance: REAL = 1.0;
        bz: Bezier;
        v0, v1, v2, v3: Vec;
        f: REAL = 0.5522848; -- 4*(sqrt2-1)/3
        a: Vec _ [f*d.x, f*d.y];
        v0 _ [b.x-d.y, b.y+d.x];
        v3 _ [b.x+d.x, b.y+d.y];
        v1 _ [v0.x+a.x, v0.y+a.y];
        v2 _ [v3.x-a.y, v3.y+a.x];
        bz.b0 _ p; bz.b1 _ map[v1]; bz.b2 _ map[v2]; bz.b3 _ map[v3];
        vertex[p]; Subdivide[bz, vertex, tolerance];
        v0 _ v3;
        v3 _ [b.x+d.y, b.y-d.x];
        v1 _ [v0.x+a.y, v0.y-a.x];
        v2 _ [v3.x+a.x, v3.y+a.y];
        bz.b0 _ bz.b3; bz.b1 _ map[v1]; bz.b2 _ map[v2]; bz.b3 _ q;
        Subdivide[bz, vertex, tolerance]; close[];
        };
      ENDCASE => ERROR;
    };
  HNormal: PROC[a, b: Vec] RETURNS[Vec] = {
    -- returns a vector of length halfWidth normal to the direction from a to b
    d: Vec _ V.Sub[b, a];
    u: Vec _ Normalize[d];
    RETURN[[-u.y*h, u.x*h]];
    };
  ComputeP: PROC[b, u: Vec] RETURNS[Vec] = INLINE { RETURN[map[[b.x+u.x, b.y+u.y]]] };
  ComputeQ: PROC[b, u: Vec] RETURNS[Vec] = INLINE { RETURN[map[[b.x-u.x, b.y-u.y]]] };
  tolerance: REAL = 1.0;
  DoCurve: PROC[b0, b1, b2, b3: Vec, p0, p3: Vec, depth: NAT _ 0] = {
    b01: Vec _ Midpoint[b0, b1];
    b12: Vec _ Midpoint[b1, b2];
    b23: Vec _ Midpoint[b2, b3];
    b11: Vec _ Midpoint[b01, b12];
    b22: Vec _ Midpoint[b12, b23];
    bm: Vec _ Midpoint[b11, b22]; -- middle point on curve (at t=1/2)
    um: Vec _ HNormal[b11, b22]; -- half width normal at bm
    pm: Vec _ map[V.Add[bm, um]]; -- point on stroke edge
    IF depth>0 -- always divide at least once --
     AND Flat[p0, pm, p3, tolerance] THEN vertex[p3]
    ELSE {
      DoCurve[b0, b01, b11, bm, p0, pm, depth+1]; -- first half
      DoCurve[bm, b22, b23, b3, pm, p3, depth+1]; -- second half
      };
    };
  EmitSegment: PROC = {
    IF NOT m0.p THEN p0 _ ComputeP[b0, u0];
    IF NOT m0.q THEN q0 _ ComputeQ[b0, u0];
    IF NOT m3.p THEN p3 _ ComputeP[b3, u3];
    IF NOT m3.q THEN q3 _ ComputeQ[b3, u3];
    vertex[p0];
    IF line THEN vertex[p3] ELSE DoCurve[b0, b1, b2, b3, p0, p3];
    vertex[q3];
    IF line THEN vertex[q0] ELSE DoCurve[b3, b2, b1, b0, q3, q0];
    close[];
    };
  EmitPreviousSegment: PROC[u4: Vec, miter: BOOLEAN _ TRUE] = {
    IF begin THEN {
      IF closed THEN {
        fb0 _ b3; fu0 _ u4;
        defer _ TRUE }
      ELSE {
        u3 _ u4;
        p3 _ ComputeP[b3, u3]; q3 _ ComputeQ[b3, u3];
        m3 _ [TRUE,TRUE];
        Cap[b3, [-u3.y, u3.x], q3, p3];
        defer _ FALSE };
      begin _ FALSE }
    ELSE {
      IF miter THEN {
        sin: REAL _ u3.x*u4.y-u3.y*u4.x;
        IF ABS[sin]<smallSin THEN {
          cos: REAL _ u3.x*u4.x+u3.y*u4.y;
          IF cos>0 THEN {
            p3 _ ComputeP[b3, u3]; q3 _ ComputeQ[b3, u3];
            m3 _ [TRUE,TRUE] }
          ELSE m3 _ [FALSE,FALSE];
          }
        ELSE {
          f: REAL _ hsq/sin;
          r: Vec _ [f*(u4.y-u3.y), f*(u3.x-u4.x)];
          IF sin<0 THEN { p3 _ map[[b3.x+r.x, b3.y+r.y]]; m3 _ [p: TRUE, q: FALSE] }
          ELSE { q3 _ map[[b3.x-r.x, b3.y-r.y]]; m3 _ [p: FALSE, q: TRUE] };
          };
        }
      ELSE m3 _ [FALSE,FALSE];
      IF defer THEN {
        fb1 _ b1; fb2 _ b2; fb3 _ b3; fline _ line;
        fu3 _ u3; fp3 _ p3; fq3 _ q3; fm3 _ m3;
        defer _ FALSE }
      ELSE EmitSegment[];
      };
    b0 _ b3; u0 _ u4; p0 _ p3; q0 _ q3; m0 _ m3;
    };

  Move: PROC[v: Vec] = { b3 _ v; begin _ TRUE };
  Line: PROC[v: Vec] = {
    u0: Vec _ HNormal[b3, v];
    EmitPreviousSegment[u0];
    b3 _ v; u3 _ u0; line _ TRUE;
    };
  Curve: PROC[v1, v2, v3: Vec] = {
    u0: Vec _ (IF v1#b3 THEN HNormal[b3, v1] ELSE HNormal[b3, v2]);
    EmitPreviousSegment[u0];
    b1 _ v1; b2 _ v2; b3 _ v3; line _ FALSE;
    u3 _ (IF b2#b3 THEN HNormal[b2, b3] ELSE HNormal[b1, b2]);
    };
  Close: PROC = {
    IF begin THEN {
      -- trajectory is a single point
      IF closed THEN { --ignore-- }
      ELSE {
        -- just pick an arbitrary direction for end caps
        u3 _ [0, h];
        p3 _ ComputeP[b3, u3];
        q3 _ ComputeQ[b3, u3];
        Cap[b3, [-h, 0], q3, p3];
        Cap[b3, [h, 0], p3, q3];
        };
      }
    ELSE IF closed THEN {
      IF b3#fb0 THEN Line[fb0]; -- if necessary, extend a line to the first point
      EmitPreviousSegment[fu0]; -- emit the segment ending at the first point
      b1 _ fb1; b2 _ fb2; b3 _ fb3; line _ fline;
      u3 _ fu3; p3 _ fp3; q3 _ fq3; m3 _ fm3;
      EmitSegment[]; -- emit the deferred first segment
      }
    ELSE { EmitPreviousSegment[u3, FALSE]; Cap[b3, [u3.y, -u3.x], p3, q3] };
    };
  CGPath.MapAndFilter[path, NIL, Move, Line, Curve, Close];
  };

}.