MINC/SoftwareDevelopment/MINC1-voxel-loop-example

this is a dump of an old doc from Peter Neelin from yonks ago

The basic idea is that you call voxel_loop with a list of input files, a list of output files and a function to call:

   /* figure out input files */
   num_outputs = 2;
   output_files = &argv[argc - num_output_files];
   num_inputs = argc - num_outputs - 1;
   input_files = &argv[1];

   /* Set up program_data pointer */
   program_data = MALLOC(sizeof(*program_data));
   program_data->field1 = junk;

   /* Set up loop options. You can just pass it NULL if the defaults
      are okay, but I put some here as examples. */
   loop_options = create_loop_options();
   set_loop_verbose(loop_options, TRUE);
   set_loop_clobber(loop_options, FALSE);
   set_loop_datatype(loop_options, NC_SHORT, TRUE, 0.0, 32000.0);
   set_loop_copy_all_header(loop_options, FALSE);
   set_loop_buffer_size(loop_options, (long) 1024 * max_buffer_size_in_kb);

   /* Do it */
   voxel_loop(num_inputs, input_files, num_outputs, output_files, 
              history_string, loop_options, 
              voxel_function, (void *) program_data);

   /* Free loop options */
   free_loop_options(loop_options);

/******************************************************************/

#define INVALID_DATA (-DBL_MAX)

void voxel_function(void *caller_data, long num_voxels, 
                    int input_num_buffers, int input_vector_length, 
                    double *input_data[],
                    int output_num_buffers, int output_vector_length, 
                    double *output_data[],
                    Loop_Info *loop_info)
/* ARGSUSED */
{
   Program_Data *program_data;
   long ivox, num_values;
   int ibuff;
   double sum1, sum2, value;

   /* Get pointer to program data */
   program_data = (Program_Data *) caller_data;
   num_values = num_voxels * input_vector_length;

   /* Do the good stuff */
   for (ivox=0; ivox < num_values; ivox++) {

      sum1 = sum2 = 0.0;
      for (ibuff=0; ibuff < input_num_buffers; ibuff++) {

         /* Do something useful with the data! */
         value = input_data[ibuff][ivox];
         if (value != INVALID_DATA) {
            sum1++;
            sum2 += value;
         }

      }

      output_data[0] = sum1;
      output_data[1] = (sum1 > 0.0 ? sum2/sum1 : INVALID_DATA);

   }

   return;
}
<syntaxhighlight>

Of course, this is a bad example, since I am loading in a buffer for
each volume and then looping over all the buffers just once. For a lot
of input files I am doing a great deal of unnecessary buffering. For
the case where you have one pass through the data, you should set the
loop option

<syntaxhighlight lang="c">
   set_loop_accumulate(loop_options, TRUE, num_extra_buffers, 
                       start_function, end_function);

This causes voxel_loop to call the voxel function with only one buffer (but once for each file). The idea is that you take the ibuff loop and put it on the outside. However, now you need to keep track of sum1 and sum2 in buffers - well, just use the output buffers for that purpose. What if you only have one output file and need two buffers? Then set num_extra_buffers to 1 - added to the one output buffer gives two available buffers (only the first buffer is written out). The functions start_function and end_function are used to initialize and do final computation on the output buffers. Remember to make sure that your final data is in the first num_outputs buffers after end_function because that is what is written out.