Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | RSS feed

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