Skip to content

Instantly share code, notes, and snippets.

@domitry

domitry/acepirin.c

Last active Nov 3, 2015
Embed
What would you like to do?
#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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment