#include <stdio.h>
#include <math.h>

#define CHUNKSIZELOG 9
#define CHUNKSIZE (1<<CHUNKSIZELOG)

typedef double ext;
  /* I used to use long double format, but memory
   * consumption is a bit excessive */
typedef struct { ext r; ext i; } cmp;
typedef cmp CHUNK[CHUNKSIZE];
typedef CHUNK IMAGE[CHUNKSIZE];
  /* All my images are square, that avoids my confusing
   * absissa and ordinate :-) */

ext invsqrt2;
CHUNK dtabchunk,rtabchunk,inchunk,outchunk;
IMAGE image,inter,spect;

typedef unsigned char PPMIMAGE[CHUNKSIZE][CHUNKSIZE][3];

PPMIMAGE inimage,outimage;

void preparefft(void)
{
  int i;
  invsqrt2 = M_SQRT1_2;
  for (i=0;i<CHUNKSIZE;i++) {
    dtabchunk[i].r = cos(-2.0L*M_PI*((ext)i)/((ext)CHUNKSIZE));
    dtabchunk[i].i = sin(-2.0L*M_PI*((ext)i)/((ext)CHUNKSIZE));
  }
  for (i=0;i<CHUNKSIZE;i++) {
    rtabchunk[i].r = cos(2.0L*M_PI*((ext)i)/((ext)CHUNKSIZE));
    rtabchunk[i].i = sin(2.0L*M_PI*((ext)i)/((ext)CHUNKSIZE));
  }
}

void fft(int loglen,cmp *input,cmp *output,cmp *table)
{
  int i,l,l2;
  CHUNK nextinput,nextoutput;
  if (loglen==0) {
    output[0] = input[0];
    return;
  }
  if (loglen==1) {
    output[0].r = invsqrt2*(input[0].r+input[1].r);
    output[0].i = invsqrt2*(input[0].i+input[1].i);
    output[1].r = invsqrt2*(input[0].r-input[1].r);
    output[1].i = invsqrt2*(input[0].i-input[1].i);
    return;
  }
  l=(1<<loglen);
  l2=l/2;
  for (i=0;i<l2;i++) {
    nextinput[i] = input[i*2];
    nextinput[i+l2] = input[i*2+1];
  }
  fft(loglen-1,nextinput,nextoutput,table);
  fft(loglen-1,nextinput+l2,nextoutput+l2,table);
  for (i=0;i<l;i++) {
    nextoutput[i].r *= invsqrt2;
    nextoutput[i].i *= invsqrt2;
  }
  for (i=0;i<l;i++) {
    output[i].r = nextoutput[i%l2].r
                + table[i<<(CHUNKSIZELOG-loglen)].r*nextoutput[(i%l2)+l2].r
                - table[i<<(CHUNKSIZELOG-loglen)].i*nextoutput[(i%l2)+l2].i;
    output[i].i = nextoutput[i%l2].i
                + table[i<<(CHUNKSIZELOG-loglen)].i*nextoutput[(i%l2)+l2].r
                + table[i<<(CHUNKSIZELOG-loglen)].r*nextoutput[(i%l2)+l2].i;
  }
}

void directfft(void)
{
  int i,j;
  fprintf(stderr,"Direct Fourier transform.\n");
  fprintf(stderr,"Horizontal part.\n");
  for (i=0;i<CHUNKSIZE;i++) {
    fft(CHUNKSIZELOG,image[i],inter[i],dtabchunk);
  }
  fprintf(stderr,"Vertical part.\n");
  for (i=0;i<CHUNKSIZE;i++) {
    for (j=0;j<CHUNKSIZE;j++) inchunk[j]=inter[j][i];
    fft(CHUNKSIZELOG,inchunk,outchunk,dtabchunk);
    for (j=0;j<CHUNKSIZE;j++) spect[j][i]=outchunk[j];
  }
}

void reversefft(void)
{
  int i,j;
  fprintf(stderr,"Reverse Fourier transform.\n");
  fprintf(stderr,"Horizontal part.\n");
  for (i=0;i<CHUNKSIZE;i++) {
    fft(CHUNKSIZELOG,spect[i],inter[i],rtabchunk);
  }
  fprintf(stderr,"Vertical part.\n");
  for (i=0;i<CHUNKSIZE;i++) {
    for (j=0;j<CHUNKSIZE;j++) inchunk[j]=inter[j][i];
    fft(CHUNKSIZELOG,inchunk,outchunk,rtabchunk);
    for (j=0;j<CHUNKSIZE;j++) image[j][i]=outchunk[j];
  }
}

