iefft.pas

·ImageEn 2.3.0 ImageEn一组用于图像处理、查看和分析的Delphi控件。能够保存几种图像格式

(*
Copyright (c) 1998-2007 by HiComponents. All rights reserved.

This software comes without express or implied warranty.
In no case shall the author be liable for any damage or unwanted behavior of any
computer hardware and/or software.

HiComponents grants you the right to include the compiled component
in your application, whether COMMERCIAL, SHAREWARE, or FREEWARE,
BUT YOU MAY NOT DISTRIBUTE THIS SOURCE CODE OR ITS COMPILED .DCU IN ANY FORM.

ImageEn, IEvolution and ImageEn ActiveX may not be included in any commercial,
shareware or freeware libraries or components.

email: [e-mail]

http://www.hicomponents.com
*)

unit iefft;

{$R-}
{$Q-}

{$I ie.inc}

interface

uses Windows, Graphics, Classes, hyiedefs, hyieutils;

type

  TIEsinglearray = array[0..16384] of single;
  PIEsinglearray = ^TIEsinglearray;

  TIECOMPLEX_IMAGE = array[0..16384] of PIEsinglearray;
  PIECOMPLEX_IMAGE = ^TIECOMPLEX_IMAGE;

  PIEdouble = ^double;
  PIEinteger = ^integer;
  PIEsingle = ^single;
  PIElongint = ^longint;

  TIEComplexColor = packed record
    // red channel
    real_Red: PIEsingle;
    imag_Red: PIEsingle;
    // blue channel
    real_Blue: PIEsingle;
    imag_Blue: PIEsingle;
    // green channel
    real_Green: PIEsingle;
    imag_Green: PIEsingle;
    // gray scale
    imag_gray: PIEsingle;
    real_gray: PIEsingle;
  end;

  TIEComplex = packed record
    real: single;
    imag: single;
  end;

{!!
<FS>TIEFtImageType

<FM>Declaration<FC>
}
  TIEFtImageType = (ieitUnknown, ieitRGB, ieitGrayscale);
{!!}

  tdwordarray = array[0..maxint div 16] of dword;
  pdwordarray = ^tdwordarray;

  tlongintarray = array[0..maxint div 16] of longint;
  plongintarray = ^tlongintarray;

  TIEFtImage = class
  private
    sintabpt: PIEdouble;
    bittabpt: PIEinteger;
    powers: PIEinteger;
    numpts: integer;
    nn: integer;
    direction: single;
    scalef: double;
    FFTN: integer;
    NORMALIZE: boolean;
    //
    fftr: PIECOMPLEX_IMAGE; // FT of R
    fftg: PIECOMPLEX_IMAGE; // FT of G
    fftb: PIECOMPLEX_IMAGE; // FT of B
    fftgray: PIECOMPLEX_IMAGE; // FT of gray scale
    fOX, fOY: integer; // origins of fOrigBitmap inside fBitmap
    fOrigBitmapWidth, fOrigBitmapHeight: integer; // original input bitmap sizes
    fOnProgress: TIEProgressEvent;
    fImageType: TIEFtImageType;
    //
    //procedure image_fft(image:TBitmap; var output:PIECOMPLEX_IMAGE; ch:integer);
    procedure image_fftoc(image: TIEBitmap; var output: PIECOMPLEX_IMAGE; ch: integer; var Progress: TProgressRec);
    procedure fftx(image: TIEBitmap; oc: integer; var output: PIECOMPLEX_IMAGE; ch: integer; var Progress: TProgressRec);
    procedure image_fftinv(image: PIECOMPLEX_IMAGE; var output: PIECOMPLEX_IMAGE; var Progress: TProgressRec);
    procedure pairsort(arr: PIEsinglearray; iarr: pdwordarray; n: integer);
    procedure fqsort(arr: PIEsinglearray; iarr: pdwordarray; l: integer; r: integer);
    //procedure image_fftinvoc(image:PIECOMPLEX_IMAGE; var output:PIECOMPLEX_IMAGE);
    procedure fftinvx(image: PIECOMPLEX_IMAGE; oc: integer; var output: PIECOMPLEX_IMAGE; var Progress: TProgressRec);
    procedure fft2d(image: PIECOMPLEX_IMAGE; direction: single; var Progress: TProgressRec);
    procedure filt_orig(xarray: PIECOMPLEX_IMAGE);
    procedure realtoint(fim: PIECOMPLEX_IMAGE; H: plongintarray);
    procedure fft(data: PIEsinglearray; dir: single);
    procedure _fft(tseries: PIEsinglearray; level: integer; chunk: integer);
    procedure fftinit(nopts: integer);
    function bitrev(bits: integer): integer;
    function newcomplex(im: TIEBitmap; ch: integer): PIECOMPLEX_IMAGE;
    function dupcomplex(im: PIECOMPLEX_IMAGE): PIECOMPLEX_IMAGE;
    procedure filt_toint(oimage: PIECOMPLEX_IMAGE; output: TIEBitmap; ch: integer);
    function GetComplexImage(x, y: integer): TIEComplexColor;
  protected
    procedure FreeAll;
  public
    constructor Create;
    destructor Destroy; override;
    // transformations
    procedure BuildFT(fOrigBitmap: TIEBitmap; ImageType: TIEFtImageType); // direct Fourier transf
    procedure BuildBitmap(Bitmap: TIEBitmap); // inverse Fourier transf
    procedure GetFTImage(Bitmap: TIEBitmap); // display an image of Fourier transf
    // processing transformed image
    procedure HiPass(radius: integer);
    procedure LoPass(radius: integer);
    procedure ClearZone(x1, y1, x2, y2: integer);
    // complex image access
    property ComplexPixel[x, y: integer]: TIEComplexColor read GetComplexImage;
    property Imagetype: TIEFtImageType read fImageType;
    property ComplexWidth: integer read FFTN;
    property ComplexHeight: integer read FFTN;
    procedure Assign(Source: TIEFtImage);
    // events
    property OnProgress: TIEProgressEvent read fOnProgress write fOnProgress;
  end;

