classify.cc

KNN algorithm, and related

  /* initialize the feature_matrix */
  for_less( j, 0, num_features) 
    for_less( i, 0, num_samples) 
      feature_matrix[i][j] = 0.0;   

  /* start traversing the tag file */
  for_less( i, 0, num_samples) {

    /* get the world coordinate from the tag file */
    wx = tags[i][0];
    wy = tags[i][1];
    wz = tags[i][2];
    
    /* convert world into voxel coordinates using first volume */
    convert_3D_world_to_voxel(in_volume[0], wx, wy, wz, &v1, &v2, &v3);
    
    /* NOTE: check to see if the training point (voxel coordinate)
     falls in the volume. Since all the volumes have the same sizes,
     take volume 0, (that is why 'voxel_is_in_volume( v1, v2, v3)'. If
     later a -world switch is supplied, each tag should be tested in
     each feature volume.  At this point, check only for volume 0 */
 
    /* if voxel is in the volume, get the feature volumes values */
    if ( voxel_is_in_volume( v1, v2, v3)) {

      for_less( j, 0, num_features) {

	GET_VALUE_3D(value, in_volume[j], ROUND(v1), ROUND(v2), ROUND(v3));
	feature_matrix[num_adj_samples][j] = value;
      }

      /* convert the label into integer */
      class_column[num_adj_samples] = atoi( labels[i] );

      if ( debug > 5 ) {
	fprintf(stdout, "tag %d (%.1f %.1f %.1f) : ", i, v1, v2, v3); 
	for_less ( j, 0, num_features ) 
	  fprintf(stdout, "%f, ", feature_matrix[num_adj_samples][j]); 
	fprintf(stdout, "-> %d.\n", class_column[num_adj_samples]);
      }

      /* keep track of highest class label to set voxel and image max */
      if ( max_class_index <  class_column[num_adj_samples] )
	max_class_index = class_column[num_adj_samples] ;
	
      num_adj_samples++;  /* adjust new sample size */

    }
    else {
      fprintf(stderr, "tag %d is not in the volume - ignoring\n", i);

    }

  } /* for_less( i, 0, num_samples) */

  /* after checking and ignoring 'out of volume tags', re-adjust the num_samples */
  num_samples = num_adj_samples;

  /*  COUNTING CLASSES COUNTING CLASSES COUNTING CLASSES COUNTING CLASSES */
  /*  COUNTING CLASSES COUNTING CLASSES COUNTING CLASSES COUNTING CLASSES */
  /*  COUNTING CLASSES COUNTING CLASSES COUNTING CLASSES COUNTING CLASSES */
  /*  COUNTING CLASSES COUNTING CLASSES COUNTING CLASSES COUNTING CLASSES */


  /* reserve space for counting classes */
  ALLOC( class_counter, num_samples ); 

  /* initialize class_counter vector */
  for_less( i, 0, num_samples) 
    class_counter[i] = 0;


  /* count for the number of classes */
  for_less( i, 0, num_samples) {

    if ( debug > 4 ) {
      
      /* dump the feature matrix and class_column here */
      fprintf(stdout, "feature[%d] : ", i); 
      for_less ( j, 0, num_features ) 
	fprintf(stdout, "%f, ", feature_matrix[i][j]); 
      fprintf(stdout, "-> %d.\n", class_column[i]);
    }

    /* increase the classes counter */
    class_counter[ class_column[i] ]++; 

  } /* for_less */


  for_less( i, 0, num_samples) {

    if ( class_counter[i] != 0) 
      num_classes++;
  }
  
  /*  COUNTING CLASSES COUNTING CLASSES COUNTING CLASSES COUNTING CLASSES */
  /*  COUNTING CLASSES COUNTING CLASSES COUNTING CLASSES COUNTING CLASSES */
  /*  COUNTING CLASSES COUNTING CLASSES COUNTING CLASSES COUNTING CLASSES */
  /*  COUNTING CLASSES COUNTING CLASSES COUNTING CLASSES COUNTING CLASSES */

  /* reserve a character pointer ( char *) for each class name */
  ALLOC( class_name, num_classes ); 

  /* reserve a array to keep count of each class  */
  ALLOC( class_count, num_classes ); 

  /* initialize the count of each class to zero, and class_name to NULL */
  for_less( i, 0, num_classes) {

    class_count[i] = 0;
    class_name[i] = NULL;
  }

  /* now that you know how many classes exist, get rid of class_counter */
  FREE(class_counter);

  /* NOTE: Now get the name of each class as it occures in
     class_column[i] vector and copy its name from labels[i] array,
     into the class_name array. At the same time, establish a count of
     each class, and put it in the array class_count.  The way this is
     done is as follows: take each sample from class_column[i] vector.
     If its name exists in the class_name array, increase its
     class_count by one, if not, take the first available spot in the
     class_name[i] array, and copy its name there from labels[i]
     array, and increase its class count by one.  */

  /* for each sample in class_column */
  for_less( i, 0, num_samples) {

    /* and for each class */
    for_less( j, 0, num_classes ) {

      /* see if there is an empty spot (NULL), OR it exists (strcmp) . */
      if (labels[i] == NULL) {
	fprintf(stderr, "Error in tag labels; error reading tag file (possibly incompatible tag file format)\n");
	fprintf(stderr, "class_name[%d]: %s labels[%d]: %s\n", 
		j, class_name[j], i, labels[i]);
	exit(EXIT_FAILURE);
      }
	
      if ( class_name[j] == NULL || !strcmp( class_name[j], labels[i] ) ) {

	/* increase the class_count of that class in that spot j */
	class_count[j]++;

	/* if the spot was empty, that it means that this is a new class name,
	   copy the name of the class in the array, by reserving memory 1st*/
	if ( class_name[j] == NULL ) {

	  /* reserve space for the class name */
	  ALLOC( class_name[j], strlen( labels[i] ) );
	  strcpy( class_name[j], labels[i] );
	}

	/* go to the next element in the class_column, don't bother with the
	   rest of the class cases, since you just increase a count */
	break;
      }
    } /* for j */
  } /* for i */

      
  if ( debug >= 1 ) {
    
    fprintf(stdout, "num_classes = %d\n", num_classes);
    fprintf(stdout, "num_features = %d\n", num_features);
    fprintf(stdout, "num_samples = %d\n", num_samples);

    for_less( j, 0, num_classes) {

      fprintf(stdout, "class_name[%d]  = %s\n", j, class_name[j]);
      fprintf(stdout, "class_count[%d] = %d\n", j, class_count[j]);
    }

  }

  /* dump the feature_matrix and class_column if requested */
  if ( dump_features ) {

    for_less ( i, 0, num_samples ) {

      for_less ( j, 0, num_features ) 
	fprintf(stdout, "%f, ", feature_matrix[i][j]); 
      fprintf(stdout, "%d.\n", class_column[i]);

    }

    /* if you are only dumping features, exit after completion */
    if ( dump_features ) 
      exit( EXIT_SUCCESS );


  }

}      


