/*
  JPEG_Quantize : 量子化用クラス
*/
class JPEG_Quantize
{
  unsigned int quality_; // 量子化レベル ( 0=量子化を行わない ; 1(最大)〜 )

public:

  // 量子化レベル quality を指定して構築
  JPEG_Quantize( unsigned int quality )
    : quality_( quality ) {}

  // 量子化ルーチン
  void quantize( Matrix< double >* mcu, const Matrix< unsigned char >& qTable );

  // 逆量子化ルーチン
  void invQuantize( Matrix< double >* mcu, const Matrix< unsigned char >& qTable );
};

/*
  MCUデータ列 mcu を量子化・逆量子化する
*/
template< class Op >
void QuantizeMain( Matrix< double >* mcu, const Matrix< unsigned char >& qTable, unsigned int quality, Op f )
{
  if ( quality == 0 ) return;

  assert( mcu != 0 && &qTable != 0 );

  Coord< int > sz = mcu->size();
  Coord< Matrix< unsigned char >::size_type > qSz = qTable.size();

  Coord< int > c;
  for ( c.y = 0 ; c.y < sz.y ; ++( c.y ) ) {
    for ( c.x = 0 ; c.x < sz.x ; ++( c.x ) ) {
      double rate = static_cast< double >( qTable( c.x % qSz.x, c.y % qSz.y ) ) / quality;
      if ( rate > 1 ) (*mcu)[c] = f( (*mcu)[c], rate );
    }
  }
}

/*
  MCUデータ列 mcu を量子化する
*/
void JPEG_Quantize::quantize( Matrix< double >* mcu, const Matrix< unsigned char >& qTable )
{ QuantizeMain( mcu, qTable, quality_, std::divides< double >() ); }

/*
  MCUデータ列 mcu を逆量子化する
*/
void JPEG_Quantize::invQuantize( Matrix< double >* mcu, const Matrix< unsigned char >& qTable )
{ QuantizeMain( mcu, qTable, quality_, std::multiplies< double >() ); }

#!/usr/bin/perl -W

my $size = 8;           # 係数行列のサイズ
my $pi = atan2(1,1) *4; # π の値

for ( my $v = 0 ; $v < $size ; $v++ ) {
  for ( my $u = 0 ; $u < $size ; $u++ ) {
    my $max = 0; # 最大値
    my $min = 0; # 最小値
    for( my $y = 0 ; $y < $size ; $y++ ) {
      for( my $x = 0 ; $x < $size ; $x++ ) {
        my $value = 255.0 * cos( ( 2 * $y + 1 ) * $u / ( $size * 2 ) * $pi ) * cos( ( 2 * $x + 1 ) * $v / ( $size * 2 ) * $pi );
        $max += $value if $value > 0;
        $min += $value if $value < 0;
      }
    }
    $max /= sqrt( 2 ) if $u == 0;
    $max /= sqrt( 2 ) if $v == 0;
    $min /= sqrt( 2 ) if $u == 0;
    $min /= sqrt( 2 ) if $v == 0;
    printf( "( %.3f - %.3f )", $min * 2 / $size, $max * 2 / $size );
  }
  print "\n";
}

// ハフマン符号の型
typedef unsigned int HuffmanCode;

// ハフマン符号長の型
typedef unsigned char HuffmanLen;

// ハフマン符号型変数の構成
struct HuffmanCodeType
{
  HuffmanCode code; // 符号
  HuffmanLen len;   // 符号長

  // コンストラクタ
  HuffmanCodeType( HuffmanCode c, HuffmanLen l )
  : code( c ), len( l ) {}
};

const HuffmanLen SSSS_SIZE = 16; // SSSS の最大値

/*
  RoundHalfUp : 浮動小数点型のデータ d を小数点以下四捨五入して整数型として返す
*/
int RoundHalfUp( double d )
{
  if ( d > 0 )
    d += 0.5;
  else
    d -= 0.5;

  return( static_cast< int >( d ) );
}

