２次元フーリエをするには？ (ID:52382)

名前

ホームページ(ブログ、Twitterなど)のURL (省略可)

本文

入力画像サイズが128×128（24bit）以上でも正しい結果が得られるようにするには、どのようにしたら良いでしょうか？256×256だとうまく結果がでません。
どなたか教えて頂けないでしょうか？？
下位は今の時点の作成プログラムです。
                        
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define PI 3.141592
#define OPT 1	//OPT=1 光学的DFT(直流分が中央)
				//OPT=0 通常のDFT(直流分が左端)

#define X_EXP 7     //X_EXP = log2(X_SIZE)
#define Y_EXP 7     //Y_EXP = log2(Y_SIZE)

#define ysize 128
#define xsize 128
unsigned char image_in[ysize][xsize][3]; 
unsigned char image_out[ysize][xsize][3];

/* BMPファイル用 */
typedef unsigned short WORD;
typedef unsigned long DWORD;

WORD  bfType;
DWORD bfSize;
WORD  bfReserved1,
      bfReserved2;
DWORD bfOffBits;
DWORD biSize, 
      biWidth, biHeight;
WORD  biPlanes,
      biBitCount;
DWORD biCompression, 
      biSizeImage,
      biXPelsPerMeter,
      biYPelsPerMeter,
      biClrUsed, 
      biClrImportant;

void readBMP(char *filename,unsigned char image[ysize][xsize][3]);
void writeBmp(char *filename,unsigned char image[ysize][xsize][3]);

void fft1core(double a_rl[], double a_im[], int length, int ex, double sin_tb1[], double cos_tb1[], double buf[]);
void cstb(int length, int inv, double sin_tb1[], double cos_tb1[]);
void birv(double a[], int length, int ex, double b[]);

int fft2(double a_rl[ysize][xsize], double a_im[ysize][xsize], int inv);
int fft(double a_rl[ysize][xsize], double a_im[ysize][xsize], int inv);

void rvmtx1(double a[ysize][xsize], double b[xsize][ysize]);
void rvmtx2(double a[xsize][ysize], double b[ysize][xsize]);

int fftimge(unsigned char image_in[ysize][xsize][3], unsigned char image_out[ysize][xsize][3]);
int fftfilter(unsigned char image_in[ysize][xsize][3], 
              unsigned char image_out[ysize][xsize][3],int a, int b);

int inv=1;

void main(void)
{
	
	char input[100],output[100];
	
	
	
	printf("入力画像ファイル名(input.bmp):");
	scanf("%s",input);
	printf("出力画像ファイル名(output.bmp):");
	scanf("%s",output);
	
	readBMP(input,image_in);

fftimge(image_in,image_out);
	writeBmp(output,image_out);

}

/*----------------------------------------------------------------------
                  BMP形式のカラー画像を読み込む関数
 image:    カラー画像配列
 xsize:    画像横サイズ
 ysize:    画像縦サイズ
 filename: ファイル名(.bmp)
-----------------------------------------------------------------------*/

void readBMP(char *filename,unsigned char image[ysize][xsize][3])
{
	long i,j,k;
	FILE *fp;

if((fp = fopen(filename,"rb")) == NULL){
		printf("readBMP: Open erorr!\n");
		exit(1);
	}

/* ヘッダー情報 */
	fread(&bfType, sizeof(bfType), 1, fp);
	fread(&bfSize, sizeof(bfSize), 1, fp);
	fread(&bfReserved1, sizeof(bfReserved1), 1, fp);
	fread(&bfReserved2, sizeof(bfReserved2), 1, fp);
	fread(&bfOffBits, sizeof(bfOffBits), 1, fp);

fread(&biSize, sizeof(biSize), 1, fp);
	fread(&biWidth, sizeof(biWidth), 1, fp);
	fread(&biHeight, sizeof(biHeight), 1, fp);
	fread(&biPlanes, sizeof(biPlanes), 1, fp);
	fread(&biBitCount, sizeof(biBitCount), 1, fp);
	fread(&biCompression, sizeof(biCompression), 1, fp);
	fread(&biSizeImage, sizeof(biSizeImage), 1, fp);
	fread(&biXPelsPerMeter, sizeof(biXPelsPerMeter), 1, fp);
	fread(&biYPelsPerMeter, sizeof(biYPelsPerMeter), 1, fp);
	fread(&biClrUsed, sizeof(biClrUsed), 1, fp);
	fread(&biClrImportant, sizeof(biClrImportant), 1, fp);

/* ビットマップデータ */
	for (i=0; i<biHeight; i++)
	{
		for (j=0; j<biWidth; j++)
			{
				for (k=0; k<3; k++) 
					fread(&image[i][j][2-k], 1, 1, fp);
			}
	}	
	fclose(fp);

}

