unit variogram;


interface

uses
  Windows, Messages, SysUtils, Math, Classes, Graphics, Controls, Forms, Dialogs,
  StdCtrls, Gauges;

type
  TBmpForm = class(TForm)
    Button1: TButton;
    Gauge1: TGauge;
    Gauge2: TGauge;
    OpenDialog1: TOpenDialog;
    Label1: TLabel;
    procedure FormDestroy(Sender: TObject);
    procedure FormPaint(Sender: TObject);
    procedure Button1Click(Sender: TObject);
    procedure FormCreate(Sender: TObject);

    private
    { Private declarations }
    Bitmap: TBitmap;
    procedure VarioBitmap;

    procedure WMEraseBkgnd(var m: TWMEraseBkgnd); message WM_ERASEBKGND;

  public
    { Public declarations }
  end;

var
  BmpForm: TBmpForm;
  Fout: Text;
  Outfile: String[80];

implementation


{$R *.DFM}


procedure TBmpForm.FormCreate(Sender: TObject);
begin
  Bitmap := TBitmap.Create;
  OpenDialog1.Execute;
  OpenDialog1.Filename;
  Bitmap.LoadFromFile( OpenDialog1.Filename );
  VertScrollBar.Range:= Bitmap.Height;
  HorzScrollBar.Range:= Bitmap.Width;

end;

procedure TBmpForm.FormDestroy(Sender: TObject);
begin
  Bitmap.Free;
end;


procedure TBmpform.WMEraseBkgnd(var m : TWMEraseBkgnd);
begin
  m.Result := LRESULT(False);
end;


procedure TBmpForm.FormPaint(Sender: TObject);
begin

      Canvas.Draw(-HorzScrollBar.Position, -VertScrollBar.Position, Bitmap);

end;

procedure TBmpForm.Button1Click(Sender: TObject);
begin
  if Button1.Caption <> 'Finished' then
  begin
    Button1.Caption:= 'Computing Variogram';
    VarioBitmap;
    Button1.Caption:= 'Finished';
    Invalidate;
  end;
end;



procedure TBmpForm.VarioBitmap;
   const
      nAngDiv = 9;
      AllInfo = TRUE;
      max_h = 367;
      Filter = 0;

   var
      n00, n90, h: Integer;
      x, y, i, j: Integer;
      sum1, sum2, semivar, var0: Extended;
      P0 : PByteArray;
      ang : Extended;
      
      Sumx2, Sumy2, Sumxy, Weight: array[1..50] of Extended;
      {
      rho, Angle, theta, Compact, Compaction: Extended;
      hBreak, hNeg: Integer;
      hHorz, hVert: Extended;
      var_x, var_y, cov_xy, Det : Extended;
      A_ellipse, B_ellipse, C_ellipse, Major, Minor, Rot_angle: Extended;
      }

      w, hTrue: Extended;

Function Atan2(top, bot :Extended ):Extended;
Const tol =1.0e-17;
Var Arctan2 : Extended;
begin
   Arctan2:= 0.0;
  if top > 0.0 then
               begin
                 if bot > tol then Arctan2:=arctan(top/bot)
                 else if bot <-tol then Arctan2 := arctan(top/bot) + pi
                 else Arctan2:=pi/2;
               end
  else if top < 0.0 then
               begin
                 if bot < -tol then Arctan2:=arctan(top/bot) + pi
                 else if bot > tol then Arctan2:=arctan(top/bot) + 2.0*pi
                 else Arctan2:=1.5*pi
               end
  else if top = 0.0 then
               begin
                 if bot >= -tol then Arctan2:=0.0
                 else Arctan2:=pi;
               end;
  if Arctan2 >= 2*pi then Arctan2:=Arctan2-2*pi;
  if Arctan2 < 0.0   then Arctan2:=2*pi-Arctan2;
  Atan2:= Arctan2;
end;

Procedure Canonical_Form( a, b, c, f: Extended;
                          Var Major, Minor, Rot_angle: Extended );
{ a*x^2 + 2*b*xy + c*y^2 = f -->  AA*X^2 + BB*Y^2 = f }
Var AA, BB: Extended;
begin
  Rot_angle:= ATAN2( 2*b, (a-c) )/2.0;
  AA:= a*SQR(COS(Rot_angle)) + b*SIN(2*Rot_angle) + c*SQR(SIN(Rot_angle));
  BB:= a*SQR(SIN(Rot_angle)) - b*SIN(2*Rot_angle) + c*SQR(COS(Rot_angle));
  if (AA < BB ) then begin
                       Major:= SQRT(1.0/AA);
                       Minor:= SQRT(1.0/BB);
                     end
                else begin
                       Major:= SQRT(1.0/BB);
                       Minor:= SQRT(1.0/AA);
                       Rot_angle:= Rot_angle + pi/2.0;
                     end;
