ex66.c

#include <stdio.h>
#include <math.h>
#include <GL/glut.h>

#define PI 3.14159265
#define WD 256
#define HT 256
#define REV 1.0
#define NON_REV -1.0
#define init_mat(mat, w, h) mat.width=w; mat.height=h

/***********************************************:::
    LOW pass & HIGH pass -- reduced order
:::***********************************************/

typedef unsigned char uchar;

typedef struct ImaginaryMatrix{
  double real[WD*HT];
  double im[WD*HT];
  int width;
  int height;
}ImMat;

uchar *im, *result;

void transpose(ImMat *in, ImMat *out){
  int w=in->width, h=in->height, x, y;
  for(y=0; y<h; y++)
    for(x=0; x<w; x++){
      int src=y*w+x, dst=x*w+y;
      out->real[dst]=in->real[src];
      out->im[dst]=in->im[src];
    }
}

void dft(double *in_real, double *in_im, double *out_real, double *out_im, int len, double sign){
  int k,n;
  for(k=0; k<len; k++){
    out_real[k] = 0.;
    out_im[k] = 0.;
    for(n=0; n<len; n++){
      double a = in_real[n], b = in_im[n];
      double x = sign*(2*PI*k*n)/len;
      double cos_ = cos(x), sin_ = sin(x);
      out_real[k] += (a*cos_ - b*sin_);
      out_im[k] += (a*sin_ + b*cos_);
    }
  }
}

void dft2(ImMat *in, ImMat *out){
  ImMat W, tW, tOut;
  int w=in->width, h=in->height, y;
  init_mat(W, w, h); init_mat(tW, w, h); init_mat(tOut, w, h);
  
  for(y=0; y<h; y++)
    dft(&(in->real[y*w]), &(in->im[y*w]), &(W.real[y*w]), &(W.im[y*w]), w, NON_REV);
  transpose(&W, &tW);

  for(y=0; y<h; y++)
    dft(&(tW.real[y*w]), &(tW.im[y*w]), &(tOut.real[y*w]), &(tOut.im[y*w]), w, NON_REV);
  transpose(&tOut, out);
}

void dft2_rev(ImMat *in, ImMat *out){
  ImMat W, tW, tOut;
  int w=in->width, h=in->height, x, y;
  init_mat(W, w, h); init_mat(tW, w, h); init_mat(tOut, w, h);
  
  for(y=0; y<h; y++)
    dft(&(in->real[y*w]), &(in->im[y*w]), &(W.real[y*w]), &(W.im[y*w]), w, REV);
  transpose(&W, &tW);
  for(y=0; y<h; y++)
    dft(&(tW.real[y*w]), &(tW.im[y*w]), &(tOut.real[y*w]), &(tOut.im[y*w]), w, REV);
  transpose(&tOut, out);

  double wh=w*h;
  for(y=0; y<h; y++)
    for(x=0; x<w; x++){
      out->real[y*w+x]/=wh;
      out->im[y*w+x]/=wh;
    }
}

uchar* image_from_bin(char *fname, int size){
  FILE *fd;
  uchar *ret = (uchar*)malloc(sizeof(uchar)*size);
  fd = fopen(fname, "rb");
  fread(ret, 1, size, fd);
  fclose(fd);
  return ret;
}

void display(){
  glDrawPixels(WD, HT, GL_LUMINANCE, GL_UNSIGNED_BYTE, result);
  glutSwapBuffers();
}

void reshape(int w, int h){
  glViewport(0,0,(GLsizei)w, (GLsizei)h);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();
}

void dispose(){
  free(im);
  free(result);
}

void keyboard(uchar key, int x, int y){
  if(key == 'q'){
    dispose();
    exit(0);
  }
}

int main(int argc, char **argv){
  ImMat image, mat;
  int i;
  
  result = (uchar*)malloc(sizeof(uchar)*WD*HT);
  im = image_from_bin("sample.bin", WD*HT);
  
  for(i=0; i<WD*HT; i++){
    image.real[i] = (double)im[i];
    image.im[i] = 0.0;
  }
  
  init_mat(image, WD, HT);
  init_mat(mat, WD, HT);

  dft2(&image, &mat);
  dft2_rev(&mat, &image);
  
  for(i=0; i<WD*HT; i++)result[i]=(uchar)image.real[i];
  
  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
  glutInitWindowSize(WD, HT);
  glutInitWindowPosition(100, 100);
  glutCreateWindow(argv[0]);
  
  glClearColor(0, 0, 0, 0);

  glutDisplayFunc(display);
  glutReshapeFunc(reshape);
  glutKeyboardFunc(keyboard);
  glutMainLoop();

  return 0;
}