tcq_quantizer.c

JPEG2000实现的源码

/*****************************************************************************/
/* Copyright 1999 Science Applications International Corporation (SAIC).     */
/* Copyright 1995 University of Arizona, Arizona Board of Regents.           */
/* All rights reserved                                                       */
/* File: "tcq_quantizer.c"                                                   */
/* Description: Trellis coded quantization                                   */
/* Author: Joe Triscari/Tom Flohr                                            */
/* Affiliation: SAIC                                                         */
/* Version: VM7.1                                                            */
/* Last Revised: 20 April, 2000                                              */
/*****************************************************************************/

/*****************************************************************************************************
 *
 *   Functions to perfrom perfrom Adaptive Trellis Coded Quantization of sequences of data.
 *
 *****************************************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include <float.h>

/******
 *
 * tcq: Trellis quantizer. Implements the new Universal TCQ.
 *      This code has been rewritten from scratch to make it more efficient.
 *      Much of the logic now resides in look up tables.  Comments have been expanded.
 *
 ******/

#include <stdlib.h>
#include <stdio.h>

#define PRIME 2999957					    /* used to calc random pixel position */
#define CHECK_CNT 10					    /* how many pixels before mse is scaled */

void ftcq (float *class_coef, long *quant_code_words, long n1, float ql_step,
		 int *scan_order)
{
  long i;
  long stage, ovr_flow;
  long group;
  long  pos;
  float datum=0, ql_step05;
  float t0, t1, t2, t3, err_buf[8];
  float err0;
  unsigned char *direct;
  int k;

  /* Given the direction to prev state (up/down) and the current state, determines next state */
  
  unsigned char next_state[2][8] = {{0, 2, 4, 6, 0, 2, 4, 6}, {1, 3, 5, 7, 1, 3, 5, 7}};
  
  /* Given the direction to the previous state and the current state determine which code book */
  
  unsigned char next_book[2][8] = {{0, 1, 2, 3, 2, 3, 0, 1}, {2, 3, 0, 1, 0, 1, 2, 3}};
  
  /* For positive ql's, determin the delta from index to make the correct ql for each code book. */
  
  signed char delta_p[8][4] = {{0, 0, 2, 2}, {0, 0, 1, 1}, {3, 0, 1, 1}, {2, -2, 0, 0},
				 {2, 2, 0, 0}, {1, 1, 0, 0}, {1, 1, 3, 0}, {0, 0, 2, -2}};
  
  /* same for negative quantisation levels */
  
  signed char delta_m[8][4] = {{0, 0, -2, -2}, {0, 0, -1, -1}, {0, -3, -1, -1}, {2, -2, 0, 0},
				 {-2, -2, 0, 0}, {-1, -1, 0, 0}, {-1, -1, 1, -3},{0, 0, 2, -2}};

  direct = (unsigned char *)(malloc(n1 * sizeof(unsigned char)));
  if(direct == 0){
    fprintf(stderr, "forward tcq -- error on malloc of buffer memory.\n");
    exit(2);
  }
  ql_step05 = 0.5F * ql_step;
  pos = ovr_flow = 0;
  
  /* Force Trellis start state to be zero */
  for(i = 0; i < 8; i++){				    /* set initial errors at each state to zero */
    err_buf[i] = FLT_MAX;
  }
  err_buf[0] = 0;

  /* Perform Quantization */
  {
    register long index;
    
    for(k = 0; k < n1; k++){                /* loop through all the pixels */
      stage = scan_order[k];
      {
	register long e0=0, e1=0, e2=0, e3=0;
#ifdef RANDOM
	datum = class_coef[pos];
	pos = (pos + PRIME) % n1;
#else
        datum = class_coef[stage];
#endif
	index = (long)(datum / ql_step05);		    /* now many quantization levels are we? */
	if(index >= 0){
	  switch(index){
	  case 0:					    /* handle +- first 4 ql as special case */
	    e0 = e1 = 0; e2 = 1; e3 = -1;
	    break;
	  case 1:
	    e0 = e1 = 0; e2 = 1; e3 = 2;
	    break;
	  case 2:
	    e0 = e1 = 0; e2 = 1; e3 = 2;
	    break;
	  case 3:
	    e0 = 3; e1 = 0; e2 = 1; e3 = 2;
	    break;
	  case 4:
	    e0 = 3; e1 = 4; e2 = 1; e3 = 2;
	    break;
	  case 5:  
	    e0 = 3; e1 = 4; e2 = 1; e3 = 2;
	    break;
	  case 6:
	    e0 = 3; e1 = 4; e2 = 5; e3 = 2;
	    break;
	  case 7:  
	    e0 = 3; e1 = 4; e2 = 5; e3 = 2;
	    break;
	    default:
	    group = index % 8;				    /* find which subgroup we are in */
	    index = (index + 1) / 2;
	    switch(group){
	    case 0:
	      e0 = index - 1; e1 = index; e2 = index + 1; e3 = index + 2; /* 4  3456 */
	      break;
	    case 1:
	      e0 = index - 2; e1 = index - 1; e2 = index; e3 = index + 1; /* 5  3456 */
	      break;
	    case 2:
	      e0 = index + 2; e1 = index - 1; e2 = index; e3 = index + 1; /* 5  7456 */
	      break;
	    case 3:
	      e0 = index + 1; e1 = index - 2; e2 = index - 1; e3 = index; /* 6  7456 */
	      break;
	    case 4:
	      e0 = index + 1; e1 = index + 2; e2 = index - 1; e3 = index; /* 6  7856 */
	      break;
	    case 5:
	      e0 = index; e1 = index + 1; e2 = index - 2; e3 = index - 1; /* 7  7856 */
	      break;
	    case 6:
	      e0 = index; e1 = index + 1; e2 = index + 2; e3 = index - 1; /* 7  7896 */
	      break;
	    case 7:
	      e0 = index - 1; e1 = index; e2 = index + 1; e3 = index - 2; /* 8  7896 */
	      break;
	    }
	    break;
	  }
	}
	else{
	  switch(index){
	  case -1:
	    e0 = e1 = 0; e2 = -2; e3 = -1;
	    break;
	  case -2:
	    e0 = e1 = 0; e2 = -2; e3 = -1;
	    break;
	  case -3:
	    e0 = 0; e1 = -3; e2 = -2; e3 = -1;
	    break;
	  case -4:
	    e0 = -4; e1 = -3; e2 = -2; e3 = -1;
	    break;
	  case -5:  
	    e0 = -4; e1 = -3; e2 = -2; e3 = -1;
	    break;
	  case -6: case -7:
	    e0 = -4; e1 = -3; e2 = -2; e3 = -5;
	    break;
	    default:
	    group = index % 8;				    /* find which subgroup we are in */
	    index = (index - 1) / 2;
	    switch(group){
	      case 0:
	      e0 = index; e1 = index + 1; e2 = index - 2; e3 = index - 1; /* -4 4365 */
	      break;
	    case -1:
	      e0 = index + 1; e1 = index + 2; e2 = index - 1; e3 = index; /* -5 4365 */
	      break;
	    case -2:
	      e0 = index + 1; e1 = index - 2; e2 = index - 1; e3 = index; /* -5 4765 */
	      break;
	    case -3:
	      e0 = index + 2; e1 = index - 1; e2 = index; e3 = index + 1; /* -6 4765 */
	      break;
	    case -4:
	      e0 = index - 2; e1 = index - 1; e2 = index; e3 = index + 1; /* -6 8765 */
	      break;
	    case -5:
	      e0 = index - 1; e1 = index; e2 = index + 1; e3 = index + 2; /* -7 8765 */
	      break;
	    case -6:
	      e0 = index - 1; e1 = index; e2 = index + 1; e3 = index - 2; /* -7 8769 */
	      break;
	    case -7:
	      e0 = index; e1 = index + 1; e2 = index + 2; e3 = index - 1; /* -8 8769 */
	      break;
	    }
	    break;
	  }
	}
	{
	  register float datum1;
	  datum1 = datum / ql_step;
	  t0 = ((float) e0) - datum1;  /* error between datum and level d0 */
	  t1 = ((float) e1) - datum1;
	  t2 = ((float) e2) - datum1;
	  t3 = ((float) e3) - datum1;
	}
      }
      {
	register float *err;
	float e0, e1, e2, e3;
	float tmp0, tmp1;
	{
	  register float err00, err02, err10, err12;
	  
	  err = &err_buf[0];
	  e0 = t0 * t0;					    /* squared error */
	  e1 = t1 * t1;
	  e2 = t2 * t2;
	  e3 = t3 * t3;

	  err00 = *err + e0;
	  err02 = *err + e2;
	  err10 = *(err + 1) + e0;
	  err12 = *(err + 1) + e2;
	  
	  if(err00 <= err12){				    /* find min SE at state 0 */
	    *err = err00;
	    group = 0;
	  }
	  else{
	    *err = err12;
	    group = 1;
	  }
	  tmp0 = *(err + 4);
	  if(err02 <= err10){				    /* find min SE at state 4 tmp0 = *(err + 4) */
	    *(err + 4) = err02;
	  }
	  else{
	    *(err + 4) = err10;
	    group |= 16;
	  }
	}
	{
	  register float err40, err42, err50, err52;
	  
	  err40 = tmp0 + e0;
	  err42 = tmp0 + e2;
	  err50 = *(err + 5) + e0;
	  err52 = *(err + 5) + e2;
	  
	  tmp0 = *(err + 2);
	  if(err42 <= err50){				    /* find min SE at state 2 tmp0 = *(err + 2) */
	    *(err + 2) = err42;
	  }
	  else{
	    *(err + 2) = err50;
	    group |= 4;
	  }
	  tmp1 = *(err + 6);
	  if(err40 <= err52){				    /* find min SE at state 6 tmp1 = *(err + 6) */
	    *(err + 6) = err40;
	  }
	  else{
	    *(err + 6) = err52;
	    group |= 64;
	  }
	}
	{
	  register float err21, err23, err31, err33;
	  
	  err21 = tmp0 + e1;
	  err23 = tmp0 + e3;
	  err31 = *(err + 3) + e1;
	  err33 = *(err + 3) + e3;
	  
	  if(err21 <= err33){				    /* find min SE at state 1 */
	    *(err + 1) = err21;				    /* came from state 2 book 1 */
	  }
	  else{
	    *(err + 1) = err33;
	    group |= 2;
	  }
	  if(err23 <= err31){				    /* find min SE at state 5 */
	    *(err + 5) = err23;
	  }
	  else{
	    *(err + 5) = err31;
	    group |= 32;
	  }
	}
	{
	  register float err61, err63, err71, err73;
	  
	  err61 = tmp1 + e1;
	  err63 = tmp1 + e3;
	  err71 = *(err + 7) + e1;
	  err73 = *(err + 7) + e3;
	  
	  if(err63 <= err71){				    /* find min SE at state 3 */
	    *(err + 3) = err63;
	  }
	  else{
	    *(err + 3) = err71;
	    group |= 8;
	  }
	  if(err61 <= err73){				    /* find min SE at state 7 */
	    *(err + 7) = err61;
	  }
	  else{
	    *(err + 7) = err73;
	    group |= 128;
	  }
	}
	direct[stage] = (unsigned char) group;
	if(ovr_flow++ >= CHECK_CNT){			    /* keep sum in bound of float */
	  ovr_flow = 0;
	  tmp0 = err_buf[0];				    /* Find Minimum Error */
	  for(i = 1; i < 8; i++){
	    if(tmp0 > err_buf[i]){
	      tmp0 = err_buf[i];
	    }
	  }
	  for(i = 0; i < 8; i++){
	    err_buf[i] -= tmp0;
	  }
	}
      }
    }
  }

  /* Backtrace trellis path and output codewords */
  
  {
    long index, code_book, cw=0;
    float datum;
    long which_state;

    which_state = 0;					    /* Find Minimum Error */
    err0 = err_buf[0];
    for(i = 1; i < 8; i++){
      if(err0 > err_buf[i]){
	err0 = err_buf[i];
	which_state = i;
      }
    }

    for(k = n1 - 1; k >= 0; k--){                   /* loop through all the pixels */
      stage = scan_order[k];
      {
	register long bit_pos, up_down;
	
	bit_pos = 1 << which_state;			    /* bit mask of current state */
	up_down = !!(bit_pos & direct[stage]);		    /* Which branch did we take? */
	code_book = next_book[up_down][which_state];	    /* What code book is the ql out of? */
	which_state = next_state[up_down][which_state];	    /* What state are we in now? */
      }
#ifdef RANDOM
      pos = (pos - PRIME) % n1;
      if(pos < 0) pos += n1;
#else
      pos = stage;
#endif
      datum = class_coef[pos];				    /* pick up class_coef value */
      index = (long)(datum / ql_step05);		    /* now many quantization levels are we? */
      group = index % 8;				    /* granularity of 8 to handle decisions */
      switch(index){					    /* that occur at center +- 1/2 ql. */
	case 0:
	  switch(code_book){
	    case 0: case 1:
	      cw = 0;
	      break;
	    case 2:
	      cw = +1;
	      break;
	    case 3:
	      cw = -1;
	      break;
	    }
	  break;
	case 1: case 2:
	  switch(code_book){
	    case 0: case 1:
	      cw = 0;
	      break;
	    case 2: case 3:
	      cw = 1;
	      break;
	    }
	  break;
	case 3: 
	  switch(code_book){
	    case 0:
	      cw = 2;
	      break;
	    case 1:
	      cw = 0;
	      break;
	     case 2: case 3:
	      cw = 1;
	      break;
	    }
	  break;
	case 4: 
	  switch(code_book){
	    case 0: case 1:
	      cw = 2;
	      break;
	    case 2: case 3:
	      cw = 1;
	      break;
	    }
	  break;
	case -1: case -2:
	  switch(code_book){
	    case 0: case 1:
	      cw = 0;
	      break;
	    case 2: case 3:
	      cw = -1;
	      break;
	    }
	  break;
	case -3:
	  switch(code_book){
	    case 0:
	      cw = 0;
	      break;
	    case 1:
	      cw = -2;
	      break;
	    case 2: case 3:
	      cw = -1;
	      break;
	    }
	  break;
	case -4:
	  switch(code_book){
	    case 0: case 1:
	      cw = -2;
	      break;
	    case 2: case 3:
	      cw = -1;
	      break;
	    }
	  break;
	default:
	  if(index >= 0){
	    index = (index + 1) / 2;
	    cw = (index + delta_p[group][code_book]) / 2;
	    break;
	  }
	  else{
	    index = (index - 1) / 2;
	    cw = (index + delta_m[-group][code_book]) / 2;
	    break;
	  }
      }

      switch(cw){					    /* Compute optomized cw value for the */
      case -2:					    /* 2 ql's on either side of zero */
	switch(code_book & 1){
	case 1:
        if(abs(cw)>1)
          quant_code_words[pos] = cw;		    /* If in superset 1 negate codeword */
        else
          quant_code_words[pos] = -cw;
	break;
      default:
	quant_code_words[pos] = cw;
	break;
      }
	break;
      case -1:
	switch(code_book & 1){
	case 1:
	  if(abs(cw)>1)
	    quant_code_words[pos] = cw;		    /* If in superset 1 negate codeword */
	  else
	    quant_code_words[pos] = -cw;
	  break;
	default:
	  quant_code_words[pos] = cw;
	  break;
	}
	break;
      case 1:
	switch(code_book & 1){
	case 1:
	  if(abs(cw)>1)
	    quant_code_words[pos] = cw;		    /* If in superset 1 negate codeword */
	  else
	    quant_code_words[pos] = -cw;
	  break;
	default:
	  quant_code_words[pos] = cw;
	  break;
	}
	break;
      case 2:
	switch(code_book & 1){
	case 1:
	  if(abs(cw)>1)
	    quant_code_words[pos] = cw;		    /* If in superset 1 negate codeword */
	  else
	    quant_code_words[pos] = -cw;
	  break;
	default:
	  quant_code_words[pos] = cw;
	  break;
	}
	break;
      default:
	switch(code_book & 1){			    /* output codewords */
	case 1:
	  if(abs(cw)>1)
	    quant_code_words[pos] = cw;		    /* If in superset 1 negate codeword */
	  else
	    quant_code_words[pos] = -cw;
	  break;					    /* make distributions similar for the */
	default:					    /* arithmetic encoder. */
	  quant_code_words[pos] = cw;
	  break;
	}
      }
    }
    
    assert(which_state == 0);
    
  }
  free(direct);
  return;
}