/* stbi-1.03 - public domain JPEG/PNG reader - http://nothings.org/stb_image.c
                      when you control the images you're loading

   QUICK NOTES:
      Primarily of interest to game developers and other people who can
          avoid problematic images and only need the trivial interface

      JPEG baseline (no JPEG progressive, no oddball channel decimations)
      PNG non-interlaced
      BMP non-1bpp, non-RLE
      TGA (not sure what subset, if a subset)
      HDR (radiance rgbE format)
      writes BMP,TGA (define STBI_NO_WRITE to remove code)
      decoded from memory or through stdio FILE (define STBI_NO_STDIO to remove code)

   TODO:
      stbi_info_*
      PSD loader

   history:
      1.03   bugfixes to STBI_NO_STDIO, STBI_NO_HDR
      1.02   support for (subset of) HDR files, float interface for preferred access to them
      1.01   fix bug: possible bug in handling right-side up bmps... not sure
             fix bug: the stbi_bmp_load() and stbi_tga_load() functions didn't work at all
      1.00   interface to zlib that skips zlib header
      0.99   correct handling of alpha in palette
      0.98   TGA loader by lonesock; dynamically add loaders (untested)
      0.97   jpeg errors on too large a file; also catch another malloc failure
      0.96   fix detection of invalid v value - particleman@mollyrocket forum
      0.95   during header scan, seek to markers in case of padding
      0.94   STBI_NO_STDIO to disable stdio usage; rename all #defines the same
      0.93   handle jpegtran output; verbose errors
      0.92   read 4,8,16,24,32-bit BMP files of several formats
      0.91   output 24-bit Windows 3.0 BMP files
      0.90   fix a few more warnings; bump version number to approach 1.0
      0.61   bugfixes due to Marc LeBlanc, Christopher Lloyd
      0.60   fix compiling as c++
      0.59   fix warnings: merge Dave Moore's -Wall fixes
      0.58   fix bug: zlib uncompressed mode len/nlen was wrong endian
      0.57   fix bug: jpg last huffman symbol before marker was >9 bits but less
                      than 16 available
      0.56   fix bug: zlib uncompressed mode len vs. nlen
      0.55   fix bug: restart_interval not initialized to 0
      0.54   allow NULL for 'int *comp'
      0.53   fix bug in png 3->4; speedup png decoding
      0.52   png handles req_comp=3,4 directly; minor cleanup; jpeg comments
      0.51   obey req_comp requests, 1-component jpegs return as 1-component,
             on 'test' only check type, not whether we support this variant
*/

#include "stb_image_aug.h"

#ifndef STBI_NO_STDIO
#include <stdio.h>
#endif
#include <stdlib.h>
#include <memory.h>
#include <assert.h>
#include <stdarg.h>

#ifndef _MSC_VER
#define __forceinline
#endif

// implementation:
typedef unsigned char uint8;
typedef unsigned short uint16;
typedef   signed short  int16;
typedef unsigned int   uint32;
typedef   signed int    int32;
typedef unsigned int   uint;

// should produce compiler error if size is wrong
typedef unsigned char validate_uint32[sizeof(uint32)==4];

#if defined(STBI_NO_STDIO) && !defined(STBI_NO_WRITE)
#define STBI_NO_WRITE
#endif

#ifndef STBI_NO_DDS
#include "stbi_DDS_aug.h"
#endif

//	I (JLD) want full messages for SOIL
#define STBI_FAILURE_USERMSG 1

//////////////////////////////////////////////////////////////////////////////
//
// Generic API that works on all image types
//

static char *failure_reason;

char *stbi_failure_reason(void)
{
   return failure_reason;
}

static int e(char *str)
{
   failure_reason = str;
   return 0;
}

#ifdef STBI_NO_FAILURE_STRINGS
   #define e(x,y)  0
#elif defined(STBI_FAILURE_USERMSG)
   #define e(x,y)  e(y)
#else
   #define e(x,y)  e(x)
#endif

#define ep(x,y)   (e(x,y)?NULL:NULL)

void stbi_image_free(unsigned char *retval_from_stbi_load)
{
   free(retval_from_stbi_load);
}

#define MAX_LOADERS  32
stbi_loader *loaders[MAX_LOADERS];
static int max_loaders = 0;

int stbi_register_loader(stbi_loader *loader)
{
   int i;
   for (i=0; i < MAX_LOADERS; ++i) {
      // already present?
      if (loaders[i] == loader)
         return 1;
      // end of the list?
      if (loaders[i] == NULL) {
         loaders[i] = loader;
         max_loaders = i+1;
         return 1;
      }
   }
   // no room for it
   return 0;
}

#ifndef STBI_NO_HDR
static float   *ldr_to_hdr(stbi_uc *data, int x, int y, int comp);
static stbi_uc *hdr_to_ldr(float   *data, int x, int y, int comp);
#endif

#ifndef STBI_NO_STDIO
unsigned char *stbi_load(char *filename, int *x, int *y, int *comp, int req_comp)
{
   FILE *f = fopen(filename, "rb");
   unsigned char *result;
   if (!f) return ep("can't fopen", "Unable to open file");
   result = stbi_load_from_file(f,x,y,comp,req_comp);
   fclose(f);
   return result;
}

unsigned char *stbi_load_from_file(FILE *f, int *x, int *y, int *comp, int req_comp)
{
   int i;
   if (stbi_jpeg_test_file(f))
      return stbi_jpeg_load_from_file(f,x,y,comp,req_comp);
   if (stbi_png_test_file(f))
      return stbi_png_load_from_file(f,x,y,comp,req_comp);
   if (stbi_bmp_test_file(f))
      return stbi_bmp_load_from_file(f,x,y,comp,req_comp);
   #ifndef STBI_NO_DDS
   if (stbi_dds_test_file(f))
      return stbi_dds_load_from_file(f,x,y,comp,req_comp);
   #endif
   #ifndef STBI_NO_HDR
   if (stbi_hdr_test_file(f)) {
      float *hdr = stbi_hdr_load_from_file(f, x,y,comp,req_comp);
      return hdr_to_ldr(hdr, *x, *y, req_comp ? req_comp : *comp);
   }
   #endif
   for (i=0; i < max_loaders; ++i)
      if (loaders[i]->test_file(f))
         return loaders[i]->load_from_file(f,x,y,comp,req_comp);
   // test tga last because it's a crappy test!
   if (stbi_tga_test_file(f))
      return stbi_tga_load_from_file(f,x,y,comp,req_comp);
   return ep("unknown image type", "Image not of any known type, or corrupt");
}
#endif

