unit Unit10;
{
Area inside convex Jordan curve
===============================
calculated numerically by grid refinement
=========================================
This software has been designed and it is
CopyLefted by Han de Bruijn:

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

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, ExtCtrls, Grafisch;

type
  TForm1 = class(TForm)
    Image1: TImage;
    procedure Toetsdruk(Sender: TObject; var Key: Char);
    procedure Scheppen(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form1: TForm1;

IMPLEMENTATION

{$R *.dfm}

type
  punt = record
    x,y : double;
  end;
  punten = array[0..3] of punt;
const
  achille : string = ' 3.2040508691540993294644483966405124836692674716427904';
var
  basis : integer; { Refinement factor }
  diepte : integer; { Depth of refinement }
{ Inner / Outer area / Perimeter: }
  Opp,dOpp,Uit,dUit,Lengte : double;
  teller : array of integer; { # of squares }
  kijken : boolean;

procedure Read_Parameters(var basis,diepte : integer);
{
  Syntax: [program] [number base] [refinement depth]
}
const
  maximum : array[2..3] of integer = (21,13);
var
  OK : boolean;
  tel : integer;
  woord : string;
begin
  OK := true;
  tel := ParamCount;
  if not (tel < 3) then OK := false;
  if not OK then Halt;
  basis := 3; diepte := maximum[3];
  if tel = 0 then Exit;

  woord := ParamStr(1);
  if (Length(woord) > 1)
  or not (woord[1] in ['2','3']) then
    Writeln(woord,' out of [2..3] range; default = 3')
  else basis := StrToInt(woord);
  if basis = 2 then diepte := maximum[2];
  if tel = 1 then Exit;

  woord := ParamStr(2);
  OK := not (Length(woord) > 2);
  OK := OK and (woord[1] in ['0'..'9']);
  if Length(woord) = 2 then
    OK := OK and (woord[2] in ['0'..'9']);
  if OK then
    diepte := StrToInt(woord)
  else
    Writeln(woord,' out of range; default chosen');
  if diepte > maximum[basis] then
  begin
    diepte := maximum[basis];
    Writeln(woord,' out of range; default chosen');
  end;
end;

procedure coordinaat(diep,getal : integer; var c0,c1 : double);
{
  Calculate normed coordinates (c) of squares
  with depth (diep) and given coding (getal)
}
var
  k,macht,g,m : integer;
  c : double;
  rij : array of integer;
begin
  SetLength(rij,diep);
  g := abs(getal);
  for k := 1 to diep do
  begin
    m := g mod basis;
    rij[k-1] := m;
    g := g div basis;
  end;
  macht := 1;
  c := 0;
  for k := diep downto 1 do
  begin
    macht := macht*basis;
    m := rij[k-1];
    c := c + m/macht;
  end;
  c0 := c; c1 := c0 + 1/macht;
end;

function F(x,y : double) : double;
{
  Sample Jordan curve
  dividing the plane
}
begin
  F := sqrt(sqr(x-2)+sqr(y-3)) + 2*sqrt(sqr(x-3)+sqr(y-1)) - 4;
end;

procedure extra(p : punten; var dA,dU,dL : double);
{
  Additional Isolines Approach
}
const
  hoe : array[0..3] of array[0..2] of integer =
        ((0,1,2),(1,0,3),(2,3,0),(3,2,1));
var
  g : array[0..3] of double;
  hoek : array[0..7] of boolean;
  boven,onder,links,rechts : boolean;
  x0,x1,y0,y1 : double;
  k,m,i,j : integer;
  vierkant : double;
begin
  for k := 0 to 3 do
    g[k] := F(p[k].x,p[k].y);

  vierkant := (p[1].x-p[0].x)*(p[2].y-p[0].y);

  hoek[0] := (g[0] > 0)  and (g[1] <= 0) and (g[2] <= 0) and (g[3] <= 0);
  hoek[1] := (g[0] <= 0) and (g[1] > 0)  and (g[2] <= 0) and (g[3] <= 0);
  hoek[2] := (g[0] <= 0) and (g[1] <= 0) and (g[2] > 0)  and (g[3] <= 0);
  hoek[3] := (g[0] <= 0) and (g[1] <= 0) and (g[2] <= 0) and (g[3] > 0) ;

  hoek[4] := (g[0] <= 0) and (g[1] > 0)  and (g[2] > 0)  and (g[3] > 0);
  hoek[5] := (g[0] > 0)  and (g[1] <= 0) and (g[2] > 0)  and (g[3] > 0);
  hoek[6] := (g[0] > 0)  and (g[1] > 0)  and (g[2] <= 0) and (g[3] > 0);
  hoek[7] := (g[0] > 0)  and (g[1] > 0)  and (g[2] > 0)  and (g[3] <= 0);

  for m := 0 to 7 do
  begin
    if not hoek[m] then Continue;
    i := hoe[m mod 4][0]; j := hoe[m mod 4][1];
    k := hoe[m mod 4][2];
    x0 := (g[j]*p[i].x-g[i]*p[j].x)/(g[j]-g[i]);
    y0 := (g[k]*p[i].y-g[i]*p[k].y)/(g[k]-g[i]);
    if kijken then
    begin
      Form1.Image1.Canvas.MoveTo(x2i(x0),y2j(p[i].y));
      Form1.Image1.Canvas.LineTo(x2i(p[i].x),y2j(y0));
    end;
    if m < 4 then
    begin
      dU := abs((x0-p[i].x)*(y0-p[i].y))/2;
      dA := vierkant - dU;
    end else begin
      dA := abs((x0-p[i].x)*(y0-p[i].y))/2;
      dU := vierkant - dA;
    end;
    dL := sqrt(sqr(x0-p[i].x)+sqr(y0-p[i].y));
  end;

  boven := (g[0] > 0)  and (g[1] > 0)  and (g[2] <= 0) and (g[3] <= 0);
  onder := (g[0] <= 0) and (g[1] <= 0) and (g[2] > 0)  and (g[3] > 0);
  if boven or onder then
  begin
    y0 := (g[2]*p[0].y-g[0]*p[2].y)/(g[2]-g[0]);
    y1 := (g[3]*p[1].y-g[1]*p[3].y)/(g[3]-g[1]);
    if kijken then
    begin
      Form1.Image1.Canvas.MoveTo(x2i(p[0].x),y2j(y0));
      Form1.Image1.Canvas.LineTo(x2i(p[1].x),y2j(y1));
    end;
    dL := sqrt(sqr(p[1].x-p[0].x)+sqr(y1-y0));
    if boven then
    begin
      dU := (p[1].x-p[0].x)*((y1+y0)/2-p[0].y);
      dA := vierkant - dU;
    end else begin
      dA := (p[1].x-p[0].x)*((y1+y0)/2-p[0].y);
      dU := vierkant - dA;
    end;
  end;

  rechts := (g[0] > 0)  and (g[1] <= 0) and (g[2] > 0)  and (g[3] <= 0);
  links  := (g[0] <= 0) and (g[1] > 0)  and (g[2] <= 0) and (g[3] > 0);
  if links or rechts then
  begin
    x0 := (g[1]*p[0].x-g[0]*p[1].x)/(g[1]-g[0]);
    x1 := (g[3]*p[2].x-g[2]*p[3].x)/(g[3]-g[2]);
    if kijken then
    begin
      Form1.Image1.Canvas.MoveTo(x2i(x0),y2j(p[0].y));
      Form1.Image1.Canvas.LineTo(x2i(x1),y2j(p[2].y));
    end;
    dL := sqrt(sqr(x1-x0)+sqr(p[2].y-p[0].y));
    if rechts then
    begin
      dU := (p[2].y-p[0].y)*((x1+x0)/2-p[0].x);
      dA := vierkant - dU;
    end else begin
      dA := (p[2].y-p[0].y)*((x1+x0)/2-p[0].x);
      dU := vierkant - dA;
    end;
  end;
end;

procedure brutaal(d,i,j : integer);
{
  Recursive Refinement
}
var
  p : punten;
  binnen,buiten : boolean;
  k,m : integer;
  c0,c1,x0,y0,x1,y1 : double;
  dA,dU,dL : double;
begin
  if d > diepte then Exit; 
  teller[d] := teller[d] + 1;

{ Calculate squares coordinates}
  coordinaat(d,i,c0,c1);
  x0 := xmin + (xmax-xmin)*c0;
  x1 := xmin + (xmax-xmin)*c1;
  coordinaat(d,j,c0,c1);
  y0 := ymin + (ymax-ymin)*c0;
  y1 := ymin + (ymax-ymin)*c1;
  p[0].x := x0; p[0].y := y0;
  p[1].x := x1; p[1].y := y0;
  p[2].x := x0; p[2].y := y1;
  p[3].x := x1; p[3].y := y1;

{ Area inside }
  binnen := true;
  for k := 0 to 3 do
    binnen := binnen and (F(p[k].x,p[k].y) <= 0);
  if binnen then
  begin
    Opp := Opp + (x1-x0)*(y1-y0);
    teller[d] := teller[d] - 1;
    Exit;
  end;
{ Area outside }
  buiten := true;
  for k := 0 to 3 do
    buiten := buiten and (F(p[k].x,p[k].y) > 0);
  if buiten then
  begin
    Uit := Uit + (x1-x0)*(y1-y0);
    teller[d] := teller[d] - 1;
    Exit;
  end;

  if d = diepte then
  begin
    if kijken then
    with Form1.Image1.Canvas do
    begin
      Pen.Color := clBlack;
      Pen.Width := 1;
      MoveTo(x2i(p[0].x),y2j(p[0].y));
      LineTo(x2i(p[1].x),y2j(p[1].y));
      LineTo(x2i(p[3].x),y2j(p[3].y));
      LineTo(x2i(p[2].x),y2j(p[2].y));
      LineTo(x2i(p[0].x),y2j(p[0].y));
      Pen.Color := clRed;
      Pen.Width := 2;
    end;
    extra(p,dA,dU,dL);
    dOpp := dOpp + dA;
    dUit := dUit + dU;
    Lengte := Lengte + dL;
  end;

{ Recursion }
  for k := 0 to basis-1 do
  begin
    for m := 0 to basis-1 do
    begin
      brutaal(d+1,basis*i+k,basis*j+m);
    end;
  end;
end;

function macht(n : integer) : integer;
{
  Used for (basis^n)
}
var
  p,k : integer;
begin
  p := 1;
  for k := 1 to n do
    p := p*basis;
  macht := p;
end;

procedure just_doit;
{
  Just do it
}
var
  k,m : integer;
  Rand,dx,dy : double;
begin
  Opp := 0; Uit := 0;
  dOpp := 0; dUit := 0;
  SetLength(teller,diepte+1);
  for k := 0 to diepte do
    teller[k] := 0;
  for k := 0 to basis-1 do
  begin
    for m := 0 to basis-1 do
    begin
      brutaal(1,k,m);
    end;
  end;
 
  for k := 0 to diepte do
    Writeln('# squares at depth = ',k,' : ',teller[k]);
  dx := (xmax-xmin)/macht(diepte);
  dy := (ymax-ymin)/macht(diepte);
  Rand := teller[diepte]*dx*dy;
  Writeln('Boundary + Inner Area + Outer Area = Area of Image, in world coordinates:');
  Writeln(Rand,' +',Opp,' +',Uit,' =');
  Writeln(Rand + Opp + Uit,' =',(xmax-xmin)*(ymax-ymin));
  Writeln('Final result with squares:');
  Uit := (xmax-xmin)*(ymax-ymin)-Uit;
  Writeln(Opp:9:6,' < Area <',Uit:9:6);
  Writeln('Isoline refinement:');
  Writeln('  ',Opp + dOpp);
  Writeln(' =',achille);
  Writeln('Perimeter:');
  Writeln(Lengte);
end;

procedure Blauwdruk;
{
  Picture to start with
}
var
  i,j : integer;
  x,y : double;
begin
  for j := 0 to Hoog-1 do
  begin
    y := j2y(j);
    for i := 0 to Wijd-1 do
    begin
      x := i2x(i);
      with Form1.Image1.Canvas do
      begin
        if F(x,y) < 0 then Pixels[i,j] := clYellow
        else Pixels[i,j] := clWhite;
      end;
    end;
  end;
end;

procedure TForm1.Toetsdruk(Sender: TObject; var Key: Char);
{
  Inactive
}
begin
  Exit;
end;

procedure TForm1.Scheppen(Sender: TObject);
{
  At moment of creation
}
const
  naam : array[0..5] of string = ('nul','een','twee','drie','vier','vijf');
begin
  Read_Parameters(basis,diepte);
  xmin := 1.5; xmax := 4.0;
  ymin := 0.2; ymax := 2.8;
  TV(Form1.Image1);
  ClearDevice;
  Blauwdruk;
{ Pictures at given depth }
  kijken := (diepte < Round(ln(Wijd)/ln(basis)));
  just_doit;
{ Pictured if possible }
  if kijken then
    Form1.Image1.Picture.SaveToFile(naam[basis]+IntToStr(diepte)+'.bmp');
end;

END.