image_bw.cpp

设计并实现了两种分层多描述视频编码器.通过对小波域的运动估计算法进行了分析和研究

// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

//                     I M A G E   C L A S S

// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
//           > > > >    C++ version 10.07 -  05/29/95   < < < <

// Amir Said - [e-mail]
// University of Campinas (UNICAMP)
// Campinas, SP 13081, Brazil

// William A. Pearlman - [e-mail]
// Rensselaer Polytechnic Institute
// Troy, NY 12180, USA

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

// Copyright (c) 1995 Amir Said & William A. Pearlman

// This program is Copyright (c) by Amir Said & William A. Pearlman.
// It may be freely redistributed in its entirety provided that this
// copyright notice is not removed. It may not be sold for profit or
// incorporated in commercial programs without the written permission
// of the copyright holders. This program is provided as is, without any
// express or implied warranty, without even the warranty of fitness
// for a particular purpose.


// - - Inclusion - - - - - - - - - - - - - - - - - - - - - - - - - - - -

#include "general.h"
#include "image_bw.h"
#include <string.h>

/*
#define RESIDUAL_IMAGE
*/

// - - Constants - - - - - - - - - - - - - - - - - - - - - - - - - - - -

static char * M_MSG = "< Image_SP >";

static char * R_MSG = "< Image_BW > cannot read from file";

static char * W_MSG = "< Image_BW > cannot write to file";

static char * L_MSG = "< Image_BW > larger than specified dimension";


// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

//  Auxiliary functions

// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

#ifdef LOSSLESS

static void SP_Transform(int m, int in[], int l[], int h[])
{
  int i, k, d1, d2, mm = m - 1;

  for (i = k = 0; i < m; i++, k += 2) {
    l[i] = (in[k] + in[k+1]) >> 1;
    h[i] =  in[k] - in[k+1]; }

  h[0] -= (d2 = l[0] - l[1]) >> 2;
  for (i = 1; i < mm; i++) {
    d1 = d2;  d2 = l[i] - l[i+1];
    h[i] -= (((d1 + d2 - h[i+1]) << 1) + d2 + 3) >> 3; }
  h[i] -= d2 >> 2;
}

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

static void SP_Recover(int m, int l[], int h[], int out[])
{
  int i, k, d1, d2, t;

  t = (h[m-1] += (d1 = l[m-2] - l[m-1]) >> 2);
  for (i = m - 2; i > 0; i--) {
    d2 = d1;  d1 = l[i-1] - l[i];
    t = (h[i] += (((d1 + d2 - t) << 1) + d2 + 3) >> 3); }
  h[0] += d1 >> 2;

  for (i = k = 0; i < m; i++, k += 2) {
    out[k] = l[i] + ((h[i] + 1) >> 1);
    out[k+1] = out[k] - h[i]; }
}

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

#else

static const int NumbTap = 4;


/*static const double T_LowPass[5] =
  { 0.852699,  0.377403, -0.110624, -0.023849, 0.037829 };

static const double T_HighPass[5] =
  { 0.788485, -0.418092, -0.040690,  0.064539, 0.0 };

static const double R_LowPass[5] =
  { 0.852699,  0.418092, -0.110624, -0.064539, 0.037829 };

static const double R_HighPass[5] =
  { 0.788485, -0.377403, -0.040690,  0.023849, 0.0 };
*/
static const double T_LowPass[5] =
  { 8.526986790094022e-01,  3.774028556126537e-01, -1.106244044184226e-01, -2.384946501937986e-02, 3.782845550699535e-02 };

static const double T_HighPass[5] =
  { 7.884856164056651e-01, -4.180922732222124e-01, -4.068941760955867e-02,  6.453888262893856e-02, 0.0 };

static const double R_LowPass[5] =
  { 8.526986790094022e-01,  4.180922732222124e-01, -1.106244044184226e-01, -6.453888262893856e-02, 3.782845550699535e-02 };

static const double R_HighPass[5] =
  { 7.884856164056651e-01, -3.774028556126537e-01, -4.068941760955867e-02,  2.384946501937986e-02, 0.0 };

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

