image_3c.c

SPIHT 源代码

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

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

// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
//           > > > >    C++ version  1.04 -  06/25/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_3c.h"
#include <string.h>


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

static char * M_MSG = "< Image_3C >";

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

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

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

static const int NumbTap = 4;

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

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

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

static const float R_HighPass[5] =
  { 0.788485, -0.377403, -0.040690,  0.023849, 0.0 };


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

//  Auxiliary functions

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

inline float Filter_L(const float * f, float * 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 float Filter_H(const float * f, float * 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(float * h, float * t)
{  
  for (int i = 1; i <= NumbTap; i++) {
    h[-i] = h[i];  t[i] = t[-i]; }
}

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

//  Functions of the class  < Image_3C >

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

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

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

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

void Image_3C::assign_mem(Image_Coord d)
{
  if ((levels >= 0) && (dim.x == d.x) && (dim.y == d.y)) return;
  free_mem();
  if ((d.x < 8) || (d.y < 8))
    Error("< Image_3C > dimension is too small or negative");
  dim = d;
  pdim.x = (d.x < 200 ? (d.x + 7) & 0x3FF8 : (d.x + 15) & 0x3FF0);
  pdim.y = (d.y < 200 ? (d.y + 7) & 0x3FF8 : (d.y + 15) & 0x3FF0);

  int c, i = NumbTap + NumbTap + Max(pdim.x, pdim.y);
  NEW_VECTOR(r_line, i, float, M_MSG);
  t_line = r_line + NumbTap;
  for (c = 0; c < 3; c++) {
    NEW_VECTOR(coeff, pdim.x, float *, M_MSG);
    for (i = 0; i < pdim.x; i++) {
      NEW_VECTOR(coeff[i], pdim.y, float, M_MSG); }
    cor[c] = coeff; }
}

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

void Image_3C::free_mem(void)
{
  if (levels >= 0) {
    int c, i;
    for (c = 2; c >= 0; c--) {
      coeff = cor[c];
      for (i = pdim.x - 1; i >= 0; i--) delete [] coeff[i];
      delete [] coeff; }
    delete [] r_line; }
  dim.x = dim.y = 0;  levels = -1;
}

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

void Image_3C::extend_image(int c)
{
  coeff = cor[c];
  int i, j;
  for (j = dim.y - 1; j < pdim.y - 1; j++) {
    coeff[0][j+1] = (coeff[0][j] + coeff[1][j]) / 2;
    coeff[dim.x-1][j+1] = (coeff[dim.x-1][j] + coeff[dim.x-2][j]) / 2;
    for (i = dim.x - 2; i > 0; i--)
      coeff[i][j+1] = (coeff[i-1][j] + coeff[i][j] + coeff[i+1][j]) / 3; }
  for (i = dim.x - 1; i < pdim.x - 1; i++) {
    coeff[i+1][0] = (coeff[i][0] + coeff[i][1]) / 2;
    coeff[i+1][pdim.y-1] = (coeff[i][pdim.y-1] + coeff[i][pdim.y-2]) / 2;
    for (j = pdim.y - 2; j > 0; j--)
      coeff[i+1][j] = (coeff[i][j-1] + coeff[i][j] + coeff[i][j+1]) / 3; }
}

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

#define ROTATION(M,i,j,k,l) g = M[i][j];  h = M[k][l];\
         M[i][j] = g-s*(h+g*tau);  M[k][l] = h+s*(g-h*tau)

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

void Image_3C::Jacobi(double m[3], double A[3][3])
{
  int i, j, k, p, q, itr = 0;
  double t = 1.0 / (double(dim.x) * double(dim.y));
  double v, tau, s, h, g, c;
  double b[3], d[3], X[3][3];

  for (i = 0; i < 3; i++) m[i] *= t;
  for (i = 0; i < 3; i++) {
    for (j = i; j < 3; j++) {
      X[i][j] = X[j][i] = 0;
      A[i][j] = A[j][i] = t * A[i][j] - m[i] * m[j]; }
    mean[i] = int(0.5 + m[i]);  m[i] = 0;
    b[i] = d[i] = A[i][i];  X[i][i] = 1; }

  for (;;) {

    for (s = i = 0; i < 2; i++)
      for (j = i + 1; j < 3; j++)
        if (s < (t = fabs(A[i][j]))) {
          s = t;  p = i;  q = j; }

    if (s == 0) break;

    if (++itr > 50)
      Error("Too many iterations in function < Image_3C::Jacobi >");

    g = 100 * s;  h = d[q] - d[p];
    if (g + fabs(h) == fabs(h))
      t = A[p][q] / h;
    else {
      v = 0.5 * h / A[p][q];  t = 1 / (fabs(v) + sqrt(1 + v * v));
      if (v < 0) t = -t; }
    c = 1 / sqrt(1 + t * t);  s = t * c;
    tau = s / (1 + c);  h = t * A[p][q];  A[p][q] = 0;
    m[p] -= h;  m[q] += h;  d[p] -= h;  d[q] += h;

    for (k = 0; k < p; k++)     { ROTATION(A,k,p,k,q); }
    for (k = p + 1; k < q; k++) { ROTATION(A,p,k,k,q); }
    for (k = q + 1; k < 3; k++) { ROTATION(A,p,k,q,k); }
    for (k = 0; k < 3; k++)     { ROTATION(X,k,p,k,q); }

    for (k = 0; k < 3; k++) {
      b[k] += m[k];  d[k] = b[k];  m[k] = 0; }
  }

  for (i = 0; i < 3; i++) {
    for (s = k = -1, j = 0; j < 3; j++)
      if (d[j] > s) s = d[k=j];
    if (k < 0) Error("Incorrect decomposition");
    for (d[k] = -2, j = 0; j < 3;j++) 
      R[i][j] = 1e-4 * Round(1e4 * X[j][k]); }
}

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

void Image_3C::averages(int M[3])
{
  for (int i = 0; i < 3; i++) M[i] = mean[i];
}

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

void Image_3C::color_transform(int T[3][3])
{
  int i, j;
  for (i = 0; i < 3; i++)
    for (j = 0; j < 3; j++)
      T[i][j] = int(Round(1e4 * (R[i][j] + 1)));
}

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

void Image_3C::set_color(int c)
{
  if ((c < 0) || (c > 3)) Error("invalid < Image_3C > color number");
  coeff = cor[c];
}

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

void Image_3C::read_pic(Image_Coord d, char * file_name)
{
  assign_mem(d);  levels = 0;

  int i, j, k, l, v, p[3];
  double m[3], S[3][3];
  FILE * in_file = fopen(file_name, "rb");
  if (in_file == NULL) Error(R_MSG);

  for (k = 0; k < 3; k++)
    for (m[k] = l = 0; l < 3; l++) S[k][l] = 0;

  for (i = 0; i < dim.x; i++)
    for (j = 0; j < dim.y; j++) {
      for (k = 0; k < 3; k++) {
        if ((v = getc(in_file)) == EOF) Error(R_MSG);
        cor[k][i][j] = p[k] = v;  m[k] += v; }
      for (k = 0; k < 3; k++)
        for (l = k; l < 3; l++) S[k][l] += p[k] * p[l]; }
  if (getc(in_file) != EOF) Error(L_MSG);
  fclose(in_file);

  Jacobi(m, S);

  for (k = 0; k < 3; k++) extend_image(k);
}

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

float Image_3C::compare(char * file_name)
{
  if (levels) Error("cannot compare < Image_3C >");

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

  double mse = 0.0;
  int i, j, k, c;
  for (i = 0; i < dim.x; i++)
    for (j = 0; j < dim.y; j++)
      for (k = 0; k < 3; k++) {
        if ((c = getc(in_file)) == EOF) Error(R_MSG);
        mse += Sqr(c - int(0.499 + cor[k][i][j])); }
  if (getc(in_file) != EOF) Error(L_MSG);
  fclose(in_file);

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

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

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

  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++)
      for (k = 0; k < 3; k++) {
      if (putc(int(0.499 + cor[k][i][j]), out_file) == EOF)
        Error(W_MSG); }

  fclose(out_file);
}

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

void Image_3C::reset(Image_Coord d)
{
  assign_mem(d);

  int i, k, j = pdim.y * sizeof(float);
  for (i = 0; i < 3; i++)
    for (mean[i] = j = 0; j < 3; j++) R[i][j] = 0;

  levels = Min(max_levels(pdim.x), max_levels(pdim.y));
//  for (k = 0; k < 3; k++)
//    for (i = 0; i < pdim.x; i++) memset(cor[k][i], 0, j);
  for (k = 0; k < 3; k++)
    for (i = 0; i < pdim.x; i++)
      for (j = 0; j < pdim.y; j++) cor[k][i][j] = 0;
}

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

void Image_3C::reset(Image_Coord d, int m[3], int T[3][3])
{
  reset(d);

  int i, j;
  for (i = 0; i < 3; i++) {
    mean[i] = m[i];
    for (j = 0; j < 3; j++) R[i][j] = 1e-4 * T[i][j] - 1.0; }
}

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

void Image_3C::transform(void)
{
  if (levels) Error("cannot transform < Image_3C >");
  if ((dim.x < 32) || (dim.y < 32))
    Error("< Image_3C > is too small to transform");

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

  levels = Min(max_levels(pdim.x), max_levels(pdim.y));
  if (levels < 2) Error("invalid < Image_3C > dimension");

  int i, j, k, l;
  float s, p[3];

  for (i = 0; i < pdim.x; i++)
    for (j = 0; j < pdim.y; j++) {
      for (k = 0; k < 3; k++) p[k] = cor[k][i][j] - mean[k];
      for (k = 0; k < 3; k++) {
        for (s = l = 0; l < 3; l++) s += R[k][l] * p[l];
        cor[k][i][j] = s; } }

  for (k = 0; k < 3; k++) w_transform(cor[k]);

  cpu_time.display(" Image transformed in");
}

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

void Image_3C::w_transform(float ** tcf)
{
  int i, j, nx, ny, mx = pdim.x, my = pdim.y;
  float * t;

  for (int lv = 0; lv < levels; lv++) {

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

    for (j = 0; j < ny; j++) {
      for (i = 0; i < nx; i++) t_line[i] = tcf[i][j];
      Reflection(t_line, t_line + nx - 1);
      for (i = 0, t = t_line; i < mx; i++) {
        tcf[i][j] = Filter_L(T_LowPass, t++);
        tcf[i+mx][j] = Filter_H(T_HighPass, t++); } }

    for (i = 0; i < nx; i++) {
      memcpy(t_line, tcf[i], ny * sizeof(float));
      Reflection(t_line, t_line + ny - 1);
      for (j = 0, t = t_line; j < my; j++) {
        tcf[i][j] = Filter_L(T_LowPass, t++);
        tcf[i][j+my] = Filter_H(T_HighPass, t++); } } }
}

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

void Image_3C::recover(void)
{
  if (levels <= 0) Error("cannot recover < Image_3C >");

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

  int i, j, k, l;
  float s, p[3];

  for (k = 0; k < 3; k++) w_recover(cor[k]);

  for (i = 0; i < pdim.x; i++)
    for (j = 0; j < pdim.y; j++) {
      for (k = 0; k < 3; k++) p[k] = cor[k][i][j];
      for (k = 0; k < 3; k++) {
        for (s = l = 0; l < 3; l++) s += R[l][k] * p[l];
        l = int(0.499 + s + mean[k]);
        if (l < 0) l = 0;
        if (l > 255) l = 255;
        cor[k][i][j] = l; } }

  levels = 0;
  cpu_time.display(" Image transformed in");
}

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

void Image_3C::w_recover(float ** tcf)
{
  float * t;
  int i, j, mx, my, nx = pdim.x >> levels, ny = pdim.y >> levels;

  for (int lv = 0; lv < levels; lv++) {

    mx = nx;  nx <<= 1;  my = ny;  ny <<= 1;

    for (i = 0; i < nx; i++) {
      for (j = 0, t = t_line; j < my; j++) {
        *(t++) = tcf[i][j];  *(t++) = tcf[i][j+my]; }
      Reflection(t_line, t_line + ny - 1);
      for (j = 0, t = t_line; j < ny;) {
        tcf[i][j++] = Filter_H(R_HighPass, t++);
        tcf[i][j++] = Filter_L(R_LowPass, t++); } }

    for (j = 0; j < ny; j++) {
      for (i = 0, t = t_line; i < mx; i++) {
        *(t++) = tcf[i][j];  *(t++) = tcf[i+mx][j]; }
      Reflection(t_line, t_line + nx - 1);
      for (i = 0, t = t_line; i < nx;) {
        tcf[i++][j] = Filter_H(R_HighPass, t++);
        tcf[i++][j] = Filter_L(R_LowPass, t++); } } }
}

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

// end of file  < Image_3C.C >