Unit Contours;
{
  Processing of BitMap Contours
  =============================
  This software has been designed and is CopyLefted
  by Han de Bruijn

   (===)
  @-O^O-@
   #/_\#
    ###

  All I ask is to be credited when it is appropriate.
}
INTERFACE

Uses SysUtils, Forms, Graphics, Generaal;

type
  schakel = record
    x,y : integer;
    lnk : integer;
  end;
  RGB = (R,G,B);
  afbeelding = array of array of double;

type
  Omtrekken = class
  private
    log : Nuttig;
    blaat : boolean;
  { Graphics Input } 
    Wijd,Hoog : integer;
  { Edges/Vertices of Contours: }
    bewaar : array of schakel;
    reserve : array of schakel;
    Geteld : integer;
  { Indices of Contours: }
    randen : array of integer;

    procedure Contouren(serieus : boolean);
    procedure Verzamelen;
    procedure Gesloten;
  public
    funktie : afbeelding;
    Scherm : TBitmap;
    constructor Create;
    procedure afmetingen(W,H : integer); { Transfer Width and Height }
    procedure Rondom(blah : boolean); { Overall processing }
    function Aantal : integer; { Number of Contours }
    function Totaal(const klaar : boolean = true) : integer; { Total Length }
    procedure Tekenen(welke : integer ; C : TColor); { Draw }
    procedure Schetsen(welke : integer ; C : TColor; W : integer); { Draw }
    function StartPunt(deze : integer) : punt; { a Vertex of Contour }
    function Kromme(deze : integer) : punten; { Retrieve Contour }
  { Reducing the number of vertices in a contour by Fan Data Compression }
    procedure FanDataCompressie(deze : integer; epsilon : double);
  end;

IMPLEMENTATION

constructor Omtrekken.Create;
begin
  log := Nuttig.Create;
end;

procedure Omtrekken.afmetingen(W,H : integer);
begin
  Wijd := W;
  Hoog := H;
  SetLength(funktie,W,H);
end;

function Omtrekken.Totaal;
begin
  if not klaar then Contouren(false);
  Totaal := Geteld;
end;

function Omtrekken.Aantal : integer;
begin
  Aantal := Length(randen)-1;
end;

procedure Omtrekken.Contouren(serieus : boolean);
{
  Find Contour Lines in a Black & White Image
  -------------------------------------------
  The gist of the algorithm is that there exists
  kind of a linear interpolation between function
  values 1 for black pixels and 0 for white pixels.

  Now consider vertices in between and search for
  function values 1/2, being the contouring level.

  The edges are oriented in such a way that if you
  walk from start- to endpoint, any black pixel is
  at your right (uphill) side and any white pixel
  is at your left (downhill) side. Contours never
  intersect.
}
var
  x,y,x1,y1,x2,y2 : integer;
  getal,tel : integer;
  extra : string;
const                   { getal =  0 1 2 3 4 5 6 7 8 9 A B C D E F }
  dx1 : array[0..$F] of integer = (9,2,1,2,1,1,9,1,0,9,1,2,0,0,1,9);
  dy1 : array[0..$F] of integer = (9,1,2,1,0,0,9,0,1,9,2,1,1,1,2,9);
  dx2 : array[0..$F] of integer = (9,1,0,0,2,1,9,0,1,9,1,1,2,1,2,9);
  dy2 : array[0..$F] of integer = (9,2,1,1,1,2,9,1,0,9,0,0,1,2,1,9);
                        { geval =  0 1 1 3 1 3 4 2 1 4 3 2 3 2 2 0 }
procedure Opslaan;
{
  Storage of Contouring Edges
  where (x1,y1) = begin / odd
        (x2,y2) = end / even.
}
begin
{ Odd counts }
  bewaar[tel-1].x := x1-2;
  bewaar[tel-1].y := y1-2;
{ Even counts }
  bewaar[tel].x := x2-2;
  bewaar[tel].y := y2-2;
end;