unsigned char *stbi_load_from_memory(stbi_uc *buffer, int len, int *x, int *y, int *comp, int req_comp)
{
   int i;
   if (stbi_jpeg_test_memory(buffer,len))
      return stbi_jpeg_load_from_memory(buffer,len,x,y,comp,req_comp);
   if (stbi_png_test_memory(buffer,len))
      return stbi_png_load_from_memory(buffer,len,x,y,comp,req_comp);
   if (stbi_bmp_test_memory(buffer,len))
      return stbi_bmp_load_from_memory(buffer,len,x,y,comp,req_comp);
   #ifndef STBI_NO_DDS
   if (stbi_dds_test_memory(buffer,len))
      return stbi_dds_load_from_memory(buffer,len,x,y,comp,req_comp);
   #endif
   #ifndef STBI_NO_HDR
   if (stbi_hdr_test_memory(buffer, len)) {
      float *hdr = stbi_hdr_load_from_memory(buffer, len,x,y,comp,req_comp);
      return hdr_to_ldr(hdr, *x, *y, req_comp ? req_comp : *comp);
   }
   #endif
   for (i=0; i < max_loaders; ++i)
      if (loaders[i]->test_memory(buffer,len))
         return loaders[i]->load_from_memory(buffer,len,x,y,comp,req_comp);
   // test tga last because it's a crappy test!
   if (stbi_tga_test_memory(buffer,len))
      return stbi_tga_load_from_memory(buffer,len,x,y,comp,req_comp);
   return ep("unknown image type", "Image not of any known type, or corrupt");
}

#ifndef STBI_NO_HDR

#ifndef STBI_NO_STDIO
float *stbi_loadf(char *filename, int *x, int *y, int *comp, int req_comp)
{
   FILE *f = fopen(filename, "rb");
   float *result;
   if (!f) return ep("can't fopen", "Unable to open file");
   result = stbi_loadf_from_file(f,x,y,comp,req_comp);
   fclose(f);
   return result;
}

float *stbi_loadf_from_file(FILE *f, int *x, int *y, int *comp, int req_comp)
{
   unsigned char *data;
   #ifndef STBI_NO_HDR
   if (stbi_hdr_test_file(f))
      return stbi_hdr_load_from_file(f,x,y,comp,req_comp);
   #endif
   data = stbi_load_from_file(f, x, y, comp, req_comp);
   if (data)
      return ldr_to_hdr(data, *x, *y, req_comp ? req_comp : *comp);
   return ep("unknown image type", "Image not of any known type, or corrupt");
}
#endif

float *stbi_loadf_from_memory(stbi_uc *buffer, int len, int *x, int *y, int *comp, int req_comp)
{
   stbi_uc *data;
   #ifndef STBI_NO_HDR
   if (stbi_hdr_test_memory(buffer, len))
      return stbi_hdr_load_from_memory(buffer, len,x,y,comp,req_comp);
   #endif
   data = stbi_load_from_memory(buffer, len, x, y, comp, req_comp);
   if (data)
      return ldr_to_hdr(data, *x, *y, req_comp ? req_comp : *comp);
   return ep("unknown image type", "Image not of any known type, or corrupt");
}
#endif

// these is-hdr-or-not is defined independent of whether STBI_NO_HDR is
// defined, for API simplicity; if STBI_NO_HDR is defined, it always
// reports false!

extern int      stbi_is_hdr_from_memory(stbi_uc *buffer, int len)
{
   #ifndef STBI_NO_HDR
   return stbi_hdr_test_memory(buffer, len);
   #else
   return 0;
   #endif
}

#ifndef STBI_NO_STDIO
extern int      stbi_is_hdr          (char *filename)
{
   FILE *f = fopen(filename, "rb");
   int result=0;
   if (f) {
      result = stbi_is_hdr_from_file(f);
      fclose(f);
   }
   return result;
}

extern int      stbi_is_hdr_from_file(FILE *f)
{
   #ifndef STBI_NO_HDR
   return stbi_hdr_test_file(f);
   #else
   return 0;
   #endif
}

#endif

// @TODO: get image dimensions & components without fully decoding
#ifndef STBI_NO_STDIO
extern int      stbi_info            (char *filename,           int *x, int *y, int *comp);
extern int      stbi_info_from_file  (FILE *f,                  int *x, int *y, int *comp);
#endif
extern int      stbi_info_from_memory(stbi_uc *buffer, int len, int *x, int *y, int *comp);

#ifndef STBI_NO_HDR
static float h2l_gamma_i=1.0f/2.2f, h2l_scale_i=1.0f;
static float l2h_gamma=2.2f, l2h_scale=1.0f;

void   stbi_hdr_to_ldr_gamma(float gamma) { h2l_gamma_i = 1/gamma; }
void   stbi_hdr_to_ldr_scale(float scale) { h2l_scale_i = 1/scale; }

void   stbi_ldr_to_hdr_gamma(float gamma) { l2h_gamma = gamma; }
void   stbi_ldr_to_hdr_scale(float scale) { l2h_scale = scale; }
#endif


//////////////////////////////////////////////////////////////////////////////
//
// Common code used by all image loaders
//

// image width, height, # components
static uint32 img_x, img_y;
static int img_n, img_out_n;

enum
{
   SCAN_load=0,
   SCAN_type,
   SCAN_header,
};

// An API for reading either from memory or file.
#ifndef STBI_NO_STDIO
static FILE  *img_file;
#endif
static uint8 *img_buffer, *img_buffer_end;

#ifndef STBI_NO_STDIO
static void start_file(FILE *f)
{
   img_file = f;
}
#endif

static void start_mem(uint8 *buffer, int len)
{
#ifndef STBI_NO_STDIO
   img_file = NULL;
#endif
   img_buffer = buffer;
   img_buffer_end = buffer+len;
}

static int get8(void)
{
#ifndef STBI_NO_STDIO
   if (img_file) {
      int c = fgetc(img_file);
      return c == EOF ? 0 : c;
   }
#endif
   if (img_buffer < img_buffer_end)
      return *img_buffer++;
   return 0;
}

static int at_eof(void)
{
#ifndef STBI_NO_STDIO
   if (img_file)
      return feof(img_file);
#endif
   return img_buffer >= img_buffer_end;
}

static uint8 get8u(void)
{
   return (uint8) get8();
}

static void skip(int n)
{
#ifndef STBI_NO_STDIO
   if (img_file)
      fseek(img_file, n, SEEK_CUR);
   else
#endif
      img_buffer += n;
}

static int get16(void)
{
   int z = get8();
   return (z << 8) + get8();
}

static uint32 get32(void)
{
   uint32 z = get16();
   return (z << 16) + get16();
}

static int get16le(void)
{
   int z = get8();
   return z + (get8() << 8);
}

static uint32 get32le(void)
{
   uint32 z = get16le();
   return z + (get16le() << 16);
}

static void getn(stbi_uc *buffer, int n)
{
#ifndef STBI_NO_STDIO
   if (img_file) {
      fread(buffer, 1, n, img_file);
      return;
   }
#endif
   memcpy(buffer, img_buffer, n);
   img_buffer += n;
}

//////////////////////////////////////////////////////////////////////////////
//
//  generic converter from built-in img_n to req_comp
//    individual types do this automatically as much as possible (e.g. jpeg
//    does all cases internally since it needs to colorspace convert anyway,
//    and it never has alpha, so very few cases ). png can automatically
//    interleave an alpha=255 channel, but falls back to this for other cases
//
//  assume data buffer is malloced, so malloc a new one and free that one
//  only failure mode is malloc failing

static uint8 compute_y(int r, int g, int b)
{
   return (uint8) (((r*77) + (g*150) +  (29*b)) >> 8);
}

