watson_metric_pgm.c

waston human visual model source code 内容概要: 人类感知模型原码！

/* Computes perceptual error. The perceptual error is computed using the
 * Watson model which takes into account three factors: contrast
 * sensitivity, luminance masking and contrast masking.
 * For a description of the Watson model see section 2.3 of the paper
 * "A comparison of image quality models and metrics based on human
 * visual sensitivity".
 */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <getopt.h>
#include <math.h>
#include <float.h>
#include <sys/param.h>
#include <dct.h>
extern "C" {
#include "pgm.h"
} 

char *progname;

#define sqr(X) ((X) * (X))

double **alloc_double(int width, int height) {
  int i;

  double **a = (double **) malloc(height * sizeof(double *));

  for (i = 0; i < height; i++)
    a[i] = (double *) malloc(width * sizeof(double));

  return a;
}

void free_double(double **a, int height) {
  int i;

  for (i = 0; i < height; i++)
    free(a[i]);

  free(a);
}

double **create_visibility_threshold(double viewing_distance, double image_width, double image_height, int pixel_width, int pixel_height, int blocksize) {
  int i, j;

  // image_width .. horizontal image size in centimeters
  // pixel_width .. no. of pixels in image_width
  // viewing_distance .. viewing distance in centimeters

  // compute pixel size in radians
  double Wx = (image_width / pixel_width) / viewing_distance;
  double Wy = (image_height / pixel_height) / viewing_distance;

  // convert pixel size in radians to size in degrees
  Wx = Wx * 180.0 / M_PI;
  Wy = Wy * 180.0 / M_PI;

  // compute mean luminance in candels per square meter according to an exponential law
  double L = 240.0; // mean luminance in gray-levels, was 128
  L = pow(10, (7.0 / 255.0 * L - 4.0));

  // Compute luminance functions i.e. Tmin, fmin and K.
  double LT = 13.45, So = 94.7, at = 0.649;
  double Lk = 300, ko = 3.125, ak = 0.0706;
  double Lf = 300, fo = 6.78, af = 0.182;

  // Minimal threshold function i.e. Tmin.
  double Tmin;
  if (L > LT)
    Tmin = L/So;
  else
    Tmin = L/So*pow(LT / L, 1 - at);

  // Parabola stepness function i.e. K.
  double parabola_stepness;
  if (L > Lk)
    parabola_stepness = ko;
  else
    parabola_stepness = ko * pow(L / Lk, ak);
  
  // minimal frequency function i.e. fmin
  double fmin;
  if (L > Lf)
    fmin = fo;
  else
    fmin = fo * pow(L / Lf, af);  

  double **thres = alloc_double(blocksize, blocksize);

  for (i = 0; i < blocksize; i++)
    for (j = 0; j < blocksize; j++) {
      double fio = i;
      double foj = j;
      fio /= (double) (2.0 * blocksize);
      foj /= (double) (2.0 * blocksize);
      fio /= Wx;
      foj /= Wy;
 
      // compute spatial frequencies
      double f = sqrt(sqr(fio) + sqr(foj));

      // compute angular parameters
      double teta = 2.0 * fio * foj;
      if (i > 0 && j > 0) {       
        teta /= sqr(f);
      }
      if (teta > 1.0) teta = 1.0;
      teta = asin(teta);
  
      // determine obliqueness effect values
      double R = 0.7;
      double r = R + (1 - R) * sqr(cos(teta));

      // Compute luminance thresholds i.e. thres.
      // f[0][0] is 0 but to don't have log(0) we make f[0][0] = 1.


      if (i == 0 && j == 0)
        f = 1;

      thres[i][j] = 
        log10(Tmin / r) + 
        parabola_stepness * pow(log10(f / fmin), 2);
      thres[i][j] = pow(10, thres[i][j]);

      // According to our experiments we can multiply the thresholds
      // until a factor of 3.5 so that they are not yet visibly after luminance
      // masking corrections.
      thres[i][j] = 2.2 * thres[i][j];

    }

  // Approximate luminance threshold for DC coefficient
  thres[0][0] = MIN(thres[0][1], thres[1][0]);

  return thres;
}

