pjpeg.c

Go to the documentation of this file.

/* Copyright (C) 2013-2016, The Regents of The University of Michigan.
All rights reserved.

This software was developed in the APRIL Robotics Lab under the
direction of Edwin Olson, ebolson@umich.edu. This software may be
available under alternative licensing terms; contact the address above.

This library 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.

This library 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 this library; if not, see <http://www.gnu.org/licenses/>.
*/

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <stdint.h>
#include <string.h>

#include "pjpeg.h"

#include "image_u8.h"
#include "image_u8x3.h"

// https://www.w3.org/Graphics/JPEG/itu-t81.pdf

void pjpeg_idct_2D_double(int32_t in[64], uint8_t *out, uint32_t outstride);
void pjpeg_idct_2D_u32(int32_t in[64], uint8_t *out, uint32_t outstride);
void pjpeg_idct_2D_nanojpeg(int32_t in[64], uint8_t *out, uint32_t outstride);

struct pjpeg_huffman_code
{
    uint8_t nbits;  // how many bits should we actually consume?
    uint8_t code;   // what is the symbol that was encoded? (not actually a DCT coefficient; see encoding)
};

struct pjpeg_decode_state
{
    int error;

    uint32_t width, height;
    uint8_t *in;
    uint32_t inlen;

    uint32_t flags;

    // to decode, we load the next 16 bits of input (generally more
    // than we need). We then look up in our code book how many bits
    // we have actually consumed. For example, if there was a code
    // whose bit sequence was "0", the first 32768 entries would all
    // be copies of {.bits=1, .value=XX}; no matter what the following
    // 15 bits are, we would get the correct decode.
    //
    // Can be up to 8 tables; computed as (ACDC * 2 + htidx)
    struct pjpeg_huffman_code huff_codes[4][65536];
    int huff_codes_present[4];

    uint8_t  qtab[4][64];

    int ncomponents;
    pjpeg_component_t *components;

    int reset_interval;
    int reset_count;
    int reset_next; // What reset marker do we expect next? (add 0xd0)

    int debug;
};