/* ----------------------------- MNI Header -----------------------------------
@NAME       : voxel_is_in_volume
@INPUT      : 
@OUTPUT     : 
@RETURNS    : 
@DESCRIPTION: check to see if a voxel is in the volume (vol 0 is same as all)
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : Sep 22, 1995 ( Vasco KOLLOKIAN)
@MODIFIED   : 
---------------------------------------------------------------------------- */
int voxel_is_in_volume( Real vox1, Real vox2, Real vox3)
{
 
  if ( vox1 < -0.5 || vox1 >= (Real) first_volume_sizes[0] - 0.5) {

    return FALSE;
  }
  
  else if ( vox2 < -0.5 || vox2 >= (Real) first_volume_sizes[1] - 0.5) {
    
    return FALSE;
  }
  
  else if ( vox3 < -0.5 || vox3 >= (Real) first_volume_sizes[2] - 0.5) {
    
    return FALSE;
  }
  else
    return TRUE;
}

/* ----------------------------- MNI Header -----------------------------------
@NAME       : convert_features_to_slice_caching(void)
@INPUT      : 
@OUTPUT     : 
@RETURNS    : 
@DESCRIPTION: convert the feature volume into slice caching
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : Oct 26, 1995 ( Vasco KOLLOKIAN)
@MODIFIED   : 
---------------------------------------------------------------------------- */
void convert_features_to_slice_caching(void) 
{

  int i;

  /* enable slice access per user selected block size def {1,-1,-1} slice*/
  block_sizes[0] = user_block_sizes[0]; 
  block_sizes[1] = user_block_sizes[1];
  block_sizes[2] = user_block_sizes[2];
  
  /* this is for the previously opened random access {1,1,1} feature volumes */
  for_less ( i, 0, num_features ) 
    set_volume_cache_block_sizes(in_volume[i], block_sizes);

  /* And set slice caching for volumes to be created */
  set_default_cache_block_sizes(block_sizes);
}


