/lib/cximage-6.0/jbig/jbig.c

Large files files are truncated, but you can click here to view the full file

/*
 *  Portable Free JBIG image compression library
 *
 *  Markus Kuhn -- http://www.cl.cam.ac.uk/~mgk25/
 *
 *  $Id: jbig.c,v 1.22 2004-06-11 15:17:06+01 mgk25 Exp $
 *
 *  This module implements a portable standard C encoder and decoder
 *  using the JBIG lossless bi-level image compression algorithm as
 *  specified in International Standard ISO 11544:1993 or equivalently
 *  as specified in ITU-T Recommendation T.82. See the file jbig.doc
 *  for usage instructions and application examples.
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program 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 General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 * 
 *  If you want to use this program under different license conditions,
 *  then contact the author for an arrangement.
 *
 *  It is possible that certain products which can be built using this
 *  software module might form inventions protected by patent rights in
 *  some countries (e.g., by patents about arithmetic coding algorithms
 *  owned by IBM and AT&T in the USA). Provision of this software by the
 *  author does NOT include any licences for any patents. In those
 *  countries where a patent licence is required for certain applications
 *  of this software module, you will have to obtain such a licence
 *  yourself.
 */

#ifdef DEBUG
#include <stdio.h>
#else
#ifndef NDEBUG
#define NDEBUG
#endif
#endif

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

#include "jbig.h"


/* optional export of arithmetic coder functions for test purposes */
#ifdef TEST_CODEC
#define ARITH
#define ARITH_INL
#else
#define ARITH      static
#ifdef __GNUC__
#define ARITH_INL  static __inline__
#else
#define ARITH_INL  static
#endif
#endif

#define MX_MAX  127    /* maximal supported mx offset for
			* adaptive template in the encoder */

#define TPB2CX  0x195  /* contexts for TP special pixels */
#define TPB3CX  0x0e5
#define TPDCX   0xc3f

/* marker codes */
#define MARKER_STUFF    0x00
#define MARKER_RESERVE  0x01
#define MARKER_SDNORM   0x02
#define MARKER_SDRST    0x03
#define MARKER_ABORT    0x04
#define MARKER_NEWLEN   0x05
#define MARKER_ATMOVE   0x06
#define MARKER_COMMENT  0x07
#define MARKER_ESC      0xff

/* loop array indices */
#define STRIPE  0
#define LAYER   1
#define PLANE   2

/* special jbg_buf pointers (instead of NULL) */
#define SDE_DONE ((struct jbg_buf *) -1)
#define SDE_TODO ((struct jbg_buf *) 0)

/* object code version id */

const char jbg_version[] = 
" JBIG-KIT " JBG_VERSION " -- Markus Kuhn -- "
"$Id: jbig.c,v 1.22 2004-06-11 15:17:06+01 mgk25 Exp $ ";

/*
 * the following array specifies for each combination of the 3
 * ordering bits, which ii[] variable represents which dimension
 * of s->sde.
 */
static const int iindex[8][3] = {
  { 2, 1, 0 },    /* no ordering bit set */
  { -1, -1, -1},  /* SMID -> illegal combination */
  { 2, 0, 1 },    /* ILEAVE */
  { 1, 0, 2 },    /* SMID + ILEAVE */
  { 0, 2, 1 },    /* SEQ */
  { 1, 2, 0 },    /* SEQ + SMID */
  { 0, 1, 2 },    /* SEQ + ILEAVE */
  { -1, -1, -1 }  /* SEQ + SMID + ILEAVE -> illegal combination */
};


/*
 * Array [language][message] with text string error messages that correspond
 * to return values from public functions in this library.
 */
#define NEMSG         9  /* number of error codes */
#define NEMSG_LANG    3  /* number of supported languages */
static const char *errmsg[NEMSG_LANG][NEMSG] = {
  /* English (JBG_EN) */
  {
    "Everything is ok",                                     /* JBG_EOK */
    "Reached specified maximum size",                       /* JBG_EOK_INTR */
    "Unexpected end of data",                               /* JBG_EAGAIN */
    "Not enough memory available",                          /* JBG_ENOMEM */
    "ABORT marker found",                                   /* JBG_EABORT */
    "Unknown marker segment encountered",                   /* JBG_EMARKER */
    "Incremental BIE does not fit to previous one",         /* JBG_ENOCONT */
    "Invalid data encountered",                             /* JBG_EINVAL */
    "Unimplemented features used"                           /* JBG_EIMPL */
  },
  /* German (JBG_DE_8859_1) */
  {
    "Kein Problem aufgetreten",                             /* JBG_EOK */
    "Angegebene maximale Bildgr\366\337e erreicht",         /* JBG_EOK_INTR */
    "Unerwartetes Ende der Daten",                          /* JBG_EAGAIN */
    "Nicht gen\374gend Speicher vorhanden",                 /* JBG_ENOMEM */
    "Es wurde eine Abbruch-Sequenz gefunden",               /* JBG_EABORT */
    "Eine unbekannte Markierungssequenz wurde gefunden",    /* JBG_EMARKER */
    "Neue Daten passen nicht zu vorangegangenen Daten",     /* JBG_ENOCONT */
    "Es wurden ung\374ltige Daten gefunden",                /* JBG_EINVAL */
    "Noch nicht implementierte Optionen wurden benutzt"     /* JBG_EIMPL */
  },
  /* German (JBG_DE_UTF_8) */
  {
    "Kein Problem aufgetreten",                             /* JBG_EOK */
    "Angegebene maximale Bildgr\303\266\303\237e erreicht", /* JBG_EOK_INTR */
    "Unerwartetes Ende der Daten",                          /* JBG_EAGAIN */
    "Nicht gen\303\274gend Speicher vorhanden",             /* JBG_ENOMEM */
    "Es wurde eine Abbruch-Sequenz gefunden",               /* JBG_EABORT */
    "Eine unbekannte Markierungssequenz wurde gefunden",    /* JBG_EMARKER */
    "Neue Daten passen nicht zu vorangegangenen Daten",     /* JBG_ENOCONT */
    "Es wurden ung\303\274ltige Daten gefunden",            /* JBG_EINVAL */
    "Noch nicht implementierte Optionen wurden benutzt"     /* JBG_EIMPL */
  }
};



/*
 * The following three functions are the only places in this code, were
 * C library memory management functions are called. The whole JBIG
 * library has been designed in order to allow multi-threaded
 * execution. No static or global variables are used, so all fuctions
 * are fully reentrant. However if you want to use this multi-thread
 * capability and your malloc, realloc and free are not reentrant,
 * then simply add the necessary semaphores or mutex primitives below.
 * In contrast to C's malloc() and realloc(), but like C's calloc(),
 * these functions take two parameters nmemb and size that are multiplied
 * before being passed on to the corresponding C function. 
 * This we can catch all overflows during a size_t multiplication a
 * a single place.
 */

#ifndef SIZE_MAX
#define SIZE_MAX ((size_t) -1)     /* largest value of size_t */
#endif

static void *checked_malloc(size_t nmemb, size_t size)
{
  void *p;

  /* Full manual exception handling is ugly here for performance
   * reasons. If an adequate handling of lack of memory is required,
   * then use C++ and throw a C++ exception instead of abort(). */

  /* assert that nmemb * size <= SIZE_MAX */
  if (size > SIZE_MAX / nmemb)
    abort();
  
  p = malloc(nmemb * size);

  if (!p)
    abort();

#if 0
  fprintf(stderr, "%p = malloc(%lu * %lu)\n", p,
	  (unsigned long) nmemb, (unsigned long) size);
#endif

  return p;
}