end;

      procedure CalcVario( angle: Extended; h: Integer; var semiVar: Extended; var wgt: Extended; var hTrue: Extended );
      var sum, sum2: Extended;
          P0, P1 : PByteArray;
          y, x, ystart, yend: Integer;
          nRows, nCols: Integer;

          hVert, hHorz: Integer;

          varSemivar, PDiff: Extended;
      begin

       hHorz:=  Round( h*Cos(angle*Pi/180) );
       hVert:=  Round( h*Sin(angle*Pi/180) );

       hTrue:= Sqrt( hHorz*hHorz+hVert*hVert );
       sum:= 0.0; sum2:= 0.0;

       nRows:= BitMap.height;
       nCols:= BitMap.width;

       if hVert < 0 then
          begin
             ystart:= -hVert;
             yend:=   nRows-1;
          end
          else
          begin
             ystart:= 0;
             yend:=   nRows-1-hVert;
          end;

        for y := ystart to yend do
        begin
          P0:= BitMap.ScanLine[y];
          P1:= BitMap.ScanLine[y+hVert];


          for x := 0 to nCols-1-hHorz  do
          begin

              {*****}
               if P0[x] < Filter then P0[x]:= 255;
               if P1[x+hHorz] < Filter then P1[x+hHorz]:= 255;
              {*****}

            PDiff:= P0[x]-P1[x+hHorz];
            sum:= sum + SQR( PDiff );
            sum2:= sum2 + SQR( PDiff )*SQR( PDiff );
          end;
        end;
        semiVar:= sum/(2*(nCols - hHorz)*(nRows-Abs(hVert)));
        varSemivar:= sum2/(4*(nCols - hHorz)*(nRows-Abs(hVert)))- SQR(semiVar);
        wgt:= varSemivar/((nCols - hHorz)*(nRows-Abs(hVert)));
      end;

   begin
      OutFile:= OpenDialog1.Filename;
      OutFile[length(OutFile)]:= 'v';
      OutFile[length(OutFile)-1]:= 's';
      OutFile[length(OutFile)-2]:= 'c';

      AssignFile(Fout,OutFile);
      Rewrite(Fout);

      n00:= Bitmap.Width;
      n90:= Bitmap.Height;


      {Variance Calculation}
        sum1:= 0.0;  sum2:= 0.0;
        for y := 0 to n90-1 do
        begin
          P0:= BitMap.ScanLine[y];

          for x := 0 to n00-1  do
          begin

          {*****}
          if P0[x] < Filter then P0[x]:= 255;
          {*****}

          sum1:= sum1 + P0[x];
          sum2:= sum2 + P0[x]*P0[x];
          end;
        end;
        var0:= sum2/(n00*n90)- Sqr( sum1/(n00*n90) );

      Writeln(Fout, 'Dataset: ,',Outfile );
      Writeln(Fout, 'Filter < :,', Filter:10 );
      Writeln(Fout, 'Image Size:,', n00:15, ',', n90:15 );
      Writeln(Fout, 'Variance:,', var0:15:5 );

      for j:= 1 to max_h Div 10 do
      begin
        Sumx2[j]:= 0.0; Sumy2[j]:= 0.0; Sumxy[j]:= 0.0; Weight[j]:= 0.0;
      end;
      {
      hBreak:= max_h;
      hNeg:= max_h;
      }
   For i := -nAngDiv to nAngDiv do
   BEGIN
      ang:= i*90/nAngDiv;
      if AllInfo then
      begin
        Writeln(Fout);
        Writeln( Fout, 'Angle: ,', -ang:6:1  );
        Writeln(Fout, 'h, SemiVar Angle, 1/weight');
        Writeln(Fout,'0,0.0,1');
      end;

      for h := 1 to max_h do
      begin
        calcVario(  ang, h, semivar, w, hTrue );
        {
        rho:= 1.0 - semivar/var0;
        }
        if AllInfo then
        writeln(Fout,  h:10, ',', semivar:15:5, ',', w:15:5 );
        {
        if (rho < 0.0)  and (h < hNeg)   then hNeg:= h;
        if (rho < 0.01) and (h < hBreak) then hBreak:= h;

        rho:= Abs( rho );
        hHorz:=  h*Cos(ang*Pi/180);
        hVert:=  h*Sin(ang*Pi/180);


        for j := 1 to max_h Div 10 do
        begin
           if h <= 10*j then
           begin
                Sumx2[j]:= Sumx2[j] + h*rho*hHorz*hHorz;
                Sumy2[j]:= Sumy2[j] + h*rho*hVert*hVert;
                Sumxy[j]:= Sumxy[j] + h*rho*hHorz*hVert;
                Weight[j]:= Weight[j] + h*rho;
           end;
        end;
    }

        if (h mod (max_h Div 10)) = 0 then
        begin
          Gauge2.Progress:= (100*h) Div max_h;
        end;
      end;
      Gauge1.Progress:= ((i+nAngDiv)*100) Div (2*nAngDiv);
      Gauge2.Progress:=  0;
   END;

{
   Writeln(Fout, 'h,% Compaction,Angle' );
   For j:= 1 to max_h Div 10 do
   begin
    if 10*j < hNeg then
    begin
      var_x:=  sumx2[j]/Weight[j];
      var_y:=  sumy2[j]/Weight[j];
      cov_xy:= sumxy[j]/Weight[j];
      Det:= var_x*var_y-SQR(cov_xy);

      A_ellipse:= var_y/Det;
      B_ellipse:= -cov_xy/Det;
      C_ellipse:= var_x/Det;

      Canonical_Form( A_ellipse, B_ellipse, C_ellipse, 1.0, Major, Minor, Rot_angle );

      Compaction:= 100.0*(1-Minor/Major);
      Angle:= -Rot_angle*180/Pi;

      If Angle >  90.0  then Angle:= Angle - 180.0;
      If Angle < -90.0 then Angle:= Angle + 180.0;

      Writeln(Fout, (j*10):6, ',', Compaction:7:3, ',', Angle:7:2 );
      if 10*j <= hBreak then
      begin
        compact:= compaction;
        theta:= Angle;
      end;
    end;
   end;

    Writeln(Fout);
    Writeln(Fout, 'hBreak,hNeg');
    Writeln(Fout, hBreak:6, ',', hNeg:6);
    CloseFile(Fout);

    Label1.Caption:= Format( '        %s        %f %% compaction at%5.1f degrees        h Negative: %d',
                     [OpenDialog1.Filename, compact, theta, hNeg] );
    Label1.Visible:= True;
}
    Label1.Caption:= 'Finished';
    Label1.Visible:= True;
    CloseFile(Fout);

   end;

end.