// from K.3.3.1 (page 158)
static uint8_t mjpeg_dht[] = { // header
    0xFF,0xC4,0x01,0xA2,

    // luminance dc coefficients.
    // DC table 0
    0x00,
    // code lengths
    0x00,0x01,0x05,0x01,0x01,0x01,0x01,0x01,0x01,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
    // values
    0x00,0x01,0x02,0x03,0x04,0x05,0x06,0x07,0x08,0x09,0x0A,0x0B,

    // chrominance DC coefficents
    // DC table 1
    0x01,
    // code lengths
    0x00,0x03,0x01,0x01,0x01,0x01,0x01,0x01,0x01,0x01,0x01,0x00,0x00,0x00,0x00,0x00,
    // values
    0x00,0x01,0x02,0x03,0x04,0x05,0x06,0x07,0x08,0x09,0x0A,0x0B,

    // luminance AC coefficients
    // AC table 0
    0x10,
    // code lengths
    0x00,0x02,0x01,0x03,0x03,0x02,0x04,0x03,0x05,0x05,0x04,0x04,0x00,0x00,0x01,0x7D,
    // codes
    0x01,0x02,0x03,0x00,0x04,0x11,0x05,0x12,0x21,0x31,0x41,0x06,0x13,0x51,0x61,
    0x07,0x22,0x71,0x14,0x32,0x81,0x91,0xA1,0x08,0x23,0x42,0xB1,0xC1,0x15,0x52,0xD1,0xF0,0x24,
    0x33,0x62,0x72,0x82,0x09,0x0A,0x16,0x17,0x18,0x19,0x1A,0x25,0x26,0x27,0x28,0x29,0x2A,0x34,
    0x35,0x36,0x37,0x38,0x39,0x3A,0x43,0x44,0x45,0x46,0x47,0x48,0x49,0x4A,0x53,0x54,0x55,0x56,
    0x57,0x58,0x59,0x5A,0x63,0x64,0x65,0x66,0x67,0x68,0x69,0x6A,0x73,0x74,0x75,0x76,0x77,0x78,
    0x79,0x7A,0x83,0x84,0x85,0x86,0x87,0x88,0x89,0x8A,0x92,0x93,0x94,0x95,0x96,0x97,0x98,0x99,
    0x9A,0xA2,0xA3,0xA4,0xA5,0xA6,0xA7,0xA8,0xA9,0xAA,0xB2,0xB3,0xB4,0xB5,0xB6,0xB7,0xB8,0xB9,
    0xBA,0xC2,0xC3,0xC4,0xC5,0xC6,0xC7,0xC8,0xC9,0xCA,0xD2,0xD3,0xD4,0xD5,0xD6,0xD7,0xD8,0xD9,
    0xDA,0xE1,0xE2,0xE3,0xE4,0xE5,0xE6,0xE7,0xE8,0xE9,0xEA,0xF1,0xF2,0xF3,0xF4,0xF5,0xF6,0xF7,
    0xF8,0xF9,0xFA,

    // chrominance DC coefficients
    // DC table 1
    0x11,
    // code lengths
    0x00,0x02,0x01,0x02,0x04,0x04,0x03,0x04,0x07,0x05,0x04,0x04,0x00,0x01,0x02,0x77,
    // values
    0x00,0x01,0x02,0x03,0x11,0x04,0x05,0x21,0x31,0x06,0x12,0x41,0x51,0x07,0x61,0x71,
    0x13,0x22,0x32,0x81,0x08,0x14,0x42,0x91,0xA1,0xB1,0xC1,0x09,0x23,0x33,0x52,0xF0,0x15,0x62,
    0x72,0xD1,0x0A,0x16,0x24,0x34,0xE1,0x25,0xF1,0x17,0x18,0x19,0x1A,0x26,0x27,0x28,0x29,0x2A,
    0x35,0x36,0x37,0x38,0x39,0x3A,0x43,0x44,0x45,0x46,0x47,0x48,0x49,0x4A,0x53,0x54,0x55,0x56,
    0x57,0x58,0x59,0x5A,0x63,0x64,0x65,0x66,0x67,0x68,0x69,0x6A,0x73,0x74,0x75,0x76,0x77,0x78,
    0x79,0x7A,0x82,0x83,0x84,0x85,0x86,0x87,0x88,0x89,0x8A,0x92,0x93,0x94,0x95,0x96,0x97,0x98,
    0x99,0x9A,0xA2,0xA3,0xA4,0xA5,0xA6,0xA7,0xA8,0xA9,0xAA,0xB2,0xB3,0xB4,0xB5,0xB6,0xB7,0xB8,
    0xB9,0xBA,0xC2,0xC3,0xC4,0xC5,0xC6,0xC7,0xC8,0xC9,0xCA,0xD2,0xD3,0xD4,0xD5,0xD6,0xD7,0xD8,
    0xD9,0xDA,0xE2,0xE3,0xE4,0xE5,0xE6,0xE7,0xE8,0xE9,0xEA,0xF2,0xF3,0xF4,0xF5,0xF6,0xF7,0xF8,
    0xF9,0xFA
};

static inline uint8_t max_u8(uint8_t a, uint8_t b)
{
    return a > b ? a : b;
}

// order of coefficients in each DC block
static const char ZZ[64] = { 0,   1,  8, 16,  9,  2,  3, 10,
                             17, 24, 32, 25, 18, 11,  4,  5,
                             12, 19, 26, 33, 40, 48, 41, 34,
                             27, 20, 13,  6,  7, 14, 21, 28,
                             35, 42, 49, 56, 57, 50, 43, 36,
                             29, 22, 15, 23, 30, 37, 44, 51,
                             58, 59, 52, 45, 38, 31, 39, 46,
                             53, 60, 61, 54, 47, 55, 62, 63 };



struct bit_decoder
{
    uint8_t *in;
    uint32_t inpos;
    uint32_t inlen;

    uint32_t bits; // the low order bits contain the next nbits_avail bits.

    int nbits_avail; // how many bits in 'bits' (left aligned) are valid?

    int error;
};

