dctfilter.c

// -----------------------------------------------------------------
// Overlapping 2D-DCT (Discrete Cosine Transform) image filter
//
// Image model: z = y + n, where z is the measured image mxn
//                               y is the noise free image
//                               n is Gaussian noise
// 
// Compile with: g++ -O3 dctfilter.c dct.c -o dctfilter
// 
// Mikael Mieskolainen, {2011, 2021}
// ----------------------------------------------------------------- 

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <cstring>

#include "dct.h" // 2D-DCT Library

// This abstracts the image matrix as a linear array
// access image matrix pixel at row i and column j, with N columns
// The first for "->" and the second for "." access
#define pixel(data,i,j)   ( (data)->image[(i)*((data)->N) + (j)] )
#define pixel_h(data,i,j) ( (data).image[(i)*((data).N) + (j)] )

// This is the image object
struct IMAGE {

  IMAGE(int m, int n) {
    M = m;
    N = n;
  };

  double *image; // Data
  int M;         // Number of Rows
  int N;         // Number of Columns
};

// Subfunction prototypes
void ones(struct IMAGE*);
void zeros(struct IMAGE*);
void printValues(struct IMAGE*, int jump=1);
void pointDiv(struct IMAGE*, struct IMAGE*, struct IMAGE*);
void pointMul(struct IMAGE*, struct IMAGE*, struct IMAGE*);
bool dctfilter(struct IMAGE*, struct IMAGE*, int step, double thresh, int window_size);
double thrfunc(double* outblock, int window_size, double thresh);


// Main program
int main(int argc, char* argv[]) {

  char* inputfile;
  char* outputfile;

  // Image dimensions fixed here (no header in the binary files)
  int M = 512;              // Rows
  int N = 512;              // Cols
  
  // Filter parameters
  int step         = 1;     // Local window hopping step
  double threshold = 40;    // Frequency domain threshold
  int bytes        = 8;     // 64 bit (double) input expected (more formats TBD!)
  bool quiet       = false; // Print output

  int opt;
  while ((opt = getopt(argc, argv, "i:o:b:m:n:s:t:q")) != -1) {  
    switch(opt) {
      case 'i':
        inputfile = strdup(optarg);
        break;
      case 'o':
        outputfile = strdup(optarg);
        break;
      case 'b':
        bytes = atoi(optarg);
        if (bytes < 1) {
          printf("-b needs to be int >= 1 \n");
          return EXIT_FAILURE;
        }
        break;
      case 'm':
        M = atoi(optarg);
        if (M < 1) {
          printf("-m needs to be integer >= 1 \n");
          return EXIT_FAILURE;
        }
        break;
      case 'n':
        N = atoi(optarg);
        if (N < 1) {
          printf("-n needs to be integer >= 1 \n");
          return EXIT_FAILURE;
        }
        break;
      case 's':
        step = atoi(optarg);
        if (step < 1) {
          printf("-s needs to be integer >= 1 \n");
          return EXIT_FAILURE;
        }
        break;
      case 't':
        threshold = atof(optarg);
        if (threshold < 0) {
          printf("-t needs to be float >= 0 \n");
          return EXIT_FAILURE;
        }
        break;
      case 'q':
        quiet = true;
        break;
      case ':':
        printf("option needs a value \n");
        break;
      case '?':
        fprintf(stderr, "Usage: %s [-i inputfile] [-o outputfile] [-m nrows] [-n ncols] [-s step] [-t threshold] [-q quiet] \n", argv[0]);
        return EXIT_FAILURE;
    }
  }

  if (!quiet) {
    printf("-i inputfile:  %s \n", inputfile);  
    printf("-o outputfile: %s \n", outputfile);
    printf("-b bytes:      %d \n", bytes);
    printf("-m nrows:      %i \n", M);
    printf("-n ncols:      %i \n", N);
    printf("-s step:       %i \n", step); 
    printf("-t threshold:  %0.1f \n", threshold);
  }

  const int window_size = 8;  // Keep it at 8 (2D-DCT function is fixed 8x8)

  // ------------------------------------------------
  // Create image matrix for the original (noisy) image
  // and denoised image

  struct IMAGE z(M,N);
  struct IMAGE y_est(M,N);

  // ------------------------------------------------
  // Allocate memory

  z.image = (double*) malloc(z.M * z.N * sizeof(double));
  y_est.image = (double*) malloc(y_est.M * y_est.N * sizeof(double));
  
  if (z.image == NULL) {
    printf("dctfilter: Out of memory with z \n");
    return EXIT_FAILURE;
  }
  if (y_est.image == NULL) {
    printf("dctfilter: Out of memory with y_est \n");
    return EXIT_FAILURE;
  }

  // Init matrices with zeros
  zeros(&z);
  zeros(&y_est);
  
  // ------------------------------------------------
  // Read the image
  
  FILE *fid;
  fid = fopen(inputfile, "rb");

  if (fid) {
    int ret = fread(z.image, bytes, z.M * z.N, fid);
    fclose(fid);
  } else {
    printf("dctfilter: Error opening inputfile <%s> ! \n", inputfile);
    return EXIT_FAILURE;
  }
  
  // ------------------------------------------------
  // Filter the image

  if (dctfilter(&z, &y_est, step, threshold, window_size)) {
    // fine
  } else {
    printf("dctfilter: Problem in filtering \n");
    return EXIT_FAILURE;
  }
  
  // ------------------------------------------------
  // Write the image
  fid = fopen(outputfile, "wb");
  if (fid) { 
    int ret = fwrite(y_est.image, bytes, y_est.M * y_est.N, fid);
    fclose(fid);
  } else {
    printf("dctfilter: Error opening outputfile <%s> ! \n", outputfile);
    return EXIT_FAILURE;
  }

  // Free the image matrices
  free(z.image);
  free(y_est.image);
  
  if (!quiet) { printf("dctfilter: Filtering done \n"); }

  return EXIT_SUCCESS;
}


