チーム21B4
課題名
c言語による画像のFFTと画像処理
研究者名
3-14-035 Sugita Motoki
3-14-041 Nishimura Masahiro
概要
チーム21C2では、画像に対してフーリエ変換を行った後の結果を配列に入れ、スペクトル画像を表示するとともに逆フーリエ変換で元の画像に戻していた。そこで、フーリエ変換後の配列の内容を直接編集することによって逆フーリエ変換後の画像を編集できるのではないかと考え実験を行った。
ソースコード
ソースコードはチーム21C2のものを引用して少し変更を加えた。
#include <stdio.h>
#include <stdlib.h>
#include <gd.h>
#include <math.h>
#include "fft.h"
int main(const int argc,const char *argv[]){
FILE *out,*out2,*in; gdImagePtr im,im_new; int width,height,i,j,color,r,g,b,pixel,x,y; double result_y_re[N_points][N_points],result_y_imm[N_points][N_points]; double result_re[N_points][N_points],result_imm[N_points][N_points]; double tmp_re[N_points][N_points],tmp_imm[N_points][N_points]; double result[N_points][N_points]; double max=0; double p,re,imm;
//第1引数で指定されたファイルを読み出し用にオープン if((in=fopen(argv[1],"r"))==NULL){ printf("file open error for %s\n",argv[1]); exit(-1); } //第2引数で指定されたファイルを書き出し用にオープン if((out=fopen(argv[2],"wb"))==NULL){ printf("file open error for %s\n",argv[2]); exit(-1); } //第3引数で指定されたファイルを書き出し用にオープン if((out2=fopen(argv[3],"wb"))==NULL){ printf("file open error for %s\n",argv[3]); exit(-1); }
//im に画像を読み込み im = gdImageCreateFromJpeg(in);
//入力画像のサイズを取得 width=gdImageSX(im); height=gdImageSY(im);
for(i=0;i<width;i++){ // y軸でFFT for(j=0;j<height;j++){ int a; //im の (i,j) におけるカラーインデックスの取得 pixel=gdImageGetPixel(im,i,j); //im の (i,j) における r,g,b の値を取得 r=gdImageRed(im,pixel); g=gdImageGreen(im,pixel); b=gdImageBlue(im,pixel); // 白黒変換 a=(r+g+b)/3; data_re[j]=a; } FFT(1,exponent,data_re,data_imm); for(j=0;j<height;j++){ result_y_re[i][j]=data_re[j]; result_y_imm[i][j]=data_imm[j]; } }
for(i=0;i<width;i++){ // x軸でFFT for(j=0;j<height;j++){ data_re[j]=result_y_re[j][i]; data_imm[j]=result_y_imm[j][i]; } FFT(1,exponent,data_re,data_imm); for(j=0;j<height;j++){ result[j][i]=mod[j]; result_re[j][i]=data_re[j]; result_imm[j][i]=data_imm[j]; if(mod[j]>max){ max=mod[j]; } } }
//新しい画像を用意 im_new= gdImageCreateTrueColor(width,height); // 周波数を描画 for(i=0;i<width;i++){ for(j=0;j<height;j++){ double a; a=result[i][j]/max*255*1000; if(a>255)a=255; color=gdImageColorExact(im_new,a,a,a); gdImageSetPixel(im_new,i,j,color); } } gdImageJpeg(im_new,out,-1);
for(i=0;i<width;i++){ for(j=0;j<height;j++){ tmp_re[i][j]=result_re[i][j]; tmp_imm[i][j]=result_imm[i][j]; } }
for(i=0;i<width;i++){ for(j=0;j<height;j++){
result_re[i][j]=tmp_re[j][i]; result_imm[i][j]=tmp_imm[j][i]; } }
//逆フーリエ変換を行いその結果を描画
p=2*M_PI/width; for(y=0;y<height;y++){ for(x=0;x<width;x++){ double re=0,imm=0,a=0; for(i=0;i<width;i++){ for(j=0;j<height;j++){ re+=(result_re[i][j]*cos(p*(x*i+y*j))-result_imm[i][j]*sin(p*(x*i+y*j))); imm+=(result_re[i][j]*sin(p*(x*i+y*j))+result_imm[i][j]*cos(p*(x*i+y*j))); } } a=sqrt(re*re+imm*imm); if(a>255)a=255; color=gdImageColorExact(im_new,a,a,a); gdImageSetPixel(im_new,x,y,color); } } gdImageJpeg(im_new,out2,-1); fclose(in); fclose(out); fclose(out2);
return 0;
}
#include "math.h" //mathematical library
#define N_points 64 //number of points
#define exponent log(64)/log(2) //log2(N_points); for N_points=64 -> exponent=6
double mod[N_points]={0}; //arrays
double data_re[N_points]={0};
double data_imm[N_points]={0};
short FFT(int dir,long m,double *x,double *y)
{
long n,i,i1,j,k,i2,l,l1,l2;
double c1,c2,tx,ty,t1,t2,u1,u2,z;
/* Calculate the number of points */
n = 1;
for (i=0;i<m;i++)
n *= 2;
/* Do the bit reversal */
i2 = n >> 1;
j = 0;
for (i=0;i<n-1;i++) {
if (i < j) {
tx = x[i];
ty = y[i];
x[i] = x[j];
y[i] = y[j];
x[j] = tx;
y[j] = ty;
}
k = i2;
while (k <= j) {
j -= k;
k >>= 1;
}
j += k;
}
/* Compute the FFT */
c1 = -1.0;
c2 = 0.0;
l2 = 1;
for (l=0;l<m;l++) {
l1 = l2;
l2 <<= 1;
u1 = 1.0;
u2 = 0.0;
for (j=0;j<l1;j++) {
for (i=j;i<n;i+=l2) {
i1 = i + l1;
t1 = u1 * x[i1] - u2 * y[i1];
t2 = u1 * y[i1] + u2 * x[i1];
x[i1] = x[i] - t1;
y[i1] = y[i] - t2;
x[i] += t1;
y[i] += t2;
}
z = u1 * c1 - u2 * c2;
u2 = u1 * c2 + u2 * c1;
u1 = z;
}
c2 = sqrt((1.0 - c1) / 2.0);
if (dir == 1)
c2 = -c2;
c1 = sqrt((1.0 + c1) / 2.0);
}
/* Scaling factor for forward transform */
if (dir == 1) {
for (i=0;i<n;i++) {
x[i] /= n;
y[i] /= n;
}
}
for(i=0;i<N_points;i++)
/* Absolute value */
mod[i]=sqrt((data_re[i]*data_re[i])+(data_imm[i]*data_imm[i]))/N_points;
}
実験結果
具体的な実験内容
フーリエ変換後の結果を保存した配列result_re、result_immの内容を、一度tmp_re、tmp_immに移した後、
様々な処理を行ってからもう一度result_re、result_immに戻す。
各画像に対する処理などは以下のとおりである。
1枚目
縦と横を入れ替えて保存する。
for (i = 0; i < width; i++) {
for (j = 0; j < height; j++) { result_re[i][j] = tmp_re[i][j]; result_imm[i][j] = tmp_imm[i][j]; }
}
2枚目
縦横それぞれ半分の長さに対し、resultの値を一律で0にする
for (i = 0; i < width / 2; i++) {
for (j = 0; j < height / 2; j++) { result_re[i][j] = 0; result_imm[i][j] = 0; }
}
3枚目
縦横それぞれ半分の長さに対し、実数部分のみ配列内を100ずらす
for (i = 0; i < width / 2; i++) {
for (j = 0; j < height / 2; j++) { result_re[i][j] = tmp_re[(i+100)%width][(j+100)%height]; result_imm[i][j] = tmp_imm[i][j] }
}
4枚目
縦横それぞれ半分の長さに対し、虚数部分のみ配列内を100ずらす
for (i = 0; i < width / 2; i++) {
for (j = 0; j < height / 2; j++) { result_re[i][j] = tmp_re[i][j]; result_imm[i][j] = tmp_imm[(i+100)%width][(j+100)%height]; }
}
5枚目
縦横それぞれ半分の長さに対し、配列の内容に10を加える
for (i = 0; i < width / 2; i++) {
for (j = 0; j < height / 2; j++) { result_re[i][j] = tmp_re[i][j]+10; result_imm[i][j] = tmp_imm[i][j]+10; }
}
6枚目
縦横それぞれ半分の長さに対し、配列の内容に50を加える
for (i = 0; i < width / 2; i++) {
for (j = 0; j < height / 2; j++) { result_re[i][j] = tmp_re[i][j]+50; result_imm[i][j] = tmp_imm[i][j]+50; }
}
7枚目
縦横それぞれ半分の長さに対し、配列の内容に100を加える
for (i = 0; i < width / 2; i++) {
for (j = 0; j < height / 2; j++) { result_re[i][j] = tmp_re[i][j]+100; result_imm[i][j] = tmp_imm[i][j]+100; }
}
5,6,7枚目から数を増やせば増やすほどぼんやりしていることが分かる。これは普通の場合は0を足していると考えることが出来る。
考察
実験開始前はフーリエ変換した後に処理を行った後で逆フーリエ変換することにより画像処理が行えると考えたが、実際は画像の通りとなった。
結果からフーリエ変換することにより周波数に変換しているので、画像と周波数ではものが違うことが要因と考える。
以上の点から周波数に対する画像処理を行うためにはローパスフィルタ、ハイパスフィルタ、バンドパスフィルタのように周波数の信号に関わるフィルターが必要であると考える。
- 最終更新:2021-07-12 17:07:29