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.

625 lines
16KB

  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/avassert.h"
  27. #include "libavutil/common.h"
  28. #include "libavutil/mem.h"
  29. #include "jpeg2000dwt.h"
  30. #include "internal.h"
  31. /* Defines for 9/7 DWT lifting parameters.
  32. * Parameters are in float. */
  33. #define F_LFTG_ALPHA 1.586134342059924f
  34. #define F_LFTG_BETA 0.052980118572961f
  35. #define F_LFTG_GAMMA 0.882911075530934f
  36. #define F_LFTG_DELTA 0.443506852043971f
  37. /* Lifting parameters in integer format.
  38. * Computed as param = (float param) * (1 << 16) */
  39. #define I_LFTG_ALPHA 103949ll
  40. #define I_LFTG_BETA 3472ll
  41. #define I_LFTG_GAMMA 57862ll
  42. #define I_LFTG_DELTA 29066ll
  43. #define I_LFTG_K 80621ll
  44. #define I_LFTG_X 53274ll
  45. #define I_PRESHIFT 8
  46. static inline void extend53(int *p, int i0, int i1)
  47. {
  48. p[i0 - 1] = p[i0 + 1];
  49. p[i1] = p[i1 - 2];
  50. p[i0 - 2] = p[i0 + 2];
  51. p[i1 + 1] = p[i1 - 3];
  52. }
  53. static inline void extend97_float(float *p, int i0, int i1)
  54. {
  55. int i;
  56. for (i = 1; i <= 4; i++) {
  57. p[i0 - i] = p[i0 + i];
  58. p[i1 + i - 1] = p[i1 - i - 1];
  59. }
  60. }
  61. static inline void extend97_int(int32_t *p, int i0, int i1)
  62. {
  63. int i;
  64. for (i = 1; i <= 4; i++) {
  65. p[i0 - i] = p[i0 + i];
  66. p[i1 + i - 1] = p[i1 - i - 1];
  67. }
  68. }
  69. static void sd_1d53(int *p, int i0, int i1)
  70. {
  71. int i;
  72. if (i1 <= i0 + 1) {
  73. if (i0 == 1)
  74. p[1] <<= 1;
  75. return;
  76. }
  77. extend53(p, i0, i1);
  78. for (i = ((i0+1)>>1) - 1; i < (i1+1)>>1; i++)
  79. p[2*i+1] -= (p[2*i] + p[2*i+2]) >> 1;
  80. for (i = ((i0+1)>>1); i < (i1+1)>>1; 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. // VER_SD
  97. l = line + mv;
  98. for (lp = 0; lp < lh; lp++) {
  99. int i, j = 0;
  100. for (i = 0; i < lv; i++)
  101. l[i] = t[w*i + lp];
  102. sd_1d53(line, mv, mv + lv);
  103. // copy back and deinterleave
  104. for (i = mv; i < lv; i+=2, j++)
  105. t[w*j + lp] = l[i];
  106. for (i = 1-mv; i < lv; i+=2, j++)
  107. t[w*j + lp] = l[i];
  108. }
  109. // HOR_SD
  110. l = line + mh;
  111. for (lp = 0; lp < lv; lp++){
  112. int i, j = 0;
  113. for (i = 0; i < lh; i++)
  114. l[i] = t[w*lp + i];
  115. sd_1d53(line, mh, mh + lh);
  116. // copy back and deinterleave
  117. for (i = mh; i < lh; i+=2, j++)
  118. t[w*lp + j] = l[i];
  119. for (i = 1-mh; i < lh; i+=2, j++)
  120. t[w*lp + j] = 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. if (i0 == 1)
  129. p[1] *= F_LFTG_X * 2;
  130. else
  131. p[0] *= F_LFTG_K;
  132. return;
  133. }
  134. extend97_float(p, i0, i1);
  135. i0++; i1++;
  136. for (i = (i0>>1) - 2; i < (i1>>1) + 1; i++)
  137. p[2*i+1] -= 1.586134 * (p[2*i] + p[2*i+2]);
  138. for (i = (i0>>1) - 1; i < (i1>>1) + 1; i++)
  139. p[2*i] -= 0.052980 * (p[2*i-1] + p[2*i+1]);
  140. for (i = (i0>>1) - 1; i < (i1>>1); i++)
  141. p[2*i+1] += 0.882911 * (p[2*i] + p[2*i+2]);
  142. for (i = (i0>>1); i < (i1>>1); i++)
  143. p[2*i] += 0.443506 * (p[2*i-1] + p[2*i+1]);
  144. }
  145. static void dwt_encode97_float(DWTContext *s, float *t)
  146. {
  147. int lev,
  148. w = s->linelen[s->ndeclevels-1][0];
  149. float *line = s->f_linebuf;
  150. line += 5;
  151. for (lev = s->ndeclevels-1; lev >= 0; lev--){
  152. int lh = s->linelen[lev][0],
  153. lv = s->linelen[lev][1],
  154. mh = s->mod[lev][0],
  155. mv = s->mod[lev][1],
  156. lp;
  157. float *l;
  158. // HOR_SD
  159. l = line + mh;
  160. for (lp = 0; lp < lv; lp++){
  161. int i, j = 0;
  162. for (i = 0; i < lh; i++)
  163. l[i] = t[w*lp + i];
  164. sd_1d97_float(line, mh, mh + lh);
  165. // copy back and deinterleave
  166. for (i = mh; i < lh; i+=2, j++)
  167. t[w*lp + j] = l[i];
  168. for (i = 1-mh; i < lh; i+=2, j++)
  169. t[w*lp + j] = l[i];
  170. }
  171. // VER_SD
  172. l = line + mv;
  173. for (lp = 0; lp < lh; lp++) {
  174. int i, j = 0;
  175. for (i = 0; i < lv; i++)
  176. l[i] = t[w*i + lp];
  177. sd_1d97_float(line, mv, mv + lv);
  178. // copy back and deinterleave
  179. for (i = mv; i < lv; i+=2, j++)
  180. t[w*j + lp] = l[i];
  181. for (i = 1-mv; i < lv; i+=2, j++)
  182. t[w*j + lp] = l[i];
  183. }
  184. }
  185. }
  186. static void sd_1d97_int(int *p, int i0, int i1)
  187. {
  188. int i;
  189. if (i1 <= i0 + 1) {
  190. if (i0 == 1)
  191. p[1] = (p[1] * I_LFTG_X + (1<<14)) >> 15;
  192. else
  193. p[0] = (p[0] * I_LFTG_K + (1<<15)) >> 16;
  194. return;
  195. }
  196. extend97_int(p, i0, i1);
  197. i0++; i1++;
  198. for (i = (i0>>1) - 2; i < (i1>>1) + 1; i++)
  199. p[2 * i + 1] -= (I_LFTG_ALPHA * (p[2 * i] + p[2 * i + 2]) + (1 << 15)) >> 16;
  200. for (i = (i0>>1) - 1; i < (i1>>1) + 1; i++)
  201. p[2 * i] -= (I_LFTG_BETA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
  202. for (i = (i0>>1) - 1; i < (i1>>1); i++)
  203. p[2 * i + 1] += (I_LFTG_GAMMA * (p[2 * i] + p[2 * i + 2]) + (1 << 15)) >> 16;
  204. for (i = (i0>>1); i < (i1>>1); i++)
  205. p[2 * i] += (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
  206. }
  207. static void dwt_encode97_int(DWTContext *s, int *t)
  208. {
  209. int lev;
  210. int w = s->linelen[s->ndeclevels-1][0];
  211. int h = s->linelen[s->ndeclevels-1][1];
  212. int i;
  213. int *line = s->i_linebuf;
  214. line += 5;
  215. for (i = 0; i < w * h; i++)
  216. t[i] <<= I_PRESHIFT;
  217. for (lev = s->ndeclevels-1; lev >= 0; lev--){
  218. int lh = s->linelen[lev][0],
  219. lv = s->linelen[lev][1],
  220. mh = s->mod[lev][0],
  221. mv = s->mod[lev][1],
  222. lp;
  223. int *l;
  224. // VER_SD
  225. l = line + mv;
  226. for (lp = 0; lp < lh; lp++) {
  227. int i, j = 0;
  228. for (i = 0; i < lv; i++)
  229. l[i] = t[w*i + lp];
  230. sd_1d97_int(line, mv, mv + lv);
  231. // copy back and deinterleave
  232. for (i = mv; i < lv; i+=2, j++)
  233. t[w*j + lp] = ((l[i] * I_LFTG_X) + (1 << 15)) >> 16;
  234. for (i = 1-mv; i < lv; i+=2, j++)
  235. t[w*j + lp] = l[i];
  236. }
  237. // HOR_SD
  238. l = line + mh;
  239. for (lp = 0; lp < lv; lp++){
  240. int i, j = 0;
  241. for (i = 0; i < lh; i++)
  242. l[i] = t[w*lp + i];
  243. sd_1d97_int(line, mh, mh + lh);
  244. // copy back and deinterleave
  245. for (i = mh; i < lh; i+=2, j++)
  246. t[w*lp + j] = ((l[i] * I_LFTG_X) + (1 << 15)) >> 16;
  247. for (i = 1-mh; i < lh; i+=2, j++)
  248. t[w*lp + j] = l[i];
  249. }
  250. }
  251. for (i = 0; i < w * h; i++)
  252. t[i] = (t[i] + ((1<<I_PRESHIFT)>>1)) >> I_PRESHIFT;
  253. }
  254. static void sr_1d53(int *p, int i0, int i1)
  255. {
  256. int i;
  257. if (i1 <= i0 + 1) {
  258. if (i0 == 1)
  259. p[1] >>= 1;
  260. return;
  261. }
  262. extend53(p, i0, i1);
  263. for (i = (i0 >> 1); i < (i1 >> 1) + 1; i++)
  264. p[2 * i] -= (p[2 * i - 1] + p[2 * i + 1] + 2) >> 2;
  265. for (i = (i0 >> 1); i < (i1 >> 1); i++)
  266. p[2 * i + 1] += (p[2 * i] + p[2 * i + 2]) >> 1;
  267. }
  268. static void dwt_decode53(DWTContext *s, int *t)
  269. {
  270. int lev;
  271. int w = s->linelen[s->ndeclevels - 1][0];
  272. int32_t *line = s->i_linebuf;
  273. line += 3;
  274. for (lev = 0; lev < s->ndeclevels; lev++) {
  275. int lh = s->linelen[lev][0],
  276. lv = s->linelen[lev][1],
  277. mh = s->mod[lev][0],
  278. mv = s->mod[lev][1],
  279. lp;
  280. int *l;
  281. // HOR_SD
  282. l = line + mh;
  283. for (lp = 0; lp < lv; lp++) {
  284. int i, j = 0;
  285. // copy with interleaving
  286. for (i = mh; i < lh; i += 2, j++)
  287. l[i] = t[w * lp + j];
  288. for (i = 1 - mh; i < lh; i += 2, j++)
  289. l[i] = t[w * lp + j];
  290. sr_1d53(line, mh, mh + lh);
  291. for (i = 0; i < lh; i++)
  292. t[w * lp + i] = l[i];
  293. }
  294. // VER_SD
  295. l = line + mv;
  296. for (lp = 0; lp < lh; lp++) {
  297. int i, j = 0;
  298. // copy with interleaving
  299. for (i = mv; i < lv; i += 2, j++)
  300. l[i] = t[w * j + lp];
  301. for (i = 1 - mv; i < lv; i += 2, j++)
  302. l[i] = t[w * j + lp];
  303. sr_1d53(line, mv, mv + lv);
  304. for (i = 0; i < lv; i++)
  305. t[w * i + lp] = l[i];
  306. }
  307. }
  308. }
  309. static void sr_1d97_float(float *p, int i0, int i1)
  310. {
  311. int i;
  312. if (i1 <= i0 + 1) {
  313. if (i0 == 1)
  314. p[1] *= F_LFTG_K/2;
  315. else
  316. p[0] *= F_LFTG_X;
  317. return;
  318. }
  319. extend97_float(p, i0, i1);
  320. for (i = (i0 >> 1) - 1; i < (i1 >> 1) + 2; i++)
  321. p[2 * i] -= F_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]);
  322. /* step 4 */
  323. for (i = (i0 >> 1) - 1; i < (i1 >> 1) + 1; i++)
  324. p[2 * i + 1] -= F_LFTG_GAMMA * (p[2 * i] + p[2 * i + 2]);
  325. /*step 5*/
  326. for (i = (i0 >> 1); i < (i1 >> 1) + 1; i++)
  327. p[2 * i] += F_LFTG_BETA * (p[2 * i - 1] + p[2 * i + 1]);
  328. /* step 6 */
  329. for (i = (i0 >> 1); i < (i1 >> 1); i++)
  330. p[2 * i + 1] += F_LFTG_ALPHA * (p[2 * i] + p[2 * i + 2]);
  331. }
  332. static void dwt_decode97_float(DWTContext *s, float *t)
  333. {
  334. int lev;
  335. int w = s->linelen[s->ndeclevels - 1][0];
  336. float *line = s->f_linebuf;
  337. float *data = t;
  338. /* position at index O of line range [0-5,w+5] cf. extend function */
  339. line += 5;
  340. for (lev = 0; lev < s->ndeclevels; lev++) {
  341. int lh = s->linelen[lev][0],
  342. lv = s->linelen[lev][1],
  343. mh = s->mod[lev][0],
  344. mv = s->mod[lev][1],
  345. lp;
  346. float *l;
  347. // HOR_SD
  348. l = line + mh;
  349. for (lp = 0; lp < lv; lp++) {
  350. int i, j = 0;
  351. // copy with interleaving
  352. for (i = mh; i < lh; i += 2, j++)
  353. l[i] = data[w * lp + j];
  354. for (i = 1 - mh; i < lh; i += 2, j++)
  355. l[i] = data[w * lp + j];
  356. sr_1d97_float(line, mh, mh + lh);
  357. for (i = 0; i < lh; i++)
  358. data[w * lp + i] = l[i];
  359. }
  360. // VER_SD
  361. l = line + mv;
  362. for (lp = 0; lp < lh; lp++) {
  363. int i, j = 0;
  364. // copy with interleaving
  365. for (i = mv; i < lv; i += 2, j++)
  366. l[i] = data[w * j + lp];
  367. for (i = 1 - mv; i < lv; i += 2, j++)
  368. l[i] = data[w * j + lp];
  369. sr_1d97_float(line, mv, mv + lv);
  370. for (i = 0; i < lv; i++)
  371. data[w * i + lp] = l[i];
  372. }
  373. }
  374. }
  375. static void sr_1d97_int(int32_t *p, int i0, int i1)
  376. {
  377. int i;
  378. if (i1 <= i0 + 1) {
  379. if (i0 == 1)
  380. p[1] = (p[1] * I_LFTG_K + (1<<16)) >> 17;
  381. else
  382. p[0] = (p[0] * I_LFTG_X + (1<<15)) >> 16;
  383. return;
  384. }
  385. extend97_int(p, i0, i1);
  386. for (i = (i0 >> 1) - 1; i < (i1 >> 1) + 2; i++)
  387. p[2 * i] -= (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
  388. /* step 4 */
  389. for (i = (i0 >> 1) - 1; i < (i1 >> 1) + 1; i++)
  390. p[2 * i + 1] -= (I_LFTG_GAMMA * (p[2 * i] + p[2 * i + 2]) + (1 << 15)) >> 16;
  391. /*step 5*/
  392. for (i = (i0 >> 1); i < (i1 >> 1) + 1; i++)
  393. p[2 * i] += (I_LFTG_BETA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
  394. /* step 6 */
  395. for (i = (i0 >> 1); i < (i1 >> 1); i++)
  396. p[2 * i + 1] += (I_LFTG_ALPHA * (p[2 * i] + p[2 * i + 2]) + (1 << 15)) >> 16;
  397. }
  398. static void dwt_decode97_int(DWTContext *s, int32_t *t)
  399. {
  400. int lev;
  401. int w = s->linelen[s->ndeclevels - 1][0];
  402. int h = s->linelen[s->ndeclevels - 1][1];
  403. int i;
  404. int32_t *line = s->i_linebuf;
  405. int32_t *data = t;
  406. /* position at index O of line range [0-5,w+5] cf. extend function */
  407. line += 5;
  408. for (i = 0; i < w * h; i++)
  409. data[i] *= 1 << I_PRESHIFT;
  410. for (lev = 0; lev < s->ndeclevels; lev++) {
  411. int lh = s->linelen[lev][0],
  412. lv = s->linelen[lev][1],
  413. mh = s->mod[lev][0],
  414. mv = s->mod[lev][1],
  415. lp;
  416. int32_t *l;
  417. // HOR_SD
  418. l = line + mh;
  419. for (lp = 0; lp < lv; lp++) {
  420. int i, j = 0;
  421. // rescale with interleaving
  422. for (i = mh; i < lh; i += 2, j++)
  423. l[i] = ((data[w * lp + j] * I_LFTG_K) + (1 << 15)) >> 16;
  424. for (i = 1 - mh; i < lh; i += 2, j++)
  425. l[i] = data[w * lp + j];
  426. sr_1d97_int(line, mh, mh + lh);
  427. for (i = 0; i < lh; i++)
  428. data[w * lp + i] = l[i];
  429. }
  430. // VER_SD
  431. l = line + mv;
  432. for (lp = 0; lp < lh; lp++) {
  433. int i, j = 0;
  434. // rescale with interleaving
  435. for (i = mv; i < lv; i += 2, j++)
  436. l[i] = ((data[w * j + lp] * I_LFTG_K) + (1 << 15)) >> 16;
  437. for (i = 1 - mv; i < lv; i += 2, j++)
  438. l[i] = data[w * j + lp];
  439. sr_1d97_int(line, mv, mv + lv);
  440. for (i = 0; i < lv; i++)
  441. data[w * i + lp] = l[i];
  442. }
  443. }
  444. for (i = 0; i < w * h; i++)
  445. data[i] = (data[i] + ((1<<I_PRESHIFT)>>1)) >> I_PRESHIFT;
  446. }
  447. int ff_jpeg2000_dwt_init(DWTContext *s, int border[2][2],
  448. int decomp_levels, int type)
  449. {
  450. int i, j, lev = decomp_levels, maxlen,
  451. b[2][2];
  452. s->ndeclevels = decomp_levels;
  453. s->type = type;
  454. for (i = 0; i < 2; i++)
  455. for (j = 0; j < 2; j++)
  456. b[i][j] = border[i][j];
  457. maxlen = FFMAX(b[0][1] - b[0][0],
  458. b[1][1] - b[1][0]);
  459. while (--lev >= 0)
  460. for (i = 0; i < 2; i++) {
  461. s->linelen[lev][i] = b[i][1] - b[i][0];
  462. s->mod[lev][i] = b[i][0] & 1;
  463. for (j = 0; j < 2; j++)
  464. b[i][j] = (b[i][j] + 1) >> 1;
  465. }
  466. switch (type) {
  467. case FF_DWT97:
  468. s->f_linebuf = av_malloc_array((maxlen + 12), sizeof(*s->f_linebuf));
  469. if (!s->f_linebuf)
  470. return AVERROR(ENOMEM);
  471. break;
  472. case FF_DWT97_INT:
  473. s->i_linebuf = av_malloc_array((maxlen + 12), sizeof(*s->i_linebuf));
  474. if (!s->i_linebuf)
  475. return AVERROR(ENOMEM);
  476. break;
  477. case FF_DWT53:
  478. s->i_linebuf = av_malloc_array((maxlen + 6), sizeof(*s->i_linebuf));
  479. if (!s->i_linebuf)
  480. return AVERROR(ENOMEM);
  481. break;
  482. default:
  483. return -1;
  484. }
  485. return 0;
  486. }
  487. int ff_dwt_encode(DWTContext *s, void *t)
  488. {
  489. if (s->ndeclevels == 0)
  490. return 0;
  491. switch(s->type){
  492. case FF_DWT97:
  493. dwt_encode97_float(s, t); break;
  494. case FF_DWT97_INT:
  495. dwt_encode97_int(s, t); break;
  496. case FF_DWT53:
  497. dwt_encode53(s, t); break;
  498. default:
  499. return -1;
  500. }
  501. return 0;
  502. }
  503. int ff_dwt_decode(DWTContext *s, void *t)
  504. {
  505. if (s->ndeclevels == 0)
  506. return 0;
  507. switch (s->type) {
  508. case FF_DWT97:
  509. dwt_decode97_float(s, t);
  510. break;
  511. case FF_DWT97_INT:
  512. dwt_decode97_int(s, t);
  513. break;
  514. case FF_DWT53:
  515. dwt_decode53(s, t);
  516. break;
  517. default:
  518. return -1;
  519. }
  520. return 0;
  521. }
  522. void ff_dwt_destroy(DWTContext *s)
  523. {
  524. av_freep(&s->f_linebuf);
  525. av_freep(&s->i_linebuf);
  526. }