acepirin.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <GL/glut.h>
#include "image.h"

#define PI 3.1415926
#define P(x,y) x+y*WD

void dft(double *A, double *re, double *im){
  double x_re[WD*HT], x_im[WD*HT];

  for(int v=0;v<HT;v++)
    for(int u=0;u<WD;u++){
      x_re[P(u,v)]=0.0;
      x_im[P(u,v)]=0.0;
      for(int x=0;x<WD;x++){
        x_re[P(u,v)]+=A[P(x,v)]*cos((-2.0*PI*u*x)/WD);
        x_im[P(u,v)]+=A[P(x,v)]*sin((-2.0*PI*u*x)/WD);
      }
    }

  for(int u=0;u<WD;u++)
    for(int v=0;v<HT;v++){
      re[P(u,v)]=0.0;
      im[P(u,v)]=0.0;
      for(int y=0;y<HT;y++){
        re[P(u,v)]+=x_re[P(u,y)]*cos((-2.0*PI*v*y)/HT)-x_im[P(u,y)]*sin((-2.0*PI*v*y)/HT);
        im[P(u,v)]+=x_re[P(u,y)]*sin((-2.0*PI*v*y)/HT)+x_im[P(u,y)]*cos((-2.0*PI*v*y)/HT);
      }
    }
}

void idft(double *re, double *im, double *A){
  double x_re[WD*HT], x_im[WD*HT];

  for(int v=0;v<HT;v++){
    for(int u=0;u<WD;u++){
      x_re[P(u,v)]=0.0;
      x_im[P(u,v)]=0.0;
      for(int x=0;x<WD;x++){
        x_re[P(u,v)]+=re[P(x,v)]*cos((2.0*PI*u*x)/WD)-im[P(x,v)]*sin((2.0*PI*u*x)/WD);
        x_im[P(u,v)]+=re[P(x,v)]*sin((2.0*PI*u*x)/WD)+im[P(x,v)]*cos((2.0*PI*u*x)/WD);
      }
    }
  }
  
  for(int u=0;u<WD;u++){
    for(int v=0;v<WD;v++){
      A[P(u,v)]=0.0;
      for(int y=0;y<HT;y++){
        A[P(u,v)]+=x_re[P(u,y)]*cos((2.0*PI*v*y)/HT)-x_im[P(u,y)]*sin((2.0*PI*v*y)/HT);
      }
    }
  }

  for(int x=0;x<WD;x++)
    for(int y=0;y<HT;y++)
      A[P(x,y)]/=(WD*HT);
}


double max_double(double *buff, int size){
  double ret=-1e30;
  int i;
  for(i=0;i<size;i++)ret=(ret < buff[i] ? buff[i] : ret);
  return ret;
}

double min_double(double *buff, int size){
  double ret=1e30;
  int i;
  for(i=0;i<size;i++)ret=(ret > buff[i] ? buff[i] : ret);
  return ret;
}

void flip(double *arr){
  int y;
  for(y=0; y<HT/2; y++){
    int x, src=y, dst=(HT-y-1);
    double tmp[WD];
    for(x=0;x<WD;x++)tmp[x]=arr[src*WD+x];
    for(x=0;x<WD;x++)arr[src*WD+x]=arr[dst*WD+x];
    for(x=0;x<WD;x++)arr[dst*WD+x]=tmp[x];
  }
}

int main(int argc, char **argv){
  double re[WD*HT], im[WD*HT], tempbuf[WD*HT], ABS[WD*HT];

  read_image("sample5.ras"); printf("read\n");

  for(int i=0;i<WD*HT;i++)tempbuf[i]=buf[i];

  flip(tempbuf);
  dft(tempbuf,re,im);
  idft(re,im,tempbuf);

  for(int i=0;i<WD*HT;i++)buf[i]=(unsigned char)tempbuf[i];
  
  /*// debug from here //
  for(i=0; i<WD*HT; i++){
    ABS[i] = log(sqrt(pow(re[i], 2) + pow(im[i], 2)));
  }

  for(i=0;i<WD*HT;i++){
    double max_ = max_double(ABS, WD*HT);
    double min_ = min_double(ABS, WD*HT);
    buf[i]=(unsigned char)((511*(ABS[i]-min_))/(max_-min_));
  }
  // debug upto here //*/

  write_image("output.ras"); printf("write\n");

  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_DOUBLE|GLUT_RGB|GLUT_DEPTH);
  glutInitWindowSize(WD, HT);
  glutInitWindowPosition(100, 100);
  glutCreateWindow(argv[0]);

  init();

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

  return 0;
}