extract_features.c File Reference

#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <math.h>
#include <string.h>
#include <unistd.h>
#include <ppm.h>
#include "gift_features.h"
#include <stdint.h>

Go to the source code of this file.

Defines
#define	numH   18
#define	numS   3
#define	numV   3
#define	numGrey   4
#define	OLDFIXED
#define	num_gabor_scales   3
#define	num_total_colour_blocks   (256+64+16+4)
Functions
void	init_feature_variables (uint32_t, uint32_t **, uint32_t ***, uint32_t ***)
void	extract_gabor_features (PPM *, uint32_t ***, uint32_t ***)
void	extract_mode_features (PPM *, uint32_t *, uint32_t, uint32_t **, byte *, uint32_t *)
enum ppm_error	write_mode_features (char *, uint32_t, uint32_t *, uint32_t ***, uint32_t ***, byte *)
int	main (int argc, char *argv[])

---

Define Documentation

#define num_gabor_scales   3

Definition at line 28 of file extract_features.c.

#define num_total_colour_blocks   (256+64+16+4)

Definition at line 30 of file extract_features.c.

#define numGrey   4

Definition at line 23 of file extract_features.c.

Referenced by main().

#define numH   18

Definition at line 20 of file extract_features.c.

Referenced by main().

#define numS   3

Definition at line 21 of file extract_features.c.

Referenced by main().

#define numV   3

Definition at line 22 of file extract_features.c.

Referenced by main().

#define OLDFIXED

Definition at line 25 of file extract_features.c.

---

Function Documentation

void extract_gabor_features	(	PPM *	,
		uint32_t ***	,
		uint32_t ***
	)

Definition at line 126 of file extract_block_features.c.

References create_filter_kernels(), gabor_block_size, gabor_filter(), gabor_ranges, _PPM::height, MAX_HEIGHT, MAX_WIDTH, num_gabor_ranges, num_gabor_scales, num_gabors_per_scale, square, _PPM::value_plane_double_reversed, and _PPM::width.

Referenced by main().

                                                                                                       {

  int i, j, x, y, k;
  int scale, orientation;
  int energy_class;
  PPM *value_plane;
  double *value_image_dbl, *filtered_image;
  double *conv;
  double *conv2;
  double gabor_mean;
  double * kernelsxy[num_gabor_scales*num_gabors_per_scale];

  /* for Gabor features, this is what I'm going to do:
    - Apply each of the 12 filters (3 scales, 4 orientations) to each
      pixel in the image.
    - Use 16x16 blocks only.
    - In each block find the rms energy of each filter.
    - quantize these into num_gabor_ranges levels, as specified in the
      "array" gabor_ranges[]
  */

  /* this is all created with one malloc, so that it can be cleared in one memset() call */
  filtered_image = (double *)malloc(im_hsv->width*im_hsv->height*sizeof(double)*3);
  conv=&filtered_image[im_hsv->width*im_hsv->height];
  conv2=&conv[im_hsv->width*im_hsv->height];

  create_filter_kernels(kernelsxy);

  /* apply each filter to the image */
  for (scale = 0; scale < num_gabor_scales; scale++) {
    for (orientation = 0; orientation < num_gabors_per_scale; orientation++) {

      /* this clears conv, conv2, and filtered_image in one call. */
      memset(filtered_image, 0, MAX_WIDTH*MAX_HEIGHT*sizeof(double)*3);

      /* filter the image */
      gabor_filter(im_hsv->value_plane_double_reversed, im_hsv->width, im_hsv->height, scale, orientation, kernelsxy, conv, conv2, filtered_image);

      /* extract the rms energy for each block */
      k = 0; /* block counter */
      for (y = 0; y < im_hsv->height; y += gabor_block_size) {
      for (x = 0; x < im_hsv->width; x += gabor_block_size) {
        gabor_mean = 0;
        for (i = 0; (i < gabor_block_size) &&  (y+i < im_hsv->height); i++) {
        for (j = 0; (j < gabor_block_size) &&  (x+j < im_hsv->height); j++) {
          gabor_mean += square((filtered_image[(y + i)*im_hsv->width + (x + j)]));
        }
        }
        gabor_mean /= square(gabor_block_size);
        gabor_mean = sqrt(gabor_mean);

        /* find the energy class for this block */
        for (energy_class = 0; energy_class < num_gabor_ranges; energy_class++) {
          if (gabor_mean < gabor_ranges(energy_class))
            break;
        }

        /* store this class in the appropriate place */
        block_gabor_class[scale][orientation][k] = energy_class;

        /* increment the histogram entry for this combination */
        gabor_histogram[scale][orientation][energy_class]++;
        
        k++; /* increment block counter */
      }
      }
    }
  }
}

