--- a/progs/Proglib/voxel_loop.c
+++ b/progs/Proglib/voxel_loop.c
@@ -7,7 +7,13 @@
 @CREATED    : January 10, 1994 (Peter Neelin)
 @MODIFIED   : 
  * $Log: voxel_loop.c,v $
- * Revision 6.1  1999-10-19 14:45:15  neelin
+ * Revision 6.2  2000-09-19 14:36:05  neelin
+ * Added ability for caller to specify functions for allocating and freeing
+ * voxel buffers used in loop. This is particularly useful for embedding
+ * the voxel_loop code in other programs, such as Python, which manage memory
+ * in their own way.
+ *
+ * Revision 6.1  1999/10/19 14:45:15  neelin
  * Fixed Log subsitutions for CVS
  *
  * Revision 6.0  1997/09/12 13:23:41  neelin
@@ -67,7 +73,7 @@
 ---------------------------------------------------------------------------- */
 
 #ifndef lint
-static char rcsid[]="$Header: /private-cvsroot/minc/progs/Proglib/Attic/voxel_loop.c,v 6.1 1999-10-19 14:45:15 neelin Exp $";
+static char rcsid[]="$Header: /private-cvsroot/minc/progs/Proglib/Attic/voxel_loop.c,v 6.2 2000-09-19 14:36:05 neelin Exp $";
 #endif
 
 #include <stdlib.h>
@@ -130,6 +136,7 @@
    void *caller_data;
    Loop_Info *loop_info;
    int is_floating_type;
+   AllocateBufferFunction allocate_buffer_function;
 };
 
 struct Loopfile_Info {
@@ -1209,36 +1216,60 @@
                  chunk_incr, &chunk_num_voxels);
 
    /* Allocate space for buffers */
