jpegdec: be less picky on padding
[ffmpeg.git] / libavcodec / j2k_dwt.c
1 /*
2  * Discrete wavelet transform
3  * Copyright (c) 2007 Kamil Nowosad
4  *
5  * This file is part of FFmpeg.
6  *
7  * FFmpeg is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public
9  * License as published by the Free Software Foundation; either
10  * version 2.1 of the License, or (at your option) any later version.
11  *
12  * FFmpeg is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15  * Lesser General Public License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with FFmpeg; if not, write to the Free Software
19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
20  */
21
22 /**
23  * Discrete wavelet transform
24  * @file
25  * @author Kamil Nowosad
26  */
27
28 #include "j2k_dwt.h"
29
30 static const float scale97[] = {1.625786, 1.230174};
31
32 static inline void extend53(int *p, int i0, int i1)
33 {
34     p[i0 - 1] = p[i0 + 1];
35     p[i1    ] = p[i1 - 2];
36     p[i0 - 2] = p[i0 + 2];
37     p[i1 + 1] = p[i1 - 3];
38 }
39
40 static inline void extend97(float *p, int i0, int i1)
41 {
42     int i;
43
44     for (i = 1; i <= 4; i++){
45         p[i0 - i] = p[i0 + i];
46         p[i1 + i - 1] = p[i1 - i - 1];
47     }
48 }
49
50 static void sd_1d53(int *p, int i0, int i1)
51 {
52     int i;
53
54     if (i1 == i0 + 1)
55         return;
56
57     extend53(p, i0, i1);
58
59     for (i = (i0+1)/2 - 1; i < (i1+1)/2; i++)
60         p[2*i+1] -= (p[2*i] + p[2*i+2]) >> 1;
61     for (i = (i0+1)/2; i < (i1+1)/2; i++)
62         p[2*i] += (p[2*i-1] + p[2*i+1] + 2) >> 2;
63 }
64
65 static void dwt_encode53(DWTContext *s, int *t)
66 {
67     int lev,
68         w = s->linelen[s->ndeclevels-1][0];
69     int *line = s->linebuf;
70     line += 3;
71
72     for (lev = s->ndeclevels-1; lev >= 0; lev--){
73         int lh = s->linelen[lev][0],
74             lv = s->linelen[lev][1],
75             mh = s->mod[lev][0],
76             mv = s->mod[lev][1],
77             lp;
78         int *l;
79
80         // HOR_SD
81         l = line + mh;
82         for (lp = 0; lp < lv; lp++){
83             int i, j = 0;
84
85             for (i = 0; i < lh; i++)
86                 l[i] = t[w*lp + i];
87
88             sd_1d53(line, mh, mh + lh);
89
90             // copy back and deinterleave
91             for (i =   mh; i < lh; i+=2, j++)
92                 t[w*lp + j] = l[i];
93             for (i = 1-mh; i < lh; i+=2, j++)
94                 t[w*lp + j] = l[i];
95         }
96
97         // VER_SD
98         l = line + mv;
99         for (lp = 0; lp < lh; lp++) {
100             int i, j = 0;
101
102             for (i = 0; i < lv; i++)
103                 l[i] = t[w*i + lp];
104
105             sd_1d53(line, mv, mv + lv);
106
107             // copy back and deinterleave
108             for (i =   mv; i < lv; i+=2, j++)
109                 t[w*j + lp] = l[i];
110             for (i = 1-mv; i < lv; i+=2, j++)
111                 t[w*j + lp] = l[i];
112         }
113     }
114 }
115
116 static void sd_1d97(float *p, int i0, int i1)
117 {
118     int i;
119
120     if (i1 == i0 + 1)
121         return;
122
123     extend97(p, i0, i1);
124     i0++; i1++;
125
126     for (i = i0/2 - 2; i < i1/2 + 1; i++)
127         p[2*i+1] -= 1.586134 * (p[2*i] + p[2*i+2]);
128     for (i = i0/2 - 1; i < i1/2 + 1; i++)
129         p[2*i] -= 0.052980 * (p[2*i-1] + p[2*i+1]);
130     for (i = i0/2 - 1; i < i1/2; i++)
131         p[2*i+1] += 0.882911 * (p[2*i] + p[2*i+2]);
132     for (i = i0/2; i < i1/2; i++)
133         p[2*i] += 0.443506 * (p[2*i-1] + p[2*i+1]);
134 }
135
136 static void dwt_encode97(DWTContext *s, int *t)
137 {
138     int lev,
139         w = s->linelen[s->ndeclevels-1][0];
140     float *line = s->linebuf;
141     line += 5;
142
143     for (lev = s->ndeclevels-1; lev >= 0; lev--){
144         int lh = s->linelen[lev][0],
145             lv = s->linelen[lev][1],
146             mh = s->mod[lev][0],
147             mv = s->mod[lev][1],
148             lp;
149         float *l;
150
151         // HOR_SD
152         l = line + mh;
153         for (lp = 0; lp < lv; lp++){
154             int i, j = 0;
155
156             for (i = 0; i < lh; i++)
157                 l[i] = t[w*lp + i];
158
159             sd_1d97(line, mh, mh + lh);
160
161             // copy back and deinterleave
162             for (i =   mh; i < lh; i+=2, j++)
163                 t[w*lp + j] = scale97[mh] * l[i] / 2;
164             for (i = 1-mh; i < lh; i+=2, j++)
165                 t[w*lp + j] = scale97[mh] * l[i] / 2;
166         }
167
168         // VER_SD
169         l = line + mv;
170         for (lp = 0; lp < lh; lp++) {
171             int i, j = 0;
172
173             for (i = 0; i < lv; i++)
174                 l[i] = t[w*i + lp];
175
176             sd_1d97(line, mv, mv + lv);
177
178             // copy back and deinterleave
179             for (i =   mv; i < lv; i+=2, j++)
180                 t[w*j + lp] = scale97[mv] * l[i] / 2;
181             for (i = 1-mv; i < lv; i+=2, j++)
182                 t[w*j + lp] = scale97[mv] * l[i] / 2;
183         }
184     }
185 }
186
187 static void sr_1d53(int *p, int i0, int i1)
188 {
189     int i;
190
191     if (i1 == i0 + 1)
192         return;
193
194     extend53(p, i0, i1);
195
196     for (i = i0/2; i < i1/2 + 1; i++)
197         p[2*i] -= (p[2*i-1] + p[2*i+1] + 2) >> 2;
198     for (i = i0/2; i < i1/2; i++)
199         p[2*i+1] += (p[2*i] + p[2*i+2]) >> 1;
200 }
201
202 static void dwt_decode53(DWTContext *s, int *t)
203 {
204     int lev,
205         w = s->linelen[s->ndeclevels-1][0];
206     int *line = s->linebuf;
207     line += 3;
208
209     for (lev = 0; lev < s->ndeclevels; lev++){
210         int lh = s->linelen[lev][0],
211             lv = s->linelen[lev][1],
212             mh = s->mod[lev][0],
213             mv = s->mod[lev][1],
214             lp;
215         int *l;
216
217         // HOR_SD
218         l = line + mh;
219         for (lp = 0; lp < lv; lp++){
220             int i, j = 0;
221             // copy with interleaving
222             for (i =   mh; i < lh; i+=2, j++)
223                 l[i] = t[w*lp + j];
224             for (i = 1-mh; i < lh; i+=2, j++)
225                 l[i] = t[w*lp + j];
226
227             sr_1d53(line, mh, mh + lh);
228
229             for (i = 0; i < lh; i++)
230                 t[w*lp + i] = l[i];
231         }
232
233         // VER_SD
234         l = line + mv;
235         for (lp = 0; lp < lh; lp++){
236             int i, j = 0;
237             // copy with interleaving
238             for (i =   mv; i < lv; i+=2, j++)
239                 l[i] = t[w*j + lp];
240             for (i = 1-mv; i < lv; i+=2, j++)
241                 l[i] = t[w*j + lp];
242
243             sr_1d53(line, mv, mv + lv);
244
245             for (i = 0; i < lv; i++)
246                 t[w*i + lp] = l[i];
247         }
248     }
249 }
250
251 static void sr_1d97(float *p, int i0, int i1)
252 {
253     int i;
254
255     if (i1 == i0 + 1)
256         return;
257
258     extend97(p, i0, i1);
259
260     for (i = i0/2 - 1; i < i1/2 + 2; i++)
261         p[2*i] -= 0.443506 * (p[2*i-1] + p[2*i+1]);
262     for (i = i0/2 - 1; i < i1/2 + 1; i++)
263         p[2*i+1] -= 0.882911 * (p[2*i] + p[2*i+2]);
264     for (i = i0/2; i < i1/2 + 1; i++)
265         p[2*i] += 0.052980 * (p[2*i-1] + p[2*i+1]);
266     for (i = i0/2; i < i1/2; i++)
267         p[2*i+1] += 1.586134 * (p[2*i] + p[2*i+2]);
268 }
269
270 static void dwt_decode97(DWTContext *s, int *t)
271 {
272     int lev,
273         w = s->linelen[s->ndeclevels-1][0];
274     float *line = s->linebuf;
275     line += 5;
276
277     for (lev = 0; lev < s->ndeclevels; lev++){
278         int lh = s->linelen[lev][0],
279             lv = s->linelen[lev][1],
280             mh = s->mod[lev][0],
281             mv = s->mod[lev][1],
282             lp;
283         float *l;
284
285         // HOR_SD
286         l = line + mh;
287         for (lp = 0; lp < lv; lp++){
288             int i, j = 0;
289             // copy with interleaving
290             for (i =   mh; i < lh; i+=2, j++)
291                 l[i] = scale97[1-mh] * t[w*lp + j];
292             for (i = 1-mh; i < lh; i+=2, j++)
293                 l[i] = scale97[1-mh] * t[w*lp + j];
294
295             sr_1d97(line, mh, mh + lh);
296
297             for (i = 0; i < lh; i++)
298                 t[w*lp + i] = l[i];
299         }
300
301         // VER_SD
302         l = line + mv;
303         for (lp = 0; lp < lh; lp++){
304             int i, j = 0;
305             // copy with interleaving
306             for (i =   mv; i < lv; i+=2, j++)
307                 l[i] = scale97[1-mv] * t[w*j + lp];
308             for (i = 1-mv; i < lv; i+=2, j++)
309                 l[i] = scale97[1-mv] * t[w*j + lp];
310
311             sr_1d97(line, mv, mv + lv);
312
313             for (i = 0; i < lv; i++)
314                 t[w*i + lp] = l[i];
315         }
316     }
317 }
318
319 int ff_j2k_dwt_init(DWTContext *s, uint16_t border[2][2], int decomp_levels, int type)
320 {
321     int i, j, lev = decomp_levels, maxlen,
322         b[2][2];
323
324     if ((unsigned)decomp_levels >= FF_DWT_MAX_DECLVLS)
325         return AVERROR_INVALIDDATA;
326     s->ndeclevels = decomp_levels;
327     s->type = type;
328
329     for (i = 0; i < 2; i++)
330         for(j = 0; j < 2; j++)
331             b[i][j] = border[i][j];
332
333     maxlen = FFMAX(b[0][1] - b[0][0],
334                    b[1][1] - b[1][0]);
335
336     while(--lev >= 0){
337         for (i = 0; i < 2; i++){
338             s->linelen[lev][i] = b[i][1] - b[i][0];
339             s->mod[lev][i] = b[i][0] & 1;
340             for (j = 0; j < 2; j++)
341                 b[i][j] = (b[i][j] + 1) >> 1;
342         }
343     }
344     if (type == FF_DWT97)
345         s->linebuf = av_malloc((maxlen + 12) * sizeof(float));
346     else if (type == FF_DWT53)
347         s->linebuf = av_malloc((maxlen + 6) * sizeof(int));
348     else
349         return -1;
350
351     if (!s->linebuf)
352         return AVERROR(ENOMEM);
353
354     return 0;
355 }
356
357 int ff_j2k_dwt_encode(DWTContext *s, int *t)
358 {
359     switch(s->type){
360         case FF_DWT97:
361             dwt_encode97(s, t); break;
362         case FF_DWT53:
363             dwt_encode53(s, t); break;
364         default:
365             return -1;
366     }
367     return 0;
368 }
369
370 int ff_j2k_dwt_decode(DWTContext *s, int *t)
371 {
372     switch(s->type){
373         case FF_DWT97:
374             dwt_decode97(s, t); break;
375         case FF_DWT53:
376             dwt_decode53(s, t); break;
377         default:
378             return -1;
379     }
380     return 0;
381 }
382
383 void ff_j2k_dwt_destroy(DWTContext *s)
384 {
385     av_freep(&s->linebuf);
386 }