static unsigned char *convert_format(unsigned char *data, int img_n, int req_comp)
{
   uint i,j;
   unsigned char *good;

   if (req_comp == img_n) return data;
   assert(req_comp >= 1 && req_comp <= 4);

   good = (unsigned char *) malloc(req_comp * img_x * img_y);
   if (good == NULL) {
      free(data);
      return ep("outofmem", "Out of memory");
   }

   for (j=0; j < img_y; ++j) {
      unsigned char *src  = data + j * img_x * img_n   ;
      unsigned char *dest = good + j * img_x * req_comp;

      #define COMBO(a,b)  ((a)*8+(b))
      #define CASE(a,b)   case COMBO(a,b): for(i=0; i < img_x; ++i, src += a, dest += b)

      // convert source image with img_n components to one with req_comp components;
      // avoid switch per pixel, so use switch per scanline and massive macros
      switch(COMBO(img_n, req_comp)) {
         CASE(1,2) dest[0]=src[0], dest[1]=255; break;
         CASE(1,3) dest[0]=dest[1]=dest[2]=src[0]; break;
         CASE(1,4) dest[0]=dest[1]=dest[2]=src[0], dest[3]=255; break;
         CASE(2,1) dest[0]=src[0]; break;
         CASE(2,3) dest[0]=dest[1]=dest[2]=src[0]; break;
         CASE(2,4) dest[0]=dest[1]=dest[2]=src[0], dest[3]=src[1]; break;
         CASE(3,4) dest[0]=src[0],dest[1]=src[1],dest[2]=src[2],dest[3]=255; break;
         CASE(3,1) dest[0]=compute_y(src[0],src[1],src[2]); break;
         CASE(3,2) dest[0]=compute_y(src[0],src[1],src[2]), dest[1] = 255; break;
         CASE(4,1) dest[0]=compute_y(src[0],src[1],src[2]); break;
         CASE(4,2) dest[0]=compute_y(src[0],src[1],src[2]), dest[1] = src[3]; break;
         CASE(4,3) dest[0]=src[0],dest[1]=src[1],dest[2]=src[2]; break;
         default: assert(0);
      }
      #undef CASE
   }

   free(data);
   img_out_n = req_comp;
   return good;
}

#ifndef STBI_NO_HDR
static float   *ldr_to_hdr(stbi_uc *data, int x, int y, int comp)
{
   int i,k,n;
   float *output = (float *) malloc(x * y * comp * sizeof(float));
   if (output == NULL) { free(data); return ep("outofmem", "Out of memory"); }
   // compute number of non-alpha components
   if (comp & 1) n = comp; else n = comp-1;
   for (i=0; i < x*y; ++i) {
      for (k=0; k < n; ++k) {
         output[i*comp + k] = (float) pow(data[i*comp+k]/255.0, l2h_gamma) * l2h_scale;
      }
      if (k < comp) output[i*comp + k] = data[i*comp+k]/255.0f;
   }
   free(data);
   return output;
}

#define float2int(x)   ((int) (x))
static stbi_uc *hdr_to_ldr(float   *data, int x, int y, int comp)
{
   int i,k,n;
   stbi_uc *output = (stbi_uc *) malloc(x * y * comp);
   if (output == NULL) { free(data); return ep("outofmem", "Out of memory"); }
   // compute number of non-alpha components
   if (comp & 1) n = comp; else n = comp-1;
   for (i=0; i < x*y; ++i) {
      for (k=0; k < n; ++k) {
         float z = (float) pow(data[i*comp+k]*h2l_scale_i, h2l_gamma_i) * 255 + 0.5f;
         if (z < 0) z = 0;
         if (z > 255) z = 255;
         output[i*comp + k] = float2int(z);
      }
      if (k < comp) {
         float z = data[i*comp+k] * 255 + 0.5f;
         if (z < 0) z = 0;
         if (z > 255) z = 255;
         output[i*comp + k] = float2int(z);
      }
   }
   free(data);
   return output;
}
#endif

//////////////////////////////////////////////////////////////////////////////
//
//  "baseline" JPEG/JFIF decoder (not actually fully baseline implementation)
//
//    simple implementation
//      - channel subsampling of at most 2 in each dimension
//      - doesn't support delayed output of y-dimension
//      - simple interface (only one output format: 8-bit interleaved RGB)
//      - doesn't try to recover corrupt jpegs
//      - doesn't allow partial loading, loading multiple at once
//      - still fast on x86 (copying globals into locals doesn't help x86)
//      - allocates lots of intermediate memory (full size of all components)
//        - non-interleaved case requires this anyway
//        - allows good upsampling (see next)
//    high-quality
//      - upsampled channels are bilinearly interpolated, even across blocks
//      - quality integer IDCT derived from IJG's 'slow'
//    performance
//      - fast huffman; reasonable integer IDCT
//      - uses a lot of intermediate memory, could cache poorly
//      - load http://nothings.org/remote/anemones.jpg 3 times on 2.8Ghz P4
//          stb_jpeg:   1.34 seconds (MSVC6, default release build)
//          stb_jpeg:   1.06 seconds (MSVC6, processor = Pentium Pro)
//          IJL11.dll:  1.08 seconds (compiled by intel)
//          IJG 1998:   0.98 seconds (MSVC6, makefile provided by IJG)
//          IJG 1998:   0.95 seconds (MSVC6, makefile + proc=PPro)

int stbi_jpeg_dc_only;

// huffman decoding acceleration
#define FAST_BITS   9  // larger handles more cases; smaller stomps less cache

typedef struct
{
   uint8  fast[1 << FAST_BITS];
   // weirdly, repacking this into AoS is a 10% speed loss, instead of a win
   uint16 code[256];
   uint8  values[256];
   uint8  size[257];
   unsigned int maxcode[18];
   int    delta[17];   // old 'firstsymbol' - old 'firstcode'
} huffman;

static huffman huff_dc[4];  // baseline is 2 tables, extended is 4
static huffman huff_ac[4];
static uint8 dequant[4][64];

static int build_huffman(huffman *h, int *count)
{
   int i,j,k=0,code;
   // build size list for each symbol (from JPEG spec)
   for (i=0; i < 16; ++i)
      for (j=0; j < count[i]; ++j)
         h->size[k++] = (uint8) (i+1);
   h->size[k] = 0;

   // compute actual symbols (from jpeg spec)
   code = 0;
   k = 0;
   for(j=1; j <= 16; ++j) {
      // compute delta to add to code to compute symbol id
      h->delta[j] = k - code;
      if (h->size[k] == j) {
         while (h->size[k] == j)
            h->code[k++] = (uint16) (code++);
         if (code-1 >= (1 << j)) return e("bad code lengths","Corrupt JPEG");
      }
      // compute largest code + 1 for this size, preshifted as needed later
      h->maxcode[j] = code << (16-j);
      code <<= 1;
   }
   h->maxcode[j] = 0xffffffff;

   // build non-spec acceleration table; 255 is flag for not-accelerated
   memset(h->fast, 255, 1 << FAST_BITS);
   for (i=0; i < k; ++i) {
      int s = h->size[i];
      if (s <= FAST_BITS) {
         int c = h->code[i] << (FAST_BITS-s);
         int m = 1 << (FAST_BITS-s);
         for (j=0; j < m; ++j) {
            h->fast[c+j] = (uint8) i;
         }
      }
   }
   return 1;
}

