    (*********************************************************************
     *  G32_WConvolution                                                 *
     *                                                                   *
     *   Adapted for G32 by Vladimir Vasilyev                            *
     *   home page : http://www.nextnow.com                              *
     *   email     : Vladimir@nextnow.com                                *
     *               W-develop@mtu-net.ru                                *
     *                                                                   *
     *   28_02_2003 added procedure DoSobel( Dst, Src: TBitmap32);       *
     *                                                                   *
     *  References:                                                      *
     *  Based on Harm's example of a 3 x 3 convolution                   *
     *  using 24-bit bitmaps and scanline                                *
     *  http://www.users.uswest.net/~sharman1/                           *
     *  sharman1@uswest.net                                              *
     *                                                                   *
     *                                                                   *
     *  The contents of this file are used with permission, subject to   *
     *  the Mozilla Public License Version 1.1 (the "License"); you may  *
     *  not use this file except in compliance with the License. You may *
     *  obtain a copy of the License at                                  *
     *  http://www.mozilla.org/MPL/MPL-1.1.html                          *
     *                                                                   *
     *  Software distributed under the License is distributed on an      *
     *  "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, either express or   *
     *  implied. See the License for the specific language governing     *
     *  rights and limitations under the License.                        *
     *                                                                   *
     *********************************************************************)

unit G32_WConvolution;

interface

Uses GR32;

Const
  HoughTettaK = 6;
  HoughW      = 1500;

  DirX   : Array[0..7] of ShortInt=( 0,1,0,-1,-1, 1,1,-1);
  DirY   : Array[0..7] of ShortInt=(-1,0,1, 0,-1,-1,1, 1);


 ColorArr : array[0..4] of TColor32=( clRed32, clYellow32, clLime32, clBlue32, clFuchsia32);


 GBSD1Kernel:array[0..25]  of Integer=( 1, 4, 7, 4,1,
                                        4,16,26,16,4,
                                        7,26,41,26,7,
                                        4,16,26,16,4,
                                        1, 4, 7, 4,1, 273);

 GBKernel3x3 : array[0..9] of Integer=( 1, 2, 1,
                                        2, 4, 2,
                                        1, 2, 1, 16);

 Lap1Kernel3x3:array[0..9] of Integer=( 0,-1, 0,
                                       -1, 4,-1,
                                        0,-1, 0, 1 );

 Lap2Kernel3x3:array[0..9] of Integer=(-1,-1,-1,
                                       -1, 8,-1,
                                       -1,-1,-1, 1 );

 SobelHx : array[0..9]     of Integer=( 1, 0,-1,
                                        2, 0,-2,
                                        1, 0,-1, 4 );

 SobelHy : array[0..9]     of Integer=( 1, 2, 1,
                                        0, 0, 0,
                                       -1,-2,-1, 4 );

Type

 TPoint2f = record
            x: Single;
            y: Single;
 end;

 TPoint3f = record
            x: Single;
            y: Single;
            z: Single;
 end;

 TPoint4f = record
            x: Single;
            y: Single;
            z: Single;
            w: Single;
 end;

 
 TPoint1fArr = Array of Single;
 TPoint2fArr = Array of TPoint2f;

 T2DPoint1fArr = Array of Array of Single;
 T2DPoint2fArr = Array of Array of TPoint2f;
 T2DPoint3fArr = Array of Array of TPoint3f;
 T2DPoint4fArr = Array of Array of TPoint4f;


 TColorFunc = function(Color32: TColor32): Integer;
 TColorMode = (cmGray,cmRed,cmGreen,cmBlue,cmColor);

 TCompareDir = (cdX,cdY);

 procedure ConvolveI(ray : array of integer; z : word; aBmp : TBitmap32);
 procedure ConvolveI5x5(ray : array of integer; z : word; aBmp : TBitmap32);
 procedure Kuwahara5x5( aBmp : TBitmap32; ColorMode : TColorMode );

 procedure Do1DConvolution(Dst, Src: TBitmap32; Const Kernel: array of Integer );
 procedure Do2DConvolutionXP( Dst, Src: TBitmap32; const Kernel: array of Integer);overload;
 procedure Do2DConvolutionXP( Src: TBitmap32; const Kernel: array of Integer);overload;

 procedure DoSobel( Dst, Src: TBitmap32);
 procedure DoSobelDirect( Direction, Dst, Src: TBitmap32; const Kernelx,Kernely: array of Integer);
 procedure DoSuppression( Dst, Direction, Grad: TBitmap32);
 procedure Color2Gray( Src: TBitmap32);
 procedure Color2GrayEX( Src: TBitmap32);
 procedure Do2DConvolutionGrayEX(Dst, Src: TBitmap32; const Kernel: array of Integer);overload;
 procedure Do2DConvolutionGrayEX( Src: TBitmap32; const Kernel: array of Integer);overload;
 procedure DoHough(Dst, Src: TBitmap32);
 procedure GrayEX2Gray(Src : TBitmap32);
 procedure GetLocalMax( Dst, Src : TBitmap32; Var AMaxArr : TPoint2fArr; Threshold    : integer);
 procedure Hough2Lines( Dst: TBitmap32; AMaxArr : TPoint2fArr );
 procedure GetCrosses(var AMaxArr : TPoint2fArr; Var ACrossPointArr : T2DPoint2fArr );



 procedure SortPoint2fArr( Var Arr : TPoint2fArr; Direction : TCompareDir);
 function  Point2f(x,y : Double) : TPoint2f;
 Procedure Circle(Dst : TBitmap32; X, Y, Radius,Color: TColor32);

 procedure SortPoint2fArrSpecial(Var Arr : TPoint2fArr; Direction : TCompareDir);
 procedure CyclHoughBitmat( Dst, Src : TBitmap32);

 procedure Subtruct(Bmp1, Bmp2: TBitmap32);

 procedure SortTetta(Var Arr : TPoint2fArr); //(r,t)