/*
  CalcSSSS : data からSSSSを求めて返す
*/
HuffmanLen CalcSSSS( int data )
{
  // int 型のビット数
  const unsigned int UINT_BIT = CHAR_BIT * sizeof( int );

  unsigned int i = UINT_BIT / 2;       // チェックするビット位置の初期値
  unsigned int gap = UINT_BIT / 4;     // 次の位置への増減分(二分探索用)

  if ( data < 0 ) data *= -1;

  // ビット位置がゼロになるまでチェック
  for ( ; ; ) {
    unsigned int upper = data >> i; // 上位ビット( i だけシフトした結果 )
    if ( upper == 0 ) {          // シフトしすぎ
      if ( i == 0 ) return( 0 ); // data = 0 なら SSSS もゼロ
      i -= gap;
    } else if ( upper > 1 ) {    // シフトが足りない
      i += gap;
    } else {                     // 最上位ビットが見つかった
      break;
    }
    if ( gap > 1 ) gap /= 2;     // 増減分を半分に
  }

  return( i + 1 );
}

/*
  DC_Encode : DC成分の符号化

  mcu : 符号化対象のMCUデータ列(個々のMCUが縦横に連続して並べられている)
  mcuSize : MCU一つのサイズ
  huffmanCode : ハフマン符号表
  encData : 符号化したデータ列を保持する配列へのポインタ
*/
void DC_Encode( const Matrix< double >& mcu, int mcuSize,
                const std::map< HuffmanLen, HuffmanCodeType >& huffmanCode,
                std::vector< HuffmanCode >* encData )
{
  int prev_dc = 0;               // 前のDC成分
  Coord< int > sz( mcu.size() ); // MCU全体のサイズ

  Coord< int > c;
  HuffmanCodeType cache( 0, 0 ); // ビット列書き込み用のキャッシュ
  for ( c.y = 0 ; c.y < sz.y ; c.y += mcuSize ) {
    for ( c.x = 0 ; c.x < sz.x ; c.x += mcuSize ) {
      int diff = RoundHalfUp( mcu[c] ) - prev_dc; // 前のDC成分との差
      HuffmanLen ssss = CalcSSSS( diff );         // diffのサイズからSSSSを取得する

      // SSSSに対応するハフマンコードの取得
      std::map< HuffmanLen, HuffmanCodeType >::const_iterator cit = huffmanCode.find( ssss );
      if ( cit == huffmanCode.end() )
        throw std::domain_error( "DC_Encode : Failed to find huffman code." );

      // 符号の書き込み
      WriteToArray( ( cit->second ).code, ( cit->second ).len, &( cache.code ), &( cache.len ), encData );

      // 差の値の符号化
      if ( diff < 0 ) diff -= 1; // 負の場合は1を引く
      if ( ssss > 0 ) {
        unsigned int buff = diff;
        buff &= ~( UINT_MAX << ssss );
        WriteToArray( buff, ssss, &( cache.code ), &( cache.len ), encData );
      }
      prev_dc = RoundHalfUp( mcu[c] );
    }
  }

  // キャッシュの残りを書き込む
  if ( cache.len > 0 )
    encData->push_back( cache.code << ( HUFFMAN_CODE_BIT - cache.len ) );
}

