extract_block_features.c File Reference

#include <stdio.h>
#include <malloc.h>
#include <math.h>
#include <stdint.h>
#include <ppm.h>
#include <gift_features.h>

Go to the source code of this file.

Defines
#define	square(x)   ((x)*(x))
#define	num_gabor_ranges   10
#define	num_gabor_scales   3
#define	num_gabors_per_scale   4
#define	OLDFIXED   1
#define	image_size   256
#define	num_total_colour_blocks   (256+64+16+4)
#define	smallest_colour_block   16
#define	gabor_block_size   16
#define	num_colour_scales   4
#define	num_blocks_at_scale(i)   (256>>i*2)
#define	gabor_ranges(i)   (i+2+99989*(((i>>3)&1)&(i&1)))
Functions
void	init_feature_variables (uint32_t colmap_size, uint32_t **col_counts, uint32_t ***block_gabor_class, uint32_t ***gabor_histogram)
void	extract_gabor_features (PPM *im_hsv, uint32_t ***block_gabor_class, uint32_t ***gabor_histogram)
void	extract_mode_features (PPM *im, uint32_t *colmap, uint32_t colmap_size, uint32_t **col_counts, byte *block_mode, uint32_t *col_histogram)
enum ppm_error	write_mode_features (char *out_fname, uint32_t colmap_size, uint32_t *col_histogram, uint32_t ***block_gabor_class, uint32_t ***gabor_histogram, byte *block_mode)
void	fts2blocks (char *fts_fname, byte *block_mode, uint32_t *col_histogram)
enum ppm_error	write_feature_descriptions (FILE *out_file, int *colmap, int colmap_size)

---

Define Documentation

#define gabor_block_size   16

Definition at line 79 of file extract_block_features.c.

Referenced by extract_gabor_features(), init_feature_variables(), write_feature_descriptions(), and write_mode_features().

#define gabor_ranges

(

i

)

(i+2+99989*(((i>>3)&1)&(i&1)))

Definition at line 98 of file extract_block_features.c.

Referenced by extract_gabor_features(), and write_feature_descriptions().

#define image_size   256

Definition at line 71 of file extract_block_features.c.

Referenced by extract_mode_features(), fts2blocks(), init_feature_variables(), write_feature_descriptions(), and write_mode_features().

#define num_blocks_at_scale

(

i

)

(256>>i*2)

Definition at line 85 of file extract_block_features.c.

Referenced by init_feature_variables(), and write_feature_descriptions().

#define num_colour_scales   4

Definition at line 84 of file extract_block_features.c.

Referenced by extract_mode_features(), fts2blocks(), init_feature_variables(), and write_feature_descriptions().

#define num_gabor_ranges   10

Definition at line 56 of file extract_block_features.c.

Referenced by extract_gabor_features(), init_feature_variables(), write_feature_descriptions(), and write_mode_features().

#define num_gabor_scales   3

Definition at line 61 of file extract_block_features.c.

Referenced by create_filter_kernels(), extract_gabor_features(), init_feature_variables(), main(), write_feature_descriptions(), and write_mode_features().

#define num_gabors_per_scale   4

Definition at line 63 of file extract_block_features.c.

Referenced by create_filter_kernels(), extract_gabor_features(), gabor_filter(), init_feature_variables(), write_feature_descriptions(), and write_mode_features().

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

Definition at line 73 of file extract_block_features.c.

Referenced by fts2blocks(), main(), and write_mode_features().

#define OLDFIXED   1

Definition at line 66 of file extract_block_features.c.

#define smallest_colour_block   16

Definition at line 75 of file extract_block_features.c.

Referenced by extract_mode_features(), fts2blocks(), and write_feature_descriptions().

#define square

(

x

)

((x)*(x))

Definition at line 48 of file extract_block_features.c.

Referenced by extract_gabor_features(), extract_mode_features(), fts2blocks(), init_feature_variables(), write_feature_descriptions(), and write_mode_features().

---

Function Documentation

void extract_gabor_features	(	PPM *	im_hsv,
		uint32_t ***	block_gabor_class,
		uint32_t ***	gabor_histogram
	)

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 *	im,
		uint32_t *	colmap,
		uint32_t	colmap_size,
		uint32_t **	col_counts,
		byte *	block_mode,
		uint32_t *	col_histogram
	)

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 fts2blocks	(	char *	fts_fname,
		byte *	block_mode,
		uint32_t *	col_histogram
	)

Definition at line 532 of file extract_block_features.c.

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

