unit filterSUSAN;
(* ***** 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
 * Site :
 *   http://filters.sourceforge.net/
 *
 * ***** END LICENSE BLOCK ***** *)

{
 edurand (filters@edurand.com)
 eburgel (filters@burgel.com)
}

{
  Algorithm : http://www.fmrib.ox.ac.uk/~steve/susan/
}

interface
uses
  filter, fparameters, image, SysUtils, filterNonMaximaSuppression;

type
  TFilterSUSAN = class(TFilter)
  public
    constructor Create; override;
    destructor Destroy; override;
    procedure Run(); override;
    procedure setParameterString( const aName,aValue: String); override ;
    procedure setParameterInteger( const aName : String;  const aValue : Int64); override;
    procedure setParameterImage( const aName: String; const aImage: PBitmap32); override ;
  private
    parameterMode : TParameterString;
    parameterMaskType : TParameterString;
    parameterPixelBrightnessDifferenceThreshold : TParameterInteger;
    parameterPixelCentroidDifferenceThreshold : TParameterInteger;
    parameterImageIn, parameterImageOut : TParameterImage;
    parameterImageMask : TParameterImage;
    imageMask : PBitmap32;
    parameterPrecision : TParameterInteger;
    imageNonMaximaSuppressionOut : PBitmap32;
    filterNonMaximaSuppression : TFilterNonMaximaSuppression;
    procedure _run();
    procedure createImageMask();
    procedure deleteImages;
  end;

implementation
uses
  imageIO, Math;

constructor TFilterSUSAN.Create;
begin
  inherited;
  parameterImageIn:=addParameterImage('inImage','image to process');
  parameterImageOut:=addParameterImage('outImage','output image');
  parameterImageMask:=addParameterImage('maskImage','image of mask');
  parameterPrecision:=addParameterInteger('precision','',0,10,3);
  parameterMode:=addParameterString('mode','EDGE/CORNER','CORNER');
  parameterMaskType:=addParameterString('maskType','ROUND/SQUARE','ROUND');
  parameterPixelBrightnessDifferenceThreshold:=addParameterInteger('pixelBrightnessDifferenceThreshold','',0,255,27);
  parameterPixelCentroidDifferenceThreshold:=addParameterInteger('pixelCentroidDifferenceThreshold','',0,5,1);
  createImageMask();
  filterNonMaximaSuppression:=TFilterNonMaximaSuppression.Create;
end;

destructor TFilterSUSAN.Destroy;
begin
  inherited;
  deleteImages;
  image.freeImage(imageMask);
  filterNonMaximaSuppression.Free;
end;

procedure TFilterSUSAN.deleteImages;
begin
  image.freeImage(imageNonMaximaSuppressionOut);
end;

procedure TFilterSUSAN.setParameterString( const aName,aValue: String);
begin
  inherited;
  if aName='maskType' then begin
    createImageMask;
  end;
end;

procedure TFilterSUSAN.setParameterInteger( const aName : String;  const aValue : Int64);
begin
  inherited;
  if aName='precision' then begin
    createImageMask;
  end;
end;

procedure TFilterSUSAN.setParameterImage( const aName: String; const aImage: PBitmap32);
begin
  if aName='inImage' then begin
    deleteImages;
    imageNonMaximaSuppressionOut:=createImageFromImage(aImage);
  end;
  inherited;
end;

procedure TFilterSUSAN.createImageMask();
var
  x,y : Integer;
  d, dMax : Integer;
  c : TColor32;
  maskSize : Integer;
begin
  image.freeImage(imageMask);
  maskSize:=parameterPrecision.Value*2+1;
  imageMask:=image.createImage(maskSize,maskSize);
  dMax:=11; // =11.56=3*3+3*3
  for y:=0 to imageMask.Height-1 do begin
    for x:=0 to imageMask.Width-1 do begin
      c:=clBlack32;
      if parameterMaskType.Value='SQUARE' then begin
        c:=clWhite32;
      end else begin
        d:=(x-3)*(x-3)+(y-3)*(y-3);
        if d<=dMax then begin
          c:=clWhite32;
        end;
      end;
      Image.setPixel(imageMask,x,y,c);
    end;
  end;
  setParameterImage('maskImage',imageMask);
end;

procedure TFilterSUSAN.Run();
begin
  if (parameterImageIn.Image<>nil) and (parameterImageOut.Image<>nil) then begin
    _run();
  end;
end;

procedure TFilterSUSAN._run();
var
  pSrc, PDest : PColor32Array;
  w,h :Integer;
  value, valueNucleus : Integer;
  tIntensity, tIntensityNucleus : Integer;
  roundImageMaskWidthDiv2 :Integer;
  roundImageMaskHeightDiv2 :Integer;
  imageSrcRow, imageSrcCol : Integer;
  imageSrcRowmax, imageSrcColmax : Integer;
  imageSrcNeighborRow, imageSrcNeighborCol : Integer;
  imageSrcNeighborRowmin, imageSrcNeighborColmin : Integer;
  imageSrcNeighborRowmax, imageSrcNeighborColmax : Integer;
  diff, t, n, nMax, g, n2 : Integer;
  gTo255 : Single;
  isDetected : boolean ;
  testCorner : boolean ;
  sumRow, sumCol : LongInt ;
  centroidRow, centroidCol : Single ;
  t2 : integer ;
  colorMask:TColor32;
  colInMask,rowInMask:Integer;
  mask : PBitmap32;
  imageDest : PBitmap32;