// sizes for components, interleaved MCUs
static int img_h_max, img_v_max;
static int img_mcu_x, img_mcu_y;
static int img_mcu_w, img_mcu_h;

// definition of jpeg image component
static struct
{
   int id;
   int h,v;
   int tq;
   int hd,ha;
   int dc_pred;

   int x,y,w2,h2;
   uint8 *data;
} img_comp[4];

static unsigned long  code_buffer; // jpeg entropy-coded buffer
static int            code_bits;   // number of valid bits
static unsigned char  marker;      // marker seen while filling entropy buffer
static int            nomore;      // flag if we saw a marker so must stop

static void grow_buffer_unsafe(void)
{
   do {
      int b = nomore ? 0 : get8();
      if (b == 0xff) {
         int c = get8();
         if (c != 0) {
            marker = (unsigned char) c;
            nomore = 1;
            return;
         }
      }
      code_buffer = (code_buffer << 8) | b;
      code_bits += 8;
   } while (code_bits <= 24);
}

// (1 << n) - 1
static unsigned long bmask[17]={0,1,3,7,15,31,63,127,255,511,1023,2047,4095,8191,16383,32767,65535};

// decode a jpeg huffman value from the bitstream
__forceinline static int decode(huffman *h)
{
   unsigned int temp;
   int c,k;

   if (code_bits < 16) grow_buffer_unsafe();

   // look at the top FAST_BITS and determine what symbol ID it is,
   // if the code is <= FAST_BITS
   c = (code_buffer >> (code_bits - FAST_BITS)) & ((1 << FAST_BITS)-1);
   k = h->fast[c];
   if (k < 255) {
      if (h->size[k] > code_bits)
         return -1;
      code_bits -= h->size[k];
      return h->values[k];
   }

   // naive test is to shift the code_buffer down so k bits are
   // valid, then test against maxcode. To speed this up, we've
   // preshifted maxcode left so that it has (16-k) 0s at the
   // end; in other words, regardless of the number of bits, it
   // wants to be compared against something shifted to have 16;
   // that way we don't need to shift inside the loop.
   if (code_bits < 16)
      temp = (code_buffer << (16 - code_bits)) & 0xffff;
   else
      temp = (code_buffer >> (code_bits - 16)) & 0xffff;
   for (k=FAST_BITS+1 ; ; ++k)
      if (temp < h->maxcode[k])
         break;
   if (k == 17) {
      // error! code not found
      code_bits -= 16;
      return -1;
   }

   if (k > code_bits)
      return -1;

   // convert the huffman code to the symbol id
   c = ((code_buffer >> (code_bits - k)) & bmask[k]) + h->delta[k];
   assert((((code_buffer) >> (code_bits - h->size[c])) & bmask[h->size[c]]) == h->code[c]);

   // convert the id to a symbol
   code_bits -= k;
   return h->values[c];
}

// combined JPEG 'receive' and JPEG 'extend', since baseline
// always extends everything it receives.
__forceinline static int extend_receive(int n)
{
   unsigned int m = 1 << (n-1);
   unsigned int k;
   if (code_bits < n) grow_buffer_unsafe();
   k = (code_buffer >> (code_bits - n)) & bmask[n];
   code_bits -= n;
   // the following test is probably a random branch that won't
   // predict well. I tried to table accelerate it but failed.
   // maybe it's compiling as a conditional move?
   if (k < m)
      return (-1 << n) + k + 1;
   else
      return k;
}

// given a value that's at position X in the zigzag stream,
// where does it appear in the 8x8 matrix coded as row-major?
static uint8 dezigzag[64+15] =
{
    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,
   // let corrupt input sample past end
   63, 63, 63, 63, 63, 63, 63, 63,
   63, 63, 63, 63, 63, 63, 63
};

// decode one 64-entry block--
static int decode_block(short data[64], huffman *hdc, huffman *hac, int b)
{
   int diff,dc,k;
   int t = decode(hdc);
   if (t < 0) return e("bad huffman code","Corrupt JPEG");

   // 0 all the ac values now so we can do it 32-bits at a time
   memset(data,0,64*sizeof(data[0]));

   diff = t ? extend_receive(t) : 0;
   dc = img_comp[b].dc_pred + diff;
   img_comp[b].dc_pred = dc;
   data[0] = (short) dc;

   // decode AC components, see JPEG spec
   k = 1;
   do {
      int r,s;
      int rs = decode(hac);
      if (rs < 0) return e("bad huffman code","Corrupt JPEG");
      s = rs & 15;
      r = rs >> 4;
      if (s == 0) {
         if (rs != 0xf0) break; // end block
         k += 16;
      } else {
         k += r;
         // decode into unzigzag'd location
         data[dezigzag[k++]] = (short) extend_receive(s);
      }
   } while (k < 64);
   return 1;
}

// take a -128..127 value and clamp it and convert to 0..255
__forceinline static uint8 clamp(int x)
{
   x += 128;
   // trick to use a single test to catch both cases
   if ((unsigned int) x > 255) {
      if (x < 0) return 0;
      if (x > 255) return 255;
   }
   return (uint8) x;
}

#define f2f(x)  (int) (((x) * 4096 + 0.5))
#define fsh(x)  ((x) << 12)

// derived from jidctint -- DCT_ISLOW
#define IDCT_1D(s0,s1,s2,s3,s4,s5,s6,s7)       \
   int t0,t1,t2,t3,p1,p2,p3,p4,p5,x0,x1,x2,x3; \
   p2 = s2;                                    \
   p3 = s6;                                    \
   p1 = (p2+p3) * f2f(0.5411961f);             \
   t2 = p1 + p3*f2f(-1.847759065f);            \
   t3 = p1 + p2*f2f( 0.765366865f);            \
   p2 = s0;                                    \
   p3 = s4;                                    \
   t0 = fsh(p2+p3);                            \
   t1 = fsh(p2-p3);                            \
   x0 = t0+t3;                                 \
   x3 = t0-t3;                                 \
   x1 = t1+t2;                                 \
   x2 = t1-t2;                                 \
   t0 = s7;                                    \
   t1 = s5;                                    \
   t2 = s3;                                    \
   t3 = s1;                                    \
   p3 = t0+t2;                                 \
   p4 = t1+t3;                                 \
   p1 = t0+t3;                                 \
   p2 = t1+t2;                                 \
   p5 = (p3+p4)*f2f( 1.175875602f);            \
   t0 = t0*f2f( 0.298631336f);                 \
   t1 = t1*f2f( 2.053119869f);                 \
   t2 = t2*f2f( 3.072711026f);                 \
   t3 = t3*f2f( 1.501321110f);                 \
   p1 = p5 + p1*f2f(-0.899976223f);            \
   p2 = p5 + p2*f2f(-2.562915447f);            \
   p3 = p3*f2f(-1.961570560f);                 \
   p4 = p4*f2f(-0.390180644f);                 \
   t3 += p1+p4;                                \
   t2 += p2+p3;                                \
   t1 += p2+p4;                                \
   t0 += p1+p3;