// 2D-filter function
bool dctfilter(struct IMAGE* z, struct IMAGE* y_est, int step, double thresh, int window_size) {
  
  // Boundary around image matrix
  const int b = window_size / 2;

  // 2D-DCT 8x8 buffers (matrix row by row as linear array)
  double inblock[(window_size*window_size)]  = {0};
  double outblock[(window_size*window_size)] = {0};

  // Allocate the aggregation matrix and weight matrix
  struct IMAGE A(z->M, z->N);
  struct IMAGE W(z->M, z->N);

  A.image = (double*) malloc(A.M * A.N * sizeof(double));
  W.image = (double*) malloc(W.M * W.N * sizeof(double));
  
  if (A.image == NULL) { printf("dctfilter: Out of memory with A \n"); return false; }
  if (W.image == NULL) { printf("dctfilter: Out of memory with W \n"); return false; }
  
  // Init matrices with zeros
  zeros(&A);
  zeros(&W);
  
  // Over all rows and columns of the image
  for (int i = b; i < z->M - b + 2; i += step) {
    for (int j = b; j < z->N - b + 2; j += step) {

      double* VP = &inblock[0]; // Current pointer

      // Accumulate the local window
      for (int m = i-b, k = 0; m < i+b; ++m) {
        for (int n = j-b; n < j+b; ++n, ++k) {

          // Check boundary
          *(VP + k) = (m < z->M && n < z->N) ? pixel(z,m,n) : 0.0;
        }
      }

      // 2D-DCT for the block
      fdct(VP, outblock);
      
      // Filtering in frequency domain & weight calculation
      const double weight = thrfunc(outblock, window_size, thresh);
      
      // Inverse 2D-DCT for the block
      idct(outblock, VP);

      // Aggregate the block and weight
      for (int m = i-b, k = 0; m < i+b; ++m) {
        for (int n = j-b; n < j+b; ++n, ++k) {

          if (m < z->M && n < z->N) { // Check boundary
            pixel_h(A,m,n) += weight * (*(VP + k));
            pixel_h(W,m,n) += weight;
          }
        }
      }
    }
  }
  
  // Weighted mean of the aggregation buffer
  for (int i = 0; i < y_est->M; ++i) {
    for (int j = 0; j < y_est->N; ++j) {
      pixel(y_est,i,j) = pixel_h(A,i,j) / pixel_h(W,i,j);
    }  
  }
  
  free(A.image);
  free(W.image);

  return true;
}

// Hard-thresholding and block weight function
double thrfunc(double* outblock, int window_size, double thresh) {

  double sum  = 0.0;
  double sum2 = 0.0;
  for (int i = 0; i < (window_size*window_size); ++i, ++outblock) {

    const double c = fabs(*outblock);

    if (c < thresh) {
      *outblock = 0;
    } else {
      sum  += c;   // Accumulate coefficient values
      sum2 += c*c;
    }
  }
  const double n   = window_size*window_size;
  const double var = (sum2 - sum*sum/n) / n;
  
  // Weight for this block
  const double W = var > 0 ? 1.0/var : 1.0;

  return W;
}

// Point wise division: y = a ./ b
void pointDiv(struct IMAGE *y, struct IMAGE *a, struct IMAGE *b) {
  
  for (int i = 0; i < y->M; ++i) {
    for (int j = 0; j < y->N; ++j) {
      pixel(y,i,j) = pixel(a,i,j) / pixel(b,i,j);
    }  
  }   
}

// Point wise multiplication: y = a .* b
void pointMul(struct IMAGE *y, struct IMAGE *a, struct IMAGE *b) {
  
  for (int i = 0; i < y->M; ++i) {
    for (int j = 0; j < y->N; ++j) {
      pixel(y,i,j) = pixel(a,i,j) * pixel(b,i,j);
    }
  }
}

// Init with ones
void ones(struct IMAGE *x) {
  
  for (int i = 0; i < x->M; ++i) {
    for (int j = 0; j < x->N; ++j) {
      pixel(x,i,j) = 1.0;
    }
  }
}

// Init with zeros
void zeros(struct IMAGE *x) {
  
  for (int i = 0; i < x->M; ++i) {
    for (int j = 0; j < x->N; ++j) {
      pixel(x,i,j) = 0.0;
    }
  }
}

void printValues(struct IMAGE *x, int jump) {
  
  for (int i = 0; i < x->M; i += jump) {
    for (int j = 0; j < x->N; j += jump) {
      printf("Value at [%d,%d] is %0.2f \n", i, j, pixel(x,i,j) );
    }
  }
}