/*
  DC_Decode : DC成分の展開

  mcu はあらかじめ必要なサイズにリサイズされているものとする

  cit : 符号化データ列への反復子のポインタ
  encEnd : 符号化データ列の末尾の次の反復子
  huffmanCode : ハフマン符号表
  mcu : 展開したDC成分を書き込む行列へのポインタ
  mcuSize : MCU一つのサイズ
*/
void DC_Decode( std::vector< HuffmanCode >::const_iterator* cit, std::vector< HuffmanCode >::const_iterator encEnd,
                const std::map< HuffmanLen, HuffmanCodeType >& huffmanCode, Matrix< double >* mcu, int mcuSize )
{
  int prev_dc = 0;     // 前のDC成分
  HuffmanLen ssss = 0; // データ列から取得するSSSS
  Coord< int > sz( mcu->size() ); // MCU全体のサイズ

  // ハフマン符号からハフマン木を生成
  std::vector< HuffmanTreeNode< HuffmanLen > > huffmanTree;
  MakeTreeFromHuffmanCode( huffmanCode, &huffmanTree );

  Coord< int > c;
  HuffmanCodeType cache( 0, 0 ); // ビット列書き込み用のキャッシュ
  for ( c.y = 0 ; c.y < sz.y ; c.y += mcuSize ) {
    for ( c.x = 0 ; c.x < sz.x ; c.x += mcuSize ) {
      // ハフマン木からSSSSを読み込む
      if ( ! ReadFromTree( huffmanTree, 0, cit, encEnd, &ssss, &cache ) == huffmanTree.max_size() )
        throw std::domain_error( "DC_Decode : Failed to read huffman tree for DC." );

      if ( ssss > 0 ) {
        // 差分をデータ列から読み込む
        unsigned int buff = 0; // ビット列を読み込むバッファ
        if ( ! ReadFromArray( cit, encEnd, &( cache.code ), &( cache.len ), &buff, ssss ) )
          throw std::domain_error( "DC_Decode : Failed to read from array for DC." );

        // 最上位ビットがゼロの場合、負の数であることを示す
        if ( ( ( buff >> ( ssss - 1 ) ) & 1 ) == 0 ) {
          buff |= UINT_MAX << ssss;
          ++buff;
        }
        prev_dc += static_cast< int >( buff );
      }
      (*mcu)[c] = prev_dc;
    }
  }
}

const HuffmanLen RUN_SIZE = 16;  // run の最大値

const HuffmanLen EOB = 0;              // EOB : run = 0, SSSS = 0
const HuffmanLen ZRL = 15 * SSSS_SIZE; // ZRL : run = 15, SSSS = 0

/*
  ZigzagRun : ジグザグスキャン

  c : 更新対象の座標へのポインタ
  d : 更新方向(増分)へのポインタ

  戻り値 すでに末尾に達していたら false を返す
*/
bool ZigzagRun( Coord< int >* c, Coord< int >* d, int mcuSize )
{
  if ( c->x == mcuSize - 1 && c->y == mcuSize - 1 )
    return( false );

  c->x += d->x;
  c->y += d->y;
  if ( c->x < 0 || c->x >= mcuSize || c->y < 0 || c->y >= mcuSize ) {
    d->x *= -1;
    d->y *= -1;
    if ( c->x < 0 ) {
      if ( c->y >= mcuSize ) {
        c->x = 1;
        c->y = mcuSize - 1;
      } else {
        c->x = 0;
      }
    } else if ( c->x >= mcuSize ) {
      c->x = mcuSize - 1;
      c->y += 2;
    } else if ( c->y < 0 ) {
      c->y = 0;
    } else {
      c->x += 2;
      c->y = mcuSize - 1;
    }
  }

  return( true );
}