static void *checked_realloc(void *ptr, size_t nmemb, size_t size)
{
  void *p;

  /* Full manual exception handling is ugly here for performance
   * reasons. If an adequate handling of lack of memory is required,
   * then use C++ and throw a C++ exception here instead of abort(). */

  /* assert that nmemb * size <= SIZE_MAX */
  if (size > SIZE_MAX / nmemb)
    abort();
  
  p = realloc(ptr, nmemb * size);

  if (!p)
    abort();

#if 0
  fprintf(stderr, "%p = realloc(%p, %lu * %lu)\n", p, ptr,
	  (unsigned long) nmemb, (unsigned long) size);
#endif

  return p;
}


static void checked_free(void *ptr)
{
  free(ptr);

#if 0
  fprintf(stderr, "free(%p)\n", ptr);
#endif

}



/*
 * The next functions implement the arithmedic encoder and decoder
 * required for JBIG. The same algorithm is also used in the arithmetic
 * variant of JPEG.
 */

#ifdef DEBUG
static long encoded_pixels = 0;
#endif

ARITH void arith_encode_init(struct jbg_arenc_state *s, int reuse_st)
{
  int i;
  
  if (!reuse_st)
    for (i = 0; i < 4096; s->st[i++] = 0);
  s->c = 0;
  s->a = 0x10000L;
  s->sc = 0;
  s->ct = 11;
  s->buffer = -1;    /* empty */
  
  return;
}


ARITH void arith_encode_flush(struct jbg_arenc_state *s)
{
  unsigned long temp;

#ifdef DEBUG
  fprintf(stderr, "  encoded pixels = %ld, a = %05lx, c = %08lx\n",
	  encoded_pixels, s->a, s->c);
#endif

  /* find the s->c in the coding interval with the largest
   * number of trailing zero bits */
  if ((temp = (s->a - 1 + s->c) & 0xffff0000L) < s->c)
    s->c = temp + 0x8000;
  else
    s->c = temp;
  /* send remaining bytes to output */
  s->c <<= s->ct;
  if (s->c & 0xf8000000L) {
    /* one final overflow has to be handled */
    if (s->buffer >= 0) {
      s->byte_out(s->buffer + 1, s->file);
      if (s->buffer + 1 == MARKER_ESC)
	s->byte_out(MARKER_STUFF, s->file);
    }
    /* output 0x00 bytes only when more non-0x00 will follow */
    if (s->c & 0x7fff800L)
      for (; s->sc; --s->sc)
	s->byte_out(0x00, s->file);
  } else {
    if (s->buffer >= 0)
      s->byte_out(s->buffer, s->file); 
    /* T.82 figure 30 says buffer+1 for the above line! Typo? */
    for (; s->sc; --s->sc) {
      s->byte_out(0xff, s->file);
      s->byte_out(MARKER_STUFF, s->file);
    }
  }
  /* output final bytes only if they are not 0x00 */
  if (s->c & 0x7fff800L) {
    s->byte_out((s->c >> 19) & 0xff, s->file);
    if (((s->c >> 19) & 0xff) == MARKER_ESC)
      s->byte_out(MARKER_STUFF, s->file);
    if (s->c & 0x7f800L) {
      s->byte_out((s->c >> 11) & 0xff, s->file);
      if (((s->c >> 11) & 0xff) == MARKER_ESC)
	s->byte_out(MARKER_STUFF, s->file);
    }
  }

  return;
}


ARITH_INL void arith_encode(struct jbg_arenc_state *s, int cx, int pix) 
{
  extern short jbg_lsz[];
  extern unsigned char jbg_nmps[], jbg_nlps[];
  register unsigned lsz, ss;
  register unsigned char *st;
  long temp;

#ifdef DEBUG
  ++encoded_pixels;
#endif

  assert(cx >= 0 && cx < 4096);
  st = s->st + cx;
  ss = *st & 0x7f;
  assert(ss < 113);
  lsz = jbg_lsz[ss];

#if 0
  fprintf(stderr, "pix = %d, cx = %d, mps = %d, st = %3d, lsz = 0x%04x, "
	  "a = 0x%05lx, c = 0x%08lx, ct = %2d, buf = 0x%02x\n",
	  pix, cx, !!(s->st[cx] & 0x80), ss, lsz, s->a, s->c, s->ct,
	  s->buffer);
#endif

  if (((pix << 7) ^ s->st[cx]) & 0x80) {
    /* encode the less probable symbol */
    if ((s->a -= lsz) >= lsz) {
      /* If the interval size (lsz) for the less probable symbol (LPS)
       * is larger than the interval size for the MPS, then exchange
       * the two symbols for coding efficiency, otherwise code the LPS
       * as usual: */
      s->c += s->a;
      s->a = lsz;
    }
    /* Check whether MPS/LPS exchange is necessary
     * and chose next probability estimator status */
    *st &= 0x80;
    *st ^= jbg_nlps[ss];
  } else {
    /* encode the more probable symbol */
    if ((s->a -= lsz) & 0xffff8000L)
      return;   /* A >= 0x8000 -> ready, no renormalization required */
    if (s->a < lsz) {
      /* If the interval size (lsz) for the less probable symbol (LPS)
       * is larger than the interval size for the MPS, then exchange
       * the two symbols for coding efficiency: */
      s->c += s->a;
      s->a = lsz;
    }
    /* chose next probability estimator status */
    *st &= 0x80;
    *st |= jbg_nmps[ss];
  }

  /* renormalization of coding interval */
  do {
    s->a <<= 1;
    s->c <<= 1;
    --s->ct;
    if (s->ct == 0) {
      /* another byte is ready for output */
      temp = s->c >> 19;
      if (temp & 0xffffff00L) {
	/* handle overflow over all buffered 0xff bytes */
	if (s->buffer >= 0) {
	  ++s->buffer;
	  s->byte_out(s->buffer, s->file);
	  if (s->buffer == MARKER_ESC)
	    s->byte_out(MARKER_STUFF, s->file);
	}
	for (; s->sc; --s->sc)
	  s->byte_out(0x00, s->file);
	s->buffer = temp & 0xff;  /* new output byte, might overflow later */
	assert(s->buffer != 0xff);
	/* can s->buffer really never become 0xff here? */
      } else if (temp == 0xff) {
	/* buffer 0xff byte (which might overflow later) */
	++s->sc;
      } else {
	/* output all buffered 0xff bytes, they will not overflow any more */
	if (s->buffer >= 0)
	  s->byte_out(s->buffer, s->file);
	for (; s->sc; --s->sc) {
	  s->byte_out(0xff, s->file);
	  s->byte_out(MARKER_STUFF, s->file);
	}
	s->buffer = temp;   /* buffer new output byte (can still overflow) */
      }
      s->c &= 0x7ffffL;
      s->ct = 8;
    }
  } while (s->a < 0x8000);
 
  return;
}


ARITH void arith_decode_init(struct jbg_ardec_state *s, int reuse_st)
{
  int i;
  
  if (!reuse_st)
    for (i = 0; i < 4096; s->st[i++] = 0);
  s->c = 0;
  s->a = 1;
  s->ct = 0;
  s->result = JBG_OK;
  s->startup = 1;
  return;
}


