Skip to content

Commit

Permalink
Fix ThresholdSauvola
Browse files Browse the repository at this point in the history
  • Loading branch information
@ollydev
ollydev committed Oct 12, 2023
1 parent 953879b commit da1b6de
Show file tree
Hide file tree
Showing 4 changed files with 160 additions and 104 deletions.
6 changes: 5 additions & 1 deletion Source/Simba.lpi
View file
Original file line number Diff line number Diff line change
Expand Up @@ -340,7 +340,7 @@
<PackageName Value="LCL"/>
</Item5>
</RequiredPackages>
<Units Count="154">
<Units Count="155">
<Unit0>
<Filename Value="Simba.lpr"/>
<IsPartOfProject Value="True"/>
Expand Down Expand Up @@ -1041,6 +1041,10 @@
<Filename Value="clearpixelaa.inc"/>
<IsPartOfProject Value="True"/>
</Unit153>
<Unit154>
<Filename Value="simba.image_integral.pas"/>
<IsPartOfProject Value="True"/>
</Unit154>
</Units>
</ProjectOptions>
<CompilerOptions>
Expand Down
149 changes: 49 additions & 100 deletions Source/image/simba.image.pas
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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;

Expand Down
99 changes: 99 additions & 0 deletions Source/image/simba.image_integral.pas
View file
Original file line number Diff line number Diff line change
@@ -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.

10 changes: 7 additions & 3 deletions Source/script/imports/simbaclasses/simba.import_class_image.pas
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

(*
Expand Down Expand Up @@ -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);
Expand Down

0 comments on commit da1b6de

Please sign in to comment.