package ip.hak;
public class DHT2D {
public static void main(String args[]) {
System.out.println("Discrete Hartley Transform");
}
public static void DHT1D(float f[], int n) {
int p = (int) (Math.log(n) / Math.log(4));
int n4 = n / 4;
int r = 4;
int j = 1,i = 0;
do {
i++;
if (i < j) {
float t = f[j - 1];
f[j - 1] = f[i - 1];
f[i - 1] = t;
}
int k = n4;
while (3 * k < j) {
j = j - 3 * k;
k = k / 4;
}
j = j + k;
} while (i < n - 1);
for (i = 0; i < n; i += 4) {
float t1 = f[i] + f[i + 1];
float t2 = f[i] - f[i + 1];
float t3 = f[i + 2] + f[i + 3];
float t4 = f[i + 2] - f[i + 3];
f[i] = (float) (t1 + t3);
f[i + 1] = (float) (t1 - t3);
f[i + 2] = (float) (t2 + t4);
f[i + 3] = (float) (t2 - t4);
}
for (int l = 2; l <= p; l++) {
int e1 = (int) Math.pow(2, l + l - 3);
int e2 = e1 + e1;
int e3 = e2 + e1;
int e4 = e3 + e1;
int e5 = e4 + e1;
int e6 = e5 + e1;
int e7 = e6 + e1;
int e8 = e7 + e1;
for (j = 0; j < n; j += e8) {
float t1 = f[j] + f[j + e2];
float t2 = f[j] - f[j + e2];
float t3 = f[j + e4] + f[j + e6];
float t4 = f[j + e4] - f[j + e6];
f[j] = (float) (t1 + t3);
f[j + e2] = (float) (t1 - t3);
f[j + e4] = (float) (t2 + t4);
f[j + e6] = (float) (t2 - t4);
t1 = f[j + e1];
t2 = f[j + e3] * r;
t3 = f[j + e5];
t4 = f[j + e7] * r;
f[j + e1] = (float) (t1 + t2 + t3);
f[j + e3] = (float) (t1 - t3 + t4);
f[j + e5] = (float) (t1 - t2 + t3);
f[j + e7] = (float) (t1 - t3 - t4);
for (int k = 1; k < e1; k++) {
int l1 = j + k;
int l2 = l1 + e2;
int l3 = l1 + e4;
int l4 = l1 + e6;
int l5 = j + e2 - k;
int l6 = l5 + e2;
int l7 = l5 + e4;
int l8 = l5 + e6;
double a1 = Math.PI * k / e4;
double a2 = a1 + a1;
double a3 = a1 + a2;
double c1 = Math.cos(a1);
double c2 = Math.cos(a2);
double c3 = Math.cos(a3);
double s1 = Math.sin(a1);
double s2 = Math.sin(a2);
double s3 = Math.sin(a3);
double te5 = f[l2] * c1 + f[l6] * s1;
double te6 = f[l3] * c2 + f[l7] * s2;
double te7 = f[l4] * c3 + f[l8] * s3;
double te8 = f[l6] * c1 - f[l2] * s1;
double te9 = f[l7] * c2 - f[l3] * s2;
double te0 = f[l8] * c3 - f[l4] * s3;
double te1 = f[l5] - te9;
double te2 = f[l5] + te9;
double te3 = -te8 - te0;
double te4 = te5 - te7;
f[l5] = (float) (te1 + te4);
f[l6] = (float) (te2 + te3);
f[l7] = (float) (te1 - te4);
f[l8] = (float) (te2 - te3);
te1 = f[l1] + te6;
te2 = f[l1] - te6;
te3 = te8 - te0;
te4 = te5 + te7;
f[l1] = (float) (te1 + te4);
f[l2] = (float) (te2 + te3);
f[l3] = (float) (te1 - te4);
f[l4] = (float) (te2 - te3);
}
}
}
}
public static float[][] DHT2D(float a[][]) {
int n = a.length;
if (n != a[0].length) {
System.out.println("Image should be square!");
return null;
}
for (int r = 0; r < n; r++)
DHT1D(a[r], n);
transpose(a);
for (int r = 0; r < n; r++)
DHT1D(a[r], n);
float h[][] = new float[n][n];
for (int x = 0; x < n; x++)
for (int y = 0; y < n; y++)
h[x][y] = (float) ((1f / 2f) * (a[x][y] + a[n - x - 1][y] + a[x][n - y - 1] + a[n - x - 1][n - y - 1]));
return h;
}
public static void transpose(float a[][]) {
float temp[][] = new float[a.length][a[0].length];
for (int x = 0; x < a.length; x++)
for (int y = 0; y < a[0].length; y++)
temp[x][y] = a[y][x];
a = temp;
}
public static void normalize(float a[][], int t) {
for (int x = 0; x < a.length; x++)
for (int y = 0; y < a[0].length; y++)
a[x][y] /= t;
}
public static float[][] forwardDHT2D(float a[][]) {
float f[][] = DHT2D(a);
normalize(f, f.length * f[0].length);
return f;
}
public static float[][] inverseDHT2D(float a[][]) {
return DHT2D(a);
}
}