// file name: sparse_2.C
//  written by  Ee-Chien Chang
#include <iostream.h>
#include "../include/sparse_mask.h"                   // to include   class item
#include "sparse_2.h"

//-------------------------
//   Revised on 10 April 1997   to support non-square,  sparse matrix
//                              with non-power-of-2    width and height
//-------------------------
//                                                         E.C.Chang
//----------------------------------------------------------------
//
//    This program defines  3  classes
//
//    1)  pixel,
//    2)  sparse_2,  and
//    3)  coeff_hier.
//
//    Only coeff_hier should be used by other program.
//    (pixel and sparse_2 are 'private')
//
//    A better way is to define sparse_2 as a sub-class of coeff_hier.
//
//-----------------------------------------------
// USAGE:
//  
//    C O E F F   H I E R
//
//     // create 5 level of sparse_2 of 3-byte per entry.
//     // Each sparse_2 is divided in pieces of size 64x64 pixel.
//     // 
//     // level 0  2002 x 525    
//     // level 1  1001 x 263
//     // level 2   501 x 132
//     // level 3   251 x  66
//     // level 4   126 x  33
//
//     coeff_hier C ( 4, 2002,525, 64,3 );
//
//
//     // set  pixel (i,j) in level 0 to white colour.
//     // please refer to S P A R S E _ 2 for more info.
//
//     C[0].s(i,j).r = 255;
//     C[0].s(i,j).g = 255;
//     C[0].s(i,j).b = 255;
//----------------------------------------------------------
//         S P A R S E   _  2
//------------------------------------------------
//
//  Example of usage:
//       // create a 1001 by 263 sparse matrix.
//       // 3 is the number of byte per pixel.
//       // currently, only 3 byte per pixel is supported.
//       //
//       // The sparse matrix is represented by  64 by 32
//       // sub matrix.
//       // the parameters 1 and 0 specify the "padding" of 
//       // zeros to be added.
//
//  sparse_2   MAT ( 1001, 263,  1,0, 64,32,3 )
//
//  MAT.s(5,6).r = 255;
//  MAT.s(5,6).g = 255;
//  MAT.s(5,6).b = 255;
//
//  cout << MAT.s(5,6) << endl;
//  // give the output     ( 255, 255, 255 )
//
//  //  however, if we are sure that at (5,6), there is a value.
//  //  That is, it was assigned a value before.  Then
//
//  cout << MAT.fs (5,6) << endl;
//  // perform the same but run faster.  (fs stands for fast s )
//
//  MAT.fs (6,6);   // give segmentation error if (6,6) is not
//                  // yet assigned a value.
//
//-------------------------------------------------------
//  NOTE:
//
//     1) the difference between  MAT.fs and MAT.s  is
//          s :  check whether memory has been allocated for (i,j)
//               if not, it will allocate it.
//         fs :  does not perform the above checking.
//
//----------------------------------------------------
 
//--------------------------------------------
//     P I X E L 
//---------------------------------------------
ostream &
operator << ( ostream& os, pixel p )
{
    os << "(" << (int) p.r << "," <<(int) p.g << "," <<(int) p.b << ")" ;
}

//--------------
//     1000   --- return -->    3
//        1   --- return -->    0
//     1010   --- return -->    4
//   100110   --- return -->    6
//        0   --- return -->    0
//--------------
int sparse_2::log( int x )
{
   int tmp_x = x;
 
   int count =0;
   int flag  =0;
 
   while (  tmp_x > 1 )
     {
          if (tmp_x%2) flag = 1;
 
          count ++;
          tmp_x >>= 1;
     }

 
   return (count+flag);
}
//---------------------------------------------------------
//      S P A R S E     -  2   
//  two level sparse matrix 
//
//   each entry pixel consists of (byte_per_pixel) number of bytes;
//
//   IMPORTANT::
//        there are two type of width here.
//        Padding width : is the width after padding is added
//                width : is the width actual widht. 
//
//       get_non_padding_width ()              
//       get_padding_width     ()  
//
//   only get_padding_width is public. 
//---------------------------------------------------------
sparse_2::~sparse_2 ()
{
     for ( int i =0; i< table_size_h * table_size_wi ; i++)
        if ( tb[i] != 0 ) delete [] tb[i];
   
     if (tb) delete [] tb;
}

