unit DirectionsUnit;

interface

uses TypesUnit, GaussUnit, StatUnit;

type
TDirections = class
  private
    Dir: Integer;
    Dis: Double;
    Sigma: Double;
    h: Double;
    Gauss: TGauss1D;
    pDisInt: Integer;
    pCosA: TFloatArray;
    pSinA: TFloatArray;
    ArrayMask: TIntField;
    LengthMask: TFloatMatrix;
    IntegralMask: TFloatMatrix;
    IntegralNorm: Double;
    procedure CountDirection(i0, i1: Integer);
    procedure MirrorDirection4(i0, i1: Integer);
    procedure MirrorDirection2(i0, i1: Integer);
    procedure MirrorDirection1(i0, i1: Integer);
  public
    constructor Create(Dr: Integer; Ds, Sigm: Double);
    property DisInt: Integer read pDisInt;
    property CosA: TFloatArray read pCosA;
    property SinA: TFloatArray read pSinA;
    function Filtering(ImBr: TByteMatrix): TByteMatrix;
    function CountWeights(ImBr: TByteMatrix; X, Y: Integer): TFloatArray;
    destructor Destroy; override;
  end;

var
  Directions: TDirections;

implementation

constructor TDirections.Create(Dr: Integer; Ds, Sigm: Double);
var
  D8, D4, D2, n: Integer;
begin
  inherited Create;

  Dir:=Dr;
  Dis:=Ds;
  Sigma:=Sigm;
  h:=2*Pi/Dir;
  Gauss:=TGauss1D.Create(Dis, Sigma);
  pDisInt:=RoundUp(Dis-0.5);

  D8:=Dir div 8;
  D4:=Dir div 4;
  D2:=Dir div 2;

  n:=2*DisInt+1;
  SetLength(pCosA, Dir);
  SetLength(pSinA, Dir);
  SetLength(ArrayMask[0], Dir, n);
  SetLength(ArrayMask[1], Dir, n);
  SetLength(LengthMask, Dir, n);
  SetLength(IntegralMask, Dir, n);

  CountDirection(0, D8);
  if Dir mod 4 = 0 then
    MirrorDirection4(D8+1, D4)
  else
    CountDirection(D8+1, D4);
  if Dir mod 2 = 0 then
    MirrorDirection2(D4+1, D2)
  else
    CountDirection(D4+1, D2);
  MirrorDirection1(D2+1, Dir-1);

  IntegralNorm:=Gauss.IntegralNorm;
end;



procedure TDirections.CountDirection(i0, i1: Integer);
var i, j, k, l, m: Integer;
    A, AbsCosA, AbsSinA, HCross, VCross, NextHCross, Leng, Bound, Sum, LinInt: Double;
begin
  A:=2*pi-i0*h;
  for i:=i0 to i1 do begin
    pCosA[i]:=Cos(A);
    pSinA[i]:=Sin(A);
    AbsCosA:=Abs(CosA[i]);
    AbsSinA:=Abs(SinA[i]);
    j:=0;
    k:=0;
    l:=0;
    m:=0;
    Leng:=0;
    Bound:=0;
    Sum:=0;
    ArrayMask[0,i,0]:=DisInt;
    ArrayMask[1,i,0]:=DisInt;

    if SinA[i] <> 0 then
      HCross:=0.5/AbsSinA
    else
      HCross:=cInfty;
    NextHCross:=HCross;

    if CosA[i] <> 0 then
      VCross:=0.5/AbsCosA
    else
      VCross:=cInfty;

    while (VCross < Dis) or (HCross < Dis) do begin
      if HCross <= VCross then begin
        Bound:=HCross;
        if SinA[i] > 0 then
          k:=k-1
        else
          k:=k+1;
        NextHCross:=(Abs(k)+0.5)/AbsSinA;
      end;

      if VCross <= HCross then begin
        Bound:=VCross;
        if CosA[i] > 0 then
          l:=l+1
        else
          l:=l-1;
        VCross:=(Abs(l)+0.5)/AbsCosA;
      end;
      HCross:=NextHCross;

      LengthMask[i,j]:=Bound-Leng;
      Leng:=Bound;
      while Gauss.x[m] < Bound do
        m:=m+1;
      LinInt:=(Gauss.x[m]-Bound)/(Gauss.x[m]-Gauss.x[m-1]);
      IntegralMask[i,j]:=Gauss.IntegralSum[m-1]*LinInt+Gauss.IntegralSum[m]*(1-LinInt)-Sum;
      Sum:=Sum+IntegralMask[i,j];

      j:=j+1;
      ArrayMask[0,i,j]:=DisInt+k;
      ArrayMask[1,i,j]:=DisInt+l;
    end;
    LengthMask[i,j]:=Dis-Leng;
    IntegralMask[i,j]:=Gauss.IntegralNorm-Sum;
    SetLength(ArrayMask[0,i], j+1);
    SetLength(ArrayMask[1,i], j+1);
    SetLength(LengthMask[i], j+1);
    SetLength(IntegralMask[i], j+1);
    A:=A-h;
  end;