inline double Filter_L(const double * f, double * v)
{
  return f[0] * v[0] +
    f[1] * (v[1] + v[-1]) + f[2] * (v[2] + v[-2]) +
    f[3] * (v[3] + v[-3]) + f[4] * (v[4] + v[-4]);
}

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

inline double Filter_H(const double * f, double * v)
{
  return f[0] * v[0] +
    f[1] * (v[1] + v[-1]) + f[2] * (v[2] + v[-2]) +
    f[3] * (v[3] + v[-3]);
}

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

static void Reflection(double * h, double * t)
{  
  for (int i = 1; i <= NumbTap; i++) {
    h[-i] = h[i];  t[i] = t[-i]; }
}

#endif

// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

//  Member functions of the class  < Image_BW >

// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

// - - Private functions - - - - - - - - - - - - - - - - - - - - - - - -

int Image_BW::max_levels(int n)
{
  int l1, l2;
  for (l1 = 0; !(n & 1); l1++) n >>= 1;
  for (l2 = l1 - 3; n; l2++) n >>= 1;
  return (l1 < l2 ? l1 : l2) - 1;
}


// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

void Image_BW::assign_mem(Image_Coord d, int b)
{
  int i, j;
  if ((b < 1) || (b > 2)) Error("Invalid number of < Image_BW > bytes");
  if ((levels >= 0) && (dim.x == d.x) && (dim.y == d.y)) return;
  free_mem();
  dim = d;
  if ((d.x < 32) || (d.x & 0x3) || (d.y < 32) || (d.y & 0x3))
    Error("invalid < Image_BW > dimension (assign_mem)");

  NEW_VECTOR(coeff, dim.x, Pel_Type *, M_MSG);
  for (i = 0; i < dim.x; i++) NEW_VECTOR(coeff[i], dim.y, Pel_Type, M_MSG);

  NEW_VECTOR(coeff1, dim.x, Pel_Type *, M_MSG);
  for (i = 0; i < dim.x; i++) NEW_VECTOR(coeff1[i], dim.y, Pel_Type, M_MSG); 
  for (i = 0; i < dim.x; i++) 
	  for (j = 0; j < dim.y; j++) coeff1[i][j] = 0;
}

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

void Image_BW::free_mem(void)
{
  if (levels >= 0) {
	  for (int i = dim.x - 1; i >= 0; i--) {
		  delete [] coeff[i]; delete [] coeff1[i]; }
	  delete [] coeff; delete [] coeff1; }
  bytes = dim.x = dim.y = 0;  levels = -1;
}

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
// - - Public functions  - - - - - - - - - - - - - - - - - - - - - - - -

void Image_BW::read_pic(Image_Coord d, char * file_name, int b)
{
  assign_mem(d, b);
  mean = levels = 0;  bytes = b;

  FILE * in_file = fopen(file_name, "rb");
  if (in_file == NULL) Error(R_MSG);

  int i, j, k, p, c;
  double s=0,mea=0;
  for (i = 0; i < dim.x; i++)
    for (j = 0; j < dim.y; j++) {
      for (p = k = 0; k < bytes; k++) {
        if ((c = getc(in_file)) == EOF) Error(R_MSG);
        p = (p << 8) | c; }
#ifdef RESIDUAL_IMAGE
      coeff[i][j] = p-128; 
#else
      coeff[i][j] = p;
#endif
			}
/*	for(i=0;i<144;i++)
		for(j=0;j<176;j++)
		   s+=coeff[i][j];
	mea=Round(s/double(144)/176);
  if (getc(in_file) != EOF) Error(L_MSG);*/
  fclose(in_file);
}


void Image_BW::P_read_pic(Image_Coord d, char * file_name, int b)
{
  assign_mem(d, b);
  mean = levels = 0;  bytes = b;
  
  int check;

  FILE * in_file = fopen(file_name, "rb");
  if (in_file == NULL) Error(R_MSG);
  
  check=fseek(in_file,(d.x*d.y)*3/2,SEEK_SET);
  if(check)
	  printf("can't open the %c",file_name);

  int i, j, k, p, c;
  for (i = 0; i < dim.x; i++)
    for (j = 0; j < dim.y; j++) {
      for (p = k = 0; k < bytes; k++) {
        if ((c = getc(in_file)) == EOF) Error(R_MSG);
        p = (p << 8) | c; }
#ifdef RESIDUAL_IMAGE
      coeff[i][j] = p-128; 
#else
      coeff[i][j] = p;
#endif
			}
//  if (getc(in_file) != EOF) Error(L_MSG);
  fclose(in_file);
}

