JPEG 2000 decoder for DCinema

Based on the 2007 GSoC project from Kamil Nowosad <k.nowosad@students.mimuw.edu.pl>
Updated to current programming standards, style and many more small
fixes by Diego Biurrun <diego@biurrun.de>.
Signed-off-by: Diego Biurrun <diego@biurrun.de>

Changelog

@@ -11,6 +11,7 @@ version 10:
   filtergraph description to be read from a file
 - uniform options syntax across all filters
 - new interlace filter
+- JPEG 2000 decoder
 
 
 version 9:

doc/general.texi

@@ -373,8 +373,8 @@ following image formats are supported:
     @tab Digital Picture Exchange
 @item JPEG         @tab X @tab X
     @tab Progressive JPEG is not supported.
-@item JPEG 2000    @tab E @tab E
-    @tab decoding supported through external library libopenjpeg
+@item JPEG 2000    @tab E @tab X
+    @tab encoding supported through external library libopenjpeg
 @item JPEG-LS      @tab X @tab X
 @item LJPEG        @tab X @tab
     @tab Lossless JPEG

libavcodec/Makefile

@@ -205,6 +205,8 @@ OBJS-$(CONFIG_INDEO4_DECODER)          += indeo4.o ivi_common.o ivi_dsp.o
 OBJS-$(CONFIG_INDEO5_DECODER)          += indeo5.o ivi_common.o ivi_dsp.o
 OBJS-$(CONFIG_INTERPLAY_DPCM_DECODER)  += dpcm.o
 OBJS-$(CONFIG_INTERPLAY_VIDEO_DECODER) += interplayvideo.o
+OBJS-$(CONFIG_JPEG2000_DECODER)        += jpeg2000dec.o jpeg2000.o      \
+                                          jpeg2000dwt.o mqcdec.o mqc.o
 OBJS-$(CONFIG_JPEGLS_DECODER)          += jpeglsdec.o jpegls.o \
                                           mjpegdec.o mjpeg.o
 OBJS-$(CONFIG_JPEGLS_ENCODER)          += jpeglsenc.o jpegls.o

libavcodec/allcodecs.c

@@ -160,6 +160,7 @@ void avcodec_register_all(void)
     REGISTER_DECODER(INDEO4,            indeo4);
     REGISTER_DECODER(INDEO5,            indeo5);
     REGISTER_DECODER(INTERPLAY_VIDEO,   interplay_video);
+    REGISTER_DECODER(JPEG2000,          jpeg2000);
     REGISTER_ENCDEC (JPEGLS,            jpegls);
     REGISTER_DECODER(JV,                jv);
     REGISTER_DECODER(KGV1,              kgv1);

libavcodec/avcodec.h

