|
/Users/lyon/j4p/src/ip/hak/DHT2D.java
|
1 package ip.hak;
2
3 public class DHT2D {
4 public static void main(String args[]) {
5 System.out.println("Discrete Hartley Transform");
6 }
7
8 public static void DHT1D(float f[], int n) {
9 int p = (int) (Math.log(n) / Math.log(4));
10 int n4 = n / 4;
11 int r = 4;
12
13 int j = 1,i = 0;
14 do {
15 i++;
16 if (i < j) {
17 float t = f[j - 1];
18 f[j - 1] = f[i - 1];
19 f[i - 1] = t;
20 }
21 int k = n4;
22 while (3 * k < j) {
23 j = j - 3 * k;
24 k = k / 4;
25 }
26 j = j + k;
27
28 } while (i < n - 1);
29
30 for (i = 0; i < n; i += 4) {
31 float t1 = f[i] + f[i + 1];
32 float t2 = f[i] - f[i + 1];
33 float t3 = f[i + 2] + f[i + 3];
34 float t4 = f[i + 2] - f[i + 3];
35 f[i] = (float) (t1 + t3);
36 f[i + 1] = (float) (t1 - t3);
37 f[i + 2] = (float) (t2 + t4);
38 f[i + 3] = (float) (t2 - t4);
39 }
40
41
42 for (int l = 2; l <= p; l++) {
43 int e1 = (int) Math.pow(2, l + l - 3);
44 int e2 = e1 + e1;
45 int e3 = e2 + e1;
46 int e4 = e3 + e1;
47 int e5 = e4 + e1;
48 int e6 = e5 + e1;
49 int e7 = e6 + e1;
50 int e8 = e7 + e1;
51 for (j = 0; j < n; j += e8) {
52 float t1 = f[j] + f[j + e2];
53 float t2 = f[j] - f[j + e2];
54 float t3 = f[j + e4] + f[j + e6];
55 float t4 = f[j + e4] - f[j + e6];
56 f[j] = (float) (t1 + t3);
57 f[j + e2] = (float) (t1 - t3);
58 f[j + e4] = (float) (t2 + t4);
59 f[j + e6] = (float) (t2 - t4);
60 t1 = f[j + e1];
61 t2 = f[j + e3] * r;
62 t3 = f[j + e5];
63 t4 = f[j + e7] * r;
64 f[j + e1] = (float) (t1 + t2 + t3);
65 f[j + e3] = (float) (t1 - t3 + t4);
66 f[j + e5] = (float) (t1 - t2 + t3);
67 f[j + e7] = (float) (t1 - t3 - t4);
68
69 for (int k = 1; k < e1; k++) {
70 int l1 = j + k;
71 int l2 = l1 + e2;
72 int l3 = l1 + e4;
73 int l4 = l1 + e6;
74 int l5 = j + e2 - k;
75 int l6 = l5 + e2;
76 int l7 = l5 + e4;
77 int l8 = l5 + e6;
78 double a1 = Math.PI * k / e4;
79 double a2 = a1 + a1;
80 double a3 = a1 + a2;
81 double c1 = Math.cos(a1);
82 double c2 = Math.cos(a2);
83 double c3 = Math.cos(a3);
84 double s1 = Math.sin(a1);
85 double s2 = Math.sin(a2);
86 double s3 = Math.sin(a3);
87 double te5 = f[l2] * c1 + f[l6] * s1;
88 double te6 = f[l3] * c2 + f[l7] * s2;
89 double te7 = f[l4] * c3 + f[l8] * s3;
90 double te8 = f[l6] * c1 - f[l2] * s1;
91 double te9 = f[l7] * c2 - f[l3] * s2;
92 double te0 = f[l8] * c3 - f[l4] * s3;
93 double te1 = f[l5] - te9;
94 double te2 = f[l5] + te9;
95 double te3 = -te8 - te0;
96 double te4 = te5 - te7;
97 f[l5] = (float) (te1 + te4);
98 f[l6] = (float) (te2 + te3);
99 f[l7] = (float) (te1 - te4);
100 f[l8] = (float) (te2 - te3);
101 te1 = f[l1] + te6;
102 te2 = f[l1] - te6;
103 te3 = te8 - te0;
104 te4 = te5 + te7;
105 f[l1] = (float) (te1 + te4);
106 f[l2] = (float) (te2 + te3);
107 f[l3] = (float) (te1 - te4);
108 f[l4] = (float) (te2 - te3);
109 }
110 }
111 }
112 }
113
114 public static float[][] DHT2D(float a[][]) {
115 int n = a.length;
116 if (n != a[0].length) {
117 System.out.println("Image should be square!");
118 return null;
119 }
120 /*
121 // change [colum][row]-> [row][colum]
122 float f[][]=new float[n][n];
123
124 for (int x=0;x<n;x++)
125 for(int y=0;y<n;y++)
126 f[x][y]=a[y][x];
127 */
128
129 for (int r = 0; r < n; r++)
130 DHT1D(a[r], n);
131
132 transpose(a);
133
134 for (int r = 0; r < n; r++)
135 DHT1D(a[r], n);
136
137 //transpose(f);
138
139 float h[][] = new float[n][n];
140
141 for (int x = 0; x < n; x++)
142 for (int y = 0; y < n; y++)
143 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]));
144 return h;
145 }
146
147 public static void transpose(float a[][]) {
148 float temp[][] = new float[a.length][a[0].length];
149
150 for (int x = 0; x < a.length; x++)
151 for (int y = 0; y < a[0].length; y++)
152 temp[x][y] = a[y][x];
153 a = temp;
154 }
155
156 public static void normalize(float a[][], int t) {
157 for (int x = 0; x < a.length; x++)
158 for (int y = 0; y < a[0].length; y++)
159 a[x][y] /= t;
160 }
161
162 public static float[][] forwardDHT2D(float a[][]) {
163 float f[][] = DHT2D(a);
164 normalize(f, f.length * f[0].length);
165
166 return f;
167 }
168
169 public static float[][] inverseDHT2D(float a[][]) {
170 return DHT2D(a);
171 }
172 }
173