void extract_mode_features	(	PPM *	,
		uint32_t *	,
		uint32_t	,
		uint32_t **	,
		byte *	,
		uint32_t *
	)

Definition at line 196 of file extract_block_features.c.

References _PPM::bytes_per_pixel, colmap2rgb_ppm(), destroy_ppm(), _PPM::height, hsv2rgb_ppm(), image_size, _PPM::max_col_comp, new_ppm(), num_colour_scales, PGM_RAW, _PPM::pixel, ppm_handle_error(), PPM_OK, PPM_RAW, smallest_colour_block, square, _PPM::type, _PPM::width, and write_ppm().

Referenced by main().

                                                                                                                                                 {

  int i, j, k, last_k, k1, m, n, r, x, y;
  byte colour;
  int scale, block_size, num_blocks, old_num_blocks;
  int max_count, mode_index;
  int b1, b2, b3, b4;
#ifdef GENERATE_BLOCK_IMAGES
  PPM *out_image, *out_image_hsv, *out_image_rgb;
  enum ppm_error the_error;
  FILE *outfile;
  char out_fname[256];
#endif

#ifdef GENERATE_BLOCK_IMAGES
    /* make space for an HSV image of the same size as that read */
    out_image = new_ppm();
    out_image->type = PGM_RAW;
    out_image->width = im->width;
    out_image->height = im->height;
    out_image->max_col_comp = im->max_col_comp;
    out_image->bytes_per_pixel = 1;
    out_image->pixel = (byte *)malloc(out_image->width*out_image->height*sizeof(byte));
#endif

  /* first extract the lowest level blocks, directly from the pixels of
  the PGM image */
  scale = 0;
  block_size = smallest_colour_block;
  num_blocks = square(image_size/block_size);
  k = last_k = 0;
  for (y = 0; y < im->height; y += block_size) {
  for (x = 0; x < im->width; x += block_size) {
    for (n = 0; n < colmap_size; n++) {
      col_counts[k][n] = 0; 
    }
    for (i = 0; (i < block_size) && (y+i < im->height); i++) {
      for (j = 0; (j < block_size) && (x+j < im->width); j++) {
        colour = im->pixel[(y + j)*im->width + (x + i)];
        col_counts[k][colour]++;
      }
    }

    /* find the mode */
    max_count = mode_index = 0;
    for (n = 0; n < colmap_size; n++) {
      if (col_counts[k][n] > max_count) {
        max_count = col_counts[k][n];
        mode_index = n;
      }
    }
    block_mode[k] = mode_index;
    k++;
  }
  }

#ifdef GENERATE_BLOCK_IMAGES
  /* generate an image with this */ 
  k1 = 0;
  for (y = 0; y < im->height; y += block_size) {
  for (x = 0; x < im->width; x += block_size) {
    for (i = 0; i < block_size; i++) {
    for (j = 0; j < block_size; j++) {
      out_image->pixel[(y + j)*im->width + (x + i)] = block_mode[k1];
    }
    }
    k1++;
  }
  }
  /* convert from the colour map to an HSV image */
  if ((the_error = colmap2rgb_ppm(out_image, colmap, colmap_size, &out_image_hsv)) != PPM_OK) {
    ppm_handle_error(the_error);
    exit(1);
  }
  /* convert to RGB */
  if ((the_error = hsv2rgb_ppm(out_image_hsv, &out_image_rgb)) != PPM_OK) {
    ppm_handle_error(the_error);
    exit(1);
  }
  /* write it */
  sprintf(out_fname, "blocks_%dx%d.ppm", block_size, block_size);
  outfile = fopen(out_fname, "wb");
  if ((the_error = write_ppm(outfile, out_image_rgb, PPM_RAW)) != PPM_OK) {
    ppm_handle_error(the_error);
    exit(1);
  }
  fclose(outfile);
  destroy_ppm(&out_image_rgb);
  destroy_ppm(&out_image_hsv);
#endif


  /* now do the other scales */
  for (scale = 1; scale < num_colour_scales; scale++) {
    block_size *= 2;
    num_blocks /= 4;
    r = (int)sqrt((double)num_blocks);
    for (i = 0; i < num_blocks; i++) {
      m = (i/r)*2*r; /* note: this first is *integer* division */
      b1 = last_k + 2*i + m;
      b2 = last_k + 2*i + m + 1;
      b3 = last_k + 2*i + m + 2*r;
      b4 = last_k + 2*i + m + 2*r + 1;
      max_count = mode_index = 0;
      for (n = 0; n < colmap_size; n++) {
        col_counts[k + i][n] = 
          col_counts[b1][n] + col_counts[b2][n] + 
          col_counts[b3][n] + col_counts[b4][n];
        if (col_counts[k + i][n] > max_count) {
          max_count = col_counts[k + i][n];
          mode_index = n;
        }
      }
      block_mode[k + i] = mode_index;
    }
#ifdef GENERATE_BLOCK_IMAGES
    /* generate an image with this */ 
    for (y = 0; y < im->height; y += block_size) {
    for (x = 0; x < im->width; x += block_size) {
      for (i = 0; i < block_size; i++) {
      for (j = 0; j < block_size; j++) {
        out_image->pixel[(y + j)*im->width + (x + i)] = block_mode[k1];
      }
      }
      k1++;
    }
    }
    /* convert from the colour map to an HSV image */
    if ((the_error = colmap2rgb_ppm(out_image, colmap, colmap_size, &out_image_hsv)) != PPM_OK) {
      ppm_handle_error(the_error);
      exit(1);
    }
    /* convert to RGB */
    if ((the_error = hsv2rgb_ppm(out_image_hsv, &out_image_rgb)) != PPM_OK) {
      ppm_handle_error(the_error);
      exit(1);
    }
    /* write it */
    sprintf(out_fname, "blocks_%dx%d.ppm", block_size, block_size);
    outfile = fopen(out_fname, "wb");
    if ((the_error = write_ppm(outfile, out_image_rgb, PPM_RAW)) != PPM_OK) {
      ppm_handle_error(the_error);
      exit(1);
    }
    fclose(outfile);
    destroy_ppm(&out_image_rgb);
    destroy_ppm(&out_image_hsv);
#endif
    last_k = k;
    k += num_blocks;
  }

  /* now finally consolidate all the colour counts to produce the
  histogram for the whole image */
  k -= num_blocks;
  for (n = 0; n < colmap_size; n++)
    col_histogram[n] = 
      col_counts[k][n] + col_counts[k + 1][n] + 
      col_counts[k + 2][n] + col_counts[k + 3][n];
}