// .344 seconds on 3*anemones.jpg
static void idct_block(uint8 *out, int out_stride, short data[64], uint8 *dequantize)
{
   int i,val[64],*v=val;
   uint8 *o,*dq = dequantize;
   short *d = data;

   if (stbi_jpeg_dc_only) {
      // ok, I don't really know why this is right, but it seems to be:
      int z = 128 + ((d[0] * dq[0]) >> 3);
      for (i=0; i < 8; ++i) {
         out[0] = out[1] = out[2] = out[3] = out[4] = out[5] = out[6] = out[7] = z;
         out += out_stride;
      }
      return;
   }

   // columns
   for (i=0; i < 8; ++i,++d,++dq, ++v) {
      // if all zeroes, shortcut -- this avoids dequantizing 0s and IDCTing
      if (d[ 8]==0 && d[16]==0 && d[24]==0 && d[32]==0
           && d[40]==0 && d[48]==0 && d[56]==0) {
         //    no shortcut                 0     seconds
         //    (1|2|3|4|5|6|7)==0          0     seconds
         //    all separate               -0.047 seconds
         //    1 && 2|3 && 4|5 && 6|7:    -0.047 seconds
         int dcterm = d[0] * dq[0] << 2;
         v[0] = v[8] = v[16] = v[24] = v[32] = v[40] = v[48] = v[56] = dcterm;
      } else {
         IDCT_1D(d[ 0]*dq[ 0],d[ 8]*dq[ 8],d[16]*dq[16],d[24]*dq[24],
                 d[32]*dq[32],d[40]*dq[40],d[48]*dq[48],d[56]*dq[56])
         // constants scaled things up by 1<<12; let's bring them back
         // down, but keep 2 extra bits of precision
         x0 += 512; x1 += 512; x2 += 512; x3 += 512;
         v[ 0] = (x0+t3) >> 10;
         v[56] = (x0-t3) >> 10;
         v[ 8] = (x1+t2) >> 10;
         v[48] = (x1-t2) >> 10;
         v[16] = (x2+t1) >> 10;
         v[40] = (x2-t1) >> 10;
         v[24] = (x3+t0) >> 10;
         v[32] = (x3-t0) >> 10;
      }
   }

   for (i=0, v=val, o=out; i < 8; ++i,v+=8,o+=out_stride) {
      // no fast case since the first 1D IDCT spread components out
      IDCT_1D(v[0],v[1],v[2],v[3],v[4],v[5],v[6],v[7])
      // constants scaled things up by 1<<12, plus we had 1<<2 from first
      // loop, plus horizontal and vertical each scale by sqrt(8) so together
      // we've got an extra 1<<3, so 1<<17 total we need to remove.
      x0 += 65536; x1 += 65536; x2 += 65536; x3 += 65536;
      o[0] = clamp((x0+t3) >> 17);
      o[7] = clamp((x0-t3) >> 17);
      o[1] = clamp((x1+t2) >> 17);
      o[6] = clamp((x1-t2) >> 17);
      o[2] = clamp((x2+t1) >> 17);
      o[5] = clamp((x2-t1) >> 17);
      o[3] = clamp((x3+t0) >> 17);
      o[4] = clamp((x3-t0) >> 17);
   }
}

#define MARKER_none  0xff
// if there's a pending marker from the entropy stream, return that
// otherwise, fetch from the stream and get a marker. if there's no
// marker, return 0xff, which is never a valid marker value
static uint8 get_marker(void)
{
   uint8 x;
   if (marker != MARKER_none) { x = marker; marker = MARKER_none; return x; }
   x = get8u();
   if (x != 0xff) return MARKER_none;
   while (x == 0xff)
      x = get8u();
   return x;
}

// in each scan, we'll have scan_n components, and the order
// of the components is specified by order[]
static int scan_n, order[4];
static int restart_interval, todo;
#define RESTART(x)     ((x) >= 0xd0 && (x) <= 0xd7)

// after a restart interval, reset the entropy decoder and
// the dc prediction
static void reset(void)
{
   code_bits = 0;
   code_buffer = 0;
   nomore = 0;
   img_comp[0].dc_pred = img_comp[1].dc_pred = img_comp[2].dc_pred = 0;
   marker = MARKER_none;
   todo = restart_interval ? restart_interval : 0x7fffffff;
   // no more than 1<<31 MCUs if no restart_interal? that's plenty safe,
   // since we don't even allow 1<<30 pixels
}

static int parse_entropy_coded_data(void)
{
   reset();
   if (scan_n == 1) {
      int i,j;
      short data[64];
      int n = order[0];
      // non-interleaved data, we just need to process one block at a time,
      // in trivial scanline order
      // number of blocks to do just depends on how many actual "pixels" this
      // component has, independent of interleaved MCU blocking and such
      int w = (img_comp[n].x+7) >> 3;
      int h = (img_comp[n].y+7) >> 3;
      for (j=0; j < h; ++j) {
         for (i=0; i < w; ++i) {
            if (!decode_block(data, huff_dc+img_comp[n].hd, huff_ac+img_comp[n].ha, n)) return 0;
            idct_block(img_comp[n].data+img_comp[n].w2*j*8+i*8, img_comp[n].w2, data, dequant[img_comp[n].tq]);
            // every data block is an MCU, so countdown the restart interval
            if (--todo <= 0) {
               if (code_bits < 24) grow_buffer_unsafe();
               // if it's NOT a restart, then just bail, so we get corrupt data
               // rather than no data
               if (!RESTART(marker)) return 1;
               reset();
            }
         }
      }
   } else { // interleaved!
      int i,j,k,x,y;
      short data[64];
      for (j=0; j < img_mcu_y; ++j) {
         for (i=0; i < img_mcu_x; ++i) {
            // scan an interleaved mcu... process scan_n components in order
            for (k=0; k < scan_n; ++k) {
               int n = order[k];
               // scan out an mcu's worth of this component; that's just determined
               // by the basic H and V specified for the component
               for (y=0; y < img_comp[n].v; ++y) {
                  for (x=0; x < img_comp[n].h; ++x) {
                     int x2 = (i*img_comp[n].h + x)*8;
                     int y2 = (j*img_comp[n].v + y)*8;
                     if (!decode_block(data, huff_dc+img_comp[n].hd, huff_ac+img_comp[n].ha, n)) return 0;
                     idct_block(img_comp[n].data+img_comp[n].w2*y2+x2, img_comp[n].w2, data, dequant[img_comp[n].tq]);
                  }
               }
            }
            // after all interleaved components, that's an interleaved MCU,
            // so now count down the restart interval
            if (--todo <= 0) {
               if (code_bits < 24) grow_buffer_unsafe();
               // if it's NOT a restart, then just bail, so we get corrupt data
               // rather than no data
               if (!RESTART(marker)) return 1;
               reset();
            }
         }
      }
   }
   return 1;
}

