w3d/wavelet_xform.c

   1 #include "mem.h"
   2 #include <assert.h>
   3 #include "wavelet.h"
   4 #include "w3dtypes.h"
   5
   6
   7 static
   8 void fwd_analyze_1 (const TYPE *x, TYPE *d, int stride, int n)
   9 {
  10    int i, k=n/2;
  11
  12    for (i=0; i<k; i++)
  13       d[i] = x[(2*i+1)*stride] - x[2*i*stride];
  14 }
  15
  16
  17 static
  18 void fwd_synthesize_1 (const TYPE *x, TYPE *s, const TYPE *d, int stride, int n)
  19 {
  20    int i, k=n/2;
  21
  22    for (i=0; i<k; i++)
  23       s[i*stride] = x[2*i*stride] + d[i] / 2;
  24    if (n & 1)
  25       s[k*stride] = x[2*k*stride] + d[k-1] / 2;
  26 }
  27
  28
  29 static
  30 void inv_analyze_1 (TYPE *x, const TYPE *d, int stride, int n)
  31 {
  32    int i, k=n/2;
  33
  34    for (i=0; i<k; i++)
  35       x[(2*i+1)*stride] = d[i] + x[2*i*stride];
  36 }
  37
  38
  39 static
  40 void inv_synthesize_1 (TYPE *x, const TYPE *s, const TYPE *d, int stride, int n)
  41 {
  42    int i, k=n/2;
  43
  44    for (i=0; i<k; i++)
  45       x[2*i*stride] = s[i] - d[i] / 2;
  46    if (n & 1)
  47       x[2*k*stride] = s[k] - d[k-1] / 2;
  48 }
  49
  50
  51
  52 static
  53 void fwd_analyze_2 (const TYPE *x, TYPE *d, int stride, int n)
  54 {
  55    int i, k=n/2;
  56
  57    if (n & 1) {
  58       for (i=0; i<k; i++)
  59          d[i] = x[(2*i+1)*stride] - (x[2*i*stride] + x[(2*i+2)*stride]) / 2;
  60    } else {
  61       for (i=0; i<k-1; i++)
  62          d[i] = x[(2*i+1)*stride] - (x[2*i*stride] + x[(2*i+2)*stride]) / 2;
  63       d[k-1] = x[(n-1)*stride] - x[(n-2)*stride];
  64    }
  65 }
  66
  67
  68
  69 static
  70 void fwd_synthesize_2 (const TYPE *x, TYPE *s, const TYPE *d, int stride, int n)
  71 {
  72    int i, k=n/2;
  73
  74    s[0] = x[0] + d[1] / 2;
  75    for (i=1; i<k; i++)
  76       s[i*stride] = x[2*i*stride] + (d[i-1] + d[i]) / 4;
  77    if (n & 1)
  78       s[k*stride] = x[2*k*stride] + d[k-1] / 2;
  79 }
  80
  81
  82 static inline
  83 void inv_analyze_2 (TYPE *x, const TYPE *d, int stride, int n)
  84 {
  85    int i, k=n/2;
  86
  87    if (n & 1) {
  88       for (i=0; i<k; i++)
  89          x[(2*i+1)*stride] = d[i] + (x[2*i*stride] + x[(2*i+2)*stride]) / 2;
  90    } else {
  91       for (i=0; i<k-1; i++)
  92          x[(2*i+1)*stride] = d[i] + (x[2*i*stride] + x[(2*i+2)*stride]) / 2;
  93       x[(n-1)*stride] = d[k-1] + x[(n-2)*stride];
  94    }
  95 }
  96
  97
  98 static inline
  99 void inv_synthesize_2 (TYPE *x, const TYPE *s, const TYPE *d, int stride, int n)
 100 {
 101    int i, k=n/2;
 102
 103    x[0] = s[0] - d[1] / 2;
 104    for (i=1; i<k; i++)
 105       x[2*i*stride] = s[i] - (d[i-1] + d[i]) / 4;
 106    if (n & 1)
 107       x[2*k*stride] = s[k] - d[k-1] / 2;
 108 }
 109
 110
 111
 112 static
 113 void fwd_analyze_4 (const TYPE *x, TYPE *d, int stride, int n)
 114 {
 115    int i, k=n/2;
 116
 117    d[0] = x[stride] - (x[0] + x[2*stride]) / 2;
 118
 119    if (n & 1) {
 120       for (i=1; i<k-1; i++)
 121          d[i] = x[(2*i+1)*stride]
 122                - ((uint32_t) 9 * (x[2*i*stride] + x[(2*i+2)*stride])
 123                                - (x[(2*i-2)*stride] + x[(2*i+4)*stride])) / 16;
 124       if (k > 1)
 125          d[k-1] = x[(2*k-1)*stride] - (x[(2*k-2)*stride] + x[2*k*stride]) / 2;
 126    } else {
 127       for (i=1; i<k-2; i++)
 128          d[i] = x[(2*i+1)*stride]
 129                - ((uint32_t) 9 * (x[2*i*stride] + x[(2*i+2)*stride])
 130                                - (x[(2*i-2)*stride] + x[(2*i+4)*stride])) / 16;
 131       if (k > 2)
 132          d[k-2] = x[(2*k-3)*stride] - (x[(2*k-4)*stride]
 133                                      + x[(2*k-2)*stride]) / 2;
 134       if (k > 1)
 135          d[k-1] = x[(n-1)*stride] - x[(n-2)*stride];
 136    }
 137 }
 138
 139
 140 static
 141 void fwd_synthesize_4 (const TYPE *x, TYPE *s, const TYPE *d, int stride, int n)
 142 {
 143    int i, k=n/2;
 144
 145    s[0] = x[0] + d[1] / 2;
 146    if (k > 1)
 147       s[stride] = x[2*stride] + (d[0] + d[1]) / 4;
 148    for (i=2; i<k-1; i++)
 149       s[i*stride] = x[2*i*stride]
 150                    + ((uint32_t) 9 * (d[i-1] + d[i]) - (d[i-2] + d[i+1])) / 32;
 151    if (k > 2)
 152       s[(k-1)*stride] = x[(2*k-2)*stride] + (d[k-2] + d[k-1]) / 4;
 153    if (n & 1)
 154       s[k*stride] = x[2*k*stride] + d[k-1] / 2;
 155 }
 156
 157
 158 static
 159 void inv_analyze_4 (TYPE *x, const TYPE *d, int stride, int n)
 160 {
 161    int i, k=n/2;
 162
 163    x[stride] = d[0] + (x[0] + x[2*stride]) / 2;
 164
 165    if (n & 1) {
 166       for (i=1; i<k-1; i++)
 167          x[(2*i+1)*stride] = d[i]
 168                + ((uint32_t) 9 * (x[2*i*stride] + x[(2*i+2)*stride])
 169                                - (x[(2*i-2)*stride] + x[(2*i+4)*stride])) / 16;
 170       if (k > 1)
 171          x[(2*k-1)*stride] = d[k-1] + (x[(2*k-2)*stride] + x[2*k*stride]) / 2;
 172    } else {
 173       for (i=1; i<k-2; i++)
 174          x[(2*i+1)*stride] = d[i]
 175                + (9 * (x[2*i*stride] + x[(2*i+2)*stride])
 176                     - (x[(2*i-2)*stride] + x[(2*i+4)*stride])) / 16;
 177       if (k > 2)
 178          x[(2*k-3)*stride] = d[k-2] + (x[(2*k-4)*stride]
 179                                      + x[(2*k-2)*stride]) / 2;
 180       if (k > 1)
 181          x[(n-1)*stride] = d[k-1] + x[(n-2)*stride];
 182    }
 183 }
 184
 185
 186 static
 187 void inv_synthesize_4 (TYPE *x, const TYPE *s, const TYPE *d, int stride, int n)
 188 {
 189    int i, k=n/2;
 190
 191    x[0] = s[0] - d[1] / 2;
 192    if (k > 1)
 193       x[2*stride] = s[1] - (d[0] + d[1]) / 4;
 194    for (i=2; i<k-1; i++)
 195       x[2*i*stride] = s[i] - ((uint32_t) 9 * (d[i-1] + d[i])
 196                                             - (d[i-2] + d[i+1])) / 32;
 197    if (k > 2)
 198       x[(2*k-2)*stride] = s[k-1] - (d[k-2] + d[k-1]) / 4;
 199    if (n & 1)
 200       x[2*k*stride] = s[k] - d[k-1] / 2;
 201 }
 202
 203
 204 static inline
 205 void copyback_d (TYPE *x, const TYPE *d, int stride, int n)
 206 {
 207    int i, j, k=n/2;
 208
 209    for (i=n-k, j=0; i<n; i++, j++)
 210       x [i*stride] = d[j];
 211 }
 212
 213
 214 static inline
 215 void copy_s_d (const TYPE *x, TYPE *s_d, int stride, int n)
 216 {
 217    int i;
 218
 219    for (i=0; i<n; i++)
 220       s_d[i] = x [i*stride];
 221 }
 222
 223
 224
 225 typedef
 226 void (*FwdSFnc) (const TYPE *x, TYPE *s, const TYPE *d, int stride, int n);
 227
 228 typedef
 229 void (*FwdAFnc) (const TYPE *x, TYPE *d, int stride, int n);
 230
 231 typedef
 232 void (*InvSFnc) (TYPE *x, const TYPE *s, const TYPE *d, int stride, int n);
 233
 234 typedef
 235 void (*InvAFnc) (TYPE *x, const TYPE *d, int stride, int n);
 236
 237
 238
 239 static FwdSFnc fwd_synthesize [] = { NULL, fwd_synthesize_1, fwd_synthesize_2,
 240                                      NULL, fwd_synthesize_4 };
 241
 242 static FwdAFnc fwd_analyze [] = { NULL, fwd_analyze_1, fwd_analyze_2,
 243                                   NULL, fwd_analyze_4 };
 244
 245 static InvSFnc inv_synthesize [] = { NULL, inv_synthesize_1, inv_synthesize_2,
 246                                      NULL, inv_synthesize_4 };
 247
 248 static InvAFnc inv_analyze [] = { NULL, inv_analyze_1, inv_analyze_2,
 249                                   NULL, inv_analyze_4 };
 250
 251
 252
 253 static inline
 254 void fwd_xform (TYPE *scratchbuf, TYPE *data, int stride, int n,
 255                 int a_moments, int s_moments)
 256 {
 257    TYPE *x = data;
 258    TYPE *d = scratchbuf;
 259    TYPE *s = data;
 260
 261    assert (a_moments == 1 || a_moments == 2 || a_moments == 4);
 262    assert (s_moments == 1 || s_moments == 2 || s_moments == 4);
 263
 264    /*  XXX FIXME: Ugly hack to work around  */
 265    /*             the short-row bug in high */
 266    /*             order xform functions     */
 267    if (n < 9)
 268       a_moments = s_moments = 2;
 269    if (n < 5)
 270       a_moments = s_moments = 1;
 271
 272    fwd_analyze [a_moments] (x, d, stride, n);
 273    fwd_synthesize [s_moments] (x, s, d, stride, n);
 274    copyback_d (x, d, stride, n);
 275 }
 276
 277
 278 static inline
 279 void inv_xform (TYPE *scratchbuf, TYPE *data, int stride, int n,
 280                 int a_moments, int s_moments)
 281 {
 282    int k=n/2;
 283    TYPE *x = data;
 284    TYPE *s = scratchbuf;
 285    TYPE *d = scratchbuf + n - k;
 286
 287    assert (a_moments == 1 || a_moments == 2 || a_moments == 4);
 288    assert (s_moments == 1 || s_moments == 2 || s_moments == 4);
 289
 290    /*  XXX FIXME: Ugly hack to work around  */
 291    /*             the short-row bug in high */
 292    /*             order xform functions     */
 293    if (n < 9)
 294       a_moments = s_moments = 2;
 295    if (n < 5)
 296       a_moments = s_moments = 1;
 297
 298    copy_s_d (data, scratchbuf, stride, n);
 299    inv_synthesize [s_moments] (x, s, d, stride, n);
 300    inv_analyze [a_moments] (x, d, stride, n);
 301 }
 302
 303
 304
 305 void wavelet_3d_buf_fwd_xform (Wavelet3DBuf* buf, int a_moments, int s_moments)
 306 {
 307    int level;
 308
 309    for (level=buf->scales-1; level>0; level--) {
 310       uint32_t w = buf->w[level];
 311       uint32_t h = buf->h[level];
 312       uint32_t f = buf->f[level];
 313
 314       if (w > 1) {
 315          int row, frame;
 316          for (frame=0; frame<f; frame++) {
 317             for (row=0; row<h; row++) {
 318                TYPE *data = buf->data + (frame * buf->height + row) * buf->width;
 319                fwd_xform (buf->scratchbuf, data, 1, w, a_moments, s_moments);
 320             }
 321          }
 322       }
 323
 324       if (h > 1) {
 325          int col, frame;
 326          for (frame=0; frame<f; frame++) {
 327             for (col=0; col<w; col++) {
 328                TYPE *data = buf->data + frame * buf->width * buf->height + col;
 329                fwd_xform (buf->scratchbuf, data, buf->width, h,
 330                           a_moments, s_moments);
 331             }
 332          }
 333       }
 334
 335       if (f > 1) {
 336          int i, j;
 337          for (j=0; j<h; j++) {
 338             for (i=0; i<w; i++) {
 339                TYPE *data = buf->data + j*buf->width + i;
 340                fwd_xform (buf->scratchbuf, data, buf->width * buf->height, f,
 341                           a_moments, s_moments);
 342             }
 343          }
 344       }
 345    }
 346 }
 347
 348
 349 void wavelet_3d_buf_inv_xform (Wavelet3DBuf* buf, int a_moments, int s_moments)
 350 {
 351    int level;
 352
 353    for (level=1; level<buf->scales; level++) {
 354       uint32_t w = buf->w[level];
 355       uint32_t h = buf->h[level];
 356       uint32_t f = buf->f[level];
 357
 358       if (f > 1) {
 359          int i, j;
 360          for (j=0; j<h; j++) {
 361             for (i=0; i<w; i++) {
 362                TYPE *data = buf->data + j*buf->width + i;
 363                inv_xform (buf->scratchbuf, data, buf->width * buf->height, f,
 364                           a_moments, s_moments);
 365             }
 366          }
 367       }
 368
 369       if (h > 1) {
 370          int col, frame;
 371          for (frame=0; frame<f; frame++) {
 372             for (col=0; col<w; col++) {
 373                TYPE *data = buf->data + frame * buf->width * buf->height + col;
 374                inv_xform (buf->scratchbuf, data, buf->width, h,
 375                           a_moments, s_moments);
 376             }
 377          }
 378       }
 379
 380       if (w > 1) {
 381          int row, frame;
 382          for (frame=0; frame<f; frame++) {
 383             for (row=0; row<h; row++) {
 384                TYPE *data = buf->data + (frame * buf->height + row) * buf->width;
 385                inv_xform (buf->scratchbuf, data, 1, w, a_moments, s_moments);
 386             }
 387          }
 388       }
 389    }
 390 }
 391