Unit BMPunit;
{
  PunchCard Reader
  ================
  This program has been designed and is CopyLefted by:

* Han de Bruijn; Systems for Research "A little bit of Physics would be  (===)
* and Education (DTO/SOO), Mekelweg 6  NO Idleness in Mathematics"(HdB) @-O^O-@
* 2628 CD Delft, The Netherlands  http://huizen.dto.tudelft.nl/deBruijn  #/_\#
* E-mail: Han.deBruijn@DTO.TUDelft.NL  Tel: +31 15 27 82751. Fax: 81722   ###

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

  Contouring of Black and White BitMap
  ------------------------------------
}
INTERFACE

Uses SysUtils, Forms, Algemeen, Graphics;

var
{ Width and Height }
  Wijd, Hoog : integer;
{ Resolution in pixels / millimeter }
  PerMM : double;

function Inlezen(naam : string) : boolean; { Read BMP file }
function LaatZien : boolean; { Process BMP Header data }
procedure BMPschoon; { Clear variable arrays }
procedure BMProutines; { Overall processing }

procedure Schrijven(naam : string); { Write BMP file }

procedure KansloosMaken(deze : integer); { Disable Contour }
function isKansloos(deze : integer) : boolean; { Contour disabled ? }
function StartPunt(deze : integer) : vektor; { Vertex of Contour }

procedure Tekenen(Scherm : TBitmap; welk : integer ; C : TColor); { Draw }

function Kromme(deze : integer) : Jordan; { Contour Definition }
function Huidige : integer; { Contour Key }
function Aantal : integer; { Number of Contours }
{
Reducing the number of vertices in a contour }
procedure FanDataCompressie(deze : integer); { Fan Data Compression }

IMPLEMENTATION

const
  Geteld : integer = 0;
  tijdelijk : integer = 0;
type
  schakel = record
    x,y : integer;
    lnk : integer;
  end;
var
{ Edges/Vertices of Contours: }
  bewaar : array of schakel;
{ Indices of Contours: }
  randen : array of integer;