void processfft(void)
{
  int i,j,ii,jj;
  ext u,v,t1,t2,t3,t4,t5,t6,k,Q;
  for (i=0;i<CHUNKSIZE;i++) for (j=0;j<CHUNKSIZE;j++) {
    inter[i][j]=spect[i][j];
  }
#if 0
  /* This one can be used to recognize large shapes */
  for (i=0;i<CHUNKSIZE;i++) for (j=0;j<CHUNKSIZE;j++) {
    u=((ext)i)/(ext)CHUNKSIZE; v=((ext)j)/(ext)CHUNKSIZE;
    if (u>0.5L) u-=1.0L; if (v>0.5L) v-=1.0L;
    t1=u*u+v*v;
    t2=25.0L*sqrt(t1);
    if (t2>1.0E-6L) t3=sin(t2)/t2; else t3=1.0L;
    t4=exp(t1*24.0L);
    k=t3*t4;
    spect[i][j].r=inter[i][j].r*k;
    spect[i][j].i=inter[i][j].i*k;
  }
#endif
#if 0
  /* This is the resistance in an RCL circuit */
  for (i=0;i<CHUNKSIZE;i++) for (j=0;j<CHUNKSIZE;j++) {
    u=((ext)i)/(ext)CHUNKSIZE; v=((ext)j)/(ext)CHUNKSIZE;
    if (u>0.5L) u-=1.0L; if (v>0.5L) v-=1.0L;
    t1=10.0L*hypot(u,v);
    k=1.0L;
    Q=1.0L;
    t2=1-t1*t1/(k*k); t3=t1/(Q*k);
    t4=t2*t2+t3*t3;
    t5=-t3*t3/t4; t6=t2*t3/t4;
    spect[i][j].r=inter[i][j].r*t5-inter[i][j].i*t6;
    spect[i][j].i=inter[i][j].i*t5+inter[i][j].r*t6;
  }
#endif
  for (i=0;i<CHUNKSIZE;i++) for (j=0;j<CHUNKSIZE;j++) {
    ii=i^419; jj=j^206;
    spect[i][j]=inter[ii][jj];
  }
}

void getaplane(int plane)
{
  int i,j;
  for (i=0;i<CHUNKSIZE;i++) for (j=0;j<CHUNKSIZE;j++) {
    image[i][j].r = ((ext)inimage[i][j][plane]-127.5L)/127.5L;
    image[i][j].i = 0.0L;
  }
}

void processplane(void)
{
  int i;
  directfft();
  processfft();
  reversefft();
}

void outputplane(int plane)
{
  int i,j;
  ext d;
  for (i=0;i<CHUNKSIZE;i++) for (j=0;j<CHUNKSIZE;j++) {
    d = image[i][j].r;
    if (d<-1.0L) d=-1.0L;
    if (d>1.0L) d=1.0L;
    outimage[i][j][plane] = (unsigned char)((d*127.5L)+127.5L);
  }
}

void readimage(void)
{
  if (fread(inimage,CHUNKSIZE*CHUNKSIZE*3,1,stdin)<1) {
    fprintf(stderr,"Error reading input file.\n");
    exit(1);
  }
}

void writeimage(void)
{
  fwrite(outimage,CHUNKSIZE*CHUNKSIZE*3,1,stdout);
}

void getheader(void)
{
  char s[301];
  int i,j;
  if (strcmp(fgets(s,300,stdin),"P6\n")) {
    fprintf(stderr,"Sorry - wrong image format.\n");
    exit(1);
  }
  while (fgets(s,300,stdin)[0]=='#');
  if (sscanf(s,"%d %d",&i,&j)<2) {
    fprintf(stderr,"Error parsing format line.\n");
    exit(1);
  }
  if ((i!=CHUNKSIZE)||(j!=CHUNKSIZE)) {
    fprintf(stderr,"Sorry - I insist on my images being %dx%d.\n",CHUNKSIZE,CHUNKSIZE);
    exit(1);
  }
  fgets(s,300,stdin);
  if (sscanf(s,"%d",&i)<1) {
    fprintf(stderr,"Error parsing color range line.\n");
    exit(1);
  }
  if (i!=255) {
    fprintf(stderr,"Sorry - I insist on my images having 255 color range.\n");
    exit(1);
  }
  printf("P6\n%d %d\n255\n",CHUNKSIZE,CHUNKSIZE);
}

void main(void)
{
  int plane;
  preparefft();
  getheader();
  readimage();
  for (plane=0;plane<3;plane++) {
    fprintf(stderr,"Processing plane %d.\n",plane);
    getaplane(plane);
    processplane();
    outputplane(plane);
  }
  writeimage();
  fflush(stdout);
}