/*
  AC_MCU_Encode : AC成分のエンコード(ひとつのMCU成分を処理)

  mcu : 符号化対象のMCUデータ列(個々のMCUが縦横に連続して並べられている)
  mcuSize : MCU一つのサイズ
  mcuCoord : 対象のMCUの符号化開始位置(左上位置)
  huffmanCode : ハフマン符号表
  encData : 符号化データ列へのポインタ
  cache : データ列書き込み用のキャッシュ
*/
void AC_MCU_Encode
( const Matrix< double >& mcu, int mcuSize, const Coord< int >& mcuCoord,
  const std::map< HuffmanLen, HuffmanCodeType >& huffmanCode,
  std::vector< HuffmanCode >* encData, HuffmanCodeType* cache )
{
  // ハフマン符号表用反復子の型
  typedef std::map< HuffmanLen, HuffmanCodeType >::const_iterator HC_CIt;

  HuffmanLen run = 0; // ゼロが連続した回数
  HuffmanLen zrl_count = 0; // ZRLを出力する回数

  HC_CIt zrl = huffmanCode.find( ZRL ); // ZRL用符号の位置

  Coord< int > c( 1, 0 );  // MCU内の座標位置
  Coord< int > d( -1, 1 ); // 座標の増分(左下向きからスタート)

  do {
    int ac = RoundHalfUp( mcu( mcuCoord.x + c.x, mcuCoord.y + c.y ) );
    if ( ac == 0 ) {
      // ゼロ成分の時は、run を 1 増加
      if ( ++run >= RUN_SIZE ) {
        // run が RUN_SIZE を越えたら ZRL を出力
        ++zrl_count;
        run -= RUN_SIZE;
      }
    } else {
      // 最初に ZRL を出力
      for ( HuffmanLen i = 0 ; i < zrl_count ; ++i ) {
        if ( zrl == huffmanCode.end() )
          throw std::domain_error( "AC_MCU_Encode : Failed to find huffman code ZRL." );
        WriteToArray( ( zrl->second ).code, ( zrl->second ).len, &( cache->code ), &( cache->len ), encData );
      }
      HuffmanLen ssss = CalcSSSS( ac ); // ac のサイズから SSSS を求める

      // runとssssの組合せを符号化
      HC_CIt cit = huffmanCode.find( run * SSSS_SIZE + ssss );
      if ( cit == huffmanCode.end() )
        throw std::domain_error( "AC_MCU_Encode : Failed to find huffman code RUN+SSSS." );
      WriteToArray( ( cit->second ).code, ( cit->second ).len, &( cache->code ), &( cache->len ), encData );

      // AC成分を符号化して書き込む
      if ( ac < 0 ) ac -= 1; // 負の場合は1を引く
      unsigned int buff = ac;
      buff &= ~( UINT_MAX << ssss );
      WriteToArray( buff, ssss, &( cache->code ), &( cache->len ), encData );

      run = zrl_count = 0;
    }
  } while ( ZigzagRun( &c, &d, mcuSize ) );

  // ゼロ成分が連続したまま最後に達した場合、EOBを出力
  if ( run != 0 || zrl_count != 0 ) {
    HC_CIt cit = huffmanCode.find( EOB );
    if ( cit == huffmanCode.end() )
      throw std::domain_error( "AC_MCU_Encode : Failed to find huffman code EOB." );
    WriteToArray( ( cit->second ).code, ( cit->second ).len, &( cache->code ), &( cache->len ), encData );
  }
}

/*
  AC_Encode : AC成分のエンコード

  mcu : 符号化対象のMCUデータ列(個々のMCUが縦横に連続して並べられている)
  mcuSize : MCU一つのサイズ
  huffmanCode : ハフマン符号表
  encData : 符号化データ列へのポインタ
*/
void AC_Encode
( const Matrix< double >& mcu, int mcuSize,
  const std::map< HuffmanLen, HuffmanCodeType >& huffmanCode,
  std::vector< HuffmanCode >* encData )
{
  // MCU全体のサイズ
  Coord< int > sz = mcu.size();
  if ( sz.x % mcuSize != 0 || sz.y % mcuSize != 0 )
    throw std::domain_error( "AC_Encode : Matrix size must be dividable by one MCU size" );

  Coord< int > c;
  HuffmanCodeType cache( 0, 0 ); // データ列書き込み用のキャッシュ
  for ( c.y = 0 ; c.y < sz.y ; c.y += mcuSize )
    for ( c.x = 0 ; c.x < sz.x ; c.x += mcuSize )
      AC_MCU_Encode( mcu, mcuSize, c, huffmanCode, encData, &cache );

  // キャッシュの残りを書き込む
  if ( cache.len > 0 )
    encData->push_back( cache.code << ( HUFFMAN_CODE_BIT - cache.len ) );
}

