unit FFT;

interface

uses Windows, Graphics, Forms, extctrls;

type Tra = record // record pro urceni posuvu
             X       : single; {posun v ose X       }
             Y       : single; {posun v ose y       }
             Quality : single; {spolehlivost urceni }
           end {record};
     MagRot = record  // record pro urceni rotace a zvetseni
                Magnification : single; {zvetseni            }
                Rotation      : single; {uhel otoceni        }
                Quality       : single; {spolehlivost urceni }
              end{record};
     TTransfType = (ttTranslation,ttMagRot,ttMagRotTransl,ttNone);

{funkce, ktera urcuje vzajemne otoceni a zmenu meritka}
function FindImgMagnRot(Img1Pointer      : pointer; // pointer na 1. obraz
                        Img1Width        : integer; // sirka 1. obrazu v pixelech
                        Img1Height       : integer; // vyska 1. obrazu v pixelech
                        Img1Channels     : integer; // kolik slozek ma pixel 1. obrazu (1 nebo 3)
                        Img1DataType     : string;  // datovy typ 1. obrazu ('BYTE' nebo 'WORD' nebo 'SMALLINT' nebo 'SINGLE' nebo 'DOUBLE');
                        Img1RoiRect      : TRect;   // obdelnik v 1. obraze pouzity k urceni transformace
                        Img2Pointer      : pointer; // pointer na 2. obraz
                        Img2Width        : integer; // sirka 2. obrazu v pixelech
                        Img2Height       : integer; // vyska 2. obrazu v pixelech
                        Img2Channels     : integer; // kolik slozek ma pixel 2. obrazu (1 nebo 3)
                        Img2DataType     : string;  // datovy typ 2. obrazu ('BYTE' nebo 'WORD' nebo 'SMALLINT' nebo 'SINGLE' nebo 'DOUBLE');
                        Img2RoiRect      : TRect;   // obdelnik v 2. obraze pouzity k urceni transformace
                        LowPassFilter    : single;  // parametr LP filtru
                        SubEpsilon       : single;  // polomer okoli pro sub-pixelove hledani maxima
                        CircleRoi        : boolean; // false ... square ROI, true ... circle ROI
                        WindEdge         : single;  // parametr pro prubeh Hanningova okna (0 .. ostry okraj, 1 .. plynuly pokles od stredu okna)
                        OnlyInvFT        : boolean; // false ... probehne cela procedura, true ... zacne se LP filtrem (pouziva se pokud uz
                                                    // procedura jednou probehla a byl zmenen jen LP filtr)
                        OnlySubPixelPeak : boolean  // true ... hleda se znovu jen jen sub-pixelove maximum
                                        ): MagRot;  // vysledek - zvetseni a otoceni
// vysledek fazove korealce (peak) lze zobrazit pomoci procedury DisplayPhaseCorrResult                                   

{funkce, ktera urcuje vzajemny posun dvou obrazu}
function FindImgTranslation(Img1Pointer      : pointer;  // pointer na 1. obraz
                            Img1Width        : integer;  // sirka 1. obrazu v pixelech
                            Img1Height       : integer;  // vyska 1. obrazu v pixelech
                            Img1Channels     : integer;  // kolik slozek ma pixel 1. obrazu (1 nebo 3)
                            Img1DataType     : string;   // datovy typ 1. obrazu ('BYTE' nebo 'WORD' nebo 'SMALLINT' nebo 'SINGLE' nebo 'DOUBLE');
                            Img1RoiRect      : TRect;    // obdelnik v 1. obraze pouzity k urceni posuvu
                            Img2Pointer      :  pointer; // pointer na 2. obraz
                            Img2Width        : integer;  // sirka 2. obrazu v pixelech
                            Img2Height       : integer;  // vyska 2. obrazu v pixelech
                            Img2Channels     : integer;  // kolik slozek ma pixel 2. obrazu (1 nebo 3)
                            Img2DataType     : string;   // datovy typ 2. obrazu ('BYTE' nebo 'WORD' nebo 'SMALLINT' nebo 'SINGLE' nebo 'DOUBLE');
                            Img2RoiRect      : TRect;    // obdelnik v 2. obraze pouzity k urceni posuvu
                            LowPassFilter    : single;   // parametr LP filtru
                            SubEpsilon       : single;   // polomer okoli pro sub-pixelove hledani maxima
                            CircleRoi        : boolean;  // false ... square ROI, true ... circle ROI
                            WindEdge         : single;   // parametr pro prubeh Hanningova okna (0 .. ostry okraj, 1 .. plynuly pokles od stredu okna)
                            OnlyInvFT        : boolean;  // false ... probehne cela procedura, true ... zacne se LP filtrem (pouziva se pokud uz
                                                         // procedura jednou probehla a byl zmenen jen LP filtr)
                            OnlySubPixelPeak : boolean   // true ... hleda se znovu jen jen sub-pixelove maximum
                                           )  : Tra;     // vysledkem je record s posuvem
// vysledek fazove korealce (peak) lze zobrazit pomoci procedury DisplayPhaseCorrResult


{procedure, ktera zobrazuje vysledek fazove korealce (peak)}
procedure DisplayPhaseCorrResult(Image     :  TImage;     // Bitmapa (pf24bit) pro zobrazeni vysledku fazove korealce
                                 ScrollBox : TScrollBox;  // procedura sama scrolluje na peak, pokud je Bitmapa ve ScrollBoxu,
                                                          // pokud neni, musi byt ScrollBox=nil
                                 Max, Min  : single;
                                 Black     : single;      // 0<=Black<White, Black=0 znamena, ze min hodnota je cerna
                                 White     : single);     // Black<White<=1, <White=1 znamena, ze max hodnota (peak) je bila


implementation

uses SysUtils, MainFormPhaseCorr, Global;

const N               = MaxROISize;
      NN              = N*N;
      Ndiv2           = N div 2;
      MaxImgSize      = 6000*6000; // pro obrazy, ze kterych se bere vyrez N*N
      SubSample       = 10;        // pro sub-pixelove hledani peaku
      NumSubPixelIter = 100;       // pocet iteraci pri hledaní teziste peaku

type  {typy pro FFT a fazovou korelaci}
      Complex = record
                 Re : single;
                 Im : single;
               end;{record}

     Spectrum = array [1..N,1..N] of complex;
     FFTImage = array [1..N,1..N] of single;

     Coordinate = record
                    X : integer;
                    Y : integer;
                  end {record};

     CoordinateS = record
                    X : single;
                    Y : single;
                   end {record};

     ImgType         = (itImage1, itImage2, itSpectrum1, itSpectrum2,
                        itTranslation, itMagRot, itNone);
     TDisplay        = (dRGB,dG,dRB);
     pSpc     = ^Spectrum;

     {typy poli pro nacteni vyrezu obrazu}
     TByteArr        = array[0..MaxImgSize] of byte;
     PByteArr        = ^TByteArr;
     TWordArr        = array[0..MaxImgSize] of word;
     PWordArr        = ^TWordArr;
     TSmallintArr    = array[0..MaxImgSize] of smallint;
     PSmallintArr    = ^TSmallintArr;
     TSingleArr      = array[0..MaxImgSize] of Single;
     PSingleArr      = ^TSingleArr;
     TDoubleArr      = array[0..MaxImgSize] of Double;
     PDoubleArr      = ^TDoubleArr;
     TByteClr        = record
                        R,G,B : Byte;
                       end;{record}
     TByteClrArr     = array[0..MaxImgSize] of TByteClr;
     PByteClrArr     = ^TByteClrArr;
     TWordClr        = record
                        R,G,B : word;
                       end;{record}
     TWordClrArr     = array[0..MaxImgSize] of TWordClr;
     PWordClrArr     = ^TWordClrArr;
     TSmallintClr    = record
                         R,G,B : Smallint;
                       end;{record}
     TSmallintClrArr = array[0..MaxImgSize] of TSmallintClr;
     PSmallintClrArr = ^TSmallintClrArr;
     TSingleClr      = record
                        R,G,B : Single;
                       end;{record}
     TSingleClrArr   = array[0..MaxImgSize] of TSingleClr;
     PSingleClrArr   = ^TSingleClrArr;
     TDoubleClr      = record
                      R,G,B : Double;
                     end;{record}
     TDoubleClrArr   = array[0..MaxImgSize] of TDoubleClr;
     PDoubleClrArr   = ^TDoubleClrArr;


