２次元フーリエをするには？
入力画像サイズが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;
}
うまく結果が出ない，というのはどういう状況なのですか？
それがちゃんと説明されないと，まともに答えをつけられる人はいないと思います。
＃プログラムの結果は，あなたの書いたコードに対する正しい振る舞いだと思います。
また，デバッガで追ってみて，異常な振る舞いをする場所は無かったのですか？
ちなみに，あまりに面倒なので全くコードは読んでいません。
コードを読んで欲しいなら抜粋して必要部分だけにして，さらに必要なコメントを付けることです。
こんな長ったらしいソースつけられても、誰も最後まで見ないと思うが...
中身はよーわからんが、とりあえず以下を直してみれば。
> #define X_EXP 7 //X_EXP = log2(X_SIZE)
> #define Y_EXP 7 //Y_EXP = log2(Y_SIZE)
> #define ysize 128
> #define xsize 128
7 -> 8
128 ->256
> #define PI 3.141592
もうちょい、精度上げたほうがいいのでは。
いろいろと御意見ありがとうございます！
まだまだ自分が未熟なのでもっと勉強します。
YuOさん、ボコノン教徒さん、ありがとうございます！
| ツイート |  |