implementation

uses SysUtils, imageenproc;

{$R-}

procedure freecomplex(var x: PIECOMPLEX_IMAGE);
begin
  if (x <> nil) then
  begin
    freemem(x[0]);
    freemem(x);
    x := nil;
  end;
end;

constructor TIEFtImage.Create;
begin
  inherited Create;
  sintabpt := nil;
  bittabpt := nil;
  powers := nil;
  scalef := 0.0;
  FFTN := 0;
  NORMALIZE := false;
  //
  fftr := nil;
  fftg := nil;
  fftb := nil;
  fftgray := nil;
  fOnProgress := nil;
  fImageType := ieitUnknown;
end;

procedure TIEFtImage.FreeAll;
begin
  if fftr <> nil then
    freecomplex(fftr);
  if fftg <> nil then
    freecomplex(fftg);
  if fftb <> nil then
    freecomplex(fftb);
  if fftgray <> nil then
    freecomplex(fftgray);
end;

destructor TIEFtImage.Destroy;
begin
  FreeAll;
  if (sintabpt <> nil) then
    freemem(sintabpt);
  if (bittabpt <> nil) then
    freemem(bittabpt);
  if (powers <> nil) then
    freemem(powers);
  //
  inherited;
end;

procedure TIEFtImage.Assign(Source: TIEFtImage);
begin
  FreeAll;
  FFTN := Source.FFTN;
  NORMALIZE := Source.NORMALIZE;
  fftinit(FFTN);
  fOX := Source.fOX;
  fOY := Source.fOY;
  fOrigBitmapWidth := Source.fOrigBitmapWidth;
  fOrigBitmapHeight := Source.fOrigBitmapHeight;
  fOnProgress := Source.fOnProgress;
  fImageType := Source.fImageType;
  if assigned(Source.fftr) then
    fftr := dupcomplex(Source.fftr);
  if assigned(Source.fftg) then
    fftg := dupcomplex(Source.fftg);
  if assigned(Source.fftb) then
    fftb := dupcomplex(Source.fftb);
  if assigned(Source.fftgray) then
    fftgray := dupcomplex(Source.fftgray);
end;

// converts from bitmap to FT
// channel:
// 0=b 1=g 2=r 3=rgb 4=grayscale

procedure TIEFtImage.BuildFT(fOrigBitmap: TIEBitmap; ImageType: TIEFtImageType);
var
  fBitmap: TIEBitmap;
  ww, i: integer;
  Progress: TProgressRec;
begin
  fImageType := ImageType;
  Progress.fOnProgress := fOnProgress;
  Progress.Sender := Self;
  if fftr <> nil then
    freecomplex(fftr);
  if fftg <> nil then
    freecomplex(fftg);
  if fftb <> nil then
    freecomplex(fftb);
  if fftgray <> nil then
    freecomplex(fftgray);
  //
    //BuildOrigHistogram;	// no best results making an histogram!
  fOrigBitmapWidth := fOrigBitmap.Width;
  fOrigBitmapHeight := fOrigBitmap.Height;
  ww := imax(fOrigBitmapWidth, fOrigBitmapHeight);
  i := 0;
  while (1 shl i) < ww do
    inc(i);
  ww := 1 shl i;
  fBitmap := TIEBitmap.Create;
  fBitmap.Allocate(ww, ww, fOrigBitmap.PixelFormat);
  fBitmap.Fill(0);
  fOX := (ww - fOrigBitmapWidth) div 2;
  fOY := (ww - fOrigBitmapHeight) div 2;
  fOrigBitmap.CopyRectTo(fBitmap, 0, 0, fOX, fOY, fOrigBitmap.Width, fOrigBitmap.Height);
  //
  normalize := true;
  case fImageType of
    ieitRGB: // RGB
      begin
        Progress.tot := 3;
        Progress.val := 0;
        image_fftoc(fBitmap, fftr, 2, Progress);
        image_fftoc(fBitmap, fftg, 1, Progress);
        image_fftoc(fBitmap, fftb, 0, Progress);
      end;
    ieitGrayScale: // GRAYSCALE
      begin
        Progress.tot := 1;
        Progress.val := 0;
        image_fftoc(fBitmap, fftgray, 3, Progress);
      end;
  end;
  normalize := false;
  FreeAndNil(fBitmap);
