画像のフーリエ変換DFT

![result.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510257/eeb89873-89bf-4acb-8404-5a7116845cbc.png)

# はじめに
Processingで画像のフーリエ変換、逆変換をするプログラムです。
画像を、三角関数の集合に変換し、それをまた逆に変換すると、もとに戻るというアレです。

OpenCVは使っていません。
複素数ライブラリ(Complex Numbers)も使っていません。


# トップ画像の説明
画像の左から、
元画像、フーリエ変換画像、その絶対値の対数表示、逆変換により戻した画像
です。

## スペクトラム表示（絶対値の対数表示）について

よく見るフーリエ変換画像は、中央が「輝く星」のような画像になります。
しかし、この画像は四隅が明るいです。
よく見るフーリエ変換画像は、上下方向に半分、左右方向に半分シフト```cv2.fftshift()```しているためです。

よって、トップ画像の場合、
中心が高周波で、周辺が低周波です。

一般的なフーリエ変換画像は
中心が低周波で、周辺が高周波です。

# 元画像
AIに描かせた、16x16の4階調の画像です。
ファミコンのマリオに習って、4階調にしています。深い意図はありません。

![girl.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510257/e42445a5-3761-4ec3-8d46-b2b83f982a19.png)

# ぼかし＝高周波成分を削除、低周波成分を残す
Jpeg圧縮は、高周波成分を削除して容量削減しているので、同じ原理と考えてもいいでしょう。
ドット絵のような画像は、エッジの部分が高周波成分であるため、画質の劣化が顕著に現れることも同様です。

![result.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510257/9c64b6ba-adc0-41b3-a355-325a2d0e5460.png)

# プログラム