begin
  h:=parameterImageIn.Image.Height;
  w:=parameterImageIn.Image.Width;
  t:=parameterPixelBrightnessDifferenceThreshold.Value;
  t2:=parameterPixelCentroidDifferenceThreshold.Value ;
  mask:=parameterImageMask.Image;
  roundImageMaskWidthDiv2:=mask.Width div 2;
  roundImageMaskHeightDiv2:=mask.Height div 2;
  nMax:=0;
  for rowInMask:=0 to mask.Height-1 do begin
    for colInMask:=0 to mask.Width-1 do begin
      colorMask:=Image.getPixel(mask,colInMask,rowInMask);
      if colorMask=clWhite32 then begin
        Inc(nMax);
      end;
    end;
  end;
  if parameterMode.Value='EDGE' then begin
    g:=(3*nMax) div 4;
    testCorner:=False;
    imageDest:=parameterImageOut.Image;
  end else begin
    g:=nMax div 2;
    testCorner:=True;
    imageDest:=imageNonMaximaSuppressionOut;
  end;
  gTo255:=255 / g;
  pSrc:=parameterImageIn.Image.Bits;
  imageSrcRowmax:=h-1;
  imageSrcColmax:=w-1;
  for imageSrcRow:=0 to imageSrcRowmax do begin
    pDest:=imageDest.Bits;
    Inc(pDest,imageSrcRow*w);
    for imageSrcCol:=0 to imageSrcColmax do begin
      n:=0;
      // nucleus
      valueNucleus:=pSrc[imageSrcRow*w+imageSrcCol];
      tIntensityNucleus:=valueNucleus and $000000FF; // Intensity(value)
      // calculate n = the USAN : Univalue Segment Assimilating Nucleus
      imageSrcNeighborRowmin:=imageSrcRow-roundImageMaskHeightDiv2; if imageSrcNeighborRowmin<0 then imageSrcNeighborRowmin:=0;
      imageSrcNeighborColmin:=imageSrcCol-roundImageMaskWidthDiv2; if imageSrcNeighborColmin<0 then imageSrcNeighborColmin:=0;
      imageSrcNeighborRowmax:=imageSrcRow+roundImageMaskHeightDiv2; if imageSrcNeighborRowmax>imageSrcRowmax then imageSrcNeighborRowmax:=imageSrcRowmax;
      imageSrcNeighborColmax:=imageSrcCol+roundImageMaskWidthDiv2; if imageSrcNeighborColmax>imageSrcColmax then imageSrcNeighborColmax:=imageSrcColmax;
      for imageSrcNeighborRow:=imageSrcNeighborRowmin to imageSrcNeighborRowmax do begin
        for imageSrcNeighborCol:=imageSrcNeighborColmin to imageSrcNeighborColmax do begin
          colInMask:=roundImageMaskWidthDiv2+imageSrcNeighborCol-imageSrcCol;
          rowInMask:=roundImageMaskHeightDiv2+imageSrcNeighborRow-imageSrcRow;
          colorMask:=Image.getPixel(mask,colInMask,rowInMask);
          if colorMask=clWhite32 then begin
            value:=pSrc[imageSrcNeighborRow*w+imageSrcNeighborCol];
            tIntensity:=value and $000000FF;  // Intensity(value)
            diff:=abs(tIntensity-tIntensityNucleus);
            if diff < t then Inc(n);
          end;
        end;
      end;
      // edge response : if n is small, then it's probably an edge
      // (the smaller the USAN area, the larger the edge response)
      if n<g then begin
        isDetected := True ;
        // if we are looking for edges, Test for false positives by finding the USAN's centroid
        if testCorner Then Begin
          n2 := 0 ;
          centroidRow := 0 ;
          centroidCol := 0 ;
          sumRow := 0 ;
          sumCol := 0 ;
          for imageSrcNeighborRow:=imageSrcNeighborRowmin to imageSrcNeighborRowmax do begin
            for imageSrcNeighborCol:=imageSrcNeighborColmin to imageSrcNeighborColmax do begin
              value:=pSrc[imageSrcNeighborRow*w+imageSrcNeighborCol];
              tIntensity:=value and $000000FF;  // Intensity(value)
              diff:=abs(tIntensity-tIntensityNucleus);
              if diff < t then Begin
                inc(n2) ;
                inc(sumRow, imageSrcNeighborRow) ;
                inc(sumCol, imageSrcNeighborCol) ;
              end ;
            end;
          end;
          if n2 > 1 Then Begin
            centroidRow := (sumRow / n2) - imageSrcRow ;
            centroidCol := (sumCol / n2) - imageSrcCol ;
            if (abs(centroidRow) + abs(centroidCol)) <= t2 then
              isDetected := false ;
          end ;
        end ;
        // at last, if pixel is interesting
        if isDetected then begin
          value:=floor((g-n)*gTo255);
        end else
          value:=0;
      end else begin
        value:=0;
      end;
      // set gray value with Color32 for optimization
      pDest^[0]:= Color32(value,value,value);
      Inc(pDest);
    end;
  end;
  // final stage for corner
  if parameterMode.Value='CORNER' then begin
    filterNonMaximaSuppression.setParameterImage('inImage',imageDest);
    filterNonMaximaSuppression.setParameterImage('outImage',parameterImageOut.Image);
    filterNonMaximaSuppression.Run;
  end;
end;


end.