iefft.pas

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

      if (x > xmax) then
        xmax := x;
      if (x < xmin) then
        xmin := x;
      image^[i]^[j] := x;
      image^[i]^[j + FFTN] := 0.0;
    end;
  end;
  if ((xmin <= 0.00001) and (xmin >= -0.00001)) then
    xmin := 0.000000;
  if (xmax - xmin) = 0 then
    xdif := 0
  else
    xdif := 255 / (xmax - xmin);
  for i := 0 to n - 1 do
  begin
    if ch < 3 then
    begin
      // fill R, G or B
      rgb := output.ScanLine[i];
      for j := 0 to n - 1 do
      begin
        rgb^[ch] := blimit(trunc((image^[i]^[j] - xmin) * xdif));
        inc(prgb(rgb));
      end
    end
    else
    begin
      // fill RGB
      pix := output.ScanLine[i];
      for j := 0 to n - 1 do
      begin
        xx := blimit(trunc((image^[i]^[j] - xmin) * xdif));
        pix^ := xx;
        inc(pix);
        pix^ := xx;
        inc(pix);
        pix^ := xx;
        inc(pix);
      end;
    end;
  end;
  freecomplex(image);
end;

// Scale a floating point image to the range 0-255

procedure TIEFtImage.realtoint(fim: PIECOMPLEX_IMAGE; H: plongintarray);
var
  data: PIEsinglearray;
  dij: pdwordarray;
  i, j, k: integer;
  ii, jj: dword;
  flast, xmax, xmin, xd: single;
  nr, nc: integer;
  N: longint;
begin
  flast := 0.0;
  nr := FFTN;
  nc := FFTN;
  N := nr * nc;
  if (H = nil) then
  begin
    xmin := abs(fim^[0]^[0]);
    xmax := xmin;
    for i := 0 to nr - 1 do
      for j := 0 to nc - 1 do
      begin
        fim^[i]^[j] := abs(fim^[i]^[j]);
        if (xmax < fim^[i]^[j]) then
          xmax := fim^[i]^[j];
        if (xmin > fim^[i]^[j]) then
          xmin := fim^[i]^[j];
      end;
    if (xmax - xmin) = 0 then
      xd := 0
    else
      xd := 255 / (xmax - xmin);
    for i := 0 to nr - 1 do
      for j := 0 to nc - 1 do
        fim^[i]^[j] := (fim^[i]^[j] - xmin) * xd;
  end
  else
  begin
    k := 0;
    getmem(data, N * sizeof(single));
    getmem(dij, N * sizeof(dword));
    for i := 0 to nr - 1 do
      for j := 0 to nc - 1 do
      begin
        data[k] := abs(fim^[i]^[j]);
        dij[k] := (i shl 10) or j;
        inc(k);
      end;
    pairsort(data, dij, N);
    j := 0;
    k := 0;
    for i := 0 to 255 do
    begin
      if (k <= 0) then
        k := 0;
      while (k < H[i]) and (j < nr * nc) do
      begin
        ii := dij[j];
        jj := ii and 1023;
        ii := (ii shr 10) and 1023;
        flast := data[j];
        fim[ii][jj] := i;
        inc(k);
        inc(j);
      end;
      k := 0;
      while (data[j] = flast) do
      begin
        ii := dij[j];
        jj := ii mod 1024;
        ii := trunc(ii / 1024);
        fim[ii][jj] := i;
        inc(j);
        inc(k);
      end;
    end;
    freemem(data);
    freemem(dij);
    for i := 0 to nr - 1 do
      for j := 0 to nc - 1 do
        fim^[i]^[j] := abs(fim^[i]^[j]);
  end;
end;

(*
function vlog2(x:integer):integer;
var
 k:integer;
begin
 result := 0;
   k := 1;
 while (k<x) do begin
  k := k*2;
      inc(result);
 end;
end;
*)

