unit LinearUnit;

interface

uses Math, TypesUnit, GaussUnit;

type TLinear = class
  private
    Dis: Double;
    Sigma: Double;
    DisInt: Integer;
    Gauss: TGauss;
    LinearMask: TFloatMatrix;
    function CountMask1_8: TFloatMatrix;
    procedure MirrorFullMask(Mask1_8: TFloatMatrix);
    procedure NormMask;
  public
    constructor Create(Ds, Sigm: Double);
    function Filtering(ImBr: TByteMatrix): TByteMatrix;
    destructor Destroy; override;
  end;

var
  Linear: TLinear;

implementation

constructor TLinear.Create(Ds, Sigm: Double);
begin
  inherited Create;
  Dis:=Ds;
  Sigma:=Sigm;
  DisInt:=RoundUp(Dis-0.5);

  Gauss:=TGauss.Create(Dis, Sigma);
  MirrorFullMask(CountMask1_8);
  NormMask;
end;



function TLinear.CountMask1_8: TFloatMatrix;
var i, j, k, Over: Integer;
    im, jm, ip, jp, Alpha, Beta, Angle, a, b, S, Area: Double;
    v: array[0..3] of Double;
begin
  SetLength(Result, DisInt+1);

  for i:=0 to DisInt do begin
    SetLength(Result[i], DisInt+1-i);
    for j:=i to DisInt do begin
      Result[i,j-i]:=0;
      im:=max(i-0.5,0);
      jm:=max(j-0.5,0);
      v[0]:=Sqrt(im*im+jm*jm);
      if v[0] >= Dis then
        Continue;
      ip:=i+0.5;
      jp:=j+0.5;
      v[1]:=Sqrt(ip*ip+jm*jm);
      v[2]:=Sqrt(im*im+jp*jp);
      v[3]:=Sqrt(ip*ip+jp*jp);

      Area:=0;
      Over:=1;
      k:=-RoundUp(-v[0]/cGaussStep-1)-1;
      repeat
        k:=k+1;
        S:=0;
        while (Gauss.x[k]>v[min(Over,3)]) and (Over<=3) do
          Over:=Over+1;
        case Over of
          1 : begin
                Alpha:=ArcSin(im/Gauss.x[k]);
                Beta:=ArcCos(jm/Gauss.x[k]);
                Angle:=Beta-Alpha;
                a:=Gauss.x[k]*Cos(Alpha)-jm;
                b:=Gauss.x[k]*Sin(Beta)-im;
                S:=(a*b + Gauss.x[k]*Gauss.x[k]*(Angle-Sin(Angle)))/2 - Area;
              end;
          2 : begin
                Alpha:=ArcSin(im/Gauss.x[k]);
                Beta:=ArcSin(ip/Gauss.x[k]);
                Angle:=Beta-Alpha;
                a:=Gauss.x[k]*Cos(Alpha)-jm;
                b:=Gauss.x[k]*Cos(Beta)-jm;
                S:=b*(ip-im) + ((a-b)*(ip-im) + Gauss.x[k]*Gauss.x[k]*(Angle-Sin(Angle)))/2 - Area;
              end;
          3 : begin
                Alpha:=ArcCos(jp/Gauss.x[k]);
                Beta:=ArcSin(ip/Gauss.x[k]);
                Angle:=Beta-Alpha;
                a:=ip-Gauss.x[k]*Sin(Alpha);
                b:=jp-Gauss.x[k]*Cos(Beta);
                S:=(ip-im)*(jp-jm) - (a*b - Gauss.x[k]*Gauss.x[k]*(Angle-Sin(Angle)))/2 - Area;
              end;
          4 : S:=(ip-im)*(jp-jm) - Area;
        end;
        Result[i,j-i]:=Result[i,j-i]+S*Gauss.y[k];
        Area:=Area+S;
      until (Gauss.x[k]>v[3]) or (Gauss.x[k]=Dis);
    end;
  end;
  Result[0,0]:=4*Result[0,0];
  for j:=1 to DisInt do
    Result[0,j]:=2*Result[0,j];
end;



procedure TLinear.MirrorFullMask(Mask1_8: TFloatMatrix);
var i, j, n: Integer;
begin
  n:=2*DisInt+1;
  SetLength(LinearMask,n,n);

  for i:=DisInt to n-1 do begin
    LinearMask[i,i]:=Mask1_8[i-DisInt,0];
    for j:=i+1 to n-1 do begin
      LinearMask[i,j]:=Mask1_8[i-DisInt,j-i];
      LinearMask[j,i]:=Mask1_8[i-DisInt,j-i];
    end;
  end;

  for i:=DisInt to n-1 do
    for j:=0 to DisInt-1 do
      LinearMask[i,j]:=LinearMask[i,n-1-j];

  for i:=0 to DisInt-1 do
    for j:=0 to n-1 do
      LinearMask[i,j]:=LinearMask[n-1-i,j]
end;



procedure TLinear.NormMask;
var i, j, n: Integer;
    Norm: Double;
begin
  n:=2*DisInt;
  Norm:=0;

  for i:=0 to n do
    for j:=0 to n do
      Norm:=Norm+LinearMask[i,j];

  for i:=0 to n do
    for j:=0 to n do
      LinearMask[i,j]:=LinearMask[i,j]/Norm;
end;



function TLinear.Filtering(ImBr: TByteMatrix): TByteMatrix;
var ImH, ImW, i, j, k, l: Integer;
    NewPix: Double;
begin
  ImH:=High(ImBr)-2*DisInt+1;
  ImW:=High(ImBr[0])-2*DisInt+1;

  SetLength(Result, ImH, ImW);

  for i:=0 to ImH-1 do
    for j:=0 to ImW-1 do begin
      NewPix:=0;
      for k:=0 to 2*DisInt do
        for l:=0 to 2*DisInt do
          NewPix:=NewPix+ImBr[i+k,j+l]*LinearMask[k,l];
      Result[i,j]:=Min(Round(NewPix),255);
    end;
end;



destructor TLinear.Destroy;
begin
  try
    Gauss.Destroy;
    LinearMask:=nil;
  finally
    inherited;
  end;
end;

end.