end;

// Bitmap must be created

procedure TIEFtImage.GetFTImage(Bitmap: TIEBitmap);
begin
  Bitmap.Allocate(FFTN, FFTN, ie24RGB);
  case fImageType of
    ieitRGB:
      begin
        // RGB
        filt_toint(fftr, Bitmap, 2);
        filt_toint(fftg, Bitmap, 1);
        filt_toint(fftb, Bitmap, 0);
      end;
    ieitGrayscale:
      // gray scale
      filt_toint(fftgray, Bitmap, 3);
  end;
end;

// converts from FT to bitmap
// Bitmap.PixleFormat will be ie24RGB

procedure TIEFtImage.BuildBitmap(Bitmap: TIEBitmap);
var
  offtr, offtg, offtb: PIECOMPLEX_IMAGE;
  x, y: integer;
  rgb: PRGB;
  Progress: TProgressRec;
  bitmapheight1, bitmapwidth1: integer;
begin
  Progress.fOnProgress := fOnProgress;
  Progress.Sender := Self;
  Bitmap.Allocate(fOrigBitmapWidth, fOrigBitmapHeight, ie24RGB);
  case fImageType of
    ieitRGB:
      begin
        // RGB
        Progress.tot := 3;
        Progress.val := 0;
        image_fftinv(fftr, offtr, Progress);
        realtoint(offtr, nil);
        image_fftinv(fftg, offtg, Progress);
        realtoint(offtg, nil);
        image_fftinv(fftb, offtb, Progress);
        realtoint(offtb, nil);
        bitmapheight1 := Bitmap.Height - 1;
        bitmapwidth1 := Bitmap.Width - 1;
        for y := 0 to BitmapHeight1 do
        begin
          rgb := Bitmap.Scanline[y];
          for x := 0 to BitmapWidth1 do
          begin
            rgb^.r := trunc(offtr^[y + fOY]^[x + fOX]);
            rgb^.g := trunc(offtg^[y + fOY]^[x + fOX]);
            rgb^.b := trunc(offtb^[y + fOY]^[x + fOX]);
            inc(rgb);
          end;
        end;
        freecomplex(offtr);
        freecomplex(offtg);
        freecomplex(offtb);
      end;
    ieitGrayscale:
      begin
        // gray scale
        Progress.tot := 1;
        Progress.val := 0;
        image_fftinv(fftgray, offtg, Progress);
        realtoint(offtg, nil);
        bitmapheight1 := Bitmap.Height - 1;
        bitmapwidth1 := Bitmap.Width - 1;
        for y := 0 to BitmapHeight1 do
        begin
          rgb := Bitmap.Scanline[y];
          for x := 0 to BitmapWidth1 do
          begin
            rgb^.r := trunc(offtg^[y + fOY]^[x + fOX]);
            rgb^.g := rgb^.r;
            rgb^.b := rgb^.r;
            inc(rgb);
          end;
        end;
        freecomplex(offtg);
      end;
  end;
end;

(*
// ch: channel of im = 0:B 1:G 2:R 3:grayscale
procedure TIEFtImage.image_fft(image:TIEBitmap; var output:PIECOMPLEX_IMAGE; ch:integer);
begin
 fftx(image, 0, output,ch);
end;
*)

// ch: channel of im = 0:B 1:G 2:R 3:grayscale

procedure TIEFtImage.image_fftoc(image: TIEBitmap; var output: PIECOMPLEX_IMAGE; ch: integer; var Progress: TProgressRec);
begin
  fftx(image, 1, output, ch, Progress);
end;

// Do a 2D FFT on a set of BYTE-type pixels
// ch: channel of im = 0:B 1:G 2:R 3:grayscale

procedure TIEFtImage.fftx(image: TIEBitmap; oc: integer; var output: PIECOMPLEX_IMAGE; ch: integer; var Progress: TProgressRec);
var
  n: integer;
  cim: PIECOMPLEX_IMAGE;
begin
  n := image.Width;
  if (FFTN <> n) then
    fftinit(n);
  FFTN := n;
  cim := newcomplex(image, ch);
  if (oc <> 0) then
    filt_orig(cim);
  fft2d(cim, -1.0, Progress);
  output := cim;