```java:Processing
int N = 16;
PImage img;
float[][] image = new float[N][N];
Complex[][] F = new Complex[N][N];
Complex[][] F_H = new Complex[N][N];
Complex[][] result = new Complex[N][N];

void setup() {
  size(1024, 256);
  noLoop();

  // 1. 画像読み込み（16x16想定）
  img = loadImage("girl.png");
  for (int y = 0; y < N; y++) {
    for (int x = 0; x < N; x++) {
      float r = red(img.get(x, y)) / 255.0;
      image[y][x] = r;
    }
  }

  // 2. DFT行列 F, F_H の作成
  for (int m = 0; m < N; m++) {
    for (int n = 0; n < N; n++) {
      float angle = -TWO_PI * m * n / N;
      F[m][n] = new Complex(cos(angle), sin(angle));
      F_H[n][m] = new Complex(cos(angle), -sin(angle)).div(N); // 1D逆変換なら1/N
    }
  }

  // 3. 行方向 DFT: rowTransformed = image × F^T
  Complex[][] rowTransformed = new Complex[N][N];
  for (int y = 0; y < N; y++) {
    for (int x = 0; x < N; x++) {
      rowTransformed[y][x] = new Complex(0, 0);
      for (int k = 0; k < N; k++) {
        rowTransformed[y][x] = rowTransformed[y][x].add(F[k][x].mul(image[y][k]));
      }
    }
  }

  // 4. 列方向 DFT: spectrum = F × rowTransformed
  Complex[][] spectrum = new Complex[N][N];
  for (int y = 0; y < N; y++) {
    for (int x = 0; x < N; x++) {
      spectrum[y][x] = new Complex(0, 0);
      for (int k = 0; k < N; k++) {
        spectrum[y][x] = spectrum[y][x].add(F[y][k].mul(rowTransformed[k][x]));
      }
    }
  }



  // 高周波除去：半径R以内を残して他は0
  int R = 6;  // ぼかし強度（大きいほど強くぼかす）

  for (int y = 0; y < N; y++) {
    for (int x = 0; x < N; x++) {
      int dy = y - N/2;
      int dx = x - N/2;
      float dist = sqrt(dx*dx + dy*dy);

      if (dist < R) {
        spectrum[y][x] = new Complex(0, 0);  // 消す（高周波カット）
      }
    }
  }

  // === 逆変換 ===

  // 5. 列方向 IDFT: colReconstructed = F_H × spectrum
  Complex[][] colReconstructed = new Complex[N][N];
  for (int y = 0; y < N; y++) {
    for (int x = 0; x < N; x++) {
      colReconstructed[y][x] = new Complex(0, 0);
      for (int k = 0; k < N; k++) {
        colReconstructed[y][x] = colReconstructed[y][x].add(F_H[y][k].mul(spectrum[k][x]));
      }
    }
  }

  // 6. 行方向 IDFT: result = colReconstructed × F_H^T
  for (int y = 0; y < N; y++) {
    for (int x = 0; x < N; x++) {
      result[y][x] = new Complex(0, 0);
      for (int k = 0; k < N; k++) {
        result[y][x] = result[y][x].add(F_H[k][x].mul(colReconstructed[y][k]));
      }
    }
  }

  // 7. 平均誤差を評価
  float totalError = 0;
  for (int y = 0; y < N; y++) {
    for (int x = 0; x < N; x++) {
      float diff = abs(image[y][x] - result[y][x].re);
      totalError += diff;
    }
  }
  println("平均誤差:", totalError / (N * N));


  // 8. 描画
  drawMatrix(image, 0, 0);
  drawMatrixComplex(spectrum, 256, 0);
  drawSpectrum(spectrum, 512, 0);
  drawMatrixComplex(result, 768, 0);
  
  //save("result.png");
}

void drawMatrix(float[][] mat, int offsetX, int offsetY) {
  float scale = 16;
  for (int y = 0; y < N; y++) {
    for (int x = 0; x < N; x++) {
      fill(mat[y][x] * 255);
      noStroke();
      rect(offsetX + x * scale, offsetY + y * scale, scale, scale);
    }
  }
}

void drawMatrixComplex(Complex[][] mat, int offsetX, int offsetY) {
  float scale = 16;

  for (int y = 0; y < N; y++) {
    for (int x = 0; x < N; x++) {
      float v = map(mat[y][x].re, 0, 1, 0, 255);
      fill(constrain(v, 0, 255));
      noStroke();
      rect(offsetX + x * scale, offsetY + y * scale, scale, scale);
    }
  }
}

// 複素数クラス
class Complex {
  float re, im;
  Complex(float r, float i) {
    re = r;
    im = i;
  }
  Complex add(Complex c) {
    return new Complex(re + c.re, im + c.im);
  }
  Complex mul(float f) {
    return new Complex(re * f, im * f);
  }
  Complex mul(Complex c) {
    return new Complex(re * c.re - im * c.im, re * c.im + im * c.re);
  }
  Complex div(float f) {
    return new Complex(re / f, im / f);
  }
}

void drawSpectrum(Complex[][] spectrum, int offsetX, int offsetY) {
  float scale = 16;

  float maxAmp = 0;
  float[][] amps = new float[N][N];

  // 振幅（絶対値）をlogスケーリング用に取得
  for (int y = 0; y < N; y++) {
    for (int x = 0; x < N; x++) {
      float amp = sqrt(spectrum[y][x].re * spectrum[y][x].re + spectrum[y][x].im * spectrum[y][x].im);
      amp = log(1 + amp);  // logスケーリング
      amps[y][x] = amp;
      if (amp > maxAmp) maxAmp = amp;
    }
  }

  // グレースケールで描画
  for (int y = 0; y < N; y++) {
    for (int x = 0; x < N; x++) {
      float norm = map(amps[y][x], 0, maxAmp, 0, 255);
      fill(norm);
      noStroke();
      rect(offsetX + x * scale, offsetY + y * scale, scale, scale);
    }
  }
}
```
# 参考サイト

このサイトは便利

https://monman53.github.io/2dfft/

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510257/f05e414d-6e9c-4279-bf8e-427e984976c9.png)

https://vigne-cla.com/25-1/