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