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
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
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
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 }