1 /*
2 * Discrete wavelet transform
3 * Copyright (c) 2007 Kamil Nowosad
4 * Copyright (c) 2013 Nicolas Bertrand <nicoinattendu@gmail.com>
5 *
6 * This file is part of FFmpeg.
7 *
8 * FFmpeg is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU Lesser General Public
10 * License as published by the Free Software Foundation; either
11 * version 2.1 of the License, or (at your option) any later version.
12 *
13 * FFmpeg is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 * Lesser General Public License for more details.
17 *
18 * You should have received a copy of the GNU Lesser General Public
19 * License along with FFmpeg; if not, write to the Free Software
20 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21 */
22
23 /**
24 * @file
25 * Discrete wavelet transform
26 */
27
28 #include "libavutil/error.h"
29 #include "libavutil/macros.h"
30 #include "libavutil/mem.h"
31 #include "jpeg2000dwt.h"
32
33 /* Defines for 9/7 DWT lifting parameters.
34 * Parameters are in float. */
35 #define F_LFTG_ALPHA 1.586134342059924f
36 #define F_LFTG_BETA 0.052980118572961f
37 #define F_LFTG_GAMMA 0.882911075530934f
38 #define F_LFTG_DELTA 0.443506852043971f
39
40 /* Lifting parameters in integer format.
41 * Computed as param = (float param) * (1 << 16) */
42 #define I_LFTG_ALPHA 103949ll
43 #define I_LFTG_BETA 3472ll
44 #define I_LFTG_GAMMA 57862ll
45 #define I_LFTG_DELTA 29066ll
46 #define I_LFTG_K 80621ll
47 #define I_LFTG_X 53274ll
48 #define I_PRESHIFT 8
49
extend53(int *p, int i0, int i1)50 static inline void extend53(int *p, int i0, int i1)
51 {
52 p[i0 - 1] = p[i0 + 1];
53 p[i1] = p[i1 - 2];
54 p[i0 - 2] = p[i0 + 2];
55 p[i1 + 1] = p[i1 - 3];
56 }
57
extend97_float(float *p, int i0, int i1)58 static inline void extend97_float(float *p, int i0, int i1)
59 {
60 int i;
61
62 for (i = 1; i <= 4; i++) {
63 p[i0 - i] = p[i0 + i];
64 p[i1 + i - 1] = p[i1 - i - 1];
65 }
66 }
67
extend97_int(int32_t *p, int i0, int i1)68 static inline void extend97_int(int32_t *p, int i0, int i1)
69 {
70 int i;
71
72 for (i = 1; i <= 4; i++) {
73 p[i0 - i] = p[i0 + i];
74 p[i1 + i - 1] = p[i1 - i - 1];
75 }
76 }
77
sd_1d53(int *p, int i0, int i1)78 static void sd_1d53(int *p, int i0, int i1)
79 {
80 int i;
81
82 if (i1 <= i0 + 1) {
83 if (i0 == 1)
84 p[1] <<= 1;
85 return;
86 }
87
88 extend53(p, i0, i1);
89
90 for (i = ((i0+1)>>1) - 1; i < (i1+1)>>1; i++)
91 p[2*i+1] -= (p[2*i] + p[2*i+2]) >> 1;
92 for (i = ((i0+1)>>1); i < (i1+1)>>1; i++)
93 p[2*i] += (p[2*i-1] + p[2*i+1] + 2) >> 2;
94 }
95
dwt_encode53(DWTContext *s, int *t)96 static void dwt_encode53(DWTContext *s, int *t)
97 {
98 int lev,
99 w = s->linelen[s->ndeclevels-1][0];
100 int *line = s->i_linebuf;
101 line += 3;
102
103 for (lev = s->ndeclevels-1; lev >= 0; lev--){
104 int lh = s->linelen[lev][0],
105 lv = s->linelen[lev][1],
106 mh = s->mod[lev][0],
107 mv = s->mod[lev][1],
108 lp;
109 int *l;
110
111 // VER_SD
112 l = line + mv;
113 for (lp = 0; lp < lh; lp++) {
114 int i, j = 0;
115
116 for (i = 0; i < lv; i++)
117 l[i] = t[w*i + lp];
118
119 sd_1d53(line, mv, mv + lv);
120
121 // copy back and deinterleave
122 for (i = mv; i < lv; i+=2, j++)
123 t[w*j + lp] = l[i];
124 for (i = 1-mv; i < lv; i+=2, j++)
125 t[w*j + lp] = l[i];
126 }
127
128 // HOR_SD
129 l = line + mh;
130 for (lp = 0; lp < lv; lp++){
131 int i, j = 0;
132
133 for (i = 0; i < lh; i++)
134 l[i] = t[w*lp + i];
135
136 sd_1d53(line, mh, mh + lh);
137
138 // copy back and deinterleave
139 for (i = mh; i < lh; i+=2, j++)
140 t[w*lp + j] = l[i];
141 for (i = 1-mh; i < lh; i+=2, j++)
142 t[w*lp + j] = l[i];
143 }
144 }
145 }
sd_1d97_float(float *p, int i0, int i1)146 static void sd_1d97_float(float *p, int i0, int i1)
147 {
148 int i;
149
150 if (i1 <= i0 + 1) {
151 if (i0 == 1)
152 p[1] *= F_LFTG_X * 2;
153 else
154 p[0] *= F_LFTG_K;
155 return;
156 }
157
158 extend97_float(p, i0, i1);
159 i0++; i1++;
160
161 for (i = (i0>>1) - 2; i < (i1>>1) + 1; i++)
162 p[2*i+1] -= 1.586134 * (p[2*i] + p[2*i+2]);
163 for (i = (i0>>1) - 1; i < (i1>>1) + 1; i++)
164 p[2*i] -= 0.052980 * (p[2*i-1] + p[2*i+1]);
165 for (i = (i0>>1) - 1; i < (i1>>1); i++)
166 p[2*i+1] += 0.882911 * (p[2*i] + p[2*i+2]);
167 for (i = (i0>>1); i < (i1>>1); i++)
168 p[2*i] += 0.443506 * (p[2*i-1] + p[2*i+1]);
169 }
170
dwt_encode97_float(DWTContext *s, float *t)171 static void dwt_encode97_float(DWTContext *s, float *t)
172 {
173 int lev,
174 w = s->linelen[s->ndeclevels-1][0];
175 float *line = s->f_linebuf;
176 line += 5;
177
178 for (lev = s->ndeclevels-1; lev >= 0; lev--){
179 int lh = s->linelen[lev][0],
180 lv = s->linelen[lev][1],
181 mh = s->mod[lev][0],
182 mv = s->mod[lev][1],
183 lp;
184 float *l;
185
186 // HOR_SD
187 l = line + mh;
188 for (lp = 0; lp < lv; lp++){
189 int i, j = 0;
190
191 for (i = 0; i < lh; i++)
192 l[i] = t[w*lp + i];
193
194 sd_1d97_float(line, mh, mh + lh);
195
196 // copy back and deinterleave
197 for (i = mh; i < lh; i+=2, j++)
198 t[w*lp + j] = l[i];
199 for (i = 1-mh; i < lh; i+=2, j++)
200 t[w*lp + j] = l[i];
201 }
202
203 // VER_SD
204 l = line + mv;
205 for (lp = 0; lp < lh; lp++) {
206 int i, j = 0;
207
208 for (i = 0; i < lv; i++)
209 l[i] = t[w*i + lp];
210
211 sd_1d97_float(line, mv, mv + lv);
212
213 // copy back and deinterleave
214 for (i = mv; i < lv; i+=2, j++)
215 t[w*j + lp] = l[i];
216 for (i = 1-mv; i < lv; i+=2, j++)
217 t[w*j + lp] = l[i];
218 }
219 }
220 }
221
sd_1d97_int(int *p, int i0, int i1)222 static void sd_1d97_int(int *p, int i0, int i1)
223 {
224 int i;
225
226 if (i1 <= i0 + 1) {
227 if (i0 == 1)
228 p[1] = (p[1] * I_LFTG_X + (1<<14)) >> 15;
229 else
230 p[0] = (p[0] * I_LFTG_K + (1<<15)) >> 16;
231 return;
232 }
233
234 extend97_int(p, i0, i1);
235 i0++; i1++;
236
237 for (i = (i0>>1) - 2; i < (i1>>1) + 1; i++)
238 p[2 * i + 1] -= (I_LFTG_ALPHA * (p[2 * i] + p[2 * i + 2]) + (1 << 15)) >> 16;
239 for (i = (i0>>1) - 1; i < (i1>>1) + 1; i++)
240 p[2 * i] -= (I_LFTG_BETA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
241 for (i = (i0>>1) - 1; i < (i1>>1); i++)
242 p[2 * i + 1] += (I_LFTG_GAMMA * (p[2 * i] + p[2 * i + 2]) + (1 << 15)) >> 16;
243 for (i = (i0>>1); i < (i1>>1); i++)
244 p[2 * i] += (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
245 }
246
dwt_encode97_int(DWTContext *s, int *t)247 static void dwt_encode97_int(DWTContext *s, int *t)
248 {
249 int lev;
250 int w = s->linelen[s->ndeclevels-1][0];
251 int h = s->linelen[s->ndeclevels-1][1];
252 int i;
253 int *line = s->i_linebuf;
254 line += 5;
255
256 for (i = 0; i < w * h; i++)
257 t[i] *= 1 << I_PRESHIFT;
258
259 for (lev = s->ndeclevels-1; lev >= 0; lev--){
260 int lh = s->linelen[lev][0],
261 lv = s->linelen[lev][1],
262 mh = s->mod[lev][0],
263 mv = s->mod[lev][1],
264 lp;
265 int *l;
266
267 // VER_SD
268 l = line + mv;
269 for (lp = 0; lp < lh; lp++) {
270 int i, j = 0;
271
272 for (i = 0; i < lv; i++)
273 l[i] = t[w*i + lp];
274
275 sd_1d97_int(line, mv, mv + lv);
276
277 // copy back and deinterleave
278 for (i = mv; i < lv; i+=2, j++)
279 t[w*j + lp] = ((l[i] * I_LFTG_X) + (1 << 15)) >> 16;
280 for (i = 1-mv; i < lv; i+=2, j++)
281 t[w*j + lp] = l[i];
282 }
283
284 // HOR_SD
285 l = line + mh;
286 for (lp = 0; lp < lv; lp++){
287 int i, j = 0;
288
289 for (i = 0; i < lh; i++)
290 l[i] = t[w*lp + i];
291
292 sd_1d97_int(line, mh, mh + lh);
293
294 // copy back and deinterleave
295 for (i = mh; i < lh; i+=2, j++)
296 t[w*lp + j] = ((l[i] * I_LFTG_X) + (1 << 15)) >> 16;
297 for (i = 1-mh; i < lh; i+=2, j++)
298 t[w*lp + j] = l[i];
299 }
300
301 }
302
303 for (i = 0; i < w * h; i++)
304 t[i] = (t[i] + ((1<<I_PRESHIFT)>>1)) >> I_PRESHIFT;
305 }
306
sr_1d53(unsigned *p, int i0, int i1)307 static void sr_1d53(unsigned *p, int i0, int i1)
308 {
309 int i;
310
311 if (i1 <= i0 + 1) {
312 if (i0 == 1)
313 p[1] = (int)p[1] >> 1;
314 return;
315 }
316
317 extend53(p, i0, i1);
318
319 for (i = (i0 >> 1); i < (i1 >> 1) + 1; i++)
320 p[2 * i] -= (int)(p[2 * i - 1] + p[2 * i + 1] + 2) >> 2;
321 for (i = (i0 >> 1); i < (i1 >> 1); i++)
322 p[2 * i + 1] += (int)(p[2 * i] + p[2 * i + 2]) >> 1;
323 }
324
dwt_decode53(DWTContext *s, int *t)325 static void dwt_decode53(DWTContext *s, int *t)
326 {
327 int lev;
328 int w = s->linelen[s->ndeclevels - 1][0];
329 int32_t *line = s->i_linebuf;
330 line += 3;
331
332 for (lev = 0; lev < s->ndeclevels; lev++) {
333 int lh = s->linelen[lev][0],
334 lv = s->linelen[lev][1],
335 mh = s->mod[lev][0],
336 mv = s->mod[lev][1],
337 lp;
338 int *l;
339
340 // HOR_SD
341 l = line + mh;
342 for (lp = 0; lp < lv; lp++) {
343 int i, j = 0;
344 // copy with interleaving
345 for (i = mh; i < lh; i += 2, j++)
346 l[i] = t[w * lp + j];
347 for (i = 1 - mh; i < lh; i += 2, j++)
348 l[i] = t[w * lp + j];
349
350 sr_1d53(line, mh, mh + lh);
351
352 for (i = 0; i < lh; i++)
353 t[w * lp + i] = l[i];
354 }
355
356 // VER_SD
357 l = line + mv;
358 for (lp = 0; lp < lh; lp++) {
359 int i, j = 0;
360 // copy with interleaving
361 for (i = mv; i < lv; i += 2, j++)
362 l[i] = t[w * j + lp];
363 for (i = 1 - mv; i < lv; i += 2, j++)
364 l[i] = t[w * j + lp];
365
366 sr_1d53(line, mv, mv + lv);
367
368 for (i = 0; i < lv; i++)
369 t[w * i + lp] = l[i];
370 }
371 }
372 }
373
sr_1d97_float(float *p, int i0, int i1)374 static void sr_1d97_float(float *p, int i0, int i1)
375 {
376 int i;
377
378 if (i1 <= i0 + 1) {
379 if (i0 == 1)
380 p[1] *= F_LFTG_K/2;
381 else
382 p[0] *= F_LFTG_X;
383 return;
384 }
385
386 extend97_float(p, i0, i1);
387
388 for (i = (i0 >> 1) - 1; i < (i1 >> 1) + 2; i++)
389 p[2 * i] -= F_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]);
390 /* step 4 */
391 for (i = (i0 >> 1) - 1; i < (i1 >> 1) + 1; i++)
392 p[2 * i + 1] -= F_LFTG_GAMMA * (p[2 * i] + p[2 * i + 2]);
393 /*step 5*/
394 for (i = (i0 >> 1); i < (i1 >> 1) + 1; i++)
395 p[2 * i] += F_LFTG_BETA * (p[2 * i - 1] + p[2 * i + 1]);
396 /* step 6 */
397 for (i = (i0 >> 1); i < (i1 >> 1); i++)
398 p[2 * i + 1] += F_LFTG_ALPHA * (p[2 * i] + p[2 * i + 2]);
399 }
400
dwt_decode97_float(DWTContext *s, float *t)401 static void dwt_decode97_float(DWTContext *s, float *t)
402 {
403 int lev;
404 int w = s->linelen[s->ndeclevels - 1][0];
405 float *line = s->f_linebuf;
406 float *data = t;
407 /* position at index O of line range [0-5,w+5] cf. extend function */
408 line += 5;
409
410 for (lev = 0; lev < s->ndeclevels; lev++) {
411 int lh = s->linelen[lev][0],
412 lv = s->linelen[lev][1],
413 mh = s->mod[lev][0],
414 mv = s->mod[lev][1],
415 lp;
416 float *l;
417 // HOR_SD
418 l = line + mh;
419 for (lp = 0; lp < lv; lp++) {
420 int i, j = 0;
421 // copy with interleaving
422 for (i = mh; i < lh; i += 2, j++)
423 l[i] = data[w * lp + j];
424 for (i = 1 - mh; i < lh; i += 2, j++)
425 l[i] = data[w * lp + j];
426
427 sr_1d97_float(line, mh, mh + lh);
428
429 for (i = 0; i < lh; i++)
430 data[w * lp + i] = l[i];
431 }
432
433 // VER_SD
434 l = line + mv;
435 for (lp = 0; lp < lh; lp++) {
436 int i, j = 0;
437 // copy with interleaving
438 for (i = mv; i < lv; i += 2, j++)
439 l[i] = data[w * j + lp];
440 for (i = 1 - mv; i < lv; i += 2, j++)
441 l[i] = data[w * j + lp];
442
443 sr_1d97_float(line, mv, mv + lv);
444
445 for (i = 0; i < lv; i++)
446 data[w * i + lp] = l[i];
447 }
448 }
449 }
450
sr_1d97_int(int32_t *p, int i0, int i1)451 static void sr_1d97_int(int32_t *p, int i0, int i1)
452 {
453 int i;
454
455 if (i1 <= i0 + 1) {
456 if (i0 == 1)
457 p[1] = (p[1] * I_LFTG_K + (1<<16)) >> 17;
458 else
459 p[0] = (p[0] * I_LFTG_X + (1<<15)) >> 16;
460 return;
461 }
462
463 extend97_int(p, i0, i1);
464
465 for (i = (i0 >> 1) - 1; i < (i1 >> 1) + 2; i++)
466 p[2 * i] -= (I_LFTG_DELTA * (p[2 * i - 1] + (int64_t)p[2 * i + 1]) + (1 << 15)) >> 16;
467 /* step 4 */
468 for (i = (i0 >> 1) - 1; i < (i1 >> 1) + 1; i++)
469 p[2 * i + 1] -= (I_LFTG_GAMMA * (p[2 * i] + (int64_t)p[2 * i + 2]) + (1 << 15)) >> 16;
470 /*step 5*/
471 for (i = (i0 >> 1); i < (i1 >> 1) + 1; i++)
472 p[2 * i] += (I_LFTG_BETA * (p[2 * i - 1] + (int64_t)p[2 * i + 1]) + (1 << 15)) >> 16;
473 /* step 6 */
474 for (i = (i0 >> 1); i < (i1 >> 1); i++)
475 p[2 * i + 1] += (I_LFTG_ALPHA * (p[2 * i] + (int64_t)p[2 * i + 2]) + (1 << 15)) >> 16;
476 }
477
dwt_decode97_int(DWTContext *s, int32_t *t)478 static void dwt_decode97_int(DWTContext *s, int32_t *t)
479 {
480 int lev;
481 int w = s->linelen[s->ndeclevels - 1][0];
482 int h = s->linelen[s->ndeclevels - 1][1];
483 int i;
484 int32_t *line = s->i_linebuf;
485 int32_t *data = t;
486 /* position at index O of line range [0-5,w+5] cf. extend function */
487 line += 5;
488
489 for (i = 0; i < w * h; i++)
490 data[i] *= 1LL << I_PRESHIFT;
491
492 for (lev = 0; lev < s->ndeclevels; lev++) {
493 int lh = s->linelen[lev][0],
494 lv = s->linelen[lev][1],
495 mh = s->mod[lev][0],
496 mv = s->mod[lev][1],
497 lp;
498 int32_t *l;
499 // HOR_SD
500 l = line + mh;
501 for (lp = 0; lp < lv; lp++) {
502 int i, j = 0;
503 // rescale with interleaving
504 for (i = mh; i < lh; i += 2, j++)
505 l[i] = ((data[w * lp + j] * I_LFTG_K) + (1 << 15)) >> 16;
506 for (i = 1 - mh; i < lh; i += 2, j++)
507 l[i] = data[w * lp + j];
508
509 sr_1d97_int(line, mh, mh + lh);
510
511 for (i = 0; i < lh; i++)
512 data[w * lp + i] = l[i];
513 }
514
515 // VER_SD
516 l = line + mv;
517 for (lp = 0; lp < lh; lp++) {
518 int i, j = 0;
519 // rescale with interleaving
520 for (i = mv; i < lv; i += 2, j++)
521 l[i] = ((data[w * j + lp] * I_LFTG_K) + (1 << 15)) >> 16;
522 for (i = 1 - mv; i < lv; i += 2, j++)
523 l[i] = data[w * j + lp];
524
525 sr_1d97_int(line, mv, mv + lv);
526
527 for (i = 0; i < lv; i++)
528 data[w * i + lp] = l[i];
529 }
530 }
531
532 for (i = 0; i < w * h; i++)
533 data[i] = (data[i] + ((1LL<<I_PRESHIFT)>>1)) >> I_PRESHIFT;
534 }
535
ff_jpeg2000_dwt_init(DWTContext *s, int border[2][2], int decomp_levels, int type)536 int ff_jpeg2000_dwt_init(DWTContext *s, int border[2][2],
537 int decomp_levels, int type)
538 {
539 int i, j, lev = decomp_levels, maxlen,
540 b[2][2];
541
542 s->ndeclevels = decomp_levels;
543 s->type = type;
544
545 for (i = 0; i < 2; i++)
546 for (j = 0; j < 2; j++)
547 b[i][j] = border[i][j];
548
549 maxlen = FFMAX(b[0][1] - b[0][0],
550 b[1][1] - b[1][0]);
551 while (--lev >= 0)
552 for (i = 0; i < 2; i++) {
553 s->linelen[lev][i] = b[i][1] - b[i][0];
554 s->mod[lev][i] = b[i][0] & 1;
555 for (j = 0; j < 2; j++)
556 b[i][j] = (b[i][j] + 1) >> 1;
557 }
558 switch (type) {
559 case FF_DWT97:
560 s->f_linebuf = av_malloc_array((maxlen + 12), sizeof(*s->f_linebuf));
561 if (!s->f_linebuf)
562 return AVERROR(ENOMEM);
563 break;
564 case FF_DWT97_INT:
565 s->i_linebuf = av_malloc_array((maxlen + 12), sizeof(*s->i_linebuf));
566 if (!s->i_linebuf)
567 return AVERROR(ENOMEM);
568 break;
569 case FF_DWT53:
570 s->i_linebuf = av_malloc_array((maxlen + 6), sizeof(*s->i_linebuf));
571 if (!s->i_linebuf)
572 return AVERROR(ENOMEM);
573 break;
574 default:
575 return -1;
576 }
577 return 0;
578 }
579
ff_dwt_encode(DWTContext *s, void *t)580 int ff_dwt_encode(DWTContext *s, void *t)
581 {
582 if (s->ndeclevels == 0)
583 return 0;
584
585 switch(s->type){
586 case FF_DWT97:
587 dwt_encode97_float(s, t); break;
588 case FF_DWT97_INT:
589 dwt_encode97_int(s, t); break;
590 case FF_DWT53:
591 dwt_encode53(s, t); break;
592 default:
593 return -1;
594 }
595 return 0;
596 }
597
ff_dwt_decode(DWTContext *s, void *t)598 int ff_dwt_decode(DWTContext *s, void *t)
599 {
600 if (s->ndeclevels == 0)
601 return 0;
602
603 switch (s->type) {
604 case FF_DWT97:
605 dwt_decode97_float(s, t);
606 break;
607 case FF_DWT97_INT:
608 dwt_decode97_int(s, t);
609 break;
610 case FF_DWT53:
611 dwt_decode53(s, t);
612 break;
613 default:
614 return -1;
615 }
616 return 0;
617 }
618
ff_dwt_destroy(DWTContext *s)619 void ff_dwt_destroy(DWTContext *s)
620 {
621 av_freep(&s->f_linebuf);
622 av_freep(&s->i_linebuf);
623 }
624