@@ -2570,6 +2570,12 @@ typedef struct AVCodecContext {
 #define FF_PROFILE_MPEG4_SIMPLE_STUDIO             14
 #define FF_PROFILE_MPEG4_ADVANCED_SIMPLE           15
 
+#define FF_PROFILE_JPEG2000_CSTREAM_RESTRICTION_0   0
+#define FF_PROFILE_JPEG2000_CSTREAM_RESTRICTION_1   1
+#define FF_PROFILE_JPEG2000_CSTREAM_NO_RESTRICTION  2
+#define FF_PROFILE_JPEG2000_DCINEMA_2K              3
+#define FF_PROFILE_JPEG2000_DCINEMA_4K              4
+
     /**
      * level
      * - encoding: Set by user.

libavcodec/jpeg2000.c

+/*
+ * JPEG 2000 encoder and decoder common functions
+ * Copyright (c) 2007 Kamil Nowosad
+ * Copyright (c) 2013 Nicolas Bertrand <nicoinattendu@gmail.com>
+ *
+ * This file is part of Libav.
+ *
+ * Libav is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * Libav is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with Libav; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+/**
+ * @file
+ * JPEG 2000 image encoder and decoder common functions
+ */
+
+#include "libavutil/common.h"
+#include "libavutil/mem.h"
+#include "avcodec.h"
+#include "jpeg2000.h"
+
+#define SHL(a, n) ((n) >= 0 ? (a) << (n) : (a) >> -(n))
+
+/* tag tree routines */
+
+/* allocate the memory for tag tree */
+static int32_t tag_tree_size(uint16_t w, uint16_t h)
+{
+    uint32_t res = 0;
+    while (w > 1 || h > 1) {
+        res += w * h;
+        if (res + 1 >= INT32_MAX)
+            return -1;
+        w = (w + 1) >> 1;
+        h = (h + 1) >> 1;
+    }
+    return (int32_t)(res + 1);
+}
+
+static Jpeg2000TgtNode *ff_jpeg2000_tag_tree_init(int w, int h)
+{
+    int pw = w, ph = h;
+    Jpeg2000TgtNode *res, *t, *t2;
+    int32_t tt_size;
+
+    tt_size = tag_tree_size(w, h);
+    if (tt_size == -1)
+        return NULL;
+
+    t = res = av_mallocz_array(tt_size, sizeof(*t));
+    if (!res)
+        return NULL;
+
+    while (w > 1 || h > 1) {
+        int i, j;
+        pw = w;
+        ph = h;
+
+        w  = (w + 1) >> 1;
+        h  = (h + 1) >> 1;
+        t2 = t + pw * ph;
+
+        for (i = 0; i < ph; i++)
+            for (j = 0; j < pw; j++)
+                t[i * pw + j].parent = &t2[(i >> 1) * w + (j >> 1)];
+
+        t = t2;
+    }
+    t[0].parent = NULL;
+    return res;
+}
+
+uint8_t ff_jpeg2000_sigctxno_lut[256][4];
+
+static int getsigctxno(int flag, int bandno)
+{
+    int h, v, d;
+
+    h = ((flag & JPEG2000_T1_SIG_E)  ? 1 : 0) +
+        ((flag & JPEG2000_T1_SIG_W)  ? 1 : 0);
+    v = ((flag & JPEG2000_T1_SIG_N)  ? 1 : 0) +
+        ((flag & JPEG2000_T1_SIG_S)  ? 1 : 0);
+    d = ((flag & JPEG2000_T1_SIG_NE) ? 1 : 0) +
+        ((flag & JPEG2000_T1_SIG_NW) ? 1 : 0) +
+        ((flag & JPEG2000_T1_SIG_SE) ? 1 : 0) +
+        ((flag & JPEG2000_T1_SIG_SW) ? 1 : 0);
+    if (bandno < 3) {
+        if (bandno == 1)
+            FFSWAP(int, h, v);
+        if (h == 2)
+            return 8;
+        if (h == 1) {
+            if (v >= 1)
+                return 7;
+            if (d >= 1)
+                return 6;
+            return 5;
+        }
+        if (v == 2)
+            return 4;
+        if (v == 1)
+            return 3;
+        if (d >= 2)
+            return 2;
+        if (d == 1)
+            return 1;
+        return 0;
+    } else {
+        if (d >= 3)
+            return 8;
+        if (d == 2) {
+            if (h + v >= 1)
+                return 7;
+            return 6;
+        }
+        if (d == 1) {
+            if (h + v >= 2)
+                return 5;
+            if (h + v == 1)
+                return 4;
+            return 3;
+        }
+        if (h + v >= 2)
+            return 2;
+        if (h + v == 1)
+            return 1;
+        return 0;
+    }
+    return 0;
+}
+
+uint8_t ff_jpeg2000_sgnctxno_lut[16][16], ff_jpeg2000_xorbit_lut[16][16];
+
+static const int contribtab[3][3] = { {  0, -1,  1 }, { -1, -1,  0 }, {  1,  0,  1 } };
+static const int  ctxlbltab[3][3] = { { 13, 12, 11 }, { 10,  9, 10 }, { 11, 12, 13 } };
+static const int  xorbittab[3][3] = { {  1,  1,  1 }, {  1,  0,  0 }, {  0,  0,  0 } };
+
+static int getsgnctxno(int flag, uint8_t *xorbit)
+{
+    int vcontrib, hcontrib;
+
+    hcontrib = contribtab[flag & JPEG2000_T1_SIG_E ? flag & JPEG2000_T1_SGN_E ? 1 : 2 : 0]
+                         [flag & JPEG2000_T1_SIG_W ? flag & JPEG2000_T1_SGN_W ? 1 : 2 : 0] + 1;
+    vcontrib = contribtab[flag & JPEG2000_T1_SIG_S ? flag & JPEG2000_T1_SGN_S ? 1 : 2 : 0]
+                         [flag & JPEG2000_T1_SIG_N ? flag & JPEG2000_T1_SGN_N ? 1 : 2 : 0] + 1;
+    *xorbit = xorbittab[hcontrib][vcontrib];
+
+    return ctxlbltab[hcontrib][vcontrib];
+}
+
+void ff_jpeg2000_init_tier1_luts(void)
+{
+    int i, j;
+    for (i = 0; i < 256; i++)
+        for (j = 0; j < 4; j++)
+            ff_jpeg2000_sigctxno_lut[i][j] = getsigctxno(i, j);
+    for (i = 0; i < 16; i++)
+        for (j = 0; j < 16; j++)
+            ff_jpeg2000_sgnctxno_lut[i][j] =
+                getsgnctxno(i + (j << 8), &ff_jpeg2000_xorbit_lut[i][j]);
+}
+
+void ff_jpeg2000_set_significance(Jpeg2000T1Context *t1, int x, int y,
+                                  int negative)
+{
+    x++;
+    y++;
+    t1->flags[y][x] |= JPEG2000_T1_SIG;
+    if (negative) {
+        t1->flags[y][x + 1] |= JPEG2000_T1_SIG_W | JPEG2000_T1_SGN_W;
+        t1->flags[y][x - 1] |= JPEG2000_T1_SIG_E | JPEG2000_T1_SGN_E;
+        t1->flags[y + 1][x] |= JPEG2000_T1_SIG_N | JPEG2000_T1_SGN_N;
+        t1->flags[y - 1][x] |= JPEG2000_T1_SIG_S | JPEG2000_T1_SGN_S;
+    } else {
+        t1->flags[y][x + 1] |= JPEG2000_T1_SIG_W;
+        t1->flags[y][x - 1] |= JPEG2000_T1_SIG_E;
+        t1->flags[y + 1][x] |= JPEG2000_T1_SIG_N;
+        t1->flags[y - 1][x] |= JPEG2000_T1_SIG_S;
+    }
+    t1->flags[y + 1][x + 1] |= JPEG2000_T1_SIG_NW;
+    t1->flags[y + 1][x - 1] |= JPEG2000_T1_SIG_NE;
+    t1->flags[y - 1][x + 1] |= JPEG2000_T1_SIG_SW;
+    t1->flags[y - 1][x - 1] |= JPEG2000_T1_SIG_SE;
+}
+
+static const uint8_t lut_gain[2][4] = { { 0, 0, 0, 0 }, { 0, 1, 1, 2 } };
+
+int ff_jpeg2000_init_component(Jpeg2000Component *comp,
+                               Jpeg2000CodingStyle *codsty,
+                               Jpeg2000QuantStyle *qntsty,
+                               int cbps, int dx, int dy,
+                               AVCodecContext *avctx)
+{
+    uint8_t log2_band_prec_width, log2_band_prec_height;
+    int reslevelno, bandno, gbandno = 0, ret, i, j;
+    uint32_t csize = 1;
+
+    if (ret = ff_jpeg2000_dwt_init(&comp->dwt, comp->coord,
+                                   codsty->nreslevels2decode - 1,
+                                   codsty->transform))
+        return ret;
+    // component size comp->coord is uint16_t so ir cannot overflow
+    csize = (comp->coord[0][1] - comp->coord[0][0]) *
+            (comp->coord[1][1] - comp->coord[1][0]);
+
+    comp->data = av_malloc_array(csize, sizeof(*comp->data));
+    if (!comp->data)
+        return AVERROR(ENOMEM);
+    comp->reslevel = av_malloc_array(codsty->nreslevels, sizeof(*comp->reslevel));
+    if (!comp->reslevel)
+        return AVERROR(ENOMEM);
+    /* LOOP on resolution levels */
+    for (reslevelno = 0; reslevelno < codsty->nreslevels; reslevelno++) {
+        int declvl = codsty->nreslevels - reslevelno;    // N_L -r see  ISO/IEC 15444-1:2002 B.5
+        Jpeg2000ResLevel *reslevel = comp->reslevel + reslevelno;
+
+        /* Compute borders for each resolution level.
+         * Computation of trx_0, trx_1, try_0 and try_1.
+         * see ISO/IEC 15444-1:2002 eq. B.5 and B-14 */
+        for (i = 0; i < 2; i++)
+            for (j = 0; j < 2; j++)
+                reslevel->coord[i][j] =
+                    ff_jpeg2000_ceildivpow2(comp->coord_o[i][j], declvl - 1);
+        // update precincts size: 2^n value
+        reslevel->log2_prec_width  = codsty->log2_prec_widths[reslevelno];
+        reslevel->log2_prec_height = codsty->log2_prec_heights[reslevelno];
+
+        /* Number of bands for each resolution level */
+        if (reslevelno == 0)
+            reslevel->nbands = 1;
+        else
+            reslevel->nbands = 3;
+
+        /* Number of precincts wich span the tile for resolution level reslevelno
+         * see B.6 in ISO/IEC 15444-1:2002 eq. B-16
+         * num_precincts_x = |- trx_1 / 2 ^ log2_prec_width) -| - (trx_0 / 2 ^ log2_prec_width)
+         * num_precincts_y = |- try_1 / 2 ^ log2_prec_width) -| - (try_0 / 2 ^ log2_prec_width)
+         * for Dcinema profiles in JPEG 2000
+         * num_precincts_x = |- trx_1 / 2 ^ log2_prec_width) -|
+         * num_precincts_y = |- try_1 / 2 ^ log2_prec_width) -| */
+        if (reslevel->coord[0][1] == reslevel->coord[0][0])
+            reslevel->num_precincts_x = 0;
+        else
+            reslevel->num_precincts_x =
+                ff_jpeg2000_ceildivpow2(reslevel->coord[0][1],
+                                        reslevel->log2_prec_width) -
+                (reslevel->coord[0][0] >> reslevel->log2_prec_width);
+
+        if (reslevel->coord[1][1] == reslevel->coord[1][0])
+            reslevel->num_precincts_y = 0;
+        else
+            reslevel->num_precincts_y =
+                ff_jpeg2000_ceildivpow2(reslevel->coord[1][1],
+                                        reslevel->log2_prec_height) -
+                (reslevel->coord[1][0] >> reslevel->log2_prec_height);
+
+        reslevel->band = av_malloc_array(reslevel->nbands, sizeof(*reslevel->band));
+        if (!reslevel->band)
+            return AVERROR(ENOMEM);
+
+        for (bandno = 0; bandno < reslevel->nbands; bandno++, gbandno++) {
+            Jpeg2000Band *band = reslevel->band + bandno;
+            int cblkno, precno;
+            int nb_precincts;
+
+            /* TODO: Implementation of quantization step not finished,
+             * see ISO/IEC 15444-1:2002 E.1 and A.6.4. */
+            switch (qntsty->quantsty) {
+                uint8_t gain;
+                int numbps;
+            case JPEG2000_QSTY_NONE:
+                /* TODO: to verify. No quantization in this case */
+                numbps = cbps +
+                         lut_gain[codsty->transform][bandno + reslevelno > 0];
+                band->stepsize = (float)SHL(2048 + qntsty->mant[gbandno],
+                                            2 + numbps - qntsty->expn[gbandno]);
+                break;
+            case JPEG2000_QSTY_SI:
+                /*TODO: Compute formula to implement. */
+                band->stepsize = (float) (1 << 13);
+                break;
+            case JPEG2000_QSTY_SE:
+                /* Exponent quantization step.
+                 * Formula:
+                 * delta_b = 2 ^ (R_b - expn_b) * (1 + (mant_b / 2 ^ 11))
+                 * R_b = R_I + log2 (gain_b )
+                 * see ISO/IEC 15444-1:2002 E.1.1 eqn. E-3 and E-4 */
+                /* TODO/WARN: value of log2 (gain_b ) not taken into account
+                 * but it works (compared to OpenJPEG). Why?
+                 * Further investigation needed. */
+                gain            = cbps;
+                band->stepsize  = pow(2.0, gain - qntsty->expn[gbandno]);
+                band->stepsize *= (float)qntsty->mant[gbandno] / 2048.0 + 1.0;
+                /* FIXME: In openjepg code stespize = stepsize * 0.5. Why?
+                 * If not set output of entropic decoder is not correct. */
+                band->stepsize *= 0.5;
+                break;
+            default:
+                band->stepsize = 0;
+                av_log(avctx, AV_LOG_ERROR, "Unknown quantization format\n");
+                break;
+            }
+            /* BITEXACT computing case --> convert to int */
+            if (avctx->flags & CODEC_FLAG_BITEXACT)
+                band->stepsize = (int32_t)(band->stepsize * (1 << 16));
+
+            /* computation of tbx_0, tbx_1, tby_0, tby_1
+             * see ISO/IEC 15444-1:2002 B.5 eq. B-15 and tbl B.1
+             * codeblock width and height is computed for
+             * DCI JPEG 2000 codeblock_width = codeblock_width = 32 = 2 ^ 5 */
+            if (reslevelno == 0) {
+                /* for reslevelno = 0, only one band, x0_b = y0_b = 0 */
+                for (i = 0; i < 2; i++)
+                    for (j = 0; j < 2; j++)
+                        band->coord[i][j] =
+                            ff_jpeg2000_ceildivpow2(comp->coord_o[i][j],
+                                                    declvl - 1);
+
+                log2_band_prec_width  = reslevel->log2_prec_width;
+                log2_band_prec_height = reslevel->log2_prec_height;
+                /* see ISO/IEC 15444-1:2002 eq. B-17 and eq. B-15 */
+                band->log2_cblk_width  = FFMIN(codsty->log2_cblk_width,
+                                               reslevel->log2_prec_width);
+                band->log2_cblk_height = FFMIN(codsty->log2_cblk_height,
+                                               reslevel->log2_prec_height);
+            } else {
+                /* 3 bands x0_b = 1 y0_b = 0; x0_b = 0 y0_b = 1; x0_b = y0_b = 1 */
+                /* x0_b and y0_b are computed with ((bandno + 1 >> i) & 1) */
+                for (i = 0; i < 2; i++)
+                    for (j = 0; j < 2; j++)
+                        /* Formula example for tbx_0 = ceildiv((tcx_0 - 2 ^ (declvl - 1) * x0_b) / declvl) */
+                        band->coord[i][j] =
+                            ff_jpeg2000_ceildivpow2(comp->coord_o[i][j] -
+                                                    (((bandno + 1 >> i) & 1) << declvl - 1),
+                                                    declvl);
+                /* TODO: Manage case of 3 band offsets here or
+                 * in coding/decoding function? */
+
+                /* see ISO/IEC 15444-1:2002 eq. B-17 and eq. B-15 */
+                band->log2_cblk_width  = FFMIN(codsty->log2_cblk_width,
+                                               reslevel->log2_prec_width - 1);
+                band->log2_cblk_height = FFMIN(codsty->log2_cblk_height,
+                                               reslevel->log2_prec_height - 1);
+
+                log2_band_prec_width  = reslevel->log2_prec_width  - 1;
+                log2_band_prec_height = reslevel->log2_prec_height - 1;
+            }
+
+            band->prec = av_malloc_array(reslevel->num_precincts_x *
+                                         reslevel->num_precincts_y,
+                                         sizeof(*band->prec));
+            if (!band->prec)
+                return AVERROR(ENOMEM);
+
+            nb_precincts = reslevel->num_precincts_x * reslevel->num_precincts_y;
+
+            for (precno = 0; precno < nb_precincts; precno++) {
+                Jpeg2000Prec *prec = band->prec + precno;
+
+                /* TODO: Explain formula for JPEG200 DCINEMA. */
+                /* TODO: Verify with previous count of codeblocks per band */
+
+                /* Compute P_x0 */
+                prec->coord[0][0] = (precno % reslevel->num_precincts_x) *
+                                    (1 << log2_band_prec_width);
+                prec->coord[0][0] = FFMAX(prec->coord[0][0], band->coord[0][0]);
+
+                /* Compute P_y0 */
+                prec->coord[1][0] = (precno / reslevel->num_precincts_x) *
+                                    (1 << log2_band_prec_height);
+                prec->coord[1][0] = FFMAX(prec->coord[1][0], band->coord[1][0]);
+
+                /* Compute P_x1 */
+                prec->coord[0][1] = prec->coord[0][0] +
+                                    (1 << log2_band_prec_width);
+                prec->coord[0][1] = FFMIN(prec->coord[0][1], band->coord[0][1]);
+
+                /* Compute P_y1 */
+                prec->coord[1][1] = prec->coord[1][0] +
+                                    (1 << log2_band_prec_height);
+                prec->coord[1][1] = FFMIN(prec->coord[1][1], band->coord[1][1]);
+
+                prec->nb_codeblocks_width =
+                    ff_jpeg2000_ceildivpow2(prec->coord[0][1] -
+                                            prec->coord[0][0],
+                                            band->log2_cblk_width);
+                prec->nb_codeblocks_height =
+                    ff_jpeg2000_ceildivpow2(prec->coord[1][1] -
+                                            prec->coord[1][0],
+                                            band->log2_cblk_height);
+
+                /* Tag trees initialization */
+                prec->cblkincl =
+                    ff_jpeg2000_tag_tree_init(prec->nb_codeblocks_width,
+                                              prec->nb_codeblocks_height);
+                if (!prec->cblkincl)
+                    return AVERROR(ENOMEM);
+
+                prec->zerobits =
+                    ff_jpeg2000_tag_tree_init(prec->nb_codeblocks_width,
+                                              prec->nb_codeblocks_height);
+                if (!prec->zerobits)
+                    return AVERROR(ENOMEM);
+
+                prec->cblk = av_malloc_array(prec->nb_codeblocks_width *
+                                             prec->nb_codeblocks_height,
+                                             sizeof(*prec->cblk));
+                if (!prec->cblk)
+                    return AVERROR(ENOMEM);
+                for (cblkno = 0; cblkno < prec->nb_codeblocks_width * prec->nb_codeblocks_height; cblkno++) {
+                    Jpeg2000Cblk *cblk = prec->cblk + cblkno;
+                    uint16_t Cx0, Cy0;
+
+                    /* Compute coordinates of codeblocks */
+                    /* Compute Cx0*/
+                    Cx0 = (prec->coord[0][0] >> band->log2_cblk_width) << band->log2_cblk_width;
+                    Cx0 = Cx0 + ((cblkno % prec->nb_codeblocks_width)  << band->log2_cblk_width);
+                    cblk->coord[0][0] = FFMAX(Cx0, prec->coord[0][0]);
+
+                    /* Compute Cy0*/
+                    Cy0 = (prec->coord[1][0] >> band->log2_cblk_height) << band->log2_cblk_height;
+                    Cy0 = Cy0 + ((cblkno / prec->nb_codeblocks_width)   << band->log2_cblk_height);
+                    cblk->coord[1][0] = FFMAX(Cy0, prec->coord[1][0]);
+
+                    /* Compute Cx1 */
+                    cblk->coord[0][1] = FFMIN(Cx0 + (1 << band->log2_cblk_width),
+                                              prec->coord[0][1]);
+
+                    /* Compute Cy1 */
+                    cblk->coord[1][1] = FFMIN(Cy0 + (1 << band->log2_cblk_height),
+                                              prec->coord[1][1]);
+                    cblk->zero      = 0;
+                    cblk->lblock    = 3;
+                    cblk->length    = 0;
+                    cblk->lengthinc = 0;
+                    cblk->npasses   = 0;
+                }
+            }
+        }
+    }
+    return 0;
+}
+
+void ff_jpeg2000_cleanup(Jpeg2000Component *comp, Jpeg2000CodingStyle *codsty)
+{
+    int reslevelno, bandno, precno;
+    for (reslevelno = 0; reslevelno < codsty->nreslevels; reslevelno++) {
+        Jpeg2000ResLevel *reslevel = comp->reslevel + reslevelno;
+
+        for (bandno = 0; bandno < reslevel->nbands; bandno++) {
+            Jpeg2000Band *band = reslevel->band + bandno;
+            for (precno = 0; precno < reslevel->num_precincts_x * reslevel->num_precincts_y; precno++) {
+                Jpeg2000Prec *prec = band->prec + precno;
+                av_freep(&prec->zerobits);
+                av_freep(&prec->cblkincl);
+                av_freep(&prec->cblk);
+            }
+
+            av_freep(&band->prec);
+        }
+        av_freep(&reslevel->band);
+    }
+
+    ff_dwt_destroy(&comp->dwt);
+    av_freep(&comp->reslevel);
+    av_freep(&comp->data);
+}

libavcodec/jpeg2000.h

+/*
+ * JPEG 2000 common defines, structures and functions
+ * Copyright (c) 2007 Kamil Nowosad
+ * Copyright (c) 2013 Nicolas Bertrand <nicoinattendu@gmail.com>
+ *
+ * This file is part of Libav.
+ *
+ * Libav is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * Libav is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with Libav; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#ifndef AVCODEC_JPEG2000_H
+#define AVCODEC_JPEG2000_H
+
+/**
+ * @file
+ * JPEG 2000 structures and defines common
+ * to encoder and decoder
+ */
+
+#include <stdint.h>
+
+#include "avcodec.h"
+#include "mqc.h"
+#include "jpeg2000dwt.h"
+
+enum Jpeg2000Markers {
+    JPEG2000_SOC = 0xff4f, // start of codestream
+    JPEG2000_SIZ = 0xff51, // image and tile size
+    JPEG2000_COD,          // coding style default
+    JPEG2000_COC,          // coding style component
+    JPEG2000_TLM = 0xff55, // packed packet headers, tile-part header
+    JPEG2000_PLM = 0xff57, // tile-part lengths
+    JPEG2000_PLT,          // packet length, main header
+    JPEG2000_QCD = 0xff5c, // quantization default
+    JPEG2000_QCC,          // quantization component
+    JPEG2000_RGN,          // region of interest
+    JPEG2000_POC,          // progression order change
+    JPEG2000_PPM,          // packet length, tile-part header
+    JPEG2000_PPT,          // packed packet headers, main header
+    JPEG2000_CRG = 0xff63, // component registration
+    JPEG2000_COM,          // comment
+    JPEG2000_SOT = 0xff90, // start of tile-part
+    JPEG2000_SOP,          // start of packet
+    JPEG2000_EPH,          // end of packet header
+    JPEG2000_SOD,          // start of data
+    JPEG2000_EOC = 0xffd9, // end of codestream
+};
+
+enum Jpeg2000Quantsty { // quantization style
+    JPEG2000_QSTY_NONE, // no quantization
+    JPEG2000_QSTY_SI,   // scalar derived
+    JPEG2000_QSTY_SE    // scalar expounded
+};
+
+#define JPEG2000_MAX_CBLKW 64
+#define JPEG2000_MAX_CBLKH 64
+
+#define JPEG2000_MAX_RESLEVELS 33
+
+// T1 flags
+// flags determining significance of neighbor coefficients
+#define JPEG2000_T1_SIG_N  0x0001
+#define JPEG2000_T1_SIG_E  0x0002
+#define JPEG2000_T1_SIG_W  0x0004
+#define JPEG2000_T1_SIG_S  0x0008
+#define JPEG2000_T1_SIG_NE 0x0010
+#define JPEG2000_T1_SIG_NW 0x0020
+#define JPEG2000_T1_SIG_SE 0x0040
+#define JPEG2000_T1_SIG_SW 0x0080
+#define JPEG2000_T1_SIG_NB (JPEG2000_T1_SIG_N  | JPEG2000_T1_SIG_E  |   \
+                            JPEG2000_T1_SIG_S  | JPEG2000_T1_SIG_W  |   \
+                            JPEG2000_T1_SIG_NE | JPEG2000_T1_SIG_NW |   \
+                            JPEG2000_T1_SIG_SE | JPEG2000_T1_SIG_SW)
+// flags determining sign bit of neighbor coefficients
+#define JPEG2000_T1_SGN_N  0x0100
+#define JPEG2000_T1_SGN_S  0x0200
+#define JPEG2000_T1_SGN_W  0x0400
+#define JPEG2000_T1_SGN_E  0x0800
+
+#define JPEG2000_T1_VIS    0x1000
+#define JPEG2000_T1_SIG    0x2000
+#define JPEG2000_T1_REF    0x4000
+
+#define JPEG2000_T1_SGN    0x8000
+
+// Codeblock coding styles
+#define JPEG2000_CBLK_BYPASS    0x01 // Selective arithmetic coding bypass
+#define JPEG2000_CBLK_RESET     0x02 // Reset context probabilities
+#define JPEG2000_CBLK_TERMALL   0x04 // Terminate after each coding pass
+#define JPEG2000_CBLK_VSC       0x08 // Vertical stripe causal context formation
+#define JPEG2000_CBLK_PREDTERM  0x10 // Predictable termination
+#define JPEG2000_CBLK_SEGSYM    0x20 // Segmentation symbols present
+
+// Coding styles
+#define JPEG2000_CSTY_PREC      0x01 // Precincts defined in coding style
+#define JPEG2000_CSTY_SOP       0x02 // SOP marker present
+#define JPEG2000_CSTY_EPH       0x04 // EPH marker present
+
+// Progression orders
+#define JPEG2000_PGOD_LRCP      0x00  // Layer-resolution level-component-position progression
+#define JPEG2000_PGOD_RLCP      0x01  // Resolution level-layer-component-position progression
+#define JPEG2000_PGOD_RPCL      0x02  // Resolution level-position-component-layer progression
+#define JPEG2000_PGOD_PCRL      0x03  // Position-component-resolution level-layer progression
+#define JPEG2000_PGOD_CPRL      0x04  // Component-position-resolution level-layer progression
+
+typedef struct Jpeg2000T1Context {
+    int data[JPEG2000_MAX_CBLKW][JPEG2000_MAX_CBLKH];
+    int flags[JPEG2000_MAX_CBLKW + 2][JPEG2000_MAX_CBLKH + 2];
+    MqcState mqc;
+} Jpeg2000T1Context;
+
+typedef struct Jpeg2000TgtNode {
+    uint8_t val;
+    uint8_t vis;
+    struct Jpeg2000TgtNode *parent;
+} Jpeg2000TgtNode;
+
+typedef struct Jpeg2000CodingStyle {
+    uint8_t nreslevels;       // number of resolution levels
+    uint8_t nreslevels2decode; // number of resolution levels to decode
+    uint8_t log2_cblk_width,
+            log2_cblk_height; // exponent of codeblock size
+    uint8_t transform;        // DWT type
+    uint8_t csty;             // coding style
+    uint8_t log2_prec_width,
+            log2_prec_height; // precinct size
+    uint8_t nlayers;          // number of layers
+    uint8_t mct;              // multiple component transformation
+    uint8_t cblk_style;       // codeblock coding style
+    uint8_t prog_order;       // progression order
+    uint8_t log2_prec_widths[JPEG2000_MAX_RESLEVELS];  // precincts size according resolution levels
+    uint8_t log2_prec_heights[JPEG2000_MAX_RESLEVELS]; // TODO: initialize prec_size array with 0?
+} Jpeg2000CodingStyle;
+
+typedef struct Jpeg2000QuantStyle {
+    uint8_t expn[32 * 3];  // quantization exponent
+    uint32_t mant[32 * 3]; // quantization mantissa
+    uint8_t quantsty;      // quantization style
+    uint8_t nguardbits;    // number of guard bits
+} Jpeg2000QuantStyle;
+
+typedef struct Jpeg2000Pass {
+    uint16_t rate;
+    int64_t disto;
+} Jpeg2000Pass;
+
+typedef struct Jpeg2000Cblk {
+    uint8_t npasses;
+    uint8_t ninclpasses; // number coding of passes included in codestream
+    uint8_t nonzerobits;
+    uint16_t length;
+    uint16_t lengthinc;
+    uint8_t lblock;
+    uint8_t zero;
+    uint8_t data[8192];
+    Jpeg2000Pass passes[100];
+    uint16_t coord[2][2]; // border coordinates {{x0, x1}, {y0, y1}}
+} Jpeg2000Cblk; // code block
+
+typedef struct Jpeg2000Prec {
+    uint16_t xi0, yi0; // codeblock indexes ([xi0, xi1))
+    uint16_t nb_codeblocks_width;
+    uint16_t nb_codeblocks_height;
+    Jpeg2000TgtNode *zerobits;
+    Jpeg2000TgtNode *cblkincl;
+    Jpeg2000Cblk *cblk;
+    uint16_t coord[2][2]; // border coordinates {{x0, x1}, {y0, y1}}
+} Jpeg2000Prec; // precinct
+
+/* TODO: stepsize can be float or integer depending on
+ * reversible or irreversible transformation. */
+typedef struct Jpeg2000Band {
+    uint16_t coord[2][2]; // border coordinates {{x0, x1}, {y0, y1}}
+    uint16_t log2_cblk_width, log2_cblk_height;
+    uint16_t cblknx, cblkny;
+    float stepsize; // quantization stepsize
+    Jpeg2000Prec *prec;
+} Jpeg2000Band; // subband
+
+typedef struct Jpeg2000ResLevel {
+    uint8_t nbands;
+    uint16_t coord[2][2]; // border coordinates {{x0, x1}, {y0, y1}}
+    uint16_t num_precincts_x, num_precincts_y; // number of precincts in x/y direction
+    uint8_t log2_prec_width, log2_prec_height; // exponent of precinct size
+    Jpeg2000Band *band;
+} Jpeg2000ResLevel; // resolution level
+
+/* TODO: data can be float of integer depending of reversible/irreversible
+ * transformation.
+ */
+typedef struct Jpeg2000Component {
+    Jpeg2000ResLevel *reslevel;
+    DWTContext dwt;
+    float *data;
+    uint16_t coord[2][2];   // border coordinates {{x0, x1}, {y0, y1}} -- can be reduced with lowres option
+    uint16_t coord_o[2][2]; // border coordinates {{x0, x1}, {y0, y1}} -- original values from jpeg2000 headers
+} Jpeg2000Component;
+
+/* misc tools */
+static inline int ff_jpeg2000_ceildivpow2(int a, int b)
+{
+    return (a + (1 << b) - 1) >> b;
+}
+
+static inline int ff_jpeg2000_ceildiv(int a, int b)
+{
+    return (a + b - 1) / b;
+}
+
+/* TIER-1 routines */
+
+/* Set up lookup tables used in TIER-1. */
+void ff_jpeg2000_init_tier1_luts(void);
+
+/* Update significance of a coefficient at current position (x,y) and
+ * for neighbors. */
+void ff_jpeg2000_set_significance(Jpeg2000T1Context *t1,
+                                  int x, int y, int negative);
+
+extern uint8_t ff_jpeg2000_sigctxno_lut[256][4];
+
+/* Get context label (number in range[0..8]) of a coefficient for significance
+ * propagation and cleanup coding passes. */
+static inline int ff_jpeg2000_getsigctxno(int flag, int bandno)
+{
+    return ff_jpeg2000_sigctxno_lut[flag & 255][bandno];
+}
+
+static const uint8_t refctxno_lut[2][2] = { { 14, 15 }, { 16, 16 } };
+
+/* Get context label (number in range[14..16]) of a coefficient for magnitude
+ * refinement pass. */
+static inline int ff_jpeg2000_getrefctxno(int flag)
+{
+    return refctxno_lut[(flag >> 14) & 1][(flag & 255) != 0];
+}
+
+extern uint8_t ff_jpeg2000_sgnctxno_lut[16][16];
+extern uint8_t ff_jpeg2000_xorbit_lut[16][16];
+
+/* Get context label (number in range[9..13]) for sign decoding. */
+static inline int ff_jpeg2000_getsgnctxno(int flag, int *xorbit)
+{
+    *xorbit = ff_jpeg2000_xorbit_lut[flag & 15][(flag >> 8) & 15];
+    return ff_jpeg2000_sgnctxno_lut[flag & 15][(flag >> 8) & 15];
+}
+
+int ff_jpeg2000_init_component(Jpeg2000Component *comp,
+                               Jpeg2000CodingStyle *codsty,
+                               Jpeg2000QuantStyle *qntsty,
+                               int cbps, int dx, int dy,
+                               AVCodecContext *ctx);
+
+void ff_jpeg2000_cleanup(Jpeg2000Component *comp, Jpeg2000CodingStyle *codsty);
+
+#endif /* AVCODEC_JPEG2000_H */