function IntPower(Base: Extended; Exponent: Integer): Extended;
asm
        mov     ecx, eax
        cdq
        fld1                      { Result := 1 }
        xor     eax, edx
        sub     eax, edx          { eax := Abs(Exponent) }
        jz      @@3
        fld     Base
        jmp     @@2
@@1:    fmul    ST, ST            { X := Base * Base }
@@2:    shr     eax,1
        jnc     @@1
        fmul    ST(1),ST          { Result := Result * X }
        jnz     @@1
        fstp    st                { pop X from FPU stack }
        cmp     ecx, 0
        jge     @@3
        fld1
        fdivrp                    { Result := 1 / Result }
@@3:
        fwait
end;

function pow(Base, Exponent: double): double;
begin
  if Exponent = 0.0 then
    Result := 1.0 { n**0 = 1 }
  else if (Base = 0.0) and (Exponent > 0.0) then
    Result := 0.0 { 0**n = 0, n > 0 }
  else if (Frac(Exponent) = 0.0) and (Abs(Exponent) <= MaxInt) then
    Result := IntPower(Base, Trunc(Exponent))
  else
    Result := Exp(Exponent * Ln(Base))
end;

procedure TIEFtImage.fft(data: PIEsinglearray; dir: single);
var
  tempval: TIEComplex;
  i: integer;
  ba: PIEinteger;
begin
  direction := dir;
  if (direction <> 1.0) then
    direction := -1.0;
  _fft(data, 1, 0);
  ba := bittabpt;
  for i := 0 to numpts - 1 do
  begin
    if (ba^ <= i) then
      with tempval do
      begin
        real := data^[i];
        imag := data^[i + FFTN];
        data^[i] := data[ba^];
        data^[i + FFTN] := data[ba^ + FFTN];
        data^[ba^] := real;
        data^[ba^ + FFTN] := imag;
      end;
    inc(ba);
  end;
  for i := 0 to numpts - 1 do
  begin
    data^[i] := data^[i] * scalef;
    data^[i + FFTN] := data^[i + FFTN] * scalef;
  end;
end;

//  _fft -- A recursive FFT routine.

procedure TIEFtImage._fft(tseries: PIEsinglearray; level: integer; chunk: integer);
var
  nodes, i: integer;
  sinindx, cosindx: integer;
  dual1, dual2: integer;
  dual1val, dual2val, dualprod, wp: TIEComplex;
begin
  nodes := (PIEinteger(integer(powers) + (nn - level) * 4))^;
  sinindx := (PIEinteger(integer(bittabpt) + (chunk div nodes) * 4))^;
  cosindx := (sinindx + numpts div 4) mod numpts;
  wp.real := (PIEdouble(integer(sintabpt) + cosindx * 8))^;
  wp.imag := direction * (PIEdouble(integer(sintabpt) + sinindx * 8))^;
  for i := 0 to nodes - 1 do
  begin
    dual1 := chunk + i;
    dual2 := dual1 + nodes;
    with dual2val do
    begin
      real := tseries^[dual2];
      imag := tseries^[dual2 + FFTN];
      dualprod.real := real * wp.real - imag * wp.imag;
      dualprod.imag := real * wp.imag + imag * wp.real;
    end;
    with dual1val do
    begin
      real := tseries^[dual1];
      imag := tseries^[dual1 + FFTN];
      tseries^[dual1] := real + dualprod.real;
      tseries^[dual1 + FFTN] := imag + dualprod.imag;
      tseries^[dual2] := real - dualprod.real;
      tseries^[dual2 + FFTN] := imag - dualprod.imag;
    end;
  end;
  if (level < nn) then
  begin
    _fft(tseries, level + 1, chunk);
    _fft(tseries, level + 1, chunk + nodes);
  end;
end;

procedure TIEFtImage.fftinit(nopts: integer);
var
  i: integer;
  po: PIEinteger;
  si: PIEdouble;
  bi: PIEinteger;