/* 
   unsigned char型配列に格納された
   24ビット-ビットマップデータをBMPファイルに出力する
  [注意] 画像の横幅は４の倍数であること！
*/
void writeBmp(char *filename,unsigned char image[ysize][xsize][3])
{
	FILE *fp;
	int i, j, k;
	
	/* ファイルオープン */
	if ((fp = fopen(filename, "wb"))==NULL) {
		printf("writeBmp: Open error!\n");
		exit(1);
	}
	printf("output file : %s\n", filename);

/* ヘッダー情報 */
	fwrite(&bfType, sizeof(bfType), 1, fp);
	fwrite(&bfSize, sizeof(bfSize), 1, fp);
	fwrite(&bfReserved1, sizeof(bfReserved1), 1, fp);
	fwrite(&bfReserved2, sizeof(bfReserved2), 1, fp);
	fwrite(&bfOffBits, sizeof(bfOffBits), 1, fp);

fwrite(&biSize, sizeof(biSize), 1, fp);
	fwrite(&biWidth, sizeof(biWidth), 1, fp);
	fwrite(&biHeight, sizeof(biHeight), 1, fp);
	fwrite(&biPlanes, sizeof(biPlanes), 1, fp);
	fwrite(&biBitCount, sizeof(biBitCount), 1, fp);
	fwrite(&biCompression, sizeof(biCompression), 1, fp);
	fwrite(&biSizeImage, sizeof(biSizeImage), 1, fp);
	fwrite(&biXPelsPerMeter, sizeof(biXPelsPerMeter), 1, fp);
	fwrite(&biYPelsPerMeter, sizeof(biYPelsPerMeter), 1, fp);
	fwrite(&biClrUsed, sizeof(biClrUsed), 1, fp);
	fwrite(&biClrImportant, sizeof(biClrImportant), 1, fp);

/* ビットマップデータ */
	for (i=0; i<biHeight; i++)
	{
		for (j=0; j<biWidth; j++) 
		{
			for (k=0; k<3; k++)
			{ 
				fwrite(&image[i][j][2-k], 1, 1, fp);
			}
		}
	}
	fclose(fp);
}

/*--------------------      1次元FFTの実行      ----------------------
     a_rl: データの次数部(入出力兼用)
     a_im: データの虚数部(入出力兼用)
     ex:   データの個数を2のべき乗で与える（データの個数＝2のex乗個）
     inv:  1:FFT, -1:I_FFT
----------------------------------------------------------------------*/
int fft1(double a_rl[], double a_im[], int ex, int inv)
{
	int i,length=1;
	double *sin_tb1;    //sin計算用テーブル
	double *cos_tb1;    //cos計算用テーブル
	double *buf;        //作業用バッファ

/*- データの個数の計算 -*/
	for(i=0;i<ex;i++){
		length*=2;
	}

sin_tb1=(double *)malloc(length*sizeof(double));
	cos_tb1=(double *)malloc(length*sizeof(double));
	buf=(double *)malloc(length*sizeof(double));

if((sin_tb1==NULL)||(cos_tb1==NULL)||(buf==NULL)){
		return -1;
	}

cstb(length,inv, sin_tb1, cos_tb1);	
	fft1core(a_rl, a_im, length, ex, sin_tb1, cos_tb1, buf);
	free(sin_tb1);
	free(cos_tb1);
	return 0;

}