implementation


Uses GR32_Filters, GR32_Blend, Dialogs, Math ;

function Point2f(x,y : Double) : TPoint2f;
begin
 Result.x:=x;
 Result.y:=y;
end;

function Set255(Clr : integer) : integer;
asm
  MOV  EAX,Clr  // store value in EAX register (32-bit register)
  CMP  EAX,254  // compare it to 254
  JG   @SETHI   // if greater than 254 then go set to 255 (max value)

  CMP  EAX,1    // if less than 255, compare to 1
  JL   @SETLO   // if less than 1 go set to 0 (min value)
  RET           // otherwise it doesn't change, just exit
@SETHI:         // Set value to 255
  MOV  EAX,255  // Move 255 into the EAX register
  RET           // Exit (result value is the EAX register value)
@SETLO:         // Set value to 0
  MOV  EAX,0    // Move 0 into EAX register
end;            // Result is in EAX


function UpSet255(Clr : integer) : integer;
asm
  MOV  EAX,Clr  // store value in EAX register (32-bit register)
  CMP  EAX,254  // compare it to 254
  JG   @SETHI   // if greater than 254 then go set to 255 (max value)
  RET
@SETHI:         // Set value to 255
  MOV  EAX,255  // Move 255 into the EAX register
  RET           // Exit (result value is the EAX register value)
end;            // Result is in EAX


procedure ConvolveI(ray : array of integer; z : word; aBmp : TBitmap32);
var
  O, T, C, B : PColor32Array; // Scanlines
  x, y       : integer;
  tBufr      : TBitmap32; // temp bitmap

  Red,Green,Blue     : Integer;
begin
  tBufr := TBitmap32.Create;
  CheckParams(tBufr,aBmp);
  tBufr.Assign(aBmp);

  for x := 1 to aBmp.Height - 2 do begin  // Walk scanlines
    O := aBmp.ScanLine[x];       // New Target (Original)
    T := tBufr.ScanLine[x-1];    //old x-1  (Top)
    C := tBufr.ScanLine[x];      //old x    (Center)
    B := tBufr.ScanLine[x+1];    //old x+1  (Bottom)
  // Now do the main piece

    for y := 1 to (tBufr.Width - 2) do begin  // Walk pixels

      Red:=Set255(
          (
          (RedComponent(T[y-1])*ray[0]) + (RedComponent(T[y])*ray[1]) + (RedComponent(T[y+1])*ray[2])+
          (RedComponent(C[y-1])*ray[3]) + (RedComponent(C[y])*ray[4]) + (RedComponent(C[y+1])*ray[5])+
          (RedComponent(B[y-1])*ray[6]) + (RedComponent(B[y])*ray[7]) + (RedComponent(B[y+1])*ray[8])
          ) div z );

      Green:=Set255(
          (
          (GreenComponent(T[y-1])*ray[0]) + (GreenComponent(T[y])*ray[1]) + (GreenComponent(T[y+1])*ray[2])+
          (GreenComponent(C[y-1])*ray[3]) + (GreenComponent(C[y])*ray[4]) + (GreenComponent(C[y+1])*ray[5])+
          (GreenComponent(B[y-1])*ray[6]) + (GreenComponent(B[y])*ray[7]) + (GreenComponent(B[y+1])*ray[8])
          ) div z );

      Blue:=Set255(
          (
          (BlueComponent(T[y-1])*ray[0]) + (BlueComponent(T[y])*ray[1]) + (BlueComponent(T[y+1])*ray[2])+
          (BlueComponent(C[y-1])*ray[3]) + (BlueComponent(C[y])*ray[4]) + (BlueComponent(C[y+1])*ray[5])+
          (BlueComponent(B[y-1])*ray[6]) + (BlueComponent(B[y])*ray[7]) + (BlueComponent(B[y+1])*ray[8])
          ) div z );

      O[y]:=Color32(Red, Green, Blue);
    end;
  end;
  tBufr.Free;
end;

procedure ConvolveI5x5(ray : array of integer; z : word; aBmp : TBitmap32);
var
  O, T,T2, C, B, B2 : PColor32Array; // Scanlines
  NS,i              : integer;
  x, y              : integer;
  tBufr             : TBitmap32; // temp bitmap
  Red,Green,Blue    : Integer;