// ensure that at least 'nbits' of data is available in the bit decoder.
static inline void bd_ensure(struct bit_decoder *bd, int nbits)
{
    while (bd->nbits_avail < nbits) {

        if (bd->inpos >= bd->inlen) {
            printf("hallucinating 1s!\n");
            // we hit end of stream hallucinate an infinite stream of 1s
            bd->bits = (bd->bits << 8) | 0xff;
            bd->nbits_avail += 8;
            continue;
        }

        uint8_t nextbyte = bd->in[bd->inpos];
        bd->inpos++;

        if (nextbyte == 0xff && bd->inpos < bd->inlen && bd->in[bd->inpos] == 0x00) {
            // a stuffed byte
            nextbyte = 0xff;
            bd->inpos++;
        }

        // it's an ordinary byte
        bd->bits = (bd->bits << 8) | nextbyte;
        bd->nbits_avail += 8;
    }
}

static inline uint32_t bd_peek_bits(struct bit_decoder *bd, int nbits)
{
    bd_ensure(bd, nbits);

    return (bd->bits >> (bd->nbits_avail - nbits)) & ((1 << nbits) - 1);
}

static inline uint32_t bd_consume_bits(struct bit_decoder *bd, int nbits)
{
    assert(nbits < 32);

    bd_ensure(bd, nbits);

    uint32_t v = (bd->bits >> (bd->nbits_avail - nbits)) & ((1 << nbits) - 1);

    bd->nbits_avail -= nbits;

    return v;
}

// discard without regard for byte stuffing!
static inline void bd_discard_bytes(struct bit_decoder *bd, int nbytes)
{
    assert(bd->nbits_avail == 0);
    bd->inpos += nbytes;
}

static inline int bd_has_more(struct bit_decoder *bd)
{
    return bd->nbits_avail > 0 || bd->inpos < bd->inlen;
}

// throw away up to 7 bits of data so that the next data returned
// began on a byte boundary.
static inline void bd_discard_to_byte_boundary(struct bit_decoder *bd)
{
    bd->nbits_avail -= (bd->nbits_avail & 7);
}

static inline uint32_t bd_get_offset(struct bit_decoder *bd)
{
    return bd->inpos - bd->nbits_avail / 8;
}