ARITH_INL int arith_decode(struct jbg_ardec_state *s, int cx)
{
  extern short jbg_lsz[];
  extern unsigned char jbg_nmps[], jbg_nlps[];
  register unsigned lsz, ss;
  register unsigned char *st;
  int pix;

  /* renormalization */
  while (s->a < 0x8000 || s->startup) {
    if (s->ct < 1 && s->result != JBG_READY) {
      /* first we have to move a new byte into s->c */
      if (s->pscd_ptr >= s->pscd_end) {
	s->result = JBG_MORE;
	return -1;
      }
      if (*s->pscd_ptr == 0xff) 
	if (s->pscd_ptr + 1 >= s->pscd_end) {
	  s->result = JBG_MARKER;
	  return -1;
	} else {
	  if (*(s->pscd_ptr + 1) == MARKER_STUFF) {
	    s->c |= 0xffL << (8 - s->ct);
	    s->ct += 8;
	    s->pscd_ptr += 2;
	    s->result = JBG_OK;
	  } else
	    s->result = JBG_READY;
	}
      else {
	s->c |= (long)*(s->pscd_ptr++) << (8 - s->ct);
	s->ct += 8;
	s->result = JBG_OK;
      }
    }
    s->c <<= 1;
    s->a <<= 1;
    --s->ct;
    if (s->a == 0x10000L)
      s->startup = 0;
  }

  st = s->st + cx;
  ss = *st & 0x7f;
  assert(ss < 113);
  lsz = jbg_lsz[ss];

#if 0
  fprintf(stderr, "cx = %d, mps = %d, st = %3d, lsz = 0x%04x, a = 0x%05lx, "
	  "c = 0x%08lx, ct = %2d\n",
	  cx, !!(s->st[cx] & 0x80), ss, lsz, s->a, s->c, s->ct);
#endif

  if ((s->c >> 16) < (s->a -= lsz))
    if (s->a & 0xffff8000L)
      return *st >> 7;
    else {
      /* MPS_EXCHANGE */
      if (s->a < lsz) {
	pix = 1 - (*st >> 7);
	/* Check whether MPS/LPS exchange is necessary
	 * and chose next probability estimator status */
	*st &= 0x80;
	*st ^= jbg_nlps[ss];
      } else {
	pix = *st >> 7;
	*st &= 0x80;
	*st |= jbg_nmps[ss];
      }
    }
  else {
    /* LPS_EXCHANGE */
    if (s->a < lsz) {
      s->c -= s->a << 16;
      s->a = lsz;
      pix = *st >> 7;
      *st &= 0x80;
      *st |= jbg_nmps[ss];
    } else {
      s->c -= s->a << 16;
      s->a = lsz;
      pix = 1 - (*st >> 7);
      /* Check whether MPS/LPS exchange is necessary
       * and chose next probability estimator status */
      *st &= 0x80;
      *st ^= jbg_nlps[ss];
    }
  }

  return pix;
}



/*
 * Memory management for buffers which are used for temporarily
 * storing SDEs by the encoder.
 *
 * The following functions manage a set of struct jbg_buf storage
 * containers were each can keep JBG_BUFSIZE bytes. The jbg_buf
 * containers can be linked to form linear double-chained lists for
 * which a number of operations are provided. Blocks which are
 * tempoarily not used any more are returned to a freelist which each
 * encoder keeps. Only the destructor of the encoder actually returns
 * the block via checked_free() to the stdlib memory management.
 */


/*
 * Allocate a new buffer block and initialize it. Try to get it from
 * the free_list, and if it is empty, call checked_malloc().
 */
static struct jbg_buf *jbg_buf_init(struct jbg_buf **free_list)
{
  struct jbg_buf *new_block;
  
  /* Test whether a block from the free list is available */
  if (*free_list) {
    new_block = *free_list;
    *free_list = new_block->next;
  } else {
    /* request a new memory block */
    new_block = (struct jbg_buf *) checked_malloc(1, sizeof(struct jbg_buf));
  }
  new_block->len = 0;
  new_block->next = NULL;
  new_block->previous = NULL;
  new_block->last = new_block;
  new_block->free_list = free_list;

  return new_block;
}


/*
 * Return an entire free_list to the memory management of stdlib.
 * This is only done by jbg_enc_free().
 */
static void jbg_buf_free(struct jbg_buf **free_list)
{
  struct jbg_buf *tmp;
  
  while (*free_list) {
    tmp = (*free_list)->next;
    checked_free(*free_list);
    *free_list = tmp;
  }
  
  return;
}


/*
 * Append a single byte to a single list that starts with the block
 * *(struct jbg_buf *) head. The type of *head is void here in order to
 * keep the interface of the arithmetic encoder gereric, which uses this
 * function as a call-back function in order to deliver single bytes
 * for a PSCD.
 */
static void jbg_buf_write(int b, void *head)
{
  struct jbg_buf *now;

  now = ((struct jbg_buf *) head)->last;
  if (now->len < JBG_BUFSIZE - 1) {
    now->d[now->len++] = b;
    return;
  }
  now->next = jbg_buf_init(((struct jbg_buf *) head)->free_list);
  now->next->previous = now;
  now->next->d[now->next->len++] = b;
  ((struct jbg_buf *) head)->last = now->next;

  return;
}


/*
 * Remove any trailing zero bytes from the end of a linked jbg_buf list,
 * however make sure that no zero byte is removed which directly
 * follows a 0xff byte (i.e., keep MARKER_ESC MARKER_STUFF sequences
 * intact). This function is used to remove any redundant final zero
 * bytes from a PSCD.
 */
static void jbg_buf_remove_zeros(struct jbg_buf *head)
{
  struct jbg_buf *last;

  while (1) {
    /* remove trailing 0x00 in last block of list until this block is empty */
    last = head->last;
    while (last->len && last->d[last->len - 1] == 0)
      last->len--;
    /* if block became really empty, remove it in case it is not the
     * only remaining block and then loop to next block */
    if (last->previous && !last->len) {
      head->last->next = *head->free_list;
      *head->free_list = head->last;
      head->last = last->previous;
      head->last->next = NULL;
    } else
      break;
  }

  /*
   * If the final non-zero byte is 0xff (MARKER_ESC), then we just have
   * removed a MARKER_STUFF and we will append it again now in order
   * to preserve PSCD status of byte stream.
   */
  if (head->last->len && head->last->d[head->last->len - 1] == MARKER_ESC)
    jbg_buf_write(MARKER_STUFF, head);
 
  return;
}


/*
 * The jbg_buf list which starts with block *new_prefix is concatenated
 * with the list which starts with block **start and *start will then point
 * to the first block of the new list.
 */
static void jbg_buf_prefix(struct jbg_buf *new_prefix, struct jbg_buf **start)
{
  new_prefix->last->next = *start;
  new_prefix->last->next->previous = new_prefix->last;
  new_prefix->last = new_prefix->last->next->last;
  *start = new_prefix;
  
  return;
}


/*
 * Send the contents of a jbg_buf list that starts with block **head to
 * the call back function data_out and return the blocks of the jbg_buf
 * list to the freelist from which these jbg_buf blocks have been taken.
 * After the call, *head == NULL.
 */
static void jbg_buf_output(struct jbg_buf **head,
			void (*data_out)(unsigned char *start,
					 size_t len, void *file),
			void *file)
{
  struct jbg_buf *tmp;
  
  while (*head) {
    data_out((*head)->d, (*head)->len, file);
    tmp = (*head)->next;
    (*head)->next = *(*head)->free_list;
    *(*head)->free_list = *head;
    *head = tmp;
  }
  
  return;
}


/*
 * Calculate y = ceil(x/2) applied n times, which is equivalent to
 * y = ceil(x/(2^n)). This function is used to
 * determine the number of pixels per row or column after n resolution
 * reductions. E.g. X[d-1] = jbg_ceil_half(X[d], 1) and X[0] =
 * jbg_ceil_half(X[d], d) as defined in clause 6.2.3 of T.82.
 */
unsigned long jbg_ceil_half(unsigned long x, int n)
{
  unsigned long mask;
  
  assert(n >= 0 && n < 32);
  mask = (1UL << n) - 1;     /* the lowest n bits are 1 here */
  return (x >> n) + ((mask & x) != 0);
}


/*
 * Set L0 (the number of lines in a stripe at lowest resolution)
 * to a default value, such that there are about 35 stripes, as
 * suggested in Annex C of ITU-T T.82, without exceeding the
 * limit 128/2^D suggested in Annex A.
 */
