unit filterLocalDeviation;
(* ***** BEGIN LICENSE BLOCK *****
 * Copyright (C) 2004 Durand Emmanuel
 * Copyright (C) 2004 Burgel Eric
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 * Contact :
 *   filters@edurand.com
 *   filters@burgel.com
 *
 * ***** END LICENSE BLOCK ***** *)

{
 eburgel (filters@burgel.com)
}

interface
uses
  filter, fparameters, image, SysUtils;

type
  TFilterLocalDeviation = class(TFilter)
  public
    constructor Create; override;
    destructor Destroy; override;
    procedure Run(); override;
  private
    parameterPercentEcarttype   : TParameterInteger;
    parameterLocalDeviationSize : TParameterInteger ;
    parameterAnalysedZoneSize   : TParameterInteger ;
    parameterOffsetSize         : TParameterInteger ;
    parameterImageIn, parameterImageOut : TParameterImage;
    parameterOutputColor : TParameterString;
    imgIn, imgOut : PBitmap32 ;
    procedure _run();
  end;

implementation
uses
  imageIO, Math;

constructor TFilterLocalDeviation.Create;
begin
  inherited;
  parameterImageIn     := addParameterImage('inImage','explication...');
  parameterImageOut    := addParameterImage('outImage','explication...');

  parameterOutputColor:=addParameterString('outputColor','UNCHANGED/WHITE','WHITE');

  parameterPercentEcarttype   := addParameterInteger('percentLocalDeviation','bornInf/Sup=moyenne-/+pourcentageEcarttype',1,500,250);
  parameterLocalDeviationSize := addParameterInteger('localDeviationSize', 'Local Deviation Size', 2, 25, 3);
  parameterAnalysedZoneSize   := addParameterInteger('analysedZoneSize', 'Analysed Zone Size', 1, 25, 2) ;
  parameterOffsetSize         := addParameterInteger('offsetSize', 'Offset Size',  0, 25, 1) ;
end;

destructor TFilterLocalDeviation.Destroy;
begin
end;

procedure TFilterLocalDeviation.Run();
begin
  imgIn  := parameterImageIn.Image ;
  imgOut := parameterImageOut.Image ;
  if (imgIn<>nil) and (imgOut<>nil) then begin
    _run();
  end;
end;

{**************************************************
 In the following code, variables beginning with :
   D : Stands for LocalDeviation
   A : Analized zone
   O : Offset
***************************************************}
procedure TFilterLocalDeviation._run();
var
  PSrc, pBeginZoneDev, pEndZoneDev, pBZD, pEZD : PColor32Array;
  PDest : PColor32Array;
  pSrcCenter, pDestCenter : PColor32Array;
  trow, tcol : Integer;
  trowmax, tcolmax : Integer;
  maxX, maxY : Integer;
  x, y : Integer;
  tIntensity, tIntensityCount, tIntensityMoyenne : Integer;
  //variance, ecarttype : Single;
  delta : Integer;
  borneInf, borneSup : Integer;
  localDeviationSize, localDeviationWidth : Integer ;
  bOutpuColorWhite : Boolean;
  average, variance, stdDeviation : Single ;
  Dx, Dy : Integer ;
  s2, sum : Integer ;
  lIntensity : Integer ;
  N : Integer ;
  percentstdDev : Single ;
  // analised zone variables
  Ax, Ay, Asize, Awidth : Integer ;