begin
  tBufr := TBitmap32.Create;
  CheckParams(tBufr,aBmp);
  tBufr.Assign(aBmp);

  for x := 2 to aBmp.Height - 3 do begin // Walk scanlines
    O := aBmp.ScanLine[x];     // New Target (Original)
    T2:= tBufr.ScanLine[x-2];  //old x-2  (Top)
    T := tBufr.ScanLine[x-1];  //old x-1  (Top)
    C := tBufr.ScanLine[x];    //old x    (Center)
    B := tBufr.ScanLine[x+1];  //old x+1  (Bottom)
    B2:= tBufr.ScanLine[x+2];  //old x+2  (Bottom)

  // Now do the main piece
    for y := 2 to (tBufr.Width - 3) do begin  // Walk pixels

    //NS:=0;
    //for i:=0 to 4 do

      Red := Set255(
         (
RedComponent(T2[y-2])*ray[0] + RedComponent(T2[y-1])*ray[1] + RedComponent(T2[y])*ray[2] + RedComponent(T2[y+1])*ray[3] + RedComponent(T2[y+2])*ray[4]+
RedComponent( T[y-2])*ray[5] + RedComponent( T[y-1])*ray[6] + RedComponent( T[y])*ray[7] + RedComponent( T[y+1])*ray[8] + RedComponent( T[y+2])*ray[9]+
RedComponent( C[y-2])*ray[10]+ RedComponent( C[y-1])*ray[11]+ RedComponent( C[y])*ray[12]+ RedComponent( C[y+1])*ray[13]+ RedComponent( C[y+2])*ray[14]+
RedComponent( B[y-2])*ray[15]+ RedComponent( B[y-1])*ray[16]+ RedComponent( B[y])*ray[17]+ RedComponent( B[y+1])*ray[18]+ RedComponent( B[y+2])*ray[19]+
RedComponent(B2[y-2])*ray[20]+ RedComponent(B2[y-1])*ray[21]+ RedComponent(B2[y])*ray[22]+ RedComponent(B2[y+1])*ray[23]+ RedComponent(B2[y+2])*ray[24]
          ) div z   );

      Blue := Set255(
         (
BlueComponent(T2[y-2])*ray[0] + BlueComponent(T2[y-1])*ray[1] + BlueComponent(T2[y])*ray[2] + BlueComponent(T2[y+1])*ray[3] + BlueComponent(T2[y+2])*ray[4]+
BlueComponent( T[y-2])*ray[5] + BlueComponent( T[y-1])*ray[6] + BlueComponent( T[y])*ray[7] + BlueComponent( T[y+1])*ray[8] + BlueComponent( T[y+2])*ray[9]+
BlueComponent( C[y-2])*ray[10]+ BlueComponent( C[y-1])*ray[11]+ BlueComponent( C[y])*ray[12]+ BlueComponent( C[y+1])*ray[13]+ BlueComponent( C[y+2])*ray[14]+
BlueComponent( B[y-2])*ray[15]+ BlueComponent( B[y-1])*ray[16]+ BlueComponent( B[y])*ray[17]+ BlueComponent( B[y+1])*ray[18]+ BlueComponent( B[y+2])*ray[19]+
BlueComponent(B2[y-2])*ray[20]+ BlueComponent(B2[y-1])*ray[21]+ BlueComponent(B2[y])*ray[22]+ BlueComponent(B2[y+1])*ray[23]+ BlueComponent(B2[y+2])*ray[24]
          ) div z );

      Green := Set255(
         (
GreenComponent(T2[y-2])*ray[0] + GreenComponent(T2[y-1])*ray[1] + GreenComponent(T2[y])*ray[2] + GreenComponent(T2[y+1])*ray[3] + GreenComponent(T2[y+2])*ray[4]+
GreenComponent( T[y-2])*ray[5] + GreenComponent( T[y-1])*ray[6] + GreenComponent( T[y])*ray[7] + GreenComponent( T[y+1])*ray[8] + GreenComponent( T[y+2])*ray[9]+
GreenComponent( C[y-2])*ray[10]+ GreenComponent( C[y-1])*ray[11]+ GreenComponent( C[y])*ray[12]+ GreenComponent( C[y+1])*ray[13]+ GreenComponent( C[y+2])*ray[14]+
GreenComponent( B[y-2])*ray[15]+ GreenComponent( B[y-1])*ray[16]+ GreenComponent( B[y])*ray[17]+ GreenComponent( B[y+1])*ray[18]+ GreenComponent( B[y+2])*ray[19]+
GreenComponent(B2[y-2])*ray[20]+ GreenComponent(B2[y-1])*ray[21]+ GreenComponent(B2[y])*ray[22]+ GreenComponent(B2[y+1])*ray[23]+ GreenComponent(B2[y+2])*ray[24]
          ) div z    );

      O[y]:=Color32(Red, Green, Blue);

    end;
  end;

  tBufr.Free;


end;


procedure Kuwahara5x5( aBmp : TBitmap32; ColorMode : TColorMode );
Var
  O, T,T2,C,B,B2  : PColor32Array; // Scanlines
  NS,i,j,k,n      : integer;
  x, y            : integer;
  tBufr           : TBitmap32; // temp bitmap
  Red,Green,Blue  : Integer;

  Region1         : array of Double;
  Region2         : array of Double;
  Region3         : array of Double;
  Region4         : array of Double;

  Mean            : array[1..4] of Extended;
  StdDev          : array[1..4] of Extended;

  minStdDev       : Extended;

  ColorFunc       : TColorFunc;
  Color           : array [0..3] of integer;
  step            : integer;
  c1,c2           : integer;
Begin

tBufr := TBitmap32.Create;
Try

  CheckParams(tBufr,aBmp);
  tBufr.Assign(aBmp);


  for x := 2 to aBmp.Height - 3 do
  begin
    O := aBmp.ScanLine[x];     // New Target (Original)

    T2:= tBufr.ScanLine[x-2];  //old x-2  (Top)
    T := tBufr.ScanLine[x-1];  //old x-1  (Top)
    C := tBufr.ScanLine[x];    //old x    (Center)
    B := tBufr.ScanLine[x+1];  //old x+1  (Bottom)
    B2:= tBufr.ScanLine[x+2];  //old x+2  (Bottom)


    //Now slide the region 5x5
    for y := 2 to (tBufr.Width - 3) do
     begin

          //fill regions array
          SetLength(Region1,6);
          SetLength(Region2,6);
          SetLength(Region3,6);
          SetLength(Region4,6);


          Case ColorMode of
             cmColor : begin
                        c1:=1;
                        c2:=3;
                       end;
             else
                        c1:=Ord(ColorMode);
                        c2:=Ord(ColorMode);
             end;

          for step:=c1 to c2 do
          begin

             case step of
               0 : ColorFunc:=Intensity;
               1 : ColorFunc:=RedComponent;
               2 : ColorFunc:=GreenComponent;
               3 : ColorFunc:=BlueComponent;
             end;

             for i:=0 to 2 do
               for j:=0 to 1 do
                begin
                   k:=j*3+i;
                   case j of
                     0 : Region1[k]:=ColorFunc(T[y-i]);
                     1 : Region1[k]:=ColorFunc(T2[y-i]);
                   end;
                end;

             for i:=0 to 2 do
               for j:=0 to 1 do
                begin
                   k:=j*3+i;
                   case j of
                     0 : Region4[k]:=ColorFunc(B[y+i]);
                     1 : Region4[k]:=ColorFunc(B2[y+i]);
                   end;
                end;

             for i:=0 to 1 do
               for j:=0 to 2 do
                begin
                   k:=j*2+i;
                   case j of
                     0 : Region2[k]:=ColorFunc(C[y+i]);
                     1 : Region2[k]:=ColorFunc(T[y+i]);
                     2 : Region2[k]:=ColorFunc(T2[y+i]);
                   end;
                end;

             for i:=0 to 1 do
               for j:=0 to 2 do
                begin
                   k:=j*2+i;
                   case j of
                     0 : Region3[k]:=ColorFunc(C[y-i]);
                     1 : Region3[k]:=ColorFunc(B[y-i]);
                     2 : Region3[k]:=ColorFunc(B2[y-i]);
                   end;
                end;

             MeanAndStdDev(Region1, Mean[1], StdDev[1]);
             MeanAndStdDev(Region2, Mean[2], StdDev[2]);
             MeanAndStdDev(Region3, Mean[3], StdDev[3]);
             MeanAndStdDev(Region4, Mean[4], StdDev[4]);

             minStdDev:=StdDev[1];
             n:=1;
             for i:=2 to 4 do
              if StdDev[i]<minStdDev then
               begin
                 minStdDev:=StdDev[i];
                 n:=i;
               end;

             Color[step]:= Round(Mean[n]);

          end;//step

          Case ColorMode of
             cmGray  : O[y]:=Gray32( Color[0] );
             cmRed   : O[y]:=Color32(Color[1],0,0);
             cmGreen : O[y]:=Color32(0,Color[2],0);
             cmBlue  : O[y]:=Color32(0,0,Color[3]);
             cmColor : O[y]:=Color32(Color[1],Color[2],Color[3]);
          end;

    end;//y
  end;//x