var
    LPfilter                   : single;
    EdgeHanning                : single;
    AutoTransl                 : Tra;
    Sp1, Sp2                   : Spectrum;
    Img1, Img2, Pom, Res       : FFTImage;
    PeakCoordinate             : Coordinate;
    SubEpsilon                 : single;       // polomer okoli pro vypocet teziste peaku

function  fftw2d_create_plan( nx, ny, d, f : integer ) : pointer;
          cdecl external 'sfftw.dll';
procedure fftwnd_one( p, i, o : pointer );
          cdecl external 'sfftw.dll';

var  b : Spectrum;

const FFTW_FORWARD  = -1;
      FFTW_BACKWARD = 1;
      FFTW_IN_PLACE = 8;
      FFTW_ESTIMATE = 0;
      FFTW_MEASURE  = 1;

procedure DisplayFFTImage(var Img : FFTImage;
                          Image   : TImage); //pf24bit
var i, j, ii, jj      : integer;
    ScrollX, ScrollY  : integer;
    NormPix, k, Range : single;
    SetMax, SetMin    : single;
    b                 : byte;
    Buff              : array [1 .. 3*N] of byte; // buffer pro kresleni obrazu
begin
  SetMin:=Img[1,1];
  SetMax:=Img[1,1];
  for i:=1 to N do
    for j:=1 to N do
      begin
        if Img[i,j]>SetMax then SetMax:=Img[i,j];
        if Img[i,j]<SetMin then SetMin:=Img[i,j];
      end;
  Range:=SetMax-SetMin;
  if Range<0.0001 then Range:=0001;
  Image.Picture.Bitmap.Width:=N;
  Image.Picture.Bitmap.Height:=N;
  Image.Width:=N;
  Image.Height:=N;
  Image.Picture.Bitmap.PixelFormat:=pf24bit;
  Image.Picture.Bitmap.PixelFormat:=pf24bit;

  for j:=1 to N do
    begin
     for i:=1 to N do
       begin
         NormPix:=(Img[i,j]-SetMin)/Range;
         if NormPix>1 then NormPix:=1;
         if NormPix<0 then NormPix:=0;
         b:=round(255.0*NormPix);
         Buff[3*i]:=b;
         Buff[3*i-1]:=b;
         Buff[3*i-2]:=b;
       end;
       Move(Buff,Image.Picture.Bitmap.ScanLine[j-1]^,SizeOf(Buff));
     end;
  Image.Visible:=true;
  Image.Repaint;
end;

function Permute(X : integer) : integer;
{permutace Fourier. spektra}
begin
 Permute:=1 + (X-1+Ndiv2) mod N;
end;

//function HiPass(i,j : integer) : single;
//{ Filter HiPass aplikovatelny na slozku [i,j] Fourier spektra }
//{ (permutovaneho) tj. nulova frekvence je v [N/2,N/2].        }
//{ dle clanku - nepouziva se                                   }
//var x, y, F : single;
//begin
//  x:=i/N - 0.5;
//  y:=j/N - 0.5;
//  F:=cos(pi*x)*cos(pi*y);
//  HiPass:=(1.0-F)*(2.0-F);
//end;



procedure LowPass2(var OldSpc, NewSpc : Spectrum; sigma : single);
{ LowPass filter pro Fourier. spektrum  (pernutovane)  }
{ tj. nulova frekvence je v [N/2,N/2].                 }
{ Sigma urcuje prorustnost                             }
var i, j, ii, jj    : integer;
    x, y, d, k, f1, f2, size : single;
begin
  f2:=0.35;
  f1:=0.01;
  size:=f2-f1;
  for i:=1 to N do
    for j:=1 to N do
      begin
         ii:=permute(i);
         jj:=permute(j);
         x:=(ii-0.5);
         y:=(jj-0.5);
         d:=2*sqrt(sqr(x-N/2)+sqr(y-N/2))/N;
       if d>f2 then k:=0 else if d>f1 then k:=0.5*(cos(pi*(d-f1)/size)+1) else k:=1;
       NewSpc[i,j].Re:=k*OldSpc[i,j].Re;
       NewSpc[i,j].Im:=k*OldSpc[i,j].Im;
       Res[ii,jj]:=k;
      end;
end;



procedure LowPass(var OldSpc, NewSpc : Spectrum; sigma : single);
{ LowPass filter pro Fourier. spektrum  (pernutovane)  }
{ (permutovane) tj. nulova frekvence je v [N/2,N/2].   }
{ Sigma urcuje prorustnost                             }
var i, j, ii, jj    : integer;
    x, y, r, k, sig : single;

begin
  if sigma>0 then sig:=1/sigma else sig:=1.0;
  for i:=1 to N do
     for j:=1 to N do
       begin
         ii:=permute(i);
         jj:=permute(j);
         x:=(ii-0.5)/N-0.5;
         y:=(jj-0.5)/N-0.5;
         r:=-1.0*(sqr(x)+sqr(y))/sig;
         if r<-30 then k:=0 else k:=exp(r);
         NewSpc[i,j].Re:=k*OldSpc[i,j].Re;
         NewSpc[i,j].Im:=k*OldSpc[i,j].Im;
         Res[ii,jj]:=k;
       end;
end;

function BiLiInterpolationImg(var x, y : single; var Img : FFTImage) : single;
{Bilinearni interpolace v obraze}
var i1, i2, j1, j2 : integer;
    f1, f2, k1, k2 : single;

begin
  if (x<1)  or (x>N) or (y<1) or (y>N) then
    begin
      BiLiInterpolationImg:=0.0;
      Exit;
    end;
  i1:=trunc(x);
  if i1<N then i2:=i1+1 else i2:=i1;
  j1:=trunc(y);
  if j1<N then j2:=j1+1 else j2:=j1;
  k1:=x-i1;
  k2:=1.0-k1;
  f1:=k2*Img[i1,j1] + k1*Img[i2,j1];
  f2:=k2*Img[i1,j2] + k1*Img[i2,j2];
  k1:=y-j1;
  k2:=1.0-k1;
  BiLiInterpolationImg:=k2*f1 + k1*f2;
end;


function BiLiInterpolationAmpSpc(var x, y : single; var Img : FFTImage) : single;
{Bilinearni interpolace v amplitudovem spektru }
var i1, i2, j1, j2 : integer;
    f1, f2, k1, k2 : single;

