PROCEDURE LocalAttributesFromPixelNeghborhood(
    READONLY n: Nbr;
    VAR p: LocalAttributes.T;
  ) =
  VAR z: ARRAY [0..Pixel.N-1] OF RealNeighborhood;
  BEGIN
    SeparateNbrIntoLayers(n, z);
    (* Compute invariants: *)
    FOR k := 0 TO Pixel.N-1 DO
      Invariants.InvariantsFromNeighborhood(z[k], p[k]);
    END;
  END LocalAttributesFromPixelNeghborhood;

PROCEDURE LocalEnergyFromNeighborhood(
    READONLY n: Nbr;        (* The full neighborhood of a texel *)
    READONLY a: Attractor;                (* The reference charges *)
    scale: LONGREAL;                      (* Energy scaling factor *)
    VAR (*OUT*) E: LONGREAL;              (* Energy associated with this texel *)
    VAR (*OUT*) dEdn: Nbr;  (* Partials "dE/dn[k,y,x]" *)
    invP: InvariantsFromDerivativesProc;  (* Maps derivatives to invariants *)
  ) = 
  VAR z: ARRAY [0..Pixel.N-1] OF RealNeighborhood;
      v: PixelDerivatives;
      dEdv: PixelDerivatives;
      dEdz: ARRAY [0..Pixel.N-1] OF RealNeighborhood;
  BEGIN
    SeparateNbrIntoLayers(n, z);

    (* Compute Derivatives for all layers: *)
    FOR k := 0 TO Pixel.N-1 DO 
      DerivativesFromNeighborhood(z[k], v[k])
    END;
    
    (* Compute forces on invariants: *)
    LocalEnergyFromDerivatives(v, a, scale, E, dEdv);
    
    (* Convert partials relative to "v" into partials relative to "z": *)
    FOR k := 0 TO Pixel.N-1 DO 
      NeighborhoodPartialsFromDerivativePartials(dEdv[k], dEdz[k])
    END;
    
    PackLayersIntoNbr(dEdz, dEdn);
  END LocalEnergyFromNeighborhood;
  
PROCEDURE LocalEnergyFromDerivatives(
    READONLY v: PixelDerivatives;       (* Spatial derivatives of texture *)
    READONLY a: Attractor;              (* The reference charges *)
    scale: LONGREAL;                    (* Energy scaling factor *)
    VAR (*OUT*) E: LONGREAL;            (* Energy associated with this texel *)
    VAR (*OUT*) dEdv: PixelDerivatives; (* Partials dE/dv[k][i]  *)
  ) = 
  VAR p: LocalAttributes.T;
      dEdp: LocalAttributes.T;
  BEGIN
    (* Compute invariants for all layers: *)
    FOR k := 0 TO Pixel.N-1 DO 
      InvariantsFromDerivatives(v[k], p[k])
    END;
    
    (* Compute forces on invariants: *)
    LocalEnergyFromInvariants(p, a, scale, E, dEdp);
    
    (* Transform them into forces on texel components: *)
    FOR k := 0 TO Pixel.N-1 DO 
      DerivativePartialsFromInvariantPartials(v[k], dEdp[k], dEdv[k])
    END;
  END LocalEnergyFromDerivatives;

      fmv: LONGREAL;           (* = FLOAT(maxVal, LONGREAL) *)
      shift, scale: LONGREAL;  (* Pixel maping: pixel = round((state + shift)*scale) *)
    lo: REAL := 0.0; hi: REAL := 1.0;
    f.h.fmv := FLOAT(f.h.maxVal, LONGREAL);
    <* ASSERT lo # hi *>
    WITH 
      llo = FLOAT(lo, LONGREAL),
      hhi = FLOAT(hi, LONGREAL)
    DO
      f.h.shift := - llo; 
      f.h.scale := f.h.fmv / (hhi - llo)
    END;
    WITH
        ss = FLOAT(s[i], LONGREAL),
        fv = MAX(0.0d0, MIN(f.h.fmv, (ss + f.h.shift) * f.h.scale)),
        v = ROUND(fv)
      DO
    f.h.fmv := FLOAT(f.h.maxVal, LONGREAL);
    <* ASSERT lo # hi *>
    WITH 
      llo = FLOAT(lo, LONGREAL),
      hhi = FLOAT(hi, LONGREAL)
    DO
      f.h.scale := f.h.fmv / (hhi - llo);
      f.h.shift := - llo
    END; 
      s[i] := FLOAT(FLOAT(v, LONGREAL)/f.h.scale - f.h.shift)
        
PROCEDURE Read(
    fileName: TEXT; 
    lo: REAL := 0.0;
    hi: REAL := 1.0;
  ): REF ARRAY OF ARRAY OF Color.T =
  VAR NX, NY, maxVal: CARDINAL;
  <* FATAL Thread.Alerted, Rd.Failure, OSError.E *>
  BEGIN 
    WITH
      rd = FileRd.Open(fileName),
      rdr = NEW(Reader).init(rd, NX, NY, maxVal, lo, hi),
      mp = NEW(REF ARRAY OF ARRAY OF Color.T, NX, NY),
      m = mp^
    DO
      FOR y := NY-1 TO 0 BY -1 DO
        WITH my = m[y] DO
          FOR x := 0 TO NX-1 DO
            my[x] := rdr.get()
          END
        END
      END;
      rdr.finish();
      Rd.Close(rd);
      RETURN mp
    END;
  END Read;

PROCEDURE Write(
    fileName: TEXT; 
    READONLY m: ARRAY OF ARRAY OF Color.T;
    lo: REAL := 0.0;
    hi: REAL := 1.0;
    maxVal: CARDINAL := 255
  ) = 
  <* FATAL Thread.Alerted, Wr.Failure, OSError.E *>
  BEGIN
    WITH
      NY = NUMBER(m),
      NX = NUMBER(m[0]),
      wr = FileWr.Open(fileName),
      wrt = NEW(Writer).init(wr, NX, NY, lo, hi, maxVal, binary := TRUE)
    DO
      FOR y := NY-1 TO 0 BY -1 DO
        WITH my = m[y] DO
          FOR x := 0 TO NX-1 DO
            wrt.put(my[x])
          END
        END
      END;
      wrt.finish();
      Wr.Close(wr)
     END
   END Write;