/*---------------     1次元FFTの計算の核になる部分     ---------------
 sin_tb1: sin計算用テーブル
 cos_tb1: cos計算用テーブル
---------------------------------------------------------------------*/
void fft1core(double a_rl[], double a_im[], int length, 
              int ex, double sin_tb1[], double cos_tb1[], double buf[])
{
	int i,j,k,w,j1,j2;
	int numb,lenb,timb;
	double xr, xi, yr, yi, nrml;

if(OPT==1){
		for(i=1;i<length;i+=2){
			a_rl[i]=-a_rl[i];
			a_im[i]=-a_im[i];
		}
	}

numb=1;
	lenb=length;
	for(i=0;i<ex;i++){
		lenb/=2;
		timb=0;
		for(j=0;j<numb;j++){
			w=0;
			for(k=0;k<lenb;k++){
				j1=timb+k;
				j2=j1+lenb;
				xr=a_rl[j1];
				xi=a_im[j1];
				yr=a_rl[j2];
				yi=a_im[j2];
				a_rl[j1]=xr+yr;
				a_im[j1]=xi+yi;
				xr=xr-yr;
				xi=xi-yi;
				a_rl[j2]=xr*cos_tb1[w]-xi*sin_tb1[w];
				a_im[j2]=xr*sin_tb1[w]+xi*cos_tb1[w];
				w+=numb;
			}
			timb+=(2*lenb);
		}
		numb*=2;
	}

birv(a_rl, length, ex, buf);     //実数データの並べ替え
	birv(a_im, length, ex, buf);     //虚数データの並べ替え

if(OPT==1){
		for(i=1;i<length;i+=2){
			a_rl[i]=-a_rl[i];
			a_im[i]=-a_im[i];
		}
	}

nrml=(double)(1.0/sqrt((double)length));
	for(i=0;i<length;i++){
		a_rl[i]*=nrml;
		a_im[i]*=nrml;
	}
	
}

/*-------------      sin・cosテーブルの作成      -------------*/
void cstb(int length, int inv, double sin_tb1[], double cos_tb1[])
{
	int i;
	double xx, arg;

xx=(double)(((-PI)*2.0)/(double)length);

if(inv<0){
		xx=-xx;
	}
	
	for(i=0;i<length;i++){
		arg=(double)i*xx;
		sin_tb1[i]=(double)sin(arg);
		cos_tb1[i]=(double)cos(arg);
	}
}

/*-----------      データの並べ換え      -----------
  a:  データの配列
  b:  作業用バッファ
--------------------------------------------------*/
void birv(double a[], int length, int ex, double b[])
{
	int i,ii,k,bit;

for(i=0;i<length;i++){
		for(k=0, ii=i, bit=0; ;bit<<=1, ii>>=1){
			bit=(ii&1)|bit;
			if(++k==ex)break;
		}
		b[i]=a[bit];
	}

for(i=0;i<length;i++){
		a[i]=b[i];
	}
}

/*------------------------------------------------------------------------------
                           2次元のFFTを実行する関数
  a_rl:  データ実数部（入出力兼用）
  a_im:  データ虚数部（入出力兼用）
  inv:   1:FFT, -1:I_FFT
------------------------------------------------------------------------------*/
int fft2(double a_rl[ysize][xsize], double a_im[ysize][xsize], int inv)
{
	int i;
	double *b_rl;       //データ転置作業用配列（実数部）
	double *b_im;       //データ転置作業用配列（虚数部）
	double *hsin_tb1;   //水平用sin計算用テーブル
	double *hcos_tb1;   //水平用cos計算用テーブル
	double *vsin_tb1;   //垂直用sin計算用テーブル
	double *vcos_tb1;   //垂直用cos計算用テーブル
	double *buf_x;      //作業用バッファ（水平方向）
	double *buf_y;      //作業用バッファ（垂直方向）

/*----------      memory allocation      ----------*/
	b_rl = (double *)calloc((size_t)xsize*ysize, sizeof(double));
	b_im = (double *)calloc((size_t)xsize*ysize, sizeof(double));
	hsin_tb1 = (double *)calloc((size_t)xsize, sizeof(double));
	hcos_tb1 = (double *)calloc((size_t)xsize, sizeof(double));
	vsin_tb1 = (double *)calloc((size_t)ysize, sizeof(double));
	vcos_tb1 = (double *)calloc((size_t)ysize, sizeof(double));
	buf_x = (double *)malloc((size_t)xsize*sizeof(double));
	buf_y = (double *)malloc((size_t)ysize*sizeof(double));

if((b_rl==NULL)||(b_im==NULL)||(hsin_tb1==NULL)||(hcos_tb1==NULL)||
		(vsin_tb1==NULL)||(vcos_tb1==NULL)||(buf_x==NULL)||(buf_y==NULL)){
		return -1;
	}

/*--- 水平用sin，cosテーブル作成 ---*/
	cstb(xsize, inv, hsin_tb1, hcos_tb1);

/*--- 垂直用sin，cosテーブル作成 ---*/
	cstb(ysize, inv, vsin_tb1, vcos_tb1);

/*--- 水平方向のFFT ---*/
	for(i=0;i<ysize;i++){
		fft1core(&a_rl[(long)i][0], &a_im[(long)i][0], 
                 xsize, X_EXP, hsin_tb1, hcos_tb1, buf_x);
	}

/*--- 2次元データの転置 ---*/
	rvmtx1((double (*)[xsize])a_rl, (double (*)[xsize])b_rl);
	rvmtx1((double (*)[xsize])a_im, (double (*)[xsize])b_im);

/*--- 垂直方向のFFT ---*/
	for(i=0;i<xsize;i++){
		fft1core(&b_rl[(long)ysize*i], &b_im[(long)ysize*i], 
                 ysize, Y_EXP, vsin_tb1, vcos_tb1, buf_y);
	}

/*--- 2次元データの転置 ---*/
	rvmtx2((double (*)[ysize])b_rl, (double (*)[ysize])a_rl);
	rvmtx2((double (*)[ysize])b_im, (double (*)[ysize])a_im);

free(b_rl);
	free(b_im);
	free(hsin_tb1);
	free(hcos_tb1);
	free(vsin_tb1);
	free(vcos_tb1);

return 0;

}