static void jbg_set_default_l0(struct jbg_enc_state *s)
{
  s->l0 = jbg_ceil_half(s->yd, s->d) / 35;   /* 35 stripes/image */
  while ((s->l0 << s->d) > 128)              /* but <= 128 lines/stripe */
    --s->l0;
  if (s->l0 < 2) s->l0 = 2;
}


/*
 * Calculate the number of stripes, as defined in clause 6.2.3 of T.82.
 */
static unsigned long jbg_stripes(unsigned long l0, unsigned long yd,
				 unsigned long d)
{
  unsigned long y0 = jbg_ceil_half(yd, d);

  return y0 / l0 + (y0 % l0 != 0);
}


/*
 * Initialize the status struct for the encoder.
 */
void jbg_enc_init(struct jbg_enc_state *s, unsigned long x, unsigned long y,
                  int planes, unsigned char **p,
                  void (*data_out)(unsigned char *start, size_t len,
				   void *file),
		  void *file)
{
  unsigned long l, lx;
  int i;

  extern char jbg_resred[], jbg_dptable[];

  s->xd = x;
  s->yd = y;
  s->yd1 = y; /* This is the hight initially announced in BIH. To provoke
                 generation of NEWLEN for T.85 compatibility tests,
                 overwrite with new value s->yd1 > s->yd  */
  s->planes = planes;
  s->data_out = data_out;
  s->file = file;

  s->d = 0;
  s->dl = 0;
  s->dh = s->d;
  jbg_set_default_l0(s);
  s->mx = 8;
  s->my = 0;
  s->order = JBG_ILEAVE | JBG_SMID;
  s->options = JBG_TPBON | JBG_TPDON | JBG_DPON;
  s->dppriv = jbg_dptable;
  s->res_tab = jbg_resred;
  
  s->highres = (int *) checked_malloc(planes, sizeof(int));
  s->lhp[0] = p;
  s->lhp[1] = (unsigned char **)
    checked_malloc(planes, sizeof(unsigned char *));
  for (i = 0; i < planes; i++) {
    s->highres[i] = 0;
    s->lhp[1][i] = (unsigned char *)
      checked_malloc(jbg_ceil_half(y, 1), jbg_ceil_half(x, 1+3));
  }
  
  s->free_list = NULL;
  s->s = (struct jbg_arenc_state *) 
    checked_malloc(s->planes, sizeof(struct jbg_arenc_state));
  s->tx = (int *) checked_malloc(s->planes, sizeof(int));
  lx = jbg_ceil_half(x, 1);
  s->tp = (char *) checked_malloc(lx, sizeof(char));
  for (l = 0; l < lx; s->tp[l++] = 2);
  s->sde = NULL;

  return;
}


/*
 * This function selects the number of differential layers based on
 * the maximum size requested for the lowest resolution layer. If
 * possible, a number of differential layers is selected, which will
 * keep the size of the lowest resolution layer below or equal to the
 * given width x and height y. However not more than 6 differential
 * resolution layers will be used. In addition, a reasonable value for
 * l0 (height of one stripe in the lowest resolution layer) is
 * selected, which obeys the recommended limitations for l0 in annex A
 * and C of the JBIG standard. The selected number of resolution layers
 * is returned. 
 */
int jbg_enc_lrlmax(struct jbg_enc_state *s, unsigned long x, 
		   unsigned long y)
{
  for (s->d = 0; s->d < 6; s->d++)
    if (jbg_ceil_half(s->xd, s->d) <= x && jbg_ceil_half(s->yd, s->d) <= y)
      break;
  s->dl = 0;
  s->dh = s->d;
  jbg_set_default_l0(s);
  return s->d;
}


/*
 * As an alternative to jbg_enc_lrlmax(), the following function allows
 * to specify the number of layers directly. The stripe height and layer
 * range is also adjusted automatically here.
 */
void jbg_enc_layers(struct jbg_enc_state *s, int d)
{
  if (d < 0 || d > 31)
    return;
  s->d  = d;
  s->dl = 0;
  s->dh = s->d;
  jbg_set_default_l0(s);
  return;
}


/*
 * Specify the highest and lowest resolution layers which will be
 * written to the output file. Call this function not before
 * jbg_enc_layers() or jbg_enc_lrlmax(), because these two functions
 * reset the lowest and highest resolution layer to default values.
 * Negative values are ignored. The total number of layers is returned.
 */
int jbg_enc_lrange(struct jbg_enc_state *s, int dl, int dh)
{
  if (dl >= 0     && dl <= s->d) s->dl = dl;
  if (dh >= s->dl && dh <= s->d) s->dh = dh;

  return s->d;
}


/*
 * The following function allows to specify the bits describing the
 * options of the format as well as the maximum AT movement window and
 * the number of layer 0 lines per stripes.
 */
void jbg_enc_options(struct jbg_enc_state *s, int order, int options,
		     unsigned long l0, int mx, int my)
{
  if (order >= 0 && order <= 0x0f) s->order = order;
  if (options >= 0) s->options = options;
  if (l0 > 0) s->l0 = l0;
  if (mx >= 0 && my < 128) s->mx = mx;
  if (my >= 0 && my < 256) s->my = my;

  return;
}


/*
 * This function actually does all the tricky work involved in producing
 * a SDE, which is stored in the appropriate s->sde[][][] element
 * for later output in the correct order.
 */