//---------------------------------------
//
//   unit_w, unit_h       :  It is not necessary for unit_w and unit_h
//                           to be power of 2.
//                           However, to save space, 
//                           it is RECOMMENDED that both  be power of 2 
//
//   IMPORTANT::
//        there are two type of width here.
//        Padding width : is the width after padding is added
//                width : is the width actual widht. 
//
//---------------------------------------

sparse_2::sparse_2( int Width    , int Height,
                    int padding_w, int padding_h,
                    int u_size_w,int u_size_h,
                    int b_per_pixel )
{

   if (b_per_pixel != 3 )
        {
           cout << "sparse_2::ERROR... not yet implemement" << endl;
           exit(0);
        }
   
   non_padding_width    = Width;
   non_padding_height   = Height;

   width                = Width +padding_w;
   height               = Height+padding_h;

   byte_per_pixel       = b_per_pixel;
   actual_unit_size_w   = u_size_w;
   actual_unit_size_h   = u_size_h;

   unit_size_w          = u_size_w;
   unit_size_h          = u_size_h; 

   log_unit_size_w      = log ( unit_size_w );
   unit_size_w          =  1 << log_unit_size_w   ;
   unit_size_mask_w     = (1 << log_unit_size_w   ) -1 ;

   log_unit_size_h      = log ( unit_size_h );
   unit_size_h          =  1 << log_unit_size_h   ;
   unit_size_mask_h     = (1 << log_unit_size_h   ) -1 ;

   table_size_wi        =  width / unit_size_w;
   table_size_h         =  height/ unit_size_h;

   if (height%unit_size_h)   table_size_h ++;
   if (width %unit_size_w)   table_size_wi++;

   log_table_size_w     =  log ( table_size_wi );
   table_size_wi        =  1 << log_table_size_w;

   tb             = new pixel * [ table_size_wi* table_size_h ];
   number_yet_rec = new int     [ table_size_wi* table_size_h ];


   for (int i =0; i< table_size_wi * table_size_h ; i++)
      {  tb [i] = 0;  number_yet_rec [i] = 0; }

   for (int i =0; i< table_size_h ; i++)
   for (int j =0; j< table_size_wi; j++)
      {
          int tw =   width - j*unit_size_w;
          int th =   height- i*unit_size_h;
 
          if (th<0) th=0;  if (th > unit_size_h) th = unit_size_h;
          if (tw<0) tw=0;  if (tw > unit_size_w) tw = unit_size_w;

          number_yet_rec[ ( i << log_table_size_w) | j  ] = tw*th;
      }

}

//-------------------
//   check the sparseness of the coeff
//-----------------------------------------------
//  Note the different betweent sparseness of pyramid and
//                              spareeness of coeff.
//
//  the size of pyramid at level 0 is double the
//      size of the  coeff.
//----

int sparse_2::sparseness ( int x1, int y1, int x2, int y2 )
{
    int total = 0;

    x1 = x1 *2;  x2 = (x2)*2;
    y1 = y1 *2;  y2 = (y2)*2;

    for (int i =   (y1 >> log_unit_size_h); i <= (y2 >> log_unit_size_h); i++ )
    for (int j =   (x1 >> log_unit_size_w); j <= (x2 >> log_unit_size_w); j++ )
      {
             total += number_yet_rec [ (i<< log_table_size_w) | j ];
      }

    if (total<0) // Something is wrong
      {
         cout << " st row " << x1 << "en row " << x2 
              << " st col " << y1 << "en col " << y2 << endl;

         cout << total << endl;
         for (int i=0 ; i< table_size_h;  i++ )
           {
              for (int j=0 ; j< table_size_wi; j++ )
                  cout << " " << number_yet_rec [ (i << log_table_size_w )|j ];
              cout << endl;
           }
         exit (0);
      }
    return (total);
}