double **contrast_masking(double **block, double **thres, int blocksize)
{
  int i, j;

  // Contrast masking degree i.e. w.
  double w = 0.7;

  // Compute the mean luminance in DCT domain.
  double mean_lum = 128 * blocksize;

  // Compute luminance masking threshods
  double lum_thres = fabs(block[0][0]);

  if (lum_thres == 0 || lum_thres == blocksize)
    return NULL;

  double **mask = (double **) alloc_double(blocksize, blocksize);

  for (i = 0; i < blocksize; i++) {
    for (j = 0; j < blocksize; j++) {
  
      // the factor 2.2 was added to increase the luminance masking
      double lumthres = 3.2 * thres[i][j] * pow(lum_thres/mean_lum, 0.65);
      if (lum_thres < 600)
        lumthres = lumthres * 2.2;  // very dark areas can be marked a bit more

      double actthres = pow(abs(block[i][j]), w) * pow(lumthres, (1 - w));
      if (i == 0 && j == 0)
        actthres = lumthres;

      // Compute contrast masking thresholds
      mask[i][j] = lumthres;
      if (lumthres < actthres)
        mask[i][j] = actthres;
    }
  }

  return mask;
}

void print_msg(char *msg) {
  fprintf(stderr, "%s", msg);
}

void print_double(double **a, int width, int height) {
  int i, j;
  
  for (i = 0; i < height; i++) {
    for (j = 0; j < width; j++) {
      fprintf(stderr, "%10.4f ", a[i][j]);
    }
    fprintf(stderr, "\n");
  }
}

void usage(void) {
  fprintf(stderr, "usage: %s [-h] [-i file] [-o file] file\n", progname);
  fprintf(stderr, "\t-h\t\tprint usage\n"); 
  fprintf(stderr, "\t-i file\t\toriginal image\n");
  fprintf(stderr, "\t-o file\t\toutput image file\n");
  exit(0);
}

