4e7bfa12b46e046cee6feb0c975f4f1ad907534f
[ffmpeg.git] / libavcodec / jpeg2000dwt.c
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 #define F_LFTG_K      1.230174104914001f
41 #define F_LFTG_X      1.625732422f
42 /* FIXME: Why use 1.625732422 instead of 1/F_LFTG_K?
43  * Incorrect value in JPEG2000 norm.
44  * see (ISO/IEC 15444:1 (version 2002) F.3.8.2 */
45
46 /* Lifting parameters in integer format.
47  * Computed as param = (float param) * (1 << 16) */
48 #define I_LFTG_ALPHA  103949
49 #define I_LFTG_BETA     3472
50 #define I_LFTG_GAMMA   57862
51 #define I_LFTG_DELTA   29066
52 #define I_LFTG_K       80621
53 #define I_LFTG_X      106544
54
55 static inline void extend53(int *p, int i0, int i1)
56 {
57     p[i0 - 1] = p[i0 + 1];
58     p[i1]     = p[i1 - 2];
59     p[i0 - 2] = p[i0 + 2];
60     p[i1 + 1] = p[i1 - 3];
61 }
62
63 static inline void extend97_float(float *p, int i0, int i1)
64 {
65     int i;
66
67     for (i = 1; i <= 4; i++) {
68         p[i0 - i]     = p[i0 + i];
69         p[i1 + i - 1] = p[i1 - i - 1];
70     }
71 }
72
73 static inline void extend97_int(int32_t *p, int i0, int i1)
74 {
75     int i;
76
77     for (i = 1; i <= 4; i++) {
78         p[i0 - i]     = p[i0 + i];
79         p[i1 + i - 1] = p[i1 - i - 1];
80     }
81 }
82
83 static void sd_1d53(int *p, int i0, int i1)
84 {
85     int i;
86
87     if (i1 <= i0 + 1) {
88         if (i0 == 1)
89             p[1] <<= 1;
90         return;
91     }
92
93     extend53(p, i0, i1);
94
95     for (i = (i0+1)/2 - 1; i < (i1+1)/2; i++)
96         p[2*i+1] -= (p[2*i] + p[2*i+2]) >> 1;
97     for (i = (i0+1)/2; i < (i1+1)/2; i++)
98         p[2*i] += (p[2*i-1] + p[2*i+1] + 2) >> 2;
99 }
100
101 static void dwt_encode53(DWTContext *s, int *t)
102 {
103     int lev,
104         w = s->linelen[s->ndeclevels-1][0];
105     int *line = s->i_linebuf;
106     line += 3;
107
108     for (lev = s->ndeclevels-1; lev >= 0; lev--){
109         int lh = s->linelen[lev][0],
110             lv = s->linelen[lev][1],
111             mh = s->mod[lev][0],
112             mv = s->mod[lev][1],
113             lp;
114         int *l;
115
116         // VER_SD
117         l = line + mv;
118         for (lp = 0; lp < lh; lp++) {
119             int i, j = 0;
120
121             for (i = 0; i < lv; i++)
122                 l[i] = t[w*i + lp];
123
124             sd_1d53(line, mv, mv + lv);
125
126             // copy back and deinterleave
127             for (i =   mv; i < lv; i+=2, j++)
128                 t[w*j + lp] = l[i];
129             for (i = 1-mv; i < lv; i+=2, j++)
130                 t[w*j + lp] = l[i];
131         }
132
133         // HOR_SD
134         l = line + mh;
135         for (lp = 0; lp < lv; lp++){
136             int i, j = 0;
137
138             for (i = 0; i < lh; i++)
139                 l[i] = t[w*lp + i];
140
141             sd_1d53(line, mh, mh + lh);
142
143             // copy back and deinterleave
144             for (i =   mh; i < lh; i+=2, j++)
145                 t[w*lp + j] = l[i];
146             for (i = 1-mh; i < lh; i+=2, j++)
147                 t[w*lp + j] = l[i];
148         }
149     }
150 }
151 static void sd_1d97_float(float *p, int i0, int i1)
152 {
153     int i;
154
155     if (i1 <= i0 + 1) {
156         if (i0 == 1)
157             p[1] *= F_LFTG_X;
158         else
159             p[0] *= F_LFTG_K;
160         return;
161     }
162
163     extend97_float(p, i0, i1);
164     i0++; i1++;
165
166     for (i = i0/2 - 2; i < i1/2 + 1; i++)
167         p[2*i+1] -= 1.586134 * (p[2*i] + p[2*i+2]);
168     for (i = i0/2 - 1; i < i1/2 + 1; i++)
169         p[2*i] -= 0.052980 * (p[2*i-1] + p[2*i+1]);
170     for (i = i0/2 - 1; i < i1/2; i++)
171         p[2*i+1] += 0.882911 * (p[2*i] + p[2*i+2]);
172     for (i = i0/2; i < i1/2; i++)
173         p[2*i] += 0.443506 * (p[2*i-1] + p[2*i+1]);
174 }
175
176 static void dwt_encode97_float(DWTContext *s, float *t)
177 {
178     int lev,
179         w = s->linelen[s->ndeclevels-1][0];
180     float *line = s->f_linebuf;
181     line += 5;
182
183     for (lev = s->ndeclevels-1; lev >= 0; lev--){
184         int lh = s->linelen[lev][0],
185             lv = s->linelen[lev][1],
186             mh = s->mod[lev][0],
187             mv = s->mod[lev][1],
188             lp;
189         float *l;
190
191         av_assert1(!mh && !mv);
192         // HOR_SD
193         l = line + mh;
194         for (lp = 0; lp < lv; lp++){
195             int i, j = 0;
196
197             for (i = 0; i < lh; i++)
198                 l[i] = t[w*lp + i];
199
200             sd_1d97_float(line, mh, mh + lh);
201
202             // copy back and deinterleave
203             for (i =   mh; i < lh; i+=2, j++)
204                 t[w*lp + j] = F_LFTG_X * l[i] / 2;
205             for (i = 1-mh; i < lh; i+=2, j++)
206                 t[w*lp + j] = F_LFTG_K * l[i] / 2;
207         }
208
209         // VER_SD
210         l = line + mv;
211         for (lp = 0; lp < lh; lp++) {
212             int i, j = 0;
213
214             for (i = 0; i < lv; i++)
215                 l[i] = t[w*i + lp];
216
217             sd_1d97_float(line, mv, mv + lv);
218
219             // copy back and deinterleave
220             for (i =   mv; i < lv; i+=2, j++)
221                 t[w*j + lp] = F_LFTG_X * l[i] / 2;
222             for (i = 1-mv; i < lv; i+=2, j++)
223                 t[w*j + lp] = F_LFTG_K * l[i] / 2;
224         }
225     }
226 }
227
228 static void sd_1d97_int(int *p, int i0, int i1)
229 {
230     int i;
231
232     if (i1 <= i0 + 1) {
233         if (i0 == 1)
234             p[1] = (p[1] * I_LFTG_X + (1<<15)) >> 16;
235         else
236             p[0] = (p[0] * I_LFTG_K + (1<<15)) >> 16;
237         return;
238     }
239
240     extend97_int(p, i0, i1);
241     i0++; i1++;
242
243     for (i = i0/2 - 2; i < i1/2 + 1; i++)
244         p[2 * i + 1] -= (I_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
245     for (i = i0/2 - 1; i < i1/2 + 1; i++)
246         p[2 * i]     -= (I_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
247     for (i = i0/2 - 1; i < i1/2; i++)
248         p[2 * i + 1] += (I_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
249     for (i = i0/2; i < i1/2; i++)
250         p[2 * i]     += (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
251 }
252
253 static void dwt_encode97_int(DWTContext *s, int *t)
254 {
255     int lev,
256         w = s->linelen[s->ndeclevels-1][0];
257     int *line = s->i_linebuf;
258     line += 5;
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         av_assert1(!mh && !mv);
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         // HOR_SD
288         l = line + mh;
289         for (lp = 0; lp < lv; lp++){
290             int i, j = 0;
291
292             for (i = 0; i < lh; i++)
293                 l[i] = t[w*lp + i];
294
295             sd_1d97_int(line, mh, mh + lh);
296
297             // copy back and deinterleave
298             for (i =   mh; i < lh; i+=2, j++)
299                 t[w*lp + j] = ((l[i] * I_LFTG_X) + (1 << 16)) >> 17;
300             for (i = 1-mh; i < lh; i+=2, j++)
301                 t[w*lp + j] = ((l[i] * I_LFTG_K) + (1 << 16)) >> 17;
302         }
303
304     }
305 }
306
307 static void sr_1d53(int *p, int i0, int i1)
308 {
309     int i;
310
311     if (i1 <= i0 + 1) {
312         if (i0 == 1)
313             p[1] >>= 1;
314         return;
315     }
316
317     extend53(p, i0, i1);
318
319     for (i = i0 / 2; i < i1 / 2 + 1; i++)
320         p[2 * i] -= (p[2 * i - 1] + p[2 * i + 1] + 2) >> 2;
321     for (i = i0 / 2; i < i1 / 2; i++)
322         p[2 * i + 1] += (p[2 * i] + p[2 * i + 2]) >> 1;
323 }
324
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
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/2;
383         return;
384     }
385
386     extend97_float(p, i0, i1);
387
388     for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
389         p[2 * i]     -= F_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]);
390     /* step 4 */
391     for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
392         p[2 * i + 1] -= F_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]);
393     /*step 5*/
394     for (i = i0 / 2; i < i1 / 2 + 1; i++)
395         p[2 * i]     += F_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]);
396     /* step 6 */
397     for (i = i0 / 2; i < i1 / 2; i++)
398         p[2 * i + 1] += F_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]);
399 }
400
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] * F_LFTG_K;
424             for (i = 1 - mh; i < lh; i += 2, j++)
425                 l[i] = data[w * lp + j] * F_LFTG_X;
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] * F_LFTG_K;
440             for (i = 1 - mv; i < lv; i += 2, j++)
441                 l[i] = data[w * j + lp] * F_LFTG_X;
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
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<<16)) >> 17;
460         return;
461     }
462
463     extend97_int(p, i0, i1);
464
465     for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
466         p[2 * i]     -= (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
467     /* step 4 */
468     for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
469         p[2 * i + 1] -= (I_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
470     /*step 5*/
471     for (i = i0 / 2; i < i1 / 2 + 1; i++)
472         p[2 * i]     += (I_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
473     /* step 6 */
474     for (i = i0 / 2; i < i1 / 2; i++)
475         p[2 * i + 1] += (I_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
476 }
477
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     int32_t *line = s->i_linebuf;
483     int32_t *data = t;
484     /* position at index O of line range [0-5,w+5] cf. extend function */
485     line += 5;
486
487     for (lev = 0; lev < s->ndeclevels; lev++) {
488         int lh = s->linelen[lev][0],
489             lv = s->linelen[lev][1],
490             mh = s->mod[lev][0],
491             mv = s->mod[lev][1],
492             lp;
493         int32_t *l;
494         // HOR_SD
495         l = line + mh;
496         for (lp = 0; lp < lv; lp++) {
497             int i, j = 0;
498             // rescale with interleaving
499             for (i = mh; i < lh; i += 2, j++)
500                 l[i] = ((data[w * lp + j] * I_LFTG_K) + (1 << 15)) >> 16;
501             for (i = 1 - mh; i < lh; i += 2, j++)
502                 l[i] = ((data[w * lp + j] * I_LFTG_X) + (1 << 15)) >> 16;
503
504             sr_1d97_int(line, mh, mh + lh);
505
506             for (i = 0; i < lh; i++)
507                 data[w * lp + i] = l[i];
508         }
509
510         // VER_SD
511         l = line + mv;
512         for (lp = 0; lp < lh; lp++) {
513             int i, j = 0;
514             // rescale with interleaving
515             for (i = mv; i < lv; i += 2, j++)
516                 l[i] = ((data[w * j + lp] * I_LFTG_K) + (1 << 15)) >> 16;
517             for (i = 1 - mv; i < lv; i += 2, j++)
518                 l[i] = ((data[w * j + lp] * I_LFTG_X) + (1 << 15)) >> 16;
519
520             sr_1d97_int(line, mv, mv + lv);
521
522             for (i = 0; i < lv; i++)
523                 data[w * i + lp] = l[i];
524         }
525     }
526 }
527
528 int ff_jpeg2000_dwt_init(DWTContext *s, uint16_t border[2][2],
529                          int decomp_levels, int type)
530 {
531     int i, j, lev = decomp_levels, maxlen,
532         b[2][2];
533
534     s->ndeclevels = decomp_levels;
535     s->type       = type;
536
537     for (i = 0; i < 2; i++)
538         for (j = 0; j < 2; j++)
539             b[i][j] = border[i][j];
540
541     maxlen = FFMAX(b[0][1] - b[0][0],
542                    b[1][1] - b[1][0]);
543     while (--lev >= 0)
544         for (i = 0; i < 2; i++) {
545             s->linelen[lev][i] = b[i][1] - b[i][0];
546             s->mod[lev][i]     = b[i][0] & 1;
547             for (j = 0; j < 2; j++)
548                 b[i][j] = (b[i][j] + 1) >> 1;
549         }
550     switch (type) {
551     case FF_DWT97:
552         s->f_linebuf = av_malloc_array((maxlen + 12), sizeof(*s->f_linebuf));
553         if (!s->f_linebuf)
554             return AVERROR(ENOMEM);
555         break;
556      case FF_DWT97_INT:
557         s->i_linebuf = av_malloc_array((maxlen + 12), sizeof(*s->i_linebuf));
558         if (!s->i_linebuf)
559             return AVERROR(ENOMEM);
560         break;
561     case FF_DWT53:
562         s->i_linebuf = av_malloc_array((maxlen +  6), sizeof(*s->i_linebuf));
563         if (!s->i_linebuf)
564             return AVERROR(ENOMEM);
565         break;
566     default:
567         return -1;
568     }
569     return 0;
570 }
571
572 int ff_dwt_encode(DWTContext *s, void *t)
573 {
574     switch(s->type){
575         case FF_DWT97:
576             dwt_encode97_float(s, t); break;
577         case FF_DWT97_INT:
578             dwt_encode97_int(s, t); break;
579         case FF_DWT53:
580             dwt_encode53(s, t); break;
581         default:
582             return -1;
583     }
584     return 0;
585 }
586
587 int ff_dwt_decode(DWTContext *s, void *t)
588 {
589     switch (s->type) {
590     case FF_DWT97:
591         dwt_decode97_float(s, t);
592         break;
593     case FF_DWT97_INT:
594         dwt_decode97_int(s, t);
595         break;
596     case FF_DWT53:
597         dwt_decode53(s, t);
598         break;
599     default:
600         return -1;
601     }
602     return 0;
603 }
604
605 void ff_dwt_destroy(DWTContext *s)
606 {
607     av_freep(&s->f_linebuf);
608     av_freep(&s->i_linebuf);
609 }
610
611 #ifdef TEST
612
613 #include "libavutil/lfg.h"
614
615 #define MAX_W 256
616
617 static int test_dwt(int *array, int *ref, uint16_t border[2][2], int decomp_levels, int type) {
618     int ret, j;
619     DWTContext s1={{{0}}}, *s= &s1;
620
621     ret = ff_jpeg2000_dwt_init(s,  border, decomp_levels, type);
622     if (ret < 0) {
623         fprintf(stderr, "ff_jpeg2000_dwt_init failed\n");
624         return 1;
625     }
626     ret = ff_dwt_encode(s, array);
627     if (ret < 0) {
628         fprintf(stderr, "ff_dwt_encode failed\n");
629         return 1;
630     }
631     ret = ff_dwt_decode(s, array);
632     if (ret < 0) {
633         fprintf(stderr, "ff_dwt_encode failed\n");
634         return 1;
635     }
636     for (j = 0; j<MAX_W * MAX_W; j++)
637         if (array[j] != ref[j]) {
638             fprintf(stderr, "missmatch at %d (%d != %d) decomp:%d border %d %d %d %d\n",
639                     j, array[j], ref[j],decomp_levels, border[0][0], border[0][1], border[1][0], border[1][1]);
640             return 2;
641         }
642     ff_dwt_destroy(s);
643
644     return 0;
645 }
646
647 int main(void) {
648     int array[MAX_W * MAX_W];
649     int ref  [MAX_W * MAX_W];
650     AVLFG prng;
651     int i,j;
652     uint16_t border[2][2];
653     int ret, decomp_levels;
654
655     av_lfg_init(&prng, 1);
656
657     for (i = 0; i<MAX_W * MAX_W; i++)
658         array[i] = ref[i] =  av_lfg_get(&prng) % 2048;
659
660     for (i = 0; i < 40; i++) {
661         for (j=0; j<4; j++)
662             border[0][j] = av_lfg_get(&prng) % MAX_W;
663         if (border[0][0] >= border[0][1] || border[1][0] >= border[1][1])
664             continue;
665         decomp_levels = av_lfg_get(&prng) % FF_DWT_MAX_DECLVLS;
666
667         ret = test_dwt(array, ref, border, decomp_levels, FF_DWT53);
668         if (ret)
669             return ret;
670     }
671
672     return 0;
673 }
674
675 #endif