6 static void unshuffle(Pix*, int, int, Pix*);
7 static void unshuffle1(Pix*, int, Pix*);
10 hinv(Pix *a, int nx, int ny)
12 int nmax, log2n, h0, hx, hy, hc, i, j, k;
13 int nxtop, nytop, nxf, nyf, c;
20 * log2n is log2 of max(nx, ny) rounded up to next power of 2
25 log2n = log(nmax)/LN2 + 0.5;
30 * get temporary storage for shuffling elements
32 tmp = (Pix*)malloc(((nmax+1)/2) * sizeof(*tmp));
34 fprint(2, "hinv: insufficient memory\n");
41 * We're indexing a as a 2-D array with dimensions (nx,ny).
49 for(k = log2n-1; k>=0; k--) {
51 * this somewhat cryptic code generates the sequence
52 * ntop[k-1] = (ntop[k]+1)/2, where ntop[log2n] = n
67 * halve divisors on last pass
73 * unshuffle in each dimension to interleave coefficients
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);
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) {
86 * Multiply h0,hx,hy,hc by 2 (1 the last time through).
88 h0 = a[s00 ] << shift;
89 hx = a[s10 ] << shift;
90 hy = a[s00+1] << shift;
91 hc = a[s10+1] << shift;
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.
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;
108 * do last element in row if row length is odd
109 * s00+1, s10+1 are off edge
111 h0 = a[s00 ] << shift;
112 hx = a[s10 ] << shift;
113 a[s10 ] = (h0 + hx + 2) >> 2;
114 a[s00 ] = (h0 - hx + 2) >> 2;
119 * do last row if column length is odd
120 * s10, s10+1 are off edge
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;
132 * do corner element if both row and column lengths are odd
133 * s00+1, s10, s10+1 are off edge
135 h0 = a[s00 ] << shift;
136 a[s00 ] = (h0 + 2) >> 2;
145 unshuffle(Pix *a, int n, int n2, Pix *tmp)
148 int nhalf, twon2, n2xnhalf;
153 n2xnhalf = n2*nhalf; /* pointer to a[i] */
156 * copy 2nd half of array to tmp
159 p1 = &a[n2xnhalf]; /* pointer to a[i] */
160 for(i=nhalf; i<n; i++) {
167 * distribute 1st half of array to even elements
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--) {
178 * now distribute 2nd half of array (in tmp) to odd elements
181 p1 = &a[n2]; /* pointer to a[i] */
182 for(i=1; i<n; i+=2) {
191 unshuffle1(Pix *a, int n, Pix *tmp)
200 * copy 2nd half of array to tmp
203 p1 = &a[nhalf]; /* pointer to a[i] */
204 for(i=nhalf; i<n; i++) {
211 * distribute 1st half of array to even elements
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--) {
222 * now distribute 2nd half of array (in tmp) to odd elements
225 p1 = &a[1]; /* pointer to a[i] */
226 for(i=1; i<n; i+=2) {