fft_2d.pas

灰度图象处理的小例子,例子中包括图象的一些处理,很有用的!

//  Transmap在 主程序中定义为TBitmap, 用之前要Create,  并将原图装入。

unit fft_2d;

interface

uses Windows, SysUtils, Classes, Graphics, Forms, Controls, Menus,
  StdCtrls, Dialogs, Buttons, Messages, ExtCtrls, ComCtrls,math,global,
  ExtDlgs, Gauges;

procedure fft2d(map:Tbitmap);
procedure ifft2d(xxx:integer);

implementation

uses Main;
var
   sequencex,sequencey:parrayofint;
   bitnumx,bitnumy:longint;

procedure fft2d(map:Tbitmap);
var
  lengthx,lengthy,length,x,y,position:longint;
  temp:Parrayofdouble;
  i,j,k,n,bak,init,step:integer;
  tmp:byte;
  mc:PByteArray;
  ltmp,tmp1_r,tmp1_i,tmp2_r,tmp2_i,tcos,tsin,max:double;
begin
   if TRANSBEGIN_FFT then
     begin
       freemem(frq_r);freemem(frq_i);
     end;
   lengthx:=map.Height;
   if log2(lengthx)>round(log2(lengthx)) then
        bitnumx:=round(log2(lengthx))+1
      else
        bitnumx:=round(log2(lengthx));
   lengthy:=map.Width;
   if log2(lengthy)>round(log2(lengthy)) then
        bitnumy:=round(log2(lengthy))+1
      else
        bitnumy:=round(log2(lengthy));

   lengthx:=round(power(2,bitnumx));
   lengthy:=round(power(2,bitnumy));
   FFT_Width:=lengthy;
   FFT_Height:=lengthx;
   length:=lengthx*lengthy;
//-----------------------------------------
   frq_r:=allocmem(length*sizeof(double));
   frq_i:=allocmem(length*sizeof(double));
   temp:=allocmem(length*sizeof(double));
   sequencex:=allocmem(lengthx*sizeof(longint));
   sequencey:=allocmem(lengthy*sizeof(longint));
   endi:=lengthx-1;endj:=bitnumx;
    for i:=0 to endi do
        begin
          bak:=i;
          sequencex[i]:=0;
          for j:=1 to endj do
            begin
             sequencex[i]:=sequencex[i]+(bak-(bak div 2)*2)*round(power(2.0,bitnumx-j));
             bak:=bak div 2;
            end;
        end;
   endi:=lengthy-1;
   endj:=bitnumy;
    for i:=0 to endi do
        begin
          bak:=i;
          sequencey[i]:=0;
          for j:=1 to endj do
            begin
             sequencey[i]:=sequencey[i]+(bak-(bak div 2)*2)*round(power(2.0,bitnumy-j));
             bak:=bak div 2;
            end;
        end;

//----------------------------------------
   endi:=map.height-1;
   endj:=map.width-1;
   for i:=0 to endi do
     begin
      mc:=sourcemap.scanline[i];
      for j:=0 to endj do
        begin
           tmp:=mc[j];
           frq_i[i*map.width+j]:=tmp;
        end;
      endj:=lengthy-1;
      for j:=map.Width to endj do
        begin
           frq_i[i*map.width+j]:=0;
        end;
     end;
   endi:=lengthy-1;endj:=lengthx-1;
   for i:=map.Height to endi do
      for j:=0 to endj do
           frq_i[i*map.width+j]:=0;
//-----------------------------------
   endi:=map.height-1;endj:=lengthy-1;
   for i:=0 to endi do
    begin
    for j:=0 to endj do
        frq_r[i*lengthy+j]:=frq_i[i*lengthy+sequencey[j]];
    for j:=0 to endj do
        frq_i[i*lengthy+j]:=0;
    end;