/* ----------------------------- MNI Header -----------------------------------
@NAME       : write_classified_volume()
@INPUT      : 
@OUTPUT     : 
@RETURNS    : 
@DESCRIPTION: write out the classified volume into a minc file.
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : February 6, 1995 ( Vasco KOLLOKIAN)
@MODIFIED   : 
---------------------------------------------------------------------------- */
void write_classified_volume(void)
{
  Real minval = (output_range[0] == -MAXDOUBLE) ? 0 : output_range[0];
  Real maxval = (output_range[1] == -MAXDOUBLE) ? max_class_index : output_range[1];

  /* write the classified volume to a file */
  if (verbose) 
    fprintf(stdout,"\nWriting classified volume %s to file ...\n", output_filename);
   

  status = output_modified_volume( output_filename, 
				   NC_BYTE, 
				   FALSE, 
				   minval, maxval,
				   classified_volume, 
				   input_filename[0],
				   history, 
				   (minc_output_options *) NULL ) ;

  /* now delete the classified volume so that cache is flushed - as per David */
  delete_volume( classified_volume );

  if ( status != OK )
     exit(EXIT_FAILURE);
} /* write_classified_volume */

/* ----------------------------- MNI Header -----------------------------------
@NAME       : decide_fuzzy_volumes
@INPUT      : 
@OUTPUT     : 
@RETURNS    : 
@DESCRIN    : initialize create_fuzzy_volume array to decide which
              fuzzy volume to initialzie and create
@METHOD     : 
@GLOBAL     :
@CALLS      : 
@CREATED    : Oct 26, 1995 (Vasco KOLLOKIAN)
@MODIFIED   : 
-------------------------------------------------------------------------- */
void decide_fuzzy_volumes(void){

  int k, num_fuzzy_vols;
  int fuzzy_all = !strcmp(fuzzy, "all");

  /* reserve some space for fuzzy volume creation flag */
  ALLOC( create_fuzzy_volume, num_classes );

  /* for each class, set fuzzy volume creation flag to false */
  for_less( k, 0, num_classes ) 
    create_fuzzy_volume[k] = 0;

  /* if a fuzzy path is not defined, make it the current directory */
  if ( fuzzy_path == NULL ) {
    
    ALLOC( fuzzy_path, 2);
    sprintf( fuzzy_path, ".");
  }
  
  /* reserve some space for fuzzy filename prefix and give default name */
  if ( !fuzzy_prefix ) {

    char fuzzy_name[] = "fuzzy_class";
    ALLOC( fuzzy_prefix, strlen(fuzzy_name));
    strcpy( fuzzy_prefix, fuzzy_name );

  }

  /* reserve memory for the fuzzy filename */
  /* the 8 accounts for '/' '_' '.' 'm' 'n' 'c' '\0' 'spare' */

  ALLOC2D( fuzzy_filename, num_classes, (  strlen(fuzzy_path)
					 + strlen(fuzzy_prefix) 
					 + strlen(class_name[0]) + 8) );

  if ( fuzzy_all || ( strlen(fuzzy) > num_classes ) )
    num_fuzzy_vols = num_classes;
  else
    num_fuzzy_vols = strlen(fuzzy);

  /* now for each fuzzy class specification that is not zero, set the flag on.
     make sure you go as far as the shorter of the strlen(fuzzy) or num_classes */
  for_less( k, 0, num_fuzzy_vols ) {

    if ( fuzzy_all || ( fuzzy[k] != '0' ) ) {

      create_fuzzy_volume[k] = 1;

      /* compose the fuzzy volume filename by adding '.mnc' to class_name */
      sprintf( fuzzy_filename[k], "%s/%s_%s.mnc", fuzzy_path, 
	                                          fuzzy_prefix, 
	                                          class_name[k]);

      /* check to see if the volume already exists */
      if ( file_exists(fuzzy_filename[k]) && !clobber_fuzzy ) {

	fprintf(stderr, "%s exists! use -clob_fuzzy to overwrite.\n",fuzzy_filename[k]);
	exit(EXIT_FAILURE);
      }

    }

  }

  if ( debug >= 1 ) {

    for_less ( k, 0, num_fuzzy_vols ) 
      fprintf( stdout, "\ncreate_fuzzy_volume[%d] = %d", k, create_fuzzy_volume[k]);
    
    fprintf( stdout, "\n");

  }



} /* decide_fuzzy_volumes */


