画像を2DFT→2IDFTして表示

フォーラム(掲示板)ルール
フォーラム(掲示板)ルールはこちら　　※コードを貼り付ける場合は [code] と [/code] で囲って下さい。詳しくはこちら
閉鎖
おちゃのみず

画像を2DFT→2IDFTして表示

#1

投稿記事 by おちゃのみず » 9年前

お世話になります。visualstudio2010でプログラム作っています。
画像→2DFT→2IDFTし、基の画像を作成するプログラムを自作しましたが、できあがったビットマップ画像は真っ黒でなにも見えません。
128×128のグレースケール256階調のビットマップ(よくあるレナの画像)から画像データの配列を取ってきて、DFT→IDFTしています。
得られた結果を256階調で表示できるように正規化しているつもりです。
ビットマップ画像の読み込みと書き出しの関数（bmp.h、bmp.cというファイル中で処理）は動作確認済みですので、DFTとIDFTに問題があると思います。

おかしなところがあれば、ご指摘いただけると大変助かります。
よろしくお願いします。

コード: 

#include <stdio.h>
#include <math.h>
#include<stdlib.h>
#include "bmp.h"


#define PI (3.1415926535897)

#define COLORTONE (256)		//画像の階調

// 画像データの高さ&幅の構造体
typedef struct tagIMAGESAIZE{
	unsigned int imageHeight;
	unsigned int imageWidth;
} IMAGESIZE;

#define COLORPALETNUM (256) // カラーパレット構造体の数
#define MAX_IMAGE_HEIGHT (2048)	// このプログラムで扱うイメージデータの最大高さ
#define MAX_IMAGE_WIDTH (4096)	// このプログラムで扱うイメージデータの最大幅

// 最大値見つける関数
double maxElem(double *ary, unsigned int arySize){
	double max = ary[0];
	for(unsigned int i = 1; i < arySize; i++){
		if(ary[i] > max){
			max = ary[i];
		}
	}
	return max;
}

// 最小値見つける関数
double minElem(double *ary, unsigned int arySize){
	double min = ary[0];
	for(unsigned int i = 1; i < arySize; i++){
		if(ary[i] < min){
			min = ary[i];
		}
	}
	return min;
}


// DFT関数
void DFT(IMAGESIZE imageSize, unsigned char* imageDataBuf, unsigned char **DFTImageDataBufP);

int main(void)
{	
	// 画像データバッファ
	unsigned char *readImageDataBuf = new unsigned char[MAX_IMAGE_HEIGHT * MAX_IMAGE_WIDTH];
	unsigned char *DFTImageDataBuf = new unsigned char[MAX_IMAGE_HEIGHT * MAX_IMAGE_WIDTH];

	int errorCheck; // 関数のエラーチェック用

	// ----  ビットマップファイルの読み込み・画像データ取得
	
	// 画像データのサイズ(高さと幅)
	IMAGESIZE ImageSize;
	ImageSize.imageHeight = 0;
	ImageSize.imageWidth = 0;


	// ファイルオープン
	FILE *rFp;
	if(fopen_s(&rFp, "Lenna128.bmp", "rb") != 0){
		printf("File open failure.\n");
		exit(1);
	}

	// 画像データの読み込み
	errorCheck = read_BMP(rFp, &ImageSize, readImageDataBuf);
	if(errorCheck){
		printf("Reading failure.");
		exit(1);
	}

	// 2次元フーリエ変換のデータ取得のためにポインタを設定
	unsigned char **DFTImageDataBufP = new unsigned char*[ImageSize.imageHeight * ImageSize.imageWidth];
	for(unsigned int i = 0; i < (ImageSize.imageHeight * ImageSize.imageWidth); i++){
		DFTImageDataBufP[i] = &DFTImageDataBuf[i];
	}

	// ---- 2次元フーリエ変換

	DFT(ImageSize, readImageDataBuf, DFTImageDataBufP);


	// ---- ビットマップファイルの作成・画像データの書き込み
	// 書き込み用pngファイルの作成
	FILE *wFp;
	if(fopen_s(&wFp, "Lenna128_After.bmp","wb")){
		printf("Writing failure.");
		exit(1);
	}

	// 画像データの書き込み
	errorCheck = write_BMP(wFp, &ImageSize, DFTImageDataBuf);
	if(errorCheck){
		printf("writing failure.");
		exit(1);
	}
	
	fclose(rFp);
	fclose(wFp);

	delete readImageDataBuf;
	delete DFTImageDataBuf;
	return 0;
}