begin
  bOutpuColorWhite:=True;
  localDeviationSize  := parameterLocalDeviationSize.Value ;
  localDeviationWidth := localDeviationSize*2 ;

  Asize  := parameterAnalysedZoneSize.Value ;
  Awidth := Asize * 2 ;

  percentstdDev := parameterPercentEcarttype.Value/100.;


  if parameterOutputColor.Value='UNCHANGED' then bOutpuColorWhite:=False;

  maxX := scanWidth(imgIn)-1 ;
  maxY := scanHeight(imgIn)-1 ;

  if (maxX>localDeviationWidth) and (maxY>localDeviationWidth) then begin
    for y:=localDeviationSize to maxY-localDeviationSize do Begin

      pSrcCenter := scanLine(imgIn, y) ;
      inc(pSrcCenter, localDeviationSize) ;
      pDestCenter := scanLine(imgOut, y) ;
      inc(pDestCenter, localDeviationSize) ;

      // Compute the local deviation for the first square of this line,
      // with the formula:
      // StdDev = square root of ( (S2 -  (SUM^2) /N) /(N-1) )  where
      //   S2 is the sum of the squares of the numbers
      //   SUM is the sum of the numbers
      //   N is the number of numbers

      sum := 0 ;
      s2 := 0 ;
      pBeginZoneDev := scanLine(imgIn, y-localDeviationSize) ;
      pEndZoneDev := pBeginZoneDev ;
      inc(pEndZoneDev, localDeviationWidth+1) ;
      for Dy := 0 to localDeviationWidth do Begin
        pSrc := scanLine(imgIn, y+Dy-localDeviationSize) ;
        for Dx:= 0 to localDeviationWidth do Begin
          lIntensity := Intensity(pSrc^[0]) ;
          inc(sum, lIntensity) ;
          inc(s2, lIntensity*lIntensity) ;
          inc(pSrc);
        end ;
      end ;

      N := localDeviationWidth+1 ;
      N := N * N ;
      for x:= localDeviationSize to maxX-localDeviationSize do Begin
        average := sum div N ;
        variance := (S2-(sum*sum)/N) / (N-1) ;
        stdDeviation := Sqrt(variance);

        borneInf := Round(average - stdDeviation * percentstdDev);
        borneSup := Round(average + stdDeviation * percentstdDev);

        // Now iterating over the local zone,
        // to detect points out of the standard deviation

        if false Then Begin
          Ay := 0 ;
          while Ay <= Awidth do Begin
          //for Ay := 0 to Awidth do Begin
            pSrc := scanLine(imgIn, y+AY-Asize) ;
            inc(pSrc, localDeviationSize-Asize+x-Asize) ;
            pDest := scanLine(imgOut, y+AY-Asize) ;
            inc(pDest, localDeviationSize-Asize+x-Asize) ;
            for Ax := 0 to Awidth do Begin
              lIntensity := Intensity(pSrc^[0]) ;
              if (lIntensity<borneInf) or (lIntensity>borneSup) then begin
                pDest[0]:=clWhite32;
              end else Begin
                //pDest[0]:=clTeal32;
              End ;
              inc(pSrc);
              inc(pDest);
            end ;
            inc(Ay,1) ;
          end ;
        end else Begin
          lIntensity := Intensity(pSrcCenter^[0]) ;
          if (lIntensity<borneInf) or (lIntensity>borneSup) then begin
            pDestCenter[0]:=clWhite32;
          end else Begin
            pDestCenter[0]:=clTeal32;
          End ;
          inc(pSrcCenter);
          inc(pDestCenter);
        end ;

        // Adjusting average and StdDev
        pBZD := pBeginZoneDev ;
        pEZD := pEndZoneDev ;
        for Dy := 0 to localDeviationWidth do Begin
          lIntensity := Intensity(pBZD^[0]) ;
          dec(sum, lIntensity) ;
          dec(s2, lIntensity*lIntensity) ;
          lIntensity := Intensity(pEZD^[0]) ;
          inc(sum, lIntensity) ;
          inc(s2, lIntensity*lIntensity) ;
          inc(pBZD, imgIn.Width) ;
          inc(pEZD, imgIn.Width) ;
        end ;
        inc(pBeginZoneDev) ;
        inc(pEndZoneDev) ;

      end ;
    end ;
  end;
  //Image.drawLine(imgOut,10,10,100,100,clRed32);
end;

end.