/*---------------------------------------------------------------------------*/
// Baseline Wavelet Transform Coder Construction Kit
//
// Geoff Davis
// gdavis@cs.dartmouth.edu
// http://www.cs.dartmouth.edu/~gdavis
//
// Copyright 1996 Geoff Davis 9/11/96
//
// Permission is granted to use this software for research purposes as
// long as this notice stays attached to this software.
//
/*---------------------------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream.h>
#include <iomanip.h>
#include "global.h"
#include "quantizer.h"

#define DISPERSE 1
#define ONE_PASS_CODING 1
#define MISSINGVAL 10000000

/*---------------------------------------------------------------------------*/
Quantizer::Quantizer (ErrorMetric *err) : err (err)
{
  data = NULL;
  nData = 0;
  max = min = sum = sumSq = mean = var = 0;
  initialDist = 0;
}

/*---------------------------------------------------------------------------*/
void Quantizer::getStats ()
{
  max = -MaxReal;
  min = MaxReal;
  sum = sumSq = 0;

  for (int i = 0; i < nData; i++) {
    if (data[i] < min)
      min = data[i];
    if (data[i] > max)
      max = data[i];
    sum += data[i];
    sumSq += square(data[i]);
  }
  mean = sum / (Real)nData;
  var = sumSq / (Real)nData - square (mean);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
UniformQuant::UniformQuant (MonoLayerCoder *entropy, int paramPrecision,
			    int zeroCenter, int nSt, ErrorMetric *err) : 
  Quantizer (err), entropy (entropy), paramPrecision (paramPrecision), 
  zeroCenter (zeroCenter)
{
  nStages = nSt;
  setPackTables();  // Initialize packetization tables
}

/*---------------------------------------------------------------------------*/
void UniformQuant::setDataEncode (Real *newData, int newNData, int subband)
{
  data = newData;
  nData = newNData;

  getStats ();
  
  if (zeroCenter) {
    max = fabs(max) > fabs(min) ? fabs(max) : fabs(min);
    min = -max;
  }
  imax = realToInt (max, paramPrecision);
  qmax = intToReal (imax, paramPrecision);
  // Make sure qmax >= max and -qmax <= min
  while (qmax < max)
    qmax = intToReal (++imax, paramPrecision);
  
  if (zeroCenter) {
    imin = -imax;
    qmin = -qmax;
  } else {
    imin = realToInt (min, paramPrecision);
    qmin = intToReal (imin, paramPrecision);
    // Make sure qmin <= min
    while (qmin > min)
      qmin = intToReal (--imin, paramPrecision);
  }

  imean = realToInt (mean, paramPrecision);
  qmean = intToReal (imean, paramPrecision);

  initialDist = 0;
  if (zeroCenter)
    for (int i = 0; i < nData; i++)
      initialDist += (*err)(data[i]);
  else
    for (int i = 0; i < nData; i++)
      initialDist += (*err)(data[i] - qmean);
}

/*---------------------------------------------------------------------------*/
void UniformQuant::setDataDecode (Real *newData, int newNData, 
				  int imax, int imin, int imean, int subband)
{
  data = newData;
  nData = newNData;

  if (imin < imax) {
    qmax = intToReal (imin, paramPrecision);
    
    if (zeroCenter) {
      qmin = -qmax;
    } else {
      qmin = intToReal (imin, paramPrecision);
    }
    qmean = intToReal (imean, paramPrecision);
  }
}

/*---------------------------------------------------------------------------*/
void UniformQuant::getRateDist (int precision, Real minStepSize, 
				Real &rate, Real &dist) 
{
  if (precision > 0) {
    const int  nSteps = (1<<precision)-1;
    const Real stepSize = (qmax-qmin)/(Real)nSteps;
    const Real recipStepSize = 1.0/stepSize;
    
    if (stepSize < minStepSize) {
      rate = MaxReal;
      dist = MaxReal;
      return;
  }
    
    entropy->setNSym (nSteps);
    rate = dist = 0;
    
    for (int i = 0; i < nData; i++) {
      int symbol = (int)((data[i]-qmin)*recipStepSize);
      assert (symbol < nSteps && symbol >= 0);
      rate += entropy->cost (symbol, TRUE);
      Real reconstruct = qmin + ((Real)symbol + 0.5) * stepSize;
      dist += (*err) (data[i]-reconstruct);
    }
  } else {
    rate = dist = 0;
    if (zeroCenter) {
      for (int i = 0; i < nData; i++) {
	dist += (*err) (data[i]);
      }
    } else {
      for (int i = 0; i < nData; i++) {
	dist += (*err) (data[i] - qmean);
      }
    }
  } 
}

/*---------------------------------------------------------------------------*/
void UniformQuant::quantize (Encoder **encoder, int precision, int subband, 
                             PacketEnc **pkt, int modulus)
{
  if (precision > 0) {
    const int  nSteps = (1<<precision)-1;
    const Real stepSize = (qmax-qmin)/(Real)nSteps;
    const Real recipStepSize = 1.0/stepSize;
    
    entropy->setNSym (nSteps);
    
    for (int i = 0; i < nData; i++) {
      int symbol = (int)((data[i]-qmin)*recipStepSize);
      assert (symbol < nSteps && symbol >= 0);
      if (!ONE_PASS_CODING)
        entropy->write (encoder[packet(subband,i)-1], symbol, TRUE);
      else if (subband == 0)
        pkt[packet(subband,i)-1]->writeLLsymbol(encoder[packet(subband,i)-1], 
                                              symbol, 0, 0);
      else
        pkt[packet(subband,i)-1]->writeHFsymbol(encoder[packet(subband,i)-1], 
                                              symbol, 0, 0);
    }
  }
}

/*---------------------------------------------------------------------------*/
void UniformQuant::dequantize (Decoder **decoder, int precision, int subband,
                               int *missing, PacketDec **pkt, int modulus)
{
  int p;
  if (precision > 0) {
    const int  nSteps = (1<<precision)-1;
    const Real stepSize = (qmax-qmin)/(Real)nSteps;
    int symbol;
    
    entropy->setNSym (nSteps);
    
    for (int i = 0; i < nData; i++) {
      p = packet(subband,i)-1;
      if (missing[p] == 0)
      {
        if (!ONE_PASS_CODING)
          symbol = entropy->read (decoder[p], TRUE);
        else if (subband == 0)
          symbol = pkt[p]->readLLsymbol(decoder[p], 0, 0);
        else
          symbol = pkt[p]->readHFsymbol(decoder[p], 0, 0);
      }
      else symbol = 0; // will have to do concealment
      assert (symbol < nSteps && symbol >= 0);
      data[i] = qmin + ((Real)symbol + 0.5) * stepSize;
    }
  } else {
    for (int i = 0; i < nData; i++) {
      data[i] = qmean;
    }
  }
}

/*---------------------------------------------------------------------------*/
void UniformQuant::writeHeader (Encoder *encoder, int precision, int subband)
{ 
  encoder->writeNonneg (precision);

  if (precision > 0) {
    encoder->writeInt (imax);
    if (!zeroCenter)
      encoder->writeInt (imin);
  } else {
    if (!zeroCenter)
      encoder->writeInt (imean);
  }
}

/*---------------------------------------------------------------------------*/
void UniformQuant::readHeader  (Decoder *decoder, int &precision, int subband)
{
  precision = decoder->readNonneg ();

  if (precision > 0) {
    imax = decoder->readInt ();
    qmax = intToReal (imax, paramPrecision);
    
    if (zeroCenter) {
      qmin = -qmax;
    } else {
      imin = decoder->readInt ();
      qmin = intToReal (imin, paramPrecision);
    }
    qmean = 0;
  } else {
    if (!zeroCenter) {
      imean = decoder->readInt ();
      qmean = intToReal (imean, paramPrecision);
    } else {
      qmean = 0;
      imean = realToInt (qmean, paramPrecision);
    }
    qmax = qmin = qmean;
  }
}

/*---------------------------------------------------------------------------*/
void UniformQuant::setParams (int newParamPrecision, Real newMax, 
			      Real newMin, Real newMean)
{
  paramPrecision = newParamPrecision;
  qmax = newMax;
  qmin = newMin;
  qmean = newMean;
}
/*---------------------------------------------------------------------------*/

// ------------------------- Added ------------------------------------------
// Set packetization tables

void UniformQuant::setPackTables2()
{ }

void UniformQuant::setPackTables()
{
  // Later extend to include subband dimensions, etc.
  
  int i, j;
  int temp[16][16] = {
      {  1,  2,  3,  4,  9, 10, 11, 12,  1,  2,  3,  4,  9, 10, 11, 12 },
      {  5,  6,  7,  8, 13, 14, 15, 16,  5,  6,  7,  8, 13, 14, 15, 16 },
      {  9, 10, 11, 12,  1,  2,  3,  4,  9, 10, 11, 12,  1,  2,  3,  4 },
      { 13, 14, 15, 16,  5,  6,  7,  8, 13, 14, 15, 16,  5,  6,  7,  8 },
      {  2,  1,  4,  3, 10,  9, 12, 11,  2,  1,  4,  3, 10,  9, 12, 11 },
      {  6,  5,  8,  7, 14, 13, 16, 15,  6,  5,  8,  7, 14, 13, 16, 15 },
      { 10,  9, 12, 11,  2,  1,  4,  3, 10,  9, 12, 11,  2,  1,  4,  3 },
      { 14, 13, 16, 15,  6,  5,  8,  7, 14, 13, 16, 15,  6,  5,  8,  7 },
      {  1,  2,  3,  4,  9, 10, 11, 12,  1,  2,  3,  4,  9, 10, 11, 12 },
      {  5,  6,  7,  8, 13, 14, 15, 16,  5,  6,  7,  8, 13, 14, 15, 16 },
      {  9, 10, 11, 12,  1,  2,  3,  4,  9, 10, 11, 12,  1,  2,  3,  4 },
      { 13, 14, 15, 16,  5,  6,  7,  8, 13, 14, 15, 16,  5,  6,  7,  8 },
      {  2,  1,  4,  3, 10,  9, 12, 11,  2,  1,  4,  3, 10,  9, 12, 11 },
      {  6,  5,  8,  7, 14, 13, 16, 15,  6,  5,  8,  7, 14, 13, 16, 15 },
      { 10,  9, 12, 11,  2,  1,  4,  3, 10,  9, 12, 11,  2,  1,  4,  3 },
      { 14, 13, 16, 15,  6,  5,  8,  7, 14, 13, 16, 15,  6,  5,  8,  7 } };
  
    
  // level 5 of decomposition    
  // subband 0 (LL) initialised as above
  for (i=0; i<16; i++)
    for(j=0; j<16; j++)
      subLL[i][j] = temp[i][j]; 
    
  // initialize other subbands at level 5 (LH, HL and HH)
  for (i=0; i<16; i++)
    for(j=0; j<16; j++)
    {
      switch(subLL[i][j])
      {
        case 1:  sub5[0][i][j] = 7;
                 sub5[1][i][j] = 16;
                 sub5[2][i][j] = 11; break;
        case 2:  sub5[0][i][j] = 15;
                 sub5[1][i][j] = 8;
                 sub5[2][i][j] = 12; break;
        case 3:  sub5[0][i][j] = 5;
                 sub5[1][i][j] = 14;
                 sub5[2][i][j] = 9; break;
        case 4:  sub5[0][i][j] = 13;
                 sub5[1][i][j] = 6;
                 sub5[2][i][j] = 10; break;
        case 5:  sub5[0][i][j] = 11;
                 sub5[1][i][j] = 4;
                 sub5[2][i][j] = 15; break;
        case 6:  sub5[0][i][j] = 3;
                 sub5[1][i][j] = 12;
                 sub5[2][i][j] = 16; break;
        case 7:  sub5[0][i][j] = 1;
                 sub5[1][i][j] = 10;
                 sub5[2][i][j] = 13; break;
        case 8:  sub5[0][i][j] = 9;
                 sub5[1][i][j] = 2;
                 sub5[2][i][j] = 14; break;
        case 9:  sub5[0][i][j] = 16;
                 sub5[1][i][j] = 7;
                 sub5[2][i][j] = 3; break;
        case 10: sub5[0][i][j] = 8;
                 sub5[1][i][j] = 15;
                 sub5[2][i][j] = 4; break;
        case 11: sub5[0][i][j] = 6;
                 sub5[1][i][j] = 13;
                 sub5[2][i][j] = 1; break;
        case 12: sub5[0][i][j] = 14;
                 sub5[1][i][j] = 5;
                 sub5[2][i][j] = 2; break;
        case 13: sub5[0][i][j] = 4;
                 sub5[1][i][j] = 11;
                 sub5[2][i][j] = 3; break;
        case 14: sub5[0][i][j] = 12;
                 sub5[1][i][j] = 3;
                 sub5[2][i][j] = 4; break;
        case 15: sub5[0][i][j] = 2;
                 sub5[1][i][j] = 9;
                 sub5[2][i][j] = 5; break;
        case 16: sub5[0][i][j] = 10;
                 sub5[1][i][j] = 1;
                 sub5[2][i][j] = 6; break;
      }
    }
  
  // Now initialize higher frequency subbands
  
  for (i=0; i<16; i++)
    for (j=0; j<16; j++)
    {
      sub4[0][i][j] = sub5[1][i][j];        // from sub5[1]
      sub4[0][15+i][j] = sub5[1][i][j];
      sub4[0][i][15+j] = sub5[1][i][j];    
      sub4[0][15+i][15+j] = sub5[1][i][j];  
      sub4[1][i][j] = sub5[0][i][j];        // from sub5[0]
      sub4[1][15+i][j] = sub5[0][i][j];
      sub4[1][i][15+j] = sub5[0][i][j];    
      sub4[1][15+i][15+j] = sub5[0][i][j];
      sub4[2][i][j] = sub5[2][i][j];        // from sub5[2]
      sub4[2][15+i][j] = sub5[2][i][j];
      sub4[2][i][15+j] = sub5[2][i][j];    
      sub4[2][15+i][15+j] = sub5[2][i][j];      
    }
    
  for (i=0; i<32; i++)
    for (j=0; j<32; j++)
    {
      sub3[0][i][j] = sub4[0][i][j];        // from sub4[0]
      sub3[0][31+i][j] = sub4[0][i][j];
      sub3[0][i][31+j] = sub4[0][i][j];    
      sub3[0][31+i][31+j] = sub4[0][i][j];  
      sub3[1][i][j] = sub4[1][i][j];        // from sub4[1]
      sub3[1][31+i][j] = sub4[1][i][j];
      sub3[1][i][31+j] = sub4[1][i][j];    
      sub3[1][31+i][31+j] = sub4[1][i][j];
      sub3[2][i][j] = sub4[2][i][j];        // from sub4[2]
      sub3[2][31+i][j] = sub4[2][i][j];
      sub3[2][i][31+j] = sub4[2][i][j];    
      sub3[2][31+i][31+j] = sub4[2][i][j];      
    }
  
  for (i=0; i<64; i++)
    for (j=0; j<64; j++)
    {
      sub2[0][i][j] = sub3[0][i][j];        // from sub3[0]
      sub2[0][63+i][j] = sub3[0][i][j];
      sub2[0][i][63+j] = sub3[0][i][j];    
      sub2[0][63+i][63+j] = sub3[0][i][j];  
      sub2[1][i][j] = sub3[1][i][j];        // from sub3[1]
      sub2[1][63+i][j] = sub3[1][i][j];
      sub2[1][i][63+j] = sub3[1][i][j];    
      sub2[1][63+i][63+j] = sub3[1][i][j];
      sub2[2][i][j] = sub3[2][i][j];        // from sub3[2]
      sub2[2][63+i][j] = sub3[2][i][j];
      sub2[2][i][63+j] = sub3[2][i][j];    
      sub2[2][63+i][63+j] = sub3[2][i][j];      
    }
  
  for (i=0; i<128; i++)
    for (j=0; j<128; j++)
    {
      sub1[0][i][j] = sub2[0][i][j];        // from sub2[0]
      sub1[0][127+i][j] = sub2[0][i][j];
      sub1[0][i][127+j] = sub2[0][i][j];    
      sub1[0][127+i][127+j] = sub2[0][i][j];  
      sub1[1][i][j] = sub2[1][i][j];        // from sub2[1]
      sub1[1][127+i][j] = sub2[1][i][j];
      sub1[1][i][127+j] = sub2[1][i][j];    
      sub1[1][127+i][127+j] = sub2[1][i][j];
      sub1[2][i][j] = sub2[2][i][j];        // from sub2[2]
      sub1[2][127+i][j] = sub2[2][i][j];
      sub1[2][i][127+j] = sub2[2][i][j];    
      sub1[2][127+i][127+j] = sub2[2][i][j];      
    }
}


// Find the buffer (packet) to which the current symbol should be written

int UniformQuant::packet (int subband, int position)
{
  // Later Extend the input values to include subband dimensions
  
  int row, col;
    
  switch (subband)
  {
    case 0:                    // LL
      // determine row and column and return subLL[row][col]...
      row = (int) floor(position/16); // rows start at 0; later try subbandHsize?
      col = position - row*16;  // columns start at 0;
      return subLL[row][col];
      break;
    case 1:
      row = (int) floor(position/16); // rows start at 0; later try subbandHsize?
      col = position - row*16;  // columns start at 0;
      return sub5[0][row][col];
      break;
    case 2:
      row = (int) floor(position/16); // rows start at 0; later try subbandHsize?
      col = position - row*16;  // columns start at 0;
      return sub5[1][row][col];
      break;
    case 3:
      row = (int) floor(position/16); // rows start at 0; later try subbandHsize?
      col = position - row*16;  // columns start at 0;
      return sub5[2][row][col];
      break;
    case 4:
      row = (int) floor(position/32); // rows start at 0; later try subbandHsize?
      col = position - row*32;  // columns start at 0;
      return sub4[0][row][col];
      break;
    case 5:
      row = (int) floor(position/32); // rows start at 0; later try subbandHsize?
      col = position - row*32;  // columns start at 0;
      return sub4[1][row][col];
      break;
    case 6:
      row = (int) floor(position/32); // rows start at 0; later try subbandHsize?
      col = position - row*32;  // columns start at 0;
      return sub4[2][row][col];
      break;
    case 7:
      row = (int) floor(position/64); // rows start at 0; later try subbandHsize?
      col = position - row*64;  // columns start at 0;
      return sub3[0][row][col];
      break;
    case 8:
      row = (int) floor(position/64); // rows start at 0; later try subbandHsize?
      col = position - row*64;  // columns start at 0;
      return sub3[1][row][col];
      break;
    case 9:
      row = (int) floor(position/64); // rows start at 0; later try subbandHsize?
      col = position - row*64;  // columns start at 0;
      return sub3[2][row][col];
      break;
    case 10:
      row = (int) floor(position/128); // rows start at 0; later try subbandHsize?
      col = position - row*128;  // columns start at 0;
      return sub2[0][row][col];
      break;
    case 11:
      row = (int) floor(position/128); // rows start at 0; later try subbandHsize?
      col = position - row*128;  // columns start at 0;
      return sub2[1][row][col];
      break;
    case 12:
      row = (int) floor(position/128); // rows start at 0; later try subbandHsize?
      col = position - row*128;  // columns start at 0;
      return sub2[2][row][col];
      break;
    case 13:
      row = (int) floor(position/256); // rows start at 0; later try subbandHsize?
      col = position - row*256;  // columns start at 0;
      return sub1[0][row][col];
      break;
    case 14:
      row = (int) floor(position/256); // rows start at 0; later try subbandHsize?
      col = position - row*256;  // columns start at 0;
      return sub1[1][row][col];
      break;
    case 15:
      row = (int) floor(position/256); // rows start at 0; later try subbandHsize?
      col = position - row*256;  // columns start at 0;
      return sub1[2][row][col];
      break;
    default:
      // error
      return 0;
      break;
  }
}
// ---------------------- End of Added ---------------------------------------






/*---------------------------------------------------------------------------*/
LayerQuant::LayerQuant  (MultiLayerCoder *entropy, int paramPrecision,
			 int signedSym, int nLayers, int nSt, ErrorMetric *err) :
  Quantizer (err), entropy (entropy), paramPrecision (paramPrecision), 
  signedSym (signedSym), nLayers (nLayers)
{
  currentLayer = -1;
  layerRate = new Real [nLayers];
  layerDist = new Real [nLayers];
  context = NULL;
  residual = NULL;
    
  nStages = nSt;
  if (DISPERSE)
    setPackTables();  // Initialize dispersive packetization tables
  else 
    setPackTables2(); // Initialize non-dispersive packetization tables
}

/*---------------------------------------------------------------------------*/
LayerQuant::~LayerQuant ()
{
  delete [] layerRate;
  delete [] layerDist;
  if (context != NULL)
    delete [] context;
  if (residual != NULL)
    delete [] residual;
    
}

/*---------------------------------------------------------------------------*/
void LayerQuant::setDataEncode  (Real *newData, int newNData, int subband)
{
  data = newData;
  nData = newNData;

  getStats ();
  
  if (signedSym) {
    max = fabs(max) > fabs(min) ? fabs(max) : fabs(min);
    min = -max;
  }
  imax = realToInt (max, paramPrecision);
  qmax = intToReal (imax, paramPrecision);
  // Make sure qmax >= max and -qmax <= min
  while (qmax < max)
    qmax = intToReal (++imax, paramPrecision);
  
  if (signedSym) {
    imin = -imax;
    qmin = -qmax;
  } else {
    imin = realToInt (min, paramPrecision);
    qmin = intToReal (imin, paramPrecision);
    // Make sure qmin <= min
    while (qmin > min)
      qmin = intToReal (--imin, paramPrecision);
  }
  if (signedSym)
    threshold = qmax/2.0;
  else
    threshold = (qmax - qmin)/2.0;
  
  currentLayer = -1;
  if (context != NULL)
    delete [] context;
  if (residual != NULL)
    delete [] residual;
  context = new int [nData];
  residual = new Real [nData];

  resetLayer ();

  initialDist = 0;
  for (int i = 0; i < nData; i++)
    initialDist += (*err)(residual[i]);

  //if (ONE_PASS_CODING || subband >= SUBBHEADERS) 05/17/01
    //entropy->reset (0);  
    // already done when LayerCoder::entropy created
   
}

/*---------------------------------------------------------------------------*/
void LayerQuant::setDataDecode  (Real *newData, int newNData, int subband, 
                                 int imax, int imin, int imean)
{
  data = newData;
  nData = newNData;

  if (imin < imax) {
    qmax = intToReal (imin, paramPrecision);
    qmean = intToReal (imean, paramPrecision);
    
    if (signedSym) {
      qmin = -qmax;
      threshold = qmax/2.0;
    } else {
      qmin = intToReal (imin, paramPrecision);
      threshold = (qmax - qmin)/2.0;
    }
    
  }

  currentLayer = -1;
  if (context != NULL)
    delete [] context;
  context = new int [nData];
  if (signedSym) {
    for (int i = 0; i < nData; i++) {
      data[i] = 0;
      context[i] = 0;
    }
  } else {
    for (int i = 0; i < nData; i++) {
      data[i] = qmean;
      context[i] = 0;
    }
  }
  
  //-----
  if (newCoeff != NULL)
    delete [] newCoeff;
  newCoeff = new int [nData];
  //-----
      
  resetLayer ();
  
  //if (ONE_PASS_CODING || subband >= SUBBHEADERS) 05/17/01
    //entropy->reset ();         
    // done in dequantize(...), may leave it out
}

/*---------------------------------------------------------------------------*/
void LayerQuant::getRateDist (int precision, Real minStepSize, 
			      Real &rate, Real &dist)
{
  assert (precision <= nLayers);

  Real currentRate = 0;
  Real currentDist = initialDist;
  
  for (int i = 0; i < precision; i++) {
    // rates & distortions have been computed for layers up to currentLayer
    if (i <= currentLayer) {
      currentRate += layerRate[i];
      currentDist += layerDist[i];
    } else {
      if (threshold > minStepSize) {
	qLayer (NULL);
	//printf("finished qLayer\n");
	currentRate += layerRate[i];
	currentDist += layerDist[i];
      } else {
	layerRate[i] = MaxReal;
	layerDist[i] = -MaxReal;
	currentRate = MaxReal;
	currentDist = MaxReal;
      }
    }
  }
  rate = currentRate;
  dist = currentDist;
  
}

/*---------------------------------------------------------------------------*/
void LayerQuant::quantize (Encoder **encoder, int precision, int subband,
                           PacketEnc **pkt, int modulus)
{
  int i;
  
  //printf("LayerQuant::quantize... ");
  
  resetLayer ();
  
  //printf("layer reset... ");
  
  if (!ONE_PASS_CODING)
  {
    if (subband < CORESUBBANDS)
      {}
    else if (subband < SUBBHEADERS)  // reset histograms in least significant layers
      entropy->reset (max(0, precision-LESSLAYHEADERS));
    else
      entropy->reset (0);  // reset histograms in all layers
    
    //entropy->reinitialize(subband); // set histograms  
  }
  else  
  {
    // one pass coding
    // assume 16 packets
    /* better compression when no reset
    for (i = 0; i < 16; i++)
    {
      pkt[i]->entropyLL->reset (0);
      pkt[i]->entropyHF->reset (0);
    }
    */
  }
  
  //printf("reset completed... ");

  for (i = 0; i < precision; i++) 
  {
    //printf("doing layer %d... ", i);
    quantizeLayer (encoder, precision, subband, pkt, modulus);
    //printf("pkt[0] rate = %d; ", (int)pkt[0]->rate);
  }
  
  //printf("done with all layers in this band\n");
  
}

/*---------------------------------------------------------------------------*/
void LayerQuant::resetLayer ()
{
  if (residual != NULL) {
    if (signedSym) {
      for (int i = 0 ;i < nData; i++) {
	residual[i] = data[i];
	context[i] = 0;
      }
    } else {
      // on the first layer we remove the mean
      for (int i = 0 ;i < nData; i++) {
	residual[i] = data[i] - 0.5*(qmax+qmin);
	context[i] = 0;
      }
    }
  } else {
    if (signedSym) {
      for (int i = 0 ;i < nData; i++) {
	context[i] = 0;
      }
    } else {
      // on the first layer we remove the mean
      for (int i = 0 ;i < nData; i++) {
	context[i] = 0;
      }
    }
  }
  currentLayer = -1;

  if (signedSym)
    threshold = qmax/2.0;
  else
    threshold = (qmax - qmin)/2.0;
}

// ---------------------- from original quantizer.cc file ------------------
void LayerQuant::qLayer (Encoder *encoder)
{
  const Real halfThreshold = 0.5 * threshold;
  const Real threeHalvesThreshold = 1.5 * threshold;
  Real deltaRate = 0, deltaDist = 0;
  int symbol;

  currentLayer++;

  //  printf ("current layer = %d, threshold = %g\n", currentLayer, threshold);
  if (signedSym) {
    for (int i = 0; i < nData; i++) {
      deltaDist -= (*err)(residual[i]);  // subtract off old error
      //      Real oldResid = residual[i];
      
      //printf("coefficent %d, layer %d\n", i, currentLayer); 
      
      if (context[i] == 0) {
	if (residual[i] > threshold) {
	  symbol = 1;
	  residual[i] -= threeHalvesThreshold;
	} else if (residual[i] < -threshold) {
	  symbol = -1;
	  residual[i] += threeHalvesThreshold;
	} else {
	  symbol = 0;
	}
      } else {
	if (residual[i] > 0) {
	  symbol = 1;
	  residual[i] -= halfThreshold;
	} else {
	  symbol = 0;
	  residual[i] += halfThreshold;
	}
      }

      //      printf ("%d (%d): %g -> %g\n", symbol, context[i], oldResid,
      //	      residual[i]);
      deltaDist += (*err)(residual[i]);  // add in new error
      deltaRate += entropy->write (encoder, symbol, TRUE, currentLayer,
				   context[i]);
      context[i] = 2*context[i] + symbol;
    }

  } else {
    for (int i = 0; i < nData; i++) {
      deltaDist -= (*err)(residual[i]);  // subtract off old error
      
      //printf("coefficent %d, layer %d\n", i, currentLayer);

      if (residual[i] > 0) {
	symbol = 1;
	residual[i] -= halfThreshold;
      } else {
	symbol = 0;
	residual[i] += halfThreshold;
      }

      deltaDist += (*err)(residual[i]);  // add in new error
      deltaRate += entropy->write (encoder, symbol, TRUE, currentLayer,
				   context[i]);
           				   
      context[i] = 2*context[i] + symbol;   
    }
  }

  //printf ("layer= %d  cost = %g\n", currentLayer, deltaRate);
  threshold *= 0.5;
  
  //printf("layer %d, rate %f, dist %f\n", currentLayer, deltaRate, deltaDist);
  
  layerRate[currentLayer] = deltaRate;
  //printf("Rate is fine\n"); 
  
  // THIS IS WHERE SEGMENTATION FAULT OCCURS with modified histograms !!!!!!!!
  layerDist[currentLayer] = deltaDist;
  
  //printf("Rate-distortion updated \n");
  
}


// ---------------------- End ----------------------------------------------



/*---------------------------------------------------------------------------*/
// deltaRate = bits required to code current layer
// deltaDist = reduction in distortion from current layer

void LayerQuant::quantizeLayer (Encoder **encoder, int precision, int subband,
                                PacketEnc **pkt, int modulus)
{
  const Real halfThreshold = 0.5 * threshold;
  const Real threeHalvesThreshold = 1.5 * threshold;
  Real deltaRate = 0, deltaDist = 0;
  int symbol, p;
  
  currentLayer++;

  //printf("signedSym = %d\n", signedSym);
  
  //printf ("current layer = %d, threshold = %g\n", currentLayer, threshold);
  if (signedSym) {
    for (int i = 0; i < nData; i++) {
      deltaDist -= (*err)(residual[i]);  // subtract off old error
      //      Real oldResid = residual[i];
      if (context[i] == 0) {
	if (residual[i] > threshold) {
	  symbol = 1;
	  residual[i] -= threeHalvesThreshold;
	} else if (residual[i] < -threshold) {
	  symbol = -1;
	  residual[i] += threeHalvesThreshold;
	} else {
	  symbol = 0;
	}
      } else {
	if (residual[i] > 0) {
	  symbol = 1;
	  residual[i] -= halfThreshold;
	} else {
	  symbol = 0;
	  residual[i] += halfThreshold;
	}
      }

      //      printf ("%d (%d): %g -> %g\n", symbol, context[i], oldResid,
      //	      residual[i]);
      deltaDist += (*err)(residual[i]);  // add in new error
      //printf("Symbol about to be written\n");
      p = (packet(subband,i)-1 + modulus)%16;
      if (p < 0 | p > 15)
         printf("Error: packet (%d) out of range sub %d, pos %d\n",p,subband,i);
      
      if (ONE_PASS_CODING)
      { deltaRate += pkt[p]->writeHFsymbol (encoder[p], symbol, currentLayer, 
                                            context[i]); }
      
      else if (subband < CORESUBBANDS)
      { deltaRate += entropy->write (encoder[p], 
                                     symbol, FALSE, currentLayer, context[i]);}
      else if ((subband < SUBBHEADERS && currentLayer >= precision-LESSLAYHEADERS)
                || (subband >= SUBBHEADERS))   
      { deltaRate += entropy->write (encoder[p], 
                                     symbol, TRUE, currentLayer, context[i]);}
      else
      { deltaRate += entropy->write (encoder[p], 
                                     symbol, FALSE, currentLayer, context[i]);}
      
      //printf("At least one symbol written\n");                                   
      context[i] = 2*context[i] + symbol;
    }

  } else {
    for (int i = 0; i < nData; i++) {
      deltaDist -= (*err)(residual[i]);  // subtract off old error

      if (residual[i] > 0) {
	symbol = 1;
	residual[i] -= halfThreshold;
      } else {
	symbol = 0;
	residual[i] += halfThreshold;
      }

      deltaDist += (*err)(residual[i]);  // add in new error
      //printf("Symbol about to be written\n");
      p = (packet(subband,i)-1 + modulus)%16;
      if (p < 0 | p > 15)
         printf("Error: packet (%d) out of range sub %d, pos %d\n",p,subband,i);
         
      if (ONE_PASS_CODING)
      { deltaRate += pkt[p]->writeLLsymbol (encoder[p], symbol, currentLayer, 
                                            context[i]); }
            
      else if (subband < CORESUBBANDS)
      { deltaRate += entropy->write (encoder[p], 
                                     symbol, FALSE, currentLayer, context[i]);}
      else if ((subband < SUBBHEADERS && currentLayer >= precision-LESSLAYHEADERS)
                || (subband >= SUBBHEADERS))   
      { deltaRate += entropy->write (encoder[p], 
                                     symbol, TRUE, currentLayer, context[i]);}
      else
      { deltaRate += entropy->write (encoder[p], 
                                     symbol, FALSE, currentLayer, context[i]);}

      //printf("At least one symbol written\n");                                   
      context[i] = 2*context[i] + symbol;   
    }
  }

  //  printf ("layer= %d  cost = %g\n", currentLayer, deltaRate);
  threshold *= 0.5;

  layerRate[currentLayer] = deltaRate;
  layerDist[currentLayer] = deltaDist;
}

/*---------------------------------------------------------------------------*/
void LayerQuant::dequantize   (Decoder **decoder, int precision, int subband,
                               int *missing, PacketDec **pkt, int modulus)
{
  int i;
  
  resetLayer ();
  
  if (!ONE_PASS_CODING)
  {
    if (subband < CORESUBBANDS)
      {}
    else if (subband < SUBBHEADERS)  // reset histograms in least signif. layers
      entropy->reset (max(0, precision-LESSLAYHEADERS));
    else
      entropy->reset (0);  // reset histograms in all layers

    //entropy->reinitialize(subband); // set histograms  
  }
  else  // one pass coding
  {
    // assume 16 packets
    /*
    for (i = 0; i < 16; i++)
    {
      pkt[i]->entropyLL->reset (0);
      pkt[i]->entropyHF->reset (0);
    }
    */
  }

  // ------------------------------
  for (i = 0; i < nData; i++)
    newCoeff[i] = 0;
  // ------------------------------
  
  
  for (i = 0; i < precision; i++) {
    dequantizeLayer (decoder, precision, subband, missing, pkt, modulus);
    //printf("%d layer of %d subband dequantized\n", i, subband);
  }
  
  if (newCoeff != NULL)
    delete [] newCoeff;
}

/*---------------------------------------------------------------------------*/
void LayerQuant::dequantizeLayer (Decoder **decoder, int &precision, 
                                  int subband, int *missing, PacketDec **pkt,
                                  int modulus)
{
  int symbol, p, defaultValue = 0;

  currentLayer++;

  //  printf ("current layer = %d, threshold = %g\n", currentLayer, threshold);
  if (signedSym) {
    for (int i = 0; i < nData; i++) 
    {
      p = (packet(subband,i)-1 + modulus)%16;
      if (p < 0  ||  p > 15)
        printf("Error in subband %d, position %d - packet %d\n",subband, i, p);
      if (missing[p] == 0)
        {
          if (ONE_PASS_CODING)
          { symbol = pkt[p]->readHFsymbol (decoder[p], currentLayer, context[i]); }
                                    
          else if (subband < CORESUBBANDS)
          { symbol = entropy->read (decoder[p], FALSE, currentLayer, 
                                    context[i], FALSE);}
          else if ((subband < SUBBHEADERS && currentLayer >= precision-LESSLAYHEADERS)
                    || (subband >= SUBBHEADERS))
          { symbol = entropy->read (decoder[p], TRUE, currentLayer, 
                                    context[i], FALSE); }
          else
          { symbol = entropy->read (decoder[p], FALSE, currentLayer, 
                                    context[i], FALSE); }
                                    
          //if (subband == 13) printf("%d ", symbol);
          
          //printf("Read %d from non-missing packet %d (%d)\n", symbol, p, missing[p]);
        }
      else 
        {
          if (ONE_PASS_CODING)
          { 
            symbol = defaultValue;  // check later
            
          }
                                    
          else if (subband < CORESUBBANDS)
          { symbol = defaultValue; } // will have to do concealment
          else if ((subband < SUBBHEADERS && currentLayer >= precision-LESSLAYHEADERS)
                    || (subband >= SUBBHEADERS))
          { symbol = entropy->read (decoder[p], TRUE, currentLayer, 
                                    context[i], TRUE); }
          else
          { symbol = defaultValue; } // will have to do concealment
          
          
          //entropy.freq[currentLayer][context[i]]->IncCount(symbol);
          //symbol = entropy->read (NULL, FALSE, currentLayer, context[i]);
          //printf("Avoided reading missing packet %d (%d)\n", p, missing[p]);
        }

      //      int oldData = data[i];
      
      if (floor(data[i]) == MISSINGVAL  &&  modulus > 0  &&  missing[p] == 0
          && currentLayer == 0)
      {
        data[i] = 0;
        newCoeff[i] = 1;
        //context[i] = 0; // added 05/30/01 - not needed, done in resetLayer()
                          //                  called from dequantize()
      }
              
      if (missing[p]  &&  modulus==0)
        data[i] = MISSINGVAL;
      else if (missing[p]==0  &&  (newCoeff[i]  ||  modulus==0))
      {
        if (context[i] == 0)
	  data[i] += 1.5*threshold * symbol;
        else 
	  data[i] += (symbol - 0.5) * threshold;
      }
      
      //      printf ("%d (%d):  %g -> %g  (%g)\n", symbol, context[i],
      //	      oldData, data[i], threshold);
	  
      context[i] = 2*context[i] + symbol;
    }
  } else {
    for (int i = 0; i < nData; i++) 
    {
      p = (packet(subband,i)-1 + modulus)%16;
      if (missing[p] == 0)
        {
          if (ONE_PASS_CODING)
          { symbol = pkt[p]->readLLsymbol (decoder[p], currentLayer, context[i]); }
                                              
          else if (subband < CORESUBBANDS)
          { symbol = entropy->read (decoder[p], FALSE, currentLayer, 
                                    context[i], FALSE);}
          else if ((subband < SUBBHEADERS && currentLayer >= precision-LESSLAYHEADERS)
                    || (subband >= SUBBHEADERS))
          { symbol = entropy->read (decoder[p], TRUE, currentLayer, 
                                    context[i], FALSE); }
          else
          { symbol = entropy->read (decoder[p], FALSE, currentLayer, 
                                    context[i], FALSE);}
          //printf("Read from non-missing packet %d (%d)\n", p, missing[p]);
        }
      else 
        {
          if (ONE_PASS_CODING)
          { symbol = defaultValue; } 
          
          else if (subband < CORESUBBANDS)
          { symbol = defaultValue; } // will have to do concealment
          else if ((subband < SUBBHEADERS && currentLayer >= precision-LESSLAYHEADERS)
                    || (subband >= SUBBHEADERS))
          { symbol = entropy->read (decoder[p], TRUE, currentLayer, 
                                    context[i], TRUE); }
          else
          { symbol = defaultValue; } // will have to do concealment
          
          //entropy.freq[currentLayer][context[i]]->IncCount(symbol);
          //symbol = entropy->read (NULL, FALSE, currentLayer, context[i]);
          //printf("Avoided reading missing packet %d (%d)\n", p, missing[p]);
        }
      
      if (floor(data[i]) == MISSINGVAL  &&  modulus > 0  &&  missing[p] == 0
          && currentLayer == 0)
      {
        //data[i] = 0.0;
        data[i] = qmean;
        newCoeff[i] = 1;
        // context[i] = 0; // added 05/30/01
      }
        
      if (missing[p]  &&  modulus==0)
        data[i] = MISSINGVAL;
      else if (missing[p]==0  &&  (newCoeff[i]  ||  modulus==0))
        data[i] += (symbol - 0.5) * threshold;
        
      context[i] = 2*context[i] + symbol;
    }
  }

  threshold *= 0.5;
}

/*---------------------------------------------------------------------------*/
void LayerQuant::writeHeader (Encoder *encoder, int precision, int subband)
{
  encoder->writeNonneg (precision);

  if (precision > 0) {
    encoder->writeInt (imax);
    if (!signedSym)
      encoder->writeInt (imin);
  } else {
    if (!signedSym)
      encoder->writeInt (imean);
  }

  //  printf ("precision = %d\n", precision);
  //  printf ("qmax = %g  qmin = %g  qmean = %g\n", qmax, qmin, qmean);
  
  // ---------------------- Added ------------------------------------------
  if (!ONE_PASS_CODING)
  {
    // write arith coder parameters; some of these are not necessary???
    if (subband < CORESUBBANDS)
      entropy->writeHead(encoder, precision);
    else if (subband < SUBBHEADERS)
      entropy->writeHead(encoder, precision - LESSLAYHEADERS);

/*
    if (subband == 0) entropy->writeHead(encoder, 2);
      else if (subband == 3) entropy->writeHead(encoder, 1);
*/   
  }  
  // ------------------- End of added -------------------------------------
}

/*---------------------------------------------------------------------------*/
void LayerQuant::readHeader  (Decoder *decoder, int &precision, int subband)
{
  precision = decoder->readNonneg ();

  //  printf ("precision = %d\n", precision);

  if (precision > 0) {
    imax = decoder->readInt ();
    qmax = intToReal (imax, paramPrecision);
    
    if (signedSym) {
      qmin = -qmax;
      qmean = 0;
    } else {
      imin = decoder->readInt ();
      qmin = intToReal (imin, paramPrecision);
      qmean = 0.5 * (qmax + qmin);
    }
  } else {
    if (!signedSym) {
      imean = decoder->readInt ();
      qmean = intToReal (imean, paramPrecision);
    } else {
      qmean = 0;
      imean = realToInt (qmean, paramPrecision);
    }
    qmax = qmin = qmean;
  }
  //  printf ("qmax = %g  qmin = %g  qmean = %g\n", qmax, qmin, qmean);


  // ---------------------- Added ------------------------------------------
  if (!ONE_PASS_CODING)
  {
    // read arith coder parameters; some of these are not necessary???
    if (subband < CORESUBBANDS)
      entropy->readHead(decoder, precision, subband);
/*  else if (subband < SUBBHEADERS)
      entropy->readHead(decoder, precision - LESSLAYHEADERS, subband); */

/*  
    if (subband == 0) entropy->readHead(decoder, 2, 0);
    else if (subband == 3) entropy->readHead(decoder, 1, 3);
*/  
  }
  // ------------------- End of added -------------------------------------
  
}

/*---------------------------------------------------------------------------*/
void LayerQuant::setParams (int newParamPrecision, Real newMax, 
			    Real newMin, Real newMean)
{
  paramPrecision = newParamPrecision;
  qmax = newMax;
  qmin = newMin;
  qmean = newMean;
}

/*---------------------------------------------------------------------------*/

// ------------------------- Added ------------------------------------------
// Set packetization tables

void LayerQuant::setPackTables2()
{
  // Later extend to include subband dimensions, etc.
  
  // THIS FUNCTION CAUSES SEGMENTATION FAULT AND/OR BUS ERRORS???
  
  int i, j, LLh, LLv;
  int temp[16][16] = {
      {  1,  1,  2,  2,  3,  3,  4,  4,  1,  1,  2,  2,  3,  3,  4,  4 },
      {  1,  1,  2,  2,  3,  3,  4,  4,  1,  1,  2,  2,  3,  3,  4,  4 },
      {  5,  5,  6,  6,  7,  7,  8,  8,  5,  5,  6,  6,  7,  7,  8,  8 },
      {  5,  5,  6,  6,  7,  7,  8,  8,  5,  5,  6,  6,  7,  7,  8,  8 },
      {  9,  9, 10, 10, 11, 11, 12, 12,  9,  9, 10, 10, 11, 11, 12, 12 },
      {  9,  9, 10, 10, 11, 11, 12, 12,  9,  9, 10, 10, 11, 11, 12, 12 },
      { 13, 13, 14, 14, 15, 15, 16, 16, 13, 13, 14, 14, 15, 15, 16, 16 },
      { 13, 13, 14, 14, 15, 15, 16, 16, 13, 13, 14, 14, 15, 15, 16, 16 },
      {  1,  1,  2,  2,  3,  3,  4,  4,  1,  1,  2,  2,  3,  3,  4,  4 },
      {  1,  1,  2,  2,  3,  3,  4,  4,  1,  1,  2,  2,  3,  3,  4,  4 },
      {  5,  5,  6,  6,  7,  7,  8,  8,  5,  5,  6,  6,  7,  7,  8,  8 },
      {  5,  5,  6,  6,  7,  7,  8,  8,  5,  5,  6,  6,  7,  7,  8,  8 },
      {  9,  9, 10, 10, 11, 11, 12, 12,  9,  9, 10, 10, 11, 11, 12, 12 },
      {  9,  9, 10, 10, 11, 11, 12, 12,  9,  9, 10, 10, 11, 11, 12, 12 },
      { 13, 13, 14, 14, 15, 15, 16, 16, 13, 13, 14, 14, 15, 15, 16, 16 },
      { 13, 13, 14, 14, 15, 15, 16, 16, 13, 13, 14, 14, 15, 15, 16, 16 } };
  
      
  // for 512x512 image
  switch (nStages)
  {
    case 5:  LLh = 16;
             LLv = 16;
             break; 
    case 4:  LLh = 32;
             LLv = 32;
             break;
    case 3:  LLh = 64;
             LLv = 64;
             break;
    case 2:  LLh = 128;
             LLv = 128;
             break;
    case 1:  LLh = 256;
             LLv = 256;
             break;         
    default: LLh = 16;
             LLv = 16;
             break;                
  }

  
  // initialize subLL
  for (i=0; i<LLv; i++)
    for(j=0; j<LLh; j++)
      subLL[i][j] = temp[i%16][j%16];
  
  
  // initialize other subbands at level 5 (LH, HL and HH)
  for (i=0; i<LLv; i++)
    for(j=0; j<LLh; j++)
    {
      sub5[0][i][j] = subLL[i][j];
      sub5[1][i][j] = subLL[i][j];
      sub5[2][i][j] = subLL[i][j];         
    }
  
  // Now initialize higher frequency subbands
  
  if (nStages >= 2)
  for (i=0; i<LLv; i++)
    for (j=0; j<LLh; j++)
    {
      sub4[0][2*i][2*j] = sub5[0][i][j];        // from sub5[0]
      sub4[0][2*i+1][2*j] = sub5[0][i][j];
      sub4[0][2*i][2*j+1] = sub5[0][i][j];    
      sub4[0][2*i+1][2*j+1] = sub5[0][i][j];  
      sub4[1][2*i][2*j] = sub5[1][i][j];        // from sub5[1]
      sub4[1][2*i+1][2*j] = sub5[1][i][j];
      sub4[1][2*i][2*j+1] = sub5[1][i][j];    
      sub4[1][2*i+1][2*j+1] = sub5[1][i][j];
      sub4[2][2*i][2*j] = sub5[2][i][j];        // from sub5[2]
      sub4[2][2*i+1][2*j] = sub5[2][i][j];
      sub4[2][2*i][2*j+1] = sub5[2][i][j];    
      sub4[2][2*i+1][2*j+1] = sub5[2][i][j];      
    }
    
  if (nStages >= 3)
  for (i=0; i<2*LLv; i++)
    for (j=0; j<2*LLh; j++)
    {
      sub3[0][2*i][2*j] = sub4[0][i][j];        // from sub4[0]
      sub3[0][2*i+1][2*j] = sub4[0][i][j];
      sub3[0][2*i][2*j+1] = sub4[0][i][j];    
      sub3[0][2*i+1][2*j+1] = sub4[0][i][j];  
      sub3[1][2*i][2*j] = sub4[1][i][j];        // from sub4[1]
      sub3[1][2*i+1][2*j] = sub4[1][i][j];
      sub3[1][2*i][2*j+1] = sub4[1][i][j];    
      sub3[1][2*i+1][2*j+1] = sub4[1][i][j];
      sub3[2][2*i][2*j] = sub4[2][i][j];        // from sub4[2]
      sub3[2][2*i+1][2*j] = sub4[2][i][j];
      sub3[2][2*i][2*j+1] = sub4[2][i][j];    
      sub3[2][2*i+1][2*j+1] = sub4[2][i][j];      
    }
  
  if (nStages >= 4)
  for (i=0; i<4*LLv; i++)
    for (j=0; j<4*LLh; j++)
    {
      sub2[0][2*i][2*j] = sub3[0][i][j];        // from sub3[0]
      sub2[0][2*i+1][2*j] = sub3[0][i][j];
      sub2[0][2*i][2*j+1] = sub3[0][i][j];    
      sub2[0][2*i+1][2*j+1] = sub3[0][i][j];  
      sub2[1][2*i][2*j] = sub3[1][i][j];        // from sub3[1]
      sub2[1][2*i+1][2*j] = sub3[1][i][j];
      sub2[1][2*i][2*j+1] = sub3[1][i][j];    
      sub2[1][2*i+1][2*j+1] = sub3[1][i][j];
      sub2[2][2*i][2*j] = sub3[2][i][j];        // from sub3[2]
      sub2[2][2*i+1][2*j] = sub3[2][i][j];
      sub2[2][2*i][2*j+1] = sub3[2][i][j];    
      sub2[2][2*i+1][2*j+1] = sub3[2][i][j];      
    }
  
  if (nStages >= 5)
  for (i=0; i<8*LLv; i++)
    for (j=0; j<8*LLh; j++)
    {
      sub1[0][2*i][2*j] = sub2[0][i][j];        // from sub2[0]
      sub1[0][2*i+1][2*j] = sub2[0][i][j];
      sub1[0][2*i][2*j+1] = sub2[0][i][j];    
      sub1[0][2*i+1][2*j+1] = sub2[0][i][j];  
      sub1[1][2*i][2*j] = sub2[1][i][j];        // from sub2[1]
      sub1[1][2*i+1][2*j] = sub2[1][i][j];
      sub1[1][2*i][2*j+1] = sub2[1][i][j];    
      sub1[1][2*i+1][2*j+1] = sub2[1][i][j];
      sub1[2][2*i][2*j] = sub2[2][i][j];        // from sub2[2]
      sub1[2][2*i+1][2*j] = sub2[2][i][j];
      sub1[2][2*i][2*j+1] = sub2[2][i][j];    
      sub1[2][2*i+1][2*j+1] = sub2[2][i][j];      
    }
}



void LayerQuant::setPackTables()
{
  // Later extend to include subband dimensions, etc.
  
  // THIS FUNCTION CAUSES SEGMENTATION FAULT AND/OR BUS ERRORS???
  
  int i, j, LLh, LLv;
  int nP = 16;             //number of packets
  
  //int temp[16][16] = {
  //    {  1,  2,  3,  4,  9, 10, 11, 12,  1,  2,  3,  4,  9, 10, 11, 12 },
  //    {  5,  6,  7,  8, 13, 14, 15, 16,  5,  6,  7,  8, 13, 14, 15, 16 },
  //    {  9, 10, 11, 12,  1,  2,  3,  4,  9, 10, 11, 12,  1,  2,  3,  4 },
  //    { 13, 14, 15, 16,  5,  6,  7,  8, 13, 14, 15, 16,  5,  6,  7,  8 },
  //    {  2,  1,  4,  3, 10,  9, 12, 11,  2,  1,  4,  3, 10,  9, 12, 11 },
  //    {  6,  5,  8,  7, 14, 13, 16, 15,  6,  5,  8,  7, 14, 13, 16, 15 },
  //    { 10,  9, 12, 11,  2,  1,  4,  3, 10,  9, 12, 11,  2,  1,  4,  3 },
  //    { 14, 13, 16, 15,  6,  5,  8,  7, 14, 13, 16, 15,  6,  5,  8,  7 },
  //    {  1,  2,  3,  4,  9, 10, 11, 12,  1,  2,  3,  4,  9, 10, 11, 12 },
  //    {  5,  6,  7,  8, 13, 14, 15, 16,  5,  6,  7,  8, 13, 14, 15, 16 },
  //    {  9, 10, 11, 12,  1,  2,  3,  4,  9, 10, 11, 12,  1,  2,  3,  4 },
  //    { 13, 14, 15, 16,  5,  6,  7,  8, 13, 14, 15, 16,  5,  6,  7,  8 },
  //    {  2,  1,  4,  3, 10,  9, 12, 11,  2,  1,  4,  3, 10,  9, 12, 11 },
  //    {  6,  5,  8,  7, 14, 13, 16, 15,  6,  5,  8,  7, 14, 13, 16, 15 },
  //    { 10,  9, 12, 11,  2,  1,  4,  3, 10,  9, 12, 11,  2,  1,  4,  3 },
  //    { 14, 13, 16, 15,  6,  5,  8,  7, 14, 13, 16, 15,  6,  5,  8,  7 } };
  
  int temp[4][4] = {
      {  1,  2,  3,  4 }, 
      {  5,  6,  7,  8 }, 
      {  9, 10, 11, 12 }, 
      { 13, 14, 15, 16 } };
  
  // for 512x512 image
  switch (nStages)
  {
    case 5:  LLh = 16;
             LLv = 16;
             break; 
    case 4:  LLh = 32;
             LLv = 32;
             break;
    case 3:  LLh = 64;
             LLv = 64;
             break;
    case 2:  LLh = 128;
             LLv = 128;
             break;
    case 1:  LLh = 256;
             LLv = 256;
             break;         
    default: LLh = 16;
             LLv = 16;
             break;                
  }

/*
  // for 256x256 image
  switch (nStages)
  {
    case 5:  LLh = 8;
             LLv = 8;
             break; 
    case 4:  LLh = 16;
             LLv = 16;
             break;
    case 3:  LLh = 32;
             LLv = 32;
             break;
    case 2:  LLh = 64;
             LLv = 64;
             break;
    case 1:  LLh = 128;
             LLv = 128;
             break;         
    default: LLh = 128;
             LLv = 128;
             break;                
  }

*/

  
  // initialize subLL
  for (i=0; i<LLv; i++)
    for(j=0; j<LLh; j++)
      //subLL[i][j] = temp[i%16][j%16];
      subLL[i][j] = temp[i%4][j%4];
  
  
  // initialize other subbands at lowest level (LH, HL and HH)
  for (i=0; i<LLv; i++)
    for(j=0; j<LLh; j++)
    {
      switch(subLL[i][j])
      // for full dispersion at the lowest level
      {
        case 1:  sub5[0][i][j] = 7;
                 sub5[1][i][j] = 16;
                 sub5[2][i][j] = 11; break;
        case 2:  sub5[0][i][j] = 15;
                 sub5[1][i][j] = 8;
                 sub5[2][i][j] = 12; break;
        case 3:  sub5[0][i][j] = 5;
                 sub5[1][i][j] = 14;
                 sub5[2][i][j] = 9; break;
        case 4:  sub5[0][i][j] = 13;
                 sub5[1][i][j] = 6;
                 sub5[2][i][j] = 10; break;
        case 5:  sub5[0][i][j] = 11;
                 sub5[1][i][j] = 4;
                 sub5[2][i][j] = 15; break;
        case 6:  sub5[0][i][j] = 3;
                 sub5[1][i][j] = 12;
                 sub5[2][i][j] = 16; break;
        case 7:  sub5[0][i][j] = 1;
                 sub5[1][i][j] = 10;
                 sub5[2][i][j] = 13; break;
        case 8:  sub5[0][i][j] = 9;
                 sub5[1][i][j] = 2;
                 sub5[2][i][j] = 14; break;
        case 9:  sub5[0][i][j] = 16;
                 sub5[1][i][j] = 7;
                 sub5[2][i][j] = 3; break;
        case 10: sub5[0][i][j] = 8;
                 sub5[1][i][j] = 15;
                 sub5[2][i][j] = 4; break;
        case 11: sub5[0][i][j] = 6;
                 sub5[1][i][j] = 13;
                 sub5[2][i][j] = 1; break;
        case 12: sub5[0][i][j] = 14;
                 sub5[1][i][j] = 5;
                 sub5[2][i][j] = 2; break;
        case 13: sub5[0][i][j] = 4;
                 sub5[1][i][j] = 11;
                 sub5[2][i][j] = 3; break;
        case 14: sub5[0][i][j] = 12;
                 sub5[1][i][j] = 3;
                 sub5[2][i][j] = 4; break;
        case 15: sub5[0][i][j] = 2;
                 sub5[1][i][j] = 9;
                 sub5[2][i][j] = 5; break;
        case 16: sub5[0][i][j] = 10;
                 sub5[1][i][j] = 1;
                 sub5[2][i][j] = 6; break;

/*
      switch(subLL[i][j])
// for partial dispersion at the lowest level: LH, HL and HH in the same pack.
// but not in the same as LL
// all higher freq data same as LH, HL and HH
      {
        case 1:  sub5[0][i][j] = 7;
                 sub5[1][i][j] = 7;
                 sub5[2][i][j] = 7; break;
        case 2:  sub5[0][i][j] = 15;
                 sub5[1][i][j] = 15;
                 sub5[2][i][j] = 15; break;
        case 3:  sub5[0][i][j] = 5;
                 sub5[1][i][j] = 5;
                 sub5[2][i][j] = 5; break;
        case 4:  sub5[0][i][j] = 13;
                 sub5[1][i][j] = 13;
                 sub5[2][i][j] = 13; break;
        case 5:  sub5[0][i][j] = 11;
                 sub5[1][i][j] = 11;
                 sub5[2][i][j] = 11; break;
        case 6:  sub5[0][i][j] = 3;
                 sub5[1][i][j] = 3;
                 sub5[2][i][j] = 3; break;
        case 7:  sub5[0][i][j] = 1;
                 sub5[1][i][j] = 1;
                 sub5[2][i][j] = 1; break;
        case 8:  sub5[0][i][j] = 9;
                 sub5[1][i][j] = 9;
                 sub5[2][i][j] = 9; break;
        case 9:  sub5[0][i][j] = 16;
                 sub5[1][i][j] = 16;
                 sub5[2][i][j] = 16; break;
        case 10: sub5[0][i][j] = 8;
                 sub5[1][i][j] = 8;
                 sub5[2][i][j] = 8; break;
        case 11: sub5[0][i][j] = 6;
                 sub5[1][i][j] = 6;
                 sub5[2][i][j] = 6; break;
        case 12: sub5[0][i][j] = 14;
                 sub5[1][i][j] = 14;
                 sub5[2][i][j] = 14; break;
        case 13: sub5[0][i][j] = 4;
                 sub5[1][i][j] = 4;
                 sub5[2][i][j] = 4; break;
        case 14: sub5[0][i][j] = 12;
                 sub5[1][i][j] = 12;
                 sub5[2][i][j] = 12; break;
        case 15: sub5[0][i][j] = 2;
                 sub5[1][i][j] = 2;
                 sub5[2][i][j] = 2; break;
        case 16: sub5[0][i][j] = 10;
                 sub5[1][i][j] = 10;
                 sub5[2][i][j] = 10; break;

*/
                
      }

    }


/* 
 // each subband in a separate packet
  for (i=0; i<LLv; i++)
    for(j=0; j<LLh; j++)
    {
       sub5[0][i][j] = 2;
       sub5[1][i][j] = 3;
       sub5[2][i][j] = 4;
    }
*/

  
  // Now initialize higher frequency subbands
  
  if (nStages >= 2)
  for (i=0; i<LLv; i++)
    for (j=0; j<LLh; j++)
    {
      sub4[0][i][j] = sub5[1][i][j];        // from sub5[1]
      sub4[0][LLv+i][j] = sub5[1][i][j];
      sub4[0][i][LLh+j] = sub5[1][i][j];    
      sub4[0][LLv+i][LLh+j] = sub5[1][i][j];  
      sub4[1][i][j] = sub5[0][i][j];        // from sub5[0]
      sub4[1][LLv+i][j] = sub5[0][i][j];
      sub4[1][i][LLh+j] = sub5[0][i][j];    
      sub4[1][LLv+i][LLh+j] = sub5[0][i][j];
      sub4[2][i][j] = sub5[2][i][j];        // from sub5[2]
      sub4[2][LLv+i][j] = sub5[2][i][j];
      sub4[2][i][LLh+j] = sub5[2][i][j];    
      sub4[2][LLv+i][LLh+j] = sub5[2][i][j];      
    }

/*
// subband per packet
  for (i=0; i<2*LLv; i++)
    for(j=0; j<2*LLh; j++)
    {
       sub4[0][i][j] = 5;
       sub4[1][i][j] = 6;
       sub4[2][i][j] = 7;
    }
*/
    
  if (nStages >= 3)
  for (i=0; i<2*LLv; i++)
    for (j=0; j<2*LLh; j++)
    {
      sub3[0][i][j] = sub4[0][i][j];        // from sub4[0]
      sub3[0][2*LLv+i][j] = sub4[0][i][j];
      sub3[0][i][2*LLh+j] = sub4[0][i][j];    
      sub3[0][2*LLv+i][2*LLh+j] = sub4[0][i][j];  
      sub3[1][i][j] = sub4[1][i][j];        // from sub4[1]
      sub3[1][2*LLv+i][j] = sub4[1][i][j];
      sub3[1][i][2*LLh+j] = sub4[1][i][j];    
      sub3[1][2*LLv+i][2*LLh+j] = sub4[1][i][j];
      sub3[2][i][j] = sub4[2][i][j];        // from sub4[2]
      sub3[2][2*LLv+i][j] = sub4[2][i][j];
      sub3[2][i][2*LLh+j] = sub4[2][i][j];    
      sub3[2][2*LLv+i][2*LLh+j] = sub4[2][i][j];      
    }

/*
  for (i=0; i<4*LLv; i++)
    for(j=0; j<4*LLh; j++)
    {
       sub3[0][i][j] = 8;
       sub3[1][i][j] = 9;
       sub3[2][i][j] = 10;
    }
*/
  
  if (nStages >= 4)
  for (i=0; i<4*LLv; i++)
    for (j=0; j<4*LLh; j++)
    {
      sub2[0][i][j] = sub3[0][i][j];        // from sub3[0]
      sub2[0][4*LLv+i][j] = sub3[0][i][j];
      sub2[0][i][4*LLh+j] = sub3[0][i][j];    
      sub2[0][4*LLv+i][4*LLh+j] = sub3[0][i][j];  
      sub2[1][i][j] = sub3[1][i][j];        // from sub3[1]
      sub2[1][4*LLv+i][j] = sub3[1][i][j];
      sub2[1][i][4*LLh+j] = sub3[1][i][j];    
      sub2[1][4*LLv+i][4*LLh+j] = sub3[1][i][j];
      sub2[2][i][j] = sub3[2][i][j];        // from sub3[2]
      sub2[2][4*LLv+i][j] = sub3[2][i][j];
      sub2[2][i][4*LLh+j] = sub3[2][i][j];    
      sub2[2][4*LLv+i][4*LLh+j] = sub3[2][i][j];      
    }

/*
  for (i=0; i<8*LLv; i++)
    for(j=0; j<8*LLh; j++)
    {
       sub2[0][i][j] = 11;
       sub2[1][i][j] = 12;
       sub2[2][i][j] = 13;
    }
*/
  
  if (nStages >= 5)
  for (i=0; i<8*LLv; i++)
    for (j=0; j<8*LLh; j++)
    {
      sub1[0][i][j] = sub2[0][i][j];        // from sub2[0]
      sub1[0][8*LLv+i][j] = sub2[0][i][j];
      sub1[0][i][8*LLh+j] = sub2[0][i][j];    
      sub1[0][8*LLv+i][8*LLh+j] = sub2[0][i][j];  
      sub1[1][i][j] = sub2[1][i][j];        // from sub2[1]
      sub1[1][8*LLv+i][j] = sub2[1][i][j];
      sub1[1][i][8*LLh+j] = sub2[1][i][j];    
      sub1[1][8*LLv+i][8*LLh+j] = sub2[1][i][j];
      sub1[2][i][j] = sub2[2][i][j];        // from sub2[2]
      sub1[2][8*LLv+i][j] = sub2[2][i][j];
      sub1[2][i][8*LLh+j] = sub2[2][i][j];    
      sub1[2][8*LLv+i][8*LLh+j] = sub2[2][i][j];      
    }

/*
  for (i=0; i<16*LLv; i++)
    for(j=0; j<16*LLh; j++)
    {
       sub1[0][i][j] = 14;
       sub1[1][i][j] = 15;
       sub1[2][i][j] = 16;
    }
*/

}


// Find the buffer (packet) to which the current symbol should be written

int LayerQuant::packet (int subband, int position)
{
  // Later Extend the input values to include subband dimensions
  
  int row, col, LLh;
  int nP = 16;        // number of packets

  
  // for 512x512 image
  switch (nStages)
  {
    case 5:  LLh = 16;
             break; 
    case 4:  LLh = 32;
             break;
    case 3:  LLh = 64;
             break;
    case 2:  LLh = 128;
             break;
    case 1:  LLh = 256;
             break;         
    default: LLh = 16;
             break;                
  }
  
/*
  // for 256x256 image
  switch (nStages)
  {
    case 5:  LLh = 8;
             break; 
    case 4:  LLh = 16;
             break;
    case 3:  LLh = 32;
             break;
    case 2:  LLh = 64;
             break;
    case 1:  LLh = 128;
             break;         
    default: LLh = 128;
             break;                
  }
*/

  //----------------------------------------------------------
  // added on 05/29/01
  row = (int) floor(position/ LLh); // rows start at 0; later try subbandHsize?
  col = position - row* LLh;  // columns start at 0;
  return (subLL[row][col] + subband - 1)%nP + 1;
  //----------------------------------------------------------
  
/* commented out on 05/29/01  
  switch (subband)
  {
    case 0:                    
      row = (int) floor(position/ LLh); 
      col = position - row* LLh; 
      return subLL[row][col];
      break;
    case 1:
      row = (int) floor(position/ LLh);
      col = position - row* LLh; 
      return sub5[0][row][col];
      break;
    case 2:
      row = (int) floor(position/ LLh); 
      col = position - row* LLh; 
      return sub5[1][row][col];
      break;
    case 3:
      row = (int) floor(position/ LLh); 
      col = position - row* LLh; 
      return sub5[2][row][col];
      break;
    case 4:
      row = (int) floor(position/(2*LLh));
      col = position - row*(2*LLh);  
      return sub4[0][row][col];
      break;
    case 5:
      row = (int) floor(position/(2*LLh)); 
      col = position - row*(2*LLh);  
      return sub4[1][row][col];
      break;
    case 6:
      row = (int) floor(position/(2*LLh));
      col = position - row*(2*LLh); 
      return sub4[2][row][col];
      break;
    case 7:
      row = (int) floor(position/(4*LLh));
      col = position - row*(4*LLh); 
      return sub3[0][row][col];
      break;
    case 8:
      row = (int) floor(position/(4*LLh));
      col = position - row*(4*LLh); 
      return sub3[1][row][col];
      break;
    case 9:
      row = (int) floor(position/(4*LLh)); 
      col = position - row*(4*LLh);  
      return sub3[2][row][col];
      break;
    case 10:
      row = (int) floor(position/(8*LLh)); 
      col = position - row*(8*LLh); 
      return sub2[0][row][col];
      break;
    case 11:
      row = (int) floor(position/(8*LLh));
      col = position - row*(8*LLh);  
      return sub2[1][row][col];
      break;
    case 12:
      row = (int) floor(position/(8*LLh));
      col = position - row*(8*LLh); 
      return sub2[2][row][col];
      break;
    case 13:
      row = (int) floor(position/(16*LLh)); 
      col = position - row*(16*LLh); 
      return sub1[0][row][col];
      break;
    case 14:
      row = (int) floor(position/(16*LLh));
      col = position - row*(16*LLh);  
      return sub1[1][row][col];
      break;
    case 15:
      row = (int) floor(position/(16*LLh)); 
      col = position - row*(16*LLh);  
      return sub1[2][row][col];
      break;
    default:
      return 0;
      break;
  }
*/

}


int LayerQuant::getCnt(int lyr, int cntxt, int cnt)
{ 
  return entropy->getCnt(lyr, cntxt, cnt); 
}

int LayerQuant::getStp(int lyr, int cntxt, int cnt)
{ 
  return entropy->getStp(lyr, cntxt, cnt); 
}

int LayerQuant::getPts(int lyr, int cntxt, int cnt)
{ 
  return entropy->getPts(lyr, cntxt, cnt); 
}

void LayerQuant::setCnt(int cntxt, int pos, int value)
{ 
  entropy->setCnt(cntxt, pos, value); 
}

void LayerQuant::setStp(int cntxt, int pos, int value)
{ 
  entropy->setStp(cntxt, pos, value); 
}

void LayerQuant::setPts(int cntxt, int pos, int value)
{ 
  entropy->setPts(cntxt, pos, value); 
}
  


// ---------------------- End of Added ---------------------------------------

/*---------------------------------------------------------------------------*/