/* ----------------------------- MNI Header -----------------------------------
@NAME       : initialize_fuzzy_volumes
@INPUT      : 
@OUTPUT     : 
@RETURNS    : 
@DESCRIN    :      create and initialize empty fuzzy  volumes
@METHOD     : 
@GLOBAL     :
@CALLS      : 
@CREATED    : Sep 14, 1995 (Vasco KOLLOKIAN)
@MODIFIED   : 
-------------------------------------------------------------------------- */
void initialize_fuzzy_volumes(void)
{

  int    k;


  ALLOC( fuzzy_volume, num_classes );

  /* create the fuzzy volume here */   
  for_less( k, 0, num_classes) {
   
    if ( create_fuzzy_volume[k] != 0 ) {

      if (verbose) { 

	fprintf(stdout, "creating fuzzy volume for class %s\n", class_name[k]);
      }

      fuzzy_volume[k] = copy_volume_definition(in_volume[0],
					       fuzzy_type,
					       FALSE,
					       fuzzy_voxel_min, 
					       fuzzy_voxel_max);

      set_volume_voxel_range(fuzzy_volume[k], fuzzy_voxel_min, fuzzy_voxel_max); 
      set_volume_real_range(fuzzy_volume[k], fuzzy_image_min, fuzzy_image_max);


      if ( cache_set ) {
	set_cache_output_volume_parameters(fuzzy_volume[k],
					   fuzzy_filename[k],
					   fuzzy_type, 
					   FALSE, 
					   fuzzy_voxel_min, 
					   fuzzy_voxel_max,
					   input_filename[0], 
					   history, 
					   (minc_output_options *) NULL ) ;



      } /* if ( cache_set ) */

    } /* if ( create_fuzzy ...) */

  } /* for_less */
            
            
} /* initialize_fuzzy_volumes */



/* ----------------------------- MNI Header -----------------------------------
@NAME       : write_fuzzy_volumes()
@INPUT      : 
@OUTPUT     : 
@RETURNS    : 
@DESCRIPTION: write out the fuzzy class confidence levels of the classified volume
              in separate volumes for each class. Volume names are that of
	      class_name[num_classes].
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : Sep 14, 1995 ( Vasco KOLLOKIAN)
@MODIFIED   : 
---------------------------------------------------------------------------- */
void write_fuzzy_volumes(void)
{

  int i;    /* counter */


  /* write the fuzzy volumes to separate files */