begin
  i1:=trunc(x);
  j1:=trunc(y);
  k1:=x-i1;
  k2:=1.0-k1;
  i2:=i1+1;
  j2:=j1+1;
  if i1=0 then i1:=N;
  if j1=0 then j1:=N;
  if i2=N+1 then i2:=1; {vyuziti periodicity}
  if j2=N+1 then j2:=1; {vyuziti periodicity}
  f1:=k2*Img[i1,j1] + k1*Img[i2,j1];
  f2:=k2*Img[i1,j2] + k1*Img[i2,j2];
  k1:=y-j1;
  k2:=1.0-k1;
  BiLiInterpolationAmpSpc:=k2*f1 + k1*f2;
end;

procedure PolarTransf(var OldImg, NewImg : FFTImage);
{ Transformuje obraz amplitudoveho spektra do polarni}
{ soustavy rouradne}
var x, y          : single;
    i, j          : integer;
    R, Fi         : single;
    StepFi, StepR : single;
begin
  StepFi:=pi/N;
  StepR:=0.5;
  for i:=1 to N do
    for j:=1 to N do
      begin
         Fi:=StepFi*i;
         R:=StepR*j;
         x:=R*cos(Fi)+n/2;
         y:=R*sin(Fi)+n/2;
         NewImg[i,j]:=BiLiInterpolationAmpSpc(x,y,OldImg);
      end;
end;

procedure PolarTransfZoom(var OldImg, NewImg : FFTImage);
{ Transformuje obraz amplitudoveho spektra do polarni}
{ soustavy rouradne}
var x, y          : single;
    i, j          : integer;
    R, Fi         : single;
    StepFi, StepR : single;
begin
  StepFi:=pi/N;
  StepR:=1/2;
  for i:=1 to N do
    for j:=1 to N do
      begin
         Fi:=StepFi*i;
         R:=StepR*j;
         x:=R*cos(Fi)+n/2;
         y:=R*sin(Fi)+n/2;
         NewImg[i,j]:=BiLiInterpolationAmpSpc(x,y,OldImg);
      end;
end;


procedure LogPolarTransf(var OldImg, NewImg : FFTImage);
{ Transformuje obraz amplitudoveho spektra do Log-polarni s.s.}
{ tj. na ose R je Log meritko, uhel fi je normalni            }
var x, y          : single;
    i, j          : integer;
    R, Fi         : single;
    StepFi, StepR : single;
    a             : single;
begin
  a:=exp((ln(N)-ln(2))/N);
  StepFi:=pi/N;
  StepR:=1.0;
  for i:=1 to N do
    for j:=1 to N do
      begin
         Fi:=StepFi*i;
         R:=StepR*j;
         R:=exp(R*ln(a));
         x:=R*cos(Fi)+n/2;
         y:=R*sin(Fi)+n/2;
         NewImg[i,j]:=BiLiInterpolationAmpSpc(x,y,OldImg);
      end;
end;

procedure LogPolarTransf4(var OldImg, NewImg : FFTImage);
{ Transformuje obraz amplitudoveho spektra do Log-polarni s.s.}
{ tj. na ose R je Log meritko, uhel fi je normalni            }
{ Pouziva 4 nasobny oversample                                }
var x, y          : single;
    sum           : single;
    i, j, k, l    : integer;
    R, Fi         : single;
    StepFi, StepR : single;
    a             : single;
begin
  a:=exp((ln(N)-ln(2))/N);
  StepFi:=0.25*(pi/N);
  StepR:=0.25;
  for i:=0 to N-1 do
    for j:=0 to N-1 do
      begin
        sum:=0.0;
        for k:=0 to 3 do
          for l:=0 to 3 do
            begin
              Fi:=StepFi*(4*i+k);
              R:=StepR*(4*j+l);
              R:=exp(R*ln(a));
              x:=1+R*cos(Fi)+n/2;
              y:=1+R*sin(Fi)+n/2;
              sum:=sum+BiLiInterpolationAmpSpc(x,y,OldImg);
            end;
        NewImg[i+1,j+1]:=sum/16.0;
      end;
end;

procedure LogPolarTransfNew(var OldImg, NewImg : FFTImage);
{ Transformuje obraz amplitudoveho spektra do Log-polarni s.s.}
{ tj. na ose R je Log meritko, uhel fi je normalni            }
{ Pouziva 4 nasobny oversample                                }
var x, y          : single;
    sum           : single;
    wei           : integer;
    i, j, k, l    : integer;
    R1, R2, r, Fi : single;
    StepFi, StepR : single;
    a             : single;
begin
  a:=exp((ln(N)-ln(2))/N);
  StepFi:=pi/N; 
  StepR:=1;
  for i:=0 to N-1 do
    begin
      R1:=0;
      Fi:=StepFi*i;
      for j:=1 to N do
        begin
          sum:=0.0;
          wei:=0;
          R2:=StepR*j;
          R2:=exp(R2*ln(a));
          R:=R1;
          repeat
            inc(wei);
            x:=1+R*cos(Fi)+n/2;
            y:=1+R*sin(Fi)+n/2;
            sum:=sum+BiLiInterpolationAmpSpc(x,y,OldImg);
            R:=R+1;
           until R>R2;
          R1:=R2;
          NewImg[i+1,j]:=sum/wei;
        end;
    end;
end;



function arg(re,im: single) : single;
{ pocita argument komplexniho cisla }
var fi : single;
begin
  if (abs(re)<0.000001) and (abs(im)<0.000001) then fi:=0 else
   begin
     if (re=0) then if (im>0) then fi:=90 else fi:=270;
     if (im=0) then if (re>0) then fi:= 0 else fi:=180;
     if (re>0) and (im>0) then fi:=180*arctan(im/re)/pi;
     if (re>0) and (im<0) then fi:=180*arctan(im/re)/pi+360;
     if (re<0) and (im>0) then fi:=180*arctan(im/re)/pi+180;
     if (re<0) and (im<0) then fi:=180*arctan(im/re)/pi+180;
   end;
  arg:=fi;
end;


procedure FourierTransform(var Img : FFTImage; var Spc : Spectrum);
var sign : longint;
    i, j : integer;
    plan : pointer;
    pSpc : pointer;
begin

  sign:=1;
  for  i:=1 to N do
    for j:=1 to N do
      begin
        Spc[i,j].Re:=Img[i,j];
        Spc[i,j].Im:=0.0;
      end;
  pSpc := @Spc;
  plan := fftw2d_create_plan( N, N, FFTW_FORWARD, FFTW_ESTIMATE or FFTW_IN_PLACE );
  fftwnd_one( plan, pSpc, nil );
end;



procedure InverseFourierTransform(var Spc : Spectrum; var Img : FFTImage);
var sign : longint;
    i, j : integer;
    plan : pointer;
    pSpc : pointer;
begin
  b:=Spc;
  pSpc := @b;
  plan := fftw2d_create_plan( N, N, FFTW_BACKWARD, FFTW_ESTIMATE or FFTW_IN_PLACE );
  fftwnd_one( plan, pSpc, nil );
  for  i:=1 to N do
    for j:=1 to N do Img[i,j]:=b[i,j].Re;
end;



procedure AmplitudeSpectrum(var Spec : Spectrum; var Amp : FFTImage);
var i, j : Integer;
{ Spec ... vstupni spektrum              }
{ Amp .... vystupni amplitudove spekrtum }
begin
  for i:=1 to N do
    for j:=1 to N do
      Amp[i,j]:=sqrt(sqr(Spec[i,j].Re)+sqr(Spec[i,j].Im));