//----------------------------------------------------------------------

void Image_BW::send_coeff(Image_Coord  d, Pel_Type ** Pcoeff)
{
//     assign_mem(d, 1);
     for(int i=0;i<d.x;i++)
		 for(int j=0;j<d.y;j++)
			 coeff[i][j]=Pcoeff[i][j];

}

//----------------------------------------------------------------------

void Image_BW::send_coeff1(Image_Coord  d, Pel_Type ** Pcoeff)
{
//     assign_mem(d, 1);
     for(int i=0;i<d.x;i++)
		 for(int j=0;j<d.y;j++)
			 coeff1[i][j]=Pcoeff[i][j];

}

//----------------------------------------------------------------------

void Image_BW::Desend_coeff(Image_Coord  d, Pel_Type ** Pcoeff)
{
//     assign_mem(d, 1);
     for(int i=0;i<d.x;i++)
		 for(int j=0;j<d.y;j++)
			 Pcoeff[i][j]=coeff[i][j];

}/**/

//----------------------------------------------------------------------

void Image_BW::Desend_coeff1(Image_Coord  d, Pel_Type ** Pcoeff)
{
//     assign_mem(d, 1);
     for(int i=0;i<d.x;i++)
		 for(int j=0;j<d.y;j++)
			 Pcoeff[i][j]=coeff1[i][j];

}

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

double Image_BW::compare(int num, char * file_name)
{
//  if (levels) Error("cannot compare < Image_BW >");

  FILE * in_file = fopen(file_name, "rb");
  if (in_file == NULL) Error(R_MSG);

  int check=fseek(in_file,(dim.x*dim.y*num),SEEK_SET);//*3/2
  if(check)
	  printf("can't open the %c",file_name);

  double mse = 0.0;
  int i, j, k, p, c, t;
  for (i = 0; i < dim.x; i++)
    for (j = 0; j < dim.y; j++) {
#ifdef RESIDUAL_IMAGE
      coeff[i][j]+=128;
#endif
#ifdef LOSSLESS
      t = coeff[i][j];
#else
      t = int(floor(0.5 + coeff[i][j]));
#endif
      if (t < 0) t = 0;
      if ((bytes == 1) && (t > 255)) t = 255;
      for (p = k = 0; k < bytes; k++) {
        if ((c = getc(in_file)) == EOF) Error(R_MSG);
        p = (p << 8) | c; }
      mse += Sqr(p - t); }
//  if (getc(in_file) != EOF) Error(L_MSG);
  fclose(in_file);

  return (mse / dim.x) / dim.y;
}

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

void Image_BW::write_pic(char * file_name)
{
  if (levels) Error("cannot write < Image_BW >");

  FILE * out_file = fopen(file_name, "wb");
  if (out_file == NULL) Error(W_MSG);

  int i, j, k;
  for (i = 0; i < dim.x; i++)
    for (j = 0; j < dim.y; j++) {
#ifdef RESIDUAL_IMAGE
      coeff[i][j]+=128;
#endif
#ifdef LOSSLESS
      k = coeff[i][j];
#else
      k = int(floor(0.5 + coeff[i][j]));
#endif
      if (k < 0) k = 0;
      if (bytes == 2) {
        if (putc(k >> 8, out_file) == EOF) Error(W_MSG); }
      else {
        if (k > 255) k = 255; }
      if (putc(k & 0xFF, out_file) == EOF) Error(W_MSG); }

  fclose(out_file);
}

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

void Image_BW::reset(Image_Coord d, long m, int b)
{
  assign_mem(d, b);
  bytes = b;  mean = m;  levels = Min(max_levels(d.x), max_levels(d.y));

  int i, j;
  for (i = 0; i < dim.x; i++)
	  for (j = 0; j < dim.y; j++) coeff[i][j] = 0; 
}

//----------------------------------------------------------------------

void Image_BW::lvset(void)
{
	levels=-1;
}

