# coding: utf-8 # --+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8 # NAME: dft.py # AUTHOR: CASPAR003 # DATE: 2013-07-27 # DESCRIPTION: 2次元DFT # --+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8 import cmath import math from os import path from PIL import Image from sys import argv, exit # --+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8 # main # --+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8 def main(infile): name, ext = path.splitext(infile) img = Image.open(infile) img = img.convert("L") pix = img.load() # pixの情報をもとに2次元配列を作る cim = [] for x in range(img.size[0]): buf = [] for y in range(img.size[1]): buf.append(complex(float(pix[x, y]) / 255.0, 0.0)) cim.append(buf) del buf # x方向のフーリエ変換 for y in range(img.size[1]): lin = [] for x in range(img.size[0]): lin.append(cim[x][y]) lin = dft(lin) for x in range(img.size[0]): cim[x][y] = lin[x] print "\rdft_x {0:5d}/{1:d} ({2:5.1f}%)".format(y+1, img.size[1], 50.0 * float(y+1) / float(img.size[1])), # y方向のフーリエ変換 for x in range(img.size[0]): lin = [] for y in range(img.size[1]): lin.append(cim[x][y]) lin = dft(lin) for y in range(img.size[1]): cim[x][y] = lin[y] print "\rdft_y {0:5d}/{1:d} ({2:5.1f}%)".format(x+1, img.size[1], 50.0 * float(x+1) / float(img.size[1]) + 50.0), print "" # スケーリング用にフーリエ変換後の最大値を見つけておく cim_amax = float(0.0) for y in range(img.size[1]): for x in range(img.size[0]): if (cim_amax < abs(cim[x][y])): cim_amax = abs(cim[x][y]) for y in range(int(math.ceil(img.size[1] / 2.0))): for x in range(int(math.ceil(img.size[0] / 2.0))): # cim第二象限 → pix第四象限 pix[int(math.floor(img.size[0] / 2.0)) + x, int(math.floor(img.size[1] / 2.0)) + y] = int(100.0 * 255.0 * (abs(cim[x][y]) / cim_amax)**2) # cim第一象限 → pix第三象限 pix[x, int(math.floor(img.size[1] / 2.0)) + y] = int(100.0 * 255.0 * (abs(cim[int(math.floor(img.size[0] / 2.0)) + x][y]) / cim_amax)**2) # cim第三象限 → pix第一象限 pix[int(math.floor(img.size[0] / 2.0)) + x, y] = int(100.0 * 255.0 * (abs(cim[x][int(math.floor(img.size[1] / 2.0)) + y]) / cim_amax)**2) # cim第四象限 → pix第二象限 pix[x, y] = int(100.0 * 255.0 * (abs(cim[int(math.floor(img.size[0] / 2.0)) + x][int(math.floor(img.size[1] / 2.0)) + y]) / cim_amax)**2) img.save(name + "_dft" + ext) # --+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8 # dft # --+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8 def dft(f): N = len(f) s = [] for j in range(N): sum = complex(0.0, 0.0) for i in range(N): tht = 2.0 * cmath.pi / float(N) * float(j) * float(i) sum += f[i] * cmath.exp(complex(0.0, -tht)) s.append(sum) return s # --+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8 # 実行ここから # --+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8 if (len(argv) < 2): print("usage: dft.py ") exit() main(argv[1])