end;

procedure AmplitudeSpectrumPerm(var Spec : Spectrum; var Amp : FFTImage);
var i, j, ii, jj, m : integer;
{ Spec ... vstupni spektrum (nepermutovane tj. ma nulovou  }
{          frekvenci ma v prvku [0,0]                      }
{ Amp .... vystupni amplitudove spekrtum permutovane       }
{          tj. ma nulovou frekvenci ma v prvku [N/2,N/2]   }
begin
  for i:=1 to N do
    for j:=1 to N do
      begin
        ii:=Permute(i);
        jj:=Permute(j);
        Amp[i,j]:=sqrt(sqr(Spec[ii,jj].Re)+sqr(Spec[ii,jj].Im));
      end;
end;


procedure LogAmplitudeSpectrum(var Spec : Spectrum; var Amp : FFTImage);
{ Spec ... vstupni spektrum                  }
{ Amp .... vystupni Log amplitudove spekrtum }
var i, j : integer;
begin
  for i:=1 to N do
    for j:=1 to N do
     begin
      Amp[i,j]:=sqrt(sqr(Spec[i,j].Re)+sqr(Spec[i,j].Im));
      if  Amp[i,j]>=0.0000001 then Amp[i,j]:=ln(Amp[i,j]);
     end;
end;


procedure LogAmplitudeSpectrumPerm(var Spec : Spectrum; var Amp : FFTImage);
{ Spec ... vstupni spektrum (nepermutovane tj. ma nulovou  }
{          frekvenci ma v prvku [0,0] }
{ Amp .... vystupni Log amplitudove spekrtum permutovane   }
{          tj. ma nulovou frekvenci ma v prvku [N/2,N/2]   }
var i, j, M   : integer;
    ii, jj    : integer;
begin
  for i:=1 to N do
    for j:=1 to N do
     begin
      ii:=Permute(i);
      jj:=Permute(j);
      Amp[i,j]:=sqrt(sqr(Spec[ii,jj].Re)+sqr(Spec[ii,jj].Im));
      if  Amp[i,j]>=0.0000001 then Amp[i,j]:=ln(Amp[i,j]);
     end;
end;

procedure LogAmplitudeSpectrumPermModif(var Spec : Spectrum; var Amp : FFTImage);
{ Spec ... vstupni spektrum (nepermutovane tj. ma nulovou  }
{          frekvenci ma v prvku [0,0] }
{ Amp .... vystupni Log amplitudove spekrtum permutovane   }
{          tj. ma nulovou frekvenci ma v prvku [N/2,N/2]   }
var i, j, M   : integer;
    ii, jj    : integer;
begin
  for i:=1 to N do
    for j:=1 to N do
     begin
      ii:=Permute(i);
      jj:=Permute(j);
      Amp[i,j]:=sqrt(sqr(Spec[ii,jj].Re)+sqr(Spec[ii,jj].Im))+1;
      Amp[i,j]:=ln(Amp[i,j])+1;
      Amp[i,j]:=ln(Amp[i,j])+1;
      //Amp[i,j]:=ln(Amp[i,j])+1;
     end;
   //Amp[1,1]:=0.0;
end;

procedure LogAmplitudeSpectrumPermHi(var Spec : Spectrum; var Amp : FFTImage);
{ Spec ... vstupni spektrum (nepermutovane tj. ma nulovou  }
{          frekvenci ma v prvku [0,0] }
{ Amp .... vystupni Log amplitudove spektrum permutovane   }
{          tj. ma nulovou frekvenci ma v prvku [N/2,N/2].  }
{          Na spektrum je aplikovan High Pass filter.      }
var i, j, M   : integer;
    ii, jj    : integer;
begin
  for i:=1 to N do
    for j:=1 to N do
     begin
      ii:=Permute(i);
      jj:=Permute(j);
      Amp[i,j]:=sqrt(sqr(Spec[ii,jj].Re)+sqr(Spec[ii,jj].Im));
      if  Amp[i,j]>=0.0000001 then Amp[i,j]:=ln(Amp[i,j]);
     end;
end;

procedure PhaseSpectrum(var Spec : Spectrum; var Amp : FFTImage);
{ Spec ... vstupni spektrum         }
{ Amp .... vystupni fazove spektrum }
var i, j : Integer;
begin
  for i:=1 to N do
    for j:=1 to N do
      Amp[i,j]:=arg(Spec[i,j].Re,Spec[i,j].Im);
end;



procedure CrossPowerSpectrum(var Spec1, Spec2, Cross : Spectrum);
{ Spec1, Spec2 ... vstupni Fourierova spektra    }
{ Cross .......... vystupni cross-power spektrum }
var i, j : Integer;
      h  : single;
begin
  for i:=1 to N do
    for j:=1 to N do
      begin
        Cross[i,j].Re:=Spec1[i,j].Re*Spec2[i,j].Re +
                       Spec1[i,j].Im*Spec2[i,j].Im;
        Cross[i,j].Im:=Spec1[i,j].Im*Spec2[i,j].Re -
                       Spec1[i,j].Re*Spec2[i,j].Im;
        h:=sqrt(sqr(Cross[i,j].Re)+sqr(Cross[i,j].Im)); //mozno konstanty
        if h<0.00001 then h:=1.0;
        Cross[i,j].Re:=Cross[i,j].Re/h;
        Cross[i,j].Im:=Cross[i,j].Im/h;
      end;
end;


procedure CrossPowerSpectrumInPlace(var Spec1, Spec2 : Spectrum);
{ Spec1, Spec2 ... vstupni Fourierova spektra    }
{ Spec1 .......... vystupni cross-power spektrum }
var i, j : Integer;
      h  : single;
      z  : complex;
begin
  for i:=1 to N do
    for j:=1 to N do
      begin
        z.Re:=Spec1[i,j].Re*Spec2[i,j].Re +
              Spec1[i,j].Im*Spec2[i,j].Im;
        z.Im:=Spec1[i,j].Im*Spec2[i,j].Re -
              Spec1[i,j].Re*Spec2[i,j].Im;
        h:=sqrt(sqr(z.Re)+sqr(z.Im)); //mozno konstanty
        if h<0.00001 then h:=1.0;
        Spec1[i,j].Re:=z.Re/h;
        Spec1[i,j].Im:=z.Im/h;
      end;
end;


procedure CrossPowerSpectrumConstInPlace(var Spec1, Spec2 : Spectrum);
{ Spec1, Spec2 ... vstupni Fourierova spektra    }
{ Spec1 .......... vystupni cross-power spektrum }
var i, j     : Integer;
      c1, c2 : single;
      k1, k2,h : single;
      z      : complex;
begin
  c1:=100000;
  c2:=100000;
  for i:=1 to N do
    for j:=1 to N do
      begin
        k1:=sqrt(sqr(Spec1[i,j].Re)+sqr(Spec1[i,j].Im))+c1;
        k2:=sqrt(sqr(Spec2[i,j].Re)+sqr(Spec2[i,j].Im))+c2;

        z.Re:=Spec1[i,j].Re*Spec2[i,j].Re +
              Spec1[i,j].Im*Spec2[i,j].Im;
        z.Im:=Spec1[i,j].Im*Spec2[i,j].Re -
              Spec1[i,j].Re*Spec2[i,j].Im;
        h:=k1*k2;
        Spec1[i,j].Re:=z.Re/h;
        Spec1[i,j].Im:=z.Im/h;
      end;
end;



procedure CrossPowerSpectrumKoef(var Spec1, Spec2, Cross : Spectrum);
{ Spec1, Spec2 ... vstupni Fourierova spektra    }
{ Cross .......... vystupni cross-power spektrum }
var i, j       : Integer;
    h          : single;
    k, q, hk   : single ;

begin
  q:=65;
  k:=N*q;
  for i:=1 to N do
    for j:=1 to N do
      begin
        Cross[i,j].Re:=Spec1[i,j].Re*Spec2[i,j].Re +
                       Spec1[i,j].Im*Spec2[i,j].Im;
        Cross[i,j].Im:=Spec1[i,j].Im*Spec2[i,j].Re -
                       Spec1[i,j].Re*Spec2[i,j].Im;
        h:=sqrt(sqr(Cross[i,j].Re)+sqr(Cross[i,j].Im));
        h:=h+k;
        Cross[i,j].Re:=Cross[i,j].Re/h;
        Cross[i,j].Im:=Cross[i,j].Im/h;
      end;
end;



procedure RotationMagnification(var InvCPS : FFTImage; var Angle,
          Magnif, Quality : single);
var max, a    : single;
    i, j      : integer;
    dx, dy    : integer;
    pom       : double;
    sum, sums : double;
    mi, miq, sigma : double;
begin
  a:=exp((ln(N)-ln(2))/N);
  sums:=0.0;
  sum :=0.0;
  max:=-10000000000.0;
  for i:=1 to N do
    for j:=1 to N do
      begin
        pom:=InvCPS[i,j];
        sum :=sum  + pom;
        sums:=sums + pom*pom;
        if InvCPS[i,j]>=max then
          begin
            max:=InvCPS[i,j];
            dx:=i;
            dy:=j;
          end;
      end;
  mi:=sum/NN;
  miq:= sums/NN;
  sigma:=sqrt(miq - sqr(mi));
  if sigma<0.00001 then sigma:=1;
  Quality:=(max-mi)/sigma;
  dec(dx);
  dec(dy);
  Angle:=180.0*(dx/N);
  Magnif:=exp(dy*ln(a));
end;


procedure PermuteSpectrum(var SpecOld, SpecNew : Spectrum);
{ SpecOld ma nulovou frekvenci v prvku [0,0]     }
{ SpecNew ma nulovou frekvenci v prvku [N/2,N/2] }
var i, j, ii, jj : integer;
    M    : integer;
begin
  M:=N div 2;
  for i:=1 to N do
    for j:=1 to N do
      begin
        ii:=Permute(i);
        jj:=Permute(j);
        SpecNew[i,j].Re:=SpecOld[ii,jj].Re;
        SpecNew[i,j].Im:=SpecOld[ii,jj].Im;
      end;
end;


function FindPeak(var Res                 : FFTImage;
                  var Min, Max, Mi, Sigma : single) : Coordinate;
var i, j      : integer;
    dx, dy    : integer;
    pom       : double;
    sum, sums : double;
    miq       : double;
begin
  sums:=0.0;
  sum :=0.0;
  Max:=-10000000000.0;
  Min:= 10000000000.0;
  for i:=1 to N do
    for j:=1 to N do
      begin
        pom:=Res[i,j];
        sum :=sum  + pom;
        sums:=sums + pom*pom;
        if Res[i,j]>Max then
          begin
            Max:=Res[i,j];
            dx:=i;
            dy:=j;
          end;
        if Res[i,j]<Min then Min:=Res[i,j];
      end;
  Mi:=sum/NN;
  miq:= sums/NN;
  sigma:=sqrt(miq - sqr(mi));
  if sigma<0.00001 then sigma:=1;
  FindPeak.X:=dx; {v matici Res}
  FindPeak.Y:=dy; {v matici Res}
end;

// transformuje posuny na souradnice v Res : FFTImage
function InvTransfCoordinate(Pix : Coordinate) : Coordinate;
begin
  if Pix.X<0
    then InvTransfCoordinate.X:=Pix.X+N+1
    else InvTransfCoordinate.X:=Pix.X+1;
  if Pix.Y<0
    then InvTransfCoordinate.Y:=Pix.Y+N+1
    else InvTransfCoordinate.Y:=Pix.Y+1;
end;

// transformuje souradnice v Res : FFTImage na posuny
function TransfCoordinate(Pix : Coordinate) : Coordinate;
var x, y : integer;
begin
  x:=Pix.X-1; {protoze matice spektra je indexovana od 1}
  y:=Pix.Y-1;
  if x<=Ndiv2
    then TransfCoordinate.X:=x
    else TransfCoordinate.X:=x-N;
  if y<=Ndiv2
    then TransfCoordinate.Y:=y
    else TransfCoordinate.Y:=y-N;
end;


// bilinearni interpolace
function GetPixelVal(x,y : single; var Res : FFTImage) : single;
var x1, x2, y1, y2 : integer;
    r1, r2         : single;
begin
  if (x<1) or (x>=N) or (y<1) or (y>=N) then
    begin
      GetPixelVal:=0;
      Exit;
    end else
    begin
      x1:=trunc(x);
      x2:=x1+1;
      y1:=trunc(y);
      y2:=y1+1;
      r1:=Res[x1,y1]*(x2-x) + Res[x2,y1]*(x-x1);
      r2:=Res[x1,y2]*(x2-x) + Res[x2,y2]*(x-x1);
      GetPixelVal:=r1*(y2-y) + r2*(y-y1);
    end;
end;


function FindSubPixelShift(var Res : FFTImage; Peak : Coordinate; epsilon : single) : CoordinateS;
var i, j, x, y, eps, k     : integer;
    M00, M10, M01, dist, w : double;
    PeakOK, P, Q           : Coordinate;
    Centroid               : CoordinateS;
begin
  PeakOK:=TransfCoordinate(Peak);
  if epsilon<1 then
    begin  // celociselne souradnice peaku
      Result.X:=PeakOK.X;
      Result.Y:=PeakOK.Y;
      Exit;
    end;
  Centroid.X:=PeakOK.X;
  Centroid.Y:=PeakOK.Y;
  eps:=round(SubSample*epsilon);

  for k:=0 to NumSubPixelIter do
    begin
      M00:=0;
      M10:=0;
      M01:=0;
      for i:=-eps to eps  do
        for j:=-eps to eps do
          begin
            dist:=sqrt(sqr(i/SubSample)+sqr(j/SubSample));
            if dist<epsilon then
              begin
                P.X:=round(Centroid.X + i/SubSample);
                P.Y:=round(Centroid.Y + j/SubSample);
                Q:=InvTransfCoordinate(P);
                w:=Res[Q.X,Q.Y];
                if w>0 then
                  begin
                    M00:=M00 + w;
                    M10:=M10 + P.X*w;
                    M01:=M01 + P.Y*w;
                  end;
              end;
          end;
        if M00>0.0001 then
          begin
            Centroid.X:=M10/M00;
            Centroid.Y:=M01/M00;
          end;  
    end;
  if M00>0.0001 then
    begin
      Result.X:=Centroid.X;
      Result.Y:=Centroid.Y;
    end else
    begin
      Result.X:=PeakOK.X;
      Result.Y:=PeakOK.Y;
    end;
end;



function GetTranslation(var Img1,Img2,Res      : FFTImage;
                        LowPassFilter          : single;
                        SubEpsilon             : single;
                        OnlyInvFT              : boolean;
                        OnlySubPixelPeak       : boolean
                        ): Tra;
{ Img1, Img2 ... porovnavane obrazy typu array[1..N,1..N] of single  }
{ Res .......... vysledny obraz, ze ktereho je detekovan posuv,      }
{                ma smysl zobrazovat jen pro kontrolu nebo v pripade,}
{                ze GetTranslation.Quality je mala nebo pro rucni    }
{                detekci                                             }
{Max, Min ...... maximalni a minimalni hodnota pixelu v Result       }
{                (nutne pro vykresleni Result)                       }
{LowPasFilter... filtr dolni propust pouzity k vyhlazeni Res.        }
{                Pouziva se, kdyz je detekovane maximum nevyrazne    }
{                vzhledem k sumu.                                    }
{                0 filter nepouzit - nejcastejsi hodnota             }
{                1,2,4,8,16,32 postupne zvysujici se ucinnost filtu. }
{                Zadani jine hodnoty je totez co 0                   }
{OnlyInvFT ....  false znamena, ze probehne cela procedura           }
{                true predpoklada, ze Sp1 je jiz korekne naplneno a  }
{                     provede se algoritmus az od Low Pass filtru    }
{FUNKCE VRACI:                                                       }
{GetTranslation.X ... posun v pixlech od Img1 k Img2 v ose X         }
{GetTranslation.Y ... posun v pixlech od Img1 k Img2 v ose Y         }
{                     (osa Y je orientovana dolu)                    }
{GetTranslation.Quality ... duverihodnost nalezene translace         }
{                           Quality=Max/Sigma, kde Sigma je sm. odch.}
{                           hodnot pixelu v Result.                  }
{                           Quality<0 detekuje nepouzitelnost        }
var SubShift            : coordinates;
    Min, Max, Mi, Sigma : single;

begin
  if not OnlySubPixelPeak then
    begin
      if not OnlyInvFT then
        begin
          FourierTransform(Img1,Sp1);
          FourierTransform(Img2,Sp2);
          CrossPowerSpectrumInPlace(Sp1,Sp2);   //vysledek v Sp1
        end;
      if LowPassFilter>1.1 then
       begin
          LowPass(Sp1,Sp2,LowPassFilter);
          InverseFourierTransform(Sp2,Res);
       end else
       begin
        InverseFourierTransform(Sp1,Res);
       end;
    end;
   PeakCoordinate:=FindPeak(Res,Min,Max,Mi,Sigma); // v souradnicich Res
   SubShift:=FindSubPixelShift(Res,PeakCoordinate,SubEpsilon);
   GetTranslation.X:=SubShift.X;
   GetTranslation.Y:=SubShift.Y;
   GetTranslation.Quality:=Max/Sigma;

  PeakX:=PeakCoordinate.X;
  PeakY:=PeakCoordinate.Y;
  PeakSubPixX:=SubShift.X;
  PeakSubPixY:=SubShift.Y;
  DisplayPhaseCorrResult(Form1.ImagePhaseCorr,Form1.ScrollBoxPhaseCorr,Max,Min,0,1);
end;


function FindMagnifRot(SubShift : coordinates) : MagRot;
var a : double;
begin
  a:=exp((ln(N)-ln(2))/N);
  FindMagnifRot.Rotation:=180*(SubShift.X/N);
  FindMagnifRot.Magnification:=exp(SubShift.Y*ln(a));
end;


function GetMagnifRot(var Img1,Img2,Res       : FFTImage;
                      LowPassFilter           : single;
                      SubEpsilon              : single;
                      OnlyInvFT               : boolean;
                      OnlySubPixelPeak        : boolean
                      ) : MagRot;
var SubShift : Coordinates;
    Min, Max, Mi, Sigma : single;                     
begin
  if not OnlySubPixelPeak then
    begin
      if not OnlyInvFT then
        begin
          FourierTransform(Img1,Sp1);
          LogAmplitudeSpectrumPerm(Sp1,Res);
          LogPolarTransfNew(Res,Pom);
          FourierTransform(Pom,Sp1);
          FourierTransform(Img2,Sp2);
          LogAmplitudeSpectrumPerm(Sp2,Res);
          LogPolarTransfNew(Res,Pom);
          FourierTransform(Pom,Sp2);
          CrossPowerSpectrumInPlace(Sp1,Sp2); //vysledek v Sp1
        end;
      if LowPassFilter>1.1 then
       begin
          LowPass(Sp1,Sp2,LowPassFilter);
          InverseFourierTransform(Sp2,Res);
       end else
       begin
        InverseFourierTransform(Sp1,Res);
       end;
      //PeakCoordinate:=FindPeak(Res,Min,Max,Mi,Sigma); // v souradnicich Res
      //PeakX:=PeakCoordinate.X;
      //PeakY:=PeakCoordinate.Y;
    end;
  PeakCoordinate:=FindPeak(Res,Min,Max,Mi,Sigma); // v souradnicich Res
  PeakX:=PeakCoordinate.X;
  PeakY:=PeakCoordinate.Y;
  SubShift:=FindSubPixelShift(Res,PeakCoordinate,SubEpsilon);
  GetMagnifRot:=FindMagnifRot(SubShift);
  GetMagnifRot.Quality:=Max/Sigma;



  {nutne puze pro vykresleni detailu a krizku na miste subpixeloveho peaku}
  PeakSubPixX:=SubShift.X;
  PeakSubPixY:=SubShift.Y;
  DisplayPhaseCorrResult(Form1.ImagePhaseCorr,Form1.ScrollBoxPhaseCorr,Max,Min,0,1);
end;

// ------------------------------ Nove --------------------------------------------------------------

procedure SetHanningCircleWindow(size : single; ImgRoiRect : TRect; var Img : FFTImage);
//size je v polomerech (napr. 0.2 znamena ze na 80% procentech polomeru je vaha 1, na 20% klesne vaha na 0)
var i, j, w    : integer;
    d, x, y, q : single;
begin
  w:=ImgRoiRect.Right-ImgRoiRect.Left;
  if w>N then w:=N;
  for i:=1 to N do
    for j:=1 to N do
      begin
       x:=i-0.5;
       y:=j-0.5;
       d:=1-2*sqrt(sqr(x-w/2)+sqr(y-w/2))/w;
       if d<0 then Img[i,j]:=0 else
         if d<size then Img[i,j]:=Img[i,j]*0.5*(sin(pi*d/size - pi/2)+1);
      end;
end;

procedure SetHanningSquareWindow(size : single; ImgRoiRect : TRect; var Img : FFTImage);
//size je v polomerech (napr. 0.2 znamena ze na 80% procentech polomeru je vaha 1, na 20% klesne vaha na 0)
var i, j, w             : integer;
    dX, dY, x, y, qX, qY : single;
begin
  if size<0.0001 then Exit;
  w:=ImgRoiRect.Right-ImgRoiRect.Left;
  if w>N then w:=N;
  for i:=1 to N do
    for j:=1 to N do
      begin
       x:=i-0.5;
       y:=j-0.5;
       dX:=1-2*abs(x-w/2)/w;
       dY:=1-2*abs(y-w/2)/w;
       if dX<size then qX:=0.5*(sin(pi*dX/size - pi/2)+1);
       if dY<size then qY:=0.5*(sin(pi*dY/size - pi/2)+1);
       Img[i,j]:=Img[i,j]*qX*qY;
      end;
end;



procedure CopyROI(var Img : FFTImage; ImgPointer : pointer; ImgWidth : integer; ImgHeight : integer;
                  ImgChannels : integer; ImgDataType : string; ImgRoiRect  : TRect);

var i, j, x1, x2, y1, y2, w, h, ww, k : integer;
    ByteArr        : PByteArr;
    WordArr        : PWordArr;
    SmallintArr    : PSmallintArr;
    SingleArr      : PSingleArr;
    DoubleArr      : PDoubleArr;
    ByteClrArr     : PByteClrArr;
    WordClrArr     : PWordClrArr;
    SmallintClrArr : PSmallintClrArr;
    SingleClrArr   : PSingleClrArr;
    DoubleClrArr   : PDoubleClrArr;
begin
  x1:=ImgRoiRect.Left;
  x2:=ImgRoiRect.Right;
  y1:=ImgRoiRect.Top;
  y2:=ImgRoiRect.Bottom;
  w:=x2-x1;
  h:=y2-y1;
  ww:=ImgWidth;
  for i:=1 to N do for j:=1 to N do Img[i,j]:=0; // nulovani
  {B&W obrazy}
  if (ImgChannels=1) and (UpperCase(ImgDataType)='BYTE') then
    begin
      ByteArr:=ImgPointer;
      for j:=0 to h-1 do for i:=0 to w-1 do Img[i+1,j+1]:=ByteArr^[x1+i + (y1+j)*ww];
    end;
  if (ImgChannels=1) and (UpperCase(ImgDataType)='WORD') then
    begin
      WordArr:=ImgPointer;
      for j:=0 to h-1 do for i:=0 to w-1 do Img[i+1,j+1]:=WordArr^[x1+i + (y1+j)*ww];
    end;
  if (ImgChannels=1) and (UpperCase(ImgDataType)='SMALLINT') then
    begin
      SmallintArr:=ImgPointer;
      for j:=0 to h-1 do for i:=0 to w-1 do Img[i+1,j+1]:=SmallintArr^[x1+i + (y1+j)*ww];
    end;
  if (ImgChannels=1) and (UpperCase(ImgDataType)='SINGLE') then
    begin
      SingleArr:=ImgPointer;
      for j:=0 to h-1 do for i:=0 to w-1 do Img[i+1,j+1]:=SingleArr^[x1+i + (y1+j)*ww];
    end;
  if (ImgChannels=1) and (UpperCase(ImgDataType)='DOUBLE') then
    begin
      DoubleArr:=ImgPointer;
      for j:=0 to h-1 do for i:=0 to w-1 do Img[i+1,j+1]:=DoubleArr^[x1+i + (y1+j)*ww];
    end;
  {barevne obrazy}
  if (ImgChannels=3) and (UpperCase(ImgDataType)='BYTE') then
    begin
      ByteClrArr:=ImgPointer;
      for j:=0 to h-1 do for i:=0 to w-1 do
        begin
          k:=x1+i + (y1+j)*ww;
          Img[i+1,j+1]:=wR*ByteClrArr^[k].R + wG*ByteClrArr^[k].G + wB*ByteClrArr^[k].B;
        end;
    end;
  if (ImgChannels=3) and (UpperCase(ImgDataType)='WORD') then
    begin
      WordClrArr:=ImgPointer;
      for j:=0 to h-1 do for i:=0 to w-1 do
        begin
          k:=x1+i + (y1+j)*ww;
          Img[i+1,j+1]:=wR*WordClrArr^[k].R + wG*WordClrArr^[k].G + wB*WordClrArr^[k].B;
        end;
    end;
  if (ImgChannels=3) and (UpperCase(ImgDataType)='SMALLINT') then
    begin
      SmallintClrArr:=ImgPointer;
      for j:=0 to h-1 do for i:=0 to w-1 do
        begin
          k:=x1+i + (y1+j)*ww;
          Img[i+1,j+1]:=wR*SmallintClrArr^[k].R + wG*SmallintClrArr^[k].G + wB*SmallintClrArr^[k].B;
        end;
    end;
  if (ImgChannels=3) and (UpperCase(ImgDataType)='SINGLE') then
    begin
      SingleClrArr:=ImgPointer;
      for j:=0 to h-1 do for i:=0 to w-1 do
        begin
          k:=x1+i + (y1+j)*ww;
          Img[i+1,j+1]:=wR*SingleClrArr^[k].R + wG*SingleClrArr^[k].G + wB*SingleClrArr^[k].B;
        end;
    end;
   if (ImgChannels=3) and (UpperCase(ImgDataType)='DOUBLE') then
    begin
      DoubleClrArr:=ImgPointer;
      for j:=0 to h-1 do for i:=0 to w-1 do
        begin
          k:=x1+i + (y1+j)*ww;
          Img[i+1,j+1]:=wR*DoubleClrArr^[k].R + wG*DoubleClrArr^[k].G + wB*DoubleClrArr^[k].B;
        end;
    end;
end;



function FindImgMagnRot(Img1Pointer      : pointer; // pointer na 1. obraz
                        Img1Width        : integer; // sirka 1. obrazu v pixelech
                        Img1Height       : integer; // vyska 1. obrazu v pixelech
                        Img1Channels     : integer; // kolik slozek ma pixel 1. obrazu (1 nebo 3)
                        Img1DataType     : string;  // datovy typ 1. obrazu ('BYTE' nebo 'WORD' nebo 'SMALLINT' nebo 'SINGLE' nebo 'DOUBLE');
                        Img1RoiRect      : TRect;   // obdelnik v 1. obraze pouzity k urceni transformace
                        Img2Pointer      : pointer; // pointer na 2. obraz
                        Img2Width        : integer; // sirka 2. obrazu v pixelech
                        Img2Height       : integer; // vyska 2. obrazu v pixelech
                        Img2Channels     : integer; // kolik slozek ma pixel 2. obrazu (1 nebo 3)
                        Img2DataType     : string;  // datovy typ 2. obrazu ('BYTE' nebo 'WORD' nebo 'SMALLINT' nebo 'SINGLE' nebo 'DOUBLE');
                        Img2RoiRect      : TRect;   // obdelnik v 2. obraze pouzity k urceni transformace
                        LowPassFilter    : single;  // parametr LP filtru
                        SubEpsilon       : single;  // polomer okoli pro sub-pixelove hledani maxima
                        CircleRoi        : boolean; // false ... square ROI, true ... circle ROI
                        WindEdge         : single;  // parametr pro prubeh Hanningova okna (0 .. ostry okraj, 1 .. plynuly pokles od stredu okna)
                        OnlyInvFT        : boolean; // false ... probehne cela procedura, true ... zacne se LP filtrem (pouziva se pokud uz
                                                    // procedura jednou probehla a byl zmenen jen LP filtr)
                        OnlySubPixelPeak : boolean  // true ... hleda se znovy jen jen sub-pixelove maximum
                                        ): MagRot;  // vysledek - zvetseni a otoceni
// vysledek fazove korealce (peak) lze zobrazit pomoci procedury DisplayPhaseCorrResult
begin
  LPfilter:=exp(LowPassFilter*ln(2));
  EdgeHanning:=WindEdge;
  CopyROI(Img1,Img1Pointer,Img1Width,Img1Height,Img1Channels,Img1DataType,Img1RoiRect);
  CopyROI(Img2,Img2Pointer,Img2Width,Img2Height,Img2Channels,Img2DataType,Img2RoiRect);
  if CircleRoi then
    begin
      SetHanningCircleWindow(EdgeHanning,Img2RoiRect,Img1);
      SetHanningCircleWindow(EdgeHanning,Img2RoiRect,Img2);
    end else
    begin
      SetHanningSquareWindow(EdgeHanning,Img1RoiRect,Img1);
      SetHanningSquareWindow(EdgeHanning,Img1RoiRect,Img2);
    end;
   FindImgMagnRot:=GetMagnifRot(Img1,Img2,Res,LPfilter,SubEpsilon,OnlyInvFT,OnlySubPixelPeak);
end;




function FindImgTranslation(Img1Pointer      : pointer;  // pointer na 1. obraz
                            Img1Width        : integer;  // sirka 1. obrazu v pixelech
                            Img1Height       : integer;  // vyska 1. obrazu v pixelech
                            Img1Channels     : integer;  // kolik slozek ma pixel 1. obrazu (1 nebo 3)
                            Img1DataType     : string;   // datovy typ 1. obrazu ('BYTE' nebo 'WORD' nebo 'SMALLINT' nebo 'SINGLE' nebo 'DOUBLE');
                            Img1RoiRect      : TRect;    // obdelnik v 1. obraze pouzity k urceni posuvu
                            Img2Pointer      :  pointer; // pointer na 2. obraz
                            Img2Width        : integer;  // sirka 2. obrazu v pixelech
                            Img2Height       : integer;  // vyska 2. obrazu v pixelech
                            Img2Channels     : integer;  // kolik slozek ma pixel 2. obrazu (1 nebo 3)
                            Img2DataType     : string;   // datovy typ 2. obrazu ('BYTE' nebo 'WORD' nebo 'SMALLINT' nebo 'SINGLE' nebo 'DOUBLE');
                            Img2RoiRect      : TRect;    // obdelnik v 2. obraze pouzity k urceni posuvu
                            LowPassFilter    : single;   // parametr LP filtru
                            SubEpsilon       : single;   // polomer okoli pro sub-pixelove hledani maxima
                            CircleRoi        : boolean;  // false ... square ROI, true ... circle ROI
                            WindEdge         : single;   // parametr pro prubeh Hanningova okna (0 .. ostry okraj, 1 .. plynuly pokles od stredu okna)
                            OnlyInvFT        : boolean;  // false ... probehne cela procedura, true ... zacne se LP filtrem (pouziva se pokud uz
                                                         // procedura jednou probehla a byl zmenen jen LP filtr)
                            OnlySubPixelPeak : boolean   // true ... hleda se znovy jen jen sub-pixelove maximum
                                          )  : Tra;      // vysledkem je record s posuvem
begin
  LPfilter:=exp(LowPassFilter*ln(2));
  EdgeHanning:=WindEdge;
  CopyROI(Img1,Img1Pointer,Img1Width,Img1Height,Img1Channels,Img1DataType,Img1RoiRect);
  CopyROI(Img2,Img2Pointer,Img2Width,Img2Height,Img2Channels,Img2DataType,Img2RoiRect);
  if CircleRoi then
    begin
      SetHanningCircleWindow(EdgeHanning,Img2RoiRect,Img1);
      SetHanningCircleWindow(EdgeHanning,Img2RoiRect,Img2);
    end else
    begin
      SetHanningSquareWindow(EdgeHanning,Img1RoiRect,Img1);
      SetHanningSquareWindow(EdgeHanning,Img1RoiRect,Img2);
    end;
 // DisplayFFTImage(Img1,Form1.ImageHidden);
  AutoTransl:=GetTranslation(Img1,Img2,Res,LPfilter,SubEpsilon,OnlyInvFT,OnlySubPixelPeak);
  FindImgTranslation.X:=AutoTransl.X+Img1RoiRect.Left-Img2RoiRect.Left;
  FindImgTranslation.Y:=AutoTransl.Y+Img1RoiRect.Top-Img2RoiRect.Top;
  FindImgTranslation.Quality:=AutoTransl.Quality;
end;

procedure DisplayPhaseCorrResult(Image     :  TImage;    // Bitmapa (pf24bit) pro zobrazeni vysledku fazove korelace
                                 ScrollBox : TScrollBox; // procedura sama scrolluje na peak, pokud je Bitmapa ve ScrollBoxu,
                                                         // pokud neni, musi byt ScrollBox=nil
                                 Max, Min  : single;
                                 Black     : single;     // 0<=Black<White, Black=0 znamena, ze min hodnota je cerna
                                 White     : single);    // Black<White<=1, <White=1 znamena, ze max hodnota (peak) je bila
// zobrazuje vysledek fazove korealce (peak)
var i, j, ii, jj      : integer;
    ScrollX, ScrollY  : integer;
    NormPix, k, Range : single;
    SetMax, SetMin    : single;
    b                 : byte;
    Buff              : array [1..3*N] of byte; // buffer pro kresleni obrazu
begin
  if Black<0 then Black:=0;
  if White>1 then White:=1;
  if Black>White then
    begin
      k:=Black;
      Black:=White;
      White:=k;
    end;
  SetMin:=Min+Black*(Max-Min);
  SetMax:=Min+White*(Max-Min);
  Range:=SetMax-SetMin;
  if Range<0.0001 then Range:=0001;
  Image.Picture.Bitmap.Width:=N;
  Image.Picture.Bitmap.Height:=N;
  Image.Width:=N;
  Image.Height:=N;
  Image.Picture.Bitmap.PixelFormat:=pf24bit;

  for j:=1 to N do
    begin
     for i:=1 to N do
       begin
         ii:=Permute(i);
         jj:=Permute(j);
         NormPix:=(Res[ii,jj]-SetMin)/Range;
         if NormPix>1 then NormPix:=1;
         if NormPix<0 then NormPix:=0;
         b:=round(255.0*NormPix);
         Buff[3*i]:=b;
         NormPix:=sqr(NormPix);
         NormPix:=sqr(NormPix);
         b:=round(255.0*NormPix);
         Buff[3*i-1]:=b;
         NormPix:=sqr(NormPix);
         NormPix:=sqr(NormPix);
         NormPix:=sqr(NormPix);
         b:=round(255.0*NormPix);
         Buff[3*i-2]  :=b;
       end;
       Move(Buff,Image.Picture.Bitmap.ScanLine[j-1]^,SizeOf(Buff));
     end;
  PeakX:=Permute(PeakCoordinate.X);
  PeakY:=Permute(PeakCoordinate.Y);
  if ScrollBox=nil then Exit;
  ScrollX:=PeakX - ScrollBox.Width div 2;
  ScrollY:=PeakY - ScrollBox.Height div 2;
  if ScrollX<0 then ScrollX:=0;
  if ScrollY<0 then ScrollY:=0;
  if ScrollX>ScrollBox.HorzScrollBar.Range then ScrollX:=ScrollBox.HorzScrollBar.Range;
  if ScrollY>ScrollBox.VertScrollBar.Range then ScrollY:=ScrollBox.VertScrollBar.Range;  ScrollBox.HorzScrollBar.Position:=ScrollX;
  ScrollBox.VertScrollBar.Position:=ScrollY;
  Image.Repaint;
end;


end.


