hinv.c (4447B)
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 }