/*
  AC_MCU_Decode : AC成分のデコード(ひとつのMCU成分を処理)

  mcu はあらかじめ必要なサイズにリサイズされているものとする

  cit : 符号化データ列への反復子のポインタ
  encEnd : 符号化データ列の末尾の次の反復子
  huffmanTree : ハフマン木
  mcu : 展開したAC成分を書き込む行列へのポインタ
  mcuSize : MCU一つのサイズ
  mcuCoord : 対象のMCUの符号化開始位置(左上位置)
  cache : データ列書き込み用のキャッシュ
*/
void AC_MCU_Decode
( std::vector< HuffmanCode >::const_iterator* cit, std::vector< HuffmanCode >::const_iterator encEnd,
  const std::vector< HuffmanTreeNode< HuffmanLen > >& huffmanTree,
  Matrix< double >* mcu, int mcuSize, const Coord< int >& mcuCoord, HuffmanCodeType* cache )
{
  Coord< int > c( 1, 0 );  // MCU内の座標位置
  Coord< int > d( -1, 1 ); // 座標の増分(左下向きからスタート)
  bool isEnd = false;      // 最後のAC成分まで読み終えたか

  HuffmanLen ssssRun; // SSSSとRUNの組み合わせ
  unsigned int buff;  // 読み込み用バッファ

  do {
    // ハフマン木からSSSSとrunを抽出する
    if ( ! ReadFromTree( huffmanTree, 0, cit, encEnd, &ssssRun, cache ) == huffmanTree.max_size() )
      throw std::domain_error( "AC_MCU_Decpde : Failed to read from huffman tree." );

    if ( ssssRun == EOB ) {
      // EOB なら以降は全てゼロを書き込む
      do {
        (*mcu)( mcuCoord.x + c.x, mcuCoord.y + c.y ) = 0;
      } while ( ZigzagRun( &c, &d, mcuSize ) );
      break;
    } else if ( ssssRun == ZRL ) {
      // ZRL なら RUN_SIZE 分のゼロを書き込む
      for ( HuffmanLen i = 0 ; i < RUN_SIZE ; i++ ) {
        (*mcu)( mcuCoord.x + c.x, mcuCoord.y + c.y ) = 0;
        isEnd = ZigzagRun( &c, &d, mcuSize );
      }
    } else {
      // run と SSSS に分解
      HuffmanLen run = ssssRun / SSSS_SIZE;
      HuffmanLen ssss = ssssRun % SSSS_SIZE;

      // run の数だけゼロを書き込む
      for ( HuffmanLen i = 0 ; i < run ; i++ ) {
        (*mcu)( mcuCoord.x + c.x, mcuCoord.y + c.y ) = 0;
        ZigzagRun( &c, &d, mcuSize );
      }

      // AC成分の読み込み
      if ( ! ReadFromArray( cit, encEnd, &( cache->code ), &( cache->len ), &buff, ssss ) )
        throw std::domain_error( "AC_MCU_Decpde : Failed to read SSSS from array." );

      // 最上位ビットがゼロの場合、負の数であることを示す
      if ( ( ( buff >> ( ssss - 1 ) ) & 1 ) == 0 ) {
        buff |= UINT_MAX << ssss;
        ++buff;
      }
      (*mcu)( mcuCoord.x + c.x, mcuCoord.y + c.y ) = static_cast< int >( buff );
      isEnd = ZigzagRun( &c, &d, mcuSize );
    }
  } while ( isEnd );
}

/*
  AC_Decode : AC成分のデコード

  mcu はあらかじめ必要なサイズにリサイズされているものとする

  cit : 符号化データ列への反復子のポインタ
  encEnd : 符号化データ列の末尾の次の反復子
  huffmanCode : ハフマン符号表
  mcu : 展開したAC成分を書き込む行列へのポインタ
  mcuSize : MCU一つのサイズ
*/
void AC_Decode( std::vector< HuffmanCode >::const_iterator* cit, std::vector< HuffmanCode >::const_iterator encEnd,
                const std::map< HuffmanLen, HuffmanCodeType >& huffmanCode, Matrix< double >* mcu, int mcuSize )
{
  Coord< int > sz( mcu->size() ); // MCU全体のサイズ
  if ( sz.x % mcuSize != 0 || sz.y % mcuSize != 0 )
    throw std::domain_error( "AC_Decode : Matrix size must be dividable by one MCU size" );

  // ハフマン符号表からハフマン木を作成
  std::vector< HuffmanTreeNode< HuffmanLen > > huffmanTree;
  MakeTreeFromHuffmanCode( huffmanCode, &huffmanTree );

  Coord< int > c;
  HuffmanCodeType cache( 0, 0 ); // データ列書き込み用のキャッシュ
  for ( c.y = 0 ; c.y < sz.y ; c.y += mcuSize )
    for ( c.x = 0 ; c.x < sz.x ; c.x += mcuSize )
      AC_MCU_Decode( cit, encEnd, huffmanTree, mcu, mcuSize, c, &cache );
}