static int pjpeg_decode_buffer(struct pjpeg_decode_state *pjd)
{
    // XXX TODO Include sanity check that this is actually a JPG

    struct bit_decoder bd;
    memset(&bd, 0, sizeof(struct bit_decoder));
    bd.in = pjd->in;
    bd.inpos = 0;
    bd.inlen = pjd->inlen;

    int marker_sync_skipped = 0;
    int marker_sync_skipped_from_offset = 0;

    while (bd_has_more(&bd)) {

        uint32_t marker_offset = bd_get_offset(&bd);

        // Look for the 0xff that signifies the beginning of a marker
        bd_discard_to_byte_boundary(&bd);

        while (bd_consume_bits(&bd, 8) != 0xff) {
            if (marker_sync_skipped == 0)
                marker_sync_skipped_from_offset = marker_offset;
            marker_sync_skipped++;
            continue;
        }

        if (marker_sync_skipped) {
            printf("%08x: skipped %04x bytes\n", marker_sync_skipped_from_offset, marker_sync_skipped);
            marker_sync_skipped = 0;
        }

        uint8_t marker = bd_consume_bits(&bd, 8);

//        printf("marker %08x : %02x\n", marker_offset, marker);

        switch (marker) {

            case 0xd8:  // start of image. Great, continue.
                continue;

                // below are the markers that A) we don't care about
                // that B) encode length as two bytes.
                //
                // Note: Other unknown fields should not be added since
                // we should be able to skip over them by looking for
                // the next marker byte.
            case 0xe0: // JFIF header.
            case 0xe1: // EXIF header (Yuck: Payload may contain 0xff 0xff!)
            case 0xe2: // ICC Profile. (Yuck: payload may contain 0xff 0xff!)
            case 0xe6: // some other common header
            case 0xfe: // Comment
            {
                uint16_t length = bd_consume_bits(&bd, 16);
                bd_discard_bytes(&bd, length - 2);
                continue;
            }

            case 0xdb: { // DQT Define Quantization Table
                uint16_t length = bd_consume_bits(&bd, 16);

                if (((length-2) % 65) != 0)
                    return PJPEG_ERR_DQT;

                // can contain multiple DQTs
                for (int offset = 0; offset < length - 2; offset += 65) {

                    // pq: quant table element precision. 0=8bit, 1=16bit.
                    // tq: quant table destination id.
                    uint8_t  pqtq = bd_consume_bits(&bd, 8);

                    if ((pqtq & 0xf0) != 0 || (pqtq & 0x0f) >= 4)
                        return PJPEG_ERR_DQT;

                    uint8_t id = pqtq & 3;

                    for (int i = 0; i < 64; i++)
                        pjd->qtab[id][i] = bd_consume_bits(&bd, 8);
                }

                break;
            }

            case 0xc0: { // SOF, non-differential, huffman, baseline
                uint16_t length = bd_consume_bits(&bd, 16);
                (void) length;

                uint8_t p = bd_consume_bits(&bd, 8); // precision
                if (p != 8)
                    return PJPEG_ERR_SOF;

                pjd->height = bd_consume_bits(&bd, 16);
                pjd->width = bd_consume_bits(&bd, 16);

//                printf("%d x %d\n", pjd->height, pjd->width);

                int nf = bd_consume_bits(&bd, 8); // # image components

                if (nf < 1 || nf > 3)
                    return PJPEG_ERR_SOF;

                pjd->ncomponents = nf;
                pjd->components = calloc(nf, sizeof(struct pjpeg_component));

                for (int i = 0; i < nf; i++) {
                    // comp. identifier
                    pjd->components[i].id = bd_consume_bits(&bd, 8);

                    // horiz/vert sampling
                    pjd->components[i].hv = bd_consume_bits(&bd, 8);
                    pjd->components[i].scaley = pjd->components[i].hv & 0x0f;
                    pjd->components[i].scalex = pjd->components[i].hv >> 4;

                    // which quant table?
                    pjd->components[i].tq = bd_consume_bits(&bd, 8);
                }
                break;
            }

            case 0xc1: // SOF, non-differential, huffman,    extended DCT
            case 0xc2: // SOF, non-differential, huffman,    progressive DCT
            case 0xc3: // SOF, non-differential, huffman,    lossless
            case 0xc5: // SOF, differential,     huffman,    baseline DCT
            case 0xc6: // SOF, differential,     huffman,    progressive
            case 0xc7: // SOF, differential,     huffman,    lossless
            case 0xc8: // reserved
            case 0xc9: // SOF, non-differential, arithmetic, extended
            case 0xca: // SOF, non-differential, arithmetic, progressive
            case 0xcb: // SOF, non-differential, arithmetic, lossless
            case 0xcd: // SOF, differential,     arithmetic, sequential
            case 0xce: // SOF, differential,     arithmetic, progressive
            case 0xcf: // SOF, differential,     arithmetic, lossless
            {
                printf("pjepg.c: unsupported JPEG type %02x\n", marker);
                return PJEPG_ERR_UNSUPPORTED;
            }

            case 0xc4: { // DHT Define Huffman Tables
                // [ED: the encoding of these tables is really quite
                // clever!]
                uint16_t length = bd_consume_bits(&bd, 16);
                length = length - 2;

                while (length > 0) {
                    uint8_t TcTh = bd_consume_bits(&bd, 8);
                    length--;
                    uint8_t Tc = (TcTh >> 4);
                    int Th = TcTh & 0x0f; // which index are we using?

                    if (Tc >= 2 || Th >= 2)
                        // Tc must be either AC=1 or DC=0.
                        // Th must be less than 2
                        return PJPEG_ERR_DHT;

                    int htidx = Tc*2 + Th;

                    uint8_t L[17]; // how many symbols of each bit length?
                    L[0] = 0;      // no 0 bit codes :)
                    for (int nbits = 1; nbits <= 16; nbits++) {
                        L[nbits] = bd_consume_bits(&bd, 8);
                        length -= L[nbits];
                    }
                    length -= 16;

                    uint32_t code_pos = 0;

                    for (int nbits = 1; nbits <= 16; nbits++) {
                        int nvalues = L[nbits];

                        // how many entries will we fill?
                        // (a 1 bit code will fill 32768, a 2 bit code 16384, ...)
                        uint32_t ncodes = (1 << (16 - nbits));

                        // consume the values...
                        for (int vi = 0; vi < nvalues; vi++) {
                            uint8_t code = bd_consume_bits(&bd, 8);

                            if (code_pos + ncodes > 0xffff)
                                return PJPEG_ERR_DHT;

                            for (int ci = 0; ci < ncodes; ci++) {
                                pjd->huff_codes[htidx][code_pos].nbits = nbits;
                                pjd->huff_codes[htidx][code_pos].code = code;
                                code_pos++;
                            }
                        }
                    }
                    pjd->huff_codes_present[htidx] = 1;
                }
                break;
            }

                // a sequentially-encoded JPG has one SOS segment. A
                // progressive JPG will have multiple SOS segments.
            case 0xda: { // Start Of Scan (SOS)

                // Note that this marker frame (and its encoded
                // length) does NOT include the bitstream that
                // follows.

                uint16_t length = bd_consume_bits(&bd, 16);
                (void) length;

                // number of components in this scan
                uint8_t ns = bd_consume_bits(&bd, 8);

                // for each component, what is the index into our pjd->components[] array?
                uint8_t comp_idx[ns];
                memset(comp_idx, 0, ns*sizeof(uint8_t));

                for (int i = 0; i < ns; i++) {
                    // component name
                    uint8_t cs = bd_consume_bits(&bd, 8);

                    int found = 0;
                    for (int j = 0; j < pjd->ncomponents; j++) {

                        if (cs == pjd->components[j].id) {
                            // which huff tables will we use for
                            // DC (high 4 bits) and AC (low 4 bits)
                            pjd->components[j].tda = bd_consume_bits(&bd, 8);
                            comp_idx[i] = j;
                            found = 1;
                            break;
                        }
                    }

                    if (!found)
                        return PJPEG_ERR_SOS;
                }

                // start of spectral selection. baseline == 0
                uint8_t ss = bd_consume_bits(&bd, 8);

                // end of spectral selection. baseline == 0x3f
                uint8_t se = bd_consume_bits(&bd, 8);

                // successive approximation bits. baseline == 0
                uint8_t Ahl = bd_consume_bits(&bd, 8);

                if (ss != 0 || se != 0x3f || Ahl != 0x00)
                    return PJPEG_ERR_SOS;

                // compute the dimensions of each MCU in pixels
                int maxmcux = 0, maxmcuy = 0;
                for (int i = 0; i < ns; i++) {
                    struct pjpeg_component *comp = &pjd->components[comp_idx[i]];

                    maxmcux = max_u8(maxmcux, comp->scalex * 8);
                    maxmcuy = max_u8(maxmcuy, comp->scaley * 8);
                }

                // how many MCU blocks are required to encode the whole image?
                int mcus_x = (pjd->width + maxmcux - 1) / maxmcux;
                int mcus_y = (pjd->height + maxmcuy - 1) / maxmcuy;

                if (0)
                    printf("Image has %d x %d MCU blocks, each %d x %d pixels\n",
                           mcus_x, mcus_y, maxmcux, maxmcuy);

                // allocate output storage
                for (int i = 0; i < ns; i++) {
                    struct pjpeg_component *comp = &pjd->components[comp_idx[i]];
                    comp->width = mcus_x * comp->scalex * 8;
                    comp->height = mcus_y * comp->scaley * 8;
                    comp->stride = comp->width;

                    int alignment = 32;
                    if ((comp->stride % alignment) != 0)
                        comp->stride += alignment - (comp->stride % alignment);

                    comp->data = calloc(comp->height * comp->stride, 1);
                }


                // each component has its own DC prediction
                int32_t dcpred[ns];
                memset(dcpred, 0, sizeof(dcpred));

                pjd->reset_count = 0;

                for (int mcu_y = 0; mcu_y < mcus_y; mcu_y++) {
                    for (int mcu_x = 0; mcu_x < mcus_x; mcu_x++) {

                        // the next two bytes in the input stream
                        // should be 0xff 0xdN, where N is the next
                        // reset counter.
                        //
                        // Our bit decoder may have already shifted
                        // these into the buffer.  Consequently, we
                        // want to use our bit decoding functions to
                        // check for the marker. But we must first
                        // discard any fractional bits left.
                        if (pjd->reset_interval > 0 && pjd->reset_count == pjd->reset_interval) {

                            // RST markers are byte-aligned, so force
                            // the bit-decoder to the next byte
                            // boundary.
                            bd_discard_to_byte_boundary(&bd);

                            while (1) {
                                int32_t value = bd_consume_bits(&bd, 8);
                                if (bd.inpos > bd.inlen)
                                    return PJPEG_ERR_EOF;
                                if (value == 0xff)
                                    break;
                                printf("RST SYNC\n");
                            }

                            int32_t marker = bd_consume_bits(&bd, 8);

//                            printf("%04x: RESET? %02x\n", *bd.inpos,  marker);
                            if (marker != (0xd0 + pjd->reset_next))
                                return PJPEG_ERR_RESET;

                            pjd->reset_count = 0;
                            pjd->reset_next = (pjd->reset_next + 1) & 0x7;

                            memset(dcpred, 0, sizeof(dcpred));
                        }

                        for (int nsidx = 0; nsidx < ns; nsidx++) {

                            struct pjpeg_component *comp = &pjd->components[comp_idx[nsidx]];

                            int32_t block[64];

                            int qtabidx = comp->tq; // which quant table?

                            for (int sby = 0; sby < comp->scaley; sby++) {
                                for (int sbx = 0; sbx < comp->scalex; sbx++) {
                                    // decode block for component nsidx
                                    memset(block, 0, sizeof(block));

                                    int dc_huff_table_idx = comp->tda >> 4;
                                    int ac_huff_table_idx = 2 + (comp->tda & 0x0f);

                                    if (!pjd->huff_codes_present[dc_huff_table_idx] ||
                                        !pjd->huff_codes_present[ac_huff_table_idx])
                                        return PJPEG_ERR_MISSING_DHT; // probably an MJPEG.


                                    if (1) {
                                        // do DC coefficient
                                        uint32_t next16 = bd_peek_bits(&bd, 16);
                                        struct pjpeg_huffman_code *huff_code = &pjd->huff_codes[dc_huff_table_idx][next16];
                                        bd_consume_bits(&bd, huff_code->nbits);

                                        int ssss = huff_code->code & 0x0f; // ssss == number of additional bits to read
                                        int32_t value = bd_consume_bits(&bd, ssss);

                                        // if high bit is clear, it's negative
                                        if ((value & (1 << (ssss-1))) == 0)
                                            value += ((-1) << ssss) + 1;

                                        dcpred[nsidx] += value;
                                        block[0] = dcpred[nsidx] * pjd->qtab[qtabidx][0];
                                    }

                                    if (1) {
                                        // do AC coefficients
                                        for (int coeff = 1; coeff < 64; coeff++) {

                                            uint32_t next16 = bd_peek_bits(&bd, 16);

                                            struct pjpeg_huffman_code *huff_code = &pjd->huff_codes[ac_huff_table_idx][next16];
                                            bd_consume_bits(&bd, huff_code->nbits);

                                            if (huff_code->code == 0) {
                                                break; // EOB
                                            }

                                            int rrrr = huff_code->code >> 4; // run length of zeros
                                            int ssss = huff_code->code & 0x0f;

                                            int32_t value = bd_consume_bits(&bd, ssss);

                                            // if high bit is clear, it's negative
                                            if ((value & (1 << (ssss-1))) == 0)
                                                value += ((-1) << ssss) + 1;

                                            coeff += rrrr;

                                            block[(int) ZZ[coeff]] = value * pjd->qtab[qtabidx][coeff];
                                        }
                                    }

                                    // do IDCT

                                    // output block's upper-left
                                    // coordinate (in pixels) is
                                    // (comp_x, comp_y).
                                    uint32_t comp_x = (mcu_x * comp->scalex + sbx) * 8;
                                    uint32_t comp_y = (mcu_y * comp->scaley + sby) * 8;
                                    uint32_t dataidx = comp_y * comp->stride + comp_x;

//                                    pjpeg_idct_2D_u32(block, &comp->data[dataidx], comp->stride);
                                    pjpeg_idct_2D_nanojpeg(block, &comp->data[dataidx], comp->stride);
                                }
                            }
                        }

                        pjd->reset_count++;
//                        printf("%04x: reset count %d / %d\n", pjd->inpos, pjd->reset_count, pjd->reset_interval);

                    }
                }

                break;
            }

            case 0xd9: { // EOI End of Image
                goto got_end_of_image;
            }

            case 0xdd: { // Define Restart Interval
                uint16_t length = bd_consume_bits(&bd, 16);
                if (length != 4)
                    return PJPEG_ERR_DRI;

                // reset interval measured in the number of MCUs
                pjd->reset_interval = bd_consume_bits(&bd, 16);

                break;
            }

            default: {
                printf("pjepg: Unknown marker %02x at offset %04x\n", marker, marker_offset);

                // try to skip it.
                uint16_t length = bd_consume_bits(&bd, 16);
                bd_discard_bytes(&bd, length - 2);
                continue;
            }
        } // switch (marker)
    } // while inpos < inlen

  got_end_of_image:

    return PJPEG_OKAY;
}

void pjpeg_destroy(pjpeg_t *pj)
{
    if (!pj)
        return;

    for (int i = 0; i < pj->ncomponents; i++)
        free(pj->components[i].data);
    free(pj->components);

    free(pj);
}


// just grab the first component.
image_u8_t *pjpeg_to_u8_baseline(pjpeg_t *pj)
{
    assert(pj->ncomponents > 0);

    pjpeg_component_t *comp = &pj->components[0];

    assert(comp->width >= pj->width && comp->height >= pj->height);

    image_u8_t *im = image_u8_create(pj->width, pj->height);
    for (int y = 0; y < im->height; y++)
        memcpy(&im->buf[y*im->stride], &comp->data[y*comp->stride], pj->width);

    return im;
}

static inline uint8_t clampd(double v)
{
    if (v < 0)
        return 0;
    if (v > 255)
        return 255;

    return (uint8_t) v;
}

static inline uint8_t clamp_u8(int32_t v)
{
    if (v < 0)
        return 0;
    if (v > 255)
        return 255;
    return v;
}

// color conversion formulas taken from JFIF spec v 1.02
image_u8x3_t *pjpeg_to_u8x3_baseline(pjpeg_t *pj)
{
    assert(pj->ncomponents == 3);

    pjpeg_component_t *Y = &pj->components[0];
    pjpeg_component_t *Cb = &pj->components[1];
    pjpeg_component_t *Cr = &pj->components[2];

    int Cb_factor_y = Y->height / Cb->height;
    int Cb_factor_x = Y->width / Cb->width;

    int Cr_factor_y = Y->height / Cr->height;
    int Cr_factor_x = Y->width / Cr->width;

    image_u8x3_t *im = image_u8x3_create(pj->width, pj->height);

    if (Cr_factor_y == 1 && Cr_factor_x == 1 && Cb_factor_y == 1 && Cb_factor_x == 1) {

        for (int y = 0; y < pj->height; y++) {
            for (int x = 0; x < pj->width; x++) {
                int32_t y_val  = Y->data[y*Y->stride + x] * 65536;
                int32_t cb_val = Cb->data[y*Cb->stride + x] - 128;
                int32_t cr_val = Cr->data[y*Cr->stride + x] - 128;

                int32_t r_val = y_val +  91881 * cr_val;
                int32_t g_val = y_val + -22554 * cb_val - 46802 * cr_val;
                int32_t b_val = y_val + 116130 * cb_val;

                im->buf[y*im->stride + 3*x + 0 ] = clamp_u8(r_val >> 16);
                im->buf[y*im->stride + 3*x + 1 ] = clamp_u8(g_val >> 16);
                im->buf[y*im->stride + 3*x + 2 ] = clamp_u8(b_val >> 16);
            }
        }
    } else if (Cb_factor_y == Cr_factor_y && Cb_factor_x == Cr_factor_x) {
        for (int by = 0; by < pj->height / Cb_factor_y; by++) {
            for (int bx = 0; bx < pj->width / Cb_factor_x; bx++) {

                int32_t cb_val = Cb->data[by*Cb->stride + bx] - 128;
                int32_t cr_val = Cr->data[by*Cr->stride + bx] - 128;

                int32_t r0 =  91881 * cr_val;
                int32_t g0 = -22554 * cb_val - 46802 * cr_val;
                int32_t b0 = 116130 * cb_val;

                for (int dy = 0; dy < Cb_factor_y; dy++) {
                    int y = by*Cb_factor_y + dy;

                    for (int dx = 0; dx < Cb_factor_x; dx++) {
                        int x = bx*Cb_factor_x + dx;

                        int32_t y_val = Y->data[y*Y->stride + x] * 65536;

                        int32_t r_val = r0 + y_val;
                        int32_t g_val = g0 + y_val;
                        int32_t b_val = b0 + y_val;

                        im->buf[y*im->stride + 3*x + 0 ] = clamp_u8(r_val >> 16);
                        im->buf[y*im->stride + 3*x + 1 ] = clamp_u8(g_val >> 16);
                        im->buf[y*im->stride + 3*x + 2 ] = clamp_u8(b_val >> 16);
                    }
                }
            }
        }
    } else {

        for (int y = 0; y < pj->height; y++) {
            for (int x = 0; x < pj->width; x++) {
                int32_t y_val  = Y->data[y*Y->stride + x];
                int32_t cb_val = Cb->data[(y / Cb_factor_y)*Cb->stride + (x / Cb_factor_x)] - 128;
                int32_t cr_val = Cr->data[(y / Cr_factor_y)*Cr->stride + (x / Cr_factor_x)] - 128;

                uint8_t r_val = clampd(y_val + 1.402 * cr_val);
                uint8_t g_val = clampd(y_val - 0.34414 * cb_val - 0.71414 * cr_val);
                uint8_t b_val = clampd(y_val + 1.772 * cb_val);

                im->buf[y*im->stride + 3*x + 0 ] = r_val;
                im->buf[y*im->stride + 3*x + 1 ] = g_val;
                im->buf[y*im->stride + 3*x + 2 ] = b_val;
            }
        }
    }

    return im;
}

// returns NULL if file loading fails.
pjpeg_t *pjpeg_create_from_file(const char *path, uint32_t flags, int *error)
{
    FILE *f = fopen(path, "r");
    if (f == NULL)
        return NULL;

    fseek(f, 0, SEEK_END);
    long buflen = ftell(f);

    uint8_t *buf = malloc(buflen);
    fseek(f, 0, SEEK_SET);
    int res = fread(buf, 1, buflen, f);
    fclose(f);
    if (res != buflen) {
        free(buf);
        if (error)
            *error = PJPEG_ERR_FILE;
        return NULL;
    }

    pjpeg_t *pj = pjpeg_create_from_buffer(buf, buflen, flags, error);

    free(buf);
    return pj;
}

pjpeg_t *pjpeg_create_from_buffer(uint8_t *buf, int buflen, uint32_t flags, int *error)
{
    struct pjpeg_decode_state pjd;
    memset(&pjd, 0, sizeof(pjd));

    if (flags & PJPEG_MJPEG) {
        pjd.in = mjpeg_dht;
        pjd.inlen = sizeof(mjpeg_dht);
        int result = pjpeg_decode_buffer(&pjd);
        assert(result == 0);
    }

    pjd.in = buf;
    pjd.inlen = buflen;
    pjd.flags = flags;

    int result = pjpeg_decode_buffer(&pjd);
    if (error)
        *error = result;

    if (result) {
        for (int i = 0; i < pjd.ncomponents; i++)
            free(pjd.components[i].data);
        free(pjd.components);

        return NULL;
    }

    pjpeg_t *pj = calloc(1, sizeof(pjpeg_t));

    pj->width = pjd.width;
    pj->height = pjd.height;
    pj->ncomponents = pjd.ncomponents;
    pj->components = pjd.components;

    return pj;
}