2 * Discrete wavelet transform
3 * Copyright (c) 2007 Kamil Nowosad
4 * Copyright (c) 2013 Nicolas Bertrand <nicoinattendu@gmail.com>
6 * This file is part of FFmpeg.
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.
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.
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
25 * Discrete wavelet transform
28 #include "libavutil/common.h"
29 #include "libavutil/mem.h"
30 #include "jpeg2000dwt.h"
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 #define F_LFTG_K 1.230174104914001f
40 #define F_LFTG_X 1.625732422f
41 /* FIXME: Why use 1.625732422 instead of 1/F_LFTG_K?
42 * Incorrect value in JPEG2000 norm.
43 * see (ISO/IEC 15444:1 (version 2002) F.3.8.2 */
45 /* Lifting parameters in integer format.
46 * Computed as param = (float param) * (1 << 16) */
47 #define I_LFTG_ALPHA 103949
48 #define I_LFTG_BETA 3472
49 #define I_LFTG_GAMMA 57862
50 #define I_LFTG_DELTA 29066
51 #define I_LFTG_K 80621
52 #define I_LFTG_X 106544
54 static inline void extend53(int *p
, int i0
, int i1
)
56 p
[i0
- 1] = p
[i0
+ 1];
58 p
[i0
- 2] = p
[i0
+ 2];
59 p
[i1
+ 1] = p
[i1
- 3];
62 static inline void extend97_float(float *p
, int i0
, int i1
)
66 for (i
= 1; i
<= 4; i
++) {
67 p
[i0
- i
] = p
[i0
+ i
];
68 p
[i1
+ i
- 1] = p
[i1
- i
- 1];
72 static inline void extend97_int(int32_t *p
, int i0
, int i1
)
76 for (i
= 1; i
<= 4; i
++) {
77 p
[i0
- i
] = p
[i0
+ i
];
78 p
[i1
+ i
- 1] = p
[i1
- i
- 1];
82 static void sd_1d53(int *p
, int i0
, int i1
)
91 for (i
= (i0
+1)/2 - 1; i
< (i1
+1)/2; i
++)
92 p
[2*i
+1] -= (p
[2*i
] + p
[2*i
+2]) >> 1;
93 for (i
= (i0
+1)/2; i
< (i1
+1)/2; i
++)
94 p
[2*i
] += (p
[2*i
-1] + p
[2*i
+1] + 2) >> 2;
97 static void dwt_encode53(DWTContext
*s
, int *t
)
100 w
= s
->linelen
[s
->ndeclevels
-1][0];
101 int *line
= s
->i_linebuf
;
104 for (lev
= s
->ndeclevels
-1; lev
>= 0; lev
--){
105 int lh
= s
->linelen
[lev
][0],
106 lv
= s
->linelen
[lev
][1],
114 for (lp
= 0; lp
< lv
; lp
++){
117 for (i
= 0; i
< lh
; i
++)
120 sd_1d53(line
, mh
, mh
+ lh
);
122 // copy back and deinterleave
123 for (i
= mh
; i
< lh
; i
+=2, j
++)
125 for (i
= 1-mh
; i
< lh
; i
+=2, j
++)
131 for (lp
= 0; lp
< lh
; lp
++) {
134 for (i
= 0; i
< lv
; i
++)
137 sd_1d53(line
, mv
, mv
+ lv
);
139 // copy back and deinterleave
140 for (i
= mv
; i
< lv
; i
+=2, j
++)
142 for (i
= 1-mv
; i
< lv
; i
+=2, j
++)
147 static void sd_1d97_float(float *p
, int i0
, int i1
)
154 extend97_float(p
, i0
, i1
);
157 for (i
= i0
/2 - 2; i
< i1
/2 + 1; i
++)
158 p
[2*i
+1] -= 1.586134 * (p
[2*i
] + p
[2*i
+2]);
159 for (i
= i0
/2 - 1; i
< i1
/2 + 1; i
++)
160 p
[2*i
] -= 0.052980 * (p
[2*i
-1] + p
[2*i
+1]);
161 for (i
= i0
/2 - 1; i
< i1
/2; i
++)
162 p
[2*i
+1] += 0.882911 * (p
[2*i
] + p
[2*i
+2]);
163 for (i
= i0
/2; i
< i1
/2; i
++)
164 p
[2*i
] += 0.443506 * (p
[2*i
-1] + p
[2*i
+1]);
167 static void dwt_encode97_float(DWTContext
*s
, float *t
)
170 w
= s
->linelen
[s
->ndeclevels
-1][0];
171 float *line
= s
->f_linebuf
;
174 for (lev
= s
->ndeclevels
-1; lev
>= 0; lev
--){
175 int lh
= s
->linelen
[lev
][0],
176 lv
= s
->linelen
[lev
][1],
184 for (lp
= 0; lp
< lv
; lp
++){
187 for (i
= 0; i
< lh
; i
++)
190 sd_1d97_float(line
, mh
, mh
+ lh
);
192 // copy back and deinterleave
193 for (i
= mh
; i
< lh
; i
+=2, j
++)
194 t
[w
*lp
+ j
] = F_LFTG_X
* l
[i
] / 2;
195 for (i
= 1-mh
; i
< lh
; i
+=2, j
++)
196 t
[w
*lp
+ j
] = F_LFTG_K
* l
[i
] / 2;
201 for (lp
= 0; lp
< lh
; lp
++) {
204 for (i
= 0; i
< lv
; i
++)
207 sd_1d97_float(line
, mv
, mv
+ lv
);
209 // copy back and deinterleave
210 for (i
= mv
; i
< lv
; i
+=2, j
++)
211 t
[w
*j
+ lp
] = F_LFTG_X
* l
[i
] / 2;
212 for (i
= 1-mv
; i
< lv
; i
+=2, j
++)
213 t
[w
*j
+ lp
] = F_LFTG_K
* l
[i
] / 2;
218 static void sd_1d97_int(int *p
, int i0
, int i1
)
225 extend97_int(p
, i0
, i1
);
228 for (i
= i0
/2 - 2; i
< i1
/2 + 1; i
++)
229 p
[2 * i
+ 1] -= (I_LFTG_ALPHA
* (p
[2 * i
] + p
[2 * i
+ 2]) + (1 << 15)) >> 16;
230 for (i
= i0
/2 - 1; i
< i1
/2 + 1; i
++)
231 p
[2 * i
] -= (I_LFTG_BETA
* (p
[2 * i
- 1] + p
[2 * i
+ 1]) + (1 << 15)) >> 16;
232 for (i
= i0
/2 - 1; i
< i1
/2; i
++)
233 p
[2 * i
+ 1] += (I_LFTG_GAMMA
* (p
[2 * i
] + p
[2 * i
+ 2]) + (1 << 15)) >> 16;
234 for (i
= i0
/2; i
< i1
/2; i
++)
235 p
[2 * i
] += (I_LFTG_DELTA
* (p
[2 * i
- 1] + p
[2 * i
+ 1]) + (1 << 15)) >> 16;
238 static void dwt_encode97_int(DWTContext
*s
, int *t
)
241 w
= s
->linelen
[s
->ndeclevels
-1][0];
242 int *line
= s
->i_linebuf
;
245 for (lev
= s
->ndeclevels
-1; lev
>= 0; lev
--){
246 int lh
= s
->linelen
[lev
][0],
247 lv
= s
->linelen
[lev
][1],
255 for (lp
= 0; lp
< lv
; lp
++){
258 for (i
= 0; i
< lh
; i
++)
261 sd_1d97_int(line
, mh
, mh
+ lh
);
263 // copy back and deinterleave
264 for (i
= mh
; i
< lh
; i
+=2, j
++)
265 t
[w
*lp
+ j
] = ((l
[i
] * I_LFTG_X
) + (1 << 16)) >> 17;
266 for (i
= 1-mh
; i
< lh
; i
+=2, j
++)
267 t
[w
*lp
+ j
] = ((l
[i
] * I_LFTG_K
) + (1 << 16)) >> 17;
272 for (lp
= 0; lp
< lh
; lp
++) {
275 for (i
= 0; i
< lv
; i
++)
278 sd_1d97_int(line
, mv
, mv
+ lv
);
280 // copy back and deinterleave
281 for (i
= mv
; i
< lv
; i
+=2, j
++)
282 t
[w
*j
+ lp
] = ((l
[i
] * I_LFTG_X
) + (1 << 16)) >> 17;
283 for (i
= 1-mv
; i
< lv
; i
+=2, j
++)
284 t
[w
*j
+ lp
] = ((l
[i
] * I_LFTG_K
) + (1 << 16)) >> 17;
289 static void sr_1d53(int *p
, int i0
, int i1
)
298 for (i
= i0
/ 2; i
< i1
/ 2 + 1; i
++)
299 p
[2 * i
] -= (p
[2 * i
- 1] + p
[2 * i
+ 1] + 2) >> 2;
300 for (i
= i0
/ 2; i
< i1
/ 2; i
++)
301 p
[2 * i
+ 1] += (p
[2 * i
] + p
[2 * i
+ 2]) >> 1;
304 static void dwt_decode53(DWTContext
*s
, int *t
)
307 int w
= s
->linelen
[s
->ndeclevels
- 1][0];
308 int32_t *line
= s
->i_linebuf
;
311 for (lev
= 0; lev
< s
->ndeclevels
; lev
++) {
312 int lh
= s
->linelen
[lev
][0],
313 lv
= s
->linelen
[lev
][1],
321 for (lp
= 0; lp
< lv
; lp
++) {
323 // copy with interleaving
324 for (i
= mh
; i
< lh
; i
+= 2, j
++)
325 l
[i
] = t
[w
* lp
+ j
];
326 for (i
= 1 - mh
; i
< lh
; i
+= 2, j
++)
327 l
[i
] = t
[w
* lp
+ j
];
329 sr_1d53(line
, mh
, mh
+ lh
);
331 for (i
= 0; i
< lh
; i
++)
332 t
[w
* lp
+ i
] = l
[i
];
337 for (lp
= 0; lp
< lh
; lp
++) {
339 // copy with interleaving
340 for (i
= mv
; i
< lv
; i
+= 2, j
++)
341 l
[i
] = t
[w
* j
+ lp
];
342 for (i
= 1 - mv
; i
< lv
; i
+= 2, j
++)
343 l
[i
] = t
[w
* j
+ lp
];
345 sr_1d53(line
, mv
, mv
+ lv
);
347 for (i
= 0; i
< lv
; i
++)
348 t
[w
* i
+ lp
] = l
[i
];
353 static void sr_1d97_float(float *p
, int i0
, int i1
)
360 extend97_float(p
, i0
, i1
);
362 for (i
= i0
/ 2 - 1; i
< i1
/ 2 + 2; i
++)
363 p
[2 * i
] -= F_LFTG_DELTA
* (p
[2 * i
- 1] + p
[2 * i
+ 1]);
365 for (i
= i0
/ 2 - 1; i
< i1
/ 2 + 1; i
++)
366 p
[2 * i
+ 1] -= F_LFTG_GAMMA
* (p
[2 * i
] + p
[2 * i
+ 2]);
368 for (i
= i0
/ 2; i
< i1
/ 2 + 1; i
++)
369 p
[2 * i
] += F_LFTG_BETA
* (p
[2 * i
- 1] + p
[2 * i
+ 1]);
371 for (i
= i0
/ 2; i
< i1
/ 2; i
++)
372 p
[2 * i
+ 1] += F_LFTG_ALPHA
* (p
[2 * i
] + p
[2 * i
+ 2]);
375 static void dwt_decode97_float(DWTContext
*s
, float *t
)
378 int w
= s
->linelen
[s
->ndeclevels
- 1][0];
379 float *line
= s
->f_linebuf
;
381 /* position at index O of line range [0-5,w+5] cf. extend function */
384 for (lev
= 0; lev
< s
->ndeclevels
; lev
++) {
385 int lh
= s
->linelen
[lev
][0],
386 lv
= s
->linelen
[lev
][1],
393 for (lp
= 0; lp
< lv
; lp
++) {
395 // copy with interleaving
396 for (i
= mh
; i
< lh
; i
+= 2, j
++)
397 l
[i
] = data
[w
* lp
+ j
] * F_LFTG_K
;
398 for (i
= 1 - mh
; i
< lh
; i
+= 2, j
++)
399 l
[i
] = data
[w
* lp
+ j
] * F_LFTG_X
;
401 sr_1d97_float(line
, mh
, mh
+ lh
);
403 for (i
= 0; i
< lh
; i
++)
404 data
[w
* lp
+ i
] = l
[i
];
409 for (lp
= 0; lp
< lh
; lp
++) {
411 // copy with interleaving
412 for (i
= mv
; i
< lv
; i
+= 2, j
++)
413 l
[i
] = data
[w
* j
+ lp
] * F_LFTG_K
;
414 for (i
= 1 - mv
; i
< lv
; i
+= 2, j
++)
415 l
[i
] = data
[w
* j
+ lp
] * F_LFTG_X
;
417 sr_1d97_float(line
, mv
, mv
+ lv
);
419 for (i
= 0; i
< lv
; i
++)
420 data
[w
* i
+ lp
] = l
[i
];
425 static void sr_1d97_int(int32_t *p
, int i0
, int i1
)
432 extend97_int(p
, i0
, i1
);
434 for (i
= i0
/ 2 - 1; i
< i1
/ 2 + 2; i
++)
435 p
[2 * i
] -= (I_LFTG_DELTA
* (p
[2 * i
- 1] + p
[2 * i
+ 1]) + (1 << 15)) >> 16;
437 for (i
= i0
/ 2 - 1; i
< i1
/ 2 + 1; i
++)
438 p
[2 * i
+ 1] -= (I_LFTG_GAMMA
* (p
[2 * i
] + p
[2 * i
+ 2]) + (1 << 15)) >> 16;
440 for (i
= i0
/ 2; i
< i1
/ 2 + 1; i
++)
441 p
[2 * i
] += (I_LFTG_BETA
* (p
[2 * i
- 1] + p
[2 * i
+ 1]) + (1 << 15)) >> 16;
443 for (i
= i0
/ 2; i
< i1
/ 2; i
++)
444 p
[2 * i
+ 1] += (I_LFTG_ALPHA
* (p
[2 * i
] + p
[2 * i
+ 2]) + (1 << 15)) >> 16;
447 static void dwt_decode97_int(DWTContext
*s
, int32_t *t
)
450 int w
= s
->linelen
[s
->ndeclevels
- 1][0];
451 int32_t *line
= s
->i_linebuf
;
453 /* position at index O of line range [0-5,w+5] cf. extend function */
456 for (lev
= 0; lev
< s
->ndeclevels
; lev
++) {
457 int lh
= s
->linelen
[lev
][0],
458 lv
= s
->linelen
[lev
][1],
465 for (lp
= 0; lp
< lv
; lp
++) {
467 // rescale with interleaving
468 for (i
= mh
; i
< lh
; i
+= 2, j
++)
469 l
[i
] = ((data
[w
* lp
+ j
] * I_LFTG_K
) + (1 << 15)) >> 16;
470 for (i
= 1 - mh
; i
< lh
; i
+= 2, j
++)
471 l
[i
] = ((data
[w
* lp
+ j
] * I_LFTG_X
) + (1 << 15)) >> 16;
473 sr_1d97_int(line
, mh
, mh
+ lh
);
475 for (i
= 0; i
< lh
; i
++)
476 data
[w
* lp
+ i
] = l
[i
];
481 for (lp
= 0; lp
< lh
; lp
++) {
483 // rescale with interleaving
484 for (i
= mv
; i
< lv
; i
+= 2, j
++)
485 l
[i
] = ((data
[w
* j
+ lp
] * I_LFTG_K
) + (1 << 15)) >> 16;
486 for (i
= 1 - mv
; i
< lv
; i
+= 2, j
++)
487 l
[i
] = ((data
[w
* j
+ lp
] * I_LFTG_X
) + (1 << 15)) >> 16;
489 sr_1d97_int(line
, mv
, mv
+ lv
);
491 for (i
= 0; i
< lv
; i
++)
492 data
[w
* i
+ lp
] = l
[i
];
497 int ff_jpeg2000_dwt_init(DWTContext
*s
, uint16_t border
[2][2],
498 int decomp_levels
, int type
)
500 int i
, j
, lev
= decomp_levels
, maxlen
,
503 s
->ndeclevels
= decomp_levels
;
506 for (i
= 0; i
< 2; i
++)
507 for (j
= 0; j
< 2; j
++)
508 b
[i
][j
] = border
[i
][j
];
510 maxlen
= FFMAX(b
[0][1] - b
[0][0],
513 for (i
= 0; i
< 2; i
++) {
514 s
->linelen
[lev
][i
] = b
[i
][1] - b
[i
][0];
515 s
->mod
[lev
][i
] = b
[i
][0] & 1;
516 for (j
= 0; j
< 2; j
++)
517 b
[i
][j
] = (b
[i
][j
] + 1) >> 1;
521 s
->f_linebuf
= av_malloc_array((maxlen
+ 12), sizeof(*s
->f_linebuf
));
523 return AVERROR(ENOMEM
);
526 s
->i_linebuf
= av_malloc_array((maxlen
+ 12), sizeof(*s
->i_linebuf
));
528 return AVERROR(ENOMEM
);
531 s
->i_linebuf
= av_malloc_array((maxlen
+ 6), sizeof(*s
->i_linebuf
));
533 return AVERROR(ENOMEM
);
541 int ff_dwt_encode(DWTContext
*s
, void *t
)
545 dwt_encode97_float(s
, t
); break;
547 dwt_encode97_int(s
, t
); break;
549 dwt_encode53(s
, t
); break;
556 int ff_dwt_decode(DWTContext
*s
, void *t
)
560 dwt_decode97_float(s
, t
);
563 dwt_decode97_int(s
, t
);
574 void ff_dwt_destroy(DWTContext
*s
)
576 av_freep(&s
->f_linebuf
);
577 av_freep(&s
->i_linebuf
);