unit filterCorrelation;
(* ***** 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 ***** *)

{
 edurand (filters@edurand.com)
}

interface
uses
  filter, fparameters, image;

type
  TFilterCorrelation = class(TFilter)
  public
    constructor Create; override;
    destructor Destroy; override;
    procedure Run(); override;
  private
    parameterImageIn1, parameterImageIn2 : TParameterImage;
    ouputParameterCorrelation : TParameterSingle;
    parameterRangeMin, parameterRangeMax : TParameterInteger;
    parameterMonitoring : TParameterBoolean;
    parameterImageOutMonitoring : TParameterImage;
    image1, image2 : PBitmap32;
    correlation : Single;
    imgOutMonitoring : PBitmap32;
    parameterImage2CenterXinImage1, parameterImage2CenterYinImage1 : TParameterInteger;
    procedure _run();
  end;

implementation
uses
  Math, imageIO;

constructor TFilterCorrelation.Create;
begin
  inherited;
  parameterImageIn1:=addParameterImage('inImage1','input image 1 : the image to scan');
  parameterImageIn2:=addParameterImage('inImage2','input image 2 : the image of the pattern to test in image 1');
  ouputParameterCorrelation := addOutputParameterSingle('correlation', 'correlation between same size image 1 and 2');
  parameterRangeMin:=addParameterInteger('rangeMin','calcul correlation on pixel of intensity is in the range',0,255,0);
  parameterRangeMax:=addParameterInteger('rangeMax','calcul correlation on pixel of intensity is in the range',0,255,255);
  parameterMonitoring:=addParameterBoolean('monitoring','monitoring',True);
  parameterImageOutMonitoring:=addOutputParameterImage('outImageMonitoring', 'image out monitoring');
  parameterImage2CenterXinImage1:=addParameterInteger('X','the x center of image 2, in image 1',-maxint+1,maxint,maxint);
  parameterImage2CenterYinImage1:=addParameterInteger('Y','the y center of image 2, in image 1',-maxint+1,maxint,maxint);
end;

destructor TFilterCorrelation.Destroy;
begin
  inherited;
  image.freeImage(imgOutMonitoring);
end;

procedure TFilterCorrelation.Run();
begin
  image1:=parameterImageIn1.Image;
  image2:=parameterImageIn2.Image;
  if (image1<>nil) and (image2<>nil) then begin
    _run();
  end;
end;

procedure TFilterCorrelation._run();
var
  center2In1 : TFPoint;
  x,xMin,xMax,y,yMin,yMax : Integer;
  x2,x2Min,y2,y2Min : Integer;
  image2WidthDiv2, image2HeightDiv2 : Integer;
  rangeMin, rangeMax : Integer;
  p1p1,p2p2,p1p2 : Extended;
  intensity1, intensity2 : Integer;
  monitoring : boolean;
  pSrc1, pSrc2, pImgOutMonitoring : PColor32Array;
  tmpExtended : Extended;
begin
  // pixel where intensity is outside this range will be ignored
  rangeMin:=parameterRangeMin.Value;
  rangeMax:=parameterRangeMax.Value;
  // image 2 is positionned on image 1 by this parameters
  if parameterImage2CenterXinImage1.Value=maxint then center2In1.x:=image1.Width div 2 else center2In1.x:=parameterImage2CenterXinImage1.Value;
  if parameterImage2CenterYinImage1.Value=maxint then center2In1.y:=image1.Height div 2 else center2In1.y:=parameterImage2CenterYinImage1.Value;
  // calcul the part of image 1 to scan
  image2WidthDiv2:=image2.Width div 2;
  image2HeightDiv2:=image2.Height div 2;
  xMin:=Floor(center2In1.x-image2WidthDiv2);
  if xMin<0 then begin
    x2Min:=-xMin;
    xMin:=0;
  end else begin
    x2Min:=0;
  end;
  xMax:=Floor(center2In1.x+image2WidthDiv2);
  if xMax>=image1.Width then begin
    xMax:=image1.Width-1;
  end;
  yMin:=Floor(center2In1.y-image2HeightDiv2);
  if yMin<0 then begin
    y2Min:=-yMin;
    yMin:=0;
  end else begin
    y2Min:=0;
  end;
  yMax:=Floor(center2In1.y+image2HeightDiv2);
  if yMax>=image1.Height then begin
    yMax:=image1.Height-1;
  end;
  // monitoring ?
  monitoring:=parameterMonitoring.Value;
  if monitoring=true then begin
    imgOutMonitoring:=image.eraseOrCreateImageLike(imgOutMonitoring,parameterImageIn1.Image);
    image.copyImageToImage(parameterImageIn1.Image,imgOutMonitoring);
  end;
  // start the scan
  p1p1:=0;
  p2p2:=0;
  p1p2:=0;
  y2:=y2Min;
  for y:=yMin to yMax do begin
    if y2>=image2.Height then begin
      break;
    end;
    x2:=x2Min;
    pSrc1:=image1.Bits;
    Inc(pSrc1,y*image1.Width);
    Inc(pSrc1,xMin);
    if monitoring=true then begin
      pImgOutMonitoring:=imgOutMonitoring.Bits;
      Inc(pImgOutMonitoring,y*imgOutMonitoring.Width);
      Inc(pImgOutMonitoring,xMin);
    end;
    pSrc2:=image2.Bits;
    Inc(pSrc2,y2*image2.Width);
    Inc(pSrc2,x2Min);
    for x:=xMin to xMax do begin
      if x2>=image2.Width then begin
        break;
      end;
      //intensity1:=image.Intensity(image.getPixel(image1,x,y));
      //intensity2:=image.Intensity(image.getPixel(image2,x2,y2));
      intensity1:=image.Intensity(pSrc1^[0]);
      intensity2:=image.Intensity(pSrc2^[0]);
      if (intensity1>=rangeMin) and (intensity1<=rangeMax) and
       (intensity2>=rangeMin) and (intensity2<=rangeMax) then begin
        p1p1:=p1p1 + intensity1*intensity1;
        p2p2:=p2p2 + intensity2*intensity2;
        p1p2:=p1p2 + intensity1*intensity2;
        if monitoring=true then begin
          //image.setPixel(imgOutMonitoring,x,y,image.Gray32(intensity2));
          pImgOutMonitoring^[0]:=image.Gray32(intensity2);
        end;
      end;
      Inc(x2);
      Inc(pSrc1);
      Inc(pSrc2);
      if monitoring=true then begin
        Inc(pImgOutMonitoring);
      end;
    end;
    Inc(y2);
  end;
  tmpExtended := p1p1 * p2p2;
  if tmpExtended>0 then begin
    correlation := p1p2 / Sqrt( tmpExtended );
  end else begin
    correlation := 1;
  end;
  // set out put parameters
  setOutputParameterSingle('correlation', correlation);
  if monitoring=true then begin
    setOutputParameterImage('outImageMonitoring',imgOutMonitoring);
  end;