void init_feature_variables	(	uint32_t	,
		uint32_t **	,
		uint32_t ***	,
		uint32_t ***
	)

Definition at line 105 of file extract_block_features.c.

References gabor_block_size, image_size, num_blocks_at_scale, num_colour_scales, num_gabor_ranges, num_gabor_scales, num_gabors_per_scale, and square.

Referenced by main().

                                                                                                                                        {

  uint32_t i, j, k=0;
  /* colour features */
  for (i = 0; i < num_colour_scales; i++) {
    for (j = 0; j < num_blocks_at_scale(i); j++) {
      col_counts[k++] = (uint32_t *)malloc(colmap_size*sizeof(uint32_t));
    }
  }

  /* Gabor features */
  for (i = 0; i < num_gabor_scales; i++) {
    block_gabor_class[i] = (uint32_t **)malloc(num_gabors_per_scale*sizeof(uint32_t *));
    gabor_histogram[i] = (uint32_t **)malloc(num_gabors_per_scale*sizeof(uint32_t *));
    for (j = 0; j < num_gabors_per_scale; j++) {
      block_gabor_class[i][j] = (uint32_t *)malloc(square(image_size/gabor_block_size)*sizeof(uint32_t));
      gabor_histogram[i][j] = (uint32_t *)calloc(num_gabor_ranges, sizeof(uint32_t));
    }
  }
}

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 38 of file extract_features.c.

