1 #include <u.h>
2 #include <libc.h>
3 #include <bio.h>
4 #include "sky.h"
5
6 static void unshuffle(Pix*, int, int, Pix*);
7 static void unshuffle1(Pix*, int, Pix*);
8
9 void
hinv(Pix * a,int nx,int ny)10 hinv(Pix *a, int nx, int ny)
11 {
12 int nmax, log2n, h0, hx, hy, hc, i, j, k;
13 int nxtop, nytop, nxf, nyf, c;
14 int oddx, oddy;
15 int shift;
16 int s10, s00;
17 Pix *tmp;
18
19 /*
20 * log2n is log2 of max(nx, ny) rounded up to next power of 2
21 */
22 nmax = ny;
23 if(nx > nmax)
24 nmax = nx;
25 log2n = log(nmax)/LN2 + 0.5;
26 if(nmax > (1<<log2n))
27 log2n++;
28
29 /*
30 * get temporary storage for shuffling elements
31 */
32 tmp = (Pix*)malloc(((nmax+1)/2) * sizeof(*tmp));
33 if(tmp == nil) {
34 fprint(2, "hinv: insufficient memory\n");
35 exits("memory");
36 }
37
38 /*
39 * do log2n expansions
40 *
41 * We're indexing a as a 2-D array with dimensions (nx,ny).
42 */
43 shift = 1;
44 nxtop = 1;
45 nytop = 1;
46 nxf = nx;
47 nyf = ny;
48 c = 1<<log2n;
49 for(k = log2n-1; k>=0; k--) {
50 /*
51 * this somewhat cryptic code generates the sequence
52 * ntop[k-1] = (ntop[k]+1)/2, where ntop[log2n] = n
53 */
54 c = c>>1;
55 nxtop = nxtop<<1;
56 nytop = nytop<<1;
57 if(nxf <= c)
58 nxtop--;
59 else
60 nxf -= c;
61 if(nyf <= c)
62 nytop--;
63 else
64 nyf -= c;
65
66 /*
67 * halve divisors on last pass
68 */
69 if(k == 0)
70 shift = 0;
71
72 /*
73 * unshuffle in each dimension to interleave coefficients
74 */
75 for(i = 0; i<nxtop; i++)
76 unshuffle1(&a[ny*i], nytop, tmp);
77 for(j = 0; j<nytop; j++)
78 unshuffle(&a[j], nxtop, ny, tmp);
79 oddx = nxtop & 1;
80 oddy = nytop & 1;
81 for(i = 0; i<nxtop-oddx; i += 2) {
82 s00 = ny*i; /* s00 is index of a[i,j] */
83 s10 = s00+ny; /* s10 is index of a[i+1,j] */
84 for(j = 0; j<nytop-oddy; j += 2) {
85 /*
86 * Multiply h0,hx,hy,hc by 2 (1 the last time through).
87 */
88 h0 = a[s00 ] << shift;
89 hx = a[s10 ] << shift;
90 hy = a[s00+1] << shift;
91 hc = a[s10+1] << shift;
92
93 /*
94 * Divide sums by 4 (shift right 2 bits).
95 * Add 1 to round -- note that these values are always
96 * positive so we don't need to do anything special
97 * for rounding negative numbers.
98 */
99 a[s10+1] = (h0 + hx + hy + hc + 2) >> 2;
100 a[s10 ] = (h0 + hx - hy - hc + 2) >> 2;
101 a[s00+1] = (h0 - hx + hy - hc + 2) >> 2;
102 a[s00 ] = (h0 - hx - hy + hc + 2) >> 2;
103 s00 += 2;
104 s10 += 2;
105 }
106 if(oddy) {
107 /*
108 * do last element in row if row length is odd
109 * s00+1, s10+1 are off edge
110 */
111 h0 = a[s00 ] << shift;
112 hx = a[s10 ] << shift;
113 a[s10 ] = (h0 + hx + 2) >> 2;
114 a[s00 ] = (h0 - hx + 2) >> 2;
115 }
116 }
117 if(oddx) {
118 /*
119 * do last row if column length is odd
120 * s10, s10+1 are off edge
121 */
122 s00 = ny*i;
123 for(j = 0; j<nytop-oddy; j += 2) {
124 h0 = a[s00 ] << shift;
125 hy = a[s00+1] << shift;
126 a[s00+1] = (h0 + hy + 2) >> 2;
127 a[s00 ] = (h0 - hy + 2) >> 2;
128 s00 += 2;
129 }
130 if(oddy) {
131 /*
132 * do corner element if both row and column lengths are odd
133 * s00+1, s10, s10+1 are off edge
134 */
135 h0 = a[s00 ] << shift;
136 a[s00 ] = (h0 + 2) >> 2;
137 }
138 }
139 }
140 free(tmp);
141 }
142
143 static
144 void
unshuffle(Pix * a,int n,int n2,Pix * tmp)145 unshuffle(Pix *a, int n, int n2, Pix *tmp)
146 {
147 int i;
148 int nhalf, twon2, n2xnhalf;
149 Pix *p1, *p2, *pt;
150
151 twon2 = n2<<1;
152 nhalf = (n+1)>>1;
153 n2xnhalf = n2*nhalf; /* pointer to a[i] */
154
155 /*
156 * copy 2nd half of array to tmp
157 */
158 pt = tmp;
159 p1 = &a[n2xnhalf]; /* pointer to a[i] */
160 for(i=nhalf; i<n; i++) {
161 *pt = *p1;
162 pt++;
163 p1 += n2;
164 }
165
166 /*
167 * distribute 1st half of array to even elements
168 */
169 p2 = &a[n2xnhalf]; /* pointer to a[i] */
170 p1 = &a[n2xnhalf<<1]; /* pointer to a[2*i] */
171 for(i=nhalf-1; i>=0; i--) {
172 p1 -= twon2;
173 p2 -= n2;
174 *p1 = *p2;
175 }
176
177 /*
178 * now distribute 2nd half of array (in tmp) to odd elements
179 */
180 pt = tmp;
181 p1 = &a[n2]; /* pointer to a[i] */
182 for(i=1; i<n; i+=2) {
183 *p1 = *pt;
184 p1 += twon2;
185 pt++;
186 }
187 }
188
189 static
190 void
unshuffle1(Pix * a,int n,Pix * tmp)191 unshuffle1(Pix *a, int n, Pix *tmp)
192 {
193 int i;
194 int nhalf;
195 Pix *p1, *p2, *pt;
196
197 nhalf = (n+1) >> 1;
198
199 /*
200 * copy 2nd half of array to tmp
201 */
202 pt = tmp;
203 p1 = &a[nhalf]; /* pointer to a[i] */
204 for(i=nhalf; i<n; i++) {
205 *pt = *p1;
206 pt++;
207 p1++;
208 }
209
210 /*
211 * distribute 1st half of array to even elements
212 */
213 p2 = &a[nhalf]; /* pointer to a[i] */
214 p1 = &a[nhalf<<1]; /* pointer to a[2*i] */
215 for(i=nhalf-1; i>=0; i--) {
216 p1 -= 2;
217 p2--;
218 *p1 = *p2;
219 }
220
221 /*
222 * now distribute 2nd half of array (in tmp) to odd elements
223 */
224 pt = tmp;
225 p1 = &a[1]; /* pointer to a[i] */
226 for(i=1; i<n; i+=2) {
227 *p1 = *pt;
228 p1 += 2;
229 pt++;
230 }
231 }
232