1 /////////////////////////////////////////////////////////////////////////// 2 // 3 // Copyright (c) 2002, Industrial Light & Magic, a division of Lucas 4 // Digital Ltd. LLC 5 // 6 // All rights reserved. 7 // 8 // Redistribution and use in source and binary forms, with or without 9 // modification, are permitted provided that the following conditions are 10 // met: 11 // * Redistributions of source code must retain the above copyright 12 // notice, this list of conditions and the following disclaimer. 13 // * Redistributions in binary form must reproduce the above 14 // copyright notice, this list of conditions and the following disclaimer 15 // in the documentation and/or other materials provided with the 16 // distribution. 17 // * Neither the name of Industrial Light & Magic nor the names of 18 // its contributors may be used to endorse or promote products derived 19 // from this software without specific prior written permission. 20 // 21 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 22 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 23 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 24 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 25 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 26 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 27 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 28 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 29 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 30 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 31 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 32 // 33 /////////////////////////////////////////////////////////////////////////// 34 35 36 37 //----------------------------------------------------------------------------- 38 // 39 // 16-bit Haar Wavelet encoding and decoding 40 // 41 // The source code in this file is derived from the encoding 42 // and decoding routines written by Christian Rouet for his 43 // PIZ image file format. 44 // 45 //----------------------------------------------------------------------------- 46 47 48 #include <ImfWav.h> 49 50 namespace Imf { 51 namespace { 52 53 54 // 55 // Wavelet basis functions without modulo arithmetic; they produce 56 // the best compression ratios when the wavelet-transformed data are 57 // Huffman-encoded, but the wavelet transform works only for 14-bit 58 // data (untransformed data values must be less than (1 << 14)). 59 // 60 61 inline void 62 wenc14 (unsigned short a, unsigned short b, 63 unsigned short &l, unsigned short &h) 64 { 65 short as = a; 66 short bs = b; 67 68 short ms = (as + bs) >> 1; 69 short ds = as - bs; 70 71 l = ms; 72 h = ds; 73 } 74 75 76 inline void 77 wdec14 (unsigned short l, unsigned short h, 78 unsigned short &a, unsigned short &b) 79 { 80 short ls = l; 81 short hs = h; 82 83 int hi = hs; 84 int ai = ls + (hi & 1) + (hi >> 1); 85 86 short as = ai; 87 short bs = ai - hi; 88 89 a = as; 90 b = bs; 91 } 92 93 94 // 95 // Wavelet basis functions with modulo arithmetic; they work with full 96 // 16-bit data, but Huffman-encoding the wavelet-transformed data doesn't 97 // compress the data quite as well. 98 // 99 100 const int NBITS = 16; 101 const int A_OFFSET = 1 << (NBITS - 1); 102 const int M_OFFSET = 1 << (NBITS - 1); 103 const int MOD_MASK = (1 << NBITS) - 1; 104 105 106 inline void 107 wenc16 (unsigned short a, unsigned short b, 108 unsigned short &l, unsigned short &h) 109 { 110 int ao = (a + A_OFFSET) & MOD_MASK; 111 int m = ((ao + b) >> 1); 112 int d = ao - b; 113 114 if (d < 0) 115 m = (m + M_OFFSET) & MOD_MASK; 116 117 d &= MOD_MASK; 118 119 l = m; 120 h = d; 121 } 122 123 124 inline void 125 wdec16 (unsigned short l, unsigned short h, 126 unsigned short &a, unsigned short &b) 127 { 128 int m = l; 129 int d = h; 130 int bb = (m - (d >> 1)) & MOD_MASK; 131 int aa = (d + bb - A_OFFSET) & MOD_MASK; 132 b = bb; 133 a = aa; 134 } 135 136 } // namespace 137 138 139 // 140 // 2D Wavelet encoding: 141 // 142 143 void 144 wav2Encode 145 (unsigned short* in, // io: values are transformed in place 146 int nx, // i : x size 147 int ox, // i : x offset 148 int ny, // i : y size 149 int oy, // i : y offset 150 unsigned short mx) // i : maximum in[x][y] value 151 { 152 bool w14 = (mx < (1 << 14)); 153 int n = (nx > ny)? ny: nx; 154 int p = 1; // == 1 << level 155 int p2 = 2; // == 1 << (level+1) 156 157 // 158 // Hierachical loop on smaller dimension n 159 // 160 161 while (p2 <= n) 162 { 163 unsigned short *py = in; 164 unsigned short *ey = in + oy * (ny - p2); 165 int oy1 = oy * p; 166 int oy2 = oy * p2; 167 int ox1 = ox * p; 168 int ox2 = ox * p2; 169 unsigned short i00,i01,i10,i11; 170 171 // 172 // Y loop 173 // 174 175 for (; py <= ey; py += oy2) 176 { 177 unsigned short *px = py; 178 unsigned short *ex = py + ox * (nx - p2); 179 180 // 181 // X loop 182 // 183 184 for (; px <= ex; px += ox2) 185 { 186 unsigned short *p01 = px + ox1; 187 unsigned short *p10 = px + oy1; 188 unsigned short *p11 = p10 + ox1; 189 190 // 191 // 2D wavelet encoding 192 // 193 194 if (w14) 195 { 196 wenc14 (*px, *p01, i00, i01); 197 wenc14 (*p10, *p11, i10, i11); 198 wenc14 (i00, i10, *px, *p10); 199 wenc14 (i01, i11, *p01, *p11); 200 } 201 else 202 { 203 wenc16 (*px, *p01, i00, i01); 204 wenc16 (*p10, *p11, i10, i11); 205 wenc16 (i00, i10, *px, *p10); 206 wenc16 (i01, i11, *p01, *p11); 207 } 208 } 209 210 // 211 // Encode (1D) odd column (still in Y loop) 212 // 213 214 if (nx & p) 215 { 216 unsigned short *p10 = px + oy1; 217 218 if (w14) 219 wenc14 (*px, *p10, i00, *p10); 220 else 221 wenc16 (*px, *p10, i00, *p10); 222 223 *px= i00; 224 } 225 } 226 227 // 228 // Encode (1D) odd line (must loop in X) 229 // 230 231 if (ny & p) 232 { 233 unsigned short *px = py; 234 unsigned short *ex = py + ox * (nx - p2); 235 236 for (; px <= ex; px += ox2) 237 { 238 unsigned short *p01 = px + ox1; 239 240 if (w14) 241 wenc14 (*px, *p01, i00, *p01); 242 else 243 wenc16 (*px, *p01, i00, *p01); 244 245 *px= i00; 246 } 247 } 248 249 // 250 // Next level 251 // 252 253 p = p2; 254 p2 <<= 1; 255 } 256 } 257 258 259 // 260 // 2D Wavelet decoding: 261 // 262 263 void 264 wav2Decode 265 (unsigned short* in, // io: values are transformed in place 266 int nx, // i : x size 267 int ox, // i : x offset 268 int ny, // i : y size 269 int oy, // i : y offset 270 unsigned short mx) // i : maximum in[x][y] value 271 { 272 bool w14 = (mx < (1 << 14)); 273 int n = (nx > ny)? ny: nx; 274 int p = 1; 275 int p2; 276 277 // 278 // Search max level 279 // 280 281 while (p <= n) 282 p <<= 1; 283 284 p >>= 1; 285 p2 = p; 286 p >>= 1; 287 288 // 289 // Hierarchical loop on smaller dimension n 290 // 291 292 while (p >= 1) 293 { 294 unsigned short *py = in; 295 unsigned short *ey = in + oy * (ny - p2); 296 int oy1 = oy * p; 297 int oy2 = oy * p2; 298 int ox1 = ox * p; 299 int ox2 = ox * p2; 300 unsigned short i00,i01,i10,i11; 301 302 // 303 // Y loop 304 // 305 306 for (; py <= ey; py += oy2) 307 { 308 unsigned short *px = py; 309 unsigned short *ex = py + ox * (nx - p2); 310 311 // 312 // X loop 313 // 314 315 for (; px <= ex; px += ox2) 316 { 317 unsigned short *p01 = px + ox1; 318 unsigned short *p10 = px + oy1; 319 unsigned short *p11 = p10 + ox1; 320 321 // 322 // 2D wavelet decoding 323 // 324 325 if (w14) 326 { 327 wdec14 (*px, *p10, i00, i10); 328 wdec14 (*p01, *p11, i01, i11); 329 wdec14 (i00, i01, *px, *p01); 330 wdec14 (i10, i11, *p10, *p11); 331 } 332 else 333 { 334 wdec16 (*px, *p10, i00, i10); 335 wdec16 (*p01, *p11, i01, i11); 336 wdec16 (i00, i01, *px, *p01); 337 wdec16 (i10, i11, *p10, *p11); 338 } 339 } 340 341 // 342 // Decode (1D) odd column (still in Y loop) 343 // 344 345 if (nx & p) 346 { 347 unsigned short *p10 = px + oy1; 348 349 if (w14) 350 wdec14 (*px, *p10, i00, *p10); 351 else 352 wdec16 (*px, *p10, i00, *p10); 353 354 *px= i00; 355 } 356 } 357 358 // 359 // Decode (1D) odd line (must loop in X) 360 // 361 362 if (ny & p) 363 { 364 unsigned short *px = py; 365 unsigned short *ex = py + ox * (nx - p2); 366 367 for (; px <= ex; px += ox2) 368 { 369 unsigned short *p01 = px + ox1; 370 371 if (w14) 372 wdec14 (*px, *p01, i00, *p01); 373 else 374 wdec16 (*px, *p01, i00, *p01); 375 376 *px= i00; 377 } 378 } 379 380 // 381 // Next level 382 // 383 384 p2 = p; 385 p >>= 1; 386 } 387 } 388 389 390 } // namespace Imf 391