end;



procedure TDirections.MirrorDirection4(i0, i1: Integer);
var i, i1mi: Integer;
begin
  for i:=0 to i1-i0 do begin
    i1mi:=i1-i;
    pCosA[i1mi]:=-SinA[i];
    pSinA[i1mi]:=-CosA[i];
    ArrayMask[0,i1mi]:=ArrayMask[1,i];
    ArrayMask[1,i1mi]:=ArrayMask[0,i];
    LengthMask[i1mi]:=LengthMask[i];
    IntegralMask[i1mi]:=IntegralMask[i];
  end;
end;



procedure TDirections.MirrorDirection2(i0, i1: Integer);
var i, j, k, i1mi: Integer;
begin
  for i:=0 to i1-i0 do begin
    i1mi:=i1-i;
    pSinA[i1mi]:=SinA[i];
    pCosA[i1mi]:=-CosA[i];
    ArrayMask[0,i1mi]:=ArrayMask[0,i];
    k:=High(ArrayMask[1,i]);
    for j:=0 to k do
      ArrayMask[1,i1mi,j]:=2*DisInt-ArrayMask[1,i,j];
    SetLength(ArrayMask[1,i1mi], k+1);
    LengthMask[i1mi]:=LengthMask[i];
    IntegralMask[i1mi]:=IntegralMask[i];
  end;
end;



procedure TDirections.MirrorDirection1(i0, i1: Integer);
var i, j, k, i1mi, ip1: Integer;
begin
  for i:=0 to i1-i0 do begin
    i1mi:=i1-i;
    ip1:=i+1;
    pSinA[i1mi]:=-SinA[ip1];
    pCosA[i1mi]:=CosA[ip1];
    k:=High(ArrayMask[0,ip1]);
    for j:=0 to k do
      ArrayMask[0,i1mi,j]:=2*DisInt-ArrayMask[0,ip1,j];
    SetLength(ArrayMask[0,i1mi], k+1);
    ArrayMask[1,i1mi]:=ArrayMask[1,ip1];
    LengthMask[i1mi]:=LengthMask[ip1];
    IntegralMask[i1mi]:=IntegralMask[ip1];
  end;
end;



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

  SetLength(Result, ImH, ImW);
  SetLength(Direction, 2*DisInt+1);

  for i:=0 to ImH-1 do
    for j:=0 to ImW-1 do begin
      NewPix:=0;
      DeSum:=0;
      for k:=0 to High(IntegralMask) do begin
        NewPixDir:=0;
        Hig:=High(IntegralMask[k]);
        SetLength(Direction, Hig+1);
        for l:=0 to Hig do begin
          Direction[l]:=ImBr[i+ArrayMask[0,k,l],j+ArrayMask[1,k,l]];
          NewPixDir:=NewPixDir+Direction[l]*IntegralMask[k,l];
        end;
        De:=1/(Variance(Direction, LengthMask[k], Dis)+cZeroVariance);
        NewPix:=NewPix+NewPixDir*De;
        DeSum:=DeSum+De;
      end;
      NewPix:=NewPix/DeSum/IntegralNorm;
      if NewPix > 255 then
        Result[i,j]:=255
      else
        Result[i,j]:=Round(NewPix);
    end;

  Direction:=nil;
end;



function TDirections.CountWeights(ImBr: TByteMatrix; X, Y: Integer): TFloatArray;
var i, j, nx, ny, Hig: Integer;
    DeMax: Double;
    Direction: TByteArray;
begin
  nx:=X-DisInt;
  ny:=Y-DisInt;

  DeMax:=0;
  SetLength(Result, Dir+1);

  for i:=0 to Dir-1 do begin
    Hig:=High(ArrayMask[0,i]);
    SetLength(Direction, Hig+1);
    for j:=0 to Hig do begin
      Direction[j]:=ImBr[ny+ArrayMask[0,i,j],nx+ArrayMask[1,i,j]];
    end;
    Result[i]:=1/(Variance(Direction, LengthMask[i], Dis)+cZeroVariance);
    if Result[i] > DeMax then
      DeMax:=Result[i];
  end;

  Result[Dir]:=DeMax;
  Direction:=nil;
end;



destructor TDirections.Destroy;
begin
  try
    Gauss.Destroy;
    ArrayMask[0]:=nil;
    ArrayMask[1]:=nil;
    LengthMask:=nil;
    IntegralMask:=nil;
  finally
    inherited;
  end;
end;

end.