void sparse_2::reduce_sparseness ( int row, int st_col, int en_col )
{
 
  row    *= 2;
  st_col *= 2;
  en_col  = 2* en_col +1;

  int tmp_w1=   st_col >> log_unit_size_w;
  int tmp_h =   row    >> log_unit_size_h;

  int tmp_w2=   en_col >> log_unit_size_w;  

  if (tmp_w1== tmp_w2)
    {
       number_yet_rec [ (tmp_h << log_table_size_w  ) | tmp_w1 ]
                  -=   (en_col - st_col +1 ) *2;
    }
 else
    {
       int count = (tmp_h << log_table_size_w) | tmp_w1;

       number_yet_rec[ count ++ ] -= (unit_size_w-(st_col&unit_size_mask_w))*2;

       for (int i= (tmp_w1+1); i< tmp_w2;i++)
        number_yet_rec[count++] -= unit_size_w *2;

       number_yet_rec[ count    ] -= ((en_col&unit_size_mask_w)+1 )*2;
    }

}

//----------------------------------------------------------
//       C O E F F   H I E R 
//
//  a hierarchy of sparse_2
//
//  Data are stored into the structure by two routines:
//
//       void ihaar  (     );      
// and   void fillin_average (  );
//
//----------------------------------------------------------

coeff_hier::coeff_hier()
{
      nos_level_plus_1  = 0;
      unit_size         = 0;
      hier              = 0;
}

void
coeff_hier::reInit (  int Nos_level,
                      int width,
                      int height,
                      int Unit_size,
                      int sample_per_pixel )
{

     for (int i =0; i <  nos_level_plus_1; i++) delete hier[i];
     if (hier) delete [] hier;


     nos_level_plus_1 = Nos_level + 1;
     byte_per_pixel   = sample_per_pixel;
     unit_size        = Unit_size;
     hier             = new sparse_2 * [ nos_level_plus_1 ];


     int w = width;
     int h = height;

     for (int i=0;i< nos_level_plus_1; i++ )
     {

        int padding_w =  w%2;
        int padding_h =  h%2;

        if (i!= (nos_level_plus_1 - 1 )  )
             hier[i] = new sparse_2 (w,h,
                                     padding_w, padding_h,
                                     unit_size, unit_size,byte_per_pixel);
        else
          {
             hier[i] = new sparse_2 (w,h,
                                     padding_w, padding_h,
                                     w        , h        ,byte_per_pixel);

          }

        w+= padding_w;
        h+= padding_h;

        w= w/2;
        h= h/2;
     }
   
}

coeff_hier::coeff_hier ( int Nos_level, int w, int h, int u_size, int spp )
{
   reInit (Nos_level, w,h,u_size,spp);
}

coeff_hier::~coeff_hier ()
{

      for (int i=0; i<nos_level_plus_1; i++)
        {
            delete  hier[i];
        }
      if (hier) delete [] hier;
}

//-------------------------------------------------------
//  average is always stored in the "highest" level.
//      (the one with the smallest size )
//--------------------------------------------------------
// recall the data stream format is
//
//      c_1_r  c_1_g  c_1_b    c_2_r  c_2_g  c_2_b
//
//   * N O T *   the following.
//      c_1_r  c_2_r  c_3_r ....    c_k_r
//      c_1_g  c_2_g  c_3_g ....    c_k_g
//      c_1_b  c_2_b  c_3_b ....    c_k_b
//
//--------------------------------------------------------
void
coeff_hier::fillin_average ( unsigned char * data )
{
    int curr = 0;

    for (int i=0; i< hier[nos_level_plus_1 -1 ] -> get_non_padding_height(); i++)
    for (int j=0; j< hier[nos_level_plus_1 -1 ] -> get_non_padding_width (); j++)
       {
          hier[ nos_level_plus_1 - 1] -> s(i,j).r = data[curr++];
          hier[ nos_level_plus_1 - 1] -> s(i,j).g = data[curr++];
          hier[ nos_level_plus_1 - 1] -> s(i,j).b = data[curr++];
       }
}