static int process_marker(int m)
{
   int L;
   switch (m) {
      case MARKER_none: // no marker found
         return e("expected marker","Corrupt JPEG");

      case 0xC2: // SOF - progressive
         return e("progressive jpeg","JPEG format not supported (progressive)");

      case 0xDD: // DRI - specify restart interval
         if (get16() != 4) return e("bad DRI len","Corrupt JPEG");
         restart_interval = get16();
         return 1;

      case 0xDB: // DQT - define quantization table
         L = get16()-2;
         while (L > 0) {
            int z = get8();
            int p = z >> 4;
            int t = z & 15,i;
            if (p != 0) return e("bad DQT type","Corrupt JPEG");
            if (t > 3) return e("bad DQT table","Corrupt JPEG");
            for (i=0; i < 64; ++i)
               dequant[t][dezigzag[i]] = get8u();
            L -= 65;
         }
         return L==0;

      case 0xC4: // DHT - define huffman table
         L = get16()-2;
         while (L > 0) {
            uint8 *v;
            int sizes[16],i,m=0;
            int z = get8();
            int tc = z >> 4;
            int th = z & 15;
            if (tc > 1 || th > 3) return e("bad DHT header","Corrupt JPEG");
            for (i=0; i < 16; ++i) {
               sizes[i] = get8();
               m += sizes[i];
            }
            L -= 17;
            if (tc == 0) {
               if (!build_huffman(huff_dc+th, sizes)) return 0;
               v = huff_dc[th].values;
            } else {
               if (!build_huffman(huff_ac+th, sizes)) return 0;
               v = huff_ac[th].values;
            }
            for (i=0; i < m; ++i)
               v[i] = get8u();
            L -= m;
         }
         return L==0;
   }
   // check for comment block or APP blocks
   if ((m >= 0xE0 && m <= 0xEF) || m == 0xFE) {
      skip(get16()-2);
      return 1;
   }
   return 0;
}

// after we see SOS
static int process_scan_header(void)
{
   int i;
   int Ls = get16();
   scan_n = get8();
   if (scan_n < 1 || scan_n > 4 || scan_n > (int) img_n) return e("bad SOS component count","Corrupt JPEG");
   if (Ls != 6+2*scan_n) return e("bad SOS len","Corrupt JPEG");
   for (i=0; i < scan_n; ++i) {
      int id = get8(), which;
      int z = get8();
      for (which = 0; which < img_n; ++which)
         if (img_comp[which].id == id)
            break;
      if (which == img_n) return 0;
      img_comp[which].hd = z >> 4;   if (img_comp[which].hd > 3) return e("bad DC huff","Corrupt JPEG");
      img_comp[which].ha = z & 15;   if (img_comp[which].ha > 3) return e("bad AC huff","Corrupt JPEG");
      order[i] = which;
   }
   if (get8() != 0) return e("bad SOS","Corrupt JPEG");
   get8(); // should be 63, but might be 0
   if (get8() != 0) return e("bad SOS","Corrupt JPEG");

   return 1;
}

static int process_frame_header(int scan)
{
   int Lf,p,i,z, h_max=1,v_max=1;
   Lf = get16();         if (Lf < 11) return e("bad SOF len","Corrupt JPEG"); // JPEG
   p  = get8();          if (p != 8) return e("only 8-bit","JPEG format not supported: 8-bit only"); // JPEG baseline
   img_y = get16();      if (img_y == 0) return e("no header height", "JPEG format not supported: delayed height"); // Legal, but we don't handle it--but neither does IJG
   img_x = get16();      if (img_x == 0) return e("0 width","Corrupt JPEG"); // JPEG requires
   img_n = get8();
   if (img_n != 3 && img_n != 1) return e("bad component count","Corrupt JPEG");    // JFIF requires

   if (Lf != 8+3*img_n) return e("bad SOF len","Corrupt JPEG");

   for (i=0; i < img_n; ++i) {
      img_comp[i].id = get8();
      if (img_comp[i].id != i+1)   // JFIF requires
         if (img_comp[i].id != i)  // jpegtran outputs non-JFIF-compliant files!
            return e("bad component ID","Corrupt JPEG");
      z = get8();
      img_comp[i].h = (z >> 4);  if (!img_comp[i].h || img_comp[i].h > 4) return e("bad H","Corrupt JPEG");
      img_comp[i].v = z & 15;    if (!img_comp[i].v || img_comp[i].v > 4) return e("bad V","Corrupt JPEG");
      img_comp[i].tq = get8();   if (img_comp[i].tq > 3) return e("bad TQ","Corrupt JPEG");
   }

   if (scan != SCAN_load) return 1;

   if ((1 << 30) / img_x / img_n < img_y) return e("too large", "Image too large to decode");

   for (i=0; i < img_n; ++i) {
      if (img_comp[i].h > h_max) h_max = img_comp[i].h;
      if (img_comp[i].v > v_max) v_max = img_comp[i].v;
   }

   // compute interleaved mcu info
   img_h_max = h_max;
   img_v_max = v_max;
   img_mcu_w = h_max * 8;
   img_mcu_h = v_max * 8;
   img_mcu_x = (img_x + img_mcu_w-1) / img_mcu_w;
   img_mcu_y = (img_y + img_mcu_h-1) / img_mcu_h;

   for (i=0; i < img_n; ++i) {
      // number of effective pixels (e.g. for non-interleaved MCU)
      img_comp[i].x = (img_x * img_comp[i].h + h_max-1) / h_max;
      img_comp[i].y = (img_y * img_comp[i].v + v_max-1) / v_max;
      // to simplify generation, we'll allocate enough memory to decode
      // the bogus oversized data from using interleaved MCUs and their
      // big blocks (e.g. a 16x16 iMCU on an image of width 33); we won't
      // discard the extra data until colorspace conversion
      img_comp[i].w2 = img_mcu_x * img_comp[i].h * 8;
      img_comp[i].h2 = img_mcu_y * img_comp[i].v * 8;
      img_comp[i].data = (uint8 *) malloc(img_comp[i].w2 * img_comp[i].h2);
      if (img_comp[i].data == NULL) {
         for(--i; i >= 0; --i)
            free(img_comp[i].data);
         return e("outofmem", "Out of memory");
      }
   }

   return 1;
}

// use comparisons since in some cases we handle more than one case (e.g. SOF)
#define DNL(x)         ((x) == 0xdc)
#define SOI(x)         ((x) == 0xd8)
#define EOI(x)         ((x) == 0xd9)
#define SOF(x)         ((x) == 0xc0 || (x) == 0xc1)
#define SOS(x)         ((x) == 0xda)

static int decode_jpeg_header(int scan)
{
   int m;
   marker = MARKER_none; // initialize cached marker to empty
   m = get_marker();
   if (!SOI(m)) return e("no SOI","Corrupt JPEG");
   if (scan == SCAN_type) return 1;
   m = get_marker();
   while (!SOF(m)) {
      if (!process_marker(m)) return 0;
      m = get_marker();
      while (m == MARKER_none) {
         // some files have extra padding after their blocks, so ok, we'll scan
         if (at_eof()) return e("no SOF", "Corrupt JPEG");
         m = get_marker();
      }
   }
   if (!process_frame_header(scan)) return 0;
   return 1;
}