Referenced by main().

                                                                              {
  
  FILE *fts_file, *outfile;
  int num_features, feature_index;
  FEATURE_DATA *feature;
  int block_features_offset;
  int block_size, num_blocks, scale;
  int *colmap;
  int colmap_size = 18*3*3 + 4;
  int x, y, i, j, k;
  PPM *hsv_image, *qimage, *rgb_image;
  enum ppm_error the_error;
  char out_fname[256];
 
  /* read the number of features */
  fts_file = fopen(fts_fname, "rb");
  fread(&num_features, sizeof(int), 1, fts_file);
  
  /* allocate space for the features */
  feature = (FEATURE_DATA *)malloc(num_features*sizeof(FEATURE_DATA));
 
  /* read the features */
  fread(feature, sizeof(FEATURE_DATA), num_features, fts_file);

  fclose(fts_file);

  /* "unpack" the features in the array */
  feature_index = 0;

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

  /* colour histogram features */
  for (i = 0; i < colmap_size; i++)
    col_histogram[i] = 0;

  for (; feature_index < num_features; feature_index++) {
    i = feature[feature_index].id - block_features_offset;
    col_histogram[i] = (byte)(feature[feature_index].frequency);
  }

  /* make a dummy image so that we can generate a colour map */
    hsv_image = new_ppm();
    hsv_image->type = PPM_RAW;
    hsv_image->width = 256;
    hsv_image->height = 256;
    hsv_image->max_col_comp = 255;
    hsv_image->bytes_per_pixel = 3;
    hsv_image->pixel = (byte *)malloc(3*256*256*sizeof(byte));
  hsv_quantize_ppm(hsv_image, &qimage, &colmap, 18, 3, 3, 4);
  destroy_ppm(&hsv_image);

  /* now write the block images */
  k = 0;
    block_size = smallest_colour_block;
    num_blocks = square(image_size/block_size);
  for (scale = 0; scale < num_colour_scales; scale++) {
    for (y = 0; y < 256; y += block_size) {
    for (x = 0; x < 256; x += block_size) {
      for (i = 0; i < block_size; i++) {
      for (j = 0; j < block_size; j++) {
        qimage->pixel[(y + j)*256 + (x + i)] = block_mode[k];
      }
      }
      k++;
    }
    }
    /* convert from the colour map to an HSV image */
    if ((the_error = colmap2rgb_ppm(qimage, colmap, colmap_size, &hsv_image)) != PPM_OK) {
      ppm_handle_error(the_error);
      exit(1);
    }
    /* convert to RGB */
    if ((the_error = hsv2rgb_ppm(hsv_image, &rgb_image)) != PPM_OK) {
      ppm_handle_error(the_error);
      exit(1);
    }
    /* write it */
    sprintf(out_fname, "fts2blocks_%dx%d.ppm", block_size, block_size);
    outfile = fopen(out_fname, "wb");
    if ((the_error = write_ppm(outfile, rgb_image, PPM_RAW)) != PPM_OK) {
      ppm_handle_error(the_error);
      exit(1);
    }
    fclose(outfile);
    destroy_ppm(&hsv_image);
    destroy_ppm(&rgb_image);
    block_size *= 2;
    num_blocks /= 4;
  }
}

void init_feature_variables	(	uint32_t	colmap_size,
		uint32_t **	col_counts,
		uint32_t ***	block_gabor_class,
		uint32_t ***	gabor_histogram
	)

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));
    }
  }
}

enum ppm_error write_feature_descriptions	(	FILE *	out_file,
		int *	colmap,
		int	colmap_size
	)

Definition at line 629 of file extract_block_features.c.

References COL_HST, COL_POS, gabor_block_size, GABOR_HST, GABOR_POS, gabor_ranges, HUE, image_size, num_blocks_at_scale, num_colour_scales, num_gabor_ranges, num_gabor_scales, num_gabors_per_scale, PPM_OK, SATURATION, smallest_colour_block, square, and VALUE.

Referenced by main().

                                                                                        {

  int feature_index = 0;
  int block_size;
  int i, j, k;
  int scale, orientation;

  /* block features */
  block_size = smallest_colour_block;
  for (i = 0; i < num_colour_scales; i++) {
    for (j = 0; j < num_blocks_at_scale(i); j++) {
      for (k = 0; k < colmap_size; k++) {
        fprintf(out_file, "%d %d COL_POS block size = %dx%d position = %d H,S,V = %d, %d, %d\n", feature_index, COL_POS, block_size, block_size, j, colmap[3*k + HUE], colmap[3*k + SATURATION], colmap[3*k + VALUE]);
        feature_index++;
      }
    }
    block_size *= 2;
  }

  /* colour histogram features */
  for (i = 0; i < colmap_size; i++) {
    fprintf(out_file, "%d %d COL_HST H,S,V = %d, %d, %d\n", feature_index, COL_HST, colmap[3*i + HUE], colmap[3*i + SATURATION], colmap[3*i + VALUE]);
    feature_index++;
  }
    
  /* 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++) {
        for (j = 0; j < num_gabor_ranges; j++) {
#if num_gabor_ranges==16 /* technically, we should be checking the type.. */
          fprintf(out_file, "%d %d GABOR_POS block size = %dx%d position = %d SCALE, ORIENTATION, ENERGY = %d, %d, %f\n", feature_index, GABOR_POS, gabor_block_size, gabor_block_size, i, scale, orientation, gabor_ranges(j));
#else
#if num_gabor_ranges==10
          fprintf(out_file, "%d %d GABOR_POS block size = %dx%d position = %d SCALE, ORIENTATION, ENERGY = %d, %d, %d.000000\n", feature_index, GABOR_POS, gabor_block_size, gabor_block_size, i, scale, orientation, gabor_ranges(j));
#endif
#endif
        feature_index++;
        }
      }
    }
  }

  /* gabor histogram features */
  for (scale = 0; scale < num_gabor_scales; scale++) {
    for (orientation = 0; orientation < num_gabors_per_scale; orientation++) {
      for (j = 0; j < num_gabor_ranges; j++) {
        fprintf(out_file, "%d %d GABOR_HST SCALE, ORIENTATION, ENERGY UPPER BOUND = %d, %d, %f\n", feature_index, GABOR_HST, scale, orientation, j);
        feature_index++;
      }
    }
  }

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

enum ppm_error write_mode_features	(	char *	out_fname,
		uint32_t	colmap_size,
		uint32_t *	col_histogram,
		uint32_t ***	block_gabor_class,
		uint32_t ***	gabor_histogram,
		byte *	block_mode
	)

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);
}