libavcodec/jpeg2000dec.c

+/*
+ * JPEG 2000 image decoder
+ * Copyright (c) 2007 Kamil Nowosad
+ * Copyright (c) 2013 Nicolas Bertrand <nicoinattendu@gmail.com>
+ *
+ * This file is part of Libav.
+ *
+ * Libav is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * Libav is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with Libav; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+/**
+ * @file
+ * JPEG 2000 image decoder
+ */
+
+#include "libavutil/common.h"
+#include "libavutil/opt.h"
+#include "avcodec.h"
+#include "bytestream.h"
+#include "internal.h"
+#include "jpeg2000.h"
+
+#define JP2_SIG_TYPE    0x6A502020
+#define JP2_SIG_VALUE   0x0D0A870A
+#define JP2_CODESTREAM  0x6A703263
+
+#define HAD_COC 0x01
+#define HAD_QCC 0x02
+
+typedef struct Jpeg2000TilePart {
+    uint16_t tp_idx;                    // Tile-part index
+    uint8_t tile_index;                 // Tile index who refers the tile-part
+    uint32_t tp_len;                    // Length of tile-part
+    const uint8_t *tp_start_bstrm;      // Start address bit stream in tile-part
+    const uint8_t *tp_end_bstrm;        // End address of the bit stream tile part
+} Jpeg2000TilePart;
+
+/* RMK: For JPEG2000 DCINEMA 3 tile-parts in a tile
+ * one per component, so tile_part elements have a size of 3 */
+typedef struct Jpeg2000Tile {
+    Jpeg2000Component   *comp;
+    uint8_t             properties[4];
+    Jpeg2000CodingStyle codsty[4];
+    Jpeg2000QuantStyle  qntsty[4];
+    Jpeg2000TilePart    tile_part[3];
+} Jpeg2000Tile;
+
+typedef struct Jpeg2000DecoderContext {
+    AVClass         *class;
+    AVCodecContext  *avctx;
+
+    int             width, height;
+    int             image_offset_x, image_offset_y;
+    int             tile_offset_x, tile_offset_y;
+    uint8_t         cbps[4];    // bits per sample in particular components
+    uint8_t         sgnd[4];    // if a component is signed
+    uint8_t         properties[4];
+    int             cdx[4], cdy[4];
+    int             precision;
+    int             ncomponents;
+    int             tile_width, tile_height;
+    int             numXtiles, numYtiles;
+    int             maxtilelen;
+
+    Jpeg2000CodingStyle codsty[4];
+    Jpeg2000QuantStyle  qntsty[4];
+
+    const uint8_t   *buf_start;
+    const uint8_t   *buf;
+    const uint8_t   *buf_end;
+    int             bit_index;
+
+    int16_t         curtileno;
+    Jpeg2000Tile    *tile;
+
+    /*options parameters*/
+    int16_t         lowres;
+    int16_t         reduction_factor;
+} Jpeg2000DecoderContext;
+
+/* get_bits functions for JPEG2000 packet bitstream
+ * It is a get_bit function with a bit-stuffing routine. If the value of the
+ * byte is 0xFF, the next byte includes an extra zero bit stuffed into the MSB.
+ * cf. ISO-15444-1:2002 / B.10.1 Bit-stuffing routine */
+static int get_bits(Jpeg2000DecoderContext *s, int n)
+{
+    int res = 0;
+    if (s->buf_end - s->buf < ((n - s->bit_index) >> 8))
+        return AVERROR(EINVAL);
+    while (--n >= 0) {
+        res <<= 1;
+        if (s->bit_index == 0) {
+            s->bit_index = 7 + (*s->buf != 0xff);
+            s->buf++;
+        }
+        s->bit_index--;
+        res |= (*s->buf >> s->bit_index) & 1;
+    }
+    return res;
+}
+
+static void jpeg2000_flush(Jpeg2000DecoderContext *s)
+{
+    if (*s->buf == 0xff)
+        s->buf++;
+    s->bit_index = 8;
+    s->buf++;
+}
+
+/* decode the value stored in node */
+static int tag_tree_decode(Jpeg2000DecoderContext *s, Jpeg2000TgtNode *node,
+                           int threshold)
+{
+    Jpeg2000TgtNode *stack[30];
+    int sp = -1, curval = 0;
+
+    while (node && !node->vis) {
+        stack[++sp] = node;
+        node        = node->parent;
+    }
+
+    if (node)
+        curval = node->val;
+    else
+        curval = stack[sp]->val;
+
+    while (curval < threshold && sp >= 0) {
+        if (curval < stack[sp]->val)
+            curval = stack[sp]->val;
+        while (curval < threshold) {
+            int ret;
+            if ((ret = get_bits(s, 1)) > 0) {
+                stack[sp]->vis++;
+                break;
+            } else if (!ret)
+                curval++;
+            else
+                return ret;
+        }
+        stack[sp]->val = curval;
+        sp--;
+    }
+    return curval;
+}
+
+/* marker segments */
+/* get sizes and offsets of image, tiles; number of components */
+static int get_siz(Jpeg2000DecoderContext *s)
+{
+    int i;
+
+    if (s->buf_end - s->buf < 36)
+        return AVERROR(EINVAL);
+
+    s->avctx->profile = bytestream_get_be16(&s->buf); // Rsiz
+    s->width          = bytestream_get_be32(&s->buf); // Width
+    s->height         = bytestream_get_be32(&s->buf); // Height
+    s->image_offset_x = bytestream_get_be32(&s->buf); // X0Siz
+    s->image_offset_y = bytestream_get_be32(&s->buf); // Y0Siz
+    s->tile_width     = bytestream_get_be32(&s->buf); // XTSiz
+    s->tile_height    = bytestream_get_be32(&s->buf); // YTSiz
+    s->tile_offset_x  = bytestream_get_be32(&s->buf); // XT0Siz
+    s->tile_offset_y  = bytestream_get_be32(&s->buf); // YT0Siz
+    s->ncomponents    = bytestream_get_be16(&s->buf); // CSiz
+
+    if (s->buf_end - s->buf < 2 * s->ncomponents)
+        return AVERROR(EINVAL);
+
+    for (i = 0; i < s->ncomponents; i++) { // Ssiz_i XRsiz_i, YRsiz_i
+        uint8_t x = bytestream_get_byte(&s->buf);
+        s->cbps[i]   = (x & 0x7f) + 1;
+        s->precision = FFMAX(s->cbps[i], s->precision);
+        s->sgnd[i]   = (x & 0x80) == 1;
+        s->cdx[i]    = bytestream_get_byte(&s->buf);
+        s->cdy[i]    = bytestream_get_byte(&s->buf);
+    }
+
+    s->numXtiles = ff_jpeg2000_ceildiv(s->width  - s->tile_offset_x, s->tile_width);
+    s->numYtiles = ff_jpeg2000_ceildiv(s->height - s->tile_offset_y, s->tile_height);
+
+    s->tile = av_mallocz(s->numXtiles * s->numYtiles * sizeof(*s->tile));
+    if (!s->tile)
+        return AVERROR(ENOMEM);
+
+    for (i = 0; i < s->numXtiles * s->numYtiles; i++) {
+        Jpeg2000Tile *tile = s->tile + i;
+
+        tile->comp = av_mallocz(s->ncomponents * sizeof(*tile->comp));
+        if (!tile->comp)
+            return AVERROR(ENOMEM);
+    }
+
+    /* compute image size with reduction factor */
+    s->avctx->width  = ff_jpeg2000_ceildivpow2(s->width  - s->image_offset_x,
+                                               s->reduction_factor);
+    s->avctx->height = ff_jpeg2000_ceildivpow2(s->height - s->image_offset_y,
+                                               s->reduction_factor);
+
+    switch (s->avctx->profile) {
+    case FF_PROFILE_JPEG2000_DCINEMA_2K:
+    case FF_PROFILE_JPEG2000_DCINEMA_4K:
+        /* XYZ color-space for digital cinema profiles */
+        s->avctx->pix_fmt = AV_PIX_FMT_XYZ12LE;
+        break;
+    default:
+        /* For other profiles selects color-space according number of
+         * components and bit depth precision. */
+        switch (s->ncomponents) {
+        case 1:
+            if (s->precision > 8)
+                s->avctx->pix_fmt = AV_PIX_FMT_GRAY16;
+            else
+                s->avctx->pix_fmt = AV_PIX_FMT_GRAY8;
+            break;
+        case 3:
+            if (s->precision > 8)
+                s->avctx->pix_fmt = AV_PIX_FMT_RGB48;
+            else
+                s->avctx->pix_fmt = AV_PIX_FMT_RGB24;
+            break;
+        case 4:
+            s->avctx->pix_fmt = AV_PIX_FMT_BGRA;
+            break;
+        default:
+            /* pixel format can not be identified */
+            s->avctx->pix_fmt = AV_PIX_FMT_NONE;
+            break;
+        }
+        break;
+    }
+    return 0;
+}
+
+/* get common part for COD and COC segments */
+static int get_cox(Jpeg2000DecoderContext *s, Jpeg2000CodingStyle *c)
+{
+    uint8_t byte;
+
+    if (s->buf_end - s->buf < 5)
+        return AVERROR(EINVAL);
+    c->nreslevels = bytestream_get_byte(&s->buf) + 1; // num of resolution levels - 1
+
+    /* compute number of resolution levels to decode */
+    if (c->nreslevels < s->reduction_factor)
+        c->nreslevels2decode = 1;
+    else
+        c->nreslevels2decode = c->nreslevels - s->reduction_factor;
+
+    c->log2_cblk_width  = bytestream_get_byte(&s->buf) + 2; // cblk width
+    c->log2_cblk_height = bytestream_get_byte(&s->buf) + 2; // cblk height
+
+    c->cblk_style = bytestream_get_byte(&s->buf);
+    if (c->cblk_style != 0) { // cblk style
+        av_log(s->avctx, AV_LOG_ERROR, "no extra cblk styles supported\n");
+        return -1;
+    }
+    c->transform = bytestream_get_byte(&s->buf); // DWT transformation type
+    /* set integer 9/7 DWT in case of BITEXACT flag */
+    if ((s->avctx->flags & CODEC_FLAG_BITEXACT) && (c->transform == FF_DWT97))
+        c->transform = FF_DWT97_INT;
+
+    if (c->csty & JPEG2000_CSTY_PREC) {
+        int i;
+        for (i = 0; i < c->nreslevels; i++) {
+            byte = bytestream_get_byte(&s->buf);
+            c->log2_prec_widths[i]  =  byte       & 0x0F;    // precinct PPx
+            c->log2_prec_heights[i] = (byte >> 4) & 0x0F;    // precinct PPy
+        }
+    }
+    return 0;
+}
+
+/* get coding parameters for a particular tile or whole image*/
+static int get_cod(Jpeg2000DecoderContext *s, Jpeg2000CodingStyle *c,
+                   uint8_t *properties)
+{
+    Jpeg2000CodingStyle tmp;
+    int compno;
+
+    if (s->buf_end - s->buf < 5)
+        return AVERROR(EINVAL);
+
+    tmp.log2_prec_width  =
+    tmp.log2_prec_height = 15;
+
+    tmp.csty = bytestream_get_byte(&s->buf);
+
+    // get progression order
+    tmp.prog_order = bytestream_get_byte(&s->buf);
+
+    tmp.nlayers = bytestream_get_be16(&s->buf);
+    tmp.mct     = bytestream_get_byte(&s->buf); // multiple component transformation
+
+    get_cox(s, &tmp);
+    for (compno = 0; compno < s->ncomponents; compno++)
+        if (!(properties[compno] & HAD_COC))
+            memcpy(c + compno, &tmp, sizeof(tmp));
+    return 0;
+}
+
+/* Get coding parameters for a component in the whole image or a
+ * particular tile. */
+static int get_coc(Jpeg2000DecoderContext *s, Jpeg2000CodingStyle *c,
+                   uint8_t *properties)
+{
+    int compno;
+
+    if (s->buf_end - s->buf < 2)
+        return AVERROR(EINVAL);
+
+    compno = bytestream_get_byte(&s->buf);
+
+    c      += compno;
+    c->csty = bytestream_get_byte(&s->buf);
+    get_cox(s, c);
+
+    properties[compno] |= HAD_COC;
+    return 0;
+}
+
+/* Get common part for QCD and QCC segments. */
+static int get_qcx(Jpeg2000DecoderContext *s, int n, Jpeg2000QuantStyle *q)
+{
+    int i, x;
+
+    if (s->buf_end - s->buf < 1)
+        return AVERROR(EINVAL);
+
+    x = bytestream_get_byte(&s->buf); // Sqcd
+
+    q->nguardbits = x >> 5;
+    q->quantsty   = x & 0x1f;
+
+    if (q->quantsty == JPEG2000_QSTY_NONE) {
+        n -= 3;
+        if (s->buf_end - s->buf < n)
+            return AVERROR(EINVAL);
+        for (i = 0; i < n; i++)
+            q->expn[i] = bytestream_get_byte(&s->buf) >> 3;
+    } else if (q->quantsty == JPEG2000_QSTY_SI) {
+        if (s->buf_end - s->buf < 2)
+            return AVERROR(EINVAL);
+        x          = bytestream_get_be16(&s->buf);
+        q->expn[0] = x >> 11;
+        q->mant[0] = x & 0x7ff;
+        for (i = 1; i < 32 * 3; i++) {
+            int curexpn = FFMAX(0, q->expn[0] - (i - 1) / 3);
+            q->expn[i] = curexpn;
+            q->mant[i] = q->mant[0];
+        }
+    } else {
+        n = (n - 3) >> 1;
+        if (s->buf_end - s->buf < n)
+            return AVERROR(EINVAL);
+        for (i = 0; i < n; i++) {
+            x          = bytestream_get_be16(&s->buf);
+            q->expn[i] = x >> 11;
+            q->mant[i] = x & 0x7ff;
+        }
+    }
+    return 0;
+}
+
+/* Get quantization parameters for a particular tile or a whole image. */
+static int get_qcd(Jpeg2000DecoderContext *s, int n, Jpeg2000QuantStyle *q,
+                   uint8_t *properties)
+{
+    Jpeg2000QuantStyle tmp;
+    int compno;
+
+    if (get_qcx(s, n, &tmp))
+        return -1;
+    for (compno = 0; compno < s->ncomponents; compno++)
+        if (!(properties[compno] & HAD_QCC))
+            memcpy(q + compno, &tmp, sizeof(tmp));
+    return 0;
+}
+
+/* Get quantization parameters for a component in the whole image
+ * on in a particular tile. */
+static int get_qcc(Jpeg2000DecoderContext *s, int n, Jpeg2000QuantStyle *q,
+                   uint8_t *properties)
+{
+    int compno;
+
+    if (s->buf_end - s->buf < 1)
+        return AVERROR(EINVAL);
+
+    compno              = bytestream_get_byte(&s->buf);
+    properties[compno] |= HAD_QCC;
+    return get_qcx(s, n - 1, q + compno);
+}
+
+/* Get start of tile segment. */
+static uint8_t get_sot(Jpeg2000DecoderContext *s, int n)
+{
+    Jpeg2000TilePart *tp;
+    uint16_t Isot;
+    uint32_t Psot;
+    uint8_t TPsot;
+
+    if (s->buf_end - s->buf < 4)
+        return AVERROR(EINVAL);
+
+    Isot = bytestream_get_be16(&s->buf);        // Isot
+    if (Isot) {
+        av_log(s->avctx, AV_LOG_ERROR,
+               "Not a DCINEMA JP2K file: more than one tile\n");
+        return -1;
+    }
+    Psot  = bytestream_get_be32(&s->buf);       // Psot
+    TPsot = bytestream_get_byte(&s->buf);       // TPsot
+
+    /* Read TNSot but not used */
+    bytestream_get_byte(&s->buf);               // TNsot
+
+    tp             = s->tile[s->curtileno].tile_part + TPsot;
+    tp->tile_index = Isot;
+    tp->tp_len     = Psot;
+    tp->tp_idx     = TPsot;
+
+    /* Start of bit stream. Pointer to SOD marker
+     * Check SOD marker is present. */
+    if (JPEG2000_SOD == bytestream_get_be16(&s->buf))
+        tp->tp_start_bstrm = s->buf;
+    else {
+        av_log(s->avctx, AV_LOG_ERROR, "SOD marker not found \n");
+        return -1;
+    }
+
+    /* End address of bit stream =
+     *     start address + (Psot - size of SOT HEADER(n)
+     *     - size of SOT MARKER(2)  - size of SOD marker(2) */
+    tp->tp_end_bstrm = s->buf + (tp->tp_len - n - 4);
+
+    // set buffer pointer to end of tile part header
+    s->buf = tp->tp_end_bstrm;
+
+    return 0;
+}
+
+/* Tile-part lengths: see ISO 15444-1:2002, section A.7.1
+ * Used to know the number of tile parts and lengths.
+ * There may be multiple TLMs in the header.
+ * TODO: The function is not used for tile-parts management, nor anywhere else.
+ * It can be useful to allocate memory for tile parts, before managing the SOT
+ * markers. Parsing the TLM header is needed to increment the input header
+ * buffer.
+ * This marker is mandatory for DCI. */
+static uint8_t get_tlm(Jpeg2000DecoderContext *s, int n)
+{
+    uint8_t Stlm, ST, SP, tile_tlm, i;
+    bytestream_get_byte(&s->buf);               /* Ztlm: skipped */
+    Stlm = bytestream_get_byte(&s->buf);
+
+    // too complex ? ST = ((Stlm >> 4) & 0x01) + ((Stlm >> 4) & 0x02);
+    ST = (Stlm >> 4) & 0x03;
+    // TODO: Manage case of ST = 0b11 --> raise error
+    SP       = (Stlm >> 6) & 0x01;
+    tile_tlm = (n - 4) / ((SP + 1) * 2 + ST);
+    for (i = 0; i < tile_tlm; i++) {
+        switch (ST) {
+        case 0:
+            break;
+        case 1:
+            bytestream_get_byte(&s->buf);
+            break;
+        case 2:
+            bytestream_get_be16(&s->buf);
+            break;
+        case 3:
+            bytestream_get_be32(&s->buf);
+            break;
+        }
+        if (SP == 0) {
+            bytestream_get_be16(&s->buf);
+        } else {
+            bytestream_get_be32(&s->buf);
+        }
+    }
+    return 0;
+}
+
+static int init_tile(Jpeg2000DecoderContext *s, int tileno)
+{
+    int compno;
+    int tilex = tileno % s->numXtiles;
+    int tiley = tileno / s->numXtiles;
+    Jpeg2000Tile *tile = s->tile + tileno;
+    Jpeg2000CodingStyle *codsty;
+    Jpeg2000QuantStyle  *qntsty;
+
+    if (!tile->comp)
+        return AVERROR(ENOMEM);
+
+    /* copy codsty, qnsty to tile. TODO: Is it the best way?
+     * codsty, qnsty is an array of 4 structs Jpeg2000CodingStyle
+     * and Jpeg2000QuantStyle */
+    memcpy(tile->codsty, s->codsty, s->ncomponents * sizeof(*codsty));
+    memcpy(tile->qntsty, s->qntsty, s->ncomponents * sizeof(*qntsty));
+
+    for (compno = 0; compno < s->ncomponents; compno++) {
+        Jpeg2000Component *comp = tile->comp + compno;
+        int ret; // global bandno
+        codsty = tile->codsty + compno;
+        qntsty = tile->qntsty + compno;
+
+        comp->coord_o[0][0] = FFMAX(tilex       * s->tile_width  + s->tile_offset_x, s->image_offset_x);
+        comp->coord_o[0][1] = FFMIN((tilex + 1) * s->tile_width  + s->tile_offset_x, s->width);
+        comp->coord_o[1][0] = FFMAX(tiley       * s->tile_height + s->tile_offset_y, s->image_offset_y);
+        comp->coord_o[1][1] = FFMIN((tiley + 1) * s->tile_height + s->tile_offset_y, s->height);
+
+        // FIXME: add a dcinema profile check ?
+        // value is guaranteed by profile (orig=0, 1 tile)
+        comp->coord[0][0] = 0;
+        comp->coord[0][1] = s->avctx->width;
+        comp->coord[1][0] = 0;
+        comp->coord[1][1] = s->avctx->height;
+
+        if (ret = ff_jpeg2000_init_component(comp, codsty, qntsty,
+                                             s->cbps[compno], s->cdx[compno],
+                                             s->cdy[compno], s->avctx))
+            return ret;
+    }
+    return 0;
+}
+
+/* Read the number of coding passes. */
+static int getnpasses(Jpeg2000DecoderContext *s)
+{
+    int num;
+    if (!get_bits(s, 1))
+        return 1;
+    if (!get_bits(s, 1))
+        return 2;
+    if ((num = get_bits(s, 2)) != 3)
+        return num < 0 ? num : 3 + num;
+    if ((num = get_bits(s, 5)) != 31)
+        return num < 0 ? num : 6 + num;
+    num = get_bits(s, 7);
+    return num < 0 ? num : 37 + num;
+}
+
+static int getlblockinc(Jpeg2000DecoderContext *s)
+{
+    int res = 0, ret;
+    while (ret = get_bits(s, 1)) {
+        if (ret < 0)
+            return ret;
+        res++;
+    }
+    return res;
+}
+
+static int jpeg2000_decode_packet(Jpeg2000DecoderContext *s,
+                                  Jpeg2000CodingStyle *codsty,
+                                  Jpeg2000ResLevel *rlevel, int precno,
+                                  int layno, uint8_t *expn, int numgbits)
+{
+    int bandno, cblkno, ret, nb_code_blocks;
+
+    if (!(ret = get_bits(s, 1))) {
+        jpeg2000_flush(s);
+        return 0;
+    } else if (ret < 0)
+        return ret;
+
+    for (bandno = 0; bandno < rlevel->nbands; bandno++) {
+        Jpeg2000Band *band = rlevel->band + bandno;
+        Jpeg2000Prec *prec = band->prec + precno;
+
+        if (band->coord[0][0] == band->coord[0][1] ||
+            band->coord[1][0] == band->coord[1][1])
+            continue;
+        prec->yi0 = 0;
+        prec->xi0 = 0;
+        nb_code_blocks =  prec->nb_codeblocks_height *
+                          prec->nb_codeblocks_width;
+        for (cblkno = 0; cblkno < nb_code_blocks; cblkno++) {
+            Jpeg2000Cblk *cblk = prec->cblk + cblkno;
+            int incl, newpasses, llen;
+
+            if (cblk->npasses)
+                incl = get_bits(s, 1);
+            else
+                incl = tag_tree_decode(s, prec->cblkincl + cblkno, layno + 1) == layno;
+            if (!incl)
+                continue;
+            else if (incl < 0)
+                return incl;
+
+            if (!cblk->npasses)
+                cblk->nonzerobits = expn[bandno] + numgbits - 1 -
+                                    tag_tree_decode(s, prec->zerobits + cblkno,
+                                                    100);
+            if ((newpasses = getnpasses(s)) < 0)
+                return newpasses;
+            if ((llen = getlblockinc(s)) < 0)
+                return llen;
+            cblk->lblock += llen;
+            if ((ret = get_bits(s, av_log2(newpasses) + cblk->lblock)) < 0)
+                return ret;
+            cblk->lengthinc = ret;
+            cblk->npasses  += newpasses;
+        }
+    }
+    jpeg2000_flush(s);
+
+    if (codsty->csty & JPEG2000_CSTY_EPH) {
+        if (AV_RB16(s->buf) == JPEG2000_EPH)
+            s->buf += 2;
+        else
+            av_log(s->avctx, AV_LOG_ERROR, "EPH marker not found.\n");
+    }
+
+    for (bandno = 0; bandno < rlevel->nbands; bandno++) {
+        Jpeg2000Band *band = rlevel->band + bandno;
+        Jpeg2000Prec *prec = band->prec + precno;
+
+        nb_code_blocks = prec->nb_codeblocks_height * prec->nb_codeblocks_width;
+        for (cblkno = 0; cblkno < nb_code_blocks; cblkno++) {
+            Jpeg2000Cblk *cblk = prec->cblk + cblkno;
+            if (s->buf_end - s->buf < cblk->lengthinc)
+                return AVERROR(EINVAL);
+            bytestream_get_buffer(&s->buf, cblk->data, cblk->lengthinc);
+            cblk->length   += cblk->lengthinc;
+            cblk->lengthinc = 0;
+        }
+    }
+    return 0;
+}
+
+static int jpeg2000_decode_packets(Jpeg2000DecoderContext *s, Jpeg2000Tile *tile)
+{
+    int layno, reslevelno, compno, precno, ok_reslevel;
+    uint8_t prog_order = tile->codsty[0].prog_order;
+    uint16_t x;
+    uint16_t y;
+
+    s->bit_index = 8;
+    switch (prog_order) {
+    case JPEG2000_PGOD_LRCP:
+        for (layno = 0; layno < tile->codsty[0].nlayers; layno++) {
+            ok_reslevel = 1;
+            for (reslevelno = 0; ok_reslevel; reslevelno++) {
+                ok_reslevel = 0;
+                for (compno = 0; compno < s->ncomponents; compno++) {
+                    Jpeg2000CodingStyle *codsty = tile->codsty + compno;
+                    Jpeg2000QuantStyle *qntsty  = tile->qntsty + compno;
+                    if (reslevelno < codsty->nreslevels) {
+                        Jpeg2000ResLevel *rlevel = tile->comp[compno].reslevel +
+                                                   reslevelno;
+                        ok_reslevel = 1;
+                        for (precno = 0; precno < rlevel->num_precincts_x * rlevel->num_precincts_y; precno++)
+                            if (jpeg2000_decode_packet(s,
+                                                       codsty, rlevel,
+                                                       precno, layno,
+                                                       qntsty->expn + (reslevelno ? 3 * (reslevelno - 1) + 1 : 0),
+                                                       qntsty->nguardbits))
+                                return -1;
+                    }
+                }
+            }
+        }
+        break;
+
+    case JPEG2000_PGOD_CPRL:
+        for (compno = 0; compno < s->ncomponents; compno++) {
+            Jpeg2000CodingStyle *codsty = tile->codsty + compno;
+            Jpeg2000QuantStyle *qntsty  = tile->qntsty + compno;
+
+            /* Set bit stream buffer address according to tile-part.
+             * For DCinema one tile-part per component, so can be
+             * indexed by component. */
+            s->buf = tile->tile_part[compno].tp_start_bstrm;
+
+            /* Position loop (y axis)
+             * TODO: Automate computing of step 256.
+             * Fixed here, but to be computed before entering here. */
+            for (y = 0; y < s->height; y += 256) {
+                /* Position loop (y axis)
+                 * TODO: automate computing of step 256.
+                 * Fixed here, but to be computed before entering here. */
+                for (x = 0; x < s->width; x += 256) {
+                    for (reslevelno = 0; reslevelno < codsty->nreslevels; reslevelno++) {
+                        uint16_t prcx, prcy;
+                        uint8_t reducedresno = codsty->nreslevels - 1 -reslevelno; //  ==> N_L - r
+                        Jpeg2000ResLevel *rlevel = tile->comp[compno].reslevel + reslevelno;
+
+                        if (!((y % (1 << (rlevel->log2_prec_height + reducedresno)) == 0) ||
+                              (y == 0))) // TODO: 2nd condition simplified as try0 always =0 for dcinema
+                            continue;
+
+                        if (!((x % (1 << (rlevel->log2_prec_width + reducedresno)) == 0) ||
+                              (x == 0))) // TODO: 2nd condition simplified as try0 always =0 for dcinema
+                            continue;
+
+                        // check if a precinct exists
+                        prcx   = ff_jpeg2000_ceildivpow2(x, reducedresno) >> rlevel->log2_prec_width;
+                        prcy   = ff_jpeg2000_ceildivpow2(y, reducedresno) >> rlevel->log2_prec_height;
+                        precno = prcx + rlevel->num_precincts_x * prcy;
+                        for (layno = 0; layno < tile->codsty[0].nlayers; layno++) {
+                            if (jpeg2000_decode_packet(s, codsty, rlevel,
+                                                       precno, layno,
+                                                       qntsty->expn + (reslevelno ? 3 * (reslevelno - 1) + 1 : 0),
+                                                       qntsty->nguardbits))
+                                return -1;
+                        }
+                    }
+                }
+            }
+        }
+        break;
+
+    default:
+        break;
+    }
+
+    /* EOC marker reached */
+    s->buf += 2;
+
+    return 0;
+}
+
+/* TIER-1 routines */
+static void decode_sigpass(Jpeg2000T1Context *t1, int width, int height,
+                           int bpno, int bandno)
+{
+    int mask = 3 << (bpno - 1), y0, x, y;
+
+    for (y0 = 0; y0 < height; y0 += 4)
+        for (x = 0; x < width; x++)
+            for (y = y0; y < height && y < y0 + 4; y++)
+                if ((t1->flags[y + 1][x + 1] & JPEG2000_T1_SIG_NB)
+                    && !(t1->flags[y + 1][x + 1] & (JPEG2000_T1_SIG | JPEG2000_T1_VIS))) {
+                    if (ff_mqc_decode(&t1->mqc,
+                                      t1->mqc.cx_states +
+                                      ff_jpeg2000_getsigctxno(t1->flags[y + 1][x + 1],
+                                                             bandno))) {
+                        int xorbit, ctxno = ff_jpeg2000_getsgnctxno(t1->flags[y + 1][x + 1],
+                                                                    &xorbit);
+
+                        t1->data[y][x] =
+                            (ff_mqc_decode(&t1->mqc,
+                                           t1->mqc.cx_states + ctxno) ^ xorbit)
+                            ? -mask : mask;
+
+                        ff_jpeg2000_set_significance(t1, x, y,
+                                                     t1->data[y][x] < 0);
+                    }
+                    t1->flags[y + 1][x + 1] |= JPEG2000_T1_VIS;
+                }
+}
+
+static void decode_refpass(Jpeg2000T1Context *t1, int width, int height,
+                           int bpno)
+{
+    int phalf, nhalf;
+    int y0, x, y;
+
+    phalf = 1 << (bpno - 1);
+    nhalf = -phalf;
+
+    for (y0 = 0; y0 < height; y0 += 4)
+        for (x = 0; x < width; x++)
+            for (y = y0; y < height && y < y0 + 4; y++)
+                if ((t1->flags[y + 1][x + 1] & (JPEG2000_T1_SIG | JPEG2000_T1_VIS)) == JPEG2000_T1_SIG) {
+                    int ctxno = ff_jpeg2000_getrefctxno(t1->flags[y + 1][x + 1]);
+                    int r     = ff_mqc_decode(&t1->mqc,
+                                              t1->mqc.cx_states + ctxno)
+                                ? phalf : nhalf;
+                    t1->data[y][x]          += t1->data[y][x] < 0 ? -r : r;
+                    t1->flags[y + 1][x + 1] |= JPEG2000_T1_REF;
+                }
+}
+
+static void decode_clnpass(Jpeg2000DecoderContext *s, Jpeg2000T1Context *t1,
+                           int width, int height, int bpno, int bandno,
+                           int seg_symbols)
+{
+    int mask = 3 << (bpno - 1), y0, x, y, runlen, dec;
+
+    for (y0 = 0; y0 < height; y0 += 4)
+        for (x = 0; x < width; x++) {
+            if (y0 + 3 < height &&
+                !((t1->flags[y0 + 1][x + 1] & (JPEG2000_T1_SIG_NB | JPEG2000_T1_VIS | JPEG2000_T1_SIG)) ||
+                  (t1->flags[y0 + 2][x + 1] & (JPEG2000_T1_SIG_NB | JPEG2000_T1_VIS | JPEG2000_T1_SIG)) ||
+                  (t1->flags[y0 + 3][x + 1] & (JPEG2000_T1_SIG_NB | JPEG2000_T1_VIS | JPEG2000_T1_SIG)) ||
+                  (t1->flags[y0 + 4][x + 1] & (JPEG2000_T1_SIG_NB | JPEG2000_T1_VIS | JPEG2000_T1_SIG)))) {
+                if (!ff_mqc_decode(&t1->mqc, t1->mqc.cx_states + MQC_CX_RL))
+                    continue;
+                runlen = ff_mqc_decode(&t1->mqc,
+                                       t1->mqc.cx_states + MQC_CX_UNI);
+                runlen = (runlen << 1) | ff_mqc_decode(&t1->mqc,
+                                                       t1->mqc.cx_states +
+                                                       MQC_CX_UNI);
+                dec = 1;
+            } else {
+                runlen = 0;
+                dec    = 0;
+            }
+
+            for (y = y0 + runlen; y < y0 + 4 && y < height; y++) {
+                if (!dec) {
+                    if (!(t1->flags[y + 1][x + 1] & (JPEG2000_T1_SIG | JPEG2000_T1_VIS)))
+                        dec = ff_mqc_decode(&t1->mqc,
+                                            t1->mqc.cx_states +
+                                            ff_jpeg2000_getsigctxno(t1->flags[y + 1][x + 1],
+                                                                   bandno));
+                }
+                if (dec) {
+                    int xorbit;
+                    int ctxno = ff_jpeg2000_getsgnctxno(t1->flags[y + 1][x + 1],
+                                                        &xorbit);
+                    t1->data[y][x] = (ff_mqc_decode(&t1->mqc,
+                                                    t1->mqc.cx_states + ctxno) ^
+                                      xorbit)
+                                     ? -mask : mask;
+                    ff_jpeg2000_set_significance(t1, x, y, t1->data[y][x] < 0);
+                }
+                dec = 0;
+                t1->flags[y + 1][x + 1] &= ~JPEG2000_T1_VIS;
+            }
+        }
+    if (seg_symbols) {
+        int val;
+        val = ff_mqc_decode(&t1->mqc, t1->mqc.cx_states + MQC_CX_UNI);
+        val = (val << 1) + ff_mqc_decode(&t1->mqc, t1->mqc.cx_states + MQC_CX_UNI);
+        val = (val << 1) + ff_mqc_decode(&t1->mqc, t1->mqc.cx_states + MQC_CX_UNI);
+        val = (val << 1) + ff_mqc_decode(&t1->mqc, t1->mqc.cx_states + MQC_CX_UNI);
+        if (val != 0xa)
+            av_log(s->avctx, AV_LOG_ERROR,
+                   "Segmentation symbol value incorrect\n");
+    }
+}
+
+static int decode_cblk(Jpeg2000DecoderContext *s, Jpeg2000CodingStyle *codsty,
+                       Jpeg2000T1Context *t1, Jpeg2000Cblk *cblk,
+                       int width, int height, int bandpos)
+{
+    int passno = cblk->npasses, pass_t = 2, bpno = cblk->nonzerobits - 1, y;
+
+    for (y = 0; y < height + 2; y++)
+        memset(t1->flags[y], 0, (width + 2) * sizeof(width));
+
+    for (y = 0; y < height; y++)
+        memset(t1->data[y], 0, width * sizeof(width));
+
+    ff_mqc_initdec(&t1->mqc, cblk->data);
+    cblk->data[cblk->length]     = 0xff;
+    cblk->data[cblk->length + 1] = 0xff;
+
+    while (passno--) {
+        switch (pass_t) {
+        case 0:
+            decode_sigpass(t1, width, height, bpno + 1, bandpos);
+            break;
+        case 1:
+            decode_refpass(t1, width, height, bpno + 1);
+            break;
+        case 2:
+            decode_clnpass(s, t1, width, height, bpno + 1, bandpos,
+                           codsty->cblk_style & JPEG2000_CBLK_SEGSYM);
+            break;
+        }
+
+        pass_t++;
+        if (pass_t == 3) {
+            bpno--;
+            pass_t = 0;
+        }
+    }
+    return 0;
+}
+
+/* TODO: Verify dequantization for lossless case
+ * comp->data can be float or int
+ * band->stepsize can be float or int
+ * depending on the type of DWT transformation.
+ * see ISO/IEC 15444-1:2002 A.6.1 */
+
+/* Float dequantization of a codeblock.*/
+static void dequantization_float(int x, int y, Jpeg2000Cblk *cblk,
+                                 Jpeg2000Component *comp,
+                                 Jpeg2000T1Context *t1, Jpeg2000Band *band)
+{
+    int i, j, idx;
+    float *datap = &comp->data[(comp->coord[0][1] - comp->coord[0][0]) * y + x];
+    for (j = 0; j < (cblk->coord[1][1] - cblk->coord[1][0]); ++j)
+        for (i = 0; i < (cblk->coord[0][1] - cblk->coord[0][0]); ++i) {
+            idx        = (comp->coord[0][1] - comp->coord[0][0]) * j + i;
+            datap[idx] = (float)(t1->data[j][i]) * ((float)band->stepsize);
+        }
+    return;
+}
+
+/* Integer dequantization of a codeblock.*/
+static void dequantization_int(int x, int y, Jpeg2000Cblk *cblk,
+                               Jpeg2000Component *comp,
+                               Jpeg2000T1Context *t1, Jpeg2000Band *band)
+{
+    int i, j, idx;
+    int32_t *datap =
+        (int32_t *) &comp->data[(comp->coord[0][1] - comp->coord[0][0]) * y + x];
+    for (j = 0; j < (cblk->coord[1][1] - cblk->coord[1][0]); ++j)
+        for (i = 0; i < (cblk->coord[0][1] - cblk->coord[0][0]); ++i) {
+            idx        = (comp->coord[0][1] - comp->coord[0][0]) * j + i;
+            datap[idx] =
+                ((int32_t)(t1->data[j][i]) * ((int32_t)band->stepsize) + (1 << 15)) >> 16;
+        }
+    return;
+}
+
+/* Inverse ICT parameters in float and integer.
+ * int value = (float value) * (1<<16) */
+static const float f_ict_params[4] = {
+    1.402f,
+    0.34413f,
+    0.71414f,
+    1.772f
+};
+static const int   i_ict_params[4] = {
+     91881,
+     22553,
+     46802,
+    116130
+};
+
+static int mct_decode(Jpeg2000DecoderContext *s, Jpeg2000Tile *tile)
+{
+    int i, csize = 1;
+    int ret = 0;
+    int32_t *src[3],  i0,  i1,  i2;
+    float   *srcf[3], i0f, i1f, i2f;
+
+    for (i = 0; i < 3; i++)
+        if (tile->codsty[0].transform == FF_DWT97)
+            srcf[i] = tile->comp[i].data;
+        else
+            src[i] = (int32_t *)tile->comp[i].data;
+
+    for (i = 0; i < 2; i++)
+        csize *= tile->comp[0].coord[i][1] - tile->comp[0].coord[i][0];
+    switch (tile->codsty[0].transform) {
+    case FF_DWT97:
+        for (i = 0; i < csize; i++) {
+            i0f = *srcf[0] + (f_ict_params[0] * *srcf[2]);
+            i1f = *srcf[0] - (f_ict_params[1] * *srcf[1])
+                           - (f_ict_params[2] * *srcf[2]);
+            i2f = *srcf[0] + (f_ict_params[3] * *srcf[1]);
+            *srcf[0]++ = i0f;
+            *srcf[1]++ = i1f;
+            *srcf[2]++ = i2f;
+        }
+        break;
+    case FF_DWT97_INT:
+        for (i = 0; i < csize; i++) {
+            i0 = *src[0] + (((i_ict_params[0] * *src[2]) + (1 << 15)) >> 16);
+            i1 = *src[0] - (((i_ict_params[1] * *src[1]) + (1 << 15)) >> 16)
+                         - (((i_ict_params[2] * *src[2]) + (1 << 15)) >> 16);
+            i2 = *src[0] + (((i_ict_params[3] * *src[1]) + (1 << 15)) >> 16);
+            *src[0]++ = i0;
+            *src[1]++ = i1;
+            *src[2]++ = i2;
+        }
+        break;
+    case FF_DWT53:
+        for (i = 0; i < csize; i++) {
+            i1 = *src[0] - (*src[2] + *src[1] >> 2);
+            i0 = i1 + *src[2];
+            i2 = i1 + *src[1];
+            *src[0]++ = i0;
+            *src[1]++ = i1;
+            *src[2]++ = i2;
+        }
+        break;
+    }
+    return ret;
+}
+
+static int jpeg2000_decode_tile(Jpeg2000DecoderContext *s, Jpeg2000Tile *tile,
+                                AVFrame *picture)
+{
+    int compno, reslevelno, bandno;
+    int x, y;
+
+    uint8_t *line;
+    Jpeg2000T1Context t1;
+    /* Loop on tile components */
+
+    for (compno = 0; compno < s->ncomponents; compno++) {
+        Jpeg2000Component *comp     = tile->comp + compno;
+        Jpeg2000CodingStyle *codsty = tile->codsty + compno;
+        /* Loop on resolution levels */
+        for (reslevelno = 0; reslevelno < codsty->nreslevels2decode; reslevelno++) {
+            Jpeg2000ResLevel *rlevel = comp->reslevel + reslevelno;
+            /* Loop on bands */
+            for (bandno = 0; bandno < rlevel->nbands; bandno++) {
+                uint16_t nb_precincts, precno;
+                Jpeg2000Band *band = rlevel->band + bandno;
+                int cblkno = 0, bandpos;
+                bandpos = bandno + (reslevelno > 0);
+
+                nb_precincts = rlevel->num_precincts_x * rlevel->num_precincts_y;
+                /* Loop on precincts */
+                for (precno = 0; precno < nb_precincts; precno++) {
+                    Jpeg2000Prec *prec = band->prec + precno;
+
+                    /* Loop on codeblocks */
+                    for (cblkno = 0; cblkno < prec->nb_codeblocks_width * prec->nb_codeblocks_height; cblkno++) {
+                        int x, y;
+                        Jpeg2000Cblk *cblk = prec->cblk + cblkno;
+                        decode_cblk(s, codsty, &t1, cblk,
+                                    cblk->coord[0][1] - cblk->coord[0][0],
+                                    cblk->coord[1][1] - cblk->coord[1][0],
+                                    bandpos);
+
+                        /* Manage band offsets */
+                        x = cblk->coord[0][0];
+                        y = cblk->coord[1][0];
+                        if ((reslevelno > 0) && ((bandno + 1) & 1)) {
+                            Jpeg2000ResLevel *pres = comp->reslevel + (reslevelno - 1);
+                            x += pres->coord[0][1] - pres->coord[0][0];
+                        }
+                        if ((reslevelno > 0) && ((bandno + 1) & 2)) {
+                            Jpeg2000ResLevel *pres = comp->reslevel + (reslevelno - 1);
+                            y += pres->coord[1][1] - pres->coord[1][0];
+                        }
+
+                        if (s->avctx->flags & CODEC_FLAG_BITEXACT)
+                            dequantization_int(x, y, cblk, comp, &t1, band);
+                        else
+                            dequantization_float(x, y, cblk, comp, &t1, band);
+                   } /* end cblk */
+                } /*end prec */
+            } /* end band */
+        } /* end reslevel */
+
+        /* inverse DWT */
+        ff_dwt_decode(&comp->dwt, comp->data);
+    } /*end comp */
+
+    /* inverse MCT transformation */
+    if (tile->codsty[0].mct)
+        mct_decode(s, tile);
+
+    if (s->avctx->pix_fmt == PIX_FMT_BGRA) // RGBA -> BGRA
+        FFSWAP(float *, tile->comp[0].data, tile->comp[2].data);
+
+    if (s->precision <= 8) {
+        for (compno = 0; compno < s->ncomponents; compno++) {
+            Jpeg2000Component *comp = tile->comp + compno;
+            int32_t *datap = (int32_t *)comp->data;
+            y    = tile->comp[compno].coord[1][0] - s->image_offset_y;
+            line = picture->data[0] + y * picture->linesize[0];
+            for (; y < tile->comp[compno].coord[1][1] - s->image_offset_y; y += s->cdy[compno]) {
+                uint8_t *dst;
+
+                x   = tile->comp[compno].coord[0][0] - s->image_offset_x;
+                dst = line + x * s->ncomponents + compno;
+
+                for (; x < tile->comp[compno].coord[0][1] - s->image_offset_x; x += s->cdx[compno]) {
+                    *datap += 1 << (s->cbps[compno] - 1);
+                    if (*datap < 0)
+                        *datap = 0;
+                    else if (*datap >= (1 << s->cbps[compno]))
+                        *datap = (1 << s->cbps[compno]) - 1;
+                    *dst = *datap++;
+                    dst += s->ncomponents;
+                }
+                line += picture->linesize[0];
+            }
+        }
+    } else {
+        for (compno = 0; compno < s->ncomponents; compno++) {
+            Jpeg2000Component *comp = tile->comp + compno;
+            float *datap = comp->data;
+            int32_t *i_datap = (int32_t *) comp->data;
+            uint16_t *linel;
+
+            y     = tile->comp[compno].coord[1][0] - s->image_offset_y;
+            linel = (uint16_t *)picture->data[0] + y * (picture->linesize[0] >> 1);
+            for (; y < tile->comp[compno].coord[1][1] - s->image_offset_y; y += s->cdy[compno]) {
+                uint16_t *dst;
+                x   = tile->comp[compno].coord[0][0] - s->image_offset_x;
+                dst = linel + (x * s->ncomponents + compno);
+                for (; x < s->avctx->width; x += s->cdx[compno]) {
+                    int16_t val;
+                    /* DC level shift and clip see ISO 15444-1:2002 G.1.2 */
+                    if (s->avctx->flags & CODEC_FLAG_BITEXACT)
+                        val = *i_datap + (1 << (s->cbps[compno] - 1));
+                    else
+                        val = lrintf(*datap) + (1 << (s->cbps[compno] - 1));
+                    val = av_clip(val, 0, (1 << s->cbps[compno]) - 1);
+                    /* align 12 bit values in little-endian mode */
+                    *dst = val << 4;
+                    datap++;
+                    i_datap++;
+                    dst += s->ncomponents;
+                }
+                linel += picture->linesize[0] >> 1;
+            }
+        }
+    }
+    return 0;
+}
+
+static void jpeg2000_dec_cleanup(Jpeg2000DecoderContext *s)
+{
+    int tileno, compno;
+    for (tileno = 0; tileno < s->numXtiles * s->numYtiles; tileno++) {
+        for (compno = 0; compno < s->ncomponents; compno++) {
+            Jpeg2000Component *comp     = s->tile[tileno].comp   + compno;
+            Jpeg2000CodingStyle *codsty = s->tile[tileno].codsty + compno;
+
+            ff_jpeg2000_cleanup(comp, codsty);
+        }
+        av_freep(&s->tile[tileno].comp);
+    }
+    av_freep(&s->tile);
+}
+
+static int jpeg2000_read_main_headers(Jpeg2000DecoderContext *s)
+{
+    Jpeg2000CodingStyle *codsty = s->codsty;
+    Jpeg2000QuantStyle *qntsty  = s->qntsty;
+    uint8_t *properties         = s->properties;
+
+    for (;;) {
+        int len, ret = 0;
+        uint16_t marker;
+        const uint8_t *oldbuf;
+
+        if (s->buf_end - s->buf < 2) {
+            av_log(s->avctx, AV_LOG_ERROR, "Missing EOC\n");
+            break;
+        }
+
+        marker = bytestream_get_be16(&s->buf);
+        oldbuf = s->buf;
+
+        if (marker == JPEG2000_EOC)
+            break;
+
+        if (s->buf_end - s->buf < 2)
+            return AVERROR(EINVAL);
+        len = bytestream_get_be16(&s->buf);
+        switch (marker) {
+        case JPEG2000_SIZ:
+            ret = get_siz(s);
+            break;
+        case JPEG2000_COC:
+            ret = get_coc(s, codsty, properties);
+            break;
+        case JPEG2000_COD:
+            ret = get_cod(s, codsty, properties);
+            break;
+        case JPEG2000_QCC:
+            ret = get_qcc(s, len, qntsty, properties);
+            break;
+        case JPEG2000_QCD:
+            ret = get_qcd(s, len, qntsty, properties);
+            break;
+        case JPEG2000_SOT:
+            ret = get_sot(s, len);
+            break;
+        case JPEG2000_COM:
+            // the comment is ignored
+            s->buf += len - 2;
+            break;
+        case JPEG2000_TLM:
+            // Tile-part lengths
+            ret = get_tlm(s, len);
+            break;
+        default:
+            av_log(s->avctx, AV_LOG_ERROR,
+                   "unsupported marker 0x%.4X at pos 0x%lX\n",
+                   marker, (uint64_t)(s->buf - s->buf_start - 4));
+            s->buf += len - 2;
+            break;
+        }
+        if (((s->buf - oldbuf != len) && (marker != JPEG2000_SOT)) || ret) {
+            av_log(s->avctx, AV_LOG_ERROR,
+                   "error during processing marker segment %.4x\n", marker);
+            return ret ? ret : -1;
+        }
+    }
+    return 0;
+}
+
+/* Read bit stream packets --> T2 operation. */
+static int jpeg2000_read_bitstream_packets(Jpeg2000DecoderContext *s)
+{
+    int ret = 0;
+    Jpeg2000Tile *tile = s->tile + s->curtileno;
+
+    if (ret = init_tile(s, s->curtileno))
+        return ret;
+    if (ret = jpeg2000_decode_packets(s, tile))
+        return ret;
+
+    return 0;
+}
+
+static int jp2_find_codestream(Jpeg2000DecoderContext *s)
+{
+    int32_t atom_size;
+    int found_codestream = 0, search_range = 10;
+
+    // Skip JPEG 2000 signature atom.
+    s->buf += 12;
+
+    while (!found_codestream && search_range) {
+        atom_size = AV_RB32(s->buf);
+        if (AV_RB32(s->buf + 4) == JP2_CODESTREAM) {
+            found_codestream = 1;
+            s->buf += 8;
+        } else {
+            s->buf += atom_size;
+            search_range--;
+        }
+    }
+
+    if (found_codestream)
+        return 1;
+    return 0;
+}
+
+static int jpeg2000_decode_frame(AVCodecContext *avctx, void *data,
+                                 int *got_frame, AVPacket *avpkt)
+{
+    Jpeg2000DecoderContext *s = avctx->priv_data;
+    AVFrame *picture = data;
+    int tileno, ret;
+
+    s->avctx     = avctx;
+    s->buf       = s->buf_start = avpkt->data;
+    s->buf_end   = s->buf_start + avpkt->size;
+    s->curtileno = 0; // TODO: only one tile in DCI JP2K. to implement for more tiles
+
+    // reduction factor, i.e number of resolution levels to skip
+    s->reduction_factor = s->lowres;
+
+    ff_jpeg2000_init_tier1_luts();
+
+    if (s->buf_end - s->buf < 2)
+        return AVERROR(EINVAL);
+
+    // check if the image is in jp2 format
+    if ((AV_RB32(s->buf) == 12) &&
+        (AV_RB32(s->buf + 4) == JP2_SIG_TYPE) &&
+        (AV_RB32(s->buf + 8) == JP2_SIG_VALUE)) {
+        if (!jp2_find_codestream(s)) {
+            av_log(avctx, AV_LOG_ERROR,
+                   "couldn't find jpeg2k codestream atom\n");
+            return -1;
+        }
+    }
+
+    if (bytestream_get_be16(&s->buf) != JPEG2000_SOC) {
+        av_log(avctx, AV_LOG_ERROR, "SOC marker not present\n");
+        return -1;
+    }
+    if (ret = jpeg2000_read_main_headers(s))
+        return ret;
+
+    /* get picture buffer */
+    if ((ret = ff_get_buffer(avctx, picture, 0)) < 0) {
+        av_log(avctx, AV_LOG_ERROR, "ff_get_buffer() failed\n");
+        return ret;
+    }
+    picture->pict_type = AV_PICTURE_TYPE_I;
+    picture->key_frame = 1;
+
+    if (ret = jpeg2000_read_bitstream_packets(s))
+        return ret;
+    for (tileno = 0; tileno < s->numXtiles * s->numYtiles; tileno++)
+        if (ret = jpeg2000_decode_tile(s, s->tile + tileno, picture))
+            return ret;
+    jpeg2000_dec_cleanup(s);
+
+    *got_frame = 1;
+
+    return s->buf - s->buf_start;
+}
+
+#define OFFSET(x) offsetof(Jpeg2000DecoderContext, x)
+#define VD AV_OPT_FLAG_VIDEO_PARAM | AV_OPT_FLAG_DECODING_PARAM
+
+static const AVOption options[] = {
+    { "lowres",  "Lower the decoding resolution by a power of two",
+        OFFSET(lowres), AV_OPT_TYPE_INT, { .i64 = 0 }, 0, JPEG2000_MAX_RESLEVELS - 1, VD },
+    { NULL },
+};
+
+static const AVProfile profiles[] = {
+    { FF_PROFILE_JPEG2000_CSTREAM_RESTRICTION_0,  "JPEG 2000 codestream restriction 0"   },
+    { FF_PROFILE_JPEG2000_CSTREAM_RESTRICTION_1,  "JPEG 2000 codestream restriction 1"   },
+    { FF_PROFILE_JPEG2000_CSTREAM_NO_RESTRICTION, "JPEG 2000 no codestream restrictions" },
+    { FF_PROFILE_JPEG2000_DCINEMA_2K,             "JPEG 2000 digital cinema 2K"          },
+    { FF_PROFILE_JPEG2000_DCINEMA_4K,             "JPEG 2000 digital cinema 4K"          },
+    { FF_PROFILE_UNKNOWN },
+};
+
+static const AVClass class = {
+    .class_name = "jpeg2000",
+    .item_name  = av_default_item_name,
+    .option     = options,
+    .version    = LIBAVUTIL_VERSION_INT,
+};
+
+AVCodec ff_jpeg2000_decoder = {
+    .name           = "jpeg2000",
+    .long_name      = NULL_IF_CONFIG_SMALL("JPEG 2000"),
+    .type           = AVMEDIA_TYPE_VIDEO,
+    .id             = AV_CODEC_ID_JPEG2000,
+    .priv_data_size = sizeof(Jpeg2000DecoderContext),
+    .decode         = jpeg2000_decode_frame,
+    .priv_class     = &class,
+    .pix_fmts       = (enum PixelFormat[]) { AV_PIX_FMT_XYZ12,
+                                             AV_PIX_FMT_GRAY8,
+                                             -1 },
+    .profiles       = NULL_IF_CONFIG_SMALL(profiles)
+};