static void encode_sde(struct jbg_enc_state *s,
		       long stripe, int layer, int plane)
{
  unsigned char *hp, *lp1, *lp2, *p0, *p1, *q1, *q2;
  unsigned long hl, ll, hx, hy, lx, ly, hbpl, lbpl;
  unsigned long line_h0 = 0, line_h1 = 0;
  unsigned long line_h2, line_h3, line_l1, line_l2, line_l3;
  struct jbg_arenc_state *se;
  unsigned long i, j, y;
  long o;
  unsigned a, p, t;
  int ltp, ltp_old, cx;
  unsigned long c_all, c[MX_MAX + 1], cmin, cmax, clmin, clmax;
  int tmax, at_determined;
  int new_tx;
  long new_tx_line = -1;
  struct jbg_buf *new_jbg_buf;

#ifdef DEBUG
  static long tp_lines, tp_exceptions, tp_pixels, dp_pixels;
  static long encoded_pixels;
#endif

  /* return immediately if this stripe has already been encoded */
  if (s->sde[stripe][layer][plane] != SDE_TODO)
    return;

#ifdef DEBUG
  if (stripe == 0)
    tp_lines = tp_exceptions = tp_pixels = dp_pixels = encoded_pixels = 0;
  fprintf(stderr, "encode_sde: s/d/p = %2ld/%2d/%2d\n",
	  stripe, layer, plane);
#endif

  /* number of lines per stripe in highres image */
  hl = s->l0 << layer;
  /* number of lines per stripe in lowres image */
  ll = hl >> 1;
  /* current line number in highres image */
  y = stripe * hl;
  /* number of pixels in highres image */
  hx = jbg_ceil_half(s->xd, s->d - layer);
  hy = jbg_ceil_half(s->yd, s->d - layer);
  /* number of pixels in lowres image */
  lx = jbg_ceil_half(hx, 1);
  ly = jbg_ceil_half(hy, 1);
  /* bytes per line in highres and lowres image */
  hbpl = jbg_ceil_half(hx, 3);
  lbpl = jbg_ceil_half(lx, 3);
  /* pointer to first image byte of highres stripe */
  hp = s->lhp[s->highres[plane]][plane] + stripe * hl * hbpl;
  lp2 = s->lhp[1 - s->highres[plane]][plane] + stripe * ll * lbpl;
  lp1 = lp2 + lbpl;
  
  /* initialize arithmetic encoder */
  se = s->s + plane;
  arith_encode_init(se, stripe != 0);
  s->sde[stripe][layer][plane] = jbg_buf_init(&s->free_list);
  se->byte_out = jbg_buf_write;
  se->file = s->sde[stripe][layer][plane];

  /* initialize adaptive template movement algorithm */
  c_all = 0;
  for (t = 0; t <= s->mx; t++)
    c[t] = 0;
  if (stripe == 0)
    s->tx[plane] = 0;
  new_tx = -1;
  at_determined = 0;  /* we haven't yet decided the template move */
  if (s->mx == 0)
    at_determined = 1;

  /* initialize typical prediction */
  ltp = 0;
  if (stripe == 0)
    ltp_old = 0;
  else {
    ltp_old = 1;
    p1 = hp - hbpl;
    if (y > 1) {
      q1 = p1 - hbpl;
      while (p1 < hp && (ltp_old = (*p1++ == *q1++)) != 0);
    } else
      while (p1 < hp && (ltp_old = (*p1++ == 0)) != 0);
  }

  if (layer == 0) {

    /*
     *  Encode lowest resolution layer
     */

    for (i = 0; i < hl && y < hy; i++, y++) {

      /* check whether it is worth to perform an ATMOVE */
      if (!at_determined && c_all > 2048) {
	cmin = clmin = 0xffffffffL;
	cmax = clmax = 0;
	tmax = 0;
	for (t = (s->options & JBG_LRLTWO) ? 5 : 3; t <= s->mx; t++) {
	  if (c[t] > cmax) cmax = c[t];
	  if (c[t] < cmin) cmin = c[t];
	  if (c[t] > c[tmax]) tmax = t;
	}
	clmin = (c[0] < cmin) ? c[0] : cmin;
	clmax = (c[0] > cmax) ? c[0] : cmax;
	if (c_all - cmax < (c_all >> 3) &&
	    cmax - c[s->tx[plane]] > c_all - cmax &&
	    cmax - c[s->tx[plane]] > (c_all >> 4) &&
	    /*                     ^ T.82 said < here, fixed in Cor.1/25 */
	    cmax - (c_all - c[s->tx[plane]]) > c_all - cmax &&
	    cmax - (c_all - c[s->tx[plane]]) > (c_all >> 4) &&
	    cmax - cmin > (c_all >> 2) &&
	    (s->tx[plane] || clmax - clmin > (c_all >> 3))) {
	  /* we have decided to perform an ATMOVE */
	  new_tx = tmax;
	  if (!(s->options & JBG_DELAY_AT)) {
	    new_tx_line = i;
	    s->tx[plane] = new_tx;
	  }
#ifdef DEBUG
	  fprintf(stderr, "ATMOVE: line=%ld, tx=%d, c_all=%ld\n",
		  i, new_tx, c_all);
#endif
	}
	at_determined = 1;
      }
      assert(s->tx[plane] >= 0); /* i.e., tx can safely be cast to unsigned */
      
      /* typical prediction */
      if (s->options & JBG_TPBON) {
	ltp = 1;
	p1 = hp;
	if (y > 0) {
	  q1 = hp - hbpl;
	  while (q1 < hp && (ltp = (*p1++ == *q1++)) != 0);
	} else
	  while (p1 < hp + hbpl && (ltp = (*p1++ == 0)) != 0);
	arith_encode(se, (s->options & JBG_LRLTWO) ? TPB2CX : TPB3CX,
		     ltp == ltp_old);
#ifdef DEBUG
	tp_lines += ltp;
#endif
	ltp_old = ltp;
	if (ltp) {
	  /* skip next line */
	  hp += hbpl;
	  continue;
	}
      }

      /*
       * Layout of the variables line_h1, line_h2, line_h3, which contain
       * as bits the neighbour pixels of the currently coded pixel X:
       *
       *          76543210765432107654321076543210     line_h3
       *          76543210765432107654321076543210     line_h2
       *  76543210765432107654321X76543210             line_h1
       */
      
      line_h1 = line_h2 = line_h3 = 0;
      if (y > 0) line_h2 = (long)*(hp - hbpl) << 8;
      if (y > 1) line_h3 = (long)*(hp - hbpl - hbpl) << 8;
      
      /* encode line */
      for (j = 0; j < hx; hp++) {
	line_h1 |= *hp;
	if (j < hbpl * 8 - 8 && y > 0) {
	  line_h2 |= *(hp - hbpl + 1);
	  if (y > 1)
	    line_h3 |= *(hp - hbpl - hbpl + 1);
	}
	if (s->options & JBG_LRLTWO) {
	  /* two line template */
	  do {
	    line_h1 <<= 1;  line_h2 <<= 1;  line_h3 <<= 1;
	    if (s->tx[plane]) {
	      if ((unsigned) s->tx[plane] > j)
		a = 0;
	      else {
		o = (j - s->tx[plane]) - (j & ~7L);
		a = (hp[o >> 3] >> (7 - (o & 7))) & 1;
		a <<= 4;
	      }
	      assert(s->tx[plane] > 23 ||
		     a == ((line_h1 >> (4 + s->tx[plane])) & 0x010));
	      arith_encode(se, (((line_h2 >> 10) & 0x3e0) | a |
				((line_h1 >>  9) & 0x00f)),
			   (line_h1 >> 8) & 1);
	    }
	    else
	      arith_encode(se, (((line_h2 >> 10) & 0x3f0) |
				((line_h1 >>  9) & 0x00f)),
			   (line_h1 >> 8) & 1);
#ifdef DEBUG
	    encoded_pixels++;
#endif
	    /* statistics for adaptive template changes */
	    if (!at_determined && j >= s->mx && j < hx-2) {
	      p = (line_h1 & 0x100) != 0; /* current pixel value */
	      c[0] += ((unsigned int)((line_h2 & 0x4000) != 0)) == p; /* default position */
	      assert(!(((line_h2 >> 6) ^ line_h1) & 0x100) ==
		     (((line_h2 & 0x4000) != 0) == p));
	      for (t = 5; t <= s->mx && t <= j; t++) {
		o = (j - t) - (j & ~7L);
		a = (hp[o >> 3] >> (7 - (o & 7))) & 1;
		assert(t > 23 ||
		       (a == p) == !(((line_h1 >> t) ^ line_h1) & 0x100));
		c[t] += a == p;
	      }
	      for (; t <= s->mx; t++) {
		c[t] += 0 == p;
	      }
	      ++c_all;
	    }
	  } while (++j & 7 && j < hx);
	} else {
	  /* three line template */
	  do {
	    line_h1 <<= 1;  line_h2 <<= 1;  line_h3 <<= 1;
	    if (s->tx[plane]) {
	      if ((unsigned) s->tx[plane] > j)
		a = 0;
	      else {
		o = (j - s->tx[plane]) - (j & ~7L);
		a = (hp[o >> 3] >> (7 - (o & 7))) & 1;
		a <<= 2;
	      }
	      assert(s->tx[plane] > 23 ||
		     a == ((line_h1 >> (6 + s->tx[plane])) & 0x004));
	      arith_encode(se, (((line_h3 >>  8) & 0x380) |
				((line_h2 >> 12) & 0x078) | a |
				((line_h1 >>  9) & 0x003)),
			   (line_h1 >> 8) & 1);
	    } else
	      arith_encode(se, (((line_h3 >>  8) & 0x380) |
				((line_h2 >> 12) & 0x07c) |
				((line_h1 >>  9) & 0x003)),
			   (line_h1 >> 8) & 1);
#ifdef DEBUG
	    encoded_pixels++;
#endif
	    /* statistics for adaptive template changes */
	    if (!at_determined && j >= s->mx && j < hx-2) {
	      p = (line_h1 & 0x100) != 0; /* current pixel value */
	      c[0] += ((unsigned int)((line_h2 & 0x4000) != 0)) == p; /* default position */
	      assert(!(((line_h2 >> 6) ^ line_h1) & 0x100) ==
		     (((line_h2 & 0x4000) != 0) == p));
	      for (t = 3; t <= s->mx && t <= j; t++) {
		o = (j - t) - (j & ~7L);
		a = (hp[o >> 3] >> (7 - (o & 7))) & 1;
		assert(t > 23 ||
		       (a == p) == !(((line_h1 >> t) ^ line_h1) & 0x100));
		c[t] += a == p;
	      }
	      for (; t <= s->mx; t++) {
		c[t] += 0 == p;
	      }
	      ++c_all;
	    }
	  } while (++j & 7 && j < hx);
	} /* if (s->options & JBG_LRLTWO) */
      } /* for (j = ...) */
    } /* for (i = ...) */

  } else {

    /*
     *  Encode differential layer
     */
    
    for (i = 0; i < hl && y < hy; i++, y++) {

      /* check whether it is worth to perform an ATMOVE */
      if (!at_determined && c_all > 2048) {
	cmin = clmin = 0xffffffffL;
	cmax = clmax = 0;
	tmax = 0;
	for (t = 3; t <= s->mx; t++) {
	  if (c[t] > cmax) cmax = c[t];
	  if (c[t] < cmin) cmin = c[t];
	  if (c[t] > c[tmax]) tmax = t;
	}
	clmin = (c[0] < cmin) ? c[0] : cmin;
	clmax = (c[0] > cmax) ? c[0] : cmax;
	if (c_all - cmax < (c_all >> 3) &&
	    cmax - c[s->tx[plane]] > c_all - cmax &&
	    cmax - c[s->tx[plane]] > (c_all >> 4) &&
	    /*                     ^ T.82 said < here, fixed in Cor.1/25 */
	    cmax - (c_all - c[s->tx[plane]]) > c_all - cmax &&
	    cmax - (c_all - c[s->tx[plane]]) > (c_all >> 4) &&
	    cmax - cmin > (c_all >> 2) &&
	    (s->tx[plane] || clmax - clmin > (c_all >> 3))) {
	  /* we have decided to perform an ATMOVE */
	  new_tx = tmax;
	  if (!(s->options & JBG_DELAY_AT)) {
	    new_tx_line = i;
	    s->tx[plane] = new_tx;
	  }
#ifdef DEBUG
	  fprintf(stderr, "ATMOVE: line=%ld, tx=%d, c_all=%ld\n",
		  i, new_tx, c_all);
#endif
	}
	at_determined = 1;
      }
      
      if ((i >> 1) >= ll - 1 || (y >> 1) >= ly - 1)
	lp1 = lp2;

      /* typical prediction */
      if (s->options & JBG_TPDON && (i & 1) == 0) {
	q1 = lp1; q2 = lp2;
	p0 = p1 = hp;
	if (i < hl - 1 && y < hy - 1)
	  p0 = hp + hbpl;
	if (y > 1)
	  line_l3 = (long)*(q2 - lbpl) << 8;
	else
	  line_l3 = 0;
	line_l2 = (long)*q2 << 8;
	line_l1 = (long)*q1 << 8;
	ltp = 1;
	for (j = 0; j < lx && ltp; q1++, q2++) {
	  if (j < lbpl * 8 - 8) {
	    if (y > 1)
	      line_l3 |= *(q2 - lbpl + 1);
	    line_l2 |= *(q2 + 1);
	    line_l1 |= *(q1 + 1);
	  }
	  do {
	    if ((j >> 2) < hbpl) {
	      line_h1 = *(p1++);
	      line_h0 = *(p0++);
	    }
	    do {
	      line_l3 <<= 1;
	      line_l2 <<= 1;
	      line_l1 <<= 1;
	      line_h1 <<= 2;
	      line_h0 <<= 2;
	      cx = (((line_l3 >> 15) & 0x007) |
		    ((line_l2 >> 12) & 0x038) |
		    ((line_l1 >> 9)  & 0x1c0));
	      if (cx == 0x000)
		if ((line_h1 & 0x300) == 0 && (line_h0 & 0x300) == 0)
		  s->tp[j] = 0;
		else {
		  ltp = 0;
#ifdef DEBUG
		  tp_exceptions++;
#endif
		}
	      else if (cx == 0x1ff)
		if ((line_h1 & 0x300) == 0x300 && (line_h0 & 0x300) == 0x300)
		  s->tp[j] = 1;
		else {
		  ltp = 0;
#ifdef DEBUG
		  tp_exceptions++;
#endif
		}
	      else
		s->tp[j] = 2;
	    } while (++j & 3 && j < lx);
	  } while (j & 7 && j < lx);
	} /* for (j = ...) */
	arith_encode(se, TPDCX, !ltp);
#ifdef DEBUG
	tp_lines += ltp;
#endif
      }


      /*
       * Layout of the variables line_h1, line_h2, line_h3, which contain
       * as bits the high resolution neighbour pixels of the currently coded
       * highres pixel X:
       *
       *            76543210 76543210 76543210 76543210     line_h3
       *            76543210 76543210 76543210 76543210     line_h2
       *   76543210 76543210 7654321X 76543210              line_h1
       *
       * Layout of the variables line_l1, line_l2, line_l3, which contain
       * the low resolution pixels near the currently coded pixel as bits.
       * The lowres pixel in which the currently coded highres pixel is
       * located is marked as Y:
       *
       *            76543210 76543210 76543210 76543210     line_l3
       *            76543210 7654321Y 76543210 76543210     line_l2
       *            76543210 76543210 76543210 76543210     line_l1
       */
      

      line_h1 = line_h2 = line_h3 = line_l1 = line_l2 = line_l3 = 0;
      if (y > 0) line_h2 = (long)*(hp - hbpl) << 8;
      if (y > 1) {
	line_h3 = (long)*(hp - hbpl - hbpl) << 8;
	line_l3 = (long)*(lp2 - lbpl) << 8;
      }
      line_l2 = (long)*lp2 << 8;
      line_l1 = (long)*lp1 << 8;
      
      /* encode line */
      for (j = 0; j < hx; lp1++, lp2++) {
	if ((j >> 1) < lbpl * 8 - 8) {
	  if (y > 1)
	    line_l3 |= *(lp2 - lbpl + 1);
	  line_l2 |= *(lp2 + 1);
	  line_l1 |= *(lp1 + 1);
	}
	do { /* ... while (j & 15 && j < hx) */

	  assert(hp - (s->lhp[s->highres[plane]][plane] +
		       (stripe * hl + i) * hbpl)
		 == (ptrdiff_t) j >> 3);

	  assert(lp2 - (s->lhp[1-s->highres[plane]][plane] +
			(stripe * ll + (i>>1)) * lbpl)
		 == (ptrdiff_t) j >> 4);

	  line_h1 |= *hp;
	  if (j < hbpl * 8 - 8) {
	    if (y > 0) {
	      line_h2 |= *(hp - hbpl + 1);
	      if (y > 1)
		line_h3 |= *(hp - hbpl - hbpl + 1);
	    }
	  }
	  do { /* ... while (j & 7 && j < hx) */
	    line_l1 <<= 1;  line_l2 <<= 1;  line_l3 <<= 1;
	    if (ltp && s->tp[j >> 1] < 2) {
	      /* pixel are typical and have not to be encoded */
	      line_h1 <<= 2;  line_h2 <<= 2;  line_h3 <<= 2;
#ifdef DEBUG
	      do {
		++tp_pixels;
	      } while (++j & 1 && j < hx);
#else
	      j += 2;
#endif
	    } else
	      do { /* ... while (++j & 1 && j < hx) */
		line_h1 <<= 1;  line_h2 <<= 1;  line_h3 <<= 1;

		/* deterministic prediction */
		if (s->options & JBG_DPON) {
		  if ((y & 1) == 0) {
		    if ((j & 1) == 0) {
		      /* phase 0 */
		      if (s->dppriv[((line_l3 >> 16) & 0x003) |
				    ((line_l2 >> 14) & 0x00c) |
				    ((line_h1 >> 5)  & 0x010) |
				    ((line_h2 >> 10) & 0x0e0)] < 2) {
#ifdef DEBUG
			++dp_pixels;
#endif
			continue;
		      }
		    } else {
		      /* phase 1 */
		      if (s->dppriv[(((line_l3 >> 16) & 0x003) |
				     ((line_l2 >> 14) & 0x00c) |
				     ((line_h1 >> 5)  & 0x030) |
				     ((line_h2 >> 10) & 0x1c0)) + 256] < 2) {
#ifdef DEBUG
			++dp_pixels;
#endif
			continue;
		      }
		    }
		  } else {
		    if ((j & 1) == 0) {
		      /* phase 2 */
		      if (s->dppriv[(((line_l3 >> 16) & 0x003) |
				     ((line_l2 >> 14) & 0x00c) |
				     ((line_h1 >> 5)  & 0x010) |
				     ((line_h2 >> 10) & 0x0e0) |
				     ((line_h3 >> 7) & 0x700)) + 768] < 2) {
#ifdef DEBUG
			++dp_pixels;
#endif
			continue;
		      }
		    } else {
		      /* phase 3 */
		      if (s->dppriv[(((line_l3 >> 16) & 0x003) |
				     ((line_l2 >> 14) & 0x00c) |
				     ((line_h1 >> 5)  & 0x030) |
				     ((line_h2 >> 10) & 0x1c0) |
				     ((line_h3 >> 7)  & 0xe00)) + 2816] < 2) {
#ifdef DEBUG
			++dp_pixels;
#endif
			continue;
		      }
		    }	
		  }	
		}

		/* determine context */
		if (s->tx[plane]) {
		  if ((unsigned) s->tx[plane] > j)
		    a = 0;
		  else {
		    o = (j - s->tx[plane]) - (j & ~7L);
		    a = (hp[o >> 3] >> (7 - (o & 7))) & 1;
		    a <<= 4;
		  }
		  assert(s->tx[plane] > 23 ||
			 a == ((line_h1 >> (4 + s->tx[plane])) & 0x010));
		  cx = (((line_h1 >> 9)  & 0x003) | a |
			((line_h2 >> 13) & 0x00c) |
			((line_h3 >> 11) & 0x020));
		} else
		  cx = (((line_h1 >> 9)  & 0x003) |
			((line_h2 >> 13) & 0x01c) |
			((line_h3 >> 11) & 0x020));
		if (j & 1)
		  cx |= (((line_l2 >> 9)  & 0x0c0) |
			 ((line_l1 >> 7)  & 0x300)) | (1UL << 10);
		else
		  cx |= (((line_l2 >> 10) & 0x0c0) |
			 ((line_l1 >> 8)  & 0x300));
		cx |= (y & 1) << 11;

		arith_encode(se, cx, (line_h1 >> 8) & 1);
#ifdef DEBUG
		encoded_pixels++;
#endif
		
		/* statistics for adaptive template changes */
		if (!at_determined && j >= s->mx) {
		  c[0] += !(((line_h2 >> 6) ^ line_h1) & 0x100);
		  for (t = 3; t <= s->mx; t++)
		    c[t] += !(((line_h1 >> t) ^ line_h1) & 0x100);
		  ++c_all;
		}
		
	      } while (++j & 1 && j < hx);
	  } while (j & 7 && j < hx);
	  hp++;
	} while (j & 15 && j < hx);
      } /* for (j = ...) */

      /* low resolution pixels are used twice */
      if ((i & 1) == 0) {
	lp1 -= lbpl;
	lp2 -= lbpl;
      }
      
    } /* for (i = ...) */
  }
  
  arith_encode_flush(se);
  jbg_buf_remove_zeros(s->sde[stripe][layer][plane]);
  jbg_buf_write(MARKER_ESC, s->sde[stripe][layer][plane]);
  jbg_buf_write(MARKER_SDNORM, s->sde[stripe][layer][plane]);

  /* add ATMOVE */
  if (new_tx != -1) {
    if (s->options & JBG_DELAY_AT) {
      /* ATMOVE will become active at the first line of the next stripe */
      s->tx[plane] = new_tx;
      jbg_buf_write(MARKER_ESC, s->sde[stripe][layer][plane]);
      jbg_buf_write(MARKER_ATMOVE, s->sde[stripe][layer][plane]);
      jbg_buf_write(0, s->sde[stripe][layer][plane]);
      jbg_buf_write(0, s->sde[stripe][layer][plane]);
      jbg_buf_write(0, s->sde[stripe][layer][plane]);
      jbg_buf_write(0, s->sde[stripe][layer][plane]);
      jbg_buf_write(s->tx[plane], s->sde[stripe][layer][plane]);
      jbg_buf_write(0, s->sde[stripe][layer][plane]);
    } else {
      /* ATMOVE has already become active during this stripe
       * => we have to prefix the SDE data with an ATMOVE marker */
      new_jbg_buf = jbg_buf_init(&s->free_list);
      jbg_buf_write(MARKER_ESC, new_jbg_buf);
      jbg_buf_write(MARKER_ATMOVE, new_jbg_buf);
      jbg_buf_write((new_tx_line >> 24) & 0xff, new_jbg_buf);
      jbg_buf_write((new_tx_line >> 16) & 0xff, new_jbg_buf);
      jbg_buf_write((new_tx_line >> 8) & 0xff, new_jbg_buf);
      jbg_buf_write(new_tx_line & 0xff, new_jbg_buf);
      jbg_buf_write(new_tx, new_jbg_buf);
      jbg_buf_write(0, new_jbg_buf);
      jbg_buf_prefix(new_jbg_buf, &s->sde[stripe][layer][plane]);
    }
  }

#if 0
  if (stripe == s->stripes - 1)
    fprintf(stderr, "tp_lines = %ld, tp_exceptions = %ld, tp_pixels = %ld, "
	    "dp_pixels = %ld, encoded_pixels = %ld\n",
	    tp_lines, tp_exceptions, tp_pixels, dp_pixels, encoded_pixels);
#endif

  return;
}


