|
HPC/並列プログラミングポータルでは、HPC（High Performance Computing）プログラミングや並列プログラミングに関する情報を集積・発信しています。 |
［記事一覧を見る］
edit/* fft1.c */
#include stdio.h
#include string.h
#include stdlib.h
#include memory.h
#include math.h
#include time.h
#include "fftw3.h"
/*
fft1 input file output_file
*/
#define FMODE_WRITE "w"
#define NORM(x,y) (sqrt(x*x+y*y))
#define POWR(x,y) (x*x+y*y)
#define c_re(x) (x[0])
#define c_im(x) (x[1])
/* prototypes */
unsigned char* load_image( char* filename,
int* pWidth,
int* pHeight,
int* pDepth );
void write_image( char* filename,
unsigned char* buff,
int width,
int height,
int depth );
unsigned char* load_image( char* filename,
int* pWidth,
int* pHeight,
int* pDepth ) {
FILE* fp_in;
unsigned char* buff;
char szBuf[512];
int height = 0;
int width = 0;
int depth = 0;
fp_in = fopen( filename, "r" );
if( fp_in == NULL ) {
fprintf( stderr, "cannot open file: %s\n", filename );
return NULL;
}
/* read PGM Header */
fgets( szBuf, 512, fp_in );
if( strlen( szBuf ) 3 ) {
fprintf( stderr, "invalid image file: %s\n", filename );
fclose( fp_in );
return NULL;
}
if( strcmp( szBuf, "P5\n" ) != 0 ) {
fprintf( stderr, "invalid image file: %s\n", filename );
fclose( fp_in );
return NULL;
}
/* read #comment */
fgets( szBuf, 512, fp_in );
while( szBuf[0] == '#' ) {
while( (szBuf[strlen(szBuf)-1] != '\n') (!feof(fp_in)) ) {
fgets( szBuf, 512, fp_in );
}
fgets( szBuf, 512, fp_in );
}
if(feof(fp_in)) {
fprintf( stderr, "invalid image file: %s\n", filename );
fclose( fp_in );
return NULL;
}
/* read image size */
if (szBuf[strlen(szBuf)-1] != '\n') {
fprintf( stderr, "image file is strange: %s\n", filename );
fclose( fp_in );
return NULL;
}
sscanf( szBuf, "%d %d", width, height );
if( (width == 0) || (height == 0) ) {
fprintf( stderr, "image file is strange: %s\n", filename );
fclose( fp_in );
return NULL;
}
fgets( szBuf, 512, fp_in );
/* read depth */
if (szBuf[strlen(szBuf)-1] != '\n') {
fprintf( stderr, "image file is strange: %s\n", filename );
fclose( fp_in );
return NULL;
}
sscanf( szBuf, "%d", depth );
if( depth == 0 ) {
fprintf( stderr, "image file is strange: %s\n", filename );
fclose( fp_in );
return NULL;
}
printf( "image size: %d x %d, depth: %d\n", width, height, depth );
/* allocate buffer */
buff = (unsigned char*)malloc( width * height * sizeof(unsigned char) );
if( buff == NULL ) {
fclose( fp_in );
return NULL;
}
fread( buff, 1, width*height, fp_in );
fclose( fp_in );
*pWidth = width;
*pHeight = height;
*pDepth = depth;
return buff;
}
void write_image( char* filename,
unsigned char* buff,
int width,
int height,
int depth ) {
FILE* fp_out;
fp_out = fopen( filename, FMODE_WRITE );
if( fp_out == NULL ) {
fprintf( stderr, "cannot open file: %s\n", filename );
return;
}
/* write PGM Header */
fprintf( fp_out, "P5\n" );
/* write image size */
fprintf( fp_out, "%d %d\n", width, height );
/* write depth */
fprintf( fp_out, "%d\n", depth );
fwrite( buff, 1, width*height, fp_out );
if( ferror(fp_out) ) {
fprintf( stderr, "file output error: %s\n", filename );
return;
}
}
void FFTWExchange(int sizex,int sizey, fftw_complex *src, fftw_complex *dst) {
int i,j,u,v;
int imgxd2=sizex/2, imgyd2=sizey/2;
for(j=0;jsizey;j++){
for(i=0;isizex;i++){
u = (iimgxd2)? imgxd2+i: i-imgxd2;
v = (jimgyd2)? imgyd2+j: j-imgyd2;
c_re(dst[j*sizex+i]) = c_re(src[u+sizey*v]);
c_im(dst[j*sizex+i]) = c_im(src[u+sizey*v]);
}
}
}
unsigned char* execute_fft( unsigned char* input_image_buf,
int width,
int height,
int depth ) {
fftw_complex *in, *out;
fftw_plan p;
unsigned char* output_image_buf;
int n;
int i,j;
double dMax, dMin, dNorm;
double* pTmp;
dMax = 0.0;
dMin = 0.0;
pTmp = (double*)malloc( width*height*sizeof(double) );
if( pTmp == NULL ) {
return NULL;
}
/* allocate memory */
in = fftw_malloc(sizeof(fftw_complex) * width * height);
out = fftw_malloc(sizeof(fftw_complex) * width * height);
/* copy data to fftw buffer */
for( n = 0; n width*height; n++ ) {
in[n][0] = (float)input_image_buf[n];
in[n][1] = 0.0;
}
/* create plan */
p = fftw_plan_dft_2d(height, width, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
/* do FFT */
fftw_execute(p);
/* allocate buffer */
output_image_buf =
(unsigned char*)malloc( width * height * sizeof(unsigned char) );
if( output_image_buf != NULL ) {
/* copy data to fftw buffer */
for( n = 0; n width*height; n++ ) {
pTmp[n] = log(POWR(out[n][0],out[n][1]) + 1);
if( pTmp[n] dMin ) {
dMin = pTmp[n];
}
if( dMax pTmp[n] ) {
dMax = pTmp[n];
}
}
for( n = 0; n width*height; n++ ) {
double dWork;
dWork = (pTmp[n] * 255.0 / (dMax - dMin));
output_image_buf[n] = (unsigned char)(dWork);
}
}
/* free memories */
fftwnd_destroy_plan(p);
free( pTmp );
fftw_free(in);
fftw_free(out);
return output_image_buf;
}
unsigned char* execute_ifft( unsigned char* input_image_buf,
int width,
int height,
int depth ) {
fftw_complex *in, *out;
fftw_plan p;
unsigned char* output_image_buf;
int n;
/* allocate memory and create plan */
in = fftw_malloc(sizeof(fftw_complex) * width * height);
out = fftw_malloc(sizeof(fftw_complex) * width * height);
p = fftw_plan_dft_2d(height, width, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
/* allocate output buffer */
/* copy data to fftw buffer */
for( n = 0; n width*height; n++ ) {
in[n][0] = input_image_buf[n]; /* Re */
in[n][1] = 0; /* Im */
}
/* do FFT */
fftw_execute(p);
/* allocate buffer */
output_image_buf =
(unsigned char*)malloc( width * height * sizeof(unsigned char) );
if( output_image_buf != NULL ) {
/* copy data to fftw buffer */
for( n = 0; n width*height; n++ ) {
output_image_buf[n] = out[n][0];
}
}
/* free memories */
/*fftwnd_destroy_plan(p);*/
fftw_free(in);
fftw_free(out);
return output_image_buf;
}
int main( int argv, char** argc ) {
char* in_file_path;
char* out_file_path;
unsigned char* input_image_buf;
unsigned char* output_image_buf;
unsigned char* image_buf;
int width;
int height;
int depth;
clock_t clock_start, clock_end;
if( argv 3 ) {
printf( "%s input_file output_file\n", argc[0] );
return -1;
}
in_file_path = argc[1];
out_file_path = argc[2];
input_image_buf = load_image( in_file_path, width, height, depth );
if( input_image_buf == NULL ) {
fprintf( stderr, "cannot load image. exit.\n" );
return -1;
}
clock_start = clock();
output_image_buf =
execute_fft( input_image_buf, width, height, depth );
clock_end = clock();
if( output_image_buf != NULL ) {
write_image( out_file_path, output_image_buf, width, height, depth );
} else {
fprintf( stderr, "FFT operation faild.\n" );
}
free( input_image_buf );
free( output_image_buf );
printf( "start: %d end: %d elapse: %d.\n", clock_start, clock_end, clock_end - clock_start );
printf( "elapse time: %lf sec.", ((double)(clock_end - clock_start)) / (double)CLOCKS_PER_SEC );
/* free( image_buf ); */
return 0;
};
[PageInfo]
LastUpdate: 2009-11-18 20:35:35, ModifiedBy: hiromichi-m
[Permissions]
view:all, edit:login users, delete/config:members