References destroy_ppm(), extract_gabor_features(), extract_mode_features(), hsv_quantize_ppm(), init_feature_variables(), num_gabor_scales, num_total_colour_blocks, numGrey, numH, numS, numV, PGM_ASC, PGM_RAW, PPM_ASC, ppm_handle_error(), PPM_OK, PPM_RAW, read_magic_no(), read_ppm(), rgb2hsv_ppm(), and write_mode_features().

                                 {

  char *in_fname, *out_fname;
  char *point_pos;
  FILE *ppm_file;
  PPM *im_rgb, *im_hsv, *im_quant;
  uint32_t *colmap, colmap_size;
  enum file_types ppm_type;
  enum ppm_error the_error;
  uint32_t ** block_gabor_class[num_gabor_scales];
  uint32_t ** gabor_histogram[num_gabor_scales];
  /* the "mode" of each color block. */
  byte block_mode[num_total_colour_blocks];
  /* the colour count of each block. for quantizing the image */
  uint32_t * col_counts[num_total_colour_blocks];
  uint32_t * col_histogram;
 
  switch(argc) {
  case 2:
    in_fname = argv[1];
    break;
  default:
    fprintf(stderr, "Usage: %s ppm_file \n\n", argv[0]);
    exit(1);
    break;
  }

  if ((ppm_file = fopen (in_fname, "r")) == NULL) {
    fprintf(stderr, "Can't open file: %s", in_fname);
    exit(1);
  }

  /* now get the filename prefix */
  out_fname = (char *)malloc((strlen(in_fname) + 10)*sizeof(char));
  if ((point_pos = strchr(in_fname, '.')) == NULL) {
    fprintf(stderr, "File %s has no ""."" - can't generate features filename\n\n", in_fname);
    exit(1);
  }
  else {
    /* must find the *last* "." in the filename */
    point_pos = &in_fname[strlen(in_fname)];
    while (*point_pos != '.') {
      point_pos--;
    }
    strncpy(out_fname, in_fname, (int)(point_pos - in_fname));
    strcat(out_fname, ".fts");
  }
#ifdef DEBUG
  fprintf(stderr, "Features will be written to file %s\n", out_fname);
#endif
  

  /* read the rgb image from we are going to extract features */
  switch(ppm_type = read_magic_no(ppm_file)) {
  case PGM_ASC: case PPM_ASC: case PGM_RAW: case PPM_RAW:
    if ((the_error = read_ppm(ppm_file, &im_rgb, ppm_type)) != PPM_OK) {
      ppm_handle_error(the_error);
      exit(1);
    }
    break;
  default:
    fprintf(stderr, "Unrecognized file type.\n");
    exit(1);
    break;
  }

  /* convert it to hsv */
  if ((the_error = rgb2hsv_ppm(im_rgb, &im_hsv)) != PPM_OK) {
    ppm_handle_error(the_error);
    exit(1);
  }

  /* quantize the image */
  colmap_size = numH*numS*numV + numGrey;
  if ((the_error = hsv_quantize_ppm(im_hsv, &im_quant, &colmap, numH, numS, numV, numGrey)) != PPM_OK) {
    ppm_handle_error(the_error);
    exit(1);
  }

  /* FIXME: move this back into init_feature_variables. */
  col_histogram = (uint32_t *)malloc(colmap_size*sizeof(uint32_t));

  /* initialise the variables required for the feature extraction */
  init_feature_variables(colmap_size, col_counts, block_gabor_class, gabor_histogram);

  /* extract the features */

  extract_mode_features(im_quant, colmap, colmap_size, col_counts, block_mode, col_histogram);

  extract_gabor_features(im_hsv, block_gabor_class, gabor_histogram);

  /* write them to the file */
  if ((the_error = write_mode_features(out_fname, colmap_size, col_histogram, block_gabor_class, gabor_histogram, block_mode)) != PPM_OK) {
    ppm_handle_error(the_error);
    exit(1);
  }
  
  /* everything is OK */
  destroy_ppm(&im_hsv);
  destroy_ppm(&im_quant);
  free(colmap);
  exit(0);
}