static int decode_jpeg_image(void)
{
   int m;
   restart_interval = 0;
   if (!decode_jpeg_header(SCAN_load)) return 0;
   m = get_marker();
   while (!EOI(m)) {
      if (SOS(m)) {
         if (!process_scan_header()) return 0;
         if (!parse_entropy_coded_data()) return 0;
      } else {
         if (!process_marker(m)) return 0;
      }
      m = get_marker();
   }
   return 1;
}

// static jfif-centered resampling with cross-block smoothing
// here by cross-block smoothing what I mean is that the resampling
// is bilerp and crosses blocks; I dunno what IJG means

#define div4(x) ((uint8) ((x) >> 2))

static void resample_v_2(uint8 *out1, uint8 *input, int w, int h, int s)
{
   // need to generate two samples vertically for every one in input
   uint8 *above;
   uint8 *below;
   uint8 *source;
   uint8 *out2;
   int i,j;
   source = input;
   out2 = out1+w;
   for (j=0; j < h; ++j) {
      above = source;
      source = input + j*s;
      below = source + s; if (j == h-1) below = source;
      for (i=0; i < w; ++i) {
         int n = source[i]*3;
         out1[i] = div4(above[i] + n);
         out2[i] = div4(below[i] + n);
      }
      out1 += w*2;
      out2 += w*2;
   }
}

static void resample_h_2(uint8 *out, uint8 *input, int w, int h, int s)
{
   // need to generate two samples horizontally for every one in input
   int i,j;
   if (w == 1) {
      for (j=0; j < h; ++j)
         out[j*2+0] = out[j*2+1] = input[j*s];
      return;
   }
   for (j=0; j < h; ++j) {
      out[0] = input[0];
      out[1] = div4(input[0]*3 + input[1]);
      for (i=1; i < w-1; ++i) {
         int n = input[i]*3;
         out[i*2-2] = div4(input[i-1] + n);
         out[i*2-1] = div4(input[i+1] + n);
      }
      out[w*2-2] = div4(input[w-2]*3 + input[w-1]);
      out[w*2-1] = input[w-1];
      out += w*2;
      input += s;
   }
}

// .172 seconds on 3*anemones.jpg
static void resample_hv_2(uint8 *out, uint8 *input, int w, int h, int s)
{
   // need to generate 2x2 samples for every one in input
   int i,j;
   int os = w*2;
   // generate edge samples... @TODO lerp them!
   for (i=0; i < w; ++i) {
      out[i*2+0] = out[i*2+1] = input[i];
      out[i*2+(2*h-1)*os+0] = out[i*2+(2*h-1)*os+1] = input[i+(h-1)*w];
   }
   for (j=0; j < h; ++j) {
      out[j*os*2+0] = out[j*os*2+os+0] = input[j*w];
      out[j*os*2+os-1] = out[j*os*2+os+os-1] = input[j*w+i-1];
   }
   // now generate interior samples; i & j point to top left of input
   for (j=0; j < h-1; ++j) {
      uint8 *in1 = input+j*s;
      uint8 *in2 = in1 + s;
      uint8 *out1 = out + (j*2+1)*os + 1;
      uint8 *out2 = out1 + os;
      for (i=0; i < w-1; ++i) {
         int p00 = in1[0], p01=in1[1], p10=in2[0], p11=in2[1];
         int p00_3 = p00*3, p01_3 = p01*3, p10_3 = p10*3, p11_3 = p11*3;

         #define div16(x)  ((uint8) ((x) >> 4))

         out1[0] = div16(p00*9 + p01_3 + p10_3 + p11);
         out1[1] = div16(p01*9 + p00_3 + p01_3 + p10);
         out2[0] = div16(p10*9 + p11_3 + p00_3 + p01);
         out2[1] = div16(p11*9 + p10_3 + p01_3 + p00);
         out1 += 2;
         out2 += 2;
         ++in1;
         ++in2;
      }
   }
}

#define float2fixed(x)  ((int) ((x) * 65536 + 0.5))

// 0.38 seconds on 3*anemones.jpg   (0.25 with processor = Pro)
// VC6 without processor=Pro is generating multiple LEAs per multiply!
static void YCbCr_to_RGB_row(uint8 *out, uint8 *y, uint8 *pcb, uint8 *pcr, int count, int step)
{
   int i;
   for (i=0; i < count; ++i) {
      int y_fixed = (y[i] << 16) + 32768; // rounding
      int r,g,b;
      int cr = pcr[i] - 128;
      int cb = pcb[i] - 128;
      r = y_fixed + cr*float2fixed(1.40200f);
      g = y_fixed - cr*float2fixed(0.71414f) - cb*float2fixed(0.34414f);
      b = y_fixed                            + cb*float2fixed(1.77200f);
      r >>= 16;
      g >>= 16;
      b >>= 16;
      if ((unsigned) r > 255) { if (r < 0) r = 0; else r = 255; }
      if ((unsigned) g > 255) { if (g < 0) g = 0; else g = 255; }
      if ((unsigned) b > 255) { if (b < 0) b = 0; else b = 255; }
      out[0] = (uint8)r;
      out[1] = (uint8)g;
      out[2] = (uint8)b;
      if (step == 4) out[3] = 255;
      out += step;
   }
}

// clean up the temporary component buffers
static void cleanup_jpeg(void)
{
   int i;
   for (i=0; i < img_n; ++i) {
      if (img_comp[i].data) {
         free(img_comp[i].data);
         img_comp[i].data = NULL;
      }
   }
}