//------------------------------------
//    showmessage('begin!');
    endx:=map.Height-1;
    for x:=0 to endx do
     begin
     //formmain.gauge1.Progress:=50 *x div endx;
     n:=2;while(n <= lengthy) do
        begin
          init:=0;while init<lengthy do
             begin
               step:=n div 2;
               i:=init;
               while i< n/2+init do
                   begin
                     position:=x*lengthy+i;
                     tcos:=cos(-2*pi/n*i);
                     tsin:=sin(-2*pi/n*i);
                     tmp2_r:=frq_r[position+step]*tcos-frq_i[position+step]*tsin;
                     tmp2_i:=frq_r[position+step]*tsin+frq_i[position+step]*tcos;
                     tmp1_r:=frq_r[position];
                     tmp1_i:=frq_i[position];
                     frq_i[position]:=tmp1_i+tmp2_i;
                     frq_r[position]:=tmp1_r+tmp2_r;
                     frq_i[position+step]:=tmp1_i-tmp2_i;
                     frq_r[position+step]:=tmp1_r-tmp2_r;
                   i:=i+1;end;
          init:=init+n; end;
     n:=n*2;end;
   end;
//------------------------------------------
   endi:=lengthy-1;endj:=lengthx-1;
   for i:=0 to endi do
    for j:=0 to endj do
        temp[i*lengthy+j]:=frq_r[i+lengthy*sequencex[j]];
   endi:=length-1;
   for i:=0 to endi do
        frq_r[i]:=temp[i];
   endi:=lengthy-1;endj:=lengthx-1;
   for i:=0 to endi do
    for j:=0 to endj do
        temp[i*lengthy+j]:=frq_i[i+lengthy*sequencex[j]];
   endi:=length-1;
   for i:=0 to endi do
        frq_i[i]:=temp[i];
//------------------------------------------
    endi:=lengthy-1;
    for x:=0 to endi do
    begin
    //mainform.gauge1.Progress:=mainform.gauge1.Progress+50 * x div endi;
    n:=2;while(n <= lengthx) do
        begin
          init:=0;while init<lengthx do
             begin
               step:=n div 2;
               i:=init;while i< n/2+init do
                   begin
                     position:=x*lengthx+i;
                     tcos:=cos(-2*pi/n*i);
                     tsin:=sin(-2*pi/n*i);
                     tmp2_r:=frq_r[position+step]*tcos-frq_i[position+step]*tsin;
                     tmp2_i:=frq_r[position+step]*tsin+frq_i[position+step]*tcos;
                     tmp1_r:=frq_r[position];
                     tmp1_i:=frq_i[position];
                     frq_i[position]:=tmp1_i+tmp2_i;
                     frq_r[position]:=tmp1_r+tmp2_r;
                     frq_i[position+step]:=tmp1_i-tmp2_i;
                     frq_r[position+step]:=tmp1_r-tmp2_r;
                     i:=i+1;
                   end;
          init:=init+n; end;
         n:=n*2;
        end;
      end;

//------------------------------------------
// showmessage('END!');
 max:=0;
 endi:=length-1;
 tmp1_r:=sqrt(length);
 for i:=0 to endi do
   begin
      temp[i]:=round(sqrt(sqr(frq_i[i])+sqr(frq_r[i])))/tmp1_r;
     if max<temp[i] then
        max:=temp[i];
      frq_i[i]:=frq_i[i]/tmp1_r;
      frq_r[i]:=frq_r[i]/tmp1_r;
   end;

 transmap.height:=lengthx;
 transmap.width:=lengthy;
 endi:=lengthx-1;endj:=lengthy-1;
 for i:=0 to endi do
   begin
      mc:=transmap.scanline[i];
      for j:=0 to endj do
        begin
              ltmp:=sqrt(sqr(frq_i[j*map.height+i])+sqr(frq_r[j*map.Height+i]))/max*255*tmp1_r;
             tmp:=round(log2(ltmp+1)*32-1);
            mc[j]:=tmp;
        end;
   end;
 freemem(temp);
 //formmain.image1.picture.assign(transmap);
end;










procedure ifft2d(xxx:integer);
var
  lengthx,lengthy,length,x,y,position:longint;
  temp:Parrayofdouble;
  i,j,k,n,bak,init,step:integer;
  tmp:byte;
  mc:PByteArray;
  ltmp,tmp1_r,tmp1_i,tmp2_r,tmp2_i,tcos,tsin,max:double;
begin
   lengthx:=FFT_Height;
   lengthy:=FFT_Width;
   length:=lengthx*lengthy;

   temp:=allocmem(length*sizeof(double));
