<* UNUSED *>     
PROCEDURE SampleCurvatureAtExtremaOnly(
    READONLY c: PZLRChain.T;
    epsilon, delta: LONGREAL;
  ): REF PZLRChain.T =
  BEGIN
    WITH
      xIndex  = FindIndexExtrema(c, epsilon, delta)^,
      NS = NUMBER(xIndex),
      rs = NEW(REF PZLRChain.T, NS),
      s = rs^
    DO
      FOR j := 0 TO NS-1 DO s[j] := c[xIndex[j]] END;
      RETURN rs
    END
  END SampleCurvatureAtExtremaOnly;

 
PROCEDURE FindIndexExtrema(
    READONLY c: PZLRChain.T;
    epsilon, delta: LONGREAL;
  ): REF IntArray  =
  VAR NX: CARDINAL := 0;
      qPrev, qThis, qNext: INTEGER;
      j, k: CARDINAL;
      iMax, iMin: CARDINAL;
  BEGIN
    WITH
      NC = NUMBER(c),
      x = NEW(REF IntArray, NC)^
    DO
      (* Quantized curvature at point 0: *)
      qThis := ROUND(Squeeze(c[0], epsilon, delta));
      iMax := 0;
      iMin := 0;
      (* Looks for the previous point with a different quantized curvature: *)
      j := NC;
      LOOP
        DEC(j);
        qPrev := ROUND(Squeeze(c[j], epsilon, delta));
        IF qPrev # qThis THEN EXIT END;
        IF c[j] > c[iMin] THEN iMin := j END;
        IF c[j] < c[iMax] THEN iMax := j END;
        IF j = 0 THEN EXIT END;
      END;
      (* Process each point: *)
      FOR i := 0 TO NC-1 DO
        k := (i+1) MOD NC;
        qNext := ROUND(Squeeze(c[k], epsilon, delta));
        (*
          "qThis" is the quantized curvature at point "i"
          "qNext" is the quantized curvature at point "k"
          "j" is the last index before "i (mod "NC") whose quantized
              curvature is different from "qThis"
          "qPrev" is the quantized curvature at point "j".
          "iMax" is the point of maximum curvature in [j..i]
          "iMin" is the point of minimum curvature in [j..i]
        *)
        IF qNext # qThis THEN
          IF (qThis > qNext AND qThis > qPrev) THEN
            x[NX] := iMax; INC(NX);
          ELSIF (qThis < qNext AND qThis < qPrev) THEN
            x[NX] := iMin; INC(NX);
          END;
          j := i;
          qPrev := qThis;
          qThis := qNext;
          iMax := k;
          iMin := k
        ELSE
          IF c[k] < c[iMin] THEN iMin := k END;
          IF c[k] > c[iMax] THEN iMax := k END;
        END
      END;
      WITH ry = NEW(REF IntArray,NX) DO
        ry^ := SUBARRAY(x, 0, NX);
        RETURN ry
      END
    END
  END FindIndexExtrema;
  