var
{ File stored in memory }
  Heel : array of byte;
  Groot : integer;
  Start : integer;
  zoveel : integer;
{ # of Black Pixels: }
  donker : integer;

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

function Inlezen(naam : string) : boolean;
{
  Read BMP file
}
const
  debug : boolean = false;
var
  ff : file;
  gelezen, L, t : integer;
  Lees : array[1..2048] of byte;
begin
  Inlezen := false;
  if not FileExists(naam) then
  begin
    Writeln(naam + ': does Not exist');
    Exit;
  end;
  Assign(ff,naam);
  Reset(ff,1);
  Groot := FileSize(ff);
  if debug then Writeln(Groot);
  SetLength(Heel, Groot);
  t := 0;
  repeat
    BlockRead(ff,Lees,SizeOf(Lees),gelezen);
    for L := 1 to gelezen do
    begin
      t := t + 1;
      Heel[t-1] := Lees[L];
    end;
  until (gelezen = 0) or (gelezen <> SizeOf(Lees));
  Close(ff);
  Inlezen := true;
end;

procedure Schrijven(naam : string);
{
  Write BMP file
}
const
  debug : boolean = false;
var
  ff : file;
  totaal, L,t,S : integer;
  Schrijf : array[1..2048] of byte;
begin
  S := SizeOf(Schrijf);
  Assign(ff,naam);
  Rewrite(ff,1);
  t := 0;
  repeat
    totaal := Groot-t;
    if totaal > S then totaal := S;
    for L := 1 to totaal do
    begin
      t := t + 1;
      Schrijf[L] := Heel[t-1];
    end;
    BlockWrite(ff,Schrijf,totaal);
  until (totaal < S);
  Close(ff);
end;

function getGetal(k1,k2 : integer) : integer;
{
  Extract BMP Number
}
var 
  lengte : integer;
  k : integer;
begin
  lengte := 0;
  for k := k2 downto (k1+1) do
  begin
    lengte := lengte*256 + Heel[k];
  end;
  getGetal := lengte;
end;

procedure setGetal(k1,k2,g : integer);
{
  Submit BMP Number
}
var 
  lengte : integer;
  k : integer;
begin
  lengte := g;
  for k := k1+1 to k2 do
  begin
    Heel[k] := (lengte mod 256);
    lengte := lengte div 256;
  end;
end;

procedure Gehakt; 
{ 
  Correcting a Resolution error
}
var
  O : integer;
begin
  Inlezen('fortran4.bmp');
  O := 14-1;
  setGetal(O+24,O+28,3937);
  setGetal(O+28,O+32,3937);
  Schrijven('fortran4.bmp');
end;

function GetPixel(x,y : integer) : boolean;
{
  Black Pixel from BMP
}
var
  waar : integer;
  beetje,bit : byte;
begin
  GetPixel := false;
  if x > Wijd-1 then Exit;
  if y > Hoog-1 then Exit;

{ localize byte in BMP file }
  waar := Start + (Hoog-1 - y) * zoveel + (x div 8);

  beetje := Heel[waar];
{ localize bit }
  bit := (beetje shr (7 - (x mod 8))) and 1;
  GetPixel := (bit = 0);
end;

function LaatZien;
{
  Show BMP Header
}
const 
  debug : boolean = false;
var
  hand : string;
  O : integer;
  zwart,wit : integer;
  OK : boolean;
  Horizontaal, Vertikaal : integer;
begin
  if debug then
    Writeln('Main Header:');
  SetLength(hand,2);
  hand[1] := char(Heel[0]);
  hand[2] := char(Heel[1]);
  if debug then
    Writeln(hand,' = BM');
  OK := (hand = 'BM');
  SetLength(hand,0);

  O := -1; { starts at 0 instead of 1 }
  if debug then
    Writeln(getGetal(O+2,O+6),' = ',Groot);
  if debug then
    Writeln(getGetal(O+6,O+10),' = 0');
  OK := OK and (getGetal(O+2,O+6) = Groot);
  if debug then
    Writeln('Offset = ',getGetal(O+10,O+14));
  Start := getGetal(O+10,O+14);
  if debug then Writeln;

  O := O + 14;
  if debug then
    Writeln('Info Header:');
  if debug then
    Writeln('40 = ',getGetal(O+0,O+4));
  OK := OK and (40 = getGetal(O+0,O+4));

  Writeln('Width = ',getGetal(O+4,O+8));
  Wijd := getGetal(O+4,O+8);
  Writeln('Height = ',getGetal(O+8,O+12));
  Hoog := getGetal(O+8,O+12);

{ # of bytes in scanline }
  zoveel := ((Wijd-1) div 8) + 1;
{ with padding up to 32bit boundary }
  zoveel := (((zoveel-1) div 4) + 1) * 4;

  if debug then
    Writeln('1 = ',getGetal(O+12,O+14));
  OK := OK and (1 = getGetal(O+12,O+14));
  if debug then
    Writeln('1 = ',getGetal(O+14,O+16));
  OK := OK and (1 = getGetal(O+14,O+16));
  if debug then
    Writeln('0 = ',getGetal(O+16,O+20));
  OK := OK and (0 = getGetal(O+16,O+20));

  Horizontaal := getGetal(O+24,O+28);
  Writeln('Horizontal Resolution = ',Horizontaal);
  Vertikaal := getGetal(O+28,O+32);
  Writeln('Vertical Resolution = ',Vertikaal);
  OK := OK and (Horizontaal = Vertikaal);
  PerMM := Horizontaal/1000;

  if debug then
  begin
    Writeln('? = ',getGetal(O+20,O+24));
    Writeln('2/0 = ',getGetal(O+32,O+36));
    Writeln('2/0 = ',getGetal(O+36,O+40));
    Writeln;
  end;

  O := O + 40;
  if debug then
  begin
    Writeln('Color Table:');
  { Black }
    Writeln('0 = ',getGetal(O+0,O+1));
    Writeln('0 = ',getGetal(O+1,O+2));
    Writeln('0 = ',getGetal(O+2,O+3));
    Writeln('0 = ',getGetal(O+3,O+4));
  { White }
    Writeln('255 = ',getGetal(O+4,O+5));
    Writeln('255 = ',getGetal(O+5,O+6));
    Writeln('255 = ',getGetal(O+6,O+7));
    Writeln('0 = ',getGetal(O+7,O+8));
  end;
  zwart := getGetal(O+0,O+4);
  wit := getGetal(O+4,O+8);
  OK := OK and (zwart = 0) and (wit = $FFFFFF);
  if not OK then
    Writeln('Not a proper B/W picture !!')
  else
    if debug then Writeln('B/W picture OK');
  LaatZien := OK; 
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;
  { Pixels beyond Screen edges assumed Black }
    zwart := ((x = 0)    and (ix[k] = 0))
          or ((x = Wijd) and (ix[k] = 1))
          or ((y = 0)    and (iy[k] = 0))
          or ((y = Hoog) and (iy[k] = 1));
    if zwart then getal := getal + 1;
    if zwart then Continue;
    zwart := GetPixel(x+ix[k]-1,y+iy[k]-1);
    if zwart then getal := getal + 1;
  end;
  if zwart then donker := donker + 1;
  MetVieren := getal;
end;

procedure 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 }
  for k := 1 to Geteld do
    bewaar[k] := bewaar[2*k]; 
end;

procedure 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-1;
  bewaar[tel-1].y := y1-1;
{ Even counts }
  bewaar[tel].x := x2-1;
  bewaar[tel].y := y2-1;
end;

begin
  TijdMeting('');
  tel := 0;
  donker := 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';
  TijdMeting('Contouring pass I' + extra);
end;

procedure 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;
{
  Just a Test
}
var
  k : integer;
  ff : TextFile;
  w : schakel;
begin
  AssignFile(ff,'bewaar.txt');
  Rewrite(ff);
  for k := 1 to Geteld do
  begin
    w := bewaar[k];
    Writeln(ff,w.x:6,w.y:6,w.lnk:6);
  end;
  CloseFile(ff);
end;

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

procedure BMProutines;
{
  All Together
}
const
  debug : boolean = false;
var
  alles,genoeg,k : integer;
begin
  Contouren(false);
  donker := donker - (Wijd + Hoog - 1);
  alles := Hoog*Wijd;
  if debug then
  begin
    Writeln('Contour Pieces = ' + IntToStr(Geteld));
    Writeln('Black Pixels # = ' + IntToStr(donker));
    Writeln('Black/White = '+ IntToStr(Round(100*donker/alles))+' %');
  end;

  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;

procedure BMPschoon;
{
  Clear structures
}
begin
  SetLength(bewaar,0);
  SetLength(Heel,0);
  SetLength(randen,0);
end;

function Kromme;
{
  Define Contour
}
var
  k,L,tel : integer;
  rond : Jordan;
begin
  tijdelijk := deze;
  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 Tekenen;
{
  Do the Drawing
  of a Contour
}
var
  k,L, x,y : integer;
begin
  k := randen[welk] ;
  x := bewaar[k].x div 2;
  y := bewaar[k].y div 2;
  Scherm.Canvas.Pen.Color := C;
  Scherm.Canvas.MoveTo(x,y);
  L := k;
  while true do
  begin
    L := bewaar[L].lnk;
    if L = k then Break;
    x := bewaar[L].x div 2;
    y := bewaar[L].y div 2;
    Scherm.Canvas.LineTo(x,y);
  end;
end;

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

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

function StartPunt;
{
  Start of Contour
}
var
  ingang : integer;
  p : vektor;
begin
  ingang := randen[deze];
  p.x := bewaar[ingang].x/2;
  p.y := bewaar[ingang].y/2;
  StartPunt := p;
end;

function Huidige;
{
  Key of Contour
}
begin
  Huidige := tijdelijk;
end;

procedure FanDataCompressie;
{
  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 (1/3 pixel size)
}
var
  L,k,m,kept,tel, memo : integer;
  pu,pv, dist, h,dh : double;
  keep,verder : boolean;
  tmp,o,p : vektor;
const
  eps : double = 1/3;
  u1 : double = 0 ; u2 : double = 0 ;
  f1 : double = 0 ; f2 : double = 0 ;
  f3 : double = 0 ;
begin
  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;

  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.