Très utile en traitement d'image et analyse d'image
Utiliser avec la classe CxImage que vous pouvez trouver au :
http://www.xdp.it/
Source / Exemple :
bool CxImage::FFT2(CxImage* srcReal, CxImage* srcImag, CxImage* dstReal, CxImage* dstImag,
long direction, bool bForceFFT, bool bMagnitude)
{
//check if there is something to convert
if (srcReal==NULL && srcImag==NULL) return false;
long w,h;
//get width and height
if (srcReal) {
w=srcReal->GetWidth();
h=srcReal->GetHeight();
} else {
w=srcImag->GetWidth();
h=srcImag->GetHeight();
}
bool bXpow2 = IsPowerof2(w);
bool bYpow2 = IsPowerof2(h);
//if bForceFFT, width AND height must be powers of 2
if (bForceFFT && !(bXpow2 && bYpow2)) {
long i;
i=0;
while((1<<i)<w) i++;
w=1<<i;
bXpow2=true;
i=0;
while((1<<i)<h) i++;
h=1<<i;
bYpow2=true;
}
// I/O images for FFT
CxImage *tmpReal,*tmpImag;
// select output
tmpReal = (dstReal) ? dstReal : srcReal;
tmpImag = (dstImag) ? dstImag : srcImag;
// src!=dst -> copy the image
if (srcReal && dstReal) tmpReal->Copy(*srcReal,true,false,false);
if (srcImag && dstImag) tmpImag->Copy(*srcImag,true,false,false);
// dst&&src are empty -> create new one, else turn to GrayScale
if (srcReal==0 && dstReal==0){
tmpReal = new CxImage(w,h,8);
tmpReal->Clear(0);
tmpReal->SetGrayPalette();
} else {
if (!tmpReal->IsGrayScale()) tmpReal->GrayScale();
}
if (srcImag==0 && dstImag==0){
tmpImag = new CxImage(w,h,8);
tmpImag->Clear(0);
tmpImag->SetGrayPalette();
} else {
if (!tmpImag->IsGrayScale()) tmpImag->GrayScale();
}
if (!(tmpReal->IsValid() && tmpImag->IsValid())){
if (srcReal==0 && dstReal==0) delete tmpReal;
if (srcImag==0 && dstImag==0) delete tmpImag;
return false;
}
//resample for FFT, if necessary
tmpReal->Resample(w,h,0);
tmpImag->Resample(w,h,0);
//ok, here we have 2 (w x h), grayscale images ready for a FFT
double* real;
double* imag;
long j,k,m;
_complex **grid;
//double mean = tmpReal->Mean();
/* Allocate memory for the grid */
grid = (_complex **)malloc(w * sizeof(_complex));
for (k=0;k<w;k++) {
grid[k] = (_complex *)malloc(h * sizeof(_complex));
}
for (j=0;j<h;j++) {
for (k=0;k<w;k++) {
grid[k][j].x = tmpReal->GetPixelIndex(k,j)-128;
grid[k][j].y = tmpImag->GetPixelIndex(k,j)-128;
}
}
//DFT buffers
double *real2,*imag2;
real2 = (double*)malloc(max(w,h) * sizeof(double));
imag2 = (double*)malloc(max(w,h) * sizeof(double));
/* Transform the rows */
real = (double *)malloc(w * sizeof(double));
imag = (double *)malloc(w * sizeof(double));
m=0;
while((1<<m)<w) m++;
for (j=0;j<h;j++) {
for (k=0;k<w;k++) {
real[k] = grid[k][j].x;
imag[k] = grid[k][j].y;
}
if (bXpow2) FFT(direction,m,real,imag);
else DFT(direction,w,real,imag,real2,imag2);
for (k=0;k<w;k++) {
grid[k][j].x = real[k];
grid[k][j].y = imag[k];
}
}
free(real);
free(imag);
/* Transform the columns */
real = (double *)malloc(h * sizeof(double));
imag = (double *)malloc(h * sizeof(double));
m=0;
while((1<<m)<h) m++;
for (k=0;k<w;k++) {
for (j=0;j<h;j++) {
real[j] = grid[k][j].x;
imag[j] = grid[k][j].y;
}
if (bYpow2) FFT(direction,m,real,imag);
else DFT(direction,h,real,imag,real2,imag2);
for (j=0;j<h;j++) {
grid[k][j].x = real[j];
grid[k][j].y = imag[j];
}
}
free(real);
free(imag);
free(real2);
free(imag2);
/* converting from double to byte, there is a HUGE loss in the dynamics
"nn" tries to keep an acceptable SNR, but 8bit=48dB: don't ask more */
double nn=pow((double)2,(double)log((double)max(w,h))/(double)log((double)2)-4);
//reversed gain for reversed transform
if (direction==-1) nn=1/nn;
//bMagnitude : just to see it on the screen
if (bMagnitude) nn*=4;
for (j=0;j<h;j++) {
for (k=0;k<w;k++) {
if (bMagnitude){
tmpReal->SetPixelIndex(k,j,(BYTE)max(0,min(255,(nn*(3+log(_cabs(grid[k][j])))))));
if (grid[k][j].x==0){
tmpImag->SetPixelIndex(k,j,(BYTE)max(0,min(255,(128+(atan(grid[k][j].y/0.0000000001)*nn)))));
} else {
tmpImag->SetPixelIndex(k,j,(BYTE)max(0,min(255,(128+(atan(grid[k][j].y/grid[k][j].x)*nn)))));
}
} else {
tmpReal->SetPixelIndex(k,j,(BYTE)max(0,min(255,(128 + grid[k][j].x*nn))));
tmpImag->SetPixelIndex(k,j,(BYTE)max(0,min(255,(128 + grid[k][j].y*nn))));
}
}
}
for (k=0;k<w;k++) free (grid[k]);
free (grid);
if (srcReal==0 && dstReal==0) delete tmpReal;
if (srcImag==0 && dstImag==0) delete tmpImag;
return true;
}
Vous n'êtes pas encore membre ?
inscrivez-vous, c'est gratuit et ça prend moins d'une minute !
Les membres obtiennent plus de réponses que les utilisateurs anonymes.
Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.
Le fait d'être membre vous permet d'avoir des options supplémentaires.