end;

procedure TIEFtImage.image_fftinv(image: PIECOMPLEX_IMAGE; var output: PIECOMPLEX_IMAGE; var Progress: TProgressRec);
begin
  fftinvx(image, 0, output, Progress);
end;

(*
procedure TIEFtImage.image_fftinvoc(image:PIECOMPLEX_IMAGE; var output:PIECOMPLEX_IMAGE);
begin
 fftinvx(image, 1, output);
end;
*)

procedure TIEFtImage.fftinvx(image: PIECOMPLEX_IMAGE; oc: integer; var output: PIECOMPLEX_IMAGE; var Progress: TProgressRec);
var
  cim: PIECOMPLEX_IMAGE;
begin
  cim := dupcomplex(image);
  if (oc <> 0) then
    filt_orig(cim);
  fft2d(cim, 1.0, Progress);
  output := cim;
end;

// Complex product: res = c1*c2
(*
procedure cprod(c1r:single; c1i:single; c2r:PIEsingle; c2i:PIEsingle);
var
 real, imag:single;
begin
 real := c1r*(c2r^) - c1i*(c2i^);
 imag := c1r*(c2i^) + (c2r^)*c1i;
 c2r^ := real;
   c2i^ := imag;
end;
*)

// abs value squared: res = |c1|**2

function pix_cnorm(cr: single; ci: single): single;
begin
  result := cr * cr + ci * ci;
end;

// res = c1/c2 ... Complex
(*
procedure cdiv(c1r:single; c1i:single; c2r:PIEsingle; c2i:PIEsingle);
var
 z,real, imag:single;
begin
   z := ((c2r^)*(c2r^) + (c2i^)*(c2i^));
 if (z <> 0.0) then begin
  real := (c1r*(c2r^) + c1i*(c2i^))/z;
  imag := ((c2r^)*c1i - c1r*(c2i^))/z;
  c2r^ := real;
  c2i^ := imag;
 end else begin
  c2r^ := 0.0;
      c2i^ := 0.0;
   end;
end;
*)

//  fft2d -- Calls `fft' to perform a 2 dimensional Fast Fourier Trans-
//           form on an N x N array of complex data. (N must be defined).

procedure TIEFtImage.fft2d(image: PIECOMPLEX_IMAGE; direction: single; var Progress: TProgressRec);
var
  temp: piesinglearray;
  i, j: integer;
begin
  Progress.per1 := 100 / (FFTN * Progress.Tot);
  getmem(temp, (FFTN * 2 + 1) * 4);
  for i := 0 to FFTN - 1 do
    fft(image^[i], direction);
  for i := 0 to FFTN - 1 do
  begin
    for j := 0 to FFTN - 1 do
    begin
      temp^[j] := image^[j]^[i];
      temp^[j + FFTN] := image^[j]^[i + FFTN];
    end;
    fft(temp, direction);
    for j := 0 to FFTN - 1 do
    begin
      image^[j]^[i] := temp^[j];
      image^[j]^[i + FFTN] := temp^[j + FFTN];
    end;
    with Progress do
      if assigned(fOnProgress) then
        fOnProgress(Sender, trunc(per1 * i + per1 * Progress.Val * FFTN));
  end;
  freemem(temp);
  inc(Progress.Val);
end;

//  Fixorig -- Modifies the input data so that after being Fourier trans-
//             formed the origin (f = 0) lies at the center of the array.

procedure TIEFtImage.filt_orig(xarray: PIECOMPLEX_IMAGE);
var
  i, j: integer;
begin
  for i := 0 to FFTN - 1 do
    for j := 0 to FFTN - 1 do
      if ((i + j) mod 2) <> 0 then
      begin
        xarray^[i]^[j] := -xarray^[i]^[j];
        xarray^[i]^[j + FFTN] := -xarray^[i]^[j + FFTN];
      end;
end;

// Convert a REAL/COMPLEX image into an integer one for display
// ch: channel to fill in output (0=B, 1=G, 2=R,  3=RGB)

procedure TIEFtImage.filt_toint(oimage: PIECOMPLEX_IMAGE; output: TIEBitmap; ch: integer);
var
  xmax, xmin, x, xdif: single;
  i, j, n: integer;
  xx: byte;
  rgb: pbytearray;
  pix: pbyte;
  image: PIECOMPLEX_IMAGE;
begin
  image := dupcomplex(oimage);
  n := FFTN;
  xmax := -1.0E20;
  xmin := -xmax;
  for i := 0 to n - 1 do
  begin
    for j := 0 to n - 1 do
    begin
      x := sqrt(pix_cnorm(image^[i]^[j], image^[i]^[j + n]));
      if (x > 0.0) then
        x := ln(sqrt(x))
      else
        x := 0.0;