/*
 * Create the next lower resolution version of an image
 */
static void resolution_reduction(struct jbg_enc_state *s, int plane,
				 int higher_layer)
{
  unsigned long hx, hy, lx, ly, hbpl, lbpl;
  unsigned char *hp1, *hp2, *hp3, *lp;
  unsigned long line_h1, line_h2, line_h3, line_l2;
  unsigned long i, j;
  int pix, k, l;

  /* number of pixels in highres image */
  hx = jbg_ceil_half(s->xd, s->d - higher_layer);
  hy = jbg_ceil_half(s->yd, s->d - higher_layer);
  /* number of pixels in lowres image */
  lx = jbg_ceil_half(hx, 1);
  ly = jbg_ceil_half(hy, 1);
  /* bytes per line in highres and lowres image */
  hbpl = jbg_ceil_half(hx, 3);
  lbpl = jbg_ceil_half(lx, 3);
  /* pointers to first image bytes */
  hp2 = s->lhp[s->highres[plane]][plane];
  hp1 = hp2 + hbpl;
  hp3 = hp2 - hbpl;
  lp = s->lhp[1 - s->highres[plane]][plane];
  
#ifdef DEBUG
  fprintf(stderr, "resolution_reduction: plane = %d, higher_layer = %d\n",
	  plane, higher_layer);
#endif

  /*
   * Layout of the variables line_h1, line_h2, line_h3, which contain
   * as bits the high resolution neighbour pixels of the currently coded
   * lowres pixel /\:
   *              \/
   *
   *   76543210 76543210 76543210 76543210     line_h3
   *   76543210 76543210 765432/\ 76543210     line_h2
   *   76543210 76543210 765432\/ 76543210     line_h1
   *
   * Layout of the variable line_l2, which contains the low resolution
   * pixels near the currently coded pixel as bits. The lowres pixel
   * which is currently coded is marked as X:
   *
   *   76543210 76543210 76543210 76543210     line_l2
   *                            X
   */

  for (i = 0; i < ly; i++) {
    if (2*i + 1 >= hy)
      hp1 = hp2;
    pix = 0;
    line_h1 = line_h2 = line_h3 = line_l2 = 0;
    for (j = 0; j < lbpl * 8; j += 8) {
      *lp = 0;
      line_l2 |= i ? *(lp-lbpl) : 0;
      for (k = 0; k < 8 && j + k < lx; k += 4) {
	if (((j + k) >> 2) < hbpl) {
	  line_h3 |= i ? *hp3 : 0;
	  ++hp3;
	  line_h2 |= *(hp2++);
	  line_h1 |= *(hp1++);
	}
	for (l = 0; l < 4 && j + k + l < lx; l++) {
	  line_h3 <<= 2;
	  line_h2 <<= 2;
	  line_h1 <<= 2;
	  line_l2 <<= 1;
	  pix = s->res_tab[((line_h1 >> 8) & 0x007) |
			   ((line_h2 >> 5) & 0x038) |
			   ((line_h3 >> 2) & 0x1c0) |
			   (pix << 9) | ((line_l2 << 2) & 0xc00)];
	  *lp = (*lp << 1) | pix;
	}
      }
      ++lp;
    }
    *(lp - 1) <<= lbpl * 8 - lx;
    hp1 += hbpl;
    hp2 += hbpl;
    hp3 += hbpl;
  }

#ifdef DEBUG
  {
    FILE *f;
    char fn[50];
    
    sprintf(fn, "dbg_d=%02d.pbm", higher_layer - 1);
    f = fopen(fn, "wb");
    fprintf(f, "P4\n%lu %lu\n", lx, ly);
    fwrite(s->lhp[1 - s->highres[plane]][plane], 1, lbpl * ly, f);
    fclose(f);
  }
#endif

  return;
}