int main (int argc, char *argv[]) {
  FILE *in = stdin;
  FILE *orig = NULL;
  FILE *out = stdout;

  char input_name[MAXPATHLEN] = "(stdin)";
  char orig_name[MAXPATHLEN] = "(none)";
  char output_name[MAXPATHLEN] = "(stdout)";

  int in_cols, in_rows, in_format;
  gray in_maxval;
  int orig_cols, orig_rows, orig_format;
  gray orig_maxval;
  int rows, cols, format;
  gray maxval;
  int row, col;
  gray **input_image, **orig_image, **output_image;

  int i, j, c, n;

  // Local perceptual error thresholds
  int nLPE1 = 0, nLPE2 = 0;
  double maxLPE = 0.0;
  double LPEThreshold1 = 0.1; 
  double LPEThreshold2 = 0.12;
  // blocksize
  int blocksize = 8;

  progname = argv[0];

  pgm_init(&argc, argv);
  
  while ((c = getopt(argc, argv, "i:h?o:")) != EOF) {
   switch (c) {
      case 'h':
      case '?':
        usage();
        break; 
      case 'i':
        if ((orig = fopen(optarg, "rb")) == NULL) {
           fprintf(stderr, "%s: unable to open original image file %s\n", progname, optarg);
           exit(1);
        }
        strcpy(orig_name, optarg);
        break;
      case 'o':
        if ((out = fopen(optarg, "wb")) == NULL) {
          fprintf(stderr, "%s: unable to open output file %s\n", progname, optarg);
          exit(1);
        }
        strcpy(output_name, optarg);
        break; 
    }
  }

  argc -= optind; 
  argv += optind;
          
  if (argc > 1) { 
    usage();
    exit(1);
  }

  if (argc == 1 && *argv[0] != '-')
    if ((in = fopen(argv[0], "rb")) == NULL) {
      fprintf(stderr, "%s: unable to open input file %s\n", progname, argv[0]);
      exit(1);
    }
    else  
      strcpy(input_name, argv[0]);

  pgm_readpgminit(in, &in_cols, &in_rows, &in_maxval, &in_format);
  pgm_readpgminit(orig, &orig_cols, &orig_rows, &orig_maxval, &orig_format);

  if (in_cols != orig_cols || in_rows != orig_rows) {
    fprintf(stderr, "%s: input image %s does not match dimensions of original image %s\n", progname, input_name, orig_name);
    exit(1);
  }

  if (cols % NJPEG) {
    fprintf(stderr, "%s: image width %d not a multiple of %d\n", progname, cols, NJPEG);
    exit(1);
  }
      
  if (rows % NJPEG) {
    fprintf(stderr, "%s: image height %d not a multiple of %d\n", progname, rows, NJPEG);
    exit(1);
  }

  cols = in_cols;
  rows = in_rows;
  format = in_format;
  maxval = in_maxval;

  input_image = pgm_allocarray(cols, rows);
  orig_image = pgm_allocarray(cols, rows);

  for (row = 0; row < rows; row++) {
    pgm_readpgmrow(in, input_image[row], cols, in_maxval, in_format);
    pgm_readpgmrow(orig, orig_image[row], cols, orig_maxval, orig_format);
  }

  fclose(in);
  fclose(orig);

  double **thres = create_visibility_threshold(72, 16, 16, 512, 512, NJPEG);

  init_dct_8x8();
  double **input_dct8x8 = alloc_coeffs_8x8();
  double **orig_dct8x8 = alloc_coeffs_8x8();

  output_image = pgm_allocarray(cols, rows);

  double globalsum = 0.0;
  for (i = 0; i < rows / NJPEG; i++) {
    for (j = 0; j < cols / NJPEG; j++) {
      fdct_block_8x8(input_image, j * NJPEG, i * NJPEG, input_dct8x8);
      fdct_block_8x8(orig_image, j * NJPEG, i * NJPEG, orig_dct8x8);

      double **mask = contrast_masking(orig_dct8x8, thres, NJPEG);
      double localsum = 0.0;
      int k, l;
      for (k = 0; k < NJPEG; k++) {
        for (l = 0; l < NJPEG; l++) {
          double error = fabs(input_dct8x8[k][l] - orig_dct8x8[k][l]);
          double perceptual_error = error / mask[k][l];
          localsum += perceptual_error;
          globalsum += perceptual_error;
        }
      }   
      localsum /= sqr(NJPEG);
      double LPE = localsum;
// fprintf(stderr, "(%d/%d) LPE: %f\n", j*NJPEG, i*NJPEG, LPE);
      if (LPE > LPEThreshold1) nLPE1++;      
      if (LPE > LPEThreshold2) nLPE2++;      
      if (LPE > maxLPE) maxLPE = LPE;

      for (k = 0; k < NJPEG; k++) {
        for (l = 0; l < NJPEG; l++) {
          output_image[i * NJPEG + k][j * NJPEG + l] = PIXELRANGE(LPE * (255 / LPEThreshold2));
        }
      }
free_double(mask, NJPEG);
    }
  }
  globalsum /= (cols * rows);
  double GPE = globalsum;

  fprintf(stderr, "#LPE1: %d\n", nLPE1);
  fprintf(stderr, "#LPE2: %d\n", nLPE2);
  fprintf(stderr, "max. LPE: %f\n", maxLPE);
  fprintf(stderr, "GPE (global perceptual error): %f\n", GPE);

  pgm_freearray(input_image, rows);
  pgm_freearray(orig_image, rows);

  pgm_writepgminit(out, cols, rows, maxval, 0);
  for (row = 0; row < rows; row++)
    pgm_writepgmrow(out, output_image[row], cols, maxval, 0);
    
  fclose(out);
        
  pgm_freearray(output_image, rows);
    
  exit(0);
}