TRANSFORMÉ DE FOURIER RAPIDE EN TRAITEMENT D'IMAGE

//---------------------------------------------------------------------------

#include <vcl.h>
#pragma hdrstop
#include<time.h>
#include "Unit1.h"
#include "Unit2.h"
#define pi M_PI
///-------------------------
#include <stdio.h>
#include <math.h>
#include <complex>
#include <iostream>
using namespace std;
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
double **r,**g,**b;
double **r1;
unsigned long debut, fin;
complex<long double> *line, *rline, *col, *rcol,**y, **RES;
long W,H;
bool charger=false;
unsigned char C_MAX=0;

//------------------------------

void remplir()
{
long double t;
for(int i=0;i<W;i++)
for(int j=0;j<H;j++)
{
 r[i][j]=GetRValue(Form1->Image1->Canvas->Pixels[i][j]);
 g[i][j]=GetGValue(Form1->Image1->Canvas->Pixels[i][j]);
 b[i][j]=GetBValue(Form1->Image1->Canvas->Pixels[i][j]);

 t=(0.2125*r[i][j]+0.7154*g[i][j]+0.0721*b[i][j]);
 if(t>C_MAX) C_MAX=t;
 y[i][j]=(t);
}

}
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
        : TForm(Owner)
{
}
//---------------------------------------------------------------------------

void __fastcall TForm1::Ouvrir1Click(TObject *Sender)
{
if (OpenPictureDialog1->Execute())
 {
  charger=true;
  Image2->Picture=NULL;
  Image1->Picture->LoadFromFile(OpenPictureDialog1->FileName);
  Image2->Width=Image1->Width; W=Image1->Width;
  Image2->Height=Image1->Height; H=Image1->Height;
//***********************************
  r=new double *[Form1->Image1->Width];
  for (int i=0;i<Form1->Image1->Width;i++)
  r[i]=new double [Form1->Image1->Height];
  g=new double *[Form1->Image1->Width];
  for (int i=0;i<Form1->Image1->Width;i++)
  g[i]=new double [Form1->Image1->Height];
  b=new double *[Form1->Image1->Width];
  for (int i=0;i<Form1->Image1->Width;i++)
  b[i]=new double [Form1->Image1->Height];
//********************************************
  y=new complex<long double> *[W];
  for (int i=0;i<W;i++) y[i]=new complex<long double> [H];

  RES= new  complex<long double> *[W];
  for (int i=0;i<W;i++) RES[i]=new complex<long double> [H];
  line=new complex<long double> [W];
  rline=new complex<long double> [W];
  col=new complex<long double> [H];
  rcol=new complex<long double> [H];
//------------------------
  remplir();
 }
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Fermer1Click(TObject *Sender)
{
Application->Terminate();
}
//---------------------------------------------------------------------------

void fft(int n, complex<long double> X[], complex<long double> RES[])
{
complex<long double> *pair, *impair, *U, *V;
int i, i1;
long double alphai, alpha;

if(n==1) RES[0]=X[0];  // fin de la pile recursive
else{

     pair = new complex<long double> [n/2];
     impair = new complex<long double> [n/2];
     // U , V deux vecteurs RES pour les appels recursive
     U = new complex<long double> [n/2];
     V = new complex<long double> [n/2];
     // remplissage des vecteurs paire et impaire
     for(i=0; i<n/2; ++i){
        pair[i]=X[2*i];
        impair[i]=X[2*i+1];
     }
     // appel recursive pour paire avec U (resp: impair avec V)
     fft( n/2, pair, U);
     fft( n/2, impair, V);
     // Calcul de la FFT
     alpha= 2*pi/n;
     for(i=0; i<n-1; ++i){
        alphai= alpha*i;
        i1=i%(n/2);
        complex<long double> tau (cosl(alphai), sinl(alphai));
        RES[i]=U[i1] + tau * V[i1];
     }

     delete pair;
     delete impair;
     delete U;
     delete V;
}
}

void __fastcall TForm1::FFT1Click(TObject *Sender)
{
int i,j;
long double C;
double g,max=0;

if (charger==true)

debut=clock();

// FFT 1D sur les lignes

for(j=0; j<H; j++){
for(i=0;i<W; i++) line[i]=y[i][j];  // extraction de la ligne

fft(W, line, rline);    // FFT sur la ligne
for(i=0;i<W; i++) RES[i][j]=rline[i]; //Ranger la ligne dans la matrice resultat
}


// FFT 1D sur les colonnes de la matrice resultat

for(i=0; i<W; i++){
for(j=0;j<H; j++) col[j]=RES[i][j];  // extraction de la colonne

fft(H, col, rcol);    // FFT sur la colonne
for(j=0;j<H; j++) RES[i][j]=rcol[j]; //Ranger la colonne dans la matrice final
}

C=255/log(1+C_MAX);

for(i=0; i<W/2; i++)
  for(j=0;j<H/2; j++){
  g=log(abs(RES[W/2+i][H/2+j])+1);
  if (g>max) max=g;
  g=log(abs(RES[i][j])+1);
    if (g>max) max=g;
}

for(i=W/2; i<W; i++)
 for(j=0;j<H/2; j++){
 g=C*log(abs(RES[i-W/2][H/2+j])+1);
 if (g>max) max=g;
  g=C*log(abs(RES[i][j])+1);
 if (g>max) max=g;
}
for(i=0; i<W/2; i++)
  for(j=0;j<H/2; j++){
  g=C*log(abs(RES[W/2+i][H/2+j])+1);
  Form1->Image2->Canvas->Pixels[i][j]=RGB((g*255)/float(max),(g*255)/float(max),(g*255)/float(max));
  g=C*log(abs(RES[i][j])+1);
  Form1->Image2->Canvas->Pixels[W/2+i][H/2+j]=RGB((g*255)/float(max),(g*255)/float(max),(g*255)/float(max));
}

for(i=W/2; i<W; i++)
  for(j=0;j<H/2; j++){
  g=C*log(abs(RES[i-W/2][H/2+j])+1);
  Form1->Image2->Canvas->Pixels[i][j]=RGB((g*255)/float(max),(g*255)/float(max),(g*255)/float(max));
  g=C*log(abs(RES[i][j])+1);
  Form1->Image2->Canvas->Pixels[i-W/2][H/2+j]=RGB((g*255)/float(max),(g*255)/float(max),(g*255)/float(max));
}
fin=clock();
ShowMessage("Temps :"+FloatToStr(((float)(fin-debut)/1000))+" sec!");
}
//---------------------------------------------------------------------------