finally
 tBufr.Free;
end;

end;



procedure Do1DConvolution(Dst, Src: TBitmap32; Const Kernel: array of Integer );
var
 i,x,y,n     : Integer;
 pDst        : PColor32;
 C, R, G, B  : Integer;
 pPixelArray : array of PColor32;
 z           : integer;
begin
 Dst.Assign(Src);

 n:=Trunc( (Length(Kernel)-1)/2 );//полуширина 0-базированная
 z:=Kernel[High(Kernel)];
 SetLength(pPixelArray, Length(Kernel)-1);

With Src do
    for y:= 0 to Height-1 do begin
        //установка указателей на пикселы
        for i:= 0 to High(pPixelArray) do pPixelArray[i]:=PixelPtr[i, y];
        pDst:= Dst.PixelPtr[n, y];

        for x:= n to Width-1-n do begin

            R:= 0;  G:= 0;  B:= 0;
            for i:= 0 to High(pPixelArray) do begin
                C:= pPixelArray[i]^;
                R:= R + Kernel[i]*((C and $00FF0000) shr 16 );
                G:= G + Kernel[i]*((C and $0000FF00) shr 8 );
                B:= B + Kernel[i]*( C and $000000FF );
            end;

            R:=Set255( R div z);
            G:=Set255( G div z);
            B:=Set255( B div z);

            pDst^:= $FF000000 + R shl 16 + G shl 8 + B;

            for i:= 0 to High(pPixelArray) do Inc(pPixelArray[i]);
            Inc(pDst);
        end;//x
    end;//y

 Dst.Changed;
end;


procedure DoHough(Dst, Src: TBitmap32);
var
 x,y         : Integer;
 C           : Integer;
 r,t         : Integer;
 TArr        : Array [0..180*HoughTettaK-1] of double;

begin
 Dst.SetSize(HoughW*2-1, Length(TArr));
 Dst.Clear($00000000);

 //LUT градусы в радианы
 for t:=Low(TArr) to High(TArr) do TArr[t]:=DegToRad(t/HoughTettaK);

 With Src do
   for x:=0 to Width-1 do
     for y:=0 to Height-1 do begin
         C:=Pixel[x,y];
         if C>250 then //threshold
             for t:=Low(TArr) to High(TArr) do begin
                 r:=Round( x*cos(TArr[t]) + y*sin(TArr[t]) );
                 Dst.Pixel[r+(HoughW-1),t]:=Dst.Pixel[r+(HoughW-1),t] + C;
             end;
     end;

end;



procedure Do2DConvolutionXP(Dst, Src: TBitmap32; const Kernel: array of Integer);
var
 i,j,x,y     : Integer;
 pPixelArray : array of PColor32;
 pDst        : PColor32;
 C, R, G, B  : Integer;
 k,n,z       : Integer;
begin
 Dst.Assign(Src);
 n:=Round(Sqrt(Length(Kernel)-1)) div 2;//полуширина 0-базированная
 z:=Kernel[High(Kernel)];
 SetLength(pPixelArray, Length(Kernel)-1);

With Src do begin
    y := 0+n;
    While y<Height-n do begin
       // -------- установка указателей на пикселы ----------
       for j:=-n to n do
        for i:=-n to n do
         begin
           k:=(j+n)*(n*2+1)+(i+n);
           pPixelArray[k]:=PixelPtr[i+n,y+j];
         end;
       pDst:= Dst.PixelPtr[n, y];

       x:=0+n;
       While x<Width-n do begin
           R := 0;  G := 0;   B := 0;

           for i := 0 to High(pPixelArray) do begin
               C := pPixelArray[i]^;
               R := R + ((C and $00FF0000) shr 16) * Kernel[i];
               G := G + ((C and $0000FF00) shr 8) * Kernel[i];
               B := B + ( C and $000000FF) * Kernel[i];
           end;


           R:=Set255( R div z);
           G:=Set255( G div z);
           B:=Set255( B div z);

           //R:=UpSet255( ABS(R div z) );
           //G:=UpSet255( ABS(G div z) );
           //B:=UpSet255( ABS(B div z) );

           pDst^:= $FF000000 + R shl 16 + G shl 8 + B;

           for i:= 0 to High(pPixelArray) do Inc(pPixelArray[i]);
           Inc(pDst);
           Inc(x);
       end;
       Inc(y);
    end;