/* 
 * This function is called inside the three loops of jbg_enc_out() in
 * order to write the next SDE. It has first to generate the required
 * SDE and all SDEs which have to be encoded before this SDE can be
 * created. The problem here is that if we want to output a lower
 * resolution layer, we have to allpy the resolution reduction
 * algorithm in order to get it. As we try to safe as much memory as
 * possible, the resolution reduction will overwrite previous higher
 * resolution bitmaps. Consequently, we have to encode and buffer SDEs
 * which depend on higher resolution layers before we can start the
 * resolution reduction. All this logic about which SDE has to be
 * encoded before resolution reduction is allowed is handled here.
 * This approach might be a little bit more complex than alternative
 * ways to do it, but it allows us to do the encoding with the minimal
 * possible amount of temporary memory.
 */
static void output_sde(struct jbg_enc_state *s,
		       unsigned long stripe, int layer, int plane)
{
  int lfcl;     /* lowest fully coded layer */
  long i;
  unsigned long u;
  
  assert(s->sde[stripe][layer][plane] != SDE_DONE);

  if (s->sde[stripe][layer][plane] != SDE_TODO) {
#ifdef DEBUG
    fprintf(stderr, "writing SDE: s/d/p = %2lu/%2d/%2d\n",
	    stripe, layer, plane);
#endif
    jbg_buf_output(&s->sde[stripe][layer][plane], s->data_out, s->file);
    s->sde[stripe][layer][plane] = SDE_DONE;
    return;
  }

  /* Determine the smallest resolution layer in this plane for which
   * not yet all stripes have been encoded into SDEs. This layer will
   * have to be completely coded, before we can apply the next
   * resolution reduction step. */
  lfcl = 0;
  for (i = s->d; i >= 0; i--)
    if (s->sde[s->stripes - 1][i][plane] == SDE_TODO) {
      lfcl = i + 1;
      break;
    }
  if (lfcl > s->d && s->d > 0 && stripe == 0) {
    /* perform the first resolution reduction */
    resolution_reduction(s, plane, s->d);
  }
  /* In case HITOLO is not used, we have to encode and store the higher
   * resolution layers first, although we do not need them right now. */
  while (lfcl - 1 > layer) {
    for (u = 0; u < s->stripes; u++)
      encode_sde(s, u, lfcl - 1, plane);
    --lfcl;
    s->highres[plane] ^= 1;
    if (lfcl > 1)
      resolution_reduction(s, plane, lfcl - 1);
  }
  
  encode_sde(s, stripe, layer, plane);

#ifdef DEBUG
  fprintf(stderr, "writing SDE: s/d/p = %2lu/%2d/%2d\n", stripe, layer, plane);
#endif
  jbg_buf_output(&s->sde[stripe][layer][plane], s->data_out, s->file);
  s->sde[stripe][layer][plane] = SDE_DONE;
  
  if (stripe == s->stripes - 1 && layer > 0 &&
      s->sde[0][layer-1][plane] == SDE_TODO) {
    s->highres[plane] ^= 1;
    if (layer > 1)
      resolution_reduction(s, plane, layer - 1);
  }
  
  return;
}