//-----------------------------------------------------
//  from a list of raw data, (in m,v,d,h), perform
//  inverse transform and fill it into the sparse matrix
//
//  the location of the data is specfied by item (a linked list)
//  The data from level i  is filled into level i-1.
//
//--------------------------------------------------------
//
//  E.g.
//
//      level:        3
//      row :         5
//      interval:     ( 10,29 )  - ( 40, 49 )  -  (60,69 ) 
//
//  Then,  for each i in the intervals,
//  (in this example, i==0), this following is performed:
//
//  From 
//            data[0], data[1], data[2], data[3]      correspond to (m,h,d,v)
//  and from level  4, (obtained by  3+1),
//            hier[4].s(10,29)                        correspond to c
//
//  perfrom the inverse  and store the result in
//
//            hier[3].s(10,20)  hier[3].s(10,21)
//            hier[3].s(11,20)  hier[3].s(11,21)
//
//----------------------------------------------------- 
//
// recall the data stream format is
//
//      m_1_r  m_2_r  m_3_r ....    m_gap_r
//      m_1_g  m_2_g  m_3_g ....    m_gap_g
//      m_1_b  m_2_b  m_3_b ....    m_gap_b
//
//      v_1_r  v_2_r  v_3_r ....    v_gap_r
//      v_1_g  v_2_g  v_3_g ....    v_gap_g
//      v_1_b  v_2_b  v_3_b ....    v_gap_b
//
//      h_1_r  h_2_r  h_3_r ....    h_gap_r
//      h_1_g  h_2_g  h_3_g ....    h_gap_g
//      h_1_b  h_2_b  h_3_b ....    h_gap_b
//
//      d_1_r  d_2_r  d_3_r ....    d_gap_r
//      d_1_g  d_2_g  d_3_g ....    d_gap_g
//      d_1_b  d_2_b  d_3_b ....    d_gap_b
//
//
// refer to coeff.c for details
//
//   Note the role of the parameter gap:
//
//          
//------------------------------------------------------
int
coeff_hier::ihaar
                ( item * interval, int gap,
                                   int lev, int row, unsigned char *curr_data )
{
   int    curr          = 0;

   int    curr_mr       = 0;
   int    curr_mg       =   gap;
   int    curr_mb       = 2*gap;

   int    curr_vr       = 3*gap;
   int    curr_vg       = 4*gap;
   int    curr_vb       = 5*gap;

   int    curr_hr       = 6*gap;
   int    curr_hg       = 7*gap;
   int    curr_hb       = 8*gap;

   int    curr_dr       = 9*gap;
   int    curr_dg       = 10*gap;
   int    curr_db       = 11*gap;


   int    up_log_w         = hier[lev+1] -> log_unit_size_w;
   int    up_mask_w        = hier[lev+1] -> unit_size_mask_w;

   int    my_log_w         = hier[lev] -> log_unit_size_w;
   int    my_mask_w        = hier[lev] -> unit_size_mask_w;

   for ( item *  curr_interval  =  interval; 
                 curr_interval !=  0;
                 curr_interval  = curr_interval -> next )
     {
          

           hier[lev]->reduce_sparseness ( row, curr_interval -> start_from,
                                               curr_interval -> end_at );



           int    up_j_r = curr_interval->start_from;
           pixel *up_tmp = & hier[lev+1]->fs(row, (up_j_r>>up_log_w)<<up_log_w);

           int    my_j_r       = 2*  curr_interval -> start_from;
           pixel *my_j_2row    = &hier[lev]->s(2*row,
                                             (my_j_r>>my_log_w)<<my_log_w);
           pixel *my_j_2row_p_1= &hier[lev]->s(2*row+1, 
                                             (my_j_r>>my_log_w)<<my_log_w);

           for (int j= curr_interval -> start_from; j<= curr_interval->end_at; j++)
             {
                  //----------------------------------------------------
                  //   compose from mixed bit
                  // 
                  //  Alternative: might move this block to stream.
                  //               i.e let the network takecare of all this bits
                  //               operation. 
                  //----------------------------------------------------
                  int mr = curr_data[curr_mr++  ];
                  int mg = curr_data[curr_mg++  ];
                  int mb = curr_data[curr_mb++  ];


                  int vr = ( (  curr_data[curr_vr++]      << 2       )  |
                             ( (mr                        >> 4 ) & 3 )    )-256;
                  int vg = ( (  curr_data[curr_vg++]      << 2       )  |
                             ( (mg                        >> 4 ) & 3 )    )-256;
                  int vb = ( (  curr_data[curr_vb++]      << 2       )  |
                             ( (mb                        >> 4 ) & 3 )    )-256; 


                  int hr = ( (  curr_data[curr_hr++]      << 2       )  |
                             ( (mr                        >> 2 ) & 3 )    )-256;
                  int hg = ( (  curr_data[curr_hg++]      << 2       )  |
                             ( (mg                        >> 2 ) & 3 )    )-256;
                  int hb = ( (  curr_data[curr_hb++]      << 2       )  |
                             ( (mb                        >> 2 ) & 3 )    )-256;             


                  int dr = ( (  curr_data[curr_dr++]      << 2       )  |
                             ( (mr                        >> 0 ) & 3 )    )-512;
                  int dg = ( (  curr_data[curr_dg++]      << 2       )  |
                             ( (mg                        >> 0 ) & 3 )    )-512;
                  int db = ( (  curr_data[curr_db++]      << 2       )  |
                             ( (mb                        >> 0 ) & 3 )    )-512;

                  up_j_r =  j & up_mask_w;  
                  if (!up_j_r) up_tmp = & hier[lev+1] -> fs(row,j); 
                   
                  int cr = ( up_tmp[up_j_r] . r << 2  )  | ( mr >> 6       );
                  int cg = ( up_tmp[up_j_r] . g << 2  )  | ( mg >> 6       );
                  int cb = ( up_tmp[up_j_r] . b << 2  )  | ( mb >> 6       );

 
                  //---------------------------------------------------
                  //  now, perform ihaar
                  //---------------------------------------------------
                  

                  int tmp1 = cr + hr; hr = cr - hr      ; cr = tmp1;
                  int tmp2 = vr + dr; dr = vr - dr      ; vr = tmp2; 
                      tmp1 = cr + vr; vr =(cr - vr) >> 2; cr = tmp1 >> 2;
                      tmp2 = hr + dr; dr =(hr - dr) >> 2; hr = tmp2 >> 2;
                    
                      tmp1 = cg + hg; hg = cg - hg      ; cg = tmp1;
                      tmp2 = vg + dg; dg = vg - dg      ; vg = tmp2; 
                      tmp1 = cg + vg; vg =(cg - vg) >> 2; cg = tmp1 >> 2;
                      tmp2 = hg + dg; dg =(hg - dg) >> 2; hg = tmp2 >> 2;

                      tmp1 = cb + hb; hb = cb - hb      ; cb = tmp1;
                      tmp2 = vb + db; db = vb - db      ; vb = tmp2; 
                      tmp1 = cb + vb; vb =(cb - vb) >> 2; cb = tmp1 >> 2;
                      tmp2 = hb + db; db =(hb - db) >> 2; hb = tmp2 >> 2;

                   //-----------------------------------------------
                   // now fill in data
                   //-----------------------------------------------

                   my_j_r = (2*j) & my_mask_w;            
                   if (!my_j_r) 
                           {
                             my_j_2row     = & hier[lev] -> s(2*row,   2*j);
                             my_j_2row_p_1 = & hier[lev] -> s(2*row+1, 2*j);
                           }

                   my_j_2row     [ my_j_r      ].r = (unsigned char) cr;
                   my_j_2row     [ my_j_r      ].g = (unsigned char) cg;
                   my_j_2row     [ my_j_r      ].b = (unsigned char) cb;

                   my_j_2row     [ my_j_r + 1  ].r = (unsigned char) vr;
                   my_j_2row     [ my_j_r + 1  ].g = (unsigned char) vg;
                   my_j_2row     [ my_j_r + 1  ].b = (unsigned char) vb;

                   my_j_2row_p_1 [ my_j_r      ].r = (unsigned char) hr;
                   my_j_2row_p_1 [ my_j_r      ].g = (unsigned char) hg;
                   my_j_2row_p_1 [ my_j_r      ].b = (unsigned char) hb;

                   my_j_2row_p_1 [ my_j_r + 1  ].r = (unsigned char) dr;
                   my_j_2row_p_1 [ my_j_r + 1  ].g = (unsigned char) dg;
                   my_j_2row_p_1 [ my_j_r + 1  ].b = (unsigned char) db;

             } 
     }

   return curr_mr;
}