//----------------------------------------------------------------------

void Image_BW::lvset(int level1)
{
	levels=level1;
}

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

void Image_BW::transform(void)
{
  if (levels) Error("cannot transform < Image_BW >");

  Chronometer cpu_time;
  cpu_time.start("\n  Starting image transformation...");

  levels = Min(max_levels(dim.x), max_levels(dim.y));
  if (levels < 1) Error("invalid < Image_BW > dimension (transform)");
  printf("Levels=%d\n",levels);

#ifdef LOSSLESS
  int i = 0, j = Max(dim.x, dim.y), k = j << 1;
#else
  int i = NumbTap, j = 0, k = Max(dim.x, dim.y) + (i << 1);
#endif
  CREATE_VECTOR(temp_line, k, Pel_Type, M_MSG);
  Pel_Type * t, * in_line = temp_line + i, * out_line = in_line + j;

// hierarchical subband or S+P transformation

  int nx, ny, mx = dim.x, my = dim.y;
  for (int lv = 0; lv < levels; lv++) {

  // shifts are halved

    nx = mx;  mx >>= 1;  ny = my;  my >>= 1;

  // transformation of columns

    for (j = 0; j < ny; j++) {
      for (i = 0; i < nx; i++) in_line[i] = coeff[i][j];
#ifdef LOSSLESS
      SP_Transform(mx, in_line, out_line, out_line + mx);
      for (i = 0; i < nx; i++) coeff[i][j] = out_line[i]; }
#else
      Reflection(in_line, in_line + nx - 1);
      for (i = 0, t = in_line; i < mx; i++) {
        coeff[i][j] = Filter_L(T_LowPass, t++);
        coeff[i+mx][j] = Filter_H(T_HighPass, t++); } }
#endif

  // transformation of rows

    for (i = 0; i < nx; i++) {
      memcpy(in_line, coeff[i], ny * sizeof(Pel_Type));
#ifdef LOSSLESS
      SP_Transform(my, in_line, coeff[i], coeff[i] + my); } }
#else
      Reflection(in_line, in_line + ny - 1);
      for (j = 0, t = in_line; j < my; j++) {
        coeff[i][j] = Filter_L(T_LowPass, t++);
        coeff[i][j+my] = Filter_H(T_HighPass, t++); } } }
#endif

// calculate and subtract mean

  Pel_Type s;
  for (s = i = 0; i < mx; i++)
    for (j = 0; j < my; j++) s += coeff[i][j];
  mean = Round((s / double(mx)) / my);
//  printf("Mean=%d, mx=%d, my=%d\n",mean,mx,my);
  for (i = 0; i < mx; i++)
    for (j = 0; j < my; j++) coeff[i][j] -= mean;

  delete [] temp_line;
  cpu_time.display(" Image transformed in");
}


void Image_BW::transform(Pel_Type ** coeff)
{
  levels =Min(max_levels(dim.x), max_levels(dim.y)) ;//3


#ifdef LOSSLESS
  int i = 0, j = Max(dim.x, dim.y), k = j << 1;
#else
  int i = NumbTap, j = 0, k = Max(dim.x, dim.y) + (i << 1);
#endif
  CREATE_VECTOR(temp_line, k, Pel_Type, M_MSG);
  Pel_Type * t, * in_line = temp_line + i, * out_line = in_line + j;

// hierarchical subband or S+P transformation

  int nx, ny, mx = dim.x, my = dim.y;
  for (int lv = 0; lv < levels; lv++) {

  // shifts are halved

    nx = mx;  mx >>= 1;  ny = my;  my >>= 1;

  // transformation of columns

    for (j = 0; j < ny; j++) {
      for (i = 0; i < nx; i++) in_line[i] = coeff[i][j];
#ifdef LOSSLESS
      SP_Transform(mx, in_line, out_line, out_line + mx);
      for (i = 0; i < nx; i++) coeff[i][j] = out_line[i]; }
#else
      Reflection(in_line, in_line + nx - 1);
      for (i = 0, t = in_line; i < mx; i++) {
        coeff[i][j] = Filter_L(T_LowPass, t++);
        coeff[i+mx][j] = Filter_H(T_HighPass, t++); } }
#endif

  // transformation of rows