You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

521 lines
14KB

  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. * @file
  24. * Discrete wavelet transform
  25. */
  26. #include "libavutil/common.h"
  27. #include "libavutil/mem.h"
  28. #include "jpeg2000dwt.h"
  29. #include "internal.h"
  30. /* Defines for 9/7 DWT lifting parameters.
  31. * Parameters are in float. */
  32. #define F_LFTG_ALPHA 1.586134342059924f
  33. #define F_LFTG_BETA 0.052980118572961f
  34. #define F_LFTG_GAMMA 0.882911075530934f
  35. #define F_LFTG_DELTA 0.443506852043971f
  36. #define F_LFTG_K 1.230174104914001f
  37. #define F_LFTG_X 1.625732422f
  38. /* FIXME: Why use 1.625732422 instead of 1/F_LFTG_K?
  39. * Incorrect value in JPEG2000 norm.
  40. * see (ISO/IEC 15444:1 (version 2002) F.3.8.2 */
  41. /* Lifting parameters in integer format.
  42. * Computed as param = (float param) * (1 << 16) */
  43. #define I_LFTG_ALPHA 103949
  44. #define I_LFTG_BETA 3472
  45. #define I_LFTG_GAMMA 57862
  46. #define I_LFTG_DELTA 29066
  47. #define I_LFTG_K 80621
  48. #define I_LFTG_X 106544
  49. static inline void extend53(int *p, int i0, int i1)
  50. {
  51. p[i0 - 1] = p[i0 + 1];
  52. p[i1] = p[i1 - 2];
  53. p[i0 - 2] = p[i0 + 2];
  54. p[i1 + 1] = p[i1 - 3];
  55. }
  56. static inline void extend97_float(float *p, int i0, int i1)
  57. {
  58. int i;
  59. for (i = 1; i <= 4; i++) {
  60. p[i0 - i] = p[i0 + i];
  61. p[i1 + i - 1] = p[i1 - i - 1];
  62. }
  63. }
  64. static inline void extend97_int(int32_t *p, int i0, int i1)
  65. {
  66. int i;
  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. static void sd_1d53(int *p, int i0, int i1)
  73. {
  74. int i;
  75. if (i1 == i0 + 1)
  76. return;
  77. extend53(p, i0, i1);
  78. for (i = (i0+1)/2 - 1; i < (i1+1)/2; i++)
  79. p[2*i+1] -= (p[2*i] + p[2*i+2]) >> 1;
  80. for (i = (i0+1)/2; i < (i1+1)/2; i++)
  81. p[2*i] += (p[2*i-1] + p[2*i+1] + 2) >> 2;
  82. }
  83. static void dwt_encode53(DWTContext *s, int *t)
  84. {
  85. int lev,
  86. w = s->linelen[s->ndeclevels-1][0];
  87. int *line = s->i_linebuf;
  88. line += 3;
  89. for (lev = s->ndeclevels-1; lev >= 0; lev--){
  90. int lh = s->linelen[lev][0],
  91. lv = s->linelen[lev][1],
  92. mh = s->mod[lev][0],
  93. mv = s->mod[lev][1],
  94. lp;
  95. int *l;
  96. // HOR_SD
  97. l = line + mh;
  98. for (lp = 0; lp < lv; lp++){
  99. int i, j = 0;
  100. for (i = 0; i < lh; i++)
  101. l[i] = t[w*lp + i];
  102. sd_1d53(line, mh, mh + lh);
  103. // copy back and deinterleave
  104. for (i = mh; i < lh; i+=2, j++)
  105. t[w*lp + j] = l[i];
  106. for (i = 1-mh; i < lh; i+=2, j++)
  107. t[w*lp + j] = l[i];
  108. }
  109. // VER_SD
  110. l = line + mv;
  111. for (lp = 0; lp < lh; lp++) {
  112. int i, j = 0;
  113. for (i = 0; i < lv; i++)
  114. l[i] = t[w*i + lp];
  115. sd_1d53(line, mv, mv + lv);
  116. // copy back and deinterleave
  117. for (i = mv; i < lv; i+=2, j++)
  118. t[w*j + lp] = l[i];
  119. for (i = 1-mv; i < lv; i+=2, j++)
  120. t[w*j + lp] = l[i];
  121. }
  122. }
  123. }
  124. static void sd_1d97_int(int *p, int i0, int i1)
  125. {
  126. int i;
  127. if (i1 == i0 + 1)
  128. return;
  129. extend97_int(p, i0, i1);
  130. i0++; i1++;
  131. for (i = i0/2 - 2; i < i1/2 + 1; i++)
  132. p[2 * i + 1] -= (I_LFTG_ALPHA * (p[2 * i] + p[2 * i + 2]) + (1 << 15)) >> 16;
  133. for (i = i0/2 - 1; i < i1/2 + 1; i++)
  134. p[2 * i] -= (I_LFTG_BETA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
  135. for (i = i0/2 - 1; i < i1/2; i++)
  136. p[2 * i + 1] += (I_LFTG_GAMMA * (p[2 * i] + p[2 * i + 2]) + (1 << 15)) >> 16;
  137. for (i = i0/2; i < i1/2; i++)
  138. p[2 * i] += (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
  139. }
  140. static void dwt_encode97_int(DWTContext *s, int *t)
  141. {
  142. int lev,
  143. w = s->linelen[s->ndeclevels-1][0];
  144. int *line = s->i_linebuf;
  145. line += 5;
  146. for (lev = s->ndeclevels-1; lev >= 0; lev--){
  147. int lh = s->linelen[lev][0],
  148. lv = s->linelen[lev][1],
  149. mh = s->mod[lev][0],
  150. mv = s->mod[lev][1],
  151. lp;
  152. int *l;
  153. // HOR_SD
  154. l = line + mh;
  155. for (lp = 0; lp < lv; lp++){
  156. int i, j = 0;
  157. for (i = 0; i < lh; i++)
  158. l[i] = t[w*lp + i];
  159. sd_1d97_int(line, mh, mh + lh);
  160. // copy back and deinterleave
  161. for (i = mh; i < lh; i+=2, j++)
  162. t[w*lp + j] = ((l[i] * I_LFTG_X) + (1 << 16)) >> 17;
  163. for (i = 1-mh; i < lh; i+=2, j++)
  164. t[w*lp + j] = ((l[i] * I_LFTG_K) + (1 << 16)) >> 17;
  165. }
  166. // VER_SD
  167. l = line + mv;
  168. for (lp = 0; lp < lh; lp++) {
  169. int i, j = 0;
  170. for (i = 0; i < lv; i++)
  171. l[i] = t[w*i + lp];
  172. sd_1d97_int(line, mv, mv + lv);
  173. // copy back and deinterleave
  174. for (i = mv; i < lv; i+=2, j++)
  175. t[w*j + lp] = ((l[i] * I_LFTG_X) + (1 << 16)) >> 17;
  176. for (i = 1-mv; i < lv; i+=2, j++)
  177. t[w*j + lp] = ((l[i] * I_LFTG_K) + (1 << 16)) >> 17;
  178. }
  179. }
  180. }
  181. static void sr_1d53(int *p, int i0, int i1)
  182. {
  183. int i;
  184. if (i1 == i0 + 1)
  185. return;
  186. extend53(p, i0, i1);
  187. for (i = i0 / 2; i < i1 / 2 + 1; i++)
  188. p[2 * i] -= (p[2 * i - 1] + p[2 * i + 1] + 2) >> 2;
  189. for (i = i0 / 2; i < i1 / 2; i++)
  190. p[2 * i + 1] += (p[2 * i] + p[2 * i + 2]) >> 1;
  191. }
  192. static void dwt_decode53(DWTContext *s, int *t)
  193. {
  194. int lev;
  195. int w = s->linelen[s->ndeclevels - 1][0];
  196. int32_t *line = s->i_linebuf;
  197. line += 3;
  198. for (lev = 0; lev < s->ndeclevels; lev++) {
  199. int lh = s->linelen[lev][0],
  200. lv = s->linelen[lev][1],
  201. mh = s->mod[lev][0],
  202. mv = s->mod[lev][1],
  203. lp;
  204. int *l;
  205. // HOR_SD
  206. l = line + mh;
  207. for (lp = 0; lp < lv; lp++) {
  208. int i, j = 0;
  209. // copy with interleaving
  210. for (i = mh; i < lh; i += 2, j++)
  211. l[i] = t[w * lp + j];
  212. for (i = 1 - mh; i < lh; i += 2, j++)
  213. l[i] = t[w * lp + j];
  214. sr_1d53(line, mh, mh + lh);
  215. for (i = 0; i < lh; i++)
  216. t[w * lp + i] = l[i];
  217. }
  218. // VER_SD
  219. l = line + mv;
  220. for (lp = 0; lp < lh; lp++) {
  221. int i, j = 0;
  222. // copy with interleaving
  223. for (i = mv; i < lv; i += 2, j++)
  224. l[i] = t[w * j + lp];
  225. for (i = 1 - mv; i < lv; i += 2, j++)
  226. l[i] = t[w * j + lp];
  227. sr_1d53(line, mv, mv + lv);
  228. for (i = 0; i < lv; i++)
  229. t[w * i + lp] = l[i];
  230. }
  231. }
  232. }
  233. static void sr_1d97_float(float *p, int i0, int i1)
  234. {
  235. int i;
  236. if (i1 == i0 + 1)
  237. return;
  238. extend97_float(p, i0, i1);
  239. /*step 1*/
  240. for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
  241. p[2 * i] *= F_LFTG_K;
  242. /* step 2*/
  243. for (i = i0 / 2 - 2; i < i1 / 2 + 2; i++)
  244. p[2 * i + 1] *= F_LFTG_X;
  245. /* step 3*/
  246. for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
  247. p[2 * i] -= F_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]);
  248. /* step 4 */
  249. for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
  250. p[2 * i + 1] -= F_LFTG_GAMMA * (p[2 * i] + p[2 * i + 2]);
  251. /*step 5*/
  252. for (i = i0 / 2; i < i1 / 2 + 1; i++)
  253. p[2 * i] += F_LFTG_BETA * (p[2 * i - 1] + p[2 * i + 1]);
  254. /* step 6 */
  255. for (i = i0 / 2; i < i1 / 2; i++)
  256. p[2 * i + 1] += F_LFTG_ALPHA * (p[2 * i] + p[2 * i + 2]);
  257. }
  258. static void dwt_decode97_float(DWTContext *s, float *t)
  259. {
  260. int lev;
  261. int w = s->linelen[s->ndeclevels - 1][0];
  262. float *line = s->f_linebuf;
  263. float *data = t;
  264. /* position at index O of line range [0-5,w+5] cf. extend function */
  265. line += 5;
  266. for (lev = 0; lev < s->ndeclevels; lev++) {
  267. int lh = s->linelen[lev][0],
  268. lv = s->linelen[lev][1],
  269. mh = s->mod[lev][0],
  270. mv = s->mod[lev][1],
  271. lp;
  272. float *l;
  273. // HOR_SD
  274. l = line + mh;
  275. for (lp = 0; lp < lv; lp++) {
  276. int i, j = 0;
  277. // copy with interleaving
  278. for (i = mh; i < lh; i += 2, j++)
  279. l[i] = data[w * lp + j];
  280. for (i = 1 - mh; i < lh; i += 2, j++)
  281. l[i] = data[w * lp + j];
  282. sr_1d97_float(line, mh, mh + lh);
  283. for (i = 0; i < lh; i++)
  284. data[w * lp + i] = l[i];
  285. }
  286. // VER_SD
  287. l = line + mv;
  288. for (lp = 0; lp < lh; lp++) {
  289. int i, j = 0;
  290. // copy with interleaving
  291. for (i = mv; i < lv; i += 2, j++)
  292. l[i] = data[w * j + lp];
  293. for (i = 1 - mv; i < lv; i += 2, j++)
  294. l[i] = data[w * j + lp];
  295. sr_1d97_float(line, mv, mv + lv);
  296. for (i = 0; i < lv; i++)
  297. data[w * i + lp] = l[i];
  298. }
  299. }
  300. }
  301. static void sr_1d97_int(int32_t *p, int i0, int i1)
  302. {
  303. int i;
  304. if (i1 == i0 + 1)
  305. return;
  306. extend97_int(p, i0, i1);
  307. /*step 1*/
  308. for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
  309. p[2 * i] = ((p[2 * i] * I_LFTG_K) + (1 << 15)) >> 16;
  310. /* step 2*/
  311. for (i = i0 / 2 - 2; i < i1 / 2 + 2; i++)
  312. p[2 * i + 1] = ((p[2 * i + 1] * I_LFTG_X) + (1 << 15)) >> 16;
  313. /* step 3*/
  314. for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
  315. p[2 * i] -= (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
  316. /* step 4 */
  317. for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
  318. p[2 * i + 1] -= (I_LFTG_GAMMA * (p[2 * i] + p[2 * i + 2]) + (1 << 15)) >> 16;
  319. /*step 5*/
  320. for (i = i0 / 2; i < i1 / 2 + 1; i++)
  321. p[2 * i] += (I_LFTG_BETA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
  322. /* step 6 */
  323. for (i = i0 / 2; i < i1 / 2; i++)
  324. p[2 * i + 1] += (I_LFTG_ALPHA * (p[2 * i] + p[2 * i + 2]) + (1 << 15)) >> 16;
  325. }
  326. static void dwt_decode97_int(DWTContext *s, int32_t *t)
  327. {
  328. int lev;
  329. int w = s->linelen[s->ndeclevels - 1][0];
  330. int32_t *line = s->i_linebuf;
  331. int32_t *data = t;
  332. /* position at index O of line range [0-5,w+5] cf. extend function */
  333. line += 5;
  334. for (lev = 0; lev < s->ndeclevels; lev++) {
  335. int lh = s->linelen[lev][0],
  336. lv = s->linelen[lev][1],
  337. mh = s->mod[lev][0],
  338. mv = s->mod[lev][1],
  339. lp;
  340. int32_t *l;
  341. // HOR_SD
  342. l = line + mh;
  343. for (lp = 0; lp < lv; lp++) {
  344. int i, j = 0;
  345. // copy with interleaving
  346. for (i = mh; i < lh; i += 2, j++)
  347. l[i] = data[w * lp + j];
  348. for (i = 1 - mh; i < lh; i += 2, j++)
  349. l[i] = data[w * lp + j];
  350. sr_1d97_int(line, mh, mh + lh);
  351. for (i = 0; i < lh; i++)
  352. data[w * lp + i] = l[i];
  353. }
  354. // VER_SD
  355. l = line + mv;
  356. for (lp = 0; lp < lh; lp++) {
  357. int i, j = 0;
  358. // copy with interleaving
  359. for (i = mv; i < lv; i += 2, j++)
  360. l[i] = data[w * j + lp];
  361. for (i = 1 - mv; i < lv; i += 2, j++)
  362. l[i] = data[w * j + lp];
  363. sr_1d97_int(line, mv, mv + lv);
  364. for (i = 0; i < lv; i++)
  365. data[w * i + lp] = l[i];
  366. }
  367. }
  368. }
  369. int ff_jpeg2000_dwt_init(DWTContext *s, uint16_t border[2][2],
  370. int decomp_levels, int type)
  371. {
  372. int i, j, lev = decomp_levels, maxlen,
  373. b[2][2];
  374. s->ndeclevels = decomp_levels;
  375. s->type = type;
  376. for (i = 0; i < 2; i++)
  377. for (j = 0; j < 2; j++)
  378. b[i][j] = border[i][j];
  379. maxlen = FFMAX(b[0][1] - b[0][0],
  380. b[1][1] - b[1][0]);
  381. while (--lev >= 0)
  382. for (i = 0; i < 2; i++) {
  383. s->linelen[lev][i] = b[i][1] - b[i][0];
  384. s->mod[lev][i] = b[i][0] & 1;
  385. for (j = 0; j < 2; j++)
  386. b[i][j] = (b[i][j] + 1) >> 1;
  387. }
  388. switch (type) {
  389. case FF_DWT97:
  390. s->f_linebuf = av_malloc((maxlen + 12) * sizeof(*s->f_linebuf));
  391. if (!s->f_linebuf)
  392. return AVERROR(ENOMEM);
  393. break;
  394. case FF_DWT97_INT:
  395. s->i_linebuf = av_malloc((maxlen + 12) * sizeof(*s->i_linebuf));
  396. if (!s->i_linebuf)
  397. return AVERROR(ENOMEM);
  398. break;
  399. case FF_DWT53:
  400. s->i_linebuf = av_malloc((maxlen + 6) * sizeof(*s->i_linebuf));
  401. if (!s->i_linebuf)
  402. return AVERROR(ENOMEM);
  403. break;
  404. default:
  405. return -1;
  406. }
  407. return 0;
  408. }
  409. int ff_dwt_encode(DWTContext *s, int *t)
  410. {
  411. switch(s->type){
  412. case FF_DWT97_INT:
  413. dwt_encode97_int(s, t); break;
  414. case FF_DWT53:
  415. dwt_encode53(s, t); break;
  416. default:
  417. return -1;
  418. }
  419. return 0;
  420. }
  421. int ff_dwt_decode(DWTContext *s, void *t)
  422. {
  423. switch (s->type) {
  424. case FF_DWT97:
  425. dwt_decode97_float(s, t);
  426. break;
  427. case FF_DWT97_INT:
  428. dwt_decode97_int(s, t);
  429. break;
  430. case FF_DWT53:
  431. dwt_decode53(s, t);
  432. break;
  433. default:
  434. return -1;
  435. }
  436. return 0;
  437. }
  438. void ff_dwt_destroy(DWTContext *s)
  439. {
  440. av_freep(&s->f_linebuf);
  441. av_freep(&s->i_linebuf);
  442. }