avcodec/jpeg2000dwt: add test for the 5/3 wavelet
[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         return;
89
90     extend53(p, i0, i1);
91
92     for (i = (i0+1)/2 - 1; i < (i1+1)/2; i++)
93         p[2*i+1] -= (p[2*i] + p[2*i+2]) >> 1;
94     for (i = (i0+1)/2; i < (i1+1)/2; i++)
95         p[2*i] += (p[2*i-1] + p[2*i+1] + 2) >> 2;
96 }
97
98 static void dwt_encode53(DWTContext *s, int *t)
99 {
100     int lev,
101         w = s->linelen[s->ndeclevels-1][0];
102     int *line = s->i_linebuf;
103     line += 3;
104
105     for (lev = s->ndeclevels-1; lev >= 0; lev--){
106         int lh = s->linelen[lev][0],
107             lv = s->linelen[lev][1],
108             mh = s->mod[lev][0],
109             mv = s->mod[lev][1],
110             lp;
111         int *l;
112
113         av_assert1(!mh && !mv);
114
115         // VER_SD
116         l = line + mv;
117         for (lp = 0; lp < lh; lp++) {
118             int i, j = 0;
119
120             for (i = 0; i < lv; i++)
121                 l[i] = t[w*i + lp];
122
123             sd_1d53(line, mv, mv + lv);
124
125             // copy back and deinterleave
126             for (i =   mv; i < lv; i+=2, j++)
127                 t[w*j + lp] = l[i];
128             for (i = 1-mv; i < lv; i+=2, j++)
129                 t[w*j + lp] = l[i];
130         }
131
132         // HOR_SD
133         l = line + mh;
134         for (lp = 0; lp < lv; lp++){
135             int i, j = 0;
136
137             for (i = 0; i < lh; i++)
138                 l[i] = t[w*lp + i];
139
140             sd_1d53(line, mh, mh + lh);
141
142             // copy back and deinterleave
143             for (i =   mh; i < lh; i+=2, j++)
144                 t[w*lp + j] = l[i];
145             for (i = 1-mh; i < lh; i+=2, j++)
146                 t[w*lp + j] = l[i];
147         }
148     }
149 }
150 static void sd_1d97_float(float *p, int i0, int i1)
151 {
152     int i;
153
154     if (i1 <= i0 + 1) {
155         if (i0 == 1)
156             p[1] *= F_LFTG_X;
157         else
158             p[0] *= F_LFTG_K;
159         return;
160     }
161
162     extend97_float(p, i0, i1);
163     i0++; i1++;
164
165     for (i = i0/2 - 2; i < i1/2 + 1; i++)
166         p[2*i+1] -= 1.586134 * (p[2*i] + p[2*i+2]);
167     for (i = i0/2 - 1; i < i1/2 + 1; i++)
168         p[2*i] -= 0.052980 * (p[2*i-1] + p[2*i+1]);
169     for (i = i0/2 - 1; i < i1/2; i++)
170         p[2*i+1] += 0.882911 * (p[2*i] + p[2*i+2]);
171     for (i = i0/2; i < i1/2; i++)
172         p[2*i] += 0.443506 * (p[2*i-1] + p[2*i+1]);
173 }
174
175 static void dwt_encode97_float(DWTContext *s, float *t)
176 {
177     int lev,
178         w = s->linelen[s->ndeclevels-1][0];
179     float *line = s->f_linebuf;
180     line += 5;
181
182     for (lev = s->ndeclevels-1; lev >= 0; lev--){
183         int lh = s->linelen[lev][0],
184             lv = s->linelen[lev][1],
185             mh = s->mod[lev][0],
186             mv = s->mod[lev][1],
187             lp;
188         float *l;
189
190         av_assert1(!mh && !mv);
191         // HOR_SD
192         l = line + mh;
193         for (lp = 0; lp < lv; lp++){
194             int i, j = 0;
195
196             for (i = 0; i < lh; i++)
197                 l[i] = t[w*lp + i];
198
199             sd_1d97_float(line, mh, mh + lh);
200
201             // copy back and deinterleave
202             for (i =   mh; i < lh; i+=2, j++)
203                 t[w*lp + j] = F_LFTG_X * l[i] / 2;
204             for (i = 1-mh; i < lh; i+=2, j++)
205                 t[w*lp + j] = F_LFTG_K * l[i] / 2;
206         }
207
208         // VER_SD
209         l = line + mv;
210         for (lp = 0; lp < lh; lp++) {
211             int i, j = 0;
212
213             for (i = 0; i < lv; i++)
214                 l[i] = t[w*i + lp];
215
216             sd_1d97_float(line, mv, mv + lv);
217
218             // copy back and deinterleave
219             for (i =   mv; i < lv; i+=2, j++)
220                 t[w*j + lp] = F_LFTG_X * l[i] / 2;
221             for (i = 1-mv; i < lv; i+=2, j++)
222                 t[w*j + lp] = F_LFTG_K * l[i] / 2;
223         }
224     }
225 }
226
227 static void sd_1d97_int(int *p, int i0, int i1)
228 {
229     int i;
230
231     if (i1 <= i0 + 1) {
232         if (i0 == 1)
233             p[1] = (p[1] * I_LFTG_X + (1<<15)) >> 16;
234         else
235             p[0] = (p[0] * I_LFTG_K + (1<<15)) >> 16;
236         return;
237     }
238
239     extend97_int(p, i0, i1);
240     i0++; i1++;
241
242     for (i = i0/2 - 2; i < i1/2 + 1; i++)
243         p[2 * i + 1] -= (I_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
244     for (i = i0/2 - 1; i < i1/2 + 1; i++)
245         p[2 * i]     -= (I_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
246     for (i = i0/2 - 1; i < i1/2; i++)
247         p[2 * i + 1] += (I_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
248     for (i = i0/2; i < i1/2; i++)
249         p[2 * i]     += (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
250 }
251
252 static void dwt_encode97_int(DWTContext *s, int *t)
253 {
254     int lev,
255         w = s->linelen[s->ndeclevels-1][0];
256     int *line = s->i_linebuf;
257     line += 5;
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         av_assert1(!mh && !mv);
268
269         // VER_SD
270         l = line + mv;
271         for (lp = 0; lp < lh; lp++) {
272             int i, j = 0;
273
274             for (i = 0; i < lv; i++)
275                 l[i] = t[w*i + lp];
276
277             sd_1d97_int(line, mv, mv + lv);
278
279             // copy back and deinterleave
280             for (i =   mv; i < lv; i+=2, j++)
281                 t[w*j + lp] = ((l[i] * I_LFTG_X) + (1 << 16)) >> 17;
282             for (i = 1-mv; i < lv; i+=2, j++)
283                 t[w*j + lp] = ((l[i] * I_LFTG_K) + (1 << 16)) >> 17;
284         }
285
286         // HOR_SD
287         l = line + mh;
288         for (lp = 0; lp < lv; lp++){
289             int i, j = 0;
290
291             for (i = 0; i < lh; i++)
292                 l[i] = t[w*lp + i];
293
294             sd_1d97_int(line, mh, mh + lh);
295
296             // copy back and deinterleave
297             for (i =   mh; i < lh; i+=2, j++)
298                 t[w*lp + j] = ((l[i] * I_LFTG_X) + (1 << 16)) >> 17;
299             for (i = 1-mh; i < lh; i+=2, j++)
300                 t[w*lp + j] = ((l[i] * I_LFTG_K) + (1 << 16)) >> 17;
301         }
302
303     }
304 }
305
306 static void sr_1d53(int *p, int i0, int i1)
307 {
308     int i;
309
310     if (i1 <= i0 + 1) {
311         if (i0 == 1)
312             p[1] >>= 1;
313         return;
314     }
315
316     extend53(p, i0, i1);
317
318     for (i = i0 / 2; i < i1 / 2 + 1; i++)
319         p[2 * i] -= (p[2 * i - 1] + p[2 * i + 1] + 2) >> 2;
320     for (i = i0 / 2; i < i1 / 2; i++)
321         p[2 * i + 1] += (p[2 * i] + p[2 * i + 2]) >> 1;
322 }
323
324 static void dwt_decode53(DWTContext *s, int *t)
325 {
326     int lev;
327     int w     = s->linelen[s->ndeclevels - 1][0];
328     int32_t *line = s->i_linebuf;
329     line += 3;
330
331     for (lev = 0; lev < s->ndeclevels; lev++) {
332         int lh = s->linelen[lev][0],
333             lv = s->linelen[lev][1],
334             mh = s->mod[lev][0],
335             mv = s->mod[lev][1],
336             lp;
337         int *l;
338
339         // HOR_SD
340         l = line + mh;
341         for (lp = 0; lp < lv; lp++) {
342             int i, j = 0;
343             // copy with interleaving
344             for (i = mh; i < lh; i += 2, j++)
345                 l[i] = t[w * lp + j];
346             for (i = 1 - mh; i < lh; i += 2, j++)
347                 l[i] = t[w * lp + j];
348
349             sr_1d53(line, mh, mh + lh);
350
351             for (i = 0; i < lh; i++)
352                 t[w * lp + i] = l[i];
353         }
354
355         // VER_SD
356         l = line + mv;
357         for (lp = 0; lp < lh; lp++) {
358             int i, j = 0;
359             // copy with interleaving
360             for (i = mv; i < lv; i += 2, j++)
361                 l[i] = t[w * j + lp];
362             for (i = 1 - mv; i < lv; i += 2, j++)
363                 l[i] = t[w * j + lp];
364
365             sr_1d53(line, mv, mv + lv);
366
367             for (i = 0; i < lv; i++)
368                 t[w * i + lp] = l[i];
369         }
370     }
371 }
372
373 static void sr_1d97_float(float *p, int i0, int i1)
374 {
375     int i;
376
377     if (i1 <= i0 + 1) {
378         if (i0 == 1)
379             p[1] *= F_LFTG_K/2;
380         else
381             p[0] *= F_LFTG_X/2;
382         return;
383     }
384
385     extend97_float(p, i0, i1);
386
387     for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
388         p[2 * i]     -= F_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]);
389     /* step 4 */
390     for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
391         p[2 * i + 1] -= F_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]);
392     /*step 5*/
393     for (i = i0 / 2; i < i1 / 2 + 1; i++)
394         p[2 * i]     += F_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]);
395     /* step 6 */
396     for (i = i0 / 2; i < i1 / 2; i++)
397         p[2 * i + 1] += F_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]);
398 }
399
400 static void dwt_decode97_float(DWTContext *s, float *t)
401 {
402     int lev;
403     int w       = s->linelen[s->ndeclevels - 1][0];
404     float *line = s->f_linebuf;
405     float *data = t;
406     /* position at index O of line range [0-5,w+5] cf. extend function */
407     line += 5;
408
409     for (lev = 0; lev < s->ndeclevels; lev++) {
410         int lh = s->linelen[lev][0],
411             lv = s->linelen[lev][1],
412             mh = s->mod[lev][0],
413             mv = s->mod[lev][1],
414             lp;
415         float *l;
416         // HOR_SD
417         l = line + mh;
418         for (lp = 0; lp < lv; lp++) {
419             int i, j = 0;
420             // copy with interleaving
421             for (i = mh; i < lh; i += 2, j++)
422                 l[i] = data[w * lp + j] * F_LFTG_K;
423             for (i = 1 - mh; i < lh; i += 2, j++)
424                 l[i] = data[w * lp + j] * F_LFTG_X;
425
426             sr_1d97_float(line, mh, mh + lh);
427
428             for (i = 0; i < lh; i++)
429                 data[w * lp + i] = l[i];
430         }
431
432         // VER_SD
433         l = line + mv;
434         for (lp = 0; lp < lh; lp++) {
435             int i, j = 0;
436             // copy with interleaving
437             for (i = mv; i < lv; i += 2, j++)
438                 l[i] = data[w * j + lp] * F_LFTG_K;
439             for (i = 1 - mv; i < lv; i += 2, j++)
440                 l[i] = data[w * j + lp] * F_LFTG_X;
441
442             sr_1d97_float(line, mv, mv + lv);
443
444             for (i = 0; i < lv; i++)
445                 data[w * i + lp] = l[i];
446         }
447     }
448 }
449
450 static void sr_1d97_int(int32_t *p, int i0, int i1)
451 {
452     int i;
453
454     if (i1 <= i0 + 1) {
455         if (i0 == 1)
456             p[1] = (p[1] * I_LFTG_K + (1<<16)) >> 17;
457         else
458             p[0] = (p[0] * I_LFTG_X + (1<<16)) >> 17;
459         return;
460     }
461
462     extend97_int(p, i0, i1);
463
464     for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
465         p[2 * i]     -= (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
466     /* step 4 */
467     for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
468         p[2 * i + 1] -= (I_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
469     /*step 5*/
470     for (i = i0 / 2; i < i1 / 2 + 1; i++)
471         p[2 * i]     += (I_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
472     /* step 6 */
473     for (i = i0 / 2; i < i1 / 2; i++)
474         p[2 * i + 1] += (I_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
475 }
476
477 static void dwt_decode97_int(DWTContext *s, int32_t *t)
478 {
479     int lev;
480     int w       = s->linelen[s->ndeclevels - 1][0];
481     int32_t *line = s->i_linebuf;
482     int32_t *data = t;
483     /* position at index O of line range [0-5,w+5] cf. extend function */
484     line += 5;
485
486     for (lev = 0; lev < s->ndeclevels; lev++) {
487         int lh = s->linelen[lev][0],
488             lv = s->linelen[lev][1],
489             mh = s->mod[lev][0],
490             mv = s->mod[lev][1],
491             lp;
492         int32_t *l;
493         // HOR_SD
494         l = line + mh;
495         for (lp = 0; lp < lv; lp++) {
496             int i, j = 0;
497             // rescale with interleaving
498             for (i = mh; i < lh; i += 2, j++)
499                 l[i] = ((data[w * lp + j] * I_LFTG_K) + (1 << 15)) >> 16;
500             for (i = 1 - mh; i < lh; i += 2, j++)
501                 l[i] = ((data[w * lp + j] * I_LFTG_X) + (1 << 15)) >> 16;
502
503             sr_1d97_int(line, mh, mh + lh);
504
505             for (i = 0; i < lh; i++)
506                 data[w * lp + i] = l[i];
507         }
508
509         // VER_SD
510         l = line + mv;
511         for (lp = 0; lp < lh; lp++) {
512             int i, j = 0;
513             // rescale with interleaving
514             for (i = mv; i < lv; i += 2, j++)
515                 l[i] = ((data[w * j + lp] * I_LFTG_K) + (1 << 15)) >> 16;
516             for (i = 1 - mv; i < lv; i += 2, j++)
517                 l[i] = ((data[w * j + lp] * I_LFTG_X) + (1 << 15)) >> 16;
518
519             sr_1d97_int(line, mv, mv + lv);
520
521             for (i = 0; i < lv; i++)
522                 data[w * i + lp] = l[i];
523         }
524     }
525 }
526
527 int ff_jpeg2000_dwt_init(DWTContext *s, uint16_t border[2][2],
528                          int decomp_levels, int type)
529 {
530     int i, j, lev = decomp_levels, maxlen,
531         b[2][2];
532
533     s->ndeclevels = decomp_levels;
534     s->type       = type;
535
536     for (i = 0; i < 2; i++)
537         for (j = 0; j < 2; j++)
538             b[i][j] = border[i][j];
539
540     maxlen = FFMAX(b[0][1] - b[0][0],
541                    b[1][1] - b[1][0]);
542     while (--lev >= 0)
543         for (i = 0; i < 2; i++) {
544             s->linelen[lev][i] = b[i][1] - b[i][0];
545             s->mod[lev][i]     = b[i][0] & 1;
546             for (j = 0; j < 2; j++)
547                 b[i][j] = (b[i][j] + 1) >> 1;
548         }
549     switch (type) {
550     case FF_DWT97:
551         s->f_linebuf = av_malloc_array((maxlen + 12), sizeof(*s->f_linebuf));
552         if (!s->f_linebuf)
553             return AVERROR(ENOMEM);
554         break;
555      case FF_DWT97_INT:
556         s->i_linebuf = av_malloc_array((maxlen + 12), sizeof(*s->i_linebuf));
557         if (!s->i_linebuf)
558             return AVERROR(ENOMEM);
559         break;
560     case FF_DWT53:
561         s->i_linebuf = av_malloc_array((maxlen +  6), sizeof(*s->i_linebuf));
562         if (!s->i_linebuf)
563             return AVERROR(ENOMEM);
564         break;
565     default:
566         return -1;
567     }
568     return 0;
569 }
570
571 int ff_dwt_encode(DWTContext *s, void *t)
572 {
573     switch(s->type){
574         case FF_DWT97:
575             dwt_encode97_float(s, t); break;
576         case FF_DWT97_INT:
577             dwt_encode97_int(s, t); break;
578         case FF_DWT53:
579             dwt_encode53(s, t); break;
580         default:
581             return -1;
582     }
583     return 0;
584 }
585
586 int ff_dwt_decode(DWTContext *s, void *t)
587 {
588     switch (s->type) {
589     case FF_DWT97:
590         dwt_decode97_float(s, t);
591         break;
592     case FF_DWT97_INT:
593         dwt_decode97_int(s, t);
594         break;
595     case FF_DWT53:
596         dwt_decode53(s, t);
597         break;
598     default:
599         return -1;
600     }
601     return 0;
602 }
603
604 void ff_dwt_destroy(DWTContext *s)
605 {
606     av_freep(&s->f_linebuf);
607     av_freep(&s->i_linebuf);
608 }
609
610 #ifdef TEST
611
612 #include "libavutil/lfg.h"
613
614 #define MAX_W 256
615
616 int main(void) {
617     int array[MAX_W * MAX_W];
618     int ref  [MAX_W * MAX_W];
619     AVLFG prng;
620     int i,j;
621     DWTContext s1={{{0}}}, *s= &s1;
622     uint16_t border[2][2];
623     int ret, decomp_levels;
624
625     av_lfg_init(&prng, 1);
626
627     for (i = 0; i<MAX_W * MAX_W; i++)
628         array[i] = ref[i] =  av_lfg_get(&prng) % 2048;
629
630     for (i = 0; i < 40; i++) {
631         for (j=0; j<4; j++)
632             border[0][j] = av_lfg_get(&prng) % MAX_W;
633         border[0][0] = border [1][0] =0;
634         if (border[0][0] >= border[0][1] || border[1][0] >= border[1][1])
635             continue;
636         decomp_levels = av_lfg_get(&prng) % FF_DWT_MAX_DECLVLS;
637         ret = ff_jpeg2000_dwt_init(s,  border, decomp_levels, FF_DWT53);
638         if (ret < 0) {
639             fprintf(stderr, "ff_jpeg2000_dwt_init failed\n");
640             return 1;
641         }
642         ret = ff_dwt_encode(s, array);
643         if (ret < 0) {
644             fprintf(stderr, "ff_dwt_encode failed\n");
645             return 1;
646         }
647         ret = ff_dwt_decode(s, array);
648         if (ret < 0) {
649             fprintf(stderr, "ff_dwt_encode failed\n");
650             return 1;
651         }
652         for (j = 0; j<MAX_W * MAX_W; j++)
653             if (array[j] != ref[j]) {
654                 fprintf(stderr, "missmatch at %d (%d != %d) decomp:%d border %d %d %d %d\n",
655                         j, array[j], ref[j],decomp_levels, border[0][0], border[0][1], border[1][0], border[1][1]);
656                 return 2;
657             }
658         ff_dwt_destroy(s);
659     }
660
661     return 0;
662 }
663
664 #endif