end;

{procedure TFilterCorrelation.Run();
var
  pSrc1, pSrc2, pImgOutMonitoring : PColor32Array;
  x, xMin, xMax : Integer;
  intensity1, intensity2 : Integer;
  p1p1,p2p2,p1p2 : Extended;
  rangeMin, rangeMax : Integer;
  monitoring : boolean;
  center1, center2In1 : TFPoint;
  y, yMin, yMax : Integer;
begin
  monitoring:=parameterMonitoring.Value;
  if monitoring=true then begin
    imgOutMonitoring:=imageIO.eraseOrCreateImageLike(imgOutMonitoring,parameterImageIn1.Image);
    imageIO.copyImageToImage(parameterImageIn1.Image,imgOutMonitoring);
  end;
  // pixel where intensity is outside this range will be ignored
  rangeMin:=parameterRangeMin.Value;
  rangeMax:=parameterRangeMax.Value;
  // the center of image 1
  center1.x:=parameterImageIn1.Image.Width div 2;
  center1.y:=parameterImageIn1.Image.Height div 2;
  // the center of image 2
  center2In1:=center1;
  // the xMin->xMax x to scan in image 1
  xMin:=Max(0,Floor(center2In1.x-parameterImageIn2.Image.Width div 2));
  xMax:=Min(parameterImageIn1.Image.Width-1,xMin+parameterImageIn2.Image.Width-1);
  // the yMin->yMax y to scan in image 1
  yMin:=Max(0,Floor(center2In1.y-parameterImageIn2.Image.Height div 2));
  yMax:=Min(parameterImageIn1.Image.Height-1,yMin+parameterImageIn2.Image.Height-1);
  // scan 
  p1p1:=0;
  p2p2:=0;
  p1p2:=0;
  correlation:=0;
  for y:=yMin to yMax do begin
    pSrc1:=parameterImageIn1.Image.Bits;
    Inc(pSrc1,y*parameterImageIn1.Image.Width);
    Inc(pSrc1,xMin);
    pImgOutMonitoring:=imgOutMonitoring.Bits;
    Inc(pImgOutMonitoring,y*imgOutMonitoring.Width);
    Inc(pImgOutMonitoring,xMin);
    pSrc2:=parameterImageIn2.Image.Bits;
    Inc(pSrc2,(y-yMin)*parameterImageIn2.Image.Width);
    for x:=xMin to xMax do begin
      intensity1:=image.Intensity(pSrc1^[0]);
      intensity2:=image.Intensity(pSrc2^[0]);
      if (intensity1>=rangeMin) and (intensity1<=rangeMax) and
       (intensity2>=rangeMin) and (intensity2<=rangeMax) then begin
        p1p1:=p1p1 + intensity1*intensity1;
        p2p2:=p2p2 + intensity2*intensity2;
        p1p2:=p1p2 + intensity1*intensity2;
        if monitoring=true then begin
          pImgOutMonitoring^[0]:=image.Gray32(intensity2);
        end;
      end;
      Inc(pSrc1);
      Inc(pSrc2);
      Inc(pImgOutMonitoring);
    end;
  end;
  correlation:= p1p2 / Sqrt( p1p1 * p2p2 );
  setOutputParameterSingle('correlation', correlation);
  if monitoring=true then begin
    setOutputParameterImage('outImageMonitoring',imgOutMonitoring);
  end;
end;}

end.