/*---------------------------------------------------------------------------*/ /* fourier.c */ /*---------------------------------------------------------------------------*/ /*\ * Copyright (C) Feb, 2000 Oliver Knill, University of Texas in Austin * posted online December 3, 2010, with Harvard example picture added * for an exhibit for the linear algebra course Math21b * Oliver Knill, Harvard University * http://www.math.harvard.edu/~knill , knill@math.harvard.edu * * Compile: gcc -o fourier -lm fourier.c * Run: ./fourier 5 in.ppm >out.ppm * Purpose: makes cos transform and reverse transform on a ppm picture * this is for educational use. No fast cos fourier transform, just to illustrate * how to take the Fourier transform of a picture and then do the inverse transform. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. \*/ #include #include #include #include int x_size, y_size, c_size; // geometry and color of picture int X_size=4; // size of Fourier picture int Y_size=4; int f(char cen) { // char into integer 0-255 int a=cen; if (a<0) a+=256; a=a-127; return(a); } int g(char cen) { int a=cen; a=4*a+127; if (a<0) a=0; if (a>255) a=255; return(a); } char *filter(char *image_in){ int i, j; int m, n; char *ptr; char *ptr_out; char *image_out; float *fourier; float *ptr_four; float red,green,blue; float coeff; int area = x_size*y_size; float pi_x = 3.141592653/(x_size); float pi_y = 3.141592653/(y_size); if((image_out=(char *) malloc(3*x_size * y_size * sizeof(char))) == NULL){ fprintf(stderr,"Can't load image (memory allocation error)\n"); exit(1); /* image = { 0 <= m < x_size, -y_size < n <= 0} */ } if((fourier=(float *) malloc(3*X_size * Y_size * sizeof(float))) == NULL){ fprintf(stderr,"(Fourier picture memory allocation error)\n"); exit(1); /* image = { 0 <= i < X_size, -Y_size < j <= 0} */ } for(j = 0; j > -Y_size; j--) { // Fourier picture for(i = 0; i < X_size; i++) { red=0; green=0; blue=0; ptr_four = fourier + 3 * ((-j) * X_size + i ); for(n = 0; n > -y_size; n--) { for(m = 0; m < x_size; m++) { ptr = image_in + 3 * ((-n) * x_size + m ); coeff = cos(pi_x*i*m)*cos(pi_y*j*n); red += f(*ptr++)*coeff; green += f(*ptr++)*coeff; blue += f(*ptr )*coeff; } } red = red/area; green = green/area; blue = blue/area; *ptr_four++ = red; *ptr_four++ = green; *ptr_four = blue; } } for(n = 0; n > -y_size; n--) { // inverse transform for(m = 0; m < x_size; m++) { ptr_out= image_out + 3 * ((-n) * x_size + m ); red=0; green=0; blue=0; for(j = 0; j > -Y_size; j--) { for(i = 0; i < X_size; i++) { ptr_four = fourier + 3 * ((-j) * X_size + i ); coeff = cos(pi_x*i*m)*cos(pi_y*j*n); if ( (i==0) ) coeff=coeff/2; if ( (j==0) ) coeff=coeff/2; red += (*ptr_four++)*coeff; green += (*ptr_four++)*coeff; blue += (*ptr_four )*coeff; } } *ptr_out++ = g(red); *ptr_out++ = g(green); *ptr_out = g(blue); } } return(image_out); } /*---------------------------------------------------------------------------*/ /* Function write_ppm_image */ /*---------------------------------------------------------------------------*/ void write_ppm_image(char *image){ int i,j; char *ptr; fprintf(stdout,"P6\n%i %i\n%i\n", x_size, y_size, c_size); for(j = 0; j > -y_size; j--){ for(i = 0 ; i < x_size; i++){ ptr = image + 3 * ((-j) * x_size + i ); fputc(*ptr++,stdout); fputc(*ptr++,stdout); fputc(*ptr ,stdout); } } fclose(stdout); } /*---------------------------------------------------------------------------*/ /* Function read_ppm_image */ /*---------------------------------------------------------------------------*/ char *read_ppm_image(char *filename) { FILE *in; int i; char buffer[1025]; char *image,*ptr; if((in = fopen(filename,"rb")) == NULL){ fprintf(stderr,"Can't open file %s\n", filename); exit(1); } fgets(buffer,1025,in); if(strncmp(buffer,"P6",2)){ fprintf(stderr,"Unsupported file format (need PPM raw)\n"); exit(1); } do fgets(buffer,1025,in); while(*buffer == '#'); x_size = atoi(strtok(buffer," ")); y_size = atoi(strtok(NULL," ")); fgets(buffer,1025,in); c_size = atoi(buffer); if((image = (char *) malloc( 3 * x_size * y_size * sizeof(char))) == NULL){ fprintf(stderr,"Can't load image (memory allocation error)\n"); exit(1); } i = 0; ptr = image; while(!feof(in) && i < 3 * x_size * y_size){ *ptr++ = fgetc(in); i++; } fclose(in); if(i < x_size * y_size){ fprintf(stderr,"Premature end image\n"); exit(1); } return(image); } /*---------------------------------------------------------------------------*/ /* Function main */ /*---------------------------------------------------------------------------*/ int main(int argc, char *argv[]){ char *image; char *filename; int M; int i; if(argc < 1){ fprintf(stderr,"Usage: filter [N] [pnmfile] \n"); exit(1); } if(argc > 3){ fprintf(stderr,"Usage: filter [N] [pnmfile] \n"); exit(1); } if(argc == 3) { M=atoi(argv[1]); filename=argv[2]; } else { filename=argv[1]; } image=read_ppm_image(filename); if(argc == 3) {X_size=M; Y_size=M; }; // else {X_size=x_size; Y_size=y_size;} write_ppm_image(filter(image)); return 0; }