static uint8 *load_jpeg_image(int *out_x, int *out_y, int *comp, int req_comp)
{
   int i, n;
   // validate req_comp
   if (req_comp < 0 || req_comp > 4) return ep("bad req_comp", "Internal error");

   // load a jpeg image from whichever source
   if (!decode_jpeg_image()) { cleanup_jpeg(); return NULL; }

   // determine actual number of components to generate
   n = req_comp ? req_comp : img_n;

   // resample components to full size... memory wasteful, but this
   // lets us bilerp across blocks while upsampling
   for (i=0; i < img_n; ++i) {
      // if we're outputting fewer than 3 components, we're grey not RGB;
      // in that case, don't bother upsampling Cb or Cr
      if (n < 3 && i) continue;

      // check if the component scale is less than max; if so it needs upsampling
      if (img_comp[i].h != img_h_max || img_comp[i].v != img_v_max) {
         int stride = img_x;
         // allocate final size; make sure it's big enough for upsampling off
         // the edges with upsample up to 4x4 (although we only support 2x2
         // currently)
         uint8 *new_data = (uint8 *) malloc((img_x+3)*(img_y+3));
         if (new_data == NULL) {
            cleanup_jpeg();
            return ep("outofmem", "Out of memory (image too large?)");
         }
         if (img_comp[i].h*2 == img_h_max && img_comp[i].v*2 == img_v_max) {
            int tx = (img_x+1)>>1;
            resample_hv_2(new_data, img_comp[i].data, tx,(img_y+1)>>1, img_comp[i].w2);
            stride = tx*2;
         } else if (img_comp[i].h == img_h_max && img_comp[i].v*2 == img_v_max) {
            resample_v_2(new_data, img_comp[i].data, img_x,(img_y+1)>>1, img_comp[i].w2);
         } else if (img_comp[i].h*2 == img_h_max && img_comp[i].v == img_v_max) {
            int tx = (img_x+1)>>1;
            resample_h_2(new_data, img_comp[i].data, tx,img_y, img_comp[i].w2);
            stride = tx*2;
         } else {
            // @TODO resample uncommon sampling pattern with nearest neighbor
            free(new_data);
            cleanup_jpeg();
            return ep("uncommon H or V", "JPEG not supported: atypical downsampling mode");
         }
         img_comp[i].w2 = stride;
         free(img_comp[i].data);
         img_comp[i].data = new_data;
      }
   }

   // now convert components to output image
   {
      uint32 i,j;
      uint8 *output = (uint8 *) malloc(n * img_x * img_y + 1);
      if (n >= 3) { // output STBI_rgb_*
         for (j=0; j < img_y; ++j) {
            uint8 *y  = img_comp[0].data + j*img_comp[0].w2;
            uint8 *out = output + n * img_x * j;
            if (img_n == 3) {
               uint8 *cb = img_comp[1].data + j*img_comp[1].w2;
               uint8 *cr = img_comp[2].data + j*img_comp[2].w2;
               YCbCr_to_RGB_row(out, y, cb, cr, img_x, n);
            } else {
               for (i=0; i < img_x; ++i) {
                  out[0] = out[1] = out[2] = y[i];
                  out[3] = 255; // not used if n == 3
                  out += n;
               }
            }
         }
      } else {      // output STBI_grey_*
         for (j=0; j < img_y; ++j) {
            uint8 *y  = img_comp[0].data + j*img_comp[0].w2;
            uint8 *out = output + n * img_x * j;
            if (n == 1)
               for (i=0; i < img_x; ++i) *out++ = *y++;
            else
               for (i=0; i < img_x; ++i) *out++ = *y++, *out++ = 255;
         }
      }
      cleanup_jpeg();
      *out_x = img_x;
      *out_y = img_y;
      if (comp) *comp  = n; // Changed JLD: report output components
      //if (comp) *comp  = img_n; // report original components, not output
      return output;
   }
}

#ifndef STBI_NO_STDIO
unsigned char *stbi_jpeg_load_from_file(FILE *f, int *x, int *y, int *comp, int req_comp)
{
   start_file(f);
   return load_jpeg_image(x,y,comp,req_comp);
}

unsigned char *stbi_jpeg_load(char *filename, int *x, int *y, int *comp, int req_comp)
{
   unsigned char *data;
   FILE *f = fopen(filename, "rb");
   if (!f) return NULL;
   data = stbi_jpeg_load_from_file(f,x,y,comp,req_comp);
   fclose(f);
   return data;
}
#endif

unsigned char *stbi_jpeg_load_from_memory(stbi_uc *buffer, int len, int *x, int *y, int *comp, int req_comp)
{
   start_mem(buffer,len);
   return load_jpeg_image(x,y,comp,req_comp);
}

#ifndef STBI_NO_STDIO
int stbi_jpeg_test_file(FILE *f)
{
   int n,r;
   n = ftell(f);
   start_file(f);
   r = decode_jpeg_header(SCAN_type);
   fseek(f,n,SEEK_SET);
   return r;
}
#endif

int stbi_jpeg_test_memory(unsigned char *buffer, int len)
{
   start_mem(buffer,len);
   return decode_jpeg_header(SCAN_type);
}

// @TODO:
#ifndef STBI_NO_STDIO
extern int      stbi_jpeg_info            (char *filename,           int *x, int *y, int *comp);
extern int      stbi_jpeg_info_from_file  (FILE *f,                  int *x, int *y, int *comp);
#endif
extern int      stbi_jpeg_info_from_memory(stbi_uc *buffer, int len, int *x, int *y, int *comp);

// public domain zlib decode    v0.2  Sean Barrett 2006-11-18
//    simple implementation
//      - all input must be provided in an upfront buffer
//      - all output is written to a single output buffer (can malloc/realloc)
//    performance
//      - fast huffman

// fast-way is faster to check than jpeg huffman, but slow way is slower
#define ZFAST_BITS  9 // accelerate all cases in default tables
#define ZFAST_MASK  ((1 << ZFAST_BITS) - 1)

// zlib-style huffman encoding
// (jpegs packs from left, zlib from right, so can't share code)
typedef struct
{
   uint16 fast[1 << ZFAST_BITS];
   uint16 firstcode[16];
   int maxcode[17];
   uint16 firstsymbol[16];
   uint8  size[288];
   uint16 value[288];
} zhuffman;

__forceinline static int bitreverse16(int n)
{
  n = ((n & 0xAAAA) >>  1) | ((n & 0x5555) << 1);
  n = ((n & 0xCCCC) >>  2) | ((n & 0x3333) << 2);
  n = ((n & 0xF0F0) >>  4) | ((n & 0x0F0F) << 4);
  n = ((n & 0xFF00) >>  8) | ((n & 0x00FF) << 8);
  return n;
}

__forceinline static int bit_reverse(int v, int bits)
{
   assert(bits <= 16);
   // to bit reverse n bits, reverse 16 and shift
   // e.g. 11 bits, bit reverse and shift away 5
   return bitreverse16(v) >> (16-bits);
}

static int zbuild_huffman(zhuffman *z, uint8 *sizelist, int num)
{
   int i,k=0;
   int code, next_code[16], sizes[17];

   // DEFLATE spec for generating codes
   memset(sizes, 0, sizeof(sizes));
   memset(z->fast, 255, sizeof(z->fast));
   for (i=0; i < num; ++i)
      ++sizes[sizelist[i]];
   sizes[0] = 0;
   for (i=1; i < 16; ++i)
      assert(sizes[i] <= (1 << i));
   code = 0;
   for (i=1; i < 16; ++i) {
      next_code[i] = code;
      z->firstcode[i] = (uint16) code;
      z->firstsymbol[i] = (uint16) k;
      code = (code + sizes[i]);
      if (sizes[i])
         if (code-1 >= (1 << i)) return e("bad codelengths","Corrupt JPEG");
      z->maxcode[i] = code << (16-i); // preshift for inner loop
      code <<= 1;
      k += sizes[i];
   }
   z->maxcode[16] = 0x10000; // sentinel
   for (i=0; i < num; ++i) {
      int s = sizelist[i];
      if (s) {
         int c = next_code[s] - z->firstcode[s] + z->firstsymbol[s];
         z->size[c] = (uint8)s;
         z->value[c] = (uint16)i;
         if (s <= ZFAST_BITS) {
            int k = bit_reverse(next_code[s],s);
            while (k < (1 << ZFAST_BITS)) {
               z->fast[k] = (uint16) c;
               k += (1 << s);
            }
         }
         ++next_code[s];
      }
   }
   return 1;
}

// zlib-from-memory implementation for PNG reading
//    because PNG allows splitting the zlib stream arbitrarily,
//    and it's annoying structurally to have PNG call ZLIB call PNG,
//    we require PNG read all the IDATs and combine them into a single
//    memory buffer

static uint8 *zbuffer, *zbuffer_end;…