<* UNUSED *>
PROCEDURE SampleCurvatureAtExtremaAndInterpolatedPoints(
    READONLY c : PZLRChain.T; 
    T : LONGREAL; 
    epsilon, delta: LONGREAL;
    idealStep : LONGREAL;
  ): REF PZLRChain.T =
  VAR 
    NS: CARDINAL := 0;
    d: LONGREAL;
  BEGIN
    WITH
      NC = NUMBER(c),
      inputStep = T / FLOAT(NC, LONGREAL),
      maxNS = NC + CEILING(T/idealStep),
      s = NEW(REF PZLRChain.T, maxNS)^,
      xIndex  = FindIndexExtrema(c, epsilon, delta)^,
      NX = NUMBER(xIndex)
    DO
      FOR j:= 0 TO NX-1 DO 
        WITH 
          i = xIndex[(j-1) MOD NX],
          k = xIndex[j MOD NX],
          ti = FLOAT(i,LONGREAL) * inputStep
        DO
          <* ASSERT  i # k *>
          IF k > i THEN 
            d := FLOAT((k-i),LONGREAL) * inputStep;
          ELSE
            d := FLOAT((k+NC-i),LONGREAL) * inputStep;
          END;
          (* Wr.PutText(stderr,"d:  " & Fmt.LongReal(d) & "\n"); *)
          IF d > idealStep THEN 
            WITH 
              r = DivideInterval(d, idealStep),
              actualStep = d /FLOAT(r,LONGREAL)
            DO
              FOR v := 1 TO r-1 DO 
                WITH
                  t = ti + FLOAT(v,LONGREAL)*FLOAT(actualStep, LONGREAL), 
                  alpha = TRUNC(t / inputStep),
                  a = alpha MOD NC,
                  ta = FLOAT(alpha,LONGREAL) * inputStep,
                  beta = alpha + 1,
                  b = beta MOD NC,
                  tb = FLOAT(beta,LONGREAL) * inputStep
                DO
                  s[NS] := Interpolate(t, ta, c[a], tb, c[b]);
                  INC(NS);
                END
              END;
            END;  
          END;
          s[NS] := c[k];  
          INC(NS);
        END;
      END;
      WITH ry = NEW(REF PZLRChain.T, NS) DO
        ry^ := SUBARRAY(s, 0, NS);
        RETURN ry;
      END;
    END;  
  END SampleCurvatureAtExtremaAndInterpolatedPoints;

PROCEDURE SampleCurvatureAtUniformlySpacedPoints( 
   READONLY c: PZLRChain.T; 
   T : LONGREAL; 
   idealStep : LONGREAL;
  ): REF PZLRChain.T =
  BEGIN 
    WITH 
      NC = NUMBER(c),
      inputStep = T / FLOAT(NC, LONGREAL),
      NS = DivideInterval(T, idealStep),
      actualStep = T / FLOAT(NS,LONGREAL),
      s = NEW(REF PZLRChain.T, NS) 
    DO
     FOR k:= 0 TO NS-1 DO 
       WITH
         t = FLOAT(k,LONGREAL) * actualStep, 
         i0 = TRUNC(t/inputStep), 
         t0 = FLOAT(i0,LONGREAL) * inputStep, 
         i1 = (i0+1) MOD NC, 
         t1 = t0 + inputStep
       DO 
         s[k] := Interpolate(t,t0,c[i0],t1,c[i1]);
       END;
     END;
     RETURN s;
   END;
 END SampleCurvatureAtUniformlySpacedPoints;

PROCEDURE DivideInterval(d, idealStep: LONGREAL): CARDINAL =
  BEGIN
    WITH
      r0 = d / idealStep,
      rLo = FLOAT(FLOOR(r0),LONGREAL),
      rHi = rLo + 1.0d0,
      logcr0 = Math.log(r0),
      wLo  = logcr0 - Math.log(rLo), wLo2 = wLo * wLo,
      wHi  = logcr0 - Math.log(rHi), wHi2 = wHi * wHi
    DO
      <* ASSERT r0 >= 0.0d0 *>
      <* ASSERT rLo >= 0.0d0 *>
      <* ASSERT rHi >= 1.0d0 *>
      IF wLo2 <= wHi2 THEN 
        RETURN ROUND(rLo)
      ELSE  
        RETURN ROUND(rHi)
      END;
    END
  END DivideInterval;