//------------------------------------------
   endi:=lengthy-1;endj:=lengthx-1;
   for i:=0 to endi do
    for j:=0 to endj do
        temp[i*lengthy+j]:=frq_r[i+lengthy*sequencey[j]];
   endi:=length-1;
   for i:=0 to endi do
        frq_r[i]:=temp[i];
   endi:=lengthy-1;endj:=lengthx-1;
   for i:=0 to endi do
    for j:=0 to endj do
        temp[i*lengthy+j]:=frq_i[i+lengthy*sequencey[j]];
   endi:=length-1;
   for i:=0 to endi do
        frq_i[i]:=temp[i];
//------------------------------------
//    showmessage('begin!');
    endx:=fft_Height-1;
    for x:=0 to endx do
     begin
     //mainform.gauge1.Progress:=50 *x div endx;
     n:=2;while(n <= lengthy) do
        begin
          init:=0;while init<lengthy do
             begin
               step:=n div 2;
               i:=init;
               while i< n/2+init do
                   begin
                     position:=x*lengthy+i;
                     tcos:=cos(-2*pi/n*i);
                     tsin:=sin(-2*pi/n*i);
                     tmp2_r:=frq_r[position+step]*tcos-frq_i[position+step]*tsin;
                     tmp2_i:=frq_r[position+step]*tsin+frq_i[position+step]*tcos;
                     tmp1_r:=frq_r[position];
                     tmp1_i:=frq_i[position];
                     frq_i[position]:=tmp1_i+tmp2_i;
                     frq_r[position]:=tmp1_r+tmp2_r;
                     frq_i[position+step]:=tmp1_i-tmp2_i;
                     frq_r[position+step]:=tmp1_r-tmp2_r;
                   i:=i+1;end;
          init:=init+n; end;
     n:=n*2;end;
   end;
//------------------------------------------
   endi:=lengthy-1;endj:=lengthx-1;
   for i:=0 to endi do
    for j:=0 to endj do
        temp[i*lengthy+j]:=frq_r[i+lengthy*sequencex[j]];
   endi:=length-1;
   for i:=0 to endi do
        frq_r[i]:=temp[i];
   endi:=lengthy-1;endj:=lengthx-1;
   for i:=0 to endi do
    for j:=0 to endj do
        temp[i*lengthy+j]:=frq_i[i+lengthy*sequencex[j]];
   endi:=length-1;
   for i:=0 to endi do
        frq_i[i]:=temp[i];
//------------------------------------------
    endi:=lengthy-1;
    for x:=0 to endi do
    begin
    //mainform.gauge1.Progress:=mainform.gauge1.Progress+50 * x div endi;
    n:=2;while(n <= lengthx) do
        begin
          init:=0;while init<lengthx do
             begin
               step:=n div 2;
               i:=init;while i< n/2+init do
                   begin
                     position:=x*lengthx+i;
                     tcos:=cos(-2*pi/n*i);
                     tsin:=sin(-2*pi/n*i);
                     tmp2_r:=frq_r[position+step]*tcos-frq_i[position+step]*tsin;
                     tmp2_i:=frq_r[position+step]*tsin+frq_i[position+step]*tcos;
                     tmp1_r:=frq_r[position];
                     tmp1_i:=frq_i[position];
                     frq_i[position]:=tmp1_i+tmp2_i;
                     frq_r[position]:=tmp1_r+tmp2_r;
                     frq_i[position+step]:=tmp1_i-tmp2_i;
                     frq_r[position+step]:=tmp1_r-tmp2_r;
                     i:=i+1;
                   end;
          init:=init+n; end;
         n:=n*2;
        end;
      end;

//------------------------------------------
// showmessage('END!');
 tmp1_r:=sqrt(length);

 transmap.height:=lengthx;transmap.width:=lengthy;
 endi:=lengthx-1;endj:=lengthy-1;
 for i:=0 to endi do
   begin
      mc:=transmap.scanline[i];
      for j:=0 to endj do
        begin
             ltmp:=sqrt(sqr(frq_i[j*fft_height+i])+sqr(frq_r[j*fft_Height+i]));
             tmp:=round(ltmp/tmp1_r);
{             tmp:=round(log2(ltmp+1)*32-1);
             ltmp:=tmp;
             tmp:=round(power(2,ltmp/32));}
            mc[j]:=tmp;
        end;
   end;
 freemem(temp);
 formmain.image1.picture.assign(transmap);
end;
end.