/*-----------------------------------------------------------------------------------
                            2次元データの転置を行う関数
  a:     2次元入力データ
  b:     2次元出力データ
  xsize: 水平データ個数
  ysize: 垂直データ個数
-----------------------------------------------------------------------------------*/
void rvmtx1(double a[ysize][xsize], double b[xsize][ysize])
{
	int i,j;

for(i=0;i<ysize;i++){
		for(j=0;j<xsize;j++){
			b[j][i] = a[i][j];
		}
	}
}

/*-----------------------------------------------------------------------------------
                            2次元データの転置を行う関数
  a:     2次元入力データ
  b:     2次元出力データ
  xsize: 水平データ個数
  ysize: 垂直データ個数
-----------------------------------------------------------------------------------*/
void rvmtx2(double a[xsize][ysize], double b[ysize][xsize])
{
	int i,j;

for(i=0;i<ysize;i++){
		for(j=0;j<xsize;j++){
			b[i][j] = a[j][i];
		}
	}
}

/*----------------------------------------------------------------------------------------
                              2次元FFTの結果を画像化するプログラム
			解説：このプログラムでは振幅特性を濃度として画像化を行っています．

/*-----------------------------------------------------------------------------------------
                            2次元FFTの結果を画像化する関数
  image_in:   入力画像配列（実世界領域）
  image_out:  出力画像配列（周波数領域）
------------------------------------------------------------------------------------------*/
int fftimge(unsigned char image_in[ysize][xsize][3], unsigned char image_out[ysize][xsize][3])
{
	long i,j;
	double *ar;   //データ実数部（入出力兼用）
	double *ai;   //データ虚数部（入出力兼用）
	double norm;  //振幅幅
	double max;
	double data;
	int k;

/*--------------     memory allocation     --------------*/
	ar = (double *)malloc((size_t)ysize*xsize*sizeof(double));
	ai = (double *)malloc((size_t)ysize*xsize*sizeof(double));

if((ar==NULL)||(ai==NULL)){
		return -1;
	}

/*--- 原画像を読み込み，2次元FFTの入力データに変換する ---*/
	for(i=0;i<ysize;i++){
		for(j=0;j<xsize;j++){
			
				ar[xsize*i+j] = (double)image_in[i][j][3];
				ai[xsize*i+j] = 0.0;
			}
		
	}

printf("aa\n");

/*--- 2次元FFTを実行する ---*/
	if(fft2((double (*)[xsize])ar, (double (*)[ysize])ai, 1)== -1){
		return -1;
	}

/*--- FFTの結果を画像化する ---*/
	max=0;
	for(i=0;i<ysize;i++){
		for(j=0;j<xsize;j++){
			norm = ar[xsize*i+j]*ar[xsize*i+j]
                  +ai[xsize*i+j]*ai[xsize*i+j];

if(norm!=0.0){
				norm = log(norm)/2.0;
			}else{
				norm =0.0;
			}

ar[xsize*i+j] = (double)norm;

if(norm>max){
				max=norm;
			}

}
	}

for(i=0;i<ysize;i++){
		for(j=0;j<xsize;j++){
			ar[xsize*i+j] = (double)(ar[xsize*i+j]*255/max);
		}
	}

/*--- FFTの結果を画像データに変換する ---*/
	for(i=0;i<ysize;i++)
	{
		for(j=0;j<xsize;j++)
		{
			for(k=0;k<3;k++)
			{
				data = ar[xsize*i+j];
				if(data > 255) data = 255;
				if(data <   0) data = 0;
				image_out[i][j][k] = (unsigned char)data;
			}
		}
	}

free(ar);
	free(ai);

return 0;

}

←解決時は質問者本人がここをチェックしてください。

戻る