2/5 (2 avis)
Snippet vu 15 177 fois - Téléchargée 32 fois
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; }
29 sept. 2011 à 08:35
26 févr. 2008 à 18:52
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.