Subversion Repositories planix.SVN

Rev

Rev 2 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 - 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 
}