MINC/Tutorials/Programming05

Working with hyperslabsEdit

The previous tutorial stepped through every voxel of a volume, added 1 to it, and wrote the output back to the same volume. Here we will improve on that example: the output will be written to a new volume, and we will use hyperslabs for the data processing.

A hyperslab is a central concept in minc. It allows you to pull an arbitrary (but contiguous) chunk of data out of a volume, do something to it, and write it back to the volume. This is a lot faster than individually getting each voxel; it does, however, require more memory.

So on to the example: this code takes a volume, creates a copy of it the same in every way except for adding 1 to each voxel. The new volume will also always use global scaling, regardless of whether the input had slice scaling or not.

#include <minc2.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <string.h>

int main(int argc, char **argv) {
  mihandle_t    minc_volume, new_volume;
  midimhandle_t dimensions[3], *dimensions_new;
  double        min, max;
  int           result, i;
  unsigned long start[3], count[3];
  unsigned int  sizes[3];
  double        *slab;

  if (argc != 3) {
    fprintf(stderr, "Usage: input.mnc output.mnc\n");
  }

  /* open the volume - first command line argument */
  result = miopen_volume(argv[1], MI2_OPEN_READ, &minc_volume);
  /* check for error on opening */
  if (result != MI_NOERROR) {
    fprintf(stderr, "Error opening input file: %d.\n", result);
    return(1);
  }

  /* get the dimension sizes */
  miget_volume_dimensions(minc_volume, MI_DIMCLASS_SPATIAL,
			  MI_DIMATTR_ALL, MI_DIMORDER_FILE,
			  3, dimensions);
  miget_dimension_sizes(dimensions, 3, sizes);
  printf("Volume sizes: %u %u %u\n", sizes[0], sizes[1], sizes[2]);

So far, so good - none of this code should be new to you. Do note the creation of two new variables for the new volume and new dimensions.

  /* allocate new dimensions */
  dimensions_new = malloc(sizeof(midimhandle_t) * 3);

  /* copy the dimensions */
  for(i=0; i < 3; i++) {
    micopy_dimension(dimensions[i], &dimensions_new[i]);
  }

Here we allocate memory for the dimensions of the new volume, and copy the dimension definitions from the old volume.

  /* allocate memory for the hyperslab - make it size of entire volume */
  slab = malloc(sizeof(double) * sizes[0] * sizes[1] * sizes[2]);

  start[0] = start[1] = start[2] = 0.0;
  if (miget_real_value_hyperslab(minc_volume, MI_TYPE_DOUBLE,
				 start, count, slab) != MI_NOERROR) {
    fprintf(stderr, "Error getting hyperslab\n");
    return(1);
  }

The memory for the hyperslab is allocated and the hyperslab obtained - in this case the hyperslab is the entire volume. The key function is miget_real_value_hyperslab, which takes five arguments: the minc volume handle, the type you want the data to be in, an array of starts (all zeroes in this case), an array of counts( the volume sizes), and the address of where to write the data.


  /* set min and max to the first voxel plus 1 
   *  (since we add one to everything anyway) */
  min = slab[0]+1;
  max = slab[0]+1;

  /* loop over all voxels */
  for (i=0; i < sizes[0] * sizes[1] * sizes[2]; i++) {
    /* increment voxel by 1 */
    slab[i] += 1;
    
    /* check if min or max should be changed */
    if (slab[i] < min) {
      min = slab[i];
    }
    else if (slab[i] > max) {
      max = slab[i];
    }
  }

And now we loop over every voxel in the hyperslab (note how it is dealt with as a one-dimensional array), making sure to keep the modified image min and max.

  /* create the new volume */
  if (micreate_volume(argv[2], 3, dimensions_new, MI_TYPE_UBYTE,
		      MI_CLASS_REAL, NULL, &new_volume) != MI_NOERROR) {
    fprintf(stderr, "Error creating new volume\n");
    return(1);
  }

We then create the actual volume using micreate_volume. It takes seven arguments: the filename, the number of dimensions, an array of the actual dimensions (created through the copy operation above), the data type (set to unsigned byte here), the data class, special volume properties (set to NULL here), and the address of the volume handle.

  /* create the data for the new volume */
  if (micreate_volume_image(new_volume) != MI_NOERROR) {
    fprintf(stderr, "Error creating volume data\n");
    return(1);
  }

Now we make sure that the data is created inside the minc volume handle using micreate_volume_image.

  /* set valid and real range */
  miset_volume_valid_range(new_volume, 255, 0);
  miset_volume_range(new_volume, max, min);

  /* write the modified hyperslab to the file */
  if (miset_real_value_hyperslab(new_volume, MI_TYPE_DOUBLE,
				 start, count, slab) != MI_NOERROR) {
    fprintf(stderr, "Error setting hyperslab\n");
    return(1);
  }

We then set the valid range - from 0 to 255 since this is unsigned byte data - and the volume range, which was determined by keeping track of the min and max when looping over voxels. The hyperslab is the placed back into the volume using miset_real_value_hyperslab, which takes the exact same arguments as miget_real_value_hyperslab.

  /* closes the volume and makes sure all data is written to file */
  miclose_volume(minc_volume);
  miclose_volume(new_volume);

  /* free memory */
  free(dimensions_new);
  free(slab);

  return(0);
}

