1 8a3b2ceb 2004-04-24 devnull #include <u.h>
2 8a3b2ceb 2004-04-24 devnull #include <libc.h>
3 8a3b2ceb 2004-04-24 devnull #include <bio.h>
4 8a3b2ceb 2004-04-24 devnull #include "sky.h"
6 8a3b2ceb 2004-04-24 devnull static void unshuffle(Pix*, int, int, Pix*);
7 8a3b2ceb 2004-04-24 devnull static void unshuffle1(Pix*, int, Pix*);
10 8a3b2ceb 2004-04-24 devnull hinv(Pix *a, int nx, int ny)
12 8a3b2ceb 2004-04-24 devnull int nmax, log2n, h0, hx, hy, hc, i, j, k;
13 8a3b2ceb 2004-04-24 devnull int nxtop, nytop, nxf, nyf, c;
14 8a3b2ceb 2004-04-24 devnull int oddx, oddy;
15 8a3b2ceb 2004-04-24 devnull int shift;
16 8a3b2ceb 2004-04-24 devnull int s10, s00;
17 8a3b2ceb 2004-04-24 devnull Pix *tmp;
20 8a3b2ceb 2004-04-24 devnull * log2n is log2 of max(nx, ny) rounded up to next power of 2
22 8a3b2ceb 2004-04-24 devnull nmax = ny;
23 8a3b2ceb 2004-04-24 devnull if(nx > nmax)
24 8a3b2ceb 2004-04-24 devnull nmax = nx;
25 8a3b2ceb 2004-04-24 devnull log2n = log(nmax)/LN2 + 0.5;
26 8a3b2ceb 2004-04-24 devnull if(nmax > (1<<log2n))
30 8a3b2ceb 2004-04-24 devnull * get temporary storage for shuffling elements
32 8a3b2ceb 2004-04-24 devnull tmp = (Pix*)malloc(((nmax+1)/2) * sizeof(*tmp));
33 8a3b2ceb 2004-04-24 devnull if(tmp == nil) {
34 8a3b2ceb 2004-04-24 devnull fprint(2, "hinv: insufficient memory\n");
35 8a3b2ceb 2004-04-24 devnull exits("memory");
39 8a3b2ceb 2004-04-24 devnull * do log2n expansions
41 8a3b2ceb 2004-04-24 devnull * We're indexing a as a 2-D array with dimensions (nx,ny).
43 8a3b2ceb 2004-04-24 devnull shift = 1;
44 8a3b2ceb 2004-04-24 devnull nxtop = 1;
45 8a3b2ceb 2004-04-24 devnull nytop = 1;
46 8a3b2ceb 2004-04-24 devnull nxf = nx;
47 8a3b2ceb 2004-04-24 devnull nyf = ny;
48 8a3b2ceb 2004-04-24 devnull c = 1<<log2n;
49 8a3b2ceb 2004-04-24 devnull for(k = log2n-1; k>=0; k--) {
51 8a3b2ceb 2004-04-24 devnull * this somewhat cryptic code generates the sequence
52 8a3b2ceb 2004-04-24 devnull * ntop[k-1] = (ntop[k]+1)/2, where ntop[log2n] = n
54 8a3b2ceb 2004-04-24 devnull c = c>>1;
55 8a3b2ceb 2004-04-24 devnull nxtop = nxtop<<1;
56 8a3b2ceb 2004-04-24 devnull nytop = nytop<<1;
57 8a3b2ceb 2004-04-24 devnull if(nxf <= c)
60 8a3b2ceb 2004-04-24 devnull nxf -= c;
61 8a3b2ceb 2004-04-24 devnull if(nyf <= c)
64 8a3b2ceb 2004-04-24 devnull nyf -= c;
67 8a3b2ceb 2004-04-24 devnull * halve divisors on last pass
69 8a3b2ceb 2004-04-24 devnull if(k == 0)
70 8a3b2ceb 2004-04-24 devnull shift = 0;
73 8a3b2ceb 2004-04-24 devnull * unshuffle in each dimension to interleave coefficients
75 8a3b2ceb 2004-04-24 devnull for(i = 0; i<nxtop; i++)
76 8a3b2ceb 2004-04-24 devnull unshuffle1(&a[ny*i], nytop, tmp);
77 8a3b2ceb 2004-04-24 devnull for(j = 0; j<nytop; j++)
78 8a3b2ceb 2004-04-24 devnull unshuffle(&a[j], nxtop, ny, tmp);
79 8a3b2ceb 2004-04-24 devnull oddx = nxtop & 1;
80 8a3b2ceb 2004-04-24 devnull oddy = nytop & 1;
81 8a3b2ceb 2004-04-24 devnull for(i = 0; i<nxtop-oddx; i += 2) {
82 8a3b2ceb 2004-04-24 devnull s00 = ny*i; /* s00 is index of a[i,j] */
83 8a3b2ceb 2004-04-24 devnull s10 = s00+ny; /* s10 is index of a[i+1,j] */
84 8a3b2ceb 2004-04-24 devnull for(j = 0; j<nytop-oddy; j += 2) {
86 8a3b2ceb 2004-04-24 devnull * Multiply h0,hx,hy,hc by 2 (1 the last time through).
88 8a3b2ceb 2004-04-24 devnull h0 = a[s00 ] << shift;
89 8a3b2ceb 2004-04-24 devnull hx = a[s10 ] << shift;
90 8a3b2ceb 2004-04-24 devnull hy = a[s00+1] << shift;
91 8a3b2ceb 2004-04-24 devnull hc = a[s10+1] << shift;
94 8a3b2ceb 2004-04-24 devnull * Divide sums by 4 (shift right 2 bits).
95 8a3b2ceb 2004-04-24 devnull * Add 1 to round -- note that these values are always
96 8a3b2ceb 2004-04-24 devnull * positive so we don't need to do anything special
97 8a3b2ceb 2004-04-24 devnull * for rounding negative numbers.
99 8a3b2ceb 2004-04-24 devnull a[s10+1] = (h0 + hx + hy + hc + 2) >> 2;
100 8a3b2ceb 2004-04-24 devnull a[s10 ] = (h0 + hx - hy - hc + 2) >> 2;
101 8a3b2ceb 2004-04-24 devnull a[s00+1] = (h0 - hx + hy - hc + 2) >> 2;
102 8a3b2ceb 2004-04-24 devnull a[s00 ] = (h0 - hx - hy + hc + 2) >> 2;
103 8a3b2ceb 2004-04-24 devnull s00 += 2;
104 8a3b2ceb 2004-04-24 devnull s10 += 2;
106 8a3b2ceb 2004-04-24 devnull if(oddy) {
108 8a3b2ceb 2004-04-24 devnull * do last element in row if row length is odd
109 8a3b2ceb 2004-04-24 devnull * s00+1, s10+1 are off edge
111 8a3b2ceb 2004-04-24 devnull h0 = a[s00 ] << shift;
112 8a3b2ceb 2004-04-24 devnull hx = a[s10 ] << shift;
113 8a3b2ceb 2004-04-24 devnull a[s10 ] = (h0 + hx + 2) >> 2;
114 8a3b2ceb 2004-04-24 devnull a[s00 ] = (h0 - hx + 2) >> 2;
117 8a3b2ceb 2004-04-24 devnull if(oddx) {
119 8a3b2ceb 2004-04-24 devnull * do last row if column length is odd
120 8a3b2ceb 2004-04-24 devnull * s10, s10+1 are off edge
122 8a3b2ceb 2004-04-24 devnull s00 = ny*i;
123 8a3b2ceb 2004-04-24 devnull for(j = 0; j<nytop-oddy; j += 2) {
124 8a3b2ceb 2004-04-24 devnull h0 = a[s00 ] << shift;
125 8a3b2ceb 2004-04-24 devnull hy = a[s00+1] << shift;
126 8a3b2ceb 2004-04-24 devnull a[s00+1] = (h0 + hy + 2) >> 2;
127 8a3b2ceb 2004-04-24 devnull a[s00 ] = (h0 - hy + 2) >> 2;
128 8a3b2ceb 2004-04-24 devnull s00 += 2;
130 8a3b2ceb 2004-04-24 devnull if(oddy) {
132 8a3b2ceb 2004-04-24 devnull * do corner element if both row and column lengths are odd
133 8a3b2ceb 2004-04-24 devnull * s00+1, s10, s10+1 are off edge
135 8a3b2ceb 2004-04-24 devnull h0 = a[s00 ] << shift;
136 8a3b2ceb 2004-04-24 devnull a[s00 ] = (h0 + 2) >> 2;
140 8a3b2ceb 2004-04-24 devnull free(tmp);
145 8a3b2ceb 2004-04-24 devnull unshuffle(Pix *a, int n, int n2, Pix *tmp)
148 8a3b2ceb 2004-04-24 devnull int nhalf, twon2, n2xnhalf;
149 8a3b2ceb 2004-04-24 devnull Pix *p1, *p2, *pt;
151 8a3b2ceb 2004-04-24 devnull twon2 = n2<<1;
152 8a3b2ceb 2004-04-24 devnull nhalf = (n+1)>>1;
153 8a3b2ceb 2004-04-24 devnull n2xnhalf = n2*nhalf; /* pointer to a[i] */
156 8a3b2ceb 2004-04-24 devnull * copy 2nd half of array to tmp
158 8a3b2ceb 2004-04-24 devnull pt = tmp;
159 8a3b2ceb 2004-04-24 devnull p1 = &a[n2xnhalf]; /* pointer to a[i] */
160 8a3b2ceb 2004-04-24 devnull for(i=nhalf; i<n; i++) {
161 8a3b2ceb 2004-04-24 devnull *pt = *p1;
163 8a3b2ceb 2004-04-24 devnull p1 += n2;
167 8a3b2ceb 2004-04-24 devnull * distribute 1st half of array to even elements
169 8a3b2ceb 2004-04-24 devnull p2 = &a[n2xnhalf]; /* pointer to a[i] */
170 8a3b2ceb 2004-04-24 devnull p1 = &a[n2xnhalf<<1]; /* pointer to a[2*i] */
171 8a3b2ceb 2004-04-24 devnull for(i=nhalf-1; i>=0; i--) {
172 8a3b2ceb 2004-04-24 devnull p1 -= twon2;
173 8a3b2ceb 2004-04-24 devnull p2 -= n2;
174 8a3b2ceb 2004-04-24 devnull *p1 = *p2;
178 8a3b2ceb 2004-04-24 devnull * now distribute 2nd half of array (in tmp) to odd elements
180 8a3b2ceb 2004-04-24 devnull pt = tmp;
181 8a3b2ceb 2004-04-24 devnull p1 = &a[n2]; /* pointer to a[i] */
182 8a3b2ceb 2004-04-24 devnull for(i=1; i<n; i+=2) {
183 8a3b2ceb 2004-04-24 devnull *p1 = *pt;
184 8a3b2ceb 2004-04-24 devnull p1 += twon2;
191 8a3b2ceb 2004-04-24 devnull unshuffle1(Pix *a, int n, Pix *tmp)
194 8a3b2ceb 2004-04-24 devnull int nhalf;
195 8a3b2ceb 2004-04-24 devnull Pix *p1, *p2, *pt;
197 8a3b2ceb 2004-04-24 devnull nhalf = (n+1) >> 1;
200 8a3b2ceb 2004-04-24 devnull * copy 2nd half of array to tmp
202 8a3b2ceb 2004-04-24 devnull pt = tmp;
203 8a3b2ceb 2004-04-24 devnull p1 = &a[nhalf]; /* pointer to a[i] */
204 8a3b2ceb 2004-04-24 devnull for(i=nhalf; i<n; i++) {
205 8a3b2ceb 2004-04-24 devnull *pt = *p1;
211 8a3b2ceb 2004-04-24 devnull * distribute 1st half of array to even elements
213 8a3b2ceb 2004-04-24 devnull p2 = &a[nhalf]; /* pointer to a[i] */
214 8a3b2ceb 2004-04-24 devnull p1 = &a[nhalf<<1]; /* pointer to a[2*i] */
215 8a3b2ceb 2004-04-24 devnull for(i=nhalf-1; i>=0; i--) {
216 8a3b2ceb 2004-04-24 devnull p1 -= 2;
218 8a3b2ceb 2004-04-24 devnull *p1 = *p2;
222 8a3b2ceb 2004-04-24 devnull * now distribute 2nd half of array (in tmp) to odd elements
224 8a3b2ceb 2004-04-24 devnull pt = tmp;
225 8a3b2ceb 2004-04-24 devnull p1 = &a[1]; /* pointer to a[i] */
226 8a3b2ceb 2004-04-24 devnull for(i=1; i<n; i+=2) {
227 8a3b2ceb 2004-04-24 devnull *p1 = *pt;
228 8a3b2ceb 2004-04-24 devnull p1 += 2;
