sys/src/cmd/scat/hinv.c

   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 }