CPUでの2次元FFTの実装は前の項目の図（図13.2「図：2次元FFTの手順」）で解説した通りとなります。
- PNG画像をバイト型で読み込み
- RGBチャネルを一つごと(Red,Green,Blue)に独立した変数に記憶
- 読み込んだ画像の走査線を行スキャン
- 各行をスキャンした配列・スライスのFFT
- 変換した2次元スライスの走査線を列スキャン
- 各列をスキャンした配列・スライスのFFT
- フィルターの適用
- 各列をスキャンした配列・スライスの逆FFT
- 各行をスキャンした配列・スライスの逆FFT
この手順のFFTでは、FFTの前にビット反転が行われるものとします。
この手順を実装すると下記のようなソースコードとなります。
FFTCPU2D.java.
package com.book.jocl.fft;
import java.awt.Graphics;
import java.awt.image.BufferedImage;
import java.awt.image.DataBufferByte;
import java.io.File;
import java.net.URL;
import java.nio.file.Paths;
import javax.imageio.ImageIO;
import javax.swing.ImageIcon;
public class FFTCPU2D {
private static final String IMAGE_PATH = "SAMPLE.png";
private static ImageIcon img;
private static int width;
private static int height;
private static double _PI = Math.PI;
private static final int _RowSize = 256;
private static final int _ColSize = 256;
private static final int size = _RowSize*_ColSize;
private static int _SIGN = 1;
private static void forward(int N, int offset, double[] data, double[] F, int sign, int step) {
int N2 = N/2;
double c, s;
double real, imaginary;
step++;
if(N2 != 2) {
forward(N2,offset,data,F,sign,step);
forward(N2,offset+N2,data,F,sign,step);
for(int i = 0; i < N2; i+=2) {
// sin -x = -sin x
// cos -x = cos x
c = Math.cos(i*_PI/N2); // (_PI * 2 * k)
s = sign*Math.sin(i*_PI/N2); // (_PI * 2 * k)
// (c+si)*(p+qi) => (cp-sq, cq+sp)
real = F[i+N2+offset]*c + F[i+N2+1+offset]*s;
imaginary = F[i+N2+1+offset]*c - F[i+N2+offset]*s;
//System.out.println("pair: "+ (i+offset)/2+"-"+ (i+N2+offset)/2+" N2: " +N2+ " offset:"+offset + " step:"+step + " c="+c+" s="+s);
F[i+N2+offset] = F[i+offset] - real;
F[i+N2+1+offset] = F[i+1+offset] - imaginary;
F[i+offset] += real;
F[i+1+offset] += imaginary;
}
} else {
F[offset] = data[offset] + data[offset+N2];
F[offset+1] = data[offset+1] + data[offset+N2+1];
F[offset+N2] = data[offset] - data[offset+N2];
F[offset+N2+1] = data[offset+1] - data[offset+N2+1];
//System.out.println("pair: "+ (offset)/2+"-"+ (N2+offset)/2+" N2: " +N2+ " offset:"+offset + " FFT-2 " +F[offset]+" "+F[offset+N2]);
}
}
private static void inverse(int N, double[] data, double[] F) {
int N2 = N/2;
int inverseSign = -1;
forward(N,0,data,F,inverseSign,0);
for(int i = 0; i < N; i+=2) {
F[i] /= N2;
F[i+1] /= N2;
}
}
private static int reverseBit(int x, int maxBits) {
int bit = 0;
while (maxBits != 0){
bit <<=1;
bit |=( x &1 );
x >>=1;
maxBits>>=1;
}
return bit;
}
private static int resizeToPowerOfTwo(int x) {
int bit = 1;
while(x > bit) {
bit <<=1;
}
return bit;
}
private static void reverseData(int N, double[] data) {
for(int i = 2; i < N*2; i+=2) {
int rev = reverseBit(i/2,N-1) << 1;
if(i < rev) {
//System.out.println("i: "+i/2+" rev:"+rev/2);
double tmpR = data[i];
double tmpI = data[i+1];
data[i] = data[rev];
data[i+1] = data[rev+1];
data[rev] = tmpR;
data[rev+1] = tmpI;
}
}
}
public static void main(String[] args) throws Exception {
URL imgResource = FFTCPU2D.class.getResource(IMAGE_PATH);
String imgPath = Paths.get(imgResource.toURI()).toFile().getAbsolutePath();
System.out.println(imgPath);
img = new ImageIcon(imgResource);
width = img.getIconWidth();
height = img.getIconHeight();
BufferedImage bimg = new BufferedImage(width, height, BufferedImage.TYPE_3BYTE_BGR);
Graphics g = bimg.createGraphics();
img.paintIcon(null, g, 0,0);
g.dispose();
int pixelLength;
if(bimg.getAlphaRaster() != null) {
pixelLength = 4;
} else {
pixelLength = 3;
}
DataBufferByte dbb = (DataBufferByte)bimg.getRaster().getDataBuffer();
byte[] byteArray = dbb.getData();
System.out.println(byteArray.length);
System.out.println(pixelLength);
//生データ、バイト
byte[] _rawBlue = new byte[width*height];
byte[] _rawRed = new byte[width*height];
byte[] _rawGreen = new byte[width*height];
//RGBチャネルとして独立した変数に記憶
for( int i = 0; i < width; i++ ) {
for( int j = 0; j < height; j++ ) {
int index = j*width*pixelLength + i*pixelLength;
_rawBlue[j*width + i] = byteArray[index++];
_rawGreen[j*width + i] = byteArray[index++];
_rawRed[j*width + i] = byteArray[index];
}
}
//色彩配列、倍精度
double[] arBlue = new double[size*2];
double[] arGreen = new double[size*2];
double[] arRed = new double[size*2];
// Single row vector：行スキャンで走査線を記憶する変数
double[] srvBlue = new double[_RowSize*2];
double[] srvGreen = new double[_RowSize*2];
double[] srvRed = new double[_RowSize*2];
// Single row output vector：行を変換した配列を記憶する変数
double[] sroBlue = new double[_RowSize*2];
double[] sroGreen = new double[_RowSize*2];
double[] sroRed = new double[_RowSize*2];
// Single column vector：列スキャンを記憶する変数
double[] sclBlue = new double[_ColSize*2];
double[] sclGreen = new double[_ColSize*2];
double[] sclRed = new double[_ColSize*2];
// Single column output vector：列を変換した配列を記憶する変数
double[] scoBlue = new double[_ColSize*2];
double[] scoGreen = new double[_ColSize*2];
double[] scoRed = new double[_ColSize*2];
// Row matrix：行をフーリエ変換した配列を連結した行列
double[] rowBlue = new double[size*2];
double[] rowGreen = new double[size*2];
double[] rowRed = new double[size*2];
// Column matrix：列をフーリエ変換した配列を連結した行列
double[] colBlue = new double[size*2];
double[] colGreen = new double[size*2];
double[] colRed = new double[size*2];
// Output matrix：出力のための行列
double[] outBlue = new double[size*2];
double[] outGreen = new double[size*2];
double[] outRed = new double[size*2];
// cast to double；倍精度にキャスト
for(int i = 0; i < _ColSize; i++) {
for(int j = 0; j < _RowSize*2; j+=2) {
arBlue[i*_RowSize*2 + j] = (double)(_rawBlue[i*width + j/2] & 0xff);
arBlue[i*_RowSize*2 + j + 1] = (double)0.0f;
arGreen[i*_RowSize*2 + j] = (double)(_rawGreen[i*width + j/2] & 0xff);
arGreen[i*_RowSize*2 + j + 1] = (double)0.0f;
arRed[i*_RowSize*2 + j] = (double)(_rawRed[i*width + j/2] & 0xff);
arRed[i*_RowSize*2 + j + 1] = (double)0.0f;
}
}
for(int i = 0; i < _ColSize; i++) {
//get row vector：行スキャン
for(int j = 0; j < _RowSize*2; j+=2) {
srvBlue[j] = arBlue[i*_RowSize*2 + j];
srvBlue[j+1] = arBlue[i*_RowSize*2 + j + 1];
srvGreen[j] = arGreen[i*_RowSize*2 + j];
srvGreen[j+1] = arGreen[i*_RowSize*2 + j + 1];
srvRed[j] = arRed[i*_RowSize*2 + j];
srvRed[j+1] = arRed[i*_RowSize*2 + j + 1];
}
//bit reversal：ビット反転
reverseData(_RowSize,srvBlue);
reverseData(_RowSize,srvGreen);
reverseData(_RowSize,srvRed);
//forward transformation：フーリエ変換
forward(_RowSize*2, 0, srvBlue, sroBlue, _SIGN, 0);
forward(_RowSize*2, 0, srvGreen, sroGreen, _SIGN, 0);
forward(_RowSize*2, 0, srvRed, sroRed, _SIGN, 0);
// form the row matrix
for(int k = 0; k < _RowSize*2; k+=2) {
rowBlue[i*_RowSize*2 + k] = sroBlue[k];
rowBlue[i*_RowSize*2 + k + 1] = sroBlue[k+1];
rowGreen[i*_RowSize*2 + k] = sroGreen[k];
rowGreen[i*_RowSize*2 + k + 1] = sroGreen[k+1];
rowRed[i*_RowSize*2 + k] = sroRed[k];
rowRed[i*_RowSize*2 + k + 1] = sroRed[k+1];
}
}
for(int i = 0; i <_RowSize*2; i+=2) {
//get column vector：列スキャンをして配列に記憶
for(int j = 0; j < _ColSize*2; j+=2) {
sclBlue[j] = rowBlue[j*_RowSize + i];
sclBlue[j+1] = rowBlue[j*_RowSize + i + 1];
sclGreen[j] = rowGreen[j*_RowSize + i];
sclGreen[j+1] = rowGreen[j*_RowSize + i + 1];
sclRed[j] = rowRed[j*_RowSize + i];
sclRed[j+1] = rowRed[j*_RowSize + i + 1];
}
//bit reversal：ビット反転
reverseData(_ColSize,sclBlue);
reverseData(_ColSize,sclGreen);
reverseData(_ColSize,sclRed);
//forward transformation：フーリエ変換
forward(_ColSize*2, 0, sclBlue, scoBlue, _SIGN, 0);
forward(_ColSize*2, 0, sclGreen, scoGreen, _SIGN, 0);
forward(_ColSize*2, 0, sclRed, scoRed, _SIGN, 0);
for(int k = 0; k < _ColSize*2; k+=2) {
colBlue[k*_RowSize + i] = scoBlue[k];
colBlue[k*_RowSize + i + 1] = scoBlue[k+1];
colGreen[k*_RowSize + i] = scoGreen[k];
colGreen[k*_RowSize + i + 1] = scoGreen[k+1];
colRed[k*_RowSize + i] = scoRed[k];
colRed[k*_RowSize + i + 1] = scoRed[k+1];
}
}
/*
*
* Apply High Pass Filter：
*
*/
for(int i = 0; i < _RowSize*2; i+=2) {
for(int k = 0; k < _ColSize; k++) {
colBlue[k*_RowSize + i] = colBlue[k*_RowSize + i] > 256*256*200 ? colBlue[k*_RowSize + i] : 0;
colBlue[k*_RowSize + i + 1] = colBlue[k*_RowSize + i + 1] > 256*256*200 ? colBlue[k*_RowSize + i + 1] : 0;
colGreen[k*_RowSize + i] = colGreen[k*_RowSize + i] > 256*256*200 ? colGreen[k*_RowSize + i] : 0;
colGreen[k*_RowSize + i + 1] = colGreen[k*_RowSize + i + 1] > 256*256*200 ? colGreen[k*_RowSize + i + 1] : 0;
colRed[k*_RowSize + i] = colRed[k*_RowSize + i] > 256*256*200 ? colRed[k*_RowSize + i] : 0;
colRed[k*_RowSize + i + 1] = colRed[k*_RowSize + i + 1] > 256*256*200 ? colRed[k*_RowSize + i + 1] : 0;
}
}
/*
*
* Inverse FFT：逆フーリエ変換
*
*/
for(int i = 0; i <_RowSize*2; i+=2) {
//get column vector：列スキャン
for(int j = 0; j < _ColSize*2; j+=2) {
sclBlue[j] = colBlue[j*_RowSize + i];
sclBlue[j+1] = colBlue[j*_RowSize + i + 1];
sclGreen[j] = colGreen[j*_RowSize + i];
sclGreen[j+1] = colGreen[j*_RowSize + i + 1];
sclRed[j] = colRed[j*_RowSize + i];
sclRed[j+1] = colRed[j*_RowSize + i + 1];
}
//bit reversal：ビット反転
reverseData(_ColSize, sclBlue);
reverseData(_ColSize, sclGreen);
reverseData(_ColSize, sclRed);
//inverse transformation：逆変換
inverse(_ColSize*2, sclBlue, scoBlue);
inverse(_ColSize*2, sclGreen, scoGreen);
inverse(_ColSize*2, sclRed, scoRed);
for(int k = 0; k < _ColSize*2; k+=2) {
colBlue[k*_RowSize + i] = scoBlue[k];
colBlue[k*_RowSize + i + 1] = scoBlue[k+1];
colGreen[k*_RowSize + i] = scoGreen[k];
colGreen[k*_RowSize + i + 1] = scoGreen[k+1];
colRed[k*_RowSize + i] = scoRed[k];
colRed[k*_RowSize + i + 1] = scoRed[k+1];
}
}
for(int i = 0; i <_ColSize; i++) {
//get row vector：行スキャン
for(int j = 0; j < _RowSize*2; j+=2) {
srvBlue[j] = colBlue[i*_RowSize*2 + j];
srvBlue[j+1] = colBlue[i*_RowSize*2 + j + 1];
srvGreen[j] = colGreen[i*_RowSize*2 + j];
srvGreen[j+1] = colGreen[i*_RowSize*2 + j + 1];
srvRed[j] = colRed[i*_RowSize*2 + j];
srvRed[j+1] = colRed[i*_RowSize*2 + j + 1];
}
//bit reversal：ビット反転
reverseData(_RowSize, srvBlue);
reverseData(_RowSize, srvGreen);
reverseData(_RowSize, srvRed);
// inverse transformation：逆変換
inverse(_RowSize*2, srvBlue, sroBlue);
inverse(_RowSize*2, srvGreen, sroGreen);
inverse(_RowSize*2, srvRed, sroRed);
for(int k = 0; k < _RowSize*2; k+=2) {
outBlue[i*_RowSize*2 + k] = sroBlue[k];
outBlue[i*_RowSize*2 + k + 1] = sroBlue[k+1];
outGreen[i*_RowSize*2 + k] = sroGreen[k];
outGreen[i*_RowSize*2 + k + 1] = sroGreen[k+1];
outRed[i*_RowSize*2 + k] = sroRed[k];
outRed[i*_RowSize*2 + k + 1] = sroRed[k+1];
}
}
// ファイルに保存するため整数型の変数をメモリ空間に割り当て
int[] output = new int[_RowSize*_ColSize];
BufferedImage outImage = new BufferedImage(_RowSize, _ColSize, BufferedImage.TYPE_INT_RGB);
// RGB色彩を整数型に集約
for( int i = 0; i < _RowSize*2; i+=2 ) {
for( int j = 0; j < _ColSize; j++ ) {
output[j*_RowSize+i/2] = 0xFF000000
| ((int) ((byte)outBlue[j*_RowSize*2 + i] & 0xff) << 16)
| ((int) ((byte)outGreen[j*_RowSize*2 + i] & 0xff) << 8)
| ((int) ((byte)outRed[j*_RowSize*2 + i] & 0xff));
}
}
// クラスパス（mavenプロジェクトの場合は/target/classes/pathToPackage）にimage.pngファイルを生成
outImage.setRGB(0, 0, _RowSize, _ColSize, output, 0, _RowSize);
URL resource = FFTCPU2D.class.getResource("./");
String path = Paths.get(resource.toURI()).toFile().getAbsolutePath();
System.out.println(path+"/image.png");
ImageIO.write(outImage, "png", new File(path+"/image.png"));
}
}
