ref: 0b5e782882152dae80fbcc6b54e187b1926e93ec
dir: /sys/src/cmd/scat/hinv.c/
#include <u.h> #include <libc.h> #include <bio.h> #include "sky.h" static void unshuffle(Pix*, int, int, Pix*); static void unshuffle1(Pix*, int, Pix*); void hinv(Pix *a, int nx, int ny) { int nmax, log2n, h0, hx, hy, hc, i, j, k; int nxtop, nytop, nxf, nyf, c; int oddx, oddy; int shift; int s10, s00; Pix *tmp; /* * log2n is log2 of max(nx, ny) rounded up to next power of 2 */ nmax = ny; if(nx > nmax) nmax = nx; log2n = log(nmax)/LN2 + 0.5; if(nmax > (1<<log2n)) log2n++; /* * get temporary storage for shuffling elements */ tmp = (Pix*)malloc(((nmax+1)/2) * sizeof(*tmp)); if(tmp == nil) { fprint(2, "hinv: insufficient memory\n"); exits("memory"); } /* * do log2n expansions * * We're indexing a as a 2-D array with dimensions (nx,ny). */ shift = 1; nxtop = 1; nytop = 1; nxf = nx; nyf = ny; c = 1<<log2n; for(k = log2n-1; k>=0; k--) { /* * this somewhat cryptic code generates the sequence * ntop[k-1] = (ntop[k]+1)/2, where ntop[log2n] = n */ c = c>>1; nxtop = nxtop<<1; nytop = nytop<<1; if(nxf <= c) nxtop--; else nxf -= c; if(nyf <= c) nytop--; else nyf -= c; /* * halve divisors on last pass */ if(k == 0) shift = 0; /* * unshuffle in each dimension to interleave coefficients */ for(i = 0; i<nxtop; i++) unshuffle1(&a[ny*i], nytop, tmp); for(j = 0; j<nytop; j++) unshuffle(&a[j], nxtop, ny, tmp); oddx = nxtop & 1; oddy = nytop & 1; for(i = 0; i<nxtop-oddx; i += 2) { s00 = ny*i; /* s00 is index of a[i,j] */ s10 = s00+ny; /* s10 is index of a[i+1,j] */ for(j = 0; j<nytop-oddy; j += 2) { /* * Multiply h0,hx,hy,hc by 2 (1 the last time through). */ h0 = a[s00 ] << shift; hx = a[s10 ] << shift; hy = a[s00+1] << shift; hc = a[s10+1] << shift; /* * Divide sums by 4 (shift right 2 bits). * Add 1 to round -- note that these values are always * positive so we don't need to do anything special * for rounding negative numbers. */ a[s10+1] = (h0 + hx + hy + hc + 2) >> 2; a[s10 ] = (h0 + hx - hy - hc + 2) >> 2; a[s00+1] = (h0 - hx + hy - hc + 2) >> 2; a[s00 ] = (h0 - hx - hy + hc + 2) >> 2; s00 += 2; s10 += 2; } if(oddy) { /* * do last element in row if row length is odd * s00+1, s10+1 are off edge */ h0 = a[s00 ] << shift; hx = a[s10 ] << shift; a[s10 ] = (h0 + hx + 2) >> 2; a[s00 ] = (h0 - hx + 2) >> 2; } } if(oddx) { /* * do last row if column length is odd * s10, s10+1 are off edge */ s00 = ny*i; for(j = 0; j<nytop-oddy; j += 2) { h0 = a[s00 ] << shift; hy = a[s00+1] << shift; a[s00+1] = (h0 + hy + 2) >> 2; a[s00 ] = (h0 - hy + 2) >> 2; s00 += 2; } if(oddy) { /* * do corner element if both row and column lengths are odd * s00+1, s10, s10+1 are off edge */ h0 = a[s00 ] << shift; a[s00 ] = (h0 + 2) >> 2; } } } free(tmp); } static void unshuffle(Pix *a, int n, int n2, Pix *tmp) { int i; int nhalf, twon2, n2xnhalf; Pix *p1, *p2, *pt; twon2 = n2<<1; nhalf = (n+1)>>1; n2xnhalf = n2*nhalf; /* pointer to a[i] */ /* * copy 2nd half of array to tmp */ pt = tmp; p1 = &a[n2xnhalf]; /* pointer to a[i] */ for(i=nhalf; i<n; i++) { *pt = *p1; pt++; p1 += n2; } /* * distribute 1st half of array to even elements */ p2 = &a[n2xnhalf]; /* pointer to a[i] */ p1 = &a[n2xnhalf<<1]; /* pointer to a[2*i] */ for(i=nhalf-1; i>=0; i--) { p1 -= twon2; p2 -= n2; *p1 = *p2; } /* * now distribute 2nd half of array (in tmp) to odd elements */ pt = tmp; p1 = &a[n2]; /* pointer to a[i] */ for(i=1; i<n; i+=2) { *p1 = *pt; p1 += twon2; pt++; } } static void unshuffle1(Pix *a, int n, Pix *tmp) { int i; int nhalf; Pix *p1, *p2, *pt; nhalf = (n+1) >> 1; /* * copy 2nd half of array to tmp */ pt = tmp; p1 = &a[nhalf]; /* pointer to a[i] */ for(i=nhalf; i<n; i++) { *pt = *p1; pt++; p1++; } /* * distribute 1st half of array to even elements */ p2 = &a[nhalf]; /* pointer to a[i] */ p1 = &a[nhalf<<1]; /* pointer to a[2*i] */ for(i=nhalf-1; i>=0; i--) { p1 -= 2; p2--; *p1 = *p2; } /* * now distribute 2nd half of array (in tmp) to odd elements */ pt = tmp; p1 = &a[1]; /* pointer to a[i] */ for(i=1; i<n; i+=2) { *p1 = *pt; p1 += 2; pt++; } }