FFmpeg
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
j2k_dwt.c
Go to the documentation of this file.
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 
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 
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 
384 {
385  av_freep(&s->linebuf);
386 }