PROCEDURE SampleLR3ChainAtUniformlySpacedPoints( 
   READONLY c: PZLR3Chain.T; 
   T : LONGREAL; 
   idealStep : LONGREAL;
  ): REF PZLR3Chain.T =
  BEGIN 
    WITH 
      NC = NUMBER(c),
      inputStep = T / FLOAT(NC, LONGREAL),
      NS = DivideInterval(T, idealStep),
      actualStep = T / FLOAT(NS,LONGREAL),
      s = NEW(REF PZLR3Chain.T, NS) 
    DO
     FOR k:= 0 TO NS-1 DO 
       WITH
         t = FLOAT(k,LONGREAL) * actualStep, 
         i0 = TRUNC(t/inputStep), 
         t0 = FLOAT(i0,LONGREAL) * inputStep, 
         i1 = (i0+1) MOD NC, 
         t1 = t0 + inputStep
       DO 
         s[k] := PZGeo.LinearInterpolate(t,t0,c[i0],t1,c[i1]);
       END;
     END;
     RETURN s;
   END;
 END SampleLR3ChainAtUniformlySpacedPoints;

      idealStep : LONGREAL; (* Ideal spacing of curvature sampling points *)

        pp.getKeyword("-idealStep");
        o.idealStep := pp.getNextLongReal();
        
         \\\n");
          Wr.PutText(stderr, "  -idealStep    <number>
 
Wr.PutText( wr, "         idealStep: " & Fmt.LongReal(o.idealStep) & "\n");
      

(* From PZInvar.Main *)

      WITH
        uFlc = SampleLR3ChainAtUniformlySpacedPoints(p, T, o.idealStep)^, 
        uFlcCmt = cmt & "\n" & "LR3Chain - sampled at uniformly spaced points",
        lu = SampleCurvatureAtUniformlySpacedPoints(cv, T, o.idealStep)^,
        luCmt = cmt & "\n" & "Curvature - sampled at uniformly spaced points",
        cc =  QuantizeInvariantChain(flu, o.epsilon, o.delta)^,
        ccCmt = cmt & "\n" & "Coded Curvature - sampled at uniformly spaced points"
      DO
        PZLR3Chain.Write(FileWr.Open(o.outName & "-u.flc"),uFlcCmt,uFlc, o.idealStep); 
        PZLRChain.Write(FileWr.Open(o.outName & "-u.fcv"), fluCmt, flu, T,o.idealStep);
        PZSymbolChain.Write(
          FileWr.Open(o.outName & "-u.cvc"), 
          cvcCmt, cvc, T, o.epsilon, o.delta
        ); 
      END;


<* UNUSED *>  
PROCEDURE NormalizedCurvature(VAR c: PZLR3Chain.T; cMax : LONGREAL;)=
  BEGIN 
    WITH 
      NC = NUMBER(c)
    DO 
      IF cMax < 1.0d0 THEN cMax := 1.0d0 END;
      FOR i:=0 TO NC-1 DO 
        c[i][2] :=  c[i][2] / cMax;  
      END;
    END; 
  END NormalizedCurvature; 

 
<* UNUSED *>
PROCEDURE Curv(READONLY p1, p2, p3: LR3.T): LONGREAL =
  BEGIN 
    WITH 
      d = Dist(p1,p3)  
    DO  
      RETURN 2.0d0 * Det(p1, p2, p3)/(d*d*d) 
    END;
  END Curv;

PROCEDURE Det(READONLY a, b, c: LR3.T): LONGREAL =
  BEGIN 
    WITH 
      xa = a[0] - b[0],
      ya = a[1] - b[1],
      xc = c[0] - b[0],
      yc = c[1] - b[1]
    DO 
      RETURN xa*yc - xc*ya
    END
   END Det;


PROCEDURE Squeeze(x:LONGREAL; alpha: LONGREAL): LONGREAL = 
  BEGIN 
    WITH  z = 2.0d0 * alpha * x  
    DO 
      RETURN ((Math.exp(z) - 1.0d0) / (Math.exp(z) + 1.0d0))
    END
  END Squeeze;


PROCEDURE CurvChainFromInvariantChain(
    READONLY  c: FloatChain; 
    epsilon, delta: LONGREAL; 
  ): REF CurvChain =
  (*
    Given the invariant chain, converts it to a coded
    curvature chain.
  *)
  BEGIN 
    WITH 
      NP = NUMBER(c),
      rcvc = NEW(REF CurvChain, NP),
      cvc = rcvc^
    DO
     FOR i:= 0 TO NP-1 DO 
        WITH 
          h =  c[i][2]* squeeze
          hs = 26.5d0 * MIN(+1.0d0, MAX(-1.0d0, Squeeze(h,1.0d0))),
          cs = ROUND(hs)
        DO
          IF cs < 0 THEN
            cvc[i] :=  VAL(ORD('a') - cs - 1, CHAR);
          ELSIF cs > 0 THEN 
            cvc[i] :=  VAL(ORD('A') + cs - 1, CHAR);
          ELSE 
            cvc[i] :=  VAL(ORD('0'),CHAR);  
          END
        END
      END;
      RETURN cvc
    END 
  END CurvChainFromInvariantChain;


                (*
                  Wr.PutText(stderr, Fmt.Int(nSamples) & 
                    "  " & Fmt.LongReal(s[nSamples,0]) & 
                    "  " & Fmt.LongReal(s[nSamples,2]) & "\n"
                  );    
                *)

      IF (c[0,2] > c[NC-1, 2] AND c[0,2] > c[1,2]) 
      OR (c[0,2] < c[NC-1, 2] AND c[0,2] < c[1,2]) 
      THEN 
        t[nSamples] := c[0];
        INC(nSamples);
      END; 
      FOR i := 1 TO NC-1 DO 
        WITH
          a = c[(i-1) MOD NC, 2],
          b = c[i, 2],
          c = c[(i+1) MOD NC, 2]
        DO
          IF (b > a AND b > c) 
          OR (b < a AND b < c) 
          THEN
            t[nSamples] := c[i];
            INC(nSamples)
          END
        END
      END;
      WITH ry = NEW(REF FloatChain, nSamples) DO
        ry^ := SUBARRAY(t, 0, nSamples);
        RETURN ry
      END
    END

PROCEDURE GenerateInverseCurvChain(READONLY c: CurvChain): REF CurvChain =
  BEGIN
    WITH 
      n = NUMBER(c),
      cInv = NEW(REF CurvChain, n)
    DO 
      FOR i := 0 TO n-1  DO 
        IF c[i] = VAL(ORD('0'),CHAR) THEN cInv[n-i-1] := VAL(ORD('0'),CHAR) 
        ELSIF   c[i] >= VAL(ORD('A'), CHAR)  AND c[i] <= VAL(ORD('Z'), CHAR)THEN
             cInv[n-i-1] :=  VAL(ORD('a') + ORD(c[i]) - ORD('A'), CHAR);
        ELSE 
          cInv[n-i-1] :=  VAL(ORD('A') + ORD(c[i]) -  ORD('a'), CHAR);
        END;
      END;
      RETURN cInv;
    END;  
  END GenerateInverseCurvChain;
  
   
PROCEDURE CurvChainFromInvariantChain(
  READONLY  p: FloatChain; 
  squeeze : LONGREAL; 
  ): REF CurvChain =
  (*
    Given the invariant chain, converts it to a coded
    curvature chain.
  *)
  BEGIN 
    WITH 
      n = NUMBER(p),
      c = NEW(REF CurvChain, n)
    DO
     FOR i:= 0 TO n-1 DO 
        WITH 
          h  = p[i][2] * squeeze,
          hs = 26.5d0 * MIN(+1.0d0, MAX(-1.0d0, Squeeze(h,1.0d0))),
          cs = ROUND(hs)
        DO
          IF cs < 0 THEN
            c[i] :=  VAL(ORD('a') - cs - 1, CHAR);
          ELSIF cs > 0 THEN 
            c[i] :=  VAL(ORD('A') + cs - 1, CHAR);
          ELSE 
            c[i] :=  VAL(ORD('0'),CHAR);  
          END
        END
      END;
      RETURN c
    END 
  END CurvChainFromInvariantChain;