libavcodec/jpeg2000dwt.c

+/*
+ * Discrete wavelet transform
+ * Copyright (c) 2007 Kamil Nowosad
+ * Copyright (c) 2013 Nicolas Bertrand <nicoinattendu@gmail.com>
+ *
+ * This file is part of Libav.
+ *
+ * Libav is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * Libav is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with Libav; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+/**
+ * @file
+ * Discrete wavelet transform
+ */
+
+#include "libavutil/common.h"
+#include "libavutil/mem.h"
+#include "jpeg2000dwt.h"
+#include "internal.h"
+
+/* Defines for 9/7 DWT lifting parameters.
+ * Parameters are in float. */
+#define F_LFTG_ALPHA  1.586134342059924f
+#define F_LFTG_BETA   0.052980118572961f
+#define F_LFTG_GAMMA  0.882911075530934f
+#define F_LFTG_DELTA  0.443506852043971f
+#define F_LFTG_K      1.230174104914001f
+#define F_LFTG_X      1.625732422f
+/* FIXME: Why use 1.625732422 instead of 1/F_LFTG_K?
+ * Incorrect value in JPEG2000 norm.
+ * see (ISO/IEC 15444:1 (version 2002) F.3.8.2 */
+
+/* Lifting parameters in integer format.
+ * Computed as param = (float param) * (1 << 16) */
+#define I_LFTG_ALPHA  103949
+#define I_LFTG_BETA     3472
+#define I_LFTG_GAMMA   57862
+#define I_LFTG_DELTA   29066
+#define I_LFTG_K       80621
+#define I_LFTG_X      106544
+
+
+static inline void extend53(int *p, int i0, int i1)
+{
+    p[i0 - 1] = p[i0 + 1];
+    p[i1]     = p[i1 - 2];
+    p[i0 - 2] = p[i0 + 2];
+    p[i1 + 1] = p[i1 - 3];
+}
+
+static inline void extend97_float(float *p, int i0, int i1)
+{
+    int i;
+
+    for (i = 1; i <= 4; i++) {
+        p[i0 - i]     = p[i0 + i];
+        p[i1 + i - 1] = p[i1 - i - 1];
+    }
+}
+
+static inline void extend97_int(int32_t *p, int i0, int i1)
+{
+    int i;
+
+    for (i = 1; i <= 4; i++) {
+        p[i0 - i]     = p[i0 + i];
+        p[i1 + i - 1] = p[i1 - i - 1];
+    }
+}
+
+static void sr_1d53(int *p, int i0, int i1)
+{
+    int i;
+
+    if (i1 == i0 + 1)
+        return;
+
+    extend53(p, i0, i1);
+
+    for (i = i0 / 2; i < i1 / 2 + 1; i++)
+        p[2 * i] -= (p[2 * i - 1] + p[2 * i + 1] + 2) >> 2;
+    for (i = i0 / 2; i < i1 / 2; i++)
+        p[2 * i + 1] += (p[2 * i] + p[2 * i + 2]) >> 1;
+}
+
+static void dwt_decode53(DWTContext *s, int *t)
+{
+    int lev;
+    int w     = s->linelen[s->ndeclevels - 1][0];
+    int32_t *line = s->i_linebuf;
+    line += 3;
+
+    for (lev = 0; lev < s->ndeclevels; lev++) {
+        int lh = s->linelen[lev][0],
+            lv = s->linelen[lev][1],
+            mh = s->mod[lev][0],
+            mv = s->mod[lev][1],
+            lp;
+        int *l;
+
+        // HOR_SD
+        l = line + mh;
+        for (lp = 0; lp < lv; lp++) {
+            int i, j = 0;
+            // copy with interleaving
+            for (i = mh; i < lh; i += 2, j++)
+                l[i] = t[w * lp + j];
+            for (i = 1 - mh; i < lh; i += 2, j++)
+                l[i] = t[w * lp + j];
+
+            sr_1d53(line, mh, mh + lh);
+
+            for (i = 0; i < lh; i++)
+                t[w * lp + i] = l[i];
+        }
+
+        // VER_SD
+        l = line + mv;
+        for (lp = 0; lp < lh; lp++) {
+            int i, j = 0;
+            // copy with interleaving
+            for (i = mv; i < lv; i += 2, j++)
+                l[i] = t[w * j + lp];
+            for (i = 1 - mv; i < lv; i += 2, j++)
+                l[i] = t[w * j + lp];
+
+            sr_1d53(line, mv, mv + lv);
+
+            for (i = 0; i < lv; i++)
+                t[w * i + lp] = l[i];
+        }
+    }
+}
+
+static void sr_1d97_float(float *p, int i0, int i1)
+{
+    int i;
+
+    if (i1 == i0 + 1)
+        return;
+
+    extend97_float(p, i0, i1);
+
+    /*step 1*/
+    for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
+        p[2 * i]     *= F_LFTG_K;
+    /* step 2*/
+    for (i = i0 / 2 - 2; i < i1 / 2 + 2; i++)
+        p[2 * i + 1] *= F_LFTG_X;
+    /* step 3*/
+    for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
+        p[2 * i]     -= F_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]);
+    /* step 4 */
+    for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
+        p[2 * i + 1] -= F_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]);
+    /*step 5*/
+    for (i = i0 / 2; i < i1 / 2 + 1; i++)
+        p[2 * i]     += F_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]);
+    /* step 6 */
+    for (i = i0 / 2; i < i1 / 2; i++)
+        p[2 * i + 1] += F_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]);
+}
+
+static void dwt_decode97_float(DWTContext *s, float *t)
+{
+    int lev;
+    int w       = s->linelen[s->ndeclevels - 1][0];
+    float *line = s->f_linebuf;
+    float *data = t;
+    /* position at index O of line range [0-5,w+5] cf. extend function */
+    line += 5;
+
+    for (lev = 0; lev < s->ndeclevels; lev++) {
+        int lh = s->linelen[lev][0],
+            lv = s->linelen[lev][1],
+            mh = s->mod[lev][0],
+            mv = s->mod[lev][1],
+            lp;
+        float *l;
+        // HOR_SD
+        l = line + mh;
+        for (lp = 0; lp < lv; lp++) {
+            int i, j = 0;
+            // copy with interleaving
+            for (i = mh; i < lh; i += 2, j++)
+                l[i] = data[w * lp + j];
+            for (i = 1 - mh; i < lh; i += 2, j++)
+                l[i] = data[w * lp + j];
+
+            sr_1d97_float(line, mh, mh + lh);
+
+            for (i = 0; i < lh; i++)
+                data[w * lp + i] = l[i];
+        }
+
+        // VER_SD
+        l = line + mv;
+        for (lp = 0; lp < lh; lp++) {
+            int i, j = 0;
+            // copy with interleaving
+            for (i = mv; i < lv; i += 2, j++)
+                l[i] = data[w * j + lp];
+            for (i = 1 - mv; i < lv; i += 2, j++)
+                l[i] = data[w * j + lp];
+
+            sr_1d97_float(line, mv, mv + lv);
+
+            for (i = 0; i < lv; i++)
+                data[w * i + lp] = l[i];
+        }
+    }
+}
+
+static void sr_1d97_int(int32_t *p, int i0, int i1)
+{
+    int i;
+
+    if (i1 == i0 + 1)
+        return;
+
+    extend97_int(p, i0, i1);
+
+    /*step 1*/
+    for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
+        p[2 * i]      = ((p[2 * i] * I_LFTG_K)     + (1 << 15)) >> 16;
+    /* step 2*/
+    for (i = i0 / 2 - 2; i < i1 / 2 + 2; i++)
+        p[2 * i + 1]  = ((p[2 * i + 1] * I_LFTG_X) + (1 << 15)) >> 16;
+    /* step 3*/
+    for (i = i0 / 2 - 1; i < i1 / 2 + 2; i++)
+        p[2 * i]     -= (I_LFTG_DELTA * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
+    /* step 4 */
+    for (i = i0 / 2 - 1; i < i1 / 2 + 1; i++)
+        p[2 * i + 1] -= (I_LFTG_GAMMA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
+    /*step 5*/
+    for (i = i0 / 2; i < i1 / 2 + 1; i++)
+        p[2 * i]     += (I_LFTG_BETA  * (p[2 * i - 1] + p[2 * i + 1]) + (1 << 15)) >> 16;
+    /* step 6 */
+    for (i = i0 / 2; i < i1 / 2; i++)
+        p[2 * i + 1] += (I_LFTG_ALPHA * (p[2 * i]     + p[2 * i + 2]) + (1 << 15)) >> 16;
+}
+
+static void dwt_decode97_int(DWTContext *s, int32_t *t)
+{
+    int lev;
+    int w       = s->linelen[s->ndeclevels - 1][0];
+    int32_t *line = s->i_linebuf;
+    int32_t *data = t;
+    /* position at index O of line range [0-5,w+5] cf. extend function */
+    line += 5;
+
+    for (lev = 0; lev < s->ndeclevels; lev++) {
+        int lh = s->linelen[lev][0],
+            lv = s->linelen[lev][1],
+            mh = s->mod[lev][0],
+            mv = s->mod[lev][1],
+            lp;
+        int32_t *l;
+        // HOR_SD
+        l = line + mh;
+        for (lp = 0; lp < lv; lp++) {
+            int i, j = 0;
+            // copy with interleaving
+            for (i = mh; i < lh; i += 2, j++)
+                l[i] = data[w * lp + j];
+            for (i = 1 - mh; i < lh; i += 2, j++)
+                l[i] = data[w * lp + j];
+
+            sr_1d97_int(line, mh, mh + lh);
+
+            for (i = 0; i < lh; i++)
+                data[w * lp + i] = l[i];
+        }
+
+        // VER_SD
+        l = line + mv;
+        for (lp = 0; lp < lh; lp++) {
+            int i, j = 0;
+            // copy with interleaving
+            for (i = mv; i < lv; i += 2, j++)
+                l[i] = data[w * j + lp];
+            for (i = 1 - mv; i < lv; i += 2, j++)
+                l[i] = data[w * j + lp];
+
+            sr_1d97_int(line, mv, mv + lv);
+
+            for (i = 0; i < lv; i++)
+                data[w * i + lp] = l[i];
+        }
+    }
+}
+
+int ff_jpeg2000_dwt_init(DWTContext *s, uint16_t border[2][2],
+                         int decomp_levels, int type)
+{
+    int i, j, lev = decomp_levels, maxlen,
+        b[2][2];
+
+    s->ndeclevels = decomp_levels;
+    s->type       = type;
+
+    for (i = 0; i < 2; i++)
+        for (j = 0; j < 2; j++)
+            b[i][j] = border[i][j];
+
+    maxlen = FFMAX(b[0][1] - b[0][0],
+                   b[1][1] - b[1][0]);
+    while (--lev >= 0)
+        for (i = 0; i < 2; i++) {
+            s->linelen[lev][i] = b[i][1] - b[i][0];
+            s->mod[lev][i]     = b[i][0] & 1;
+            for (j = 0; j < 2; j++)
+                b[i][j] = (b[i][j] + 1) >> 1;
+        }
+    switch (type) {
+    case FF_DWT97:
+        s->f_linebuf = av_malloc((maxlen + 12) * sizeof(*s->f_linebuf));
+        if (!s->f_linebuf)
+            return AVERROR(ENOMEM);
+        break;
+     case FF_DWT97_INT:
+        s->i_linebuf = av_malloc((maxlen + 12) * sizeof(*s->i_linebuf));
+        if (!s->i_linebuf)
+            return AVERROR(ENOMEM);
+        break;
+    case FF_DWT53:
+        s->i_linebuf = av_malloc((maxlen +  6) * sizeof(*s->i_linebuf));
+        if (!s->i_linebuf)
+            return AVERROR(ENOMEM);
+        break;
+    default:
+        return -1;
+    }
+    return 0;
+}
+
+int ff_dwt_decode(DWTContext *s, void *t)
+{
+    switch (s->type) {
+    case FF_DWT97:
+        dwt_decode97_float(s, t);
+        break;
+    case FF_DWT97_INT:
+        dwt_decode97_int(s, t);
+        break;
+    case FF_DWT53:
+        dwt_decode53(s, t);
+        break;
+    default:
+        return -1;
+    }
+    return 0;
+}
+
+void ff_dwt_destroy(DWTContext *s)
+{
+    av_freep(&s->f_linebuf);
+    av_freep(&s->i_linebuf);
+}

libavcodec/jpeg2000dwt.h

+/*
+ * Discrete wavelet transform
+ * Copyright (c) 2007 Kamil Nowosad
+ *
+ * This file is part of Libav.
+ *
+ * Libav is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * Libav is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with Libav; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#ifndef AVCODEC_JPEG2000DWT_H
+#define AVCODEC_JPEG2000DWT_H
+
+/**
+ * @file
+ * Discrete wavelet transform
+ */
+
+#include <stdint.h>
+
+#define FF_DWT_MAX_DECLVLS 32 ///< max number of decomposition levels
+
+enum DWTType {
+    FF_DWT97,
+    FF_DWT53,
+    FF_DWT97_INT
+};
+
+typedef struct DWTContext {
+    /// line lengths { horizontal, vertical } in consecutive decomposition levels
+    uint16_t linelen[FF_DWT_MAX_DECLVLS][2];
+    uint8_t mod[FF_DWT_MAX_DECLVLS][2];  ///< coordinates (x0, y0) of decomp. levels mod 2
+    uint8_t ndeclevels;                  ///< number of decomposition levels
+    uint8_t type;                        ///< 0 for 9/7; 1 for 5/3
+    int32_t *i_linebuf;                  ///< int buffer used by transform
+    float   *f_linebuf;                  ///< float buffer used by transform
+} DWTContext;
+
+/**
+ * Initialize DWT.
+ * @param s                 DWT context
+ * @param border            coordinates of transformed region {{x0, x1}, {y0, y1}}
+ * @param decomp_levels     number of decomposition levels
+ * @param type              0 for DWT 9/7; 1 for DWT 5/3
+ */
+int ff_jpeg2000_dwt_init(DWTContext *s, uint16_t border[2][2],
+                         int decomp_levels, int type);
+
+int ff_dwt_decode(DWTContext *s, void *t);
+
+void ff_dwt_destroy(DWTContext *s);
+
+#endif /* AVCODEC_JPEG2000DWT_H */

libavcodec/mqc.c

+/*
+ * MQ-coder encoder and decoder common functions
+ * Copyright (c) 2007 Kamil Nowosad
+ *
+ * This file is part of Libav.
+ *
+ * Libav is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * Libav is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with Libav; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+/**
+ * MQ-coder common (decoder/encoder) functions
+ * @file
+ * @author Kamil Nowosad
+ */
+
+#include <string.h>
+#include <stdint.h>
+
+#include "mqc.h"
+
+/* MQ coder context state structure */
+typedef struct MqcCxState {
+    uint16_t qe;
+    uint8_t  nmps;
+    uint8_t  nlps;
+    uint8_t  sw;
+} MqcCxState;
+
+static const MqcCxState cx_states[47] = {
+    { 0x5601,  1,  1, 1 },
+    { 0x3401,  2,  6, 0 },
+    { 0x1801,  3,  9, 0 },
+    { 0x0AC1,  4, 12, 0 },
+    { 0x0521,  5, 29, 0 },
+    { 0x0221, 38, 33, 0 },
+    { 0x5601,  7,  6, 1 },
+    { 0x5401,  8, 14, 0 },
+    { 0x4801,  9, 14, 0 },
+    { 0x3801, 10, 14, 0 },
+    { 0x3001, 11, 17, 0 },
+    { 0x2401, 12, 18, 0 },
+    { 0x1C01, 13, 20, 0 },
+    { 0x1601, 29, 21, 0 },
+    { 0x5601, 15, 14, 1 },
+    { 0x5401, 16, 14, 0 },
+    { 0x5101, 17, 15, 0 },
+    { 0x4801, 18, 16, 0 },
+    { 0x3801, 19, 17, 0 },
+    { 0x3401, 20, 18, 0 },
+    { 0x3001, 21, 19, 0 },
+    { 0x2801, 22, 19, 0 },
+    { 0x2401, 23, 20, 0 },
+    { 0x2201, 24, 21, 0 },
+    { 0x1C01, 25, 22, 0 },
+    { 0x1801, 26, 23, 0 },
+    { 0x1601, 27, 24, 0 },
+    { 0x1401, 28, 25, 0 },
+    { 0x1201, 29, 26, 0 },
+    { 0x1101, 30, 27, 0 },
+    { 0x0AC1, 31, 28, 0 },
+    { 0x09C1, 32, 29, 0 },
+    { 0x08A1, 33, 30, 0 },
+    { 0x0521, 34, 31, 0 },
+    { 0x0441, 35, 32, 0 },
+    { 0x02A1, 36, 33, 0 },
+    { 0x0221, 37, 34, 0 },
+    { 0x0141, 38, 35, 0 },
+    { 0x0111, 39, 36, 0 },
+    { 0x0085, 40, 37, 0 },
+    { 0x0049, 41, 38, 0 },
+    { 0x0025, 42, 39, 0 },
+    { 0x0015, 43, 40, 0 },
+    { 0x0009, 44, 41, 0 },
+    { 0x0005, 45, 42, 0 },
+    { 0x0001, 45, 43, 0 },
+    { 0x5601, 46, 46, 0 }
+};
+
+uint16_t ff_mqc_qe [2 * 47];
+uint8_t ff_mqc_nlps[2 * 47];
+uint8_t ff_mqc_nmps[2 * 47];
+
+void ff_mqc_init_contexts(MqcState *mqc)
+{
+    int i;
+    memset(mqc->cx_states, 0, sizeof(mqc->cx_states));
+    mqc->cx_states[MQC_CX_UNI] = 2 * 46;
+    mqc->cx_states[MQC_CX_RL]  = 2 * 3;
+    mqc->cx_states[0]          = 2 * 4;
+
+    for (i = 0; i < 47; i++) {
+        ff_mqc_qe[2 * i]     =
+        ff_mqc_qe[2 * i + 1] = cx_states[i].qe;
+
+        ff_mqc_nlps[2 * i]     = 2 * cx_states[i].nlps + cx_states[i].sw;
+        ff_mqc_nlps[2 * i + 1] = 2 * cx_states[i].nlps + 1 - cx_states[i].sw;
+        ff_mqc_nmps[2 * i]     = 2 * cx_states[i].nmps;
+        ff_mqc_nmps[2 * i + 1] = 2 * cx_states[i].nmps + 1;
+    }
+}

libavcodec/mqc.h

+/*
+ * MQ-coder: structures, common and decoder functions
+ * Copyright (c) 2007 Kamil Nowosad
+ *
+ * This file is part of Libav.
+ *
+ * Libav is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * Libav is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with Libav; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#ifndef AVCODEC_MQC_H
+#define AVCODEC_MQC_H
+
+/**
+ * MQ-coder
+ * @file
+ * @author Kamil Nowosad
+ */
+
+#include <stdint.h>
+
+#define MQC_CX_UNI 17
+#define MQC_CX_RL  18
+
+extern uint16_t ff_mqc_qe[2 * 47];
+extern uint8_t  ff_mqc_nlps[2 * 47];
+extern uint8_t  ff_mqc_nmps[2 * 47];
+
+typedef struct MqcState {
+    uint8_t *bp, *bpstart;
+    unsigned int a;
+    unsigned int c;
+    unsigned int ct;
+    uint8_t cx_states[19];
+} MqcState;
+
+/* decoder */
+
+/**
+ * Initialize MQ-decoder.
+ * @param mqc   MQ decoder state
+ * @param bp    byte poiter
+ */
+void ff_mqc_initdec(MqcState *mqc, uint8_t *bp);
+
+/**
+ * MQ decoder.
+ * @param mqc       MQ decoder state
+ * @param cxstate   Context
+ * @return          Decision (0 ot 1)
+ */
+int ff_mqc_decode(MqcState *mqc, uint8_t *cxstate);
+
+/* common */
+
+/**
+ * MQ-coder context initialisations.
+ * @param mqc       MQ-coder context
+ */
+void ff_mqc_init_contexts(MqcState *mqc);
+
+#endif /* AVCODEC_MQC_H */

libavcodec/mqcdec.c

+/*
+ * MQ-coder decoder
+ * Copyright (c) 2007 Kamil Nowosad
+ *
+ * This file is part of Libav.
+ *
+ * Libav is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * Libav is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with Libav; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+/**
+ * MQ-coder decoder
+ * @file
+ * @author Kamil Nowosad
+ */
+
+#include "mqc.h"
+
+static void bytein(MqcState *mqc)
+{
+    if (*mqc->bp == 0xff) {
+        if (*(mqc->bp + 1) > 0x8f)
+            mqc->c++;
+        else {
+            mqc->bp++;
+            mqc->c += 2 + 0xfe00 - (*mqc->bp << 9);
+        }
+    } else {
+        mqc->bp++;
+        mqc->c += 1 + 0xff00 - (*mqc->bp << 8);
+    }
+}
+
+static int exchange(MqcState *mqc, uint8_t *cxstate, int lps)
+{
+    int d;
+    if ((mqc->a < ff_mqc_qe[*cxstate]) ^ (!lps)) {
+        if (lps)
+            mqc->a = ff_mqc_qe[*cxstate];
+        d = *cxstate & 1;
+        *cxstate = ff_mqc_nmps[*cxstate];
+    } else {
+        if (lps)
+            mqc->a = ff_mqc_qe[*cxstate];
+        d = 1 - (*cxstate & 1);
+        *cxstate = ff_mqc_nlps[*cxstate];
+    }
+    // do RENORMD: see ISO/IEC 15444-1:2002 §C.3.3
+    do {
+        if (!(mqc->c & 0xff)) {
+            mqc->c -= 0x100;
+            bytein(mqc);
+        }
+        mqc->a += mqc->a;
+        mqc->c += mqc->c;
+    } while (!(mqc->a & 0x8000));
+    return d;
+}
+
+void ff_mqc_initdec(MqcState *mqc, uint8_t *bp)
+{
+    ff_mqc_init_contexts(mqc);
+    mqc->bp = bp;
+    mqc->c  = (*mqc->bp ^ 0xff) << 16;
+    bytein(mqc);
+    mqc->c = mqc->c << 7;
+    mqc->a = 0x8000;
+}
+
+int ff_mqc_decode(MqcState *mqc, uint8_t *cxstate)
+{
+    mqc->a -= ff_mqc_qe[*cxstate];
+    if ((mqc->c >> 16) < mqc->a) {
+        if (mqc->a & 0x8000)
+            return *cxstate & 1;
+        else
+            return exchange(mqc, cxstate, 0);
+    } else {
+        mqc->c -= mqc->a << 16;
+        return exchange(mqc, cxstate, 1);
+    }
+}

tests/fate/video.mak

@@ -151,6 +151,9 @@ fate-interplay-mve-8bit: CMD = framecrc -i $(SAMPLES)/interplay-mve/interplay-lo
 FATE_SAMPLES_AVCONV-$(call DEMDEC, IPMOVIE, INTERPLAY_VIDEO) += fate-interplay-mve-16bit
 fate-interplay-mve-16bit: CMD = framecrc -i $(SAMPLES)/interplay-mve/descent3-level5-16bit-partial.mve -pix_fmt rgb24 -an
 
+FATE_SAMPLES_AVCONV-$(call DEMDEC, MXF, JPEG2000) += fate-jpeg2000-dcinema
+fate-jpeg2000-dcinema: CMD = framecrc -flags +bitexact -i $(SAMPLES)/jpeg2000/chiens_dcinema2K.mxf
+
 FATE_SAMPLES_AVCONV-$(call DEMDEC, AVI, KGV1) += fate-kgv1
 fate-kgv1: CMD = framecrc -i $(SAMPLES)/kega/kgv1.avi -pix_fmt rgb555le -an
 

tests/ref/fate/jpeg2000-dcinema

+#tb 0: 1/24
+0,          0,          0,        1, 12441600, 0xf0de508b
+0,          1,          1,        1, 12441600, 0x8e50c249