Imported Debian version 2.4.3~trusty1
[deb_ffmpeg.git] / ffmpeg / libavcodec / jpeg2000dwt.c
CommitLineData
2ba45a60
DM
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/common.h"
29#include "libavutil/mem.h"
30#include "jpeg2000dwt.h"
31#include "internal.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#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 */
44
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
53
54static inline void extend53(int *p, int i0, int i1)
55{
56 p[i0 - 1] = p[i0 + 1];
57 p[i1] = p[i1 - 2];
58 p[i0 - 2] = p[i0 + 2];
59 p[i1 + 1] = p[i1 - 3];
60}
61
62static inline void extend97_float(float *p, int i0, int i1)
63{
64 int i;
65
66 for (i = 1; i <= 4; i++) {
67 p[i0 - i] = p[i0 + i];
68 p[i1 + i - 1] = p[i1 - i - 1];
69 }
70}
71
72static inline void extend97_int(int32_t *p, int i0, int i1)
73{
74 int i;
75
76 for (i = 1; i <= 4; i++) {
77 p[i0 - i] = p[i0 + i];
78 p[i1 + i - 1] = p[i1 - i - 1];
79 }
80}
81
82static void sd_1d53(int *p, int i0, int i1)
83{
84 int i;
85
86 if (i1 == i0 + 1)
87 return;
88
89 extend53(p, i0, i1);
90
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;
95}
96
97static 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 // HOR_SD
113 l = line + mh;
114 for (lp = 0; lp < lv; lp++){
115 int i, j = 0;
116
117 for (i = 0; i < lh; i++)
118 l[i] = t[w*lp + i];
119
120 sd_1d53(line, mh, mh + lh);
121
122 // copy back and deinterleave
123 for (i = mh; i < lh; i+=2, j++)
124 t[w*lp + j] = l[i];
125 for (i = 1-mh; i < lh; i+=2, j++)
126 t[w*lp + j] = l[i];
127 }
128
129 // VER_SD
130 l = line + mv;
131 for (lp = 0; lp < lh; lp++) {
132 int i, j = 0;
133
134 for (i = 0; i < lv; i++)
135 l[i] = t[w*i + lp];
136
137 sd_1d53(line, mv, mv + lv);
138
139 // copy back and deinterleave
140 for (i = mv; i < lv; i+=2, j++)
141 t[w*j + lp] = l[i];
142 for (i = 1-mv; i < lv; i+=2, j++)
143 t[w*j + lp] = l[i];
144 }
145 }
146}
147static void sd_1d97_float(float *p, int i0, int i1)
148{
149 int i;
150
151 if (i1 == i0 + 1)
152 return;
153
154 extend97_float(p, i0, i1);
155 i0++; i1++;
156
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]);
165}
166
167static void dwt_encode97_float(DWTContext *s, float *t)
168{
169 int lev,
170 w = s->linelen[s->ndeclevels-1][0];
171 float *line = s->f_linebuf;
172 line += 5;
173
174 for (lev = s->ndeclevels-1; lev >= 0; lev--){
175 int lh = s->linelen[lev][0],
176 lv = s->linelen[lev][1],
177 mh = s->mod[lev][0],
178 mv = s->mod[lev][1],
179 lp;
180 float *l;
181
182 // HOR_SD
183 l = line + mh;
184 for (lp = 0; lp < lv; lp++){
185 int i, j = 0;
186
187 for (i = 0; i < lh; i++)
188 l[i] = t[w*lp + i];
189
190 sd_1d97_float(line, mh, mh + lh);
191
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;
197 }
198
199 // VER_SD
200 l = line + mv;
201 for (lp = 0; lp < lh; lp++) {
202 int i, j = 0;
203
204 for (i = 0; i < lv; i++)
205 l[i] = t[w*i + lp];
206
207 sd_1d97_float(line, mv, mv + lv);
208
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;
214 }
215 }
216}
217
218static void sd_1d97_int(int *p, int i0, int i1)
219{
220 int i;
221
222 if (i1 == i0 + 1)
223 return;
224
225 extend97_int(p, i0, i1);
226 i0++; i1++;
227
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;
236}
237
238static void dwt_encode97_int(DWTContext *s, int *t)
239{
240 int lev,
241 w = s->linelen[s->ndeclevels-1][0];
242 int *line = s->i_linebuf;
243 line += 5;
244
245 for (lev = s->ndeclevels-1; lev >= 0; lev--){
246 int lh = s->linelen[lev][0],
247 lv = s->linelen[lev][1],
248 mh = s->mod[lev][0],
249 mv = s->mod[lev][1],
250 lp;
251 int *l;
252
253 // HOR_SD
254 l = line + mh;
255 for (lp = 0; lp < lv; lp++){
256 int i, j = 0;
257
258 for (i = 0; i < lh; i++)
259 l[i] = t[w*lp + i];
260
261 sd_1d97_int(line, mh, mh + lh);
262
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;
268 }
269
270 // VER_SD
271 l = line + mv;
272 for (lp = 0; lp < lh; lp++) {
273 int i, j = 0;
274
275 for (i = 0; i < lv; i++)
276 l[i] = t[w*i + lp];
277
278 sd_1d97_int(line, mv, mv + lv);
279
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;
285 }
286 }
287}
288
289static void sr_1d53(int *p, int i0, int i1)
290{
291 int i;
292
293 if (i1 == i0 + 1)
294 return;
295
296 extend53(p, i0, i1);
297
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;
302}
303
304static void dwt_decode53(DWTContext *s, int *t)
305{
306 int lev;
307 int w = s->linelen[s->ndeclevels - 1][0];
308 int32_t *line = s->i_linebuf;
309 line += 3;
310
311 for (lev = 0; lev < s->ndeclevels; lev++) {
312 int lh = s->linelen[lev][0],
313 lv = s->linelen[lev][1],
314 mh = s->mod[lev][0],
315 mv = s->mod[lev][1],
316 lp;
317 int *l;
318
319 // HOR_SD
320 l = line + mh;
321 for (lp = 0; lp < lv; lp++) {
322 int i, j = 0;
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];
328
329 sr_1d53(line, mh, mh + lh);
330
331 for (i = 0; i < lh; i++)
332 t[w * lp + i] = l[i];
333 }
334
335 // VER_SD
336 l = line + mv;
337 for (lp = 0; lp < lh; lp++) {
338 int i, j = 0;
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];
344
345 sr_1d53(line, mv, mv + lv);
346
347 for (i = 0; i < lv; i++)
348 t[w * i + lp] = l[i];
349 }
350 }
351}
352
353static void sr_1d97_float(float *p, int i0, int i1)
354{
355 int i;
356
357 if (i1 == i0 + 1)
358 return;
359
360 extend97_float(p, i0, i1);
361
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]);
364 /* step 4 */
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]);
367 /*step 5*/
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]);
370 /* step 6 */
371 for (i = i0 / 2; i < i1 / 2; i++)
372 p[2 * i + 1] += F_LFTG_ALPHA * (p[2 * i] + p[2 * i + 2]);
373}
374
375static void dwt_decode97_float(DWTContext *s, float *t)
376{
377 int lev;
378 int w = s->linelen[s->ndeclevels - 1][0];
379 float *line = s->f_linebuf;
380 float *data = t;
381 /* position at index O of line range [0-5,w+5] cf. extend function */
382 line += 5;
383
384 for (lev = 0; lev < s->ndeclevels; lev++) {
385 int lh = s->linelen[lev][0],
386 lv = s->linelen[lev][1],
387 mh = s->mod[lev][0],
388 mv = s->mod[lev][1],
389 lp;
390 float *l;
391 // HOR_SD
392 l = line + mh;
393 for (lp = 0; lp < lv; lp++) {
394 int i, j = 0;
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;
400
401 sr_1d97_float(line, mh, mh + lh);
402
403 for (i = 0; i < lh; i++)
404 data[w * lp + i] = l[i];
405 }
406
407 // VER_SD
408 l = line + mv;
409 for (lp = 0; lp < lh; lp++) {
410 int i, j = 0;
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;
416
417 sr_1d97_float(line, mv, mv + lv);
418
419 for (i = 0; i < lv; i++)
420 data[w * i + lp] = l[i];
421 }
422 }
423}
424
425static void sr_1d97_int(int32_t *p, int i0, int i1)
426{
427 int i;
428
429 if (i1 == i0 + 1)
430 return;
431
432 extend97_int(p, i0, i1);
433
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;
436 /* step 4 */
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;
439 /*step 5*/
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;
442 /* step 6 */
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;
445}
446
447static void dwt_decode97_int(DWTContext *s, int32_t *t)
448{
449 int lev;
450 int w = s->linelen[s->ndeclevels - 1][0];
451 int32_t *line = s->i_linebuf;
452 int32_t *data = t;
453 /* position at index O of line range [0-5,w+5] cf. extend function */
454 line += 5;
455
456 for (lev = 0; lev < s->ndeclevels; lev++) {
457 int lh = s->linelen[lev][0],
458 lv = s->linelen[lev][1],
459 mh = s->mod[lev][0],
460 mv = s->mod[lev][1],
461 lp;
462 int32_t *l;
463 // HOR_SD
464 l = line + mh;
465 for (lp = 0; lp < lv; lp++) {
466 int i, j = 0;
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;
472
473 sr_1d97_int(line, mh, mh + lh);
474
475 for (i = 0; i < lh; i++)
476 data[w * lp + i] = l[i];
477 }
478
479 // VER_SD
480 l = line + mv;
481 for (lp = 0; lp < lh; lp++) {
482 int i, j = 0;
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;
488
489 sr_1d97_int(line, mv, mv + lv);
490
491 for (i = 0; i < lv; i++)
492 data[w * i + lp] = l[i];
493 }
494 }
495}
496
497int ff_jpeg2000_dwt_init(DWTContext *s, uint16_t border[2][2],
498 int decomp_levels, int type)
499{
500 int i, j, lev = decomp_levels, maxlen,
501 b[2][2];
502
503 s->ndeclevels = decomp_levels;
504 s->type = type;
505
506 for (i = 0; i < 2; i++)
507 for (j = 0; j < 2; j++)
508 b[i][j] = border[i][j];
509
510 maxlen = FFMAX(b[0][1] - b[0][0],
511 b[1][1] - b[1][0]);
512 while (--lev >= 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;
518 }
519 switch (type) {
520 case FF_DWT97:
521 s->f_linebuf = av_malloc_array((maxlen + 12), sizeof(*s->f_linebuf));
522 if (!s->f_linebuf)
523 return AVERROR(ENOMEM);
524 break;
525 case FF_DWT97_INT:
526 s->i_linebuf = av_malloc_array((maxlen + 12), sizeof(*s->i_linebuf));
527 if (!s->i_linebuf)
528 return AVERROR(ENOMEM);
529 break;
530 case FF_DWT53:
531 s->i_linebuf = av_malloc_array((maxlen + 6), sizeof(*s->i_linebuf));
532 if (!s->i_linebuf)
533 return AVERROR(ENOMEM);
534 break;
535 default:
536 return -1;
537 }
538 return 0;
539}
540
541int ff_dwt_encode(DWTContext *s, void *t)
542{
543 switch(s->type){
544 case FF_DWT97:
545 dwt_encode97_float(s, t); break;
546 case FF_DWT97_INT:
547 dwt_encode97_int(s, t); break;
548 case FF_DWT53:
549 dwt_encode53(s, t); break;
550 default:
551 return -1;
552 }
553 return 0;
554}
555
556int ff_dwt_decode(DWTContext *s, void *t)
557{
558 switch (s->type) {
559 case FF_DWT97:
560 dwt_decode97_float(s, t);
561 break;
562 case FF_DWT97_INT:
563 dwt_decode97_int(s, t);
564 break;
565 case FF_DWT53:
566 dwt_decode53(s, t);
567 break;
568 default:
569 return -1;
570 }
571 return 0;
572}
573
574void ff_dwt_destroy(DWTContext *s)
575{
576 av_freep(&s->f_linebuf);
577 av_freep(&s->i_linebuf);
578}