void DFT(IMAGESIZE imageSize, unsigned char* imageDataBuf, unsigned char **DFTImageDataBufP){
	
	// フーリエ変換用データバッファ
	double *ReBuf = new double[imageSize.imageHeight * imageSize.imageWidth];
	double *ImBuf = new double[imageSize.imageHeight * imageSize.imageWidth];
	double *fspectl = new double[imageSize.imageHeight * imageSize.imageWidth];
	double *DFTBUF = new double[imageSize.imageHeight * imageSize.imageWidth];

	unsigned char** DFTDataPtr = new unsigned char* [imageSize.imageHeight * imageSize.imageWidth];
	for(unsigned int i = 0; i < (imageSize.imageHeight * imageSize.imageWidth); i++){
		DFTDataPtr[i] = DFTImageDataBufP[i];
	}
	
	// フーリエ変換
	for( unsigned int v = 0; v < imageSize.imageHeight; v++){		// y方向の空間周波数についてのループ
		for( unsigned int u = 0; u < imageSize.imageWidth; u++){	// x方向の空間周波数についてのループ
			ReBuf[v * imageSize.imageWidth + u] = 0.0;				// 周波数領域要素の初期化
			ImBuf[v * imageSize.imageWidth + u] = 0.0;				// 周波数領域要素の初期化

			for( unsigned int y = 0; y < imageSize.imageHeight; y++){	
				for( unsigned int x = 0; x < imageSize.imageWidth; x++){
					
					ReBuf[v * imageSize.imageWidth + u] += (double)imageDataBuf[y * imageSize.imageWidth + x] * cos(2.0 * PI * (( u / (double)imageSize.imageWidth) + (v / (double)imageSize.imageHeight)));
					ImBuf[v * imageSize.imageWidth + u] += (double)imageDataBuf[y * imageSize.imageWidth + x] * sin(2.0 * PI * (( u / (double)imageSize.imageWidth) + (v / (double)imageSize.imageHeight)));
				}
			}
			ReBuf[v * imageSize.imageWidth + u] /= ((double)imageSize.imageWidth * (double)imageSize.imageHeight);
			ImBuf[v * imageSize.imageWidth + u] /= ((double)imageSize.imageWidth * (double)imageSize.imageHeight);
		}
	}

		
	// 逆フーリエ変換
	double *InvReBuf = new double[imageSize.imageHeight * imageSize.imageWidth];
	double *InvImBuf = new double[imageSize.imageHeight * imageSize.imageWidth];
	double *InvDFTBUF = new double[imageSize.imageHeight * imageSize.imageWidth];


	for( unsigned int y = 0; y < imageSize.imageHeight; y++){		// y方向の空間周波数についてのループ
		for( unsigned int x = 0; x < imageSize.imageWidth; x++){	// x方向の空間周波数についてのループ
			InvReBuf[y * imageSize.imageWidth + x] = 0.0;				// 周波数領域要素の初期化
			InvImBuf[y * imageSize.imageWidth + x] = 0.0;				// 周波数領域要素の初期化

			for( unsigned int v = 0; v < imageSize.imageHeight; v++){	
				for( unsigned int u = 0; u < imageSize.imageWidth; u++){
					
					InvReBuf[y * imageSize.imageWidth + x] += ReBuf[v * imageSize.imageWidth + u] * cos(2.0 * PI * (( u / (double)imageSize.imageWidth) + (v / (double)imageSize.imageHeight)))
															- ImBuf[v * imageSize.imageWidth + u] * sin(2.0 * PI * (( u / (double)imageSize.imageWidth) + (v / (double)imageSize.imageHeight)));

					InvImBuf[y * imageSize.imageWidth + x] += ReBuf[v * imageSize.imageWidth + u] * sin(2.0 * PI * (( u / (double)imageSize.imageWidth) + (v / (double)imageSize.imageHeight)))
															+ ImBuf[v * imageSize.imageWidth + u] * cos(2.0 * PI * (( u / (double)imageSize.imageWidth) + (v / (double)imageSize.imageHeight)));

					
					InvDFTBUF[y * imageSize.imageWidth + x] = sqrt((InvReBuf[y * imageSize.imageWidth + x] * InvReBuf[y * imageSize.imageWidth + x])
																+(InvImBuf[y * imageSize.imageWidth + x] * InvImBuf[y * imageSize.imageWidth + x]));
				}
			}
		}
	}

	double *logInvDFTBUF = new double[imageSize.imageHeight * imageSize.imageWidth];
	for(unsigned int i = 0; i < (imageSize.imageHeight* imageSize.imageWidth); i++){		// 0～255の間に値を正規化
		logInvDFTBUF[i] = log10(InvDFTBUF[i]);
	}
	
	
	// 256のグレースケールで表現するための正規化作業
	double maxInvDFTBUFElem = maxElem(&logInvDFTBUF[0], (imageSize.imageHeight* imageSize.imageWidth));
	double minInvDFTBUFElem = minElem(&logInvDFTBUF[0], (imageSize.imageHeight* imageSize.imageWidth));

	double *standInvDFTBUF = new double[imageSize.imageHeight * imageSize.imageWidth];

	for(unsigned int i = 0; i < (imageSize.imageHeight* imageSize.imageWidth); i++){		// 0～255の間に値を正規化
		standInvDFTBUF[i] = (fabs(logInvDFTBUF[i] - minInvDFTBUFElem) * COLORTONE) / ( maxInvDFTBUFElem - minInvDFTBUFElem);
		*(DFTDataPtr[i]) = (unsigned char)standInvDFTBUF[i];
	}
		
	delete(ReBuf);
	delete(ImBuf);
	delete(DFTBUF);
	delete (standInvDFTBUF);

	delete(DFTDataPtr);


}
ページトップ
アバター
h2so5
副管理人
記事: 2212
登録日時: 13年前
住所: 東京
連絡を取る:
h2so5 に連絡する

