[MINC-users] Re: [geeks] voxel loop
Maxime Descoteaux
mdesco@cim.mcgill.ca
Fri, 28 Nov 2003 09:52:24 -0500 (EST)
Thanks Peter,
My problem is that I run out of memory on full brain mr images. I need
help setting things up? I normally have "Volume" data structures lying
around (too many)! Why?
Well, I do a multiscale vessel extraction. So, at all time and for every
scale, I need six volumes to compute the symmetric Hessian matrix of 2nd
derivatives at every voxel. Ixx, Ixy, Ixz, Iyy, Iyz, Izz of dimensions X x
Y x Z. I do e-value analysis on the Hessian to compute a vesselness
measure.
So, I also have 5 other volumes to save the output of the analysis
1 for the vesselness measure
1 for the x component of the smallest e-vector of the Hessian
1 for the y component of the smallest e-vector of the Hessian
1 for the z component of the smallest e-vector of the Hessian
1 for the scale at which vessel is detected
This is a lot of volumes lying around.... Is voxel_loop a solution to my
problem? I am not quite familiar with the buffer idea and I am not quite
sure on how to set things up...
rough idea:
1) Declare one Volume data structure and compute Ixx. Now, write it out
to "Ixx.mnc". Compute Ixy and write it out to "Ixy.mnc" ....... to
"Izz.mnc".
2) Set my number of input volumes and array of names variables properly.
3) Initialize output volumes names properly and write them out with dummy
initialization
4) Use voxel_loop on the six volumes and define my Voxelfunction to be the
Hessian computation routine...
5) If things are done right, my "vesselness.mnc", "e_x.mnc", ... ,
"scale.mnc" will be filled with the data
Is this really more efficient than to just have volumes lying around? One
last thing, I know the voxel_function can be applied to each voxel one by
one but can a voxel_function be used on a neighborhood of voxels (like a
convolution operation)?
Thanks a lot. This really helps
Max
I don't know if this is clear?
DO you have any suggestion
Thank you very much
Max
Can you help me? DO you have advice or information on this matter?
Thank you very much
> Here's an example with comments that I sent to Andrew in the hopes that it
> could form the basis for some voxel_loop documentation and end up in the
> minc documentation somewhere - I don't know if it has.
>
> Peter
> ----
> Peter Neelin (neelin@bic.mni.mcgill.ca)
>
>
> =======================================================================
>
> 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;
> }
>
>
> 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
>
> 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.
>
> =======================================================================
>
>
>
> _______________________________________________
> MINC-users@bic.mni.mcgill.ca
> http://www.bic.mni.mcgill.ca/mailman/listinfo/minc-users
>
--
Maxime Descoteaux phone: (514) 989-8670
Centre for Intelligent Machines phone: (514) 398-5606
Masters in Computer Science
McGill University
http://www.cim.mcgill.ca/~mdesco