gif.c

#ifndef JO_GIF_H
#define JO_GIF_H

#include <stdio.h>

typedef struct jo_gif_kd_node_s {
  unsigned char pnt[3];
  unsigned char lhs, rhs;
  unsigned char idx;
} jo_gif_kd_node_t;

typedef struct {
  /* in */
  FILE* fp;
  void (*gamma_f)(unsigned char*, int);

  /* internal */
  unsigned char *buf[2];
  int frame;
  short w, h;

  short col_cnt;
  short repeat;
  unsigned char pal_siz;
  unsigned char dither;

  unsigned char pal[0x300]; 
  jo_gif_kd_node_t kd[0x100];
  short kdroot;
} jo_gif_t;

/* width/height | the same for every frame */
/* repeat       | 0 = loop forever, 1 = loop once, etc... */
/* palSize      | must be power of 2 - 1. so, 255 not 256. */
extern int
jo_gif_start(jo_gif_t*, const char* path, short w, short h, short repeat,
             int pal_siz);

/* gif        | the state (returned from jo_gif_start) */
/* img        | the pixels */
/* delay_s    | amount of time in between frames (in centiseconds) */
/* gen_pal    | true if you want a unique palette generated for this frame */
/* (does not effect future frames) */
extern void
jo_gif_frame(jo_gif_t*, unsigned char* img, short delay_cs, int gen_pal);

/* gif | the state (returned from jo_gif_start) */
extern void
jo_gif_end(jo_gif_t* gif);

#endif

#include <assert.h>
#include <memory.h>
#include <stdlib.h>

#define JO_GIF_HASH_SIZE 5003

#define jo__gif_swapi(a, b) ((a) ^= (b), (b) ^= (a), (a) ^= (b))
#define jo__gif_clamp(a, b, c) (((a) < (b)) ? (b) : (((a) > (c)) ? (c) : (a)))

static int
jo__gif_ilog2(int n) {
#ifdef _MSC_VER
  unsigned long msbp;
  _BitScanReverse(&msbp, (unsigned long)n);
  return (int)msbp;
#elif defined(__GNUC__) || defined(__clang__)
  return (int)sizeof(unsigned long) * 8 - 1 - __builtin_clzl((unsigned long)n);
#else
  #define lt(n) n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n
  static const char tbl[256] = {
    0,0,1,1,2,2,2,2,3,3,3,3,3,3,3,3,
    lt(4), lt(5), lt(5), lt(6), lt(6), lt(6), lt(6),
    lt(7), lt(7), lt(7), lt(7), lt(7), lt(7), lt(7), lt(7)};
  return tbl[n];
  #undef lt
#endif
}
static void
jo__kd_swap(jo_gif_kd_node_t* n, int x, int y) {
  if (x == y) return;
  jo__gif_swapi(n[x].pnt[0], n[y].pnt[0]);
  jo__gif_swapi(n[x].pnt[1], n[y].pnt[1]);
  jo__gif_swapi(n[x].pnt[2], n[y].pnt[2]);
  jo__gif_swapi(n[x].idx, n[y].idx);
}
static int
jo__gif_kd_median(jo_gif_kd_node_t* n, int start, int end, int idx) {
  int pivot, md, p, store;
  if (end <= start) {
    return -1;
  } else if (end == start + 1) {
    return start;
  }
  md = start + (end - start) / 2;
  while (1) {
    pivot = n[md].pnt[idx];
    jo__kd_swap(n, md, (end - 1));
    for (store = p = start; p < end; p++) {
      if (n[p].pnt[idx] < pivot) {
        if (p != store) {
          jo__kd_swap(n, p, store);
        }
        store++;
      }
    }
    jo__kd_swap(n, store, end - 1);
    if (n[store].pnt[idx] == n[md].pnt[idx]) {
      return md;
    }
    if (store > md) {
      end = store;
    } else {
      start = store;
    }
  }
}
static int
jo__gif_kd_build(jo_gif_kd_node_t* n, int t, int cnt, int i) {
  int res = 0;
  if (!cnt) {
    return 0;
  }
  res = jo__gif_kd_median(n, t, t + cnt, i);
  if (res >= 0) {
    i = (i + 1) % 3;
    n->lhs = jo__gif_kd_build(n, t, res - t, i);
    n->rhs = jo__gif_kd_build(n, res + 1, t + cnt - (res + 1), i);
  }
  return res;
}

typedef struct jo_gif_kd_res_s {
  int node;
  int dist;
} jo_gif_kd_near_t;

