チーム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;
}


実験結果

21B4_1.jpg21B4_2.jpg
21B4_3.jpg21B4_4.jpg
21B4_7.jpg21B4_6.jpg21B4_5.jpg

具体的な実験内容

フーリエ変換後の結果を保存した配列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

このWIKIを編集するにはパスワード入力が必要です

認証パスワード