end;
Dst.Changed;
end;


procedure Do2DConvolutionGrayEX( Src: TBitmap32; const Kernel: array of Integer);overload;
Var
 Dst  : TBitmap32; // temp bitmap
begin
 Dst:= TBitmap32.Create;
 try
   Do2DConvolutionGrayEX(Dst,Src,Kernel);
   Src.Assign(Dst);
 finally
   Dst.free;
 end;
end;


procedure Do2DConvolutionGrayEX(Dst, Src: TBitmap32; const Kernel: array of Integer);
var
 i,j,x,y     : Integer;
 pPixelArray : array of PColor32;
 pDst        : PColor32;
 k,n,z       : Integer;
 Res         : Integer;
begin
 Dst.Assign(Src);
 n:=Round(Sqrt(Length(Kernel)-1)) div 2;//полуширина 0-базированная
 z:=Kernel[High(Kernel)];
 SetLength(pPixelArray, Length(Kernel)-1);

With Src do begin
    y := 0+n;
    While y<Height-n do begin

       for j:=-n to n do //-------- установка указателей на пикселы ----------
        for i:=-n to n do
         begin
           k:=(j+n)*(n*2+1)+(i+n);
           pPixelArray[k]:=PixelPtr[i+n,y+j];
         end;
       pDst:= Dst.PixelPtr[n,y];

       x:=0+n;
       While x<Width-n do begin
           Res:=0;
           for i := 0 to High(pPixelArray) do
               Res:=Res + pPixelArray[i]^ * Kernel[i];

           Res:=Res div z;//не нормирую на 255!!!
           if Res<0 then Res:=0;
           pDst^:= Res;

           for i:= 0 to High(pPixelArray) do Inc(pPixelArray[i]);
           Inc(pDst);
           Inc(x);
       end;
       Inc(y);
    end;
end;
Dst.Changed;
end;


procedure Do2DConvolutionXP( Src: TBitmap32; const Kernel: array of Integer);
Var
 Dst  : TBitmap32; // temp bitmap
begin
 Dst:= TBitmap32.Create;
 try
   Do2DConvolutionXP(Dst,Src,Kernel);
   Src.Assign(Dst);
 finally
   Dst.free;
 end;
end;

procedure Color2Gray( Src: TBitmap32);
Var
 Dst : TBitmap32; // temp bitmap
begin
 Dst:=TBitmap32.Create;
 try
   ColorToGrayscale(Dst, Src);
   Src.Assign(Dst);
 finally
   Dst.free;
 end;
end;

procedure Color2GrayEX( Src: TBitmap32);
var
  i: Integer;
  S: PColor32;
begin
  S := @Src.Bits[0];
  i:=Src.Width*Src.Height;
  While i>0 do begin
     S^ := Intensity(S^);
     Inc(S); Dec(i);
  end;
  Src.Changed;
end;


procedure Hough2Lines( Dst: TBitmap32; AMaxArr : TPoint2fArr );
Var
 x,y         : Integer;
 r, tetta    : single;
 i           : integer;
begin
 for i:=0 to High(AMaxArr) do begin
    r:=AMaxArr[i].x;
    tetta:=AMaxArr[i].y;
    if Abs( sin(tetta) ) > 0.7 then  //~45 град
       for x:=0 to Dst.Width-1 do begin
           y:=Round( ( r-x*cos(tetta) )/sin(tetta) );
           if (y>=960)or(y<0) then continue;
           Dst.Pixel[x,y]:=clLime32;
       end
    else
       for y:=0 to Dst.Height-1 do begin
           x:=Round( ( r-y*sin(tetta) )/cos(tetta) );
           if (x>=1280)or(x<0) then continue;
           Dst.Pixel[x,y]:=clLime32;
       end;
 end;

end;


procedure GetCrosses(Var AMaxArr : TPoint2fArr; Var  ACrossPointArr : T2DPoint2fArr );
Var
 x,y      : Integer;
 r, tetta : single;
 tetta0   : single;

 i,j,k,n  : integer;

 Arr1     : TPoint2fArr;
 Arr2     : TPoint2fArr;
 TempArr  : TPoint2fArr;

 dt       : single;
 Len      : integer;

 down     : double;
begin
 Arr1:=nil;
 Arr2:=nil;

 dt:=0.2;


 SortTetta(AMaxArr); // в Arr1 всегда попадут "горизонтальные" линии

 tetta0:=AMaxArr[0].y;
 SetLength(Arr1,1);
 Arr1[0]:=AMaxArr[0];

 //сортировка по углам  t на 2 множества
 for n:=1 to High(AMaxArr) do begin
     tetta:=AMaxArr[n].y;
     if ABS( ABS(cos(tetta)) - ABS(cos(tetta0)) )<dt then
        begin
           Len:=Length(Arr1);
           SetLength(Arr1,Len+1);
           Arr1[Len]:=AMaxArr[n];
        end
     else
        begin
           Len:=Length(Arr2);
           SetLength(Arr2,Len+1);
           Arr2[Len]:=AMaxArr[n];
        end;
 end;


 SortPoint2fArrSpecial(Arr1,cdY);

 SortPoint2fArrSpecial(Arr2,cdX);

 //SortPoint2fArr(Arr1,cdX);
 //SortPoint2fArr(Arr2,cdX);


 SetLength(ACrossPointArr, Length(Arr1), Length(Arr2));
 for i:=0 to High(ACrossPointArr) do //arr1
  for j:=0 to High(ACrossPointArr[i]) do begin //arr2
      down:=cos(Arr1[i].y)*sin(Arr2[j].y)-cos(Arr2[j].y)*sin(Arr1[i].y);

      ACrossPointArr[i,j].x:=( Arr1[i].x*sin(Arr2[j].y)-Arr2[j].x*sin(Arr1[i].y) )/down;
      ACrossPointArr[i,j].y:=( Arr2[j].x*cos(Arr1[i].y)-Arr1[i].x*cos(Arr2[j].y) )/down;
  end;

