unit Unit0;
{
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;
var
  basis : integer; { Refinement factor }
  diepte : integer; { Depth of refinement }
  Opp,Uit : double; { Inner / Outer area }
  teller : array of integer; { # of squares }

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 brutaal(d,i,j : integer);
{
  Recursive Refinement
}
var
  p : array[0..3] of punt;
  binnen,buiten : boolean;
  k,m : integer;
  c0,c1,x0,y0,x1,y1 : 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;

{ Picture of squares at given depth }
  if diepte < Round(ln(Wijd)/ln(basis)) then
{ Refinement pictured }
  if d = diepte then
  begin
    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));
    end;
  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;
  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:');
  Uit := (xmax-xmin)*(ymax-ymin)-Uit;
  Writeln(Opp:9:6,' < Area <',Uit:9:6);
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;
  just_doit;
{ Pictured if possible }
  if diepte < Round(ln(Wijd)/ln(basis)) then
    Form1.Image1.Picture.SaveToFile(naam[basis]+IntToStr(diepte)+'.bmp');
end;

END.