function getPixel(x,y : integer) : boolean;
{
  Dark Pixel
}
begin
  getPixel := false;
  if x < 0 then Exit;
  if x > Wijd-1 then Exit;
  if y < 0 then Exit;
  if y > Hoog-1 then Exit;
  getPixel := (funktie[x,y] > 0);
end;

function MetVieren(x,y : integer) : integer;
{
  Take four pixels at a time and    0 1
  caracterize their B/W structure   2 3
  by a number called 'MetVieren'.
}
var
  getal, k : integer;
  zwart : boolean;
const
  ix : array[1..4] of integer = ( 0, 1, 0, 1);
  iy : array[1..4] of integer = ( 0, 0, 1, 1);
begin
  getal := 0;
  for k := 1 to 4 do
  begin
    getal := getal shl 1;
    zwart := getPixel(x+ix[k]-1,y+iy[k]-1);
    if zwart then getal := getal + 1;
  end;
  MetVieren := getal;
end;
{ ======================== }
{ MAIN routine starts here }
{ ======================== }
begin
  if blaat then log.TijdMeting('');
  tel := 0;
  for y := 0 to Hoog do
  begin
    Application.ProcessMessages;

    for x := 0 to Wijd do
    begin
    { Take four pixels at a time and    0 1
      caracterize their B/W structure   2 3
      by a number called 'getal': }

      getal := MetVieren(x,y);

    { All coordinates are multipied by a factor 2     0 1
      in order to preserve their "integer" nature.    2 3
      In comments below  O = White , X = Black : }

      Case getal of

        $0,$F: Continue; { Nothing to do = geval 0 }

      { One black pixel of four gives an edge which joins   0 1
        the two points midway between the one black pixel   2 3
        and two white pixels: }

        $1,$2,$4,$8,  { geval 1 }
        { X O   O X   O O   O O }
        { O O   O O   X O   O X }

      { Three black pixels of four gives an edge joining    0 1
        the two points midway between the one white pixel   2 3
        and two black pixels: }

        $D,$E,$B,$7,  { geval 2 }
        { X O   O X   X X   X X }
        { X X   X X   O X   X O }

        { Two black pixels of four gives an edge which joins  0 1
          the two points midway between the black pixels and  2 3
          the white pixels: }

        $3,$C,$A,$5:  { geval 3 }
        { X X   O O   O X   X O }
        { O O   X X   O X   X O }

        begin
          tel := tel + 2;
          if serieus then
          begin
            x1 := 2*x + dx1[getal];
            y1 := 2*y + dy1[getal];
            x2 := 2*x + dx2[getal];
            y2 := 2*y + dy2[getal];
            Opslaan;
          end;
        end;

        { Two black pixels of four gives two edges joining
          four points midway between the black pixels and    0 1
          the white pixels in such a way that both edges     2 3
          lie inside a white area. Meaning that the other
          possibility is deliberately excluded: it's Black
          on White, it's Not White on Black here ! }

        $9: { geval 4 }
        begin
          tel := tel + 2;
          if serieus then
          begin
            x1 := 2*x + 2 ; y1 := 2*y + 1 ;  { O X }
            x2 := 2*x + 1 ; y2 := 2*y + 0 ;  { X O }
            Opslaan;
          end;
          tel := tel + 2;
          if serieus then
          begin
            x1 := 2*x + 0 ; y1 := 2*y + 1 ;
            x2 := 2*x + 1 ; y2 := 2*y + 2 ;
            Opslaan;
          end;
        end;
        $6: { geval 4 }
        begin
          tel := tel + 2;
          if serieus then
          begin
            x1 := 2*x + 1 ; y1 := 2*y + 0 ;  { X O }
            x2 := 2*x + 0 ; y2 := 2*y + 1 ;  { O X }
            Opslaan;
          end;
          tel := tel + 2;
          if serieus then
          begin
            x1 := 2*x + 1 ; y1 := 2*y + 2 ;
            x2 := 2*x + 2 ; y2 := 2*y + 1 ;
            Opslaan;
          end;
        end;
      end;
    end;
  end;

  Geteld := tel div 2;
  extra := '';
  if serieus then extra := 'I';
  if blaat then log.TijdMeting('Contouring pass I' + extra);
end;

procedure Omtrekken.Verzamelen;
{
  Collect Contouring Edges
  to form Closed Curves
}
var
  start,einde : schakel;
  p,q, k : integer;
  pq,qp : boolean;
begin
  p := 0;
  while p < Geteld do
  begin
    Application.ProcessMessages;
    p := p + 1;
    if (bewaar[2*p-1].lnk > 0)
     and (bewaar[2*p].lnk > 0) then Continue;
    q := p;
    while (q <= Geteld) do
    begin
      q := q + 1;
      einde := bewaar[2*p];
      start := bewaar[2*q-1];
      pq := (start.x = einde.x)
        and (start.y = einde.y);
      if pq then
      begin
        bewaar[2*p].lnk := q;
        bewaar[2*q-1].lnk := p;
      end;
      einde := bewaar[2*q];
      start := bewaar[2*p-1];
      qp := (start.x = einde.x)
        and (start.y = einde.y);
      if qp then
      begin
        bewaar[2*q].lnk := p;
        bewaar[2*p-1].lnk := q;
      end;
      if (bewaar[2*p-1].lnk = 0)
        or (bewaar[2*p].lnk = 0) then Continue;
      if (pq or qp) then Break;
    end;
  end;

{ Keep only half of the data }
  SetLength(reserve,2*Geteld);
  for k := 1 to Geteld do
  begin
    reserve[2*k-2] := bewaar[2*k-1]; 
    reserve[2*k-1] := bewaar[2*k]; 
  end;
  for k := 1 to Geteld do
    bewaar[k] := bewaar[2*k]; 
  SetLength(bewaar,Geteld+1);
end;

procedure Omtrekken.Gesloten;
{
  Closure of Contours
}
var
  done : array of boolean;
  tel, k,L : integer;
begin
  SetLength(done,Geteld+1);
  for k := 1 to Geteld do
    done[k] := false;

  tel := 0;
  for k := 1 to Geteld do
  begin
    if done[k] then Continue;
    tel := tel + 1;
    L := k;
    while not done[L] do
    begin
      done[L] := true;
      L := bewaar[L].lnk;
    end;
  end;

  SetLength(randen,tel+1);
  for k := 1 to Geteld do
    done[k] := false;

  tel := 0;
  for k := 1 to Geteld do
  begin
    if done[k] then Continue;
    tel := tel + 1;
    randen[tel] := k;
    L := k;
    while not done[L] do
    begin
      done[L] := true;
      L := bewaar[L].lnk;
    end;
  end;

  SetLength(done,0);
end;

procedure ZamelTest(OM : Omtrekken);
{
  Just a Test
}
var
  k : integer;
  ff : TextFile;
  w : schakel;
begin
  AssignFile(ff,'bewaar.txt');
  Rewrite(ff);
  for k := 1 to OM.Geteld do
  begin
    w := OM.bewaar[k];
    Writeln(ff,w.x:6,w.y:6,w.lnk:6);
  end;
  CloseFile(ff);
end;

procedure SluitTest(OM : Omtrekken);
{
  Just a Test
}
var
  k,aantal : integer;
  ff : TextFile;
begin
  AssignFile(ff,'gesloten.txt');
  Rewrite(ff);
  aantal := Length(OM.randen)-1;
  for k := 1 to aantal do
    Writeln(ff,OM.randen[k]:5);
  CloseFile(ff);
end;

procedure Omtrekken.Rondom(blah : boolean);
{
  All Together
}
var
  genoeg,k : integer;
begin
  blaat := blah;
  Contouren(false);
  if blah then
    Writeln('# contour segments = ' + IntToStr(Geteld));

  genoeg := 2*Geteld;
  SetLength(bewaar, genoeg+1);
  for k := 1 to genoeg do
  begin
    bewaar[k].x := 0;
    bewaar[k].y := 0;
    bewaar[k].lnk := 0;
  end;

  Contouren(true);
  Verzamelen;
  Gesloten;
end;

function Omtrekken.Kromme(deze : integer) : punten;
{
  Define Contour
}
var
  k,L,tel : integer;
  rond : punten;
begin
  k := randen[deze];
  L := k;
  tel := 0;
  while true do
  begin
    L := bewaar[L].lnk;
    tel := tel + 1;
    if k = L then Break;
  end;
{ Number of Vertices }
  SetLength(rond,tel);

  L := k;
  tel := 0;
  while true do
  begin
    rond[tel].x := bewaar[L].x/2;
    rond[tel].y := bewaar[L].y/2;
    L := bewaar[L].lnk;
    tel := tel + 1;
    if k = L then Break;
  end;
  Kromme := rond;
  SetLength(rond,0);
end;

procedure Omtrekken.Tekenen(welke : integer; C : TColor);
{
  Do the Drawing of a
  fine grained Contour
}
var
  L,k, x,y, xe,ye : integer;
  x1,y1,x2,y2 : double;

  function links(xp,yp : integer) : boolean;
  {
    At the Left of an Edge (1) => (2)
  }
  var
    xd,yd : double;
  begin
    xd := xp; yd := yp;
    links := (y2-y1)*(xd-x1)-(x2-x1)*(yd-y1) < 0;
  end;

begin
  k := randen[welke];
  x2 := bewaar[k].x/2 ; y2 := bewaar[k].y/2;
  L := k;
  while true do
  begin
    L := bewaar[L].lnk;
    x1 := x2 ; y1 := y2 ;
    x2 := bewaar[L].x/2 ; y2 := bewaar[L].y/2;
    x := bewaar[L].x ; y := bewaar[L].y ;
    if (abs(x mod 2) = 1) and ((y mod 2) = 0) then
    begin
      xe := (x-1) div 2 ; ye := y div 2 ;
      With Scherm.Canvas do
      if links(xe,ye) then Pixels[xe,ye] := C;
      xe := (x+1) div 2 ; ye := y div 2 ;
      With Scherm.Canvas do
      if links(xe,ye) then Pixels[xe,ye] := C;
    end;
    if ((x mod 2) = 0) and (abs(y mod 2) = 1) then
    begin
      xe := x div 2 ; ye := (y-1) div 2 ;
      With Scherm.Canvas do
      if links(xe,ye) then Pixels[xe,ye] := C;
      xe := x div 2 ; ye := (y+1) div 2 ;
      With Scherm.Canvas do
      if links(xe,ye) then Pixels[xe,ye] := C;
    end;
    if L = k then Break;
  end;
end;

procedure Omtrekken.Schetsen(welke : integer; C : TColor; W : integer);
{
  Do the Drawing of a
  Compressed Contour
}
var
  L,k, x,y : integer;
begin
  k := randen[welke];
  x := bewaar[k].x div 2;
  y := bewaar[k].y div 2;
  Scherm.Canvas.MoveTo(x,y);
  Scherm.Canvas.Pen.Color := C;
  Scherm.Canvas.Pen.Width := W;
  L := k;
  while true do
  begin
    L := bewaar[L].lnk;
    x := bewaar[L].x div 2;
    y := bewaar[L].y div 2;
    Scherm.Canvas.LineTo(x,y);
    if L = k then Break;
  end;
end;

procedure KansloosMaken(OM : Omtrekken; deze : integer);
{
  Disable Contour
}
var
  ingang : integer;
begin
  ingang := OM.randen[deze];
  OM.bewaar[ingang].lnk := ingang;
end;

function isKansloos(OM : Omtrekken; deze : integer) : boolean;
{
  Disabled Contour
}
var
  ingang : integer;
begin
  ingang := OM.randen[deze];
  isKansloos := (OM.bewaar[ingang].lnk = ingang);
end;

function Omtrekken.StartPunt(deze : integer) : punt;
{
  Pixel near Start of Contour
}
var
  ingang : integer;
  x1,y1,x2,y2 : double;
  L,x,y,xe,ye : integer;
  p : punt;

  function rechts(xp,yp : integer) : boolean;
  {
    At the Right of an Edge (1) => (2)
  }
  var
    xd,yd : double;
  begin
    xd := xp; yd := yp;
    rechts := (y2-y1)*(xd-x1)-(x2-x1)*(yd-y1) > 0;
  end;

begin
  ingang := randen[deze];
  x1 := bewaar[ingang].x/2;
  y1 := bewaar[ingang].y/2;
  L := bewaar[ingang].lnk;
  x2 := bewaar[L].x/2;
  y2 := bewaar[L].y/2;
  x := bewaar[L].x ; y := bewaar[L].y ;
  if (abs(x mod 2) = 1) and ((y mod 2) = 0) then
  begin
    xe := (x-1) div 2 ; ye := y div 2 ;
    if rechts(xe,ye) then
    begin
      p.x := xe; p.y := ye;
    end;
    xe := (x+1) div 2 ; ye := y div 2 ;
    if rechts(xe,ye) then
    begin
      p.x := xe; p.y := ye;
    end;
  end;
  if ((x mod 2) = 0) and (abs(y mod 2) = 1) then
  begin
    xe := x div 2 ; ye := (y-1) div 2 ;
    if rechts(xe,ye) then
    begin
      p.x := xe; p.y := ye;
    end;
    xe := x div 2 ; ye := (y+1) div 2 ;
    if rechts(xe,ye) then
    begin
      p.x := xe; p.y := ye;
    end;
  end;
  StartPunt := p;
end;

procedure Omtrekken.FanDataCompressie(deze : integer; epsilon : double);
{
  Fan Compression algorithm
  -------------------------
  IEEE Computer Graphics & Applications
  "Faster Plots by Fan Data-Compression"
  Richard A. Fowell and David D. McNeil
  March 1989.

  eps = desired accuracy
}
var
  L,k,m,kept,tel, memo : integer;
  pu,pv, dist, h,dh, eps : double;
  keep,verder : boolean;
  tmp,o,p : punt;
  u1,u2,f1,f2,f3 : double;
begin
  eps := epsilon;
  kept := 0;
  tel := 0;
  k := randen[deze];
  L := k;
  o.x := bewaar[L].x/2;
  o.y := bewaar[L].y/2;
  verder := true;
  memo := L;
  u1 := 0 ; u2 := 0 ;
  f1 := 0 ; f2 := 0 ; f3 := 0 ;

  while true do
  begin
    L := bewaar[L].lnk;
    if L = k then Break;
    tel := tel + 1;
  { Local coordinates }
    if verder then
    begin
      p.x := bewaar[L].x/2 - o.x;
      p.y := bewaar[L].y/2 - o.y;
      dist := sqrt(p.x*p.x+p.y*p.y);
      if dist > eps then
      begin
        verder := false;
        u1 := p.x/dist;
        u2 := p.y/dist;
        f1 := dist;
        f2 := eps/dist;
        f3 := -eps/dist;
      end;
    end;

  { Elimination process }
    if not verder then
    begin
      keep := true;
      m := bewaar[L].lnk;
      tmp.x := bewaar[m].x/2;
      tmp.y := bewaar[m].y/2;
      pu :=  (tmp.x-o.x)*u1+(tmp.y-o.y)*u2;
      pv := -(tmp.x-o.x)*u2+(tmp.y-o.y)*u1;
      if pu >= f1 then
      begin
        h := pv/pu;
        if (h <= f2) and (h >= f3) then
        begin
          dh := eps/pu;
          keep := false;
          f1 := pu;
          if h+dh < f2 then f2 := h+dh;
          if h-dh > f3 then f3 := h-dh;
        end;
      end;
      if keep then
      begin
        kept := kept + 1;
        o.x := bewaar[L].x/2;
        o.y := bewaar[L].y/2;
        bewaar[memo].lnk := L;
        memo := L;
        verder := true;
      end;
    end;
  end;
  bewaar[memo].lnk := k;
  o.x := bewaar[k].x/2;
  o.y := bewaar[k].y/2;
  Writeln('Left after Compression = ',Round(100*kept/tel),' %');
end;

END.