-   input_buffers = MALLOC(sizeof(*input_buffers) * num_input_buffers);
-   for (ibuff=0; ibuff < num_input_buffers; ibuff++) {
-      input_buffers[ibuff] = MALLOC(sizeof(double) * chunk_num_voxels *
-                                    input_vector_length);
+
+   if (loop_options->allocate_buffer_function != NULL) {
+      loop_options->allocate_buffer_function
+         (loop_options->caller_data, TRUE, 
+          num_input_buffers, chunk_num_voxels, input_vector_length, 
+          &input_buffers, 
+          num_output_files, block_num_voxels, output_vector_length, 
+          &output_buffers, 
+          num_extra_buffers, chunk_num_voxels, output_vector_length, 
+          &extra_buffers, 
+          loop_options->loop_info);
+      
    }
+   else {
+
+      /* Allocate input buffers */
+      input_buffers = MALLOC(sizeof(*input_buffers) * num_input_buffers);
+      for (ibuff=0; ibuff < num_input_buffers; ibuff++) {
+         input_buffers[ibuff] = MALLOC(sizeof(double) * chunk_num_voxels *
+                                       input_vector_length);
+      }
+
+      /* Allocate output buffers */
+      if (num_output_files > 0) {
+         output_buffers = MALLOC(sizeof(*output_buffers) * num_output_files);
+         for (ibuff=0; ibuff < num_output_files; ibuff++) {
+            output_buffers[ibuff] = MALLOC(sizeof(double) * block_num_voxels *
+                                           output_vector_length);
+         }
+      }
+
+      /* Allocate extra buffers */
+      if (num_extra_buffers > 0) {
+         extra_buffers = MALLOC(sizeof(*extra_buffers) * num_extra_buffers);
+         for (ibuff=0; ibuff < num_extra_buffers; ibuff++) {
+            extra_buffers[ibuff] = MALLOC(sizeof(double) * chunk_num_voxels *
+                                          output_vector_length);
+         }
+      }
+
+   }
+
+   /* Set up the results pointers */
    if (num_output_buffers > 0) {
       results_buffers = MALLOC(sizeof(*results_buffers) * num_output_buffers);
-   }
-   else
-      results_buffers = NULL;
-   if (num_output_files > 0) {
-      output_buffers = MALLOC(sizeof(output_buffers) * num_output_files);
-      for (ibuff=0; ibuff < num_output_files; ibuff++) {
-         output_buffers[ibuff] = MALLOC(sizeof(double) * block_num_voxels *
-                                        output_vector_length);
-         results_buffers[ibuff] = output_buffers[ibuff];
+      for (ibuff=0; ibuff < num_output_buffers; ibuff++) {
+         if (ibuff < num_output_files) {
+            results_buffers[ibuff] = output_buffers[ibuff];
+         }
+         else {
+            results_buffers[ibuff] = extra_buffers[ibuff-num_output_files];
+         }
       }
    }
-   else
-      output_buffers = NULL;
-   if (num_extra_buffers > 0) {
-      extra_buffers = MALLOC(sizeof(extra_buffers) * num_extra_buffers);
-      for (ibuff=0; ibuff < num_extra_buffers; ibuff++) {
-         extra_buffers[ibuff] = MALLOC(sizeof(double) * chunk_num_voxels *
-                                       output_vector_length);
-         results_buffers[ibuff+num_output_files] = extra_buffers[ibuff];
-      }
-   }
-   else
-      extra_buffers = NULL;
 
    /* Initialize global min and max */
    if (num_output_files > 0) {
@@ -1507,26 +1538,51 @@
       (void) fflush(stderr);
    }
 
-   /* Free buffers */
-   for (ibuff=0; ibuff < num_input_buffers; ibuff++) {
-      FREE(input_buffers[ibuff]);
-   }
-   FREE(input_buffers);
-   if (num_output_files > 0) {
-      for (ibuff=0; ibuff < num_output_files; ibuff++) {
-         FREE(output_buffers[ibuff]);
-      }
-      FREE(output_buffers);
-   }
-   if (num_extra_buffers > 0) {
-      for (ibuff=0; ibuff < num_extra_buffers; ibuff++) {
-         FREE(extra_buffers[ibuff]);
-      }
-      FREE(extra_buffers);
-   }
+   /* Free results pointer array, but not its buffers, since these
+      were allocate as output_buffers and extra_buffers */
    if (num_output_buffers > 0) {
       FREE(results_buffers);
    }
+
+   /* Free the buffers */
+   if (loop_options->allocate_buffer_function != NULL) {
+      loop_options->allocate_buffer_function
+         (loop_options->caller_data, FALSE, 
+          num_input_buffers, chunk_num_voxels, input_vector_length, 
+          &input_buffers, 
+          num_output_files, block_num_voxels, output_vector_length, 
+          &output_buffers, 
+          num_extra_buffers, chunk_num_voxels, output_vector_length, 
+          &extra_buffers, 
+          loop_options->loop_info);
+   }
+   else {
+
+      /* Free input buffers */
+      for (ibuff=0; ibuff < num_input_buffers; ibuff++) {
+         FREE(input_buffers[ibuff]);
+      }
+      FREE(input_buffers);
+
+      /* Free output buffers */
+      if (num_output_files > 0) {
+         for (ibuff=0; ibuff < num_output_files; ibuff++) {
+            FREE(output_buffers[ibuff]);
+         }
+         FREE(output_buffers);
+      }
+
+      /* Free extra buffers */
+      if (num_extra_buffers > 0) {
+         for (ibuff=0; ibuff < num_extra_buffers; ibuff++) {
+            FREE(extra_buffers[ibuff]);
+         }
+         FREE(extra_buffers);
+      }
+
+   }
+
+   /* Free max and min arrays */
    if (num_output_files > 0) {
       FREE(global_minimum);
       FREE(global_maximum);
@@ -2549,6 +2605,8 @@
    loop_options->caller_data = NULL;
    loop_options->loop_info = create_loop_info();
 
+   loop_options->allocate_buffer_function = NULL;
+
    /* Return the structure pointer */
    return loop_options;
 }
@@ -2932,6 +2990,27 @@
    
 }
 
+/* ----------------------------- MNI Header -----------------------------------
+@NAME       : set_loop_allocate_buffer_function
+@INPUT      : loop_options - user options for looping
+              allocate_buffer_function - user function that allocates and frees
+                 input_data
+@OUTPUT     : (none)
+@RETURNS    : (nothing)
+@DESCRIPTION: Routine to set the function that allocates and frees input_data
+@METHOD     : 
+@GLOBALS    : 
+@CALLS      : 
+@CREATED    : Aug 29, 2000 (J. Taylor)
+@MODIFIED   : Sept. 19, 2000 (P. Neelin)
+---------------------------------------------------------------------------- */
+public void set_loop_allocate_buffer_function(Loop_Options *loop_options, 
+                         AllocateBufferFunction allocate_buffer_function)
+
+{
+   loop_options->allocate_buffer_function = allocate_buffer_function;
+}
+
 /* ------------ Routines to set and get loop info ------------ */
 
 /* ----------------------------- MNI Header -----------------------------------

--- a/progs/Proglib/voxel_loop.h
+++ b/progs/Proglib/voxel_loop.h
@@ -5,7 +5,13 @@
 @CREATED    : January 10, 1994 (Peter Neelin)
 @MODIFIED   : 
  * $Log: voxel_loop.h,v $
- * Revision 6.1  1999-10-19 14:45:16  neelin
+ * Revision 6.2  2000-09-19 14:36:05  neelin
+ * Added ability for caller to specify functions for allocating and freeing
+ * voxel buffers used in loop. This is particularly useful for embedding
+ * the voxel_loop code in other programs, such as Python, which manage memory
+ * in their own way.
+ *
+ * Revision 6.1  1999/10/19 14:45:16  neelin
  * Fixed Log subsitutions for CVS
  *
  * Revision 6.0  1997/09/12 13:23:41  neelin
@@ -78,13 +84,13 @@
               num_voxels - number of voxels to process. Note that the
                  total number of input values is 
                  num_voxels * input_vector_length.
-              input_num_buffers - number of input buffers to handle
+              num_input_buffers - number of input buffers to handle
                  on this call - either total number of input files or 1 
                  (for accumulating over files).
               input_vector_length - length of input vector.
               input_data - array of pointers to input buffers (1 for
                  each input file, unless we are accumulating).
-              output_num_buffers - number of output buffers to handle
+              num_output_buffers - number of output buffers to handle
                  on this call - will be the total number of output files
                  unless we are accumulating over files (see 
                  set_loop_accumulate).
@@ -102,8 +108,8 @@
 ---------------------------------------------------------------------------- */
 typedef void (*VoxelFunction) 
      (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[],
+      int num_input_buffers, int input_vector_length, double *input_data[],
+      int num_output_buffers, int output_vector_length, double *output_data[],
       Loop_Info *loop_info);
 
 /* ----------------------------- MNI Header -----------------------------------
@@ -188,6 +194,43 @@
       int output_num_buffers, int output_vector_length, double *output_data[],
       Loop_Info *loop_info);
 
+/* ----------------------------- MNI Header -----------------------------------
+@NAME       : AllocateBufferFunction
+@INPUT      : caller_data - pointer to client data.
+              do_allocations - if TRUE, allocate data, if FALSE, free data
+              num_input_buffers - number of input buffers to allocate
+              num_input_voxels - number of voxels in each input buffer
+              input_vector_length - number of values per input voxel
+              num_output_buffers - number of output buffers to allocate
+              num_output_voxels - number of voxels in each output buffer
+              output_vector_length - number of values per output voxel
+              num_extra_buffers - number of working buffers to allocate
+              num_extra_voxels - number of voxels in each extra buffer
+              extra_vector_length - number of values per extra voxel
+              loop_info - pointer that can be passed to functions returning
+                 looping information
+              input_buffers, output_buffers, extra_buffers - if do_allocation 
+                 is FALSE, then arrays of pointers to buffers that need
+                 to be freed. The pointer arrays should also be freed.
+@OUTPUT     : input_buffers, output_buffers, extra_buffers - if do_allocation 
+                 is TRUE, then these arrays of buffers should be allocated
+                 (both the array of pointers and the buffers). The pointer 
+                 array should have length num_xxx_buffers and each buffer
+                 should have length num_xxx_voxels*xxx_vector_length and be
+                 of type double.
+@RETURNS    : (nothing)
+@DESCRIPTION: Typedef for function called by voxel_loop to allocate and
+              free buffers.
+---------------------------------------------------------------------------- */
+typedef void (*AllocateBufferFunction) 
+     (void *caller_data, int do_allocation,
+      int num_input_buffers, int num_input_voxels, int input_vector_length, 
+      double ***input_buffers,
+      int num_output_buffers, int num_output_voxels, int output_vector_length, 
+      double ***output_buffers,
+      int num_extra_buffers, int num_extra_voxels, int extra_vector_length, 
+      double ***extra_buffers,
+      Loop_Info *loop_info);
 
 /* Function declarations */
 public void voxel_loop(int num_input_files, char *input_files[], 
@@ -231,6 +274,9 @@
                                 int num_extra_buffers,
                                 VoxelStartFunction start_function,
                                 VoxelFinishFunction finish_function);
+public void set_loop_allocate_buffer_function(Loop_Options *loop_options, 
+                         AllocateBufferFunction allocate_buffer_function)
+;
 public void get_info_shape(Loop_Info *loop_info, int ndims,
                            long start[], long count[]);
 public int get_info_current_file(Loop_Info *loop_info);