end;



procedure DoSobel( Dst, Src: TBitmap32);
var
 i          : integer;
 BmpX, BmpY : TBitmap32;
 D, Sx , Sy     : PColor32;
 Cx, Rx, Gx, Bx : Integer;
 Cy, Ry, Gy, By : Integer;
 R, G, B        : Integer;
begin
 BmpX:= TBitmap32.Create;
 BmpY:= TBitmap32.Create;

 try
   Do2DConvolutionXP(BmpX,Src,SobelHx);
   Do2DConvolutionXP(BmpY,Src,SobelHy);

   Sx := @BmpX.Bits[0];
   Sy := @BmpY.Bits[0];
   D  := @Dst.Bits[0];

   i:=Src.Width*Src.Height;
   While i>0 do begin

       Cx := Sx^;
       Rx :=(Cx and $00FF0000)shr 16;
       Gx :=(Cx and $0000FF00)shr 8;
       Bx :=(Cx and $000000FF);

       Cy := Sy^;
       Ry :=(Cy and $00FF0000)shr 16;
       Gy :=(Cy and $0000FF00)shr 8;
       By :=(Cy and $000000FF);

       R:=UpSet255( Round(Hypot(Rx,Ry)) );
       G:=UpSet255( Round(Hypot(Gx,Gy)) );
       B:=UpSet255( Round(Hypot(Bx,By)) );

       D^:= $FF000000 + R shl 16 + G shl 8 + B;

    Inc(Sx); Inc(Sy); Inc(D); Dec(i);
   end;

 Dst.Changed;

 finally
   BmpX.free;
   BmpY.free;
 end;

end;

procedure DoSuppression( Dst, Direction, Grad: TBitmap32);
var
 x,y,n,id : integer;
 B        : byte;
 isMax    : boolean;
 C        : integer;
begin

 CheckParams(Dst,Grad);
 Dst.Clear($00000000);

 n:=1;
 With Grad do
 for x:= n to Width-1-n do
 for y:= n to Height-1-n do begin
     C:=Pixel[x,y];
     id:=Direction.Pixel[x,y];

     Case id of
        1 : isMax:=( C>=Pixel[x,  y+1] ) and ( C>=Pixel[x,  y-1] );
        2 : isMax:=( C>=Pixel[x+1,y-1] ) and ( C>=Pixel[x-1,y+1] );
        3 : isMax:=( C>=Pixel[x+1,y  ] ) and ( C>=Pixel[x-1,  y] );
        4 : isMax:=( C>=Pixel[x+1,y+1] ) and ( C>=Pixel[x-1,y-1] );
     end;

     if not isMax
      then Dst.Pixel[x,y]:=$00000000
       else  Dst.Pixel[x,y]:=C;
 end;

end;


procedure GetLocalMax( Dst, Src : TBitmap32; Var AMaxArr : TPoint2fArr; Threshold    : integer);
var
 x,y,n,id,i,j : integer;
 x_,y_        : integer;
 B            : byte;
 isMax        : boolean;
 C            : TColor32;
 Tmp          : TBitmap32;
 Len          : integer;
 TempMaxArr   : TPoint2fArr;
 r,t          : single;

begin

 //threshold - сделать автоматическое определение величины

 AMaxArr:=nil;
 TempMaxArr:=nil;

 n:=15;// полуширина области поиска

 Dst.SetSize(Src.Width, Src.Height);
 Dst.Clear($00000000);

 With Src do
    for x:=0+ n to Width-1 -n do begin

        for y:=0 to Height-1 do begin
            C:=Pixel[x,y];
            if C<Threshold then continue;
            isMax:=true;

            // проверяем на локальный макс
            for i:=-n to n do begin
                for j:=-n to n do begin
                    X_:=x+i;

                    Y_:=y+j;
                    if Y_<0 then begin
                       Y_:=Height+Y_;
                       X_:=2*(HoughW-1) - X_;
                    end
                    else
                    if Y_> Height-1 then begin
                       Y_:=Y_-Height;
                       X_:=2*(HoughW-1) - X_;
                    end;

                    if C<Pixel[X_,Y_] then begin
                       isMax:=false;
                       break;//j
                    end;
                end;
                if not isMax then break;//i
            end;

            if not isMax then
               Dst.Pixel[x,y]:=$00000000
            else
               begin
                  Dst.Pixel[x,y]:=clWhite32;//C;
                  Len:=Length(TempMaxArr);
                  SetLength(TempMaxArr,Len+1);
                  //TempMaxArr[Len]:=Point2f(x, DegToRad(y/HoughTettaK) );
                  TempMaxArr[Len]:=Point2f(x,y);
               end;
        end;
    end;

 for i:=0 to High(TempMaxArr) do begin
     // r,t conwertion
     TempMaxArr[i].x :=TempMaxArr[i].x -(HoughW-1);
     TempMaxArr[i].y :=DegToRad(TempMaxArr[i].y/HoughTettaK);
 end;

 AMaxArr:=Copy(TempMaxArr);

end;


procedure CyclHoughBitmat( Dst, Src : TBitmap32);
var
 x,y,n : integer;
 x_,y_ : integer;
 C     : TColor32;
begin
 n:=15;// полуширина области поиска

 Dst.SetSize(Src.Width, Src.Height+n*2);
 Dst.Clear($00000000);
 Dst.Draw(0,n,Src);

 With Src do
    for x:=0 to Width-1  do begin

        for y:=n-1 downto 0 do begin

            X_:=x;
            Y_:=y - n;;
            if Y_<0 then begin
               Y_:=Height+Y_;
               X_:=2*(HoughW-1) - X_;
            end;

            C:=Pixel[X_,Y_];
            Dst.Pixel[x,y]:=C;

        end;

        for y:=0 to n-1 do begin

            X_:=x;
            Y_:=Height + y;
            if Y_> Height-1 then begin
               Y_:=Y_-Height;
               X_:=2*(HoughW-1) - X_;
            end;
            C:=Pixel[X_,Y_];
            Dst.Pixel[x,n+Height+y]:=C;
        end;

    end;