static int
jo__gif_kd_dist(const jo_gif_kd_node_t* a, const unsigned char* pnt) {
  int i, d = 0;
  for(i = 0; i < 3; ++i) {
    int t = a->pnt[i] - pnt[i];
    d += t * t;
  }
  return d;
}
static void
jo__gif_kd_nearest(jo_gif_kd_near_t* res, const jo_gif_kd_node_t* n, int root,
                   const unsigned char* pnt, int i) {
  int d, dx;
  if (!root) {
    return;
  }
  d = jo__gif_kd_dist(n + root, pnt);
  dx = n[root].pnt[i] - pnt[i];
  if (res->node < 0 || d < res->dist) {
    res->node = root;
    res->dist = d;
  }
  if (res->dist >= 16) {
    int dx2 = dx * dx;
    i = (i + 1) % 3;
    jo__gif_kd_nearest(res, n, dx > 0 ? n[root].lhs : n[root].rhs, pnt, i);
    if (dx2 < res->dist) {
      jo__gif_kd_nearest(res, n, dx > 0 ? n[root].rhs : n[root].lhs, pnt, i);
    }
  }
}
static void
jo__gif_quantize(unsigned char* rgba, int rgbaSize, int sample,
                 unsigned char* map, int col_cnt) {
  /* Based on NeuQuant algorithm */
  /* defs for freq and bias */
  const int intbiasshift = 16; /* bias for fractions */
  const int intbias = (((int)1) << intbiasshift);
  const int gammashift = 10; /* gamma = 1024 */
  const int betashift = 10;
  const int beta = (intbias >> betashift); /* beta = 1/1024 */
  const int betagamma = (intbias << (gammashift - betashift));

  /* defs for decreasing radius factor */
  const int radiusbiasshift = 6; /* at 32.0 biased by 6 bits */
  const int radiusbias = (((int)1) << radiusbiasshift);
  const int radiusdec = 30; /* factor of 1/30 each cycle */

  /* defs for decreasing alpha factor */
  const int alphabiasshift = 10; /* alpha starts at 1.0 */
  const int initalpha = (((int)1) << alphabiasshift);

  /* radbias and alpharadbias used for rad_pow calculation */
  const int radbiasshift = 8;
  const int radbias = (((int)1) << radbiasshift);
  const int alpharadbshift = (alphabiasshift + radbiasshift);
  const int alpharadbias = (((int)1) << alpharadbshift);

  int i, k, jj, m;
  int net[256][3];
  int pix = 0, bias[256] = {0}, freq[256];
  sample = sample < 1 ? 1 : sample > 30 ? 30 : sample;
  for (i = 0; i < col_cnt; ++i) {
    /* Put nurons evenly through the luminance spectrum. */
    net[i][0] = net[i][1] = net[i][2] = (i << 12) / col_cnt;
    freq[i] = intbias / col_cnt;
  }
  {
    /* Learn */
    int step = 4;
    static const int primes[5] = {499, 491, 487, 503};
    for (i = 0; i < 4; ++i) {
      /* TODO/Error? primes[i]*4? */
      if (rgbaSize > primes[i] * 4 && (rgbaSize % primes[i])) {  
        step = primes[i] * 4;
      }
    }
    sample = step == 4 ? 1 : sample;

    {int alphadec = 30 + ((sample - 1) / 3);
    int samplepixels = rgbaSize / (4 * sample);
    int d = samplepixels / 100;
    int alpha = initalpha;
    int delta = d == 0 ? 1 : d;

    int radius = (col_cnt >> 3) * radiusbias;
    int rads = radius >> radiusbiasshift;
    int rad = rads <= 1 ? 0 : rads;
    int radSq = rad * rad;
    int rad_pow[32];
    for (i = 0; i < rad; i++) {
      rad_pow[i] = alpha * (((radSq - i * i) * radbias) / radSq);
    }
    /* Randomly walk through the pixels and relax neurons to the "optimal" target. */
    for (i = 0; i < samplepixels;) {
      int r = rgba[pix + 0] << 4;
      int g = rgba[pix + 1] << 4;
      int b = rgba[pix + 2] << 4;
      int j = -1;
      {
        /* finds closest neuron (min dist) and updates freq */
        /* finds best neuron (min dist-bias) and returns position */
        /* for frequently chosen neurons, freq[k] is high and bias[k] is */
        /* negative bias[k] = gamma*((1/col_cnt)-freq[k]) */
        int bestd = 0x7FFFFFFF, bestbiasd = 0x7FFFFFFF, bestpos = -1;
        for (k = 0; k < col_cnt; k++) {
          int* n = net[k];
          int biasdist, betafreq;
          int dist = abs(n[0] - r) + abs(n[1] - g) + abs(n[2] - b);
          if (dist < bestd) {
            bestd = dist;
            bestpos = k;
          }
          biasdist = dist - ((bias[k]) >> (intbiasshift - 4));
          if (biasdist < bestbiasd) {
            bestbiasd = biasdist;
            j = k;
          }
          betafreq = freq[k] >> betashift;
          freq[k] -= betafreq;
          bias[k] += betafreq << gammashift;
        }
        freq[bestpos] += beta;
        bias[bestpos] -= betagamma;
      }
      /* Move neuron j towards biased (b,g,r) by factor alpha */
      net[j][0] -= (net[j][0] - r) * alpha / initalpha;
      net[j][1] -= (net[j][1] - g) * alpha / initalpha;
      net[j][2] -= (net[j][2] - b) * alpha / initalpha;
      if (rad != 0) {
        /* Move adjacent neurons by precomputed alpha*(1-((i-j)^2/[r]^2)) in */
        /* rad_pow[|i-j|] */
        int lo = j - rad < -1 ? -1 : j - rad;
        int hi = j + rad;

        hi = hi > col_cnt ? col_cnt : hi;
        for (jj = j + 1, m = 1; jj < hi; ++jj) {
          int a = rad_pow[m++];
          net[jj][0] -= (net[jj][0] - r) * a / alpharadbias;
          net[jj][1] -= (net[jj][1] - g) * a / alpharadbias;
          net[jj][2] -= (net[jj][2] - b) * a / alpharadbias;
        }
        for (k = j - 1, m = 1; k > lo; --k) {
          int a = rad_pow[m++];
          net[k][0] -= (net[k][0] - r) * a / alpharadbias;
          net[k][1] -= (net[k][1] - g) * a / alpharadbias;
          net[k][2] -= (net[k][2] - b) * a / alpharadbias;
        }
      }
      pix += step;
      pix = pix >= rgbaSize ? pix - rgbaSize : pix;

      /* every 1% of the image, move less over the following iterations. */
      if (++i % delta == 0) {
        alpha -= alpha / alphadec;
        radius -= radius / radiusdec;
        rad = radius >> radiusbiasshift;
        rad = rad <= 1 ? 0 : rad;
        radSq = rad * rad;
        for (j = 0; j < rad; j++) {
          rad_pow[j] = alpha * ((radSq - j * j) * radbias / radSq);
        }
      }}
    }
  }
  /* Unbias net to give byte values 0..255 */
  for (i = 0; i < col_cnt; i++) {
    map[i * 3 + 0] = net[i][0] >>= 4;
    map[i * 3 + 1] = net[i][1] >>= 4;
    map[i * 3 + 2] = net[i][2] >>= 4;
  }
}

