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.

582 lines
15KB

  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_float(float *p, int i0, int i1)
  125. {
  126. int i;
  127. if (i1 == i0 + 1)
  128. return;
  129. extend97_float(p, i0, i1);
  130. i0++; i1++;
  131. for (i = i0/2 - 2; i < i1/2 + 1; i++)
  132. p[2*i+1] -= 1.586134 * (p[2*i] + p[2*i+2]);
  133. for (i = i0/2 - 1; i < i1/2 + 1; i++)
  134. p[2*i] -= 0.052980 * (p[2*i-1] + p[2*i+1]);
  135. for (i = i0/2 - 1; i < i1/2; i++)
  136. p[2*i+1] += 0.882911 * (p[2*i] + p[2*i+2]);
  137. for (i = i0/2; i < i1/2; i++)
  138. p[2*i] += 0.443506 * (p[2*i-1] + p[2*i+1]);
  139. }
  140. static void dwt_encode97_float(DWTContext *s, float *t)
  141. {
  142. int lev,
  143. w = s->linelen[s->ndeclevels-1][0];
  144. float *line = s->f_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. float *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_float(line, mh, mh + lh);
  160. // copy back and deinterleave
  161. for (i = mh; i < lh; i+=2, j++)
  162. t[w*lp + j] = F_LFTG_X * l[i] / 2;
  163. for (i = 1-mh; i < lh; i+=2, j++)
  164. t[w*lp + j] = F_LFTG_K * l[i] / 2;
  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_float(line, mv, mv + lv);
  173. // copy back and deinterleave
  174. for (i = mv; i < lv; i+=2, j++)
  175. t[w*j + lp] = F_LFTG_X * l[i] / 2;
  176. for (i = 1-mv; i < lv; i+=2, j++)
  177. t[w*j + lp] = F_LFTG_K * l[i] / 2;
  178. }
  179. }
  180. }
  181. static void sd_1d97_int(int *p, int i0, int i1)
  182. {
  183. int i;
  184. if (i1 == i0 + 1)
  185. return;
  186. extend97_int(p, i0, i1);
  187. i0++; i1++;
  188. for (i = i0/2 - 2; i < i1/2 + 1; i++)
  189. p[2 * i + 1] -= (I_LFTG_ALPHA * (p[2 * i] + p[2 * i + 2]) + (1 << 15)) >> 16;
  190. for (i = i0/2 - 1; i < i1/2 + 1; i++)
  191. p[2 * i] -= (I_LFTG_BETA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
  192. for (i = i0/2 - 1; i < i1/2; i++)
  193. p[2 * i + 1] += (I_LFTG_GAMMA * (p[2 * i] + p[2 * i + 2]) + (1 << 15)) >> 16;
  194. for (i = i0/2; i < i1/2; i++)
  195. p[2 * i] += (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
  196. }
  197. static void dwt_encode97_int(DWTContext *s, int *t)
  198. {
  199. int lev,
  200. w = s->linelen[s->ndeclevels-1][0];
  201. int *line = s->i_linebuf;
  202. line += 5;
  203. for (lev = s->ndeclevels-1; lev >= 0; lev--){
  204. int lh = s->linelen[lev][0],
  205. lv = s->linelen[lev][1],
  206. mh = s->mod[lev][0],
  207. mv = s->mod[lev][1],
  208. lp;
  209. int *l;
  210. // HOR_SD
  211. l = line + mh;
  212. for (lp = 0; lp < lv; lp++){
  213. int i, j = 0;
  214. for (i = 0; i < lh; i++)
  215. l[i] = t[w*lp + i];
  216. sd_1d97_int(line, mh, mh + lh);
  217. // copy back and deinterleave
  218. for (i = mh; i < lh; i+=2, j++)
  219. t[w*lp + j] = ((l[i] * I_LFTG_X) + (1 << 16)) >> 17;
  220. for (i = 1-mh; i < lh; i+=2, j++)
  221. t[w*lp + j] = ((l[i] * I_LFTG_K) + (1 << 16)) >> 17;
  222. }
  223. // VER_SD
  224. l = line + mv;
  225. for (lp = 0; lp < lh; lp++) {
  226. int i, j = 0;
  227. for (i = 0; i < lv; i++)
  228. l[i] = t[w*i + lp];
  229. sd_1d97_int(line, mv, mv + lv);
  230. // copy back and deinterleave
  231. for (i = mv; i < lv; i+=2, j++)
  232. t[w*j + lp] = ((l[i] * I_LFTG_X) + (1 << 16)) >> 17;
  233. for (i = 1-mv; i < lv; i+=2, j++)
  234. t[w*j + lp] = ((l[i] * I_LFTG_K) + (1 << 16)) >> 17;
  235. }
  236. }
  237. }
  238. static void sr_1d53(int *p, int i0, int i1)
  239. {
  240. int i;
  241. if (i1 == i0 + 1)
  242. return;
  243. extend53(p, i0, i1);
  244. for (i = i0 / 2; i < i1 / 2 + 1; i++)
  245. p[2 * i] -= (p[2 * i - 1] + p[2 * i + 1] + 2) >> 2;
  246. for (i = i0 / 2; i < i1 / 2; i++)
  247. p[2 * i + 1] += (p[2 * i] + p[2 * i + 2]) >> 1;
  248. }
  249. static void dwt_decode53(DWTContext *s, int *t)
  250. {
  251. int lev;
  252. int w = s->linelen[s->ndeclevels - 1][0];
  253. int32_t *line = s->i_linebuf;
  254. line += 3;
  255. for (lev = 0; lev < s->ndeclevels; lev++) {
  256. int lh = s->linelen[lev][0],
  257. lv = s->linelen[lev][1],
  258. mh = s->mod[lev][0],
  259. mv = s->mod[lev][1],
  260. lp;
  261. int *l;
  262. // HOR_SD
  263. l = line + mh;
  264. for (lp = 0; lp < lv; lp++) {
  265. int i, j = 0;
  266. // copy with interleaving
  267. for (i = mh; i < lh; i += 2, j++)
  268. l[i] = t[w * lp + j];
  269. for (i = 1 - mh; i < lh; i += 2, j++)
  270. l[i] = t[w * lp + j];
  271. sr_1d53(line, mh, mh + lh);
  272. for (i = 0; i < lh; i++)
  273. t[w * lp + i] = l[i];
  274. }
  275. // VER_SD
  276. l = line + mv;
  277. for (lp = 0; lp < lh; lp++) {
  278. int i, j = 0;
  279. // copy with interleaving
  280. for (i = mv; i < lv; i += 2, j++)
  281. l[i] = t[w * j + lp];
  282. for (i = 1 - mv; i < lv; i += 2, j++)
  283. l[i] = t[w * j + lp];
  284. sr_1d53(line, mv, mv + lv);
  285. for (i = 0; i < lv; i++)
  286. t[w * i + lp] = l[i];
  287. }
  288. }
  289. }
  290. static void sr_1d97_float(float *p, int i0, int i1)
  291. {
  292. int i;
  293. if (i1 == i0 + 1)
  294. return;
  295. extend97_float(p, i0, i1);
  296. for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
  297. p[2 * i] -= F_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]);
  298. /* step 4 */
  299. for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
  300. p[2 * i + 1] -= F_LFTG_GAMMA * (p[2 * i] + p[2 * i + 2]);
  301. /*step 5*/
  302. for (i = i0 / 2; i < i1 / 2 + 1; i++)
  303. p[2 * i] += F_LFTG_BETA * (p[2 * i - 1] + p[2 * i + 1]);
  304. /* step 6 */
  305. for (i = i0 / 2; i < i1 / 2; i++)
  306. p[2 * i + 1] += F_LFTG_ALPHA * (p[2 * i] + p[2 * i + 2]);
  307. }
  308. static void dwt_decode97_float(DWTContext *s, float *t)
  309. {
  310. int lev;
  311. int w = s->linelen[s->ndeclevels - 1][0];
  312. float *line = s->f_linebuf;
  313. float *data = t;
  314. /* position at index O of line range [0-5,w+5] cf. extend function */
  315. line += 5;
  316. for (lev = 0; lev < s->ndeclevels; lev++) {
  317. int lh = s->linelen[lev][0],
  318. lv = s->linelen[lev][1],
  319. mh = s->mod[lev][0],
  320. mv = s->mod[lev][1],
  321. lp;
  322. float *l;
  323. // HOR_SD
  324. l = line + mh;
  325. for (lp = 0; lp < lv; lp++) {
  326. int i, j = 0;
  327. // copy with interleaving
  328. for (i = mh; i < lh; i += 2, j++)
  329. l[i] = data[w * lp + j] * F_LFTG_K;
  330. for (i = 1 - mh; i < lh; i += 2, j++)
  331. l[i] = data[w * lp + j] * F_LFTG_X;
  332. sr_1d97_float(line, mh, mh + lh);
  333. for (i = 0; i < lh; i++)
  334. data[w * lp + i] = l[i];
  335. }
  336. // VER_SD
  337. l = line + mv;
  338. for (lp = 0; lp < lh; lp++) {
  339. int i, j = 0;
  340. // copy with interleaving
  341. for (i = mv; i < lv; i += 2, j++)
  342. l[i] = data[w * j + lp] * F_LFTG_K;
  343. for (i = 1 - mv; i < lv; i += 2, j++)
  344. l[i] = data[w * j + lp] * F_LFTG_X;
  345. sr_1d97_float(line, mv, mv + lv);
  346. for (i = 0; i < lv; i++)
  347. data[w * i + lp] = l[i];
  348. }
  349. }
  350. }
  351. static void sr_1d97_int(int32_t *p, int i0, int i1)
  352. {
  353. int i;
  354. if (i1 == i0 + 1)
  355. return;
  356. extend97_int(p, i0, i1);
  357. for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
  358. p[2 * i] -= (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
  359. /* step 4 */
  360. for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
  361. p[2 * i + 1] -= (I_LFTG_GAMMA * (p[2 * i] + p[2 * i + 2]) + (1 << 15)) >> 16;
  362. /*step 5*/
  363. for (i = i0 / 2; i < i1 / 2 + 1; i++)
  364. p[2 * i] += (I_LFTG_BETA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
  365. /* step 6 */
  366. for (i = i0 / 2; i < i1 / 2; i++)
  367. p[2 * i + 1] += (I_LFTG_ALPHA * (p[2 * i] + p[2 * i + 2]) + (1 << 15)) >> 16;
  368. }
  369. static void dwt_decode97_int(DWTContext *s, int32_t *t)
  370. {
  371. int lev;
  372. int w = s->linelen[s->ndeclevels - 1][0];
  373. int32_t *line = s->i_linebuf;
  374. int32_t *data = t;
  375. /* position at index O of line range [0-5,w+5] cf. extend function */
  376. line += 5;
  377. for (lev = 0; lev < s->ndeclevels; lev++) {
  378. int lh = s->linelen[lev][0],
  379. lv = s->linelen[lev][1],
  380. mh = s->mod[lev][0],
  381. mv = s->mod[lev][1],
  382. lp;
  383. int32_t *l;
  384. // HOR_SD
  385. l = line + mh;
  386. for (lp = 0; lp < lv; lp++) {
  387. int i, j = 0;
  388. // rescale with interleaving
  389. for (i = mh; i < lh; i += 2, j++)
  390. l[i] = ((data[w * lp + j] * I_LFTG_K) + (1 << 15)) >> 16;
  391. for (i = 1 - mh; i < lh; i += 2, j++)
  392. l[i] = ((data[w * lp + j] * I_LFTG_X) + (1 << 15)) >> 16;
  393. sr_1d97_int(line, mh, mh + lh);
  394. for (i = 0; i < lh; i++)
  395. data[w * lp + i] = l[i];
  396. }
  397. // VER_SD
  398. l = line + mv;
  399. for (lp = 0; lp < lh; lp++) {
  400. int i, j = 0;
  401. // rescale with interleaving
  402. for (i = mv; i < lv; i += 2, j++)
  403. l[i] = ((data[w * j + lp] * I_LFTG_K) + (1 << 15)) >> 16;
  404. for (i = 1 - mv; i < lv; i += 2, j++)
  405. l[i] = ((data[w * j + lp] * I_LFTG_X) + (1 << 15)) >> 16;
  406. sr_1d97_int(line, mv, mv + lv);
  407. for (i = 0; i < lv; i++)
  408. data[w * i + lp] = l[i];
  409. }
  410. }
  411. }
  412. int ff_jpeg2000_dwt_init(DWTContext *s, uint16_t border[2][2],
  413. int decomp_levels, int type)
  414. {
  415. int i, j, lev = decomp_levels, maxlen,
  416. b[2][2];
  417. s->ndeclevels = decomp_levels;
  418. s->type = type;
  419. for (i = 0; i < 2; i++)
  420. for (j = 0; j < 2; j++)
  421. b[i][j] = border[i][j];
  422. maxlen = FFMAX(b[0][1] - b[0][0],
  423. b[1][1] - b[1][0]);
  424. while (--lev >= 0)
  425. for (i = 0; i < 2; i++) {
  426. s->linelen[lev][i] = b[i][1] - b[i][0];
  427. s->mod[lev][i] = b[i][0] & 1;
  428. for (j = 0; j < 2; j++)
  429. b[i][j] = (b[i][j] + 1) >> 1;
  430. }
  431. switch (type) {
  432. case FF_DWT97:
  433. s->f_linebuf = av_malloc_array((maxlen + 12), sizeof(*s->f_linebuf));
  434. if (!s->f_linebuf)
  435. return AVERROR(ENOMEM);
  436. break;
  437. case FF_DWT97_INT:
  438. s->i_linebuf = av_malloc_array((maxlen + 12), sizeof(*s->i_linebuf));
  439. if (!s->i_linebuf)
  440. return AVERROR(ENOMEM);
  441. break;
  442. case FF_DWT53:
  443. s->i_linebuf = av_malloc_array((maxlen + 6), sizeof(*s->i_linebuf));
  444. if (!s->i_linebuf)
  445. return AVERROR(ENOMEM);
  446. break;
  447. default:
  448. return -1;
  449. }
  450. return 0;
  451. }
  452. int ff_dwt_encode(DWTContext *s, void *t)
  453. {
  454. switch(s->type){
  455. case FF_DWT97:
  456. dwt_encode97_float(s, t); break;
  457. case FF_DWT97_INT:
  458. dwt_encode97_int(s, t); break;
  459. case FF_DWT53:
  460. dwt_encode53(s, t); break;
  461. default:
  462. return -1;
  463. }
  464. return 0;
  465. }
  466. int ff_dwt_decode(DWTContext *s, void *t)
  467. {
  468. if (s->ndeclevels == 0)
  469. return 0;
  470. switch (s->type) {
  471. case FF_DWT97:
  472. dwt_decode97_float(s, t);
  473. break;
  474. case FF_DWT97_INT:
  475. dwt_decode97_int(s, t);
  476. break;
  477. case FF_DWT53:
  478. dwt_decode53(s, t);
  479. break;
  480. default:
  481. return -1;
  482. }
  483. return 0;
  484. }
  485. void ff_dwt_destroy(DWTContext *s)
  486. {
  487. av_freep(&s->f_linebuf);
  488. av_freep(&s->i_linebuf);
  489. }