package transforms;
import math.Mat;
import ip.vs.ColorUtils;
import ip.vs.ImageUtils;
import utils.Timer;
public class FFT {
private final int FORWARD_FFT = -1;
private final int REVERSE_FFT = 1;
private float direction = (float) FORWARD_FFT;
static final float twoPI = (float) (2 * Math.PI);
private int N;
private int numBits;
private int width, height;
private float minPSD = 9999999;
private float maxPSD = -9999999;
ColorUtils CU = new ColorUtils();
ImageUtils iUtils = new ImageUtils();
public boolean DisplayLogPSD = false;
public float cR_r[][]; public float cG_r[][]; public float cB_r[][]; public float cR_i[][]; public float cG_i[][]; public float cB_i[][];
float r_data[] = null;
float i_data[] = null;
public void printStats() {
Mat.printStats("cR_r", cR_r);
Mat.printStats("cG_r", cG_r);
Mat.printStats("cB_r", cB_r);
Mat.printStats("cR_i", cR_i);
Mat.printStats("cG_i", cG_i);
Mat.printStats("cB_i", cB_i);
}
public void normalize() {
float sf = 1.0f / (cR_r.length * cR_r[0].length);
Mat.scale(cR_r, sf);
Mat.scale(cG_r, sf);
Mat.scale(cB_r, sf);
Mat.scale(cR_i, sf);
Mat.scale(cG_i, sf);
Mat.scale(cB_i, sf);
}
public void complexMult(FFT fft) {
for (int x = 0; x < cR_r.length; x++)
for (int y = 0; y < cR_r[0].length; y++) {
cR_r[x][y] = cR_r[x][y] * fft.cR_r[x][y];
cG_r[x][y] = cG_r[x][y] * fft.cG_r[x][y];
cB_r[x][y] = cB_r[x][y] * fft.cB_r[x][y];
cR_i[x][y] = cR_i[x][y] * fft.cR_i[x][y];
cG_i[x][y] = cG_i[x][y] * fft.cG_i[x][y];
cB_i[x][y] = cB_i[x][y] * fft.cB_i[x][y];
}
}
private void init() {
minPSD = 9999999;
maxPSD = -9999999;
}
public int[] getPsd() {
float[] magnitudeR = magnitudeSpectrum(cR_r, cR_i);
float[] magnitudeG = magnitudeSpectrum(cG_r, cG_i);
float[] magnitudeB = magnitudeSpectrum(cB_r, cB_i);
double scaleFactor = 100;
System.out.println("Max psd = " + maxPSD);
scaleFactor = 255.0 / (Math.log(1 + maxPSD));
for (int i = 0; i < N; i++) {
magnitudeR[i] = (float)
(scaleFactor * Math.log(1 + magnitudeR[i]));
magnitudeG[i] = (float)
(scaleFactor * Math.log(1 + magnitudeG[i]));
magnitudeB[i] = (float)
(scaleFactor * Math.log(1 + magnitudeB[i]));
}
return (CU.imagetoInt(magnitudeR, magnitudeG, magnitudeB));
}
public int[] getPhaseImage() {
float[] magnitudeR = getPhaseImage(cR_r, cR_i);
float[] magnitudeG = getPhaseImage(cG_r, cG_i);
float[] magnitudeB = getPhaseImage(cB_r, cB_i);
double scaleFactor = 100;
System.out.println("Max psd = " + maxPSD);
scaleFactor = 255;
System.out.println("Scalefactor = " + scaleFactor);
for (int i = 0; i < N; i++) {
magnitudeR[i] = (float)
(scaleFactor * Math.log(1 + magnitudeR[i]));
magnitudeG[i] = (float)
(scaleFactor * Math.log(1 + magnitudeG[i]));
magnitudeB[i] = (float)
(scaleFactor * Math.log(1 + magnitudeB[i]));
}
System.out.println("Minimum PSD value: " + minPSD);
return (CU.imagetoInt(magnitudeR, magnitudeG, magnitudeB));
}
public int[] forward2dFFT(float[] imageData_R, float[] imageData_G,
float[] imageData_B, int imageWidth, int imageHeight) {
init();
width = imageWidth;
height = imageHeight;
cR_r = new float[height][width];
cR_i = new float[height][width];
cG_r = new float[height][width];
cG_i = new float[height][width];
cB_r = new float[height][width];
cB_i = new float[height][width];
N = width * height;
numBits = log2(width);
copyFloatToComplex(cR_r, cR_i, imageData_R);
copyFloatToComplex(cG_r, cG_i, imageData_G);
copyFloatToComplex(cB_r, cB_i, imageData_B);
for (int i = 0; i < height; i++) {
forwardFFT(cR_r[i], cR_i[i]);
forwardFFT(cG_r[i], cG_i[i]);
forwardFFT(cB_r[i], cB_i[i]);
}
cR_r = Rotate90(cR_r);
cR_i = Rotate90(cR_i);
cG_r = Rotate90(cG_r);
cG_i = Rotate90(cG_i);
cB_r = Rotate90(cB_r);
cB_i = Rotate90(cB_i);
for (int i = 0; i < height; i++) {
forwardFFT(cR_r[i], cR_i[i]);
forwardFFT(cG_r[i], cG_i[i]);
forwardFFT(cB_r[i], cB_i[i]);
}
return getPsd();
}
public int[] forward2dFFT(short[] imageData_R, short[] imageData_G,
short[] imageData_B, int imageWidth, int imageHeight) {
init();
width = imageWidth;
height = imageHeight;
cR_r = new float[height][width];
cR_i = new float[height][width];
cG_r = new float[height][width];
cG_i = new float[height][width];
cB_r = new float[height][width];
cB_i = new float[height][width];
N = width * height;
numBits = log2(width);
copyShortToComplex(cR_r, cR_i, imageData_R);
copyShortToComplex(cG_r, cG_i, imageData_G);
copyShortToComplex(cB_r, cB_i, imageData_B);
for (int i = 0; i < height; i++) {
forwardFFT(cR_r[i], cR_i[i]);
forwardFFT(cG_r[i], cG_i[i]);
forwardFFT(cB_r[i], cB_i[i]);
}
cR_r = Rotate90(cR_r);
cR_i = Rotate90(cR_i);
cG_r = Rotate90(cG_r);
cG_i = Rotate90(cG_i);
cB_r = Rotate90(cB_r);
cB_i = Rotate90(cB_i);
for (int i = 0; i < height; i++) {
forwardFFT(cR_r[i], cR_i[i]);
forwardFFT(cG_r[i], cG_i[i]);
forwardFFT(cB_r[i], cB_i[i]);
}
return getPsd();
}
int log2(double n) {
return (int) (Math.log(n) / Math.log(2));
}
public int[] reverse2dFFT() {
init();
N = width * height;
numBits = log2(width);
for (int i = 0; i < height; i++) {
reverseFFT(cR_r[i], cR_i[i]);
reverseFFT(cG_r[i], cG_i[i]);
reverseFFT(cB_r[i], cB_i[i]);
}
cR_r = Rotate90(cR_r);
cR_i = Rotate90(cR_i);
cG_r = Rotate90(cG_r);
cG_i = Rotate90(cG_i);
cB_r = Rotate90(cB_r);
cB_i = Rotate90(cB_i);
for (int i = 0; i < height; i++) {
reverseFFT(cR_r[i], cR_i[i]);
reverseFFT(cG_r[i], cG_i[i]);
reverseFFT(cB_r[i], cB_i[i]);
}
rotateInPlace180(cR_r);
rotateInPlace180(cR_i);
rotateInPlace180(cG_r);
rotateInPlace180(cG_i);
rotateInPlace180(cB_r);
rotateInPlace180(cB_i);
float[] realR = magnitudeSpectrum(cR_r, cR_i);
float[] realG = magnitudeSpectrum(cG_r, cG_i);
float[] realB = magnitudeSpectrum(cB_r, cB_i);
return (CU.imagetoInt(realR, realG, realB));
}
public void TestFFT() {
init();
System.out.println("Testing FFT engine in 2D class...");
numBits = 3;
N = 8;
float[] in_r = new float[N];
float[] in_i = new float[N];
float[] holder_r = new float[N];
float[] holder_i = new float[N];
for (int i = 0; i < N; i++) {
in_r[i] = (float) i / N;
in_i[i] = (float) 0;
}
forwardFFT(in_r, in_i);
for (int i = 0; i < N; i++) {
holder_r[i] = in_r[i];
holder_i[i] = in_i[i];
}
for (int i = 0; i < N; i++) {
in_r[i] = in_r[i];
in_i[i] = in_i[i];
}
reverseFFT(in_r, in_i);
float[] mag = magnitudeSpectrum(in_r, in_i);
System.out.println("i x1[i] re[i] im[i] tv[i]");
for (int i = 0; i < N; i++) {
System.out.println(i + "\t" + i + "\t" + holder_r[i] + "\t\t" +
holder_i[i] + "\t\t" + mag[i]);
}
in_r = new float[262144];
in_i = new float[262144];
N = 262144;
numBits = 18;
for (int i = 0; i < N; i++) {
in_r[i] = (float) i / N;
in_i[i] = (float) 0;
}
Timer t1 = new Timer();
t1.clear();
t1.mark();
forwardFFT(in_r, in_i);
System.out.println("256K point 1D FFT (using float):");
t1.record();
t1.report();
}
public void swap(int i) {
float tempr;
int j = bitr(i);
String js = Integer.toBinaryString(j);
String is = Integer.toBinaryString(i);
tempr = r_data[j];
r_data[j] = r_data[i];
r_data[i] = tempr;
tempr = i_data[j];
i_data[j] = i_data[i];
i_data[i] = tempr;
}
public void swapInt(int i, int j) {
float tempr;
int ti;
int tj;
ti = i - 1;
tj = j - 1;
tempr = r_data[tj];
r_data[tj] = r_data[ti];
r_data[ti] = tempr;
tempr = i_data[tj];
i_data[tj] = i_data[ti];
i_data[ti] = tempr;
}
double getMaxValue(double in[]) {
double max;
max = -0.99e30;
for (int i = 0; i < in.length; i++)
if (in[i] > max)
max = in[i];
return max;
}
void bitReverse2() {
int n = r_data.length;
int j = 1;
int k;
for (int i = 1; i < n; i++) {
if (i < j) {
swapInt(i, j);
}
k = n / 2;
while (k >= 1 && k < j) {
j = j - k;
k = k / 2;
}
j = j + k;
} }
void bitReverse() {
int n = r_data.length;
int j = 1;
int k;
for (int i = 1; i < n; i++) {
if (i < j) {
swap(i);
} j = bitr(i);
} }
int bitr(int j) {
int ans = 0;
for (int i = 0; i < numBits; i++) {
ans = (ans << 1) + (j & 1);
j = j >> 1;
}
return ans;
}
public void forwardFFT(float[] in_r, float[] in_i) {
direction = FORWARD_FFT;
fft(in_r, in_i);
}
public void reverseFFT(float[] in_r, float[] in_i) {
direction = REVERSE_FFT;
fft(in_r, in_i);
}
public void fft(float in_r[], float in_i[]) {
int id;
int localN;
float wtemp, Wjk_r, Wjk_i, Wj_r, Wj_i;
float theta, tempr, tempi;
int ti, tj;
int length = 1 << numBits; int n = length;
r_data = in_r;
i_data = in_i;
bitReverse2();
for (int m = 1; m <= numBits; m++) {
localN = 1 << m;
Wjk_r = 1;
Wjk_i = 0;
theta = twoPI / localN;
Wj_r = (float) Math.cos(theta);
Wj_i = (float) (direction * Math.sin(theta));
int nby2 = localN / 2;
for (int j = 0; j < nby2; j++) {
for (int k = j; k < n; k += localN) {
id = k + nby2;
tempr = Wjk_r * r_data[id] - Wjk_i * i_data[id];
tempi = Wjk_r * i_data[id] + Wjk_i * r_data[id];
r_data[id] = r_data[k] - tempr;
i_data[id] = i_data[k] - tempi;
r_data[k] += tempr;
i_data[k] += tempi;
}
wtemp = Wjk_r;
Wjk_r = Wj_r * Wjk_r - Wj_i * Wjk_i;
Wjk_i = Wj_r * Wjk_i + Wj_i * wtemp;
}
}
}
public float[][] getRedReal() {
return (cR_r);
}
public float[][] getRedImaginary() {
return (cR_i);
}
public float[][] getGreenReal() {
return (cG_r);
}
public float[][] getGreenImaginary() {
return (cG_i);
}
public float[][] getBlueReal() {
return (cB_r);
}
public float[][] getBlueImaginary() {
return (cB_i);
}
private void copy2dArray(float[][] dst, float[][] src) {
for (int i = 0; i < height; i++) {
for (int j = 0; j < width; j++) {
dst[i][j] = src[i][j];
}
}
}
private void copyFloatToComplex(float[][] dst_r, float[][] dst_i,
float[] imageData) {
int k = 0;
float alternateSign = 1;
for (int i = 0; i < height; i++) {
for (int j = 0; j < width; j++) {
alternateSign = ((i + j) % 2 == 0) ? -1 : 1;
dst_r[i][j] = (float) (imageData[k++] * alternateSign / N);
dst_i[i][j] = (float) 0.0;
}
}
}
private void copyShortToComplex(float[][] dst_r, float[][] dst_i,
short[] imageData) {
int k = 0;
float alternateSign = 1;
for (int i = 0; i < height; i++)
for (int j = 0; j < width; j++) {
alternateSign = ((i + j) % 2 == 0) ? -1 : 1;
dst_r[i][j] = (float) (imageData[k++] * alternateSign / N);
dst_i[i][j] = (float) 0.0;
}
}
private boolean isNotInRange(int x, int low, int high) {
if (x < low) return true;
if (x >= high) return true;
return false;
}
private float[][] Rotate90(float[][] in) {
float[][] out = new float[height][width];
for (int i = 0; i < height; i++) {
for (int j = 0; j < width; j++) {
int x = height - j - 1;
int y = i;
if (isNotInRange(x, 0, in.length)) continue;
if (isNotInRange(y, 0, in[0].length)) continue;
out[i][j] = in[x][y];
}
}
return (out);
}
private void rotateInPlace180(float[][] in) {
float temp;
for (int i = 0; i < height / 2; i++) {
for (int j = 0; j < width; j++) {
temp = in[i][j];
in[i][j] = in[height - i - 1][width - j - 1];
in[height - i - 1][width - j - 1] = temp;
}
}
}
private float[] copyRealToFloat(float[][] in_r) {
float[] f_data = new float[N];
int k = 0;
for (int i = 0; i < height; i++) {
for (int j = 0; j < width; j++) {
f_data[k++] = in_r[i][j];
}
}
return (f_data);
}
private float[] magnitudeSpectrum(float[][] in_r, float[][] in_i) {
float[] mag = new float[N];
int k = 0;
for (int i = 0; i < height; i++) {
for (int j = 0; j < width; j++) {
mag[k] = (float) Math.sqrt(in_r[i][j] * in_r[i][j] +
in_i[i][j] * in_i[i][j]);
if (minPSD > mag[k])
minPSD = mag[k];
if (maxPSD < mag[k])
maxPSD = mag[k];
k++;
}
}
return (mag);
}
private float[] getPhaseImage(float[][] in_r, float[][] in_i) {
float phase[] = new float[N];
int k = 0;
for (int i = 0; i < height; i++)
for (int j = 0; j < width; j++) {
phase[k] = (float) (in_r[i][j] / in_i[i][j]);
if (minPSD > phase[k])
minPSD = phase[k];
if (maxPSD < phase[k])
maxPSD = phase[k];
k++;
}
return (phase);
}
public float[] magnitudeSpectrum(float[] in_r, float[] in_i) {
N = in_r.length;
float[] mag = new float[N];
for (int i = 0; i < N; i++) {
mag[i] = (float) Math.sqrt(in_r[i] * in_r[i] +
in_i[i] * in_i[i]);
if (minPSD > mag[i]) {
minPSD = mag[i];
}
}
return (mag);
}
}