/* make_fv_dir_all.c */
/* 各数字の方向分布特徴ベクトルの生成 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mypgm.h"

#define BLACK 0
#define WHITE 255
#define RADIUS 25.0
#define HEI 120
#define WID 80
#define SIDE 20
#define BLOCKS 24
#define DIRECTION 4

/* prototype declaration */
void bilinear_normal_image(double);
int dircode(int, int);

/* global variables */
double fv[BLOCKS][DIRECTION];


main( )
{
  FILE *fp;
  char filename[MAX_FILENAME];
  
  static int l_pat[10] = {2535, 1691, 1658, 1736, 1188,
			  1996, 1977, 1410, 2366, 1428};

  static int t_pat[10] = {1485, 2392, 1939, 2043, 2403,
			  1641, 1831, 2056, 986, 1140}; 
  
  int digit, serial;
  int code;
  int i, j, k;
  int sample = 0;
  
  
  /* fv generation of learning data */
  for (digit = 0; digit < 10; digit++) {
    for (serial = 1; serial <= l_pat[digit]; serial++) {
      sprintf(filename, "./cdrom1b/ldg%d_%04d.pgm", digit, serial);
      
      /* input of image data */ 
      load_image_file(filename);
      bilinear_normal_image(RADIUS);
      printf("sample  = %d\n", ++sample);
      
      for (j = 0; j < y_size2; j++) {
	for (i = 0; i < x_size2; i++) {
	  if (image2[j][i] == BLACK)
	    image1[j][i] = 1;
	  else
	    image1[j][i] = 0;
	}
      }
      
      for (i = 0; i < BLOCKS; i++)
	for (j = 0; j < DIRECTION; j++)
	  fv[i][j] = 0.0;

      for (j = 0; j < y_size2 - 1; j++) {
	for (i = 0; i < x_size2 - 1; i++) {
	  code = dircode(j, i);
	  if (code != -1) {
	    k = j/SIDE*(WID/SIDE) + i/SIDE;
	    fv[k][code]++;
	  }
	}
      }

      for (i = 0; i < BLOCKS; i++)
	for (j = 0; j < DIRECTION; j++)
	  fv[i][j] /= (double)SIDE * SIDE;

      /* output fv */
      sprintf(filename, "./fv_dir/lfv%d_%04d", digit, serial);
      fp = fopen(filename, "wb");
      fwrite(fv, sizeof(fv), 1, fp);
      fclose(fp);

    }
  }

  /* fv generation of test data */
  for (digit = 0; digit < 10; digit++) {
    for (serial = 1; serial <= t_pat[digit]; serial++) {
      sprintf(filename, "./cdrom1b/tdg%d_%04d.pgm", digit, serial);
      
      /* input of image data */ 
      load_image_file(filename);
      bilinear_normal_image(RADIUS);
      printf("sample  = %d\n", ++sample);
      
      for (j = 0; j < y_size2; j++) {
	for (i = 0; i < x_size2; i++) {
	  if (image2[j][i] == BLACK)
	    image1[j][i] = 1;
	  else
	    image1[j][i] = 0;
	}
      }
      
      for (i = 0; i < BLOCKS; i++)
	for (j = 0; j < DIRECTION; j++)
	  fv[i][j] = 0.0;

      for (j = 0; j < y_size2 - 1; j++) {
	for (i = 0; i < x_size2 - 1; i++) {
	  code = dircode(j, i);
	  if (code != -1) {
	    k = j/SIDE*(WID/SIDE) + i/SIDE;
	    fv[k][code]++;
	  }
	}
      }

      for (i = 0; i < BLOCKS; i++)
	for (j = 0; j < DIRECTION; j++)
	  fv[i][j] /= (double)SIDE * SIDE;

      /* output fv */
      sprintf(filename, "./fv_dir/tfv%d_%04d", digit, serial);
      fp = fopen(filename, "wb");
      fwrite(fv, sizeof(fv), 1, fp);
      fclose(fp);

    }
  }
}


void bilinear_normal_image(double norm_radius)
{
  int x, y;
  int count;
  double gx, gy;
  double rad, ratio;
  double x_new, y_new, x_frac, y_frac;
  double gray_new;
  int gv_new;
  int m, n;
  
  gx = gy = 0.0;
  count = 0;
  for (y = 0; y < y_size1; y++) {
    for (x = 0; x < x_size1; x++) {
      if (image1[y][x] == BLACK) {
	gx += (double)x;
	gy += (double)y;
	count++;
      }
    }
  }
  gx /= (double)count;
  gy /= (double)count;
  rad = 0.0;
  for (y = 0; y < y_size1; y++) {
    for (x = 0; x < x_size1; x++) {
      if (image1[y][x] == BLACK) {
	rad += sqrt(pow((double)x - gx, 2) + pow((double)y - gy, 2));
      }
    }
  }
  rad /= (double)count;
  if (rad == 0.0) exit(1);
  ratio = norm_radius / rad;
  
  x_size2 = x_size1;
  y_size2 = y_size1;
  for (y = 0; y < y_size2; y++) {
    for (x = 0; x < x_size2; x++) {
      x_new = (x - x_size2/2) / ratio + gx;
      y_new = (y - y_size2/2) / ratio + gy;
      m = (int)floor(x_new);
      n = (int)floor(y_new);
      x_frac = x_new - m;
      y_frac = y_new - n;
      
      if (m >= 0 && m+1 < x_size2 && n >= 0 && n+1 < y_size2) {
	gray_new = (1.0 - y_frac)*((1.0 - x_frac)*image1[n][m] +
				   x_frac *image1[n][m+1]) +
	  y_frac *((1.0 - x_frac)*image1[n+1][m] +
		   x_frac *image1[n+1][m+1]);
	gv_new = gray_new >= MAX_BRIGHTNESS/2.0 ? MAX_BRIGHTNESS : 0;
      } else {
	gv_new = MAX_BRIGHTNESS;
      }
      image2[y][x] = (unsigned char)gv_new;
    }
  }
}


int dircode(int y, int x)
{
  int peri[4], sum;

  peri[0] = image1[y][x];
  peri[1] = image1[y][x+1];
  peri[2] = image1[y+1][x];
  peri[3] = image1[y+1][x+1];
  sum = peri[0] + peri[1] + peri[2] + peri[3];
  if (sum == 2) {
    if (peri[0] * peri[1] == 1 || peri[2] * peri[3] == 1) {
      return 0;  /* Horizontal */
    }
    if (peri[0] * peri[2] == 1 || peri[1] * peri[3] == 1) {
      return 2;  /* Vertical */
    }
	return -1;
  } else if (sum == 1) {
    if (peri[0] == 1 || peri[3] == 1) {
      return 1;  /* Right-diagonal */
    }
    if (peri[1] == 1 || peri[2] == 1) {
      return 3;  /* Left-diagonal */
    }
	return -1;
  } else {
    return -1;
  }
}