(* FOR PZCandidate.m3: *)
  
PROCEDURE CandidateOverlap(
    READONLY a, b: PZCandidate.T;
    minSteps: NAT;
    maxAdjust: NAT;
    VAR (*WORK*) ra, rb: REF PZMatch.T;   (* Work areas *)
  ): BOOL =
  BEGIN
    FOR j := 0 TO 1 DO
      WITH
        aj = a.seg[j], bj = b.seg[j],
        ov = PZSegment.Overlap(aj, bj)
      DO
        Wr.PutText(stderr, " " & FI(ov, 4));
        IF ov < minSteps THEN 
          Wr.PutText(stderr, "!");
          RETURN FALSE
        END;
      END
    END;
    WITH
      ov = PackedMatchOverlap(a.pm, b.pm, a.seg, b.seg, maxAdjust, ra,rb)
    DO
      Wr.PutText(stderr, " " & FI(ov, 4));
      IF ov < minSteps THEN 
        Wr.PutText(stderr, "!");
        RETURN FALSE
      END;
    END;
    RETURN TRUE
  END CandidateOverlap;

PROCEDURE ReallocMatch(VAR r: REF PZMatch.T; np: NAT) =
  (*
    Reallocates "r^" if needed to store at least "np" elements.
    May discard the old contents. *)
  BEGIN
    IF r = NIL THEN
      r := NEW(REF PZMatch.T, np)
    ELSIF NUMBER(r^) < np THEN
      r := NEW(REF PZMatch.T, MAX(2*NUMBER(r^), np))
    END
  END ReallocMatch;

PROCEDURE PackedMatchOverlap(
    a, b: REF PZMatch.PackedT;
    READONLY aSeg, bSeg: PZSegment.Pair;
    maxAdjust: NAT;
    VAR (*WORK*) ra, rb: REF PZMatch.T;   (* Work area *)
  ): NAT =
  VAR na, nb: NAT;
  BEGIN
    IF a = NIL THEN na := MAX(aSeg[0].ns, aSeg[1].ns) ELSE na := a.np END;
    ReallocMatch(ra, na);
    IF b = NIL THEN nb := MAX(bSeg[0].ns, bSeg[1].ns) ELSE nb := b.np END;
    ReallocMatch(rb, nb);
    WITH
      am = SUBARRAY(ra^, 0, na),
      bm = SUBARRAY(rb^, 0, nb)
    DO
      IF a = NIL THEN
        PZMatch.DoMostPerfect(aSeg[0].ns, aSeg[1].ns, am)
      ELSE
        PZMatch.DoUnpack(a^, am)
      END;
      IF b = NIL THEN
        PZMatch.DoMostPerfect(bSeg[0].ns, bSeg[1].ns, bm)
      ELSE
        PZMatch.DoUnpack(b^, bm)
      END;
      RETURN MatchingOverlap(am, bm, aSeg, bSeg, maxAdjust)
    END
  END PackedMatchOverlap;


(* FOR PZMatch.m3 *)

TYPE Pair = PZMatch.Pair;

PROCEDURE MatchingOverlap(
    READONLY a, b: PZMatch.T;
    READONLY aSeg, bSeg: PZSegment.Pair;
    maxAdjust: NAT;
  ): NAT =
  
  VAR N: NAT := 0;
      tot: Pair;  (* Total samples in each contour. *)
  
  PROCEDURE Tally(
      ai: NAT; READONLY aiP: Pair;
      af: NAT; READONLY afP: Pair;
      bi: NAT; READONLY biP: Pair;
      bf: NAT; READONLY bfP: Pair;
    ) =
    (*
      Adds TO "N" the number of sample pairs "a[ak]" that 
      are within "maxAdjust" some sample pair "b[bk]", where
      "ak IN [ai..af]" and "bk IN [bi..bf]". 
      
      Assumes that "aiP", "afP", "biP", "bfP" are the sample index
      pairs corresponding to "a[ai]", "a[af]", etc., already mapped to
      absolute sample indices.
    *)
    BEGIN
      Wr.PutText(stderr, "[");
      <* ASSERT ai < af *>
      <* ASSERT bi <= bf *>
      IF NOT BoxesOverlap(aiP, afP, biP, bfP) THEN 
        (* Do nothing. *)
      ELSIF af - ai >= bf - bi AND af - ai >= 2 THEN 
        (* The "a" matching is longer and has at least two steps; split it: *)
        WITH
          am = (ai + af) DIV 2,
          amP = PZMatch.AbsPairClip(a[am], aSeg)
        DO
          <* ASSERT ai < am AND am < af *>
          Tally(ai, aiP, am, amP, bi, biP, bf, bfP);
          Tally(am, amP, af, afP, bi, biP, bf, bfP);
        END
      ELSIF bf - bi >= af - ai AND bf - bi >= 2 THEN 
        (* The "b" matching is longer and has at least two steps; split it: *)
        WITH
          bm = (bi + bf) DIV 2,
          bmP = PZMatch.AbsPairClip(b[bm], bSeg)
        DO
          <* ASSERT bi < bm AND bm < bf *>
          Tally(ai, aiP, af, afP, bi, biP, bm, bmP);
          Tally(ai, aiP, af, afP, bm, bmP, bf, bfP);
        END
      ELSE
        (*
          Both "a" and "b" are reduced to a single step.  Let "B" be
          the step "b" bloated by the square ball with radius
          "maxAdjust".  Since "B" is convex, it suffices to check
          whether both endpoints of "a" are in "B".  Moreover the only
          integer points in "B" are also in the bloating of the two
          endpoints.
        *)
        IF  MIN(TorDist(aiP,biP), TorDist(aiP,bfP)) <= maxAdjust
        AND MIN(TorDist(afP,biP), TorDist(afP,bfP)) <= maxAdjust THEN
          INC(N)
        END
      END;
      Wr.PutText(stderr, "]");
    END Tally;
    
  PROCEDURE BoxesOverlap(READONLY aiP, afP, biP, bfP: Pair): BOOL =
    (*
      Returns true if the bounding box of "a" is at most
      "maxAdjust" away from the bounding box of "b". *)
    BEGIN
      FOR j := 0 TO 1 DO
        Wr.PutText(stderr, " " & FI(aiP[j],1) & "-" & FI(afP[j],1) & ":" & FI(biP[j],1) & "-" & FI(bfP[j],1));
        IF NOT aSeg[j].rev THEN
          WITH
            m = tot[j],
            na = (afP[j] - aiP[j]) MOD m + 1,
            dab = (biP[j] - aiP[j]) MOD m,
            nb = (bfP[j] - biP[j]) MOD m + 1,
            dba = (aiP[j] - biP[j]) MOD m
          DO
            IF dab > na + maxAdjust AND dba > nb + maxAdjust THEN RETURN FALSE END
          END
        ELSE
          WITH
            m = tot[j],
            na = (aiP[j] - afP[j]) MOD m + 1,
            dab = (bfP[j] - afP[j]) MOD m,
            nb = (biP[j] - bfP[j]) MOD m + 1,
            dba = (afP[j] - bfP[j]) MOD m
          DO
            IF dab > na + maxAdjust AND dba > nb + maxAdjust THEN RETURN FALSE END
          END
        END
      END;
      RETURN TRUE
    END BoxesOverlap;
    
  PROCEDURE TorDist(READONLY ap, bp: Pair): NAT =
    (*
      Minimum distance between the pairs "ap" and "bp" on the 
      toroidal grid defined by the two contours. *)
    BEGIN
      WITH
        d0 = MIN((ap[0] - bp[0]) MOD tot[0], (bp[0] - ap[0]) MOD tot[0]), 
        d1 = MIN((ap[1] - bp[1]) MOD tot[1], (bp[1] - ap[1]) MOD tot[1])
      DO
        RETURN MAX(d0, d1)
      END
    END TorDist;
    
  BEGIN
    IF aSeg[0].cvx # bSeg[0].cvx OR aSeg[1].cvx # bSeg[1].cvx THEN RETURN 0 END;
    IF aSeg[0].rev # bSeg[0].rev OR aSeg[1].rev # bSeg[1].rev THEN RETURN 0 END;
    <* ASSERT aSeg[0].tot = bSeg[0].tot AND aSeg[1].tot = bSeg[1].tot *>
    tot[0] := aSeg[0].tot; tot[1] := aSeg[1].tot;
    N := 0;
    WITH
      ai = 0,       aiP = PZMatch.AbsPairClip(a[ai], aSeg),
      af = LAST(a), afP = PZMatch.AbsPairClip(a[af], aSeg),
      bi = 0,       biP = PZMatch.AbsPairClip(b[bi], bSeg),
      bf = LAST(b), bfP = PZMatch.AbsPairClip(b[bf], bSeg)
    DO
      Tally(ai, aiP, af, afP, bi, biP, bf, bfP);
      RETURN N
    END
  END MatchingOverlap;


PROCEDURE OrVecs(READONLY a, b: BOOLS): REF BOOLS =
  BEGIN
    WITH
      na = NUMBER(a),
      nb = NUMBER(b),
      n = MAX(na, nb),
      rc = NEW(REF BOOLS, n), c = rc^
    DO
      FOR i := 0 TO n-1 DO 
        c[i] := FALSE; 
        IF i < na AND a[i] THEN c[i] := TRUE END;
        IF i < nb AND b[i] THEN c[i] := TRUE END;
      END;
      RETURN rc
    END
  END OrVecs;
      


PROCEDURE ReadAllCVCFiles(
    chainPrefix : TEXT;
    nChains, band: CARDINAL; 
    readChains: BOOLEAN;
  ): REF ARRAY OF PZSymbolChain.ReadData =
  <* FATAL Rd.Failure *>
  BEGIN 
   WITH 
     cvcChain = NEW(REF ARRAY OF PZSymbolChain.ReadData, nChains)
   DO 
     FOR i := 0 TO nChains-1 DO  
        WITH 
          chainTag  = Fmt.Pad(Fmt.Int(i), 4, '0'),
          bandTag = Fmt.Pad(Fmt.Int(band), 3, '0'),
          fileName = 
            chainTag & "/" & chainPrefix & "-" & chainTag & "-f" & bandTag & "-c-u.cvc", 
          rd = FileRd.Open(fileName),
          data = PZSymbolChain.Read(rd, headerOnly := NOT readChains)
        DO
          cvcChain[i] := data;
          Rd.Close(rd);
        END;
      END;
      RETURN cvcChain;
    END;  
  END ReadAllCVCFiles;
 
 
(* Last edited on 1999-11-15 12:24:36 by hcgl *)