enum ppm_error write_mode_features	(	char *	,
		uint32_t	,
		uint32_t *	,
		uint32_t ***	,
		uint32_t ***	,
		byte *
	)

HACK: convert all this features to floats.

WARNING: we write here floats into an integer field of a structure

Why is this necessary? The original version rounds at many places So we cannot just set FREQ_MAX=1. (WM)

Definition at line 357 of file extract_block_features.c.

References FILE_CLOSE_ERROR, FILE_OPEN_ERROR, FILE_WRITE_ERROR, FREQ_MAX, _FEATURE_DATA::frequency, gabor_block_size, _FEATURE_DATA::id, image_size, num_gabor_ranges, num_gabor_scales, num_gabors_per_scale, num_total_colour_blocks, PPM_OK, and square.

Referenced by main().

                                                                                                                                                                                     {

  FILE *out_file = NULL;
  int feature_index;
  int block_features_offset;
  int scale, orientation, energy_class;
  int num_features;
  FEATURE_DATA *feature;
  int i;

  /***** count the features *****/

  /* we know that there will be a mode color for each block */
  num_features = num_total_colour_blocks;

  /* and there will be one for each non-zero histogram entry. Convert
  histogram entries to range [0, FREQ_MAX] while we're at it */
  for (i = 0; i < colmap_size; i++) {
    col_histogram[i] = (int)(rint(FREQ_MAX*(double)col_histogram[i]/(double)(square(image_size))));
    if (col_histogram[i] != 0)
      num_features++;
  }

#ifdef NO_FLAT_FEATURES
  /* count the block entries which are not of class zero (which indicates
  a filter response too low to write), and the non-zero gabor histogram
  entries */
  for (scale = 0; scale < num_gabor_scales; scale++) {
    for (orientation = 0; orientation < num_gabors_per_scale; orientation++) {
      /* blocks */
      for (i = 0; i < square((image_size/gabor_block_size)); i++) {
        if (block_gabor_class[scale][orientation][i] != 0)
          num_features++;
      }

      /* histogram */
      for (energy_class = 0; energy_class < num_gabor_ranges; energy_class++) {
        if (gabor_histogram[scale][orientation][energy_class] != 0)
          num_features++;
      }
    }
  }
#else
  /* here we ARE storing features for filter responses in the lowest
   * energy band, so the number of block features is fixed (and v.
   * large!). We also count the non-zero gabor histogram entries */

  /* blocks */
  num_features += num_gabor_scales*num_gabors_per_scale*square((image_size/gabor_block_size));

  /* histogram */
  for (scale = 0; scale < num_gabor_scales; scale++) {
    for (orientation = 0; orientation < num_gabors_per_scale; orientation++) {
      for (energy_class = 0; energy_class < num_gabor_ranges; energy_class++) {
        if (gabor_histogram[scale][orientation][energy_class] != 0)
          num_features++;
      }
    }
  }
#endif

#ifdef DEBUG
  fprintf(stderr, "Image contains %d features.\n", num_features);
#endif

  /* allocate space to store the features */
  feature = (FEATURE_DATA *)malloc(num_features*sizeof(FEATURE_DATA));

  /***** store the features in the array *****/
  feature_index = 0;

  /* colour block features */
  for (i = 0; i < num_total_colour_blocks; i++) {
    /* note that each block can have colmap_size values */
    feature[feature_index].id = i*colmap_size + (int)block_mode[i];
    feature[feature_index].frequency = FREQ_MAX; /* binary features */
    feature_index++;
  }
  block_features_offset = num_total_colour_blocks*colmap_size;

  /* colour histogram features */
  for (i = 0; i < colmap_size; i++) {
    if (col_histogram[i] != 0) {
      feature[feature_index].id = block_features_offset + i;
      feature[feature_index].frequency = (freq_type)col_histogram[i];
      feature_index++;
    }
  }
  block_features_offset += colmap_size;

  /* gabor block features */
  for (scale = 0; scale < num_gabor_scales; scale++) {
    for (orientation = 0; orientation < num_gabors_per_scale; orientation++) {
      /* blocks */
      for (i = 0; i < square((image_size/gabor_block_size)); i++) {
#ifdef NO_FLAT_FEATURES
        if (block_gabor_class[scale][orientation][i] != 0) {
          feature[feature_index].id = block_features_offset + block_gabor_class[scale][orientation][i];
          feature[feature_index].frequency = FREQ_MAX;
          feature_index++;
        }
#else
        feature[feature_index].id = block_features_offset + block_gabor_class[scale][orientation][i];
        feature[feature_index].frequency = FREQ_MAX;
        feature_index++;
#endif
        block_features_offset += num_gabor_ranges;
      }
    }
  }

  /* gabor histogram features */
  for (scale = 0; scale < num_gabor_scales; scale++) {
    for (orientation = 0; orientation < num_gabors_per_scale; orientation++) {
      for (energy_class = 0; energy_class < num_gabor_ranges; energy_class++) {
        if (gabor_histogram[scale][orientation][energy_class] != 0) {
          feature[feature_index].id = block_features_offset;
          feature[feature_index].frequency = (freq_type)rint(FREQ_MAX*(double)gabor_histogram[scale][orientation][energy_class]/(double)square((image_size/gabor_block_size)));
          feature_index++;
        }
        block_features_offset++;
      }
    }
  }

#ifdef DEBUG
  fprintf(stderr, "%d features found.\n", feature_index);
  fprintf(stderr, "block_features_offset = %d.\n", block_features_offset);
#endif  

    
  /***** now write the file *****/

  /* open it */
  if ((out_file = fopen(out_fname, "wb")) == NULL) {
    fprintf(stderr, "Error opening file %s for writing.\n\n", out_fname);
    return(FILE_OPEN_ERROR);
  }

  /* write the number of features */
  if (fwrite(&num_features, sizeof(int), 1, out_file) != 1) {
    fprintf(stderr, "Error writing file %s.\n\n", out_fname);
    return(FILE_WRITE_ERROR);
  }

  for(i=0;
      i<num_features;
      i++){
    *((float*)&(feature[i].frequency))=((float)feature[i].frequency/(float)FREQ_MAX);//(WM HACK)
  };

  /* write the features, as a block */
  if (fwrite(feature, sizeof(FEATURE_DATA), num_features, out_file) != num_features) {
    fprintf(stderr, "Error writing file %s.\n\n", out_fname);
    return(FILE_WRITE_ERROR);
  }

  if (fclose(out_file) == EOF) {
    fprintf(stderr, "Error closing file %s.\n\n", out_fname);
    return(FILE_CLOSE_ERROR);
  }

  /* everything is OK if we got this far */
  return(PPM_OK);
}