MINC/Tutorials/Programming02

From Wikibooks, open books for an open world
Jump to navigation Jump to search

The two MINC coordinate systems.[edit | edit source]

MINC has two coordinate systems: a voxel and a world coordinate system. Voxel coordinates are like arrays, with the origin fixed in one corner of the image and incrementing to the length - 1 of that dimension. World coordinates, on the other hand, reflect the position in real world units of each voxel. Within a particular space, such as Talairach space, each particular coordinate thus has a meaning and is often shared across volumes.

The previous example just obtained a voxel using voxel coordinates - below is the same code modified to obtain a voxel at the origin of the world coordinate system.

#include <minc2.h>
#include <stdio.h>

int main(int argc, char **argv) {
  mihandle_t    minc_volume;
  double        voxel;
  int           result, i;
  double        world_location[3];
  double        dvoxel_location[3];
  unsigned long voxel_location[3];
  
  /* open the volume - first and only 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);
  }

So far, so good - this is almost the same code as in the first example with the addition of some extra variables.

  /* convert world coordinates to voxel coordinates */
  world_location[0] = world_location[1] = world_location[2] = 0.0;
  miconvert_world_to_voxel(minc_volume, world_location, dvoxel_location);

This is where the actual conversion happens, taken care of by the miconvert_world_to_voxel function. There are three arguments: the volume struct (initialized by miopen_volume), an array of doubles containing the world coordinates that are to be converted (initialized to be zeros across the board in the line above), and an array of doubles which will hold the corresponding voxel coordinates.

  /* miconvert_world_to_voxel needs an array of doubles     *
   * but miget_real_value needs unsigned longs - so we cast */
  for (i=0; i<3; i++) { 
    voxel_location[i] = (unsigned long) dvoxel_location[i];
  }

  printf("Voxel location of xyz 0 0 0: %lu %lu %lu\n", 
         voxel_location[0], voxel_location[1], voxel_location[2]);

  /* print the value at that location */
  miget_real_value(minc_volume, voxel_location, 3, &voxel);
  printf("Voxel at xyz 0 0 0 was: %f\n", voxel);

  return(0);
}

The rest of the code should once again seem fairly familiar. There is one slightly awkward cast in there, as miconvert_world_to_voxel assigns the converted voxel coordinates to an array of doubles, but miget_real_value needs an array of unsigned longs. We could do some real interpolation, but we'll be lazy and just cast the double back to the unsigned long.

code[edit | edit source]

#include <minc2.h>
#include <stdio.h>

int main(int argc, char **argv) {
  mihandle_t    minc_volume;
  double        voxel;
  int           result, i;
  double        world_location[3];
  double        dvoxel_location[3];
  unsigned long voxel_location[3];
  
  /* open the volume - first and only 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);
  }

  /* convert world coordinates to voxel coordinates */
  world_location[0] = world_location[1] = world_location[2] = 0.0;
  miconvert_world_to_voxel(minc_volume, world_location, dvoxel_location);

  /* miconvert_world_to_voxel needs an array of doubles     *
   * but miget_real_value needs unsigned longs - so we cast */
  for (i=0; i<3; i++) { 
    voxel_location[i] = (unsigned long) dvoxel_location[i];
  }

  printf("Voxel location of xyz 0 0 0: %lu %lu %lu\n", 
	 voxel_location[0], voxel_location[1], voxel_location[2]);

  /* print the value at that location */
  miget_real_value(minc_volume, voxel_location, 3, &voxel);
  printf("Voxel at xyz 0 0 0 was: %f\n", voxel);

  return(0);
}