end;




procedure DoSobelDirect( Direction, Dst, Src: TBitmap32; const Kernelx,Kernely: array of Integer);
var
 Rx, Gx, Bx  : Integer;
 Ry, Gy, By  : Integer;
 R, G, B     : Integer;
 C           : TColor32;
 i,j,x,y     : Integer;
 pPixelArray : array of PColor32;

 pDst        : PColor32;
 pDir        : PColor32;

 k,n,z       : Integer;
 Int,IntX,IntY : Integer;
 angle       : single;
 angleId     : integer;
begin
 CheckParams(Dst,Src);
 Dst.Clear($00000000);

 CheckParams(Direction,Src);
 Direction.Clear($00000000);

 n:=Round(Sqrt(Length(Kernelx)-1)) div 2;//полуширина 0-базированная
 z:=Kernelx[High(Kernelx)];
 SetLength(pPixelArray, Length(Kernelx)-1);

With Src do begin
    y := 0+n;
    While y<Height-n do begin

       for j:=-n to n do // -------- установка указателей на пикселы ----------
        for i:=-n to n do
         begin
           k:=(j+n)*(n*2+1)+(i+n);
           pPixelArray[k]:=PixelPtr[i+n,y+j];
         end;
       pDst:= Dst.PixelPtr[n, y];
       pDir:= Direction.PixelPtr[n, y];

       x:=0+n;
       While x<Width-n do begin
           //Rx:=0; Gx:=0; Bx:=0;
           //Ry:=0; Gy:=0; By:=0;
           IntX:=0; IntY:=0;

           for i := 0 to High(pPixelArray) do begin
               C := pPixelArray[i]^;

               //R:=(C and $00FF0000)shr 16;
               //G:=(C and $0000FF00)shr 8;
               //B:=(C and $000000FF);
               //Rx:= Rx+R*Kernelx[i]; Gx:= Gx+G*Kernelx[i]; Bx:= Bx+B*Kernelx[i];
               //Ry:= Ry+R*Kernely[i]; Gy:= Gy+G*Kernely[i]; By:= By+B*Kernely[i];

               //Int:=Intensity(C);
               Int:=C;//GrayInt

               IntX:=IntX+Int*Kernelx[i];
               IntY:=IntY+Int*Kernely[i];

           end;
           //R:=UpSet255( Round( Hypot(Rx,Ry)/z ) );
           //G:=UpSet255( Round( Hypot(Gx,Gy)/z ) );
           //B:=UpSet255( Round( Hypot(Bx,By)/z ) );
           //pDst^:= $FF000000 + R shl 16 + G shl 8 + B;

           C:= Round( Hypot(IntX,IntY) );// не нормирую на 255 и вес Z
           //C:=ABS(IntX) + ABS(IntY);
           pDst^:=C;

           angle:=ABS(ArcTan2(IntY,IntX));
           angleId:=0;
           if angle<=Pi*1/8 then angleId:=3 else
           if angle<=Pi*3/8 then angleId:=2 else
           if angle<=Pi*5/8 then angleId:=1 else
           if angle<=Pi*7/8 then angleId:=4 else
           if angle >Pi*7/8 then angleId:=3 else
           ShowMessage('Error at angle computation');
           //pDir^:=ColorArr[angleId];
           pDir^:=angleId;

           for i:= 0 to High(pPixelArray) do Inc(pPixelArray[i]);
           Inc(pDst); Inc(pDir);
           Inc(x);
       end;
       Inc(y);
    end;
end;
Dst.Changed;
Direction.Changed;
end;


procedure GrayEX2Gray(Src : TBitmap32);
Var
 i          : integer;
 S          : PColor32;
 C          : Cardinal;
 max        : Cardinal;
 k          : double;

begin
   max:=0;

   S := @Src.Bits[0];
   i:=Src.Width*Src.Height;
   While i>0 do begin
       if S^ > max then max:=S^;
       Inc(S); Dec(i);
   end;

   k:=max/255;

   S := @Src.Bits[0];
   i:=Src.Width*Src.Height;
   While i>0 do begin
         C:=Set255( Round(S^/k) );
         S^:= $00000000 + C shl 16 + C shl 8 + C;
         Inc(S); Dec(i);
   end;

   Src.Changed;

end;


//--упорядочивание точек по TCompareDir = (cdX,cdY);
procedure SortPoint2fArr(Var Arr : TPoint2fArr; Direction : TCompareDir);
{
Type
 TCompareFunc = function(i, j : Integer) : Integer;
Var
 Compare : TCompareFunc;
}
    function CompareX(i, j : Integer) : Integer;
    begin
       if Arr[i].x>Arr[j].x then Result:=1
       else
       if Arr[i].x<Arr[j].x then Result:=-1
       else
       Result:=0;
    end;

    function CompareY(i, j : Integer) : Integer;
    begin
       if Arr[i].y>Arr[j].y then Result:=1
       else
       if Arr[i].y<Arr[j].y then Result:=-1
       else
       Result:=0;
    end;

    procedure Swap(i, j : Integer);
    Var
     P : TPoint2f;
    begin
       P:=Arr[i];
       Arr[i]:=Arr[j];
       Arr[j]:=P;
    end;

    procedure DoSort(From, Count: Integer);
          procedure   Sort(iL, iR: Integer);
          var
           L, R, M : Integer;
          begin
              repeat
                  L := iL;
                  R := iR;
                  M := ( L + R ) div 2;
                  repeat
                      //while Compare(From + L, From + M) < 0 do Inc(L);
                      //while Compare(From + M, From + R) < 0 do Dec(R);
                      Case Direction of
                       cdX : begin
                                while CompareX(From + L, From + M) < 0 do Inc(L);
                                while CompareX(From + M, From + R) < 0 do Dec(R);
                             end;
                       cdY : begin
                                while CompareY(From + L, From + M) < 0 do Inc(L);
                                while CompareY(From + M, From + R) < 0 do Dec(R);
                             end;
                      end;

                      if L <= R then begin
                          Swap(From + L, From + R);
                          if M = L then M := R
                          else
                          if M = R then M := L;
                          Inc(L);
                          Dec(R);
                      end;
                  until L > R;

                  if (R - iL) > (iR - L) then begin
                      if L < iR then Sort(L, iR);
                      iR := R;
                  end
                  else begin
                      if iL < R then Sort(iL, R);
                      iL := L;
                  end;

              until iL >= iR;
          end;
    begin
        if Count > 1 then Sort(0, Count - 1);
    end;

begin
{
  Case Direction of
   cdX : Compare:=CompareX;
   cdY : Compare:=CompareY;
  end;
}
  DoSort(0,Length(Arr));
end;


procedure SortPoint2fArrSpecial(Var Arr : TPoint2fArr; Direction : TCompareDir);

    function CompareX(i, j : Integer) : Integer;
    begin
       if Arr[i].x/cos(Arr[i].y) > Arr[j].x/cos(Arr[j].y) then Result:=1
       else
       if Arr[i].x/cos(Arr[i].y) < Arr[j].x/cos(Arr[j].y) then Result:=-1
       else
       Result:=0;
    end;

    function CompareY(i, j : Integer) : Integer;
    begin
       if Arr[i].x/sin(Arr[i].y) > Arr[j].x/sin(Arr[j].y) then Result:=1
       else
       if Arr[i].x/sin(Arr[i].y) < Arr[j].x/sin(Arr[j].y) then Result:=-1
       else
       Result:=0;
    end;

    procedure Swap(i, j : Integer);
    Var
     P : TPoint2f;
    begin
       P:=Arr[i];
       Arr[i]:=Arr[j];
       Arr[j]:=P;
    end;

    procedure DoSort(From, Count: Integer);
          procedure   Sort(iL, iR: Integer);
          var
           L, R, M : Integer;
          begin
              repeat
                  L := iL;
                  R := iR;
                  M := ( L + R ) div 2;
                  repeat
                      Case Direction of
                       cdX : begin
                                while CompareX(From + L, From + M) < 0 do Inc(L);
                                while CompareX(From + M, From + R) < 0 do Dec(R);
                             end;
                       cdY : begin
                                while CompareY(From + L, From + M) < 0 do Inc(L);
                                while CompareY(From + M, From + R) < 0 do Dec(R);
                             end;
                      end;

                      if L <= R then begin
                          Swap(From + L, From + R);
                          if M = L then M := R
                          else
                          if M = R then M := L;
                          Inc(L);
                          Dec(R);
                      end;
                  until L > R;

                  if (R - iL) > (iR - L) then begin
                      if L < iR then Sort(L, iR);
                      iR := R;
                  end
                  else begin
                      if iL < R then Sort(iL, R);
                      iL := L;
                  end;

              until iL >= iR;
          end;
    begin
        if Count > 1 then Sort(0, Count - 1);
    end;

begin
  DoSort(0,Length(Arr));
end;




procedure SortTetta(Var Arr : TPoint2fArr); //(r,t)

    function Compare(i, j : Integer) : Integer;
    begin
       if ABS(cos(Arr[i].y)) > ABS(cos(Arr[j].y)) then Result:=1
       else
       if ABS(cos(Arr[i].y)) < ABS(cos(Arr[j].y)) then Result:=-1
       else
       Result:=0;
    end;

    procedure Swap(i, j : Integer);
    Var
     P : TPoint2f;
    begin
       P:=Arr[i];
       Arr[i]:=Arr[j];
       Arr[j]:=P;
    end;

    procedure DoSort(From, Count: Integer);
          procedure   Sort(iL, iR: Integer);
          var
           L, R, M : Integer;
          begin
              repeat
                  L := iL;
                  R := iR;
                  M := ( L + R ) div 2;
                  repeat
                      while Compare(From + L, From + M) < 0 do Inc(L);
                      while Compare(From + M, From + R) < 0 do Dec(R);

                      if L <= R then begin
                          Swap(From + L, From + R);
                          if M = L then M := R
                          else
                          if M = R then M := L;
                          Inc(L);
                          Dec(R);
                      end;
                  until L > R;

                  if (R - iL) > (iR - L) then begin
                      if L < iR then Sort(L, iR);
                      iR := R;
                  end
                  else begin
                      if iL < R then Sort(iL, R);
                      iL := L;
                  end;

              until iL >= iR;
          end;
    begin
        if Count > 1 then Sort(0, Count - 1);
    end;

begin
  DoSort(0,Length(Arr));
end;




Procedure Circle(Dst : TBitmap32; X, Y, Radius,Color: TColor32);
Var
 a, b, r2   : Integer;
 i : integer;
Begin
 R2 := Sqr(Radius);
 a:=0;
 While a<=Radius do begin
       b:= Round(Sqrt(R2 - Sqr(a)));
       for i:=-b to b do begin
           Dst.PixelS[x+i,y-a]:=color;
           Dst.PixelS[x+i,y+a]:=color;
       end;
       a:=a+1;
 end;
end;


procedure Subtruct(Bmp1, Bmp2: TBitmap32);
var
  I: Integer;
  p1, p2: PColor32;
begin
  p1 := @Bmp1.Bits[0];
  p2 := @Bmp2.Bits[0];
  for I := 0 to Bmp1.Width * Bmp1.Height - 1 do
  begin
    p1^ := ColorSub(P1^,P2^);
    Inc(P1); Inc(P2);
  end;
  Bmp1.Changed;
end;


End.
