1235368Sgnn/////////////////////////////////////////////////////////////////////////// 2235368Sgnn// 3235368Sgnn// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas 4235368Sgnn// Digital Ltd. LLC 5235368Sgnn// 6235368Sgnn// All rights reserved. 7235368Sgnn// 8235368Sgnn// Redistribution and use in source and binary forms, with or without 9235368Sgnn// modification, are permitted provided that the following conditions are 10235368Sgnn// met: 11235368Sgnn// * Redistributions of source code must retain the above copyright 12235368Sgnn// notice, this list of conditions and the following disclaimer. 13235368Sgnn// * Redistributions in binary form must reproduce the above 14235368Sgnn// copyright notice, this list of conditions and the following disclaimer 15235368Sgnn// in the documentation and/or other materials provided with the 16235368Sgnn// distribution. 17235368Sgnn// * Neither the name of Industrial Light & Magic nor the names of 18235368Sgnn// its contributors may be used to endorse or promote products derived 19235368Sgnn// from this software without specific prior written permission. 20235368Sgnn// 21235368Sgnn// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 22235368Sgnn// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 23235368Sgnn// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 24235368Sgnn// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 25235368Sgnn// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 26235368Sgnn// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 27235368Sgnn// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 28235368Sgnn// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 29235368Sgnn// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 30235368Sgnn// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 31235368Sgnn// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 32235368Sgnn// 33235368Sgnn/////////////////////////////////////////////////////////////////////////// 34235368Sgnn 35235368Sgnn 36235368Sgnn 37235368Sgnn//----------------------------------------------------------------------------- 38235368Sgnn// 39235368Sgnn// 16-bit Haar Wavelet encoding and decoding 40235368Sgnn// 41235368Sgnn// The source code in this file is derived from the encoding 42235368Sgnn// and decoding routines written by Christian Rouet for his 43235368Sgnn// PIZ image file format. 44235368Sgnn// 45235368Sgnn//----------------------------------------------------------------------------- 46235368Sgnn 47235368Sgnn 48235368Sgnn#include <ImfWav.h> 49235368Sgnn 50235368Sgnnnamespace Imf { 51235368Sgnnnamespace { 52235368Sgnn 53235368Sgnn 54235368Sgnn// 55235368Sgnn// Wavelet basis functions without modulo arithmetic; they produce 56235368Sgnn// the best compression ratios when the wavelet-transformed data are 57235368Sgnn// Huffman-encoded, but the wavelet transform works only for 14-bit 58235368Sgnn// data (untransformed data values must be less than (1 << 14)). 59235368Sgnn// 60235368Sgnn 61235368Sgnninline void 62235368Sgnnwenc14 (unsigned short a, unsigned short b, 63235368Sgnn unsigned short &l, unsigned short &h) 64235368Sgnn{ 65235368Sgnn short as = a; 66235368Sgnn short bs = b; 67235368Sgnn 68235368Sgnn short ms = (as + bs) >> 1; 69235368Sgnn short ds = as - bs; 70235368Sgnn 71235368Sgnn l = ms; 72235368Sgnn h = ds; 73235368Sgnn} 74235368Sgnn 75235368Sgnn 76235368Sgnninline void 77235368Sgnnwdec14 (unsigned short l, unsigned short h, 78235368Sgnn unsigned short &a, unsigned short &b) 79235368Sgnn{ 80235368Sgnn short ls = l; 81235368Sgnn short hs = h; 82235368Sgnn 83235368Sgnn int hi = hs; 84235368Sgnn int ai = ls + (hi & 1) + (hi >> 1); 85235368Sgnn 86235368Sgnn short as = ai; 87235368Sgnn short bs = ai - hi; 88235368Sgnn 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 100const int NBITS = 16; 101const int A_OFFSET = 1 << (NBITS - 1); 102const int M_OFFSET = 1 << (NBITS - 1); 103const int MOD_MASK = (1 << NBITS) - 1; 104 105 106inline void 107wenc16 (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 124inline void 125wdec16 (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 143void 144wav2Encode 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 263void 264wav2Decode 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