/* FFT2fltr_powB.c */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mypgm.h"
#include "FFT1.h"
#include "FFT2.h"

void bandpass_powB_filtering( )
{
  int i, j;
  int sum;
  double rad, rad_L, rad_H, rad_G;
  double norm, value, mean, var, sdv;
  double th_pow;
  
  x_size2 = num_of_data;
  y_size2 = num_of_data;
  
  printf("\nFFT image is generated.\n\n");
  
  /* FFT bandpass filtering */
  rad_G = 8.0;
  rad_L = 24.0;
  rad_H = 80.0;
  
  /* mean power calculation */
  mean = 0.0;
  var = 0.0;
  sum = 0;
  for (j = 0; j < num_of_data; j++) {
    for (i = 0; i < num_of_data; i++) {
      rad = sqrt(pow(j - num_of_data/2, 2) + pow(i - num_of_data/2, 2));
      if (rad <= rad_H && (rad >= rad_L || rad <= rad_G)) {
	value = pow(data[j][i], 2) + pow(jdata[j][i], 2);
	norm = log(1.0 + value);
	mean += norm;
	var += norm * norm;
	sum++;
      }
    }
  }
  mean /= (double)sum;
  var = (var - mean * mean * (double)sum) / (double)(sum - 1);
  sdv = sqrt(var);
  th_pow = mean;
  
  /* frequency and power filtering */
  for (j = 0; j < num_of_data; j++) {
    for (i = 0; i < num_of_data; i++) {
      rad = sqrt(pow(j - num_of_data/2, 2) + pow(i - num_of_data/2, 2));
      value = pow(data[j][i], 2) + pow(jdata[j][i], 2);
      norm = log(1.0 + value);
      if (rad > rad_H || (rad < rad_L && rad > rad_G)
	  || norm <= th_pow) {
	data[j][i] = 0.0;
	jdata[j][i] = 0.0;
      } else {
	data[j][i] *= norm;
	jdata[j][i] *= norm;
      }
    }
  }
}


main( )
{
  int i, j;
  int minsize, cutsize, x_offset, y_offset;
  
  load_image_data( );      /* Input image load      */
  
  cutsize = 1;
  minsize = x_size1 > y_size1 ? y_size1 : x_size1;
  while (1) {
    if (2*cutsize > minsize) break;
    cutsize *= 2;
  }
  
  printf("cut size = %d\n", cutsize);
  
  x_offset = (x_size1 - cutsize)/2;
  y_offset = (y_size1 - cutsize)/2;
  
  /* image size change to the power of two */
  x_size1 = y_size1 = cutsize;
  
  for (j = 0; j < y_size1; j++) {
    for (i = 0; i < x_size1; i++) {
      image2[j][i] = image1[j+y_offset][i+x_offset];
    }
  }
  for (j = 0; j < y_size1; j++) {
    for (i = 0; i < x_size1; i++) {
      image1[j][i] = image2[j][i];
    }
  }
  
  make_original_data( );   /* FFT arrays generation   */
  FFT2(1);                 /* 2-d FFT  */
  origin_centered_FFT();   /* FFT rearrangement */
  bandpass_powB_filtering( );   /* bandpass and power filtering */
  origin_centered_FFT();   /* FFT rearrangement */
  FFT2(-1);                /* 2-d inverse FFT */
  make_output_image_enhanced( );    /* Image generation after IFFT */
  save_image_data( );      /* Output image save */
  return 0;
}