begin
  numpts := nopts;
  nn := trunc((ln(numpts) / ln(2.0)) + 0.5);
  scalef := 1.0 / sqrt(numpts);
  if (sintabpt <> nil) then
    freemem(sintabpt);
  if (bittabpt <> nil) then
    freemem(bittabpt);
  if (powers <> nil) then
    freemem(powers);
  getmem(sintabpt, nopts * sizeof(double));
  getmem(bittabpt, nopts * sizeof(integer));
  getmem(powers, (nn + 1) * sizeof(integer));
  po := powers;
  for i := 0 to nn do
  begin
    po^ := trunc(pow(2.0, i) + 0.5);
    inc(po);
  end;
  si := sintabpt;
  bi := bittabpt;
  for i := 0 to numpts - 1 do
  begin
    si^ := sin(6.283185307179587 * i / numpts);
    bi^ := bitrev(i);
    inc(si);
    inc(bi);
  end;
end;

function TIEFtImage.bitrev(bits: integer): integer;
var
  i, lookmask, setmask, tempbit: integer;
begin
  lookmask := 1;
  setmask := numpts;
  setmask := setmask shr 1;
  tempbit := 0;
  for i := 0 to nn - 1 do
  begin
    if ((bits and lookmask) = lookmask) then
      tempbit := tempbit or setmask;
    lookmask := lookmask shl 1;
    setmask := setmask shr 1;
  end;
  result := tempbit;
end;

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

function TIEFtImage.newcomplex(im: TIEBitmap; ch: integer): PIECOMPLEX_IMAGE;
var
  i, j: integer;
  x: PIECOMPLEX_IMAGE;
  y: PIEsinglearray;
  xmax: single;
  pix: pbytearray;
begin
  getmem(x, (im.Height * sizeof(PIEsinglearray)));
  getmem(y, (sizeof(single) * im.Height * im.Width * 2));
  for i := 0 to im.Height - 1 do
    x^[i] := @(y^[i * FFTN * 2]);
  xmax := 1;
  for i := 0 to im.Height - 1 do
  begin
    pix := im.Scanline[i];
    if ch = 3 then
    begin
      // convert to grayscale
      for j := 0 to im.Width - 1 do
      begin
        with PRGB(pix)^ do
          x^[i]^[j] := (r * gRedToGrayCoef + g * gGreenToGrayCoef + b * gBlueToGrayCoef) div 100;
        x^[i]^[j + im.Width] := 0.0;
        if (x^[i]^[j] > xmax) then
          xmax := x^[i]^[j];
        inc(PRGB(pix));
      end;
    end
    else
    begin
      // rgb
      for j := 0 to im.Width - 1 do
      begin
        x^[i]^[j] := pix^[ch];
        x^[i]^[j + im.Width] := 0.0;
        if (x^[i]^[j] > xmax) then
          xmax := x^[i]^[j];
        inc(pbyte(pix), 3);
      end;
    end;
  end;
  if NORMALIZE then
    for i := 0 to im.Height - 1 do
      for j := 0 to im.Width - 1 do
        x^[i]^[j] := x^[i]^[j] / xmax;
  result := x;
end;

function TIEFtImage.dupcomplex(im: PIECOMPLEX_IMAGE): PIECOMPLEX_IMAGE;
var
  i: integer;
  x: PIECOMPLEX_IMAGE;
  y: PIEsinglearray;
begin
  getmem(x, (FFTN * sizeof(PIEsinglearray)));
  getmem(y, (sizeof(single) * FFTN * FFTN * 2));
  x^[0] := y;
  for i := 1 to FFTN - 1 do
    x^[i] := @(y^[i * 2 * FFTN]);
  for i := 0 to FFTN - 1 do
    move(im^[i]^[0], x^[i]^[0], 2 * FFTN * sizeof(single));
  result := x;
end;

procedure TIEFtImage.pairsort(arr: PIEsinglearray; iarr: pdwordarray; n: integer);
begin
  fqsort(arr, iarr, 0, n - 1);
end;

procedure TIEFtImage.fqsort(arr: PIEsinglearray; iarr: pdwordarray; l: integer; r: integer);
var
  i, j: integer;
  k: dword;
  x, w: single;
begin
  i := l;
  j := r;
  x := arr^[(l + r) div 2];
  repeat
    while (arr^[i] < x) do
      inc(i);
    while (x < arr^[j]) do
      dec(j);
    if (i <= j) then
    begin
      w := arr^[i];
      arr^[i] := arr^[j];
      arr^[j] := w;
      k := iarr^[i];
      iarr^[i] := iarr^[j];
      iarr^[j] := k;
      inc(i);
      dec(j);
    end;
  until not (i <= j);
  if (l < j) then
    fqsort(arr, iarr, l, j);
  if (i < r) then
    fqsort(arr, iarr, i, r);
end;

function TIEFtImage.GetComplexImage(x, y: integer): TIEComplexColor;
begin
  zeromemory(@result, sizeof(TIEComplexColor));
  with result do
  begin
    if assigned(fftr) then
    begin
      real_Red := @fftr^[y]^[x];
      imag_Red := @fftr^[y]^[x + FFTN];
    end;
    if assigned(fftg) then
    begin
      real_Green := @fftg^[y]^[x];
      imag_Green := @fftg^[y]^[x + FFTN];
    end;
    if assigned(fftb) then
    begin
      real_Blue := @fftb^[y]^[x];
      imag_Blue := @fftb^[y]^[x + FFTN];
    end;
    if assigned(fftgray) then
    begin
      real_Gray := @fftgray^[y]^[x];
      imag_Gray := @fftgray^[y]^[x + FFTN];
    end;
  end;
end;

procedure TIEFtImage.HiPass(radius: integer);
var
  cc, i, j, dist, c1: integer;
begin
  cc := FFTN div 2;
  for i := 0 to FFTN - 1 do
  begin
    c1 := (i - cc) * (i - cc);
    for j := 0 to FFTN - 1 do
    begin
      dist := trunc(sqrt((c1 + (j - cc) * (j - cc))));
      if (dist <= radius) then
        with ComplexPixel[j, i] do
          case fImageType of
            ieitRGB:
              begin
                real_Red^ := 0;
                imag_Red^ := 0;
                real_Green^ := 0;
                imag_Green^ := 0;
                real_Blue^ := 0;
                imag_Blue^ := 0;
              end;
            ieitGrayscale:
              begin
                real_Gray^ := 0;
                imag_Gray^ := 0;
              end;
          end;
    end;
  end;
end;

procedure TIEFtImage.LoPass(radius: integer);
var
  cc, i, j, dist, c1: integer;
begin
  cc := FFTN div 2;
  for i := 0 to FFTN - 1 do
  begin
    c1 := (i - cc) * (i - cc);
    for j := 0 to FFTN - 1 do
    begin
      dist := trunc(sqrt((c1 + (j - cc) * (j - cc))));
      if (dist > radius) then
        with ComplexPixel[j, i] do
          case fImageType of
            ieitRGB:
              begin
                real_Red^ := 0;
                imag_Red^ := 0;
                real_Green^ := 0;
                imag_Green^ := 0;
                real_Blue^ := 0;
                imag_Blue^ := 0;
              end;
            ieitGrayscale:
              begin
                real_Gray^ := 0;
                imag_Gray^ := 0;
              end;
          end;
    end;
  end;
end;

// set to 0 the complex region x1,y1...x2,y2

procedure TIEFtImage.ClearZone(x1, y1, x2, y2: integer);
var
  x, y: integer;
begin
  for x := x1 to x2 do
    for y := y1 to y2 do
      with ComplexPixel[x, y] do
        case fImageType of
          ieitRGB:
            begin
              real_Red^ := 0;
              imag_Red^ := 0;
              real_Green^ := 0;
              imag_Green^ := 0;
              real_Blue^ := 0;
              imag_Blue^ := 0;
            end;
          ieitGrayscale:
            begin
              real_Gray^ := 0;
              imag_Gray^ := 0;
            end;
        end;
end;

end.