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.

728 lines
20KB

  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. #define F_LFTG_K 1.230174104914001f
  38. #define F_LFTG_X 1.625732422f
  39. /* FIXME: Why use 1.625732422 instead of 1/F_LFTG_K?
  40. * Incorrect value in JPEG2000 norm.
  41. * see (ISO/IEC 15444:1 (version 2002) F.3.8.2 */
  42. /* Lifting parameters in integer format.
  43. * Computed as param = (float param) * (1 << 16) */
  44. #define I_LFTG_ALPHA 103949
  45. #define I_LFTG_BETA 3472
  46. #define I_LFTG_GAMMA 57862
  47. #define I_LFTG_DELTA 29066
  48. #define I_LFTG_K 80621
  49. #define I_LFTG_X 106544
  50. static inline void extend53(int *p, int i0, int i1)
  51. {
  52. p[i0 - 1] = p[i0 + 1];
  53. p[i1] = p[i1 - 2];
  54. p[i0 - 2] = p[i0 + 2];
  55. p[i1 + 1] = p[i1 - 3];
  56. }
  57. static inline void extend97_float(float *p, int i0, int i1)
  58. {
  59. int i;
  60. for (i = 1; i <= 4; i++) {
  61. p[i0 - i] = p[i0 + i];
  62. p[i1 + i - 1] = p[i1 - i - 1];
  63. }
  64. }
  65. static inline void extend97_int(int32_t *p, int i0, int i1)
  66. {
  67. int i;
  68. for (i = 1; i <= 4; i++) {
  69. p[i0 - i] = p[i0 + i];
  70. p[i1 + i - 1] = p[i1 - i - 1];
  71. }
  72. }
  73. static void sd_1d53(int *p, int i0, int i1)
  74. {
  75. int i;
  76. if (i1 <= i0 + 1) {
  77. if (i0 == 1)
  78. p[1] <<= 1;
  79. return;
  80. }
  81. extend53(p, i0, i1);
  82. for (i = (i0+1)/2 - 1; i < (i1+1)/2; i++)
  83. p[2*i+1] -= (p[2*i] + p[2*i+2]) >> 1;
  84. for (i = (i0+1)/2; i < (i1+1)/2; i++)
  85. p[2*i] += (p[2*i-1] + p[2*i+1] + 2) >> 2;
  86. }
  87. static void dwt_encode53(DWTContext *s, int *t)
  88. {
  89. int lev,
  90. w = s->linelen[s->ndeclevels-1][0];
  91. int *line = s->i_linebuf;
  92. line += 3;
  93. for (lev = s->ndeclevels-1; lev >= 0; lev--){
  94. int lh = s->linelen[lev][0],
  95. lv = s->linelen[lev][1],
  96. mh = s->mod[lev][0],
  97. mv = s->mod[lev][1],
  98. lp;
  99. int *l;
  100. // VER_SD
  101. l = line + mv;
  102. for (lp = 0; lp < lh; lp++) {
  103. int i, j = 0;
  104. for (i = 0; i < lv; i++)
  105. l[i] = t[w*i + lp];
  106. sd_1d53(line, mv, mv + lv);
  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. // HOR_SD
  114. l = line + mh;
  115. for (lp = 0; lp < lv; lp++){
  116. int i, j = 0;
  117. for (i = 0; i < lh; i++)
  118. l[i] = t[w*lp + i];
  119. sd_1d53(line, mh, mh + lh);
  120. // copy back and deinterleave
  121. for (i = mh; i < lh; i+=2, j++)
  122. t[w*lp + j] = l[i];
  123. for (i = 1-mh; i < lh; i+=2, j++)
  124. t[w*lp + j] = l[i];
  125. }
  126. }
  127. }
  128. static void sd_1d97_float(float *p, int i0, int i1)
  129. {
  130. int i;
  131. if (i1 <= i0 + 1) {
  132. if (i0 == 1)
  133. p[1] *= F_LFTG_X;
  134. else
  135. p[0] *= F_LFTG_K;
  136. return;
  137. }
  138. extend97_float(p, i0, i1);
  139. i0++; i1++;
  140. for (i = i0/2 - 2; i < i1/2 + 1; i++)
  141. p[2*i+1] -= 1.586134 * (p[2*i] + p[2*i+2]);
  142. for (i = i0/2 - 1; i < i1/2 + 1; i++)
  143. p[2*i] -= 0.052980 * (p[2*i-1] + p[2*i+1]);
  144. for (i = i0/2 - 1; i < i1/2; i++)
  145. p[2*i+1] += 0.882911 * (p[2*i] + p[2*i+2]);
  146. for (i = i0/2; i < i1/2; i++)
  147. p[2*i] += 0.443506 * (p[2*i-1] + p[2*i+1]);
  148. }
  149. static void dwt_encode97_float(DWTContext *s, float *t)
  150. {
  151. int lev,
  152. w = s->linelen[s->ndeclevels-1][0];
  153. float *line = s->f_linebuf;
  154. line += 5;
  155. for (lev = s->ndeclevels-1; lev >= 0; lev--){
  156. int lh = s->linelen[lev][0],
  157. lv = s->linelen[lev][1],
  158. mh = s->mod[lev][0],
  159. mv = s->mod[lev][1],
  160. lp;
  161. float *l;
  162. // HOR_SD
  163. l = line + mh;
  164. for (lp = 0; lp < lv; lp++){
  165. int i, j = 0;
  166. for (i = 0; i < lh; i++)
  167. l[i] = t[w*lp + i];
  168. sd_1d97_float(line, mh, mh + lh);
  169. // copy back and deinterleave
  170. for (i = mh; i < lh; i+=2, j++)
  171. t[w*lp + j] = F_LFTG_X * l[i] / 2;
  172. for (i = 1-mh; i < lh; i+=2, j++)
  173. t[w*lp + j] = F_LFTG_K * l[i] / 2;
  174. }
  175. // VER_SD
  176. l = line + mv;
  177. for (lp = 0; lp < lh; lp++) {
  178. int i, j = 0;
  179. for (i = 0; i < lv; i++)
  180. l[i] = t[w*i + lp];
  181. sd_1d97_float(line, mv, mv + lv);
  182. // copy back and deinterleave
  183. for (i = mv; i < lv; i+=2, j++)
  184. t[w*j + lp] = F_LFTG_X * l[i] / 2;
  185. for (i = 1-mv; i < lv; i+=2, j++)
  186. t[w*j + lp] = F_LFTG_K * l[i] / 2;
  187. }
  188. }
  189. }
  190. static void sd_1d97_int(int *p, int i0, int i1)
  191. {
  192. int i;
  193. if (i1 <= i0 + 1) {
  194. if (i0 == 1)
  195. p[1] = (p[1] * I_LFTG_X + (1<<15)) >> 16;
  196. else
  197. p[0] = (p[0] * I_LFTG_K + (1<<15)) >> 16;
  198. return;
  199. }
  200. extend97_int(p, i0, i1);
  201. i0++; i1++;
  202. for (i = i0/2 - 2; i < i1/2 + 1; i++)
  203. p[2 * i + 1] -= (I_LFTG_ALPHA * (p[2 * i] + p[2 * i + 2]) + (1 << 15)) >> 16;
  204. for (i = i0/2 - 1; i < i1/2 + 1; i++)
  205. p[2 * i] -= (I_LFTG_BETA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
  206. for (i = i0/2 - 1; i < i1/2; i++)
  207. p[2 * i + 1] += (I_LFTG_GAMMA * (p[2 * i] + p[2 * i + 2]) + (1 << 15)) >> 16;
  208. for (i = i0/2; i < i1/2; i++)
  209. p[2 * i] += (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
  210. }
  211. static void dwt_encode97_int(DWTContext *s, int *t)
  212. {
  213. int lev,
  214. w = s->linelen[s->ndeclevels-1][0];
  215. int *line = s->i_linebuf;
  216. line += 5;
  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 << 16)) >> 17;
  234. for (i = 1-mv; i < lv; i+=2, j++)
  235. t[w*j + lp] = ((l[i] * I_LFTG_K) + (1 << 16)) >> 17;
  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 << 16)) >> 17;
  247. for (i = 1-mh; i < lh; i+=2, j++)
  248. t[w*lp + j] = ((l[i] * I_LFTG_K) + (1 << 16)) >> 17;
  249. }
  250. }
  251. }
  252. static void sr_1d53(int *p, int i0, int i1)
  253. {
  254. int i;
  255. if (i1 <= i0 + 1) {
  256. if (i0 == 1)
  257. p[1] >>= 1;
  258. return;
  259. }
  260. extend53(p, i0, i1);
  261. for (i = i0 / 2; i < i1 / 2 + 1; i++)
  262. p[2 * i] -= (p[2 * i - 1] + p[2 * i + 1] + 2) >> 2;
  263. for (i = i0 / 2; i < i1 / 2; i++)
  264. p[2 * i + 1] += (p[2 * i] + p[2 * i + 2]) >> 1;
  265. }
  266. static void dwt_decode53(DWTContext *s, int *t)
  267. {
  268. int lev;
  269. int w = s->linelen[s->ndeclevels - 1][0];
  270. int32_t *line = s->i_linebuf;
  271. line += 3;
  272. for (lev = 0; lev < s->ndeclevels; lev++) {
  273. int lh = s->linelen[lev][0],
  274. lv = s->linelen[lev][1],
  275. mh = s->mod[lev][0],
  276. mv = s->mod[lev][1],
  277. lp;
  278. int *l;
  279. // HOR_SD
  280. l = line + mh;
  281. for (lp = 0; lp < lv; lp++) {
  282. int i, j = 0;
  283. // copy with interleaving
  284. for (i = mh; i < lh; i += 2, j++)
  285. l[i] = t[w * lp + j];
  286. for (i = 1 - mh; i < lh; i += 2, j++)
  287. l[i] = t[w * lp + j];
  288. sr_1d53(line, mh, mh + lh);
  289. for (i = 0; i < lh; i++)
  290. t[w * lp + i] = l[i];
  291. }
  292. // VER_SD
  293. l = line + mv;
  294. for (lp = 0; lp < lh; lp++) {
  295. int i, j = 0;
  296. // copy with interleaving
  297. for (i = mv; i < lv; i += 2, j++)
  298. l[i] = t[w * j + lp];
  299. for (i = 1 - mv; i < lv; i += 2, j++)
  300. l[i] = t[w * j + lp];
  301. sr_1d53(line, mv, mv + lv);
  302. for (i = 0; i < lv; i++)
  303. t[w * i + lp] = l[i];
  304. }
  305. }
  306. }
  307. static void sr_1d97_float(float *p, int i0, int i1)
  308. {
  309. int i;
  310. if (i1 <= i0 + 1) {
  311. if (i0 == 1)
  312. p[1] *= F_LFTG_K/2;
  313. else
  314. p[0] *= F_LFTG_X/2;
  315. return;
  316. }
  317. extend97_float(p, i0, i1);
  318. for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
  319. p[2 * i] -= F_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]);
  320. /* step 4 */
  321. for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
  322. p[2 * i + 1] -= F_LFTG_GAMMA * (p[2 * i] + p[2 * i + 2]);
  323. /*step 5*/
  324. for (i = i0 / 2; i < i1 / 2 + 1; i++)
  325. p[2 * i] += F_LFTG_BETA * (p[2 * i - 1] + p[2 * i + 1]);
  326. /* step 6 */
  327. for (i = i0 / 2; i < i1 / 2; i++)
  328. p[2 * i + 1] += F_LFTG_ALPHA * (p[2 * i] + p[2 * i + 2]);
  329. }
  330. static void dwt_decode97_float(DWTContext *s, float *t)
  331. {
  332. int lev;
  333. int w = s->linelen[s->ndeclevels - 1][0];
  334. float *line = s->f_linebuf;
  335. float *data = t;
  336. /* position at index O of line range [0-5,w+5] cf. extend function */
  337. line += 5;
  338. for (lev = 0; lev < s->ndeclevels; lev++) {
  339. int lh = s->linelen[lev][0],
  340. lv = s->linelen[lev][1],
  341. mh = s->mod[lev][0],
  342. mv = s->mod[lev][1],
  343. lp;
  344. float *l;
  345. // HOR_SD
  346. l = line + mh;
  347. for (lp = 0; lp < lv; lp++) {
  348. int i, j = 0;
  349. // copy with interleaving
  350. for (i = mh; i < lh; i += 2, j++)
  351. l[i] = data[w * lp + j] * F_LFTG_K;
  352. for (i = 1 - mh; i < lh; i += 2, j++)
  353. l[i] = data[w * lp + j] * F_LFTG_X;
  354. sr_1d97_float(line, mh, mh + lh);
  355. for (i = 0; i < lh; i++)
  356. data[w * lp + i] = l[i];
  357. }
  358. // VER_SD
  359. l = line + mv;
  360. for (lp = 0; lp < lh; lp++) {
  361. int i, j = 0;
  362. // copy with interleaving
  363. for (i = mv; i < lv; i += 2, j++)
  364. l[i] = data[w * j + lp] * F_LFTG_K;
  365. for (i = 1 - mv; i < lv; i += 2, j++)
  366. l[i] = data[w * j + lp] * F_LFTG_X;
  367. sr_1d97_float(line, mv, mv + lv);
  368. for (i = 0; i < lv; i++)
  369. data[w * i + lp] = l[i];
  370. }
  371. }
  372. }
  373. static void sr_1d97_int(int32_t *p, int i0, int i1)
  374. {
  375. int i;
  376. if (i1 <= i0 + 1) {
  377. if (i0 == 1)
  378. p[1] = (p[1] * I_LFTG_K + (1<<16)) >> 17;
  379. else
  380. p[0] = (p[0] * I_LFTG_X + (1<<16)) >> 17;
  381. return;
  382. }
  383. extend97_int(p, i0, i1);
  384. for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
  385. p[2 * i] -= (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
  386. /* step 4 */
  387. for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
  388. p[2 * i + 1] -= (I_LFTG_GAMMA * (p[2 * i] + p[2 * i + 2]) + (1 << 15)) >> 16;
  389. /*step 5*/
  390. for (i = i0 / 2; i < i1 / 2 + 1; i++)
  391. p[2 * i] += (I_LFTG_BETA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
  392. /* step 6 */
  393. for (i = i0 / 2; i < i1 / 2; i++)
  394. p[2 * i + 1] += (I_LFTG_ALPHA * (p[2 * i] + p[2 * i + 2]) + (1 << 15)) >> 16;
  395. }
  396. static void dwt_decode97_int(DWTContext *s, int32_t *t)
  397. {
  398. int lev;
  399. int w = s->linelen[s->ndeclevels - 1][0];
  400. int32_t *line = s->i_linebuf;
  401. int32_t *data = t;
  402. /* position at index O of line range [0-5,w+5] cf. extend function */
  403. line += 5;
  404. for (lev = 0; lev < s->ndeclevels; lev++) {
  405. int lh = s->linelen[lev][0],
  406. lv = s->linelen[lev][1],
  407. mh = s->mod[lev][0],
  408. mv = s->mod[lev][1],
  409. lp;
  410. int32_t *l;
  411. // HOR_SD
  412. l = line + mh;
  413. for (lp = 0; lp < lv; lp++) {
  414. int i, j = 0;
  415. // rescale with interleaving
  416. for (i = mh; i < lh; i += 2, j++)
  417. l[i] = ((data[w * lp + j] * I_LFTG_K) + (1 << 15)) >> 16;
  418. for (i = 1 - mh; i < lh; i += 2, j++)
  419. l[i] = ((data[w * lp + j] * I_LFTG_X) + (1 << 15)) >> 16;
  420. sr_1d97_int(line, mh, mh + lh);
  421. for (i = 0; i < lh; i++)
  422. data[w * lp + i] = l[i];
  423. }
  424. // VER_SD
  425. l = line + mv;
  426. for (lp = 0; lp < lh; lp++) {
  427. int i, j = 0;
  428. // rescale with interleaving
  429. for (i = mv; i < lv; i += 2, j++)
  430. l[i] = ((data[w * j + lp] * I_LFTG_K) + (1 << 15)) >> 16;
  431. for (i = 1 - mv; i < lv; i += 2, j++)
  432. l[i] = ((data[w * j + lp] * I_LFTG_X) + (1 << 15)) >> 16;
  433. sr_1d97_int(line, mv, mv + lv);
  434. for (i = 0; i < lv; i++)
  435. data[w * i + lp] = l[i];
  436. }
  437. }
  438. }
  439. int ff_jpeg2000_dwt_init(DWTContext *s, uint16_t border[2][2],
  440. int decomp_levels, int type)
  441. {
  442. int i, j, lev = decomp_levels, maxlen,
  443. b[2][2];
  444. s->ndeclevels = decomp_levels;
  445. s->type = type;
  446. for (i = 0; i < 2; i++)
  447. for (j = 0; j < 2; j++)
  448. b[i][j] = border[i][j];
  449. maxlen = FFMAX(b[0][1] - b[0][0],
  450. b[1][1] - b[1][0]);
  451. while (--lev >= 0)
  452. for (i = 0; i < 2; i++) {
  453. s->linelen[lev][i] = b[i][1] - b[i][0];
  454. s->mod[lev][i] = b[i][0] & 1;
  455. for (j = 0; j < 2; j++)
  456. b[i][j] = (b[i][j] + 1) >> 1;
  457. }
  458. switch (type) {
  459. case FF_DWT97:
  460. s->f_linebuf = av_malloc_array((maxlen + 12), sizeof(*s->f_linebuf));
  461. if (!s->f_linebuf)
  462. return AVERROR(ENOMEM);
  463. break;
  464. case FF_DWT97_INT:
  465. s->i_linebuf = av_malloc_array((maxlen + 12), sizeof(*s->i_linebuf));
  466. if (!s->i_linebuf)
  467. return AVERROR(ENOMEM);
  468. break;
  469. case FF_DWT53:
  470. s->i_linebuf = av_malloc_array((maxlen + 6), sizeof(*s->i_linebuf));
  471. if (!s->i_linebuf)
  472. return AVERROR(ENOMEM);
  473. break;
  474. default:
  475. return -1;
  476. }
  477. return 0;
  478. }
  479. int ff_dwt_encode(DWTContext *s, void *t)
  480. {
  481. switch(s->type){
  482. case FF_DWT97:
  483. dwt_encode97_float(s, t); break;
  484. case FF_DWT97_INT:
  485. dwt_encode97_int(s, t); break;
  486. case FF_DWT53:
  487. dwt_encode53(s, t); break;
  488. default:
  489. return -1;
  490. }
  491. return 0;
  492. }
  493. int ff_dwt_decode(DWTContext *s, void *t)
  494. {
  495. switch (s->type) {
  496. case FF_DWT97:
  497. dwt_decode97_float(s, t);
  498. break;
  499. case FF_DWT97_INT:
  500. dwt_decode97_int(s, t);
  501. break;
  502. case FF_DWT53:
  503. dwt_decode53(s, t);
  504. break;
  505. default:
  506. return -1;
  507. }
  508. return 0;
  509. }
  510. void ff_dwt_destroy(DWTContext *s)
  511. {
  512. av_freep(&s->f_linebuf);
  513. av_freep(&s->i_linebuf);
  514. }
  515. #ifdef TEST
  516. #include "libavutil/lfg.h"
  517. #define MAX_W 256
  518. static int test_dwt(int *array, int *ref, uint16_t border[2][2], int decomp_levels, int type, int max_diff) {
  519. int ret, j;
  520. DWTContext s1={{{0}}}, *s= &s1;
  521. int64_t err2 = 0;
  522. ret = ff_jpeg2000_dwt_init(s, border, decomp_levels, type);
  523. if (ret < 0) {
  524. fprintf(stderr, "ff_jpeg2000_dwt_init failed\n");
  525. return 1;
  526. }
  527. ret = ff_dwt_encode(s, array);
  528. if (ret < 0) {
  529. fprintf(stderr, "ff_dwt_encode failed\n");
  530. return 1;
  531. }
  532. ret = ff_dwt_decode(s, array);
  533. if (ret < 0) {
  534. fprintf(stderr, "ff_dwt_encode failed\n");
  535. return 1;
  536. }
  537. for (j = 0; j<MAX_W * MAX_W; j++) {
  538. if (FFABS(array[j] - ref[j]) > max_diff) {
  539. fprintf(stderr, "missmatch at %d (%d != %d) decomp:%d border %d %d %d %d\n",
  540. j, array[j], ref[j],decomp_levels, border[0][0], border[0][1], border[1][0], border[1][1]);
  541. return 2;
  542. }
  543. err2 += (array[j] - ref[j]) * (array[j] - ref[j]);
  544. array[j] = ref[j];
  545. }
  546. ff_dwt_destroy(s);
  547. printf("%s, decomp:%2d border %3d %3d %3d %3d milli-err2:%9"PRId64"\n",
  548. type == FF_DWT53 ? "5/3i" : "9/7i",
  549. decomp_levels, border[0][0], border[0][1], border[1][0], border[1][1],
  550. 1000*err2 / ((border[0][1] - border[0][0])*(border[1][1] - border[1][0])));
  551. return 0;
  552. }
  553. static int test_dwtf(float *array, float *ref, uint16_t border[2][2], int decomp_levels, float max_diff) {
  554. int ret, j;
  555. DWTContext s1={{{0}}}, *s= &s1;
  556. double err2 = 0;
  557. ret = ff_jpeg2000_dwt_init(s, border, decomp_levels, FF_DWT97);
  558. if (ret < 0) {
  559. fprintf(stderr, "ff_jpeg2000_dwt_init failed\n");
  560. return 1;
  561. }
  562. ret = ff_dwt_encode(s, array);
  563. if (ret < 0) {
  564. fprintf(stderr, "ff_dwt_encode failed\n");
  565. return 1;
  566. }
  567. ret = ff_dwt_decode(s, array);
  568. if (ret < 0) {
  569. fprintf(stderr, "ff_dwt_encode failed\n");
  570. return 1;
  571. }
  572. for (j = 0; j<MAX_W * MAX_W; j++) {
  573. if (FFABS(array[j] - ref[j]) > max_diff) {
  574. fprintf(stderr, "missmatch at %d (%f != %f) decomp:%d border %d %d %d %d\n",
  575. j, array[j], ref[j],decomp_levels, border[0][0], border[0][1], border[1][0], border[1][1]);
  576. return 2;
  577. }
  578. err2 += (array[j] - ref[j]) * (array[j] - ref[j]);
  579. array[j] = ref[j];
  580. }
  581. ff_dwt_destroy(s);
  582. printf("9/7f, decomp:%2d border %3d %3d %3d %3d err2:%20.4f\n",
  583. decomp_levels, border[0][0], border[0][1], border[1][0], border[1][1],
  584. err2 / ((border[0][1] - border[0][0])*(border[1][1] - border[1][0])));
  585. return 0;
  586. }
  587. int main(void) {
  588. int array[MAX_W * MAX_W];
  589. int ref [MAX_W * MAX_W];
  590. float arrayf[MAX_W * MAX_W];
  591. float reff [MAX_W * MAX_W];
  592. AVLFG prng;
  593. int i,j;
  594. uint16_t border[2][2];
  595. int ret, decomp_levels;
  596. av_lfg_init(&prng, 1);
  597. for (i = 0; i<MAX_W * MAX_W; i++)
  598. arrayf[i] = reff[i] = array[i] = ref[i] = av_lfg_get(&prng) % 2048;
  599. for (i = 0; i < 100; i++) {
  600. for (j=0; j<4; j++)
  601. border[j>>1][j&1] = av_lfg_get(&prng) % MAX_W;
  602. if (border[0][0] >= border[0][1] || border[1][0] >= border[1][1])
  603. continue;
  604. decomp_levels = av_lfg_get(&prng) % FF_DWT_MAX_DECLVLS;
  605. ret = test_dwt(array, ref, border, decomp_levels, FF_DWT53, 0);
  606. if (ret)
  607. return ret;
  608. ret = test_dwt(array, ref, border, decomp_levels, FF_DWT97_INT, FFMIN(7+5*decomp_levels, 15+3*decomp_levels));
  609. if (ret)
  610. return ret;
  611. ret = test_dwtf(arrayf, reff, border, decomp_levels, 1.0);
  612. if (ret)
  613. return ret;
  614. }
  615. return 0;
  616. }
  617. #endif