/*
 * Convert the table which controls the deterministic prediction
 * process from the internal format into the representation required
 * for the 1728 byte long DPTABLE element of a BIH.
 *
 * The bit order of the DPTABLE format (see also ITU-T T.82 figure 13) is
 *
 * high res:   4  5  6     low res:  0  1
 *             7  8  9               2  3
 *            10 11 12
 *
 * were 4 table entries are packed into one byte, while we here use
 * internally an unpacked 6912 byte long table indexed by the following
 * bit order:
 *
 * high res:   7  6  5     high res:   8  7  6     low res:  1  0
 * (phase 0)   4  .  .     (phase 1)   5  4  .               3  2
 *             .  .  .                 .  .  .
 *
 * high res:  10  9  8     high res:  11 10  9
 * (phase 2)   7  6  5     (phase 3)   8  7  6
 *             4  .  .                 5  4  .
 */
void jbg_int2dppriv(unsigned char *dptable, const char *internal)
{
  int i, j, k;
  int trans0[ 8] = { 1, 0, 3, 2, 7, 6, 5, 4 };
  int trans1[ 9] = { 1, 0, 3, 2, 8, 7, 6, 5, 4 };
  int trans2[11] = { 1, 0, 3, 2, 10, 9, 8, 7, 6, 5, 4 };
  int trans3[12] = { 1, 0, 3, 2, 11, 10, 9, 8, 7, 6, 5, 4 };
  
  for (i = 0; i < 1728; dptable[i++] = 0);

#define FILL_TABLE1(offset, len, trans) \
  for (i = 0; i < len; i++) { \
    k = 0; \
    for (j = 0; j < 8; j++) \
      k |= ((i >> j) & 1) << trans[j]; \
    dptable[(i + offset) >> 2] |= \
      (internal[k + offset] & 3) << ((3 - (i&3)) << 1); \
  }

  FILL_TABLE1(   0,  256, trans0);
  FILL_TABLE1( 256,  512, trans1);
  FILL_TABLE1( 768, 2048, trans2);
  FILL_TABLE1(2816, 4096, trans3);

  return;
}


/*
 * Convert the table which controls the deterministic prediction
 * process from the 1728 byte long DPTABLE format into the 6912 byte long
 * internal format.
 */
void jbg_dppriv2int(char *internal, const unsigned char *dptable)
{
  int i, j, k;
  int trans0[ 8] = { 1, 0, 3, 2, 7, 6, 5, 4 };
  int trans1[ 9] = { 1, 0, 3, 2, 8, 7, 6, 5, 4 };
  int trans2[11] = { 1, 0, 3, 2, 10, 9, 8, 7, 6, 5, …