Re: 画像を2DFT→2IDFTして表示

#2

投稿記事 by h2so5 » 9年前

ぱっと見て気になったところは

・165行目でInvDFTBUFの同じ座標に何度も上書きしている（ループの最後の値しか使っていない）
・配列の解放にdeleteを使っている（delete[] を使うのが正しい）
・不要な配列が多い（一時的にしか使わない値をわざわざ配列に保存する意味がない、 182行目など）
・引数として渡した配列を書き換えるのにダブルポインタを使う必要がない（配列は実質的に参照渡しなので）
ページトップ
sleep

Re: 画像を2DFT→2IDFTして表示

#3

投稿記事 by sleep » 9年前

私もさらっとしか見てませんが
wikipediaの離散フーリエ変換(2次元での変換)と見比べた感じだと
計算式に xとy および DEFの sin cos に －（マイナス）記号が足りない様に見えます。
（ただし、－（マイナス）記号に関しては DEFかIDEFどちらか片方にさえ付いていれば良いのかな？そのへんあまり詳しくないですが）
とりあえず、wikipedia通り 変換側（おちゃのみずさんのコードの場合逆変換側）に －（マイナス）記号を付けてます。
おちゃのみず さんが書きました:

コード: 

ReBuf[v * imageSize.imageWidth + u] += (double)imageDataBuf[y * imageSize.imageWidth + x] * cos(2.0 * PI * (( u / (double)imageSize.imageWidth) + (v / (double)imageSize.imageHeight)));
ImBuf[v * imageSize.imageWidth + u] += (double)imageDataBuf[y * imageSize.imageWidth + x] * sin(2.0 * PI * (( u / (double)imageSize.imageWidth) + (v / (double)imageSize.imageHeight)));

コード: 

ReBuf[v * imageSize.imageWidth + u] += (double)imageDataBuf[y * imageSize.imageWidth + x] * cos( 2.0 * PI * (( u * x / (double)imageSize.imageWidth) + ( v * y / (double)imageSize.imageHeight)));
ImBuf[v * imageSize.imageWidth + u] += (double)imageDataBuf[y * imageSize.imageWidth + x] * sin( 2.0 * PI * (( u * x / (double)imageSize.imageWidth) + ( v * y / (double)imageSize.imageHeight)));
おちゃのみず さんが書きました:

コード: 

InvReBuf[y * imageSize.imageWidth + x] += ReBuf[v * imageSize.imageWidth + u] * cos(2.0 * PI * (( u / (double)imageSize.imageWidth) + (v / (double)imageSize.imageHeight)))
										- ImBuf[v * imageSize.imageWidth + u] * sin(2.0 * PI * (( u / (double)imageSize.imageWidth) + (v / (double)imageSize.imageHeight)));

コード: 

InvReBuf[y * imageSize.imageWidth + x] += ReBuf[v * imageSize.imageWidth + u] * cos( -2.0 * PI * (( u * x / (double)imageSize.imageWidth) + ( v * y / (double)imageSize.imageHeight)))
                                        - ImBuf[v * imageSize.imageWidth + u] * sin( -2.0 * PI * (( u * x / (double)imageSize.imageWidth) + ( v * y / (double)imageSize.imageHeight)));
ページトップ
おちゃのみず

Re: 画像を2DFT→2IDFTして表示

#4

投稿記事 by おちゃのみず » 9年前

ご指摘たくさんありがとうございます！おちゃのみずです。
いろいろおかしな点が気づけて大変助かりました。
確かに、ダブルポインタは使わなくてもいいですね＾＾；


もう少し悩むところがあってまだ解決しておりませんが、
取り急ぎ、お礼申し上げます。

また何かしら質問させていただいた際は、よろしくお願いします。
ページトップ
閉鎖

“C言語何でも質問掲示板” へ戻る