typedef struct {
  FILE* fp;
  unsigned char buf[256];
  int idx;

  int bit_cnt;
  int out_bits;
  int cur_bits;
} jo_gif_lzw_t;

static void
jo__gif_lzw_write(jo_gif_lzw_t* s, int code) {
  s->out_bits |= code << s->cur_bits;
  s->cur_bits += s->bit_cnt;

  while (s->cur_bits >= 8) {
    s->buf[s->idx++] = s->out_bits & 255;
    s->out_bits >>= 8;
    s->cur_bits -= 8;

    if (s->idx >= 255) {
      putc(s->idx, s->fp);
      fwrite(s->buf, s->idx, 1, s->fp);
      s->idx = 0;
    }
  }
}
static void
jo__gif_lzw_encode(unsigned char* in, int len, FILE* fp) {
  /* Note: 30k stack space for dictionary =| */
  short code_tbl[JO_GIF_HASH_SIZE];
  int hash_tbl[JO_GIF_HASH_SIZE];
  int maxcode = 511;
  int free_ent = 0x102;
  int ent = *in++;

  jo_gif_lzw_t state;
  state.bit_cnt = 9;
  state.fp = fp;

  memset(hash_tbl, 0xFF, sizeof(hash_tbl));
  jo__gif_lzw_write(&state, 0x100);

l0:
  while (--len) {
    int c = *in++;
    int fcode = (c << 12) + ent;
    int key = (c << 4) ^ ent;  /* xor hashing */
    while (hash_tbl[key] >= 0) {
      if (hash_tbl[key] == fcode) {
        ent = code_tbl[key];
        goto l0;
      }
      ++key;
      key = key >= JO_GIF_HASH_SIZE ? key - JO_GIF_HASH_SIZE : key;
    }
    jo__gif_lzw_write(&state, ent);
    ent = c;
    if (free_ent < 4096) {
      if (free_ent > maxcode) {
        ++state.bit_cnt;
        if (state.bit_cnt == 12) {
          maxcode = 4096;
        } else {
          maxcode = (1 << state.bit_cnt) - 1;
        }
      }
      code_tbl[key] = free_ent++;
      hash_tbl[key] = fcode;
    } else {
      memset(hash_tbl, 0xFF, sizeof(hash_tbl));
      free_ent = 0x102;
      jo__gif_lzw_write(&state, 0x100);
      state.bit_cnt = 9;
      maxcode = 511;
    }
  }
  jo__gif_lzw_write(&state, ent);
  jo__gif_lzw_write(&state, 0x101);
  jo__gif_lzw_write(&state, 0);

  if (state.idx) {
    putc(state.idx, fp);
    fwrite(state.buf, state.idx, 1, fp);
  }
}
extern int
jo_gif_start(jo_gif_t* gif, const char* filename, short w, short h,
             short repeat, int col_cnt) {
  gif->col_cnt = col_cnt > 256 ? 256 : col_cnt < 2 ? 2 : col_cnt;
  gif->pal_siz = jo__gif_ilog2(gif->col_cnt - 1);
  gif->w = w, gif->h = h;
  gif->repeat = repeat;
  if (filename) {
    gif->fp = fopen(filename, "wb");
    if (!gif->fp) {
      printf("Error: Could not open gif output file %s\n", filename);
      return -1;
    }
  } else {
    assert(gif->fp);
    if (!gif->fp) {
      printf("Missing gif output file pointer!\n");
      return -2;
    }
  }
  /* frames for delta buffer */
  gif->buf[0] = malloc((size_t)(w * h));
  gif->buf[1] = (gif->col_cnt != 256) ? malloc((size_t)(w * h)) : 0;
  fwrite("GIF89a", 6, 1, gif->fp);
  /* Logical Screen Descriptor */
  fwrite(&gif->w, 2, 1, gif->fp);
  fwrite(&gif->h, 2, 1, gif->fp);
  putc(0xF0 | gif->pal_siz, gif->fp);
  fwrite("\x00\x00", 2, 1, gif->fp);  /* bg color index, aspect ratio */
  return 0;
}
static void
jo__gif_dither(unsigned char* dith_pix, int w, const unsigned char* pal_pix,
               const unsigned char* pal, int size, int k) {
  /* Floyd-Steinberg Error Diffusion */
  /* TODO: Use something better -- http://caca.zoy.org/study/part3.html */
  int i, diff[3];
  diff[0] = dith_pix[k + 0] - pal[pal_pix[k / 4] * 3 + 0];
  diff[1] = dith_pix[k + 1] - pal[pal_pix[k / 4] * 3 + 1];
  diff[2] = dith_pix[k + 2] - pal[pal_pix[k / 4] * 3 + 2];

  if (k + 4 < size * 4) {
    dith_pix[k + 4 + 0] = (unsigned char)jo__gif_clamp(
        dith_pix[k + 4 + 0] + (diff[0] * 7 / 16), 0, 255);
    dith_pix[k + 4 + 1] = (unsigned char)jo__gif_clamp(
        dith_pix[k + 4 + 1] + (diff[1] * 7 / 16), 0, 255);
    dith_pix[k + 4 + 2] = (unsigned char)jo__gif_clamp(
        dith_pix[k + 4 + 2] + (diff[2] * 7 / 16), 0, 255);
  }
  if (k + w * 4 + 4 < size * 4) {
    for (i = 0; i < 3; ++i) {
      dith_pix[k - 4 + w * 4 + i] = (unsigned char)jo__gif_clamp(
          dith_pix[k - 4 + w * 4 + i] + (diff[i] * 3 / 16), 0, 255);
      dith_pix[k + w * 4 + i] = (unsigned char)jo__gif_clamp(
          dith_pix[k + w * 4 + i] + (diff[i] * 5 / 16), 0, 255);
      dith_pix[k + w * 4 + 4 + i] = (unsigned char)jo__gif_clamp(
          dith_pix[k + w * 4 + 4 + i] + (diff[i] * 1 / 16), 0, 255);
    }
  }
}
extern void
jo_gif_frame(jo_gif_t* gif, unsigned char* rgba, short delayCsec, int gen_pal) {
  int k, i;
  int size = gif->w * gif->h;
  unsigned char local_pal_tbl[0x300];
  unsigned char* pal = gif->frame == 0 || !gen_pal ? gif->pal : local_pal_tbl;
  int no_green = gif->col_cnt == 256 || gif->frame == 0 || gen_pal;
  unsigned char *pal_pix = gif->buf[0], *oldPixels = gif->buf[1];
  if (!gif->fp) {
    return;
  }
  if (gen_pal) {
    unsigned char *map = pal + 3 * (gif->col_cnt != 256);
    jo__gif_quantize(rgba, size * 4, 1, map, gif->col_cnt);
  }
  if (gif->frame == 0 || gen_pal) {
    int n = gif->col_cnt + (gif->col_cnt != 256);
    for (i = (gif->col_cnt != 256); i < n; i++) {
      gif->kd[i].pnt[0] = pal[i * 3 + 0];
      gif->kd[i].pnt[1] = pal[i * 3 + 1];
      gif->kd[i].pnt[2] = pal[i * 3 + 2];
      gif->kd[i].idx = (unsigned char)i;
    }
    gif->kdroot =  jo__gif_kd_build(gif->kd, (gif->col_cnt != 256), gif->col_cnt, 0);
  }
  if (no_green && oldPixels) {
    /* if using all colors, on the first frame, or new palette, can't delta-frame */
    memset(oldPixels, 0, (size_t)size);
  }
  {
    unsigned char* dith_pix = (unsigned char*)malloc((size_t)size * 4);
    memcpy(dith_pix, rgba, (size_t)size * 4);
    for (k = 0; k < size * 4; k += 4) {
      jo_gif_kd_near_t at = {0};
      at.dist = 0x7FFFFFFF;
      at.node = -1;
      jo__gif_kd_nearest(&at, gif->kd, gif->kdroot, &dith_pix[k], 0);
      
      pal_pix[k / 4] = gif->kd[at.node].idx;
      if (gif->dither) {
        jo__gif_dither(dith_pix, gif->w, pal_pix, pal, size, k);
      }
    }
    free(dith_pix);
  }
  if (gif->gamma_f) {
    gif->gamma_f(pal + 3 * (gif->col_cnt != 256), gif->col_cnt * 3);
  }
  if (!gif->frame) {
    /* Global Color Table */
    fwrite(pal, 3 * (1 << (gif->pal_siz + 1)), 1, gif->fp);
    if (gif->repeat >= 0) {
      /* Netscape Extension */
      fwrite("\x21\xff\x0bNETSCAPE2.0\x03\x01", 16, 1, gif->fp);
      /* loop count (extra iterations, 0=repeat forever) */
      fwrite(&gif->repeat, 2, 1, gif->fp);
      putc(0, gif->fp);  /* block terminator */
    }
  }
  /* calculate delta-frame */
  for (i = 0; i < size * !no_green; i++) {
    if (pal_pix[i] == oldPixels[i]) {
      pal_pix[i] = 0;
    } else {
      oldPixels[i] = pal_pix[i];
    }
  }
  /* Graphic Control Extension */
  /* last byte here tells whether to use a transparent colour and delta frame */
  fwrite(no_green ? "\x21\xf9\x04\x00" : "\x21\xf9\x04\x05", 4, 1, gif->fp);
  fwrite(&delayCsec, 2, 1, gif->fp);  /* delayCsec x 1/100 sec */
  fwrite("\x00\x00", 2, 1, gif->fp);  /* transparent color index (first byte) */
  /* Image Descriptor */
  fwrite("\x2c\x00\x00\x00\x00", 5, 1, gif->fp);  /* header, x,y */
  fwrite(&gif->w, 2, 1, gif->fp);
  fwrite(&gif->h, 2, 1, gif->fp);
  if (!gif->frame || !gen_pal) {
    putc(0, gif->fp);
  } else {
    putc(0x80 | gif->pal_siz, gif->fp);
    fwrite(pal, 3 * (1 << (gif->pal_siz + 1)), 1, gif->fp);
  }
  putc(8, gif->fp);  /* block terminator */
  jo__gif_lzw_encode(pal_pix, size, gif->fp);
  putc(0, gif->fp);  /* block terminator */
  ++gif->frame;
}
extern void
jo_gif_end(jo_gif_t* gif) {
  free(gif->buf[0]);
  free(gif->buf[1]);
  if (gif->fp) {
    /* gif trailer */
    putc(0x3b, gif->fp);
    fclose(gif->fp);
  }
}