And finally we close the volumes, free the allocated memory, and exit.

codeEdit

#include <minc2.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <string.h>

int main(int argc, char **argv) {
  mihandle_t    minc_volume, new_volume;
  midimhandle_t dimensions[3], *dimensions_new;
  double        min, max;
  int           result, i;
  unsigned long start[3], count[3];
  unsigned int  sizes[3];
  double        *slab;

  if (argc != 3) {
    fprintf(stderr, "Usage: input.mnc output.mnc\n");
  }

  /* open the volume - first command line argument */
  result = miopen_volume(argv[1], MI2_OPEN_READ, &minc_volume);
  /* check for error on opening */
  if (result != MI_NOERROR) {
    fprintf(stderr, "Error opening input file: %d.\n", result);
    return(1);
  }

  /* get the dimension sizes */
  miget_volume_dimensions(minc_volume, MI_DIMCLASS_SPATIAL,
			  MI_DIMATTR_ALL, MI_DIMORDER_FILE,
			  3, dimensions);
  miget_dimension_sizes(dimensions, 3, sizes);
  printf("Volume sizes: %u %u %u\n", sizes[0], sizes[1], sizes[2]);

  /* allocate new dimensions */
  dimensions_new = malloc(sizeof(midimhandle_t) * 3);

  /* copy the dimensions */
  for(i=0; i < 3; i++) {
    micopy_dimension(dimensions[i], &dimensions_new[i]);
  }

  start[0] = start[1] = start[2] = 0;

  for (i=0; i < 3; i++) {
    count[i] = (unsigned long) sizes[i];
  }

  /* allocate memory for the hyperslab - make it size of entire volume */
  slab = malloc(sizeof(double) * sizes[0] * sizes[1] * sizes[2]);

  if (miget_real_value_hyperslab(minc_volume, MI_TYPE_DOUBLE,
				 start, count, slab) != MI_NOERROR) {
    fprintf(stderr, "Error getting hyperslab\n");
    return(1);
  }
  
  /* set min and max to the first voxel plus 1 
   *  (since we add one to everything anyway) */
  min = slab[i]+1;
  max = slab[i]+1;

  /* loop over all voxels */
  for (i=0; i < sizes[0] * sizes[1] * sizes[2]; i++) {
    /* increment voxel by 1 */
    slab[i] += 1;
    
    /* check if min or max should be changed */
    if (slab[i] < min) {
      min = slab[i];
    }
    else if (slab[i] > max) {
      max = slab[i];
    }
  }

  /* create the new volume */
  if (micreate_volume(argv[2], 3, dimensions_new, MI_TYPE_UBYTE,
		      MI_CLASS_REAL, NULL, &new_volume) != MI_NOERROR) {
    fprintf(stderr, "Error creating new volume\n");
    return(1);
  }
  
  /* create the data for the new volume */
  if (micreate_volume_image(new_volume) != MI_NOERROR) {
    fprintf(stderr, "Error creating volume data\n");
    return(1);
  }

  /* set valid and real range */
  miset_volume_valid_range(new_volume, 255, 0);
  miset_volume_range(new_volume, max, min);

  /* write the modified hyperslab to the file */
  if (miset_real_value_hyperslab(new_volume, MI_TYPE_DOUBLE,
				 start, count, slab) != MI_NOERROR) {
    fprintf(stderr, "Error setting hyperslab\n");
    return(1);
  }
  
  /* closes the volume and makes sure all data is written to file */
  miclose_volume(minc_volume);
  miclose_volume(new_volume);

  /* free memory */
  free(dimensions_new);
  free(slab);

  return(0);
}