From da1b6de8898d4aeb6ac78f653ea78459d26bc045 Mon Sep 17 00:00:00 2001 From: Olly Date: Fri, 13 Oct 2023 00:16:01 +0100 Subject: [PATCH] Fix ThresholdSauvola --- Source/Simba.lpi | 6 +- Source/image/simba.image.pas | 149 ++++++------------ Source/image/simba.image_integral.pas | 99 ++++++++++++ .../simbaclasses/simba.import_class_image.pas | 10 +- 4 files changed, 160 insertions(+), 104 deletions(-) create mode 100644 Source/image/simba.image_integral.pas diff --git a/Source/Simba.lpi b/Source/Simba.lpi index ba0d59a1f..c124a7324 100644 --- a/Source/Simba.lpi +++ b/Source/Simba.lpi @@ -340,7 +340,7 @@ - + @@ -1041,6 +1041,10 @@ + + + + diff --git a/Source/image/simba.image.pas b/Source/image/simba.image.pas index 5020a4044..e959af7e0 100644 --- a/Source/image/simba.image.pas +++ b/Source/image/simba.image.pas @@ -192,7 +192,7 @@ TSimbaImage = class(TSimbaBaseClass) function ToMatrix(X1, Y1, X2, Y2: Integer): TIntegerMatrix; overload; function ThresholdAdaptive(Alpha, Beta: Byte; AInvert: Boolean; Method: ESimbaImageThreshMethod; K: Integer): TSimbaImage; - function ThresholdSauvola(Window: Integer; AInvert: Boolean; K: Single): TSimbaImage; + function ThresholdSauvola(Radius: Integer; R: Single = 128; K: Single = 0.5): TSimbaImage; function RowPtrs: TSimbaImageRowPtrs; @@ -229,7 +229,7 @@ implementation simba.arraybuffer, simba.geometry, simba.tpa, simba.encoding, simba.compress, simba.nativeinterface, simba.singlematrix, - simba.image_lazbridge, simba.rgbsumtable; + simba.image_lazbridge, simba.rgbsumtable, simba.image_integral; function GetDistinctColor(const Color, Index: Integer): Integer; inline; const @@ -2216,115 +2216,64 @@ function TSimbaImage.ThresholdAdaptive(Alpha, Beta: Byte; AInvert: Boolean; Meth end; end; -function TSimbaImage.ThresholdSauvola(Window: Integer; AInvert: Boolean; K: Single): TSimbaImage; -var - SumTable: TDoubleMatrix; - SumTableSquared: TDoubleMatrix; - - function CreateSumTable(Matrix: TByteMatrix): TDoubleMatrix; - var - W, H, X, Y: Integer; - begin - SetLength(Result, Matrix.Height, Matrix.Width); - - W := Result.Width - 1; - H := Result.Height - 1; - - for Y := 0 to H do - Result[0, Y] := Matrix[0, Y]; - - for Y := 1 to H do - for X := 0 to W do - Result[Y, X] := Matrix[Y,X] + Result[Y-1, X]; - - for Y := 0 to H do - for X := 1 to W do - Result[Y, X] += Result[Y, X-1]; - end; - - function CreateSumTableSquared(Matrix: TByteMatrix): TDoubleMatrix; - var - W, H, X, Y: Integer; - begin - SetLength(Result, Matrix.Height, Matrix.Width); - - W := Result.Width - 1; - H := Result.Height - 1; - - for Y := 0 to H do - Result[0, Y] := Sqr(Matrix[0, Y]); - - for Y := 1 to H do - for X := 0 to W do - Result[Y, X] := Sqr(Matrix[Y,X]) + Result[Y-1, X]; - - for Y := 0 to H do - for X := 1 to W do - Result[Y, X] += Result[Y, X-1]; - end; - - function QuerySumTable(Table: TDoubleMatrix; X1, Y1, X2, Y2: Integer): Double; - begin - Result := Table[Y2, X2]; - - if (Y1 > 0) then - Result := Result - Table[Y1 - 1, X2]; - if (X1 > 0) then - Result := Result - Table[Y2, X1 - 1]; - if (Y1 > 0) and (X1 > 0) then - Result := Result + Table[Y1 - 1, X1 - 1]; - end; - - function GetMean(B: TBox): Single; - begin - Result := QuerySumTable(SumTable, B.X1, B.Y1, B.X2, B.Y2) / B.Area(); - end; - - function GetStdDev(B: TBox; Mean: Double): Single; - begin - Result := QuerySumTable(SumTableSquared, B.X1, B.Y1, B.X2, B.Y2) - (Sqr(Mean) * B.Area); - Result := Sqrt(Result / B.Area); - end; - +{ + Radius = Window size + R = dynamic range of standard deviation (default = 128) + K = constant value in range 0.2..0.5 (default = 0.5) +} +function TSimbaImage.ThresholdSauvola(Radius: Integer; R: Single = 128; K: Single = 0.5): TSimbaImage; +const + HIT: TColorBGRA = (B: 255; G: 255; R: 255; A: 0); + MISS: TColorBGRA = (B: 0; G: 0; R: 0; A: 0); var - X, Y, W, H: Integer; - Mean, StdDev, Thresh: Single; - B: TBox; - Matrix: TByteMatrix; - Background, Foreground: TColorBGRA; + Mat: TByteMatrix; + Integral: TSimbaIntegralImageF; + X,Y,W,H: Integer; + Left, Right, Top, Bottom: Integer; + Count: Integer; + Sum, SumSquares: Double; + Mean, Stdev, Threshold: Double; begin Result := TSimbaImage.Create(FWidth, FHeight); - Background := Default(TColorBGRA); - Foreground := Default(TColorBGRA); - Foreground.R := 255; - if AInvert then - Swap(Integer(Background), Integer(Foreground)); + Mat := Self.ToGreyMatrix(); + Mat.GetSize(W, H); - Matrix := ToGreyMatrix(); + Dec(W); + Dec(H); - SumTable := CreateSumTable(Matrix); - SumTableSquared := CreateSumTableSquared(Matrix); - - W := FWidth - 1; - H := FHeight - 1; + Radius := (Radius - 1) div 2; + Integral := TSimbaIntegralImageF.Create(Mat); for Y := 0 to H do for X := 0 to W do begin - B.X1 := Max(X-Window, 0); - B.Y1 := Max(Y-Window, 0); - B.X2 := Min(X+Window, W); - B.Y2 := Min(Y+Window, H); - - Mean := GetMean(B); - StdDev := GetStdDev(B, Mean); - Thresh := Mean * (1.0 + K * ((StdDev / 128.0) - 1.0)); - - if Matrix[Y, X] < Thresh then - Result.Data[Y*FWidth+X] := Foreground + Left := Max(X-W, 0); + Right := Min(X+W, W); + Top := Max(Y-W, 0); + Bottom := Min(Y+W, H); + + //Sum := 0; + //SumSquares := 0; + // + //for y := top to bottom do + // for x := left to right do + // begin + // Sum += mat[y,x]; + // SumSquares += mat[y,x]*mat[y,x]; + // end; + + Integral.Query(Left, Top, Right, Bottom, Sum, SumSquares); + + Count := (Bottom - Top + 1) * (Right - Left + 1); + Mean := Sum / Count; + Stdev := Sqrt((SumSquares / Count) - Sqr(Mean)); + Threshold := Mean * (1.0 + K * ((Stdev / R) - 1.0)); + + if Mat[Y, X] < Threshold then + Result.FData[Y * FWidth + X] := MISS else - Result.Data[Y*FWidth+X] := Background; + Result.FData[Y * FWidth + X] := HIT; end; end; diff --git a/Source/image/simba.image_integral.pas b/Source/image/simba.image_integral.pas new file mode 100644 index 000000000..7e2618592 --- /dev/null +++ b/Source/image/simba.image_integral.pas @@ -0,0 +1,99 @@ +{ + Author: Raymond van Venetië and Merlijn Wajer + Project: Simba (https://github.com/MerlijnWajer/Simba) + License: GNU General Public License (https://www.gnu.org/licenses/gpl-3.0) + + -------------------------------------------------------------------------- + + An integral image, also known as a summed area table, is a data structure + used for quick and efficient calculation of the sum of values in a + rectangular subset of an image. The value at each pixel in the integral image + is the sum of all the pixels above and to the left of it in the original + image. +} +unit simba.image_integral; + +{$i simba.inc} + +interface + +uses + Classes, SysUtils, + simba.mufasatypes, simba.image; + +type + TIntegralImageDataF = record Value: Double; ValueSqr: Double; end; + + TSimbaIntegralImageF = record + public + type + TData = record Value: Double; ValueSqr: Double; end; + TDataMatrix = array of array of TData; + public + Data: TDataMatrix; + + class function Create(const From: TByteMatrix): TSimbaIntegralImageF; static; + function Query(const Left, Top, Right, Bottom: Integer; out Value, ValueSqr: Double): TSimbaIntegralImageF; + end; + +implementation + +class function TSimbaIntegralImageF.Create(const From: TByteMatrix): TSimbaIntegralImageF; +var + X, Y: Integer; +begin + SetLength(Result.Data, From.Height, From.Width); + + // Compute the first row of the integral image + for X := 0 to From.Width - 1 do + begin + Result.Data[0, X].Value := From[0, X]; + Result.Data[0, X].ValueSqr := Sqr(From[0, X]); + if (X > 0) then + begin + Result.Data[0, X].Value += Result.Data[0, X - 1].Value; + Result.Data[0, X].ValueSqr += Result.Data[0, X - 1].ValueSqr; + end; + end; + + // Compute the first column of the integral image + for Y := 1 to From.Height - 1 do + begin + Result.Data[Y, 0].Value := From[Y, 0] + Result.Data[Y - 1, 0].Value; + Result.Data[Y, 0].ValueSqr := Sqr(From[Y, 0]) + Result.Data[Y - 1, 0].ValueSqr; + end; + + // Compute the rest of the integral image + for Y := 1 to From.Height - 1 do + for X := 1 to From.Width - 1 do + begin + Result.Data[Y, X].Value := From[Y, X] + Result.Data[Y - 1, X].Value + Result.Data[Y, X - 1].Value - Result.Data[Y - 1, X - 1].Value; + Result.Data[Y, X].ValueSqr := Sqr(From[Y, X]) + Result.Data[Y - 1, X].ValueSqr + Result.Data[Y, X - 1].ValueSqr - Result.Data[Y - 1, X - 1].ValueSqr; + end; +end; + +function TSimbaIntegralImageF.Query(const Left, Top, Right, Bottom: Integer; out Value, ValueSqr: Double): TSimbaIntegralImageF; +const + ZERO: TData = (Value: 0; ValueSqr: 0); +var + A,B,C,D: TData; +begin + A := ZERO; + B := ZERO; + C := ZERO; + D := ZERO; + + if (Left - 1 >= 0) and (Top - 1 >= 0) then + A := Data[Top - 1, Left - 1]; + if (Top - 1 >= 0) then + B := Data[Top - 1, Right]; + if (Left - 1 >= 0) then + C := Data[Bottom, Left - 1]; + D := Data[Bottom, Right]; + + Value := D.Value - B.Value - C.Value + A.Value; + ValueSqr := D.ValueSqr - B.ValueSqr - C.ValueSqr + A.ValueSqr; +end; + +end. + diff --git a/Source/script/imports/simbaclasses/simba.import_class_image.pas b/Source/script/imports/simbaclasses/simba.import_class_image.pas index f8e765a71..2f626ca47 100644 --- a/Source/script/imports/simbaclasses/simba.import_class_image.pas +++ b/Source/script/imports/simbaclasses/simba.import_class_image.pas @@ -412,11 +412,15 @@ procedure _LapeImage_ThresholdAdaptive(const Params: PParamArray; const Result: (* TImage.ThresholdSauvola ~~~~~~~~~~~~~~~~~~~~~~~ -> function TImage.ThresholdSauvola(Radius: Integer; AInvert: Boolean; k: Single): TImage; +> function TImage.ThresholdSauvola(Radius: Integer; R: Single = 128; K: Single = 0.5): TImage; + + Radius = Window size + R = dynamic range of standard deviation (default = 128) + K = constant value in range 0.2..0.5 (default = 0.5) *) procedure _LapeImage_ThresholdSauvola(const Params: PParamArray; const Result: Pointer); LAPE_WRAPPER_CALLING_CONV begin - PSimbaImage(Result)^ := PSimbaImage(Params^[0])^.ThresholdSauvola(PInteger(Params^[1])^, PBoolean(Params^[2])^, PSingle(Params^[3])^); + PSimbaImage(Result)^ := PSimbaImage(Params^[0])^.ThresholdSauvola(PInteger(Params^[1])^, PSingle(Params^[2])^, PSingle(Params^[3])^); end; (* @@ -1402,7 +1406,7 @@ procedure ImportSimbaImage(Compiler: TSimbaScript_Compiler); addGlobalFunc('function TImage.ToMatrix: TIntegerMatrix; overload', @_LapeImage_ToMatrix); addGlobalFunc('function TImage.ToMatrix(X1, Y1, X2, Y2: Integer): TIntegerMatrix; overload', @_LapeImage_ToMatrixEx); addGlobalFunc('function TImage.ThresholdAdaptive(Alpha, Beta: Byte; AInvert: Boolean; Method: EImageThreshMethod; k: Integer): TImage', @_LapeImage_ThresholdAdaptive); - addGlobalFunc('function TImage.ThresholdSauvola(Radius: Integer; AInvert: Boolean; k: Single): TImage', @_LapeImage_ThresholdSauvola); + addGlobalFunc('function TImage.ThresholdSauvola(Radius: Integer; R: Single = 128; K: Single = 0.5): TImage', @_LapeImage_ThresholdSauvola); addGlobalFunc('procedure TImage.Pad(Amount: Integer)', @_LapeImage_Pad); addGlobalFunc('procedure TImage.LoadFromFile(FileName: String); overload', @_LapeImage_LoadFromFile);