MODULE PhiFeatures2D;

IMPORT Wr, Thread;
IMPORT FeatureMapper2D AS Mapper;
FROM Mapper IMPORT
  NAT, LONG, Nhood, Nhoods, CoNhoods, Features, CoFeatures,
  EvalSpaceDerivatives;

REVEAL
  T = Public BRANDED OBJECT
    OVERRIDES
      init := Init;
      print := Print;
      eval := Eval;
      diff := Diff;
      defined := Defined;
    END; 

PROCEDURE Init(s: T; NL: NAT): T =
  BEGIN
    s.NL := NL;
    s.NF := 5*NL;
    RETURN s
  END Init;

PROCEDURE Print(<*UNUSED*> s: T; wr: Wr.T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr, "Phi: the phi invariants\n");
    Wr.Flush(wr)
  END Print;

PROCEDURE Eval(s: T; READONLY n: Nhoods; VAR t: Features) =
  BEGIN
    FOR i := 0 TO s.NL-1 DO
      WITH
        n = n[i],
        t = SUBARRAY(t, 5*i, 5),
        D = EvalSpaceDerivatives(n),
        fx2 = D.fx*D.fx,
        fy2 = D.fy*D.fy,
        g = FLOAT(n[+1,00] + n[-1,00] - n[00,+1] - n[00,-1], LONG),
        h = 0.5d0 * FLOAT(n[+1,+1] + n[-1,-1] - n[+1,-1] - n[-1,+1], LONG)
      DO
        t[0] := FLOAT(n[0,0], LONG);
        t[1] := fx2 + fy2;
        t[2] := D.fxx + D.fyy;
        t[3] := D.fxx*fx2 + h*D.fx*D.fy + D.fyy*fy2;
        t[4] := h*h + g*g
      END
    END
  END Eval;

PROCEDURE Diff(
    s: T;
    READONLY dEdt: CoFeatures;
    READONLY t: Features;
    READONLY n: Nhoods;
    VAR dEdn: CoNhoods;
  ) =
  BEGIN
    FOR i := 0 TO s.NL-1 DO
      WITH
        dEdt = SUBARRAY(dEdt, 5*i, 5),
        n = n[i],
        dEdn = dEdn[i],
        t = SUBARRAY(t, 5*i, 5),
        D = EvalSpaceDerivatives(n),
        g = FLOAT(n[+1,00] + n[-1,00] - n[00,+1] - n[00,-1], LONG),
        h = 0.5d0 * FLOAT(n[+1,+1] + n[-1,-1] - n[+1,-1] - n[-1,+1], LONG),

        hfx2 = 0.5d0 * h * D.fx,
        hfy2 = 0.5d0 * h * D.fy,
        hfxfy = 0.5d0 * D.fx * D.fy,
        fxdEdt1 = D.fx * dEdt[1],
        fydEdt1 = D.fy * dEdt[1],
        hfxfydEdt3 = hfxfy * dEdt[3],
        g2dEdt4 = 2.0d0 * g * dEdt[4],
        hdEdt4 = h * dEdt[4]
      DO
        dEdn[00,00] := FLOAT(dEdt[0] - 4.0d0 * dEdt[2] - 2.0d0 * t[1]*dEdt[3]);
        dEdn[-1,00] := FLOAT(- fydEdt1 + dEdt[2] + (D.fy*(D.fy-D.fyy)-hfx2)*dEdt[3] + g2dEdt4);
        dEdn[+1,00] := FLOAT(+ fydEdt1 + dEdt[2] + (D.fy*(D.fy+D.fyy)+hfx2)*dEdt[3] + g2dEdt4);
        dEdn[00,-1] := FLOAT(- fxdEdt1 + dEdt[2] + (D.fx*(D.fx-D.fxx)-hfy2)*dEdt[3] - g2dEdt4);
        dEdn[00,+1] := FLOAT(+ fxdEdt1 + dEdt[2] + (D.fx*(D.fx+D.fxx)+hfy2)*dEdt[3] - g2dEdt4);
        dEdn[-1,-1] := FLOAT(+ hfxfydEdt3 + hdEdt4);
        dEdn[+1,-1] := FLOAT(- hfxfydEdt3 - hdEdt4);
        dEdn[-1,+1] := FLOAT(- hfxfydEdt3 - hdEdt4);
        dEdn[+1,+1] := FLOAT(+ hfxfydEdt3 + hdEdt4);
      END
    END
  END Diff;
  
PROCEDURE Defined(
    <*UNUSED*> s: T;
    READONLY m: Nhood
  ): BOOLEAN =
  BEGIN
    RETURN
      m[-1,-1] > 0.5 AND
      m[-1,00] > 0.5 AND
      m[-1,+1] > 0.5 AND
      m[00,-1] > 0.5 AND
      m[00,00] > 0.5 AND
      m[00,+1] > 0.5 AND
      m[+1,-1] > 0.5 AND
      m[+1,00] > 0.5 AND
      m[+1,+1] > 0.5
  END Defined;

BEGIN
END PhiFeatures2D.