Mercurial > hg > minc-tools
changeset 2637:0b9a5b816213
* added mincmorph
* updated README and COPYING
| author | Andrew L Janke <a.janke@gmail.com> |
|---|---|
| date | Wed, 14 Mar 2012 02:57:00 +1000 |
| parents | b0e5228b9c74 |
| children | 9a50bb928cad |
| files | COPYING ChangeLog Makefile.am README progs/mincmorph/ChangeLog progs/mincmorph/README progs/mincmorph/kernel_io.c progs/mincmorph/kernel_io.h progs/mincmorph/kernel_ops.c progs/mincmorph/kernel_ops.h progs/mincmorph/kernels/3x3_4-conn.kern progs/mincmorph/kernels/3x3_8-conn.kern progs/mincmorph/kernels/3x3x3_26-conn.kern progs/mincmorph/kernels/3x3x3_6-conn.kern progs/mincmorph/mincmorph.c |
| diffstat | 15 files changed, 2193 insertions(+), 2 deletions(-) [+] |
line wrap: on
line diff
--- a/COPYING +++ b/COPYING @@ -1,5 +1,5 @@ -Copyright Peter Neelin and David MacDonald, McConnell Brain -Imaging Centre, Montreal Neurological Institute, McGill University. +Copyright the MINC developers, McConnell Brain Imaging Centre, +Montreal Neurological Institute, McGill University. Permission to use, copy, modify, and distribute this software and its documentation for any purpose and without fee is hereby granted,
--- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,6 @@ + +2012-03-14 Andrew L Janke <a.janke@gmail.com> + * added progs/mincmorph -- Release of minc 2.1.10 --
--- a/Makefile.am +++ b/Makefile.am @@ -85,6 +85,8 @@ progs/minccalc/errx.h \ progs/minccalc/node.h \ progs/minccalc/gram.h \ + progs/mincmorph/kernel_io.h \ + progs/mincmorph/kernel_ops.h \ progs/mincresample/mincresample.h \ progs/mincreshape/mincreshape.h \ libsrc2/minc2_private.h \ @@ -202,6 +204,7 @@ mincmakescalar \ mincmakevector \ mincmath \ + mincmorph \ mincresample \ mincreshape \ mincstats \ @@ -324,6 +327,12 @@ mincmath_SOURCES = progs/mincmath/mincmath.c +mincmorph_CFLAGS = -Iprogs/mincmorph -I$(srcdir)/progs/mincmorph +mincmorph_SOURCES = \ + progs/mincmorph/mincmorph.c \ + progs/mincmorph/kernel_io.c \ + progs/mincmorph/kernel_ops.c + mincresample_SOURCES = \ progs/mincresample/mincresample.c \ progs/mincresample/resample_volumes.c \
--- a/README +++ b/README @@ -39,6 +39,10 @@ MINC requires that NetCDF and HDF5 are installed. +MINC development is performed on github + + http://github.com/BIC-MNI/minc + SUPPORT FOR HDF5 "MINC 2.0" FORMAT FILES ----------------------------------------
new file mode 100644 --- /dev/null +++ b/progs/mincmorph/ChangeLog @@ -0,0 +1,20 @@ +2009-11-10 Andrew L Janke <a.janke@gmail.com> + * Added local correlation option (I[]) + * changed to VIO_ code for clean MINC2 build + +2005-09-23 Bert Vincent <bert@bic.mni.mcgill.ca> + * Fix reading of kernel files + * Force calc_volume_range() to return with a max value strictly + greater than the min. + +2004-10-08 Andrew L. Janke <a.janke@gmail.com> + * Added 4 inbuilt kernels + +2004-05-25 Andrew L. Janke <a.janke@gmail.com> + * Added -median_dilate option for mark + +2004-04-14 Andrew L. Janke <a.janke@gmail.com> + * Fixed bug in get_string_from_string() mallocing + +2003-12-17 Andrew L. Janke <a.janke@gmail.com> + * Added support for grey level dilation and erosion
new file mode 100644 --- /dev/null +++ b/progs/mincmorph/README @@ -0,0 +1,16 @@ +mincmorph - morphological operations 101 + +Comments to: +Andrew Janke - a.janke@gmail.com + +Much of the theory behind what this does can be found at these web sites.. + http://www.mmorph.com/faq.htm + http://www.dai.ed.ac.uk/HIPR2/skeleton.htm + + SNF - http://www.umiacs.umd.edu/research/EXPAR/papers/3449/node7.html#SECTION00041000000000000000 + + CC8 - http://www.dai.ed.ac.uk/CVonline/LOCAL_COPIES/OWENS/LECT3/node2.html + - http://www.dai.ed.ac.uk/HIPR2/label.htm + +HOW TO USE IT: + mincmorph -help for more information
new file mode 100644 --- /dev/null +++ b/progs/mincmorph/kernel_io.c @@ -0,0 +1,409 @@ +/* kernel_io.c - reads kernel files */ + +#include <volume_io.h> +#include "kernel_io.h" +#define MAX_KERNEL_ELEMS 1000 + +extern int verbose; + +static const STRING KERNEL_FILE_HEADER = "MNI Morphology Kernel File"; +static const STRING KERNEL_TYPE = "Kernel_Type"; +static const STRING NORMAL_KERNEL = "Normal_Kernel"; +static const STRING KERNEL = "Kernel"; + +/* inbuilt kernels */ +int n_inbuilt_kern = 4; + +/* returns a new Kernel struct (pointer) */ +Kernel *new_kernel(int nelems) +{ + int i, j; + Kernel *tmp; + + ALLOC_VAR_SIZED_STRUCT(tmp, Real, 10); + + /* allocate and initialise the Kernel Array */ + SET_ARRAY_SIZE(tmp->K, 0, nelems, 10); + for(i = 0; i < nelems; i++){ + ALLOC(tmp->K[i], KERNEL_DIMS + 1); + + for(j = 0; j < KERNEL_DIMS; j++){ + tmp->K[i][j] = 0.0; + } + tmp->K[i][KERNEL_DIMS] = 1.0; + + } + tmp->nelems = nelems; + + return tmp; + } + +/* reads in a Kernel from a file */ +Status input_kernel(const char *kernel_file, Kernel * kernel) +{ + int i, j; + + STRING line; + STRING type_name; + STRING str; + Real tmp_real; + FILE *file; + + /* parameter checking */ + if(kernel_file == NULL){ + print_error("input_kernel(): passed NULL FILE.\n"); + return (ERROR); + } + + file = fopen(kernel_file, "r"); + if(file == NULL){ + print_error("input_kernel(): error opening Kernel file.\n"); + return (ERROR); + } + + /* okay read the header */ + if(mni_input_string(file, &line, (char)0, (char)0) != OK){ + delete_string(line); + print_error("input_kernel(): could not read header in file.\n"); + return (ERROR); + } + + if(!equal_strings(line, KERNEL_FILE_HEADER)){ + delete_string(line); + print_error("input_kernel(): invalid header in file.\n"); + return (ERROR); + } + + /* --- read the type of Kernel */ + if(mni_input_keyword_and_equal_sign(file, KERNEL_TYPE, FALSE) != OK){ + return (ERROR); + } + + if(mni_input_string(file, &type_name, (char)';', (char)0) != OK){ + print_error("input_kernel(): missing kernel type.\n"); + return (ERROR); + } + + if(mni_skip_expected_character(file, (char)';') != OK){ + return (ERROR); + } + + if(!equal_strings(type_name, NORMAL_KERNEL)){ + print_error("input_kernel(): invalid kernel type.\n"); + delete_string(type_name); + return (ERROR); + } + delete_string(type_name); + + /* --- read the next string */ + if(mni_input_string(file, &str, (char)'=', (char)0) != OK) + return (ERROR); + + if(!equal_strings(str, KERNEL)){ + print_error("Expected %s =\n", KERNEL); + delete_string(str); + return (ERROR); + } + delete_string(str); + + if(mni_skip_expected_character(file, (char)'=') != OK){ + return (ERROR); + } + + /* now read the elements (lines) of the kernel */ + if(verbose){ + fprintf(stderr, "Reading [%s]", kernel_file); + } + for(i = 0; i < MAX_KERNEL_ELEMS; i++){ + + /* allocate a bit of memory */ + SET_ARRAY_SIZE(kernel->K, kernel->nelems, kernel->nelems + 1, 10); + ALLOC(kernel->K[i], KERNEL_DIMS + 1); + + /* get the 5 dimension vectors and the coefficient */ + for(j = 0; j < 6; j++){ + if(mni_input_real(file, &tmp_real) == OK){ + kernel->K[i][j] = tmp_real; + } + else { + /* check for end */ + if(mni_skip_expected_character(file, (char)';') == OK){ + kernel->nelems = i; + if(verbose){ + fprintf(stderr, " %dx%d Kernel elements read\n", i, kernel->nelems); + } + return (OK); + } + else { + print_error("input_kernel(): error reading kernel [%d,%d]\n", i + 1, + j + 1); + return (ERROR); + } + } + } + kernel->nelems++; + + if(verbose){ + fprintf(stderr, "."); + fflush(stderr); + } + } + + /* SHOLDN'T BE REACHED */ + print_error("input_kernel(): Glark! Something is amiss in the State of Kansas\n"); + return (ERROR); + } + +/* pretty print a kernel */ +int print_kernel(Kernel * kernel) +{ + int i, j; + + fprintf(stderr, " x y z t v coeff\n"); + fprintf(stderr, " -----------------------------------------------\n"); + for(i = 0; i < kernel->nelems; i++){ + fprintf(stderr, "[%02d]", i); + for(j = 0; j < KERNEL_DIMS + 1; j++){ + fprintf(stderr, "%8.02f", kernel->K[i][j]); + } + fprintf(stderr, "\n"); + } + + return (TRUE); + } + +int setup_pad_values(Kernel * kernel) +{ + int c, n; + + /* init padding values */ + for(n = 0; n < KERNEL_DIMS; n++){ + kernel->pre_pad[n] = 0; + kernel->post_pad[n] = 0; + } + + /* find the padding sizes */ + for(c = 0; c < kernel->nelems; c++){ + for(n = 0; n < KERNEL_DIMS; n++){ + if(kernel->K[c][n] < kernel->pre_pad[n]){ + kernel->pre_pad[n] = kernel->K[c][n]; + } + + if(kernel->K[c][n] > kernel->post_pad[n]){ + kernel->post_pad[n] = kernel->K[c][n]; + } + } + } + + return (TRUE); + } + +/* 2D 4 connectivity kernel */ +/* x y z t v coeff */ +/* ----------------------------------------------- */ +/* [00] 1.00 0.00 0.00 0.00 0.00 1.00 */ +/* [01] -1.00 0.00 0.00 0.00 0.00 1.00 */ +/* [02] 0.00 1.00 0.00 0.00 0.00 1.00 */ +/* [03] 0.00 -1.00 0.00 0.00 0.00 1.00 */ +Kernel *get_2D04_kernel(void) +{ + Kernel *K = new_kernel(4); + + K->K[0][0] = 1.0; + K->K[1][0] = -1.0; + + K->K[2][1] = 1.0; + K->K[3][1] = -1.0; + + return K; + } + +/* 2D 8 connectivity kernel */ +/* x y z t v coeff */ +/* ----------------------------------------------- */ +/* [00] 1.00 1.00 0.00 0.00 0.00 1.00 */ +/* [01] 1.00 0.00 0.00 0.00 0.00 1.00 */ +/* [02] 1.00 -1.00 0.00 0.00 0.00 1.00 */ +/* */ +/* [03] 0.00 1.00 0.00 0.00 0.00 1.00 */ +/* [04] 0.00 -1.00 0.00 0.00 0.00 1.00 */ +/* */ +/* [05] -1.00 1.00 0.00 0.00 0.00 1.00 */ +/* [06] -1.00 0.00 0.00 0.00 0.00 1.00 */ +/* [07] -1.00 -1.00 0.00 0.00 0.00 1.00 */ +Kernel *get_2D08_kernel(void) +{ + Kernel *K = new_kernel(8); + + K->K[0][0] = 1.0; + K->K[0][1] = 1.0; + + K->K[1][0] = 1.0; + + K->K[2][0] = 1.0; + K->K[2][1] = -1.0; + + K->K[3][1] = 1.0; + + K->K[4][1] = -1.0; + + K->K[5][0] = -1.0; + K->K[5][1] = 1.0; + + K->K[6][0] = -1.0; + + K->K[7][0] = -1.0; + K->K[7][1] = -1.0; + + return K; + } + +/* 3D 6 connectivity kernel */ +/* x y z t v coeff */ +/* ----------------------------------------------- */ +/* [00] 1.00 0.00 0.00 0.00 0.00 1.00 */ +/* [01] -1.00 0.00 0.00 0.00 0.00 1.00 */ +/* [02] 0.00 1.00 0.00 0.00 0.00 1.00 */ +/* [03] 0.00 -1.00 0.00 0.00 0.00 1.00 */ +/* [04] 0.00 0.00 1.00 0.00 0.00 1.00 */ +/* [05] 0.00 0.00 -1.00 0.00 0.00 1.00 */ +Kernel *get_3D06_kernel(void) +{ + Kernel *K = new_kernel(6); + + K->K[0][0] = 1.0; + K->K[1][0] = -1.0; + + K->K[2][1] = 1.0; + K->K[3][1] = -1.0; + + K->K[4][2] = 1.0; + K->K[5][2] = -1.0; + + return K; + } + +/* 3D 26 connectivity kernel */ +/* x y z t v coeff */ +/* ----------------------------------------------- */ +/* [00] 1.00 1.00 1.00 0.00 0.00 1.00 */ +/* [01] 1.00 1.00 0.00 0.00 0.00 1.00 */ +/* [02] 1.00 1.00 -1.00 0.00 0.00 1.00 */ + +/* [03] 1.00 0.00 1.00 0.00 0.00 1.00 */ +/* [04] 1.00 0.00 0.00 0.00 0.00 1.00 */ +/* [05] 1.00 0.00 -1.00 0.00 0.00 1.00 */ + +/* [06] 1.00 -1.00 1.00 0.00 0.00 1.00 */ +/* [07] 1.00 -1.00 0.00 0.00 0.00 1.00 */ +/* [08] 1.00 -1.00 -1.00 0.00 0.00 1.00 */ + +/* [09] 0.00 1.00 1.00 0.00 0.00 1.00 */ +/* [10] 0.00 1.00 0.00 0.00 0.00 1.00 */ +/* [11] 0.00 1.00 -1.00 0.00 0.00 1.00 */ + +/* [12] 0.00 0.00 1.00 0.00 0.00 1.00 */ +/* ---- 0.00 0.00 0.00 0.00 0.00 ---- */ +/* [13] 0.00 0.00 -1.00 0.00 0.00 1.00 */ + +/* [14] 0.00 -1.00 1.00 0.00 0.00 1.00 */ +/* [15] 0.00 -1.00 0.00 0.00 0.00 1.00 */ +/* [16] 0.00 -1.00 -1.00 0.00 0.00 1.00 */ + +/* [17] -1.00 1.00 1.00 0.00 0.00 1.00 */ +/* [18] -1.00 1.00 0.00 0.00 0.00 1.00 */ +/* [19] -1.00 1.00 -1.00 0.00 0.00 1.00 */ + +/* [20] -1.00 0.00 1.00 0.00 0.00 1.00 */ +/* [21] -1.00 0.00 0.00 0.00 0.00 1.00 */ +/* [22] -1.00 0.00 -1.00 0.00 0.00 1.00 */ + +/* [23] -1.00 -1.00 1.00 0.00 0.00 1.00 */ +/* [24] -1.00 -1.00 0.00 0.00 0.00 1.00 */ +/* [25] -1.00 -1.00 -1.00 0.00 0.00 1.00 */ +Kernel *get_3D26_kernel(void) +{ + Kernel *K = new_kernel(26); + + K->K[0][0] = 1.0; + K->K[0][1] = 1.0; + K->K[0][2] = 1.0; + + K->K[1][0] = 1.0; + K->K[1][1] = 1.0; + + K->K[2][0] = 1.0; + K->K[2][1] = 1.0; + K->K[2][2] = -1.0; + + K->K[3][0] = 1.0; + K->K[3][2] = 1.0; + + K->K[4][0] = 1.0; + + K->K[5][0] = 1.0; + K->K[5][2] = -1.0; + + K->K[6][0] = 1.0; + K->K[6][1] = -1.0; + K->K[6][2] = 1.0; + + K->K[7][0] = 1.0; + K->K[7][1] = -1.0; + + K->K[8][0] = 1.0; + K->K[8][1] = -1.0; + K->K[8][2] = -1.0; + + K->K[9][1] = 1.0; + K->K[9][2] = 1.0; + + K->K[10][1] = 1.0; + + K->K[11][1] = 1.0; + K->K[11][2] = -1.0; + + K->K[12][2] = 1.0; + + K->K[13][2] = -1.0; + + K->K[14][1] = -1.0; + K->K[14][2] = 1.0; + + K->K[15][1] = -1.0; + + K->K[16][1] = -1.0; + K->K[16][2] = -1.0; + + K->K[17][0] = -1.0; + K->K[17][1] = 1.0; + K->K[17][2] = 1.0; + + K->K[18][0] = -1.0; + K->K[18][1] = 1.0; + + K->K[19][0] = -1.0; + K->K[19][1] = 1.0; + K->K[19][2] = -1.0; + + K->K[20][0] = -1.0; + K->K[20][2] = 1.0; + + K->K[21][0] = -1.0; + + K->K[22][0] = -1.0; + K->K[22][2] = -1.0; + + K->K[23][0] = -1.0; + K->K[23][1] = -1.0; + K->K[23][2] = 1.0; + + K->K[24][0] = -1.0; + K->K[24][1] = -1.0; + + K->K[25][0] = -1.0; + K->K[25][1] = -1.0; + K->K[25][2] = -1.0; + + return K; + }
new file mode 100644 --- /dev/null +++ b/progs/mincmorph/kernel_io.h @@ -0,0 +1,42 @@ +/* kernel_io.h */ + +#ifndef KERNEL_IO +#define KERNEL_IO + +#define KERNEL_DIMS 5 + +/* inbuilt kernels */ +int n_inbuilt_kern; + +typedef enum { + K_NULL = 0, + K_2D04, K_2D08, K_3D06, K_3D26 + } kern_types; + +/* Structure for Kernel information */ +typedef struct { + int nelems; + int pre_pad[KERNEL_DIMS]; + int post_pad[KERNEL_DIMS]; + VIO_Real **K; + } Kernel; + +/* returns a new B_Matrix struct (pointer) */ +Kernel *new_kernel(int nelems); + +/* reads in a B_Matrix from a file (pointer) */ +VIO_Status input_kernel(const char *kernel_file, Kernel * kernel); + +/* pretty print a kernel */ +int print_kernel(Kernel * kernel); + +/* calculate start and step offsets for this kernel */ +int setup_pad_values(Kernel * kernel); + +/* return the default kernel(s) */ +Kernel *get_2D04_kernel(void); +Kernel *get_2D08_kernel(void); +Kernel *get_3D06_kernel(void); +Kernel *get_3D26_kernel(void); + +#endif
new file mode 100644 --- /dev/null +++ b/progs/mincmorph/kernel_ops.c @@ -0,0 +1,797 @@ +/* kernel_ops.c */ + +#include <float.h> +#include <limits.h> +#include "kernel_ops.h" + +extern int verbose; + +/* function prototypes */ +void split_kernel(Kernel * K, Kernel * k1, Kernel * k2); +int compare_ints(const void *a, const void *b); +int compare_groups(const void *a, const void *b); + +/* structure for group information */ +typedef struct { + unsigned int orig_label; + unsigned int count; + } group_info_struct; + +typedef group_info_struct *Group_info; + +int compare_ints(const void *a, const void *b) +{ + return (*(int *)a - *(int *)b); + } + +int compare_groups(const void *a, const void *b) +{ + return (*(Group_info *) b)->count - (*(Group_info *) a)->count; + } + +void split_kernel(Kernel * K, Kernel * k1, Kernel * k2) +{ + int c, k1c, k2c; + + /* fill the two sub kernels */ + k1c = k2c = 0; + for(c = 0; c < K->nelems; c++){ + if((K->K[c][2] < 0) || + (K->K[c][1] < 0 && K->K[c][2] <= 0) || + (K->K[c][0] < 0 && K->K[c][1] <= 0 && K->K[c][2] <= 0)){ + + k1->K[k1c] = K->K[c]; + k1c++; + } + else { + k2->K[k2c] = K->K[c]; + k2c++; + } + } + k1->nelems = k1c; + k2->nelems = k2c; + } + +/* binarise a volume between a range */ +Volume *binarise(Volume * vol, double floor, double ceil, double fg, double bg) +{ + int x, y, z; + int sizes[MAX_VAR_DIMS]; + double value; + progress_struct progress; + + if(verbose){ + fprintf(stdout, "Binarising, range: [%g:%g] fg/bg: [%g:%g]\n", floor, ceil, fg, bg); + } + get_volume_sizes(*vol, sizes); + + initialize_progress_report(&progress, FALSE, sizes[2], "Binarise"); + for(z = sizes[0]; z--;){ + for(y = sizes[1]; y--;){ + for(x = sizes[2]; x--;){ + + value = get_volume_voxel_value(*vol, z, y, x, 0, 0); + if((value >= floor) && (value <= ceil)){ + set_volume_voxel_value(*vol, z, y, x, 0, 0, fg); + } + else { + set_volume_voxel_value(*vol, z, y, x, 0, 0, bg); + } + + } + } + update_progress_report(&progress, z + 1); + } + terminate_progress_report(&progress); + + return (vol); + } + +/* clamp a volume between a range */ +Volume *clamp(Volume * vol, double floor, double ceil, double bg) +{ + int x, y, z; + int sizes[MAX_VAR_DIMS]; + double value; + progress_struct progress; + + if(verbose){ + fprintf(stdout, "Clamping, range: [%g:%g] bg: %g\n", floor, ceil, bg); + } + + get_volume_sizes(*vol, sizes); + + initialize_progress_report(&progress, FALSE, sizes[2], "Clamping"); + for(z = sizes[0]; z--;){ + for(y = sizes[1]; y--;){ + for(x = sizes[2]; x--;){ + value = get_volume_voxel_value(*vol, z, y, x, 0, 0); + if((value < floor) || (value > ceil)){ + set_volume_voxel_value(*vol, z, y, x, 0, 0, bg); + } + } + } + update_progress_report(&progress, z + 1); + } + + terminate_progress_report(&progress); + return (vol); + } + +/* pad a volume using the background value */ +Volume *pad(Kernel * K, Volume * vol, double bg) +{ + int x, y, z; + int sizes[MAX_VAR_DIMS]; + + get_volume_sizes(*vol, sizes); + + /* z */ + for(y = 0; y < sizes[1]; y++){ + for(x = 0; x < sizes[2]; x++){ + for(z = 0; z < -K->pre_pad[2]; z++){ + set_volume_real_value(*vol, z, y, x, 0, 0, bg); + } + for(z = sizes[0] - K->post_pad[2]; z < sizes[0]; z++){ + set_volume_real_value(*vol, z, y, x, 0, 0, bg); + } + } + } + + /* y */ + for(z = 0; z < sizes[0]; z++){ + for(x = 0; x < sizes[2]; x++){ + for(y = 0; y < -K->pre_pad[1]; y++){ + set_volume_real_value(*vol, z, y, x, 0, 0, bg); + } + for(y = sizes[1] - K->post_pad[1]; y < sizes[1]; y++){ + set_volume_real_value(*vol, z, y, x, 0, 0, bg); + } + } + } + + /* x */ + for(z = 0; z < sizes[0]; z++){ + for(y = 0; y < sizes[1]; y++){ + for(x = 0; x < -K->pre_pad[0]; x++){ + set_volume_real_value(*vol, z, y, x, 0, 0, bg); + } + for(x = sizes[2] - K->post_pad[0]; x < sizes[2]; x++){ + set_volume_real_value(*vol, z, y, x, 0, 0, bg); + } + } + } + + return (vol); + } + +/* perform a dilation on a volume */ +Volume *dilation_kernel(Kernel * K, Volume * vol) +{ + int x, y, z, c; + double value; + int sizes[MAX_VAR_DIMS]; + progress_struct progress; + Volume tmp_vol; + + if(verbose){ + fprintf(stdout, "Dilation kernel\n"); + } + get_volume_sizes(*vol, sizes); + initialize_progress_report(&progress, FALSE, sizes[2], "Dilation"); + + /* copy the volume */ + tmp_vol = copy_volume(*vol); + + for(z = -K->pre_pad[2]; z < sizes[0] - K->post_pad[2]; z++){ + for(y = -K->pre_pad[1]; y < sizes[1] - K->post_pad[1]; y++){ + for(x = -K->pre_pad[0]; x < sizes[2] - K->post_pad[0]; x++){ + + value = get_volume_real_value(tmp_vol, z, y, x, 0, 0); + for(c = 0; c < K->nelems; c++){ + if(get_volume_real_value(*vol, + z + K->K[c][2], + y + K->K[c][1], + x + K->K[c][0], + 0 + K->K[c][3], 0 + K->K[c][4]) < value){ + set_volume_real_value(*vol, + z + K->K[c][2], + y + K->K[c][1], + x + K->K[c][0], + 0 + K->K[c][3], + 0 + K->K[c][4], value * K->K[c][5]); + } + } + } + } + + update_progress_report(&progress, z + 1); + } + + delete_volume(tmp_vol); + terminate_progress_report(&progress); + return (vol); + } + +/* perform a median kernel operation on a volume */ +Volume *median_dilation_kernel(Kernel * K, Volume * vol) +{ + int x, y, z, c, i; + int sizes[MAX_VAR_DIMS]; + progress_struct progress; + Volume tmp_vol; + double value; + + unsigned int kvalue; + unsigned int neighbours[K->nelems]; + + if(verbose){ + fprintf(stdout, "Median Dilation kernel\n"); + } + get_volume_sizes(*vol, sizes); + initialize_progress_report(&progress, FALSE, sizes[2], "Median Dilation"); + + /* copy the volume */ + tmp_vol = copy_volume(*vol); + + for(z = -K->pre_pad[2]; z < sizes[0] - K->post_pad[2]; z++){ + for(y = -K->pre_pad[1]; y < sizes[1] - K->post_pad[1]; y++){ + for(x = -K->pre_pad[0]; x < sizes[2] - K->post_pad[0]; x++){ + + /* only modify background voxels */ + value = get_volume_voxel_value(tmp_vol, z, y, x, 0, 0); + if(value == 0.0){ + + i = 0; + for(c = 0; c < K->nelems; c++){ + + kvalue = (unsigned int)get_volume_voxel_value(tmp_vol, + z + K->K[c][2], + y + K->K[c][1], + x + K->K[c][0], + 0 + K->K[c][3], + 0 + K->K[c][4]); + if(kvalue != 0){ + neighbours[i] = kvalue; + i++; + } + } + + /* only run this for adjacent voxels */ + if(i > 0){ + + /* find the median of our little array */ + qsort(&neighbours[0], (size_t) i, sizeof(unsigned int), &compare_ints); + + /* store the median value */ + set_volume_voxel_value(*vol, z, y, x, 0, 0, + (double)neighbours[(int)floor((i - 1) / 2)]); + } + } + + /* else just copy the original value over */ + else { + set_volume_voxel_value(*vol, z, y, x, 0, 0, value); + } + } + } + + update_progress_report(&progress, z + 1); + } + + delete_volume(tmp_vol); + terminate_progress_report(&progress); + return (vol); + } + +/* perform an erosion on a volume */ +Volume *erosion_kernel(Kernel * K, Volume * vol) +{ + int x, y, z, c; + double value; + int sizes[MAX_VAR_DIMS]; + progress_struct progress; + Volume tmp_vol; + + if(verbose){ + fprintf(stdout, "Erosion kernel\n"); + } + get_volume_sizes(*vol, sizes); + initialize_progress_report(&progress, FALSE, sizes[2], "Erosion"); + + /* copy the volume */ + tmp_vol = copy_volume(*vol); + + for(z = -K->pre_pad[2]; z < sizes[0] - K->post_pad[2]; z++){ + for(y = -K->pre_pad[1]; y < sizes[1] - K->post_pad[1]; y++){ + for(x = -K->pre_pad[0]; x < sizes[2] - K->post_pad[0]; x++){ + + value = get_volume_real_value(tmp_vol, z, y, x, 0, 0); + for(c = 0; c < K->nelems; c++){ + if(get_volume_real_value(*vol, + z + K->K[c][2], + y + K->K[c][1], + x + K->K[c][0], + 0 + K->K[c][3], 0 + K->K[c][4]) > value){ + set_volume_real_value(*vol, + z + K->K[c][2], + y + K->K[c][1], + x + K->K[c][0], + 0 + K->K[c][3], + 0 + K->K[c][4], value * K->K[c][5]); + } + } + + value = get_volume_real_value(tmp_vol, z, y, x, 0, 0); + } + } + update_progress_report(&progress, z + 1); + } + + delete_volume(tmp_vol); + terminate_progress_report(&progress); + return (vol); + } + +/* convolve a volume with a input kernel */ +Volume *convolve_kernel(Kernel * K, Volume * vol) +{ + int x, y, z, c; + double value; + int sizes[MAX_VAR_DIMS]; + progress_struct progress; + Volume tmp_vol; + + if(verbose){ + fprintf(stdout, "Convolve kernel\n"); + } + get_volume_sizes(*vol, sizes); + initialize_progress_report(&progress, FALSE, sizes[2], "Convolve"); + + /* copy the volume */ + tmp_vol = copy_volume(*vol); + + for(z = -K->pre_pad[2]; z < sizes[0] - K->post_pad[2]; z++){ + for(y = -K->pre_pad[1]; y < sizes[1] - K->post_pad[1]; y++){ + for(x = -K->pre_pad[0]; x < sizes[2] - K->post_pad[0]; x++){ + + value = 0; + for(c = 0; c < K->nelems; c++){ + value += get_volume_real_value(tmp_vol, + z + K->K[c][2], + y + K->K[c][1], + x + K->K[c][0], 0 + K->K[c][3], + 0 + K->K[c][4]) * K->K[c][5]; + } + set_volume_real_value(*vol, z, y, x, 0, 0, value); + } + } + + update_progress_report(&progress, z + 1); + } + + delete_volume(tmp_vol); + terminate_progress_report(&progress); + return (vol); + } + +/* should really only work on binary images */ +/* from the original 2 pass Borgefors alg */ +Volume *distance_kernel(Kernel * K, Volume * vol, double bg) +{ + int x, y, z, c; + double value, min; + int sizes[MAX_VAR_DIMS]; + progress_struct progress; + Kernel *k1, *k2; + + /* split the Kernel */ + k1 = new_kernel(K->nelems); + k2 = new_kernel(K->nelems); + split_kernel(K, k1, k2); + + setup_pad_values(k1); + setup_pad_values(k2); + + if(verbose){ + fprintf(stdout, "Distance kernel - background %g\n", bg); + fprintf(stdout, "forward direction kernel:\n"); + print_kernel(k1); + fprintf(stdout, "\nreverse direction kernel:\n"); + print_kernel(k2); + } + + get_volume_sizes(*vol, sizes); + initialize_progress_report(&progress, FALSE, sizes[2] * 2, "Distance"); + + /* forward raster direction */ + for(z = -K->pre_pad[2]; z < sizes[0] - K->post_pad[2]; z++){ + for(y = -K->pre_pad[1]; y < sizes[1] - K->post_pad[1]; y++){ + for(x = -K->pre_pad[0]; x < sizes[2] - K->post_pad[0]; x++){ + + if(get_volume_real_value(*vol, z, y, x, 0, 0) != bg){ + + /* find the minimum */ + min = DBL_MAX; + for(c = 0; c < k1->nelems; c++){ + value = get_volume_real_value(*vol, + z + k1->K[c][2], + y + k1->K[c][1], + x + k1->K[c][0], + 0 + k1->K[c][3], 0 + k1->K[c][4]) + 1; + if(value < min){ + min = value; + } + } + + set_volume_real_value(*vol, z, y, x, 0, 0, min); + } + } + } + update_progress_report(&progress, z + 1); + } + + /* reverse raster direction */ + for(z = sizes[0] - k2->post_pad[2] - 1; z >= -k2->pre_pad[2]; z--){ + for(y = sizes[1] - k2->post_pad[1] - 1; y >= -k2->pre_pad[1]; y--){ + for(x = sizes[2] - k2->post_pad[0] - 1; x >= -k2->pre_pad[0]; x--){ + + min = get_volume_real_value(*vol, z, y, x, 0, 0); + if(min != bg){ + + /* find the minimum distance to bg in the neighbouring vectors */ + for(c = 0; c < k2->nelems; c++){ + value = get_volume_real_value(*vol, + z + k2->K[c][2], + y + k2->K[c][1], + x + k2->K[c][0], + 0 + k2->K[c][3], 0 + k2->K[c][4]) + 1; + if(value < min){ + min = value; + } + } + + set_volume_real_value(*vol, z, y, x, 0, 0, min); + } + } + } + update_progress_report(&progress, sizes[2] + z + 1); + } + + free(k1); + free(k2); + terminate_progress_report(&progress); + return (vol); + } + +/* do connected components labelling on a volume */ +/* resulting groups are sorted WRT size */ +Volume *group_kernel(Kernel * K, Volume * vol, double bg) +{ + int x, y, z; + int sizes[MAX_VAR_DIMS]; + progress_struct progress; + Volume tmp_vol; + Kernel *k1, *k2; + + unsigned int *equiv; + unsigned int *counts; + unsigned int *trans; + unsigned int neighbours[K->nelems]; + + /* counters */ + unsigned int c; + unsigned int value; + unsigned int group_idx; /* label for the next group */ + unsigned int num_groups; + + unsigned int min_label; + unsigned int curr_label; + unsigned int prev_label; + unsigned int num_matches; + + /* structure for group data */ + Group_info *group_data; + + /* split the Kernel into forward and backwards kernels */ + k1 = new_kernel(K->nelems); + k2 = new_kernel(K->nelems); + split_kernel(K, k1, k2); + + setup_pad_values(k1); + setup_pad_values(k2); + + if(verbose){ + fprintf(stdout, "Group kernel - background %g\n", bg); + fprintf(stdout, "forward direction kernel:\n"); + print_kernel(k1); + fprintf(stdout, "\nreverse direction kernel:\n"); + print_kernel(k2); + } + + get_volume_sizes(*vol, sizes); + initialize_progress_report(&progress, FALSE, sizes[2], "Groups"); + + /* copy and then zero out the original volume */ + tmp_vol = copy_volume(*vol); + for(z = sizes[0]; z--;){ + for(y = sizes[1]; y--;){ + for(x = sizes[2]; x--;){ + set_volume_voxel_value(*vol, z, y, x, 0, 0, 0); + } + } + } + + /* pass 1 - forward direction (we assume a symmetric kernel) */ + + /* our first group is given the label 1 */ + group_idx = 1; + + /* initialise the equiv and counts arrays */ + SET_ARRAY_SIZE(equiv, 0, group_idx, 500); + equiv[0] = 0; + + SET_ARRAY_SIZE(counts, 0, group_idx, 500); + counts[0] = 0; + + for(z = -k1->pre_pad[2]; z < sizes[0] - k1->post_pad[2]; z++){ + for(y = -k1->pre_pad[1]; y < sizes[1] - k1->post_pad[1]; y++){ + for(x = -k1->pre_pad[0]; x < sizes[2] - k1->post_pad[0]; x++){ + + if(get_volume_voxel_value(tmp_vol, z, y, x, 0, 0) != bg){ + + /* search this voxels neighbours */ + num_matches = 0; + min_label = INT_MAX; + + for(c = 0; c < k1->nelems; c++){ + value = (unsigned int)get_volume_voxel_value(*vol, + z + k1->K[c][2], + y + k1->K[c][1], + x + k1->K[c][0], + 0 + k1->K[c][3], + 0 + k1->K[c][4]); + if(value != 0){ + if(value < min_label){ + min_label = value; + } + neighbours[num_matches] = value; + num_matches++; + } + } + + switch (num_matches){ + case 0: + /* no neighbours, make a new label and increment */ + set_volume_voxel_value(*vol, z, y, x, 0, 0, (Real) group_idx); + + SET_ARRAY_SIZE(equiv, group_idx, group_idx + 1, 500); + equiv[group_idx] = group_idx; + + SET_ARRAY_SIZE(counts, group_idx, group_idx + 1, 500); + counts[group_idx] = 1; + + group_idx++; + break; + + case 1: + /* only one neighbour, no equivalences needed */ + set_volume_voxel_value(*vol, z, y, x, 0, 0, (Real) min_label); + counts[min_label]++; + break; + + default: + /* more than one neighbour */ + + /* first sort the neighbours array */ + qsort(&neighbours[0], (size_t) num_matches, sizeof(unsigned int), + &compare_ints); + + /* find the minimum possible label for this voxel, */ + /* this is done by descending through each neighbours */ + /* equivalences until an equivalence equal to itself */ + /* is found */ + prev_label = -1; + for(c = 0; c < num_matches; c++){ + curr_label = neighbours[c]; + + /* recurse this label if we haven't yet */ + if(curr_label != prev_label){ + while(equiv[curr_label] != equiv[equiv[curr_label]]){ + curr_label = equiv[curr_label]; + } + + /* check against the current minimum value */ + if(equiv[curr_label] < min_label){ + min_label = equiv[curr_label]; + } + } + + prev_label = neighbours[c]; + } + + /* repeat, setting equivalences to the min_label */ + prev_label = -1; + for(c = 0; c < num_matches; c++){ + curr_label = neighbours[c]; + + if(curr_label != prev_label){ + while(equiv[curr_label] != equiv[equiv[curr_label]]){ + curr_label = equiv[curr_label]; + + equiv[curr_label] = min_label; + } + + /* set the label itself */ + if(equiv[neighbours[c]] != min_label){ + equiv[neighbours[c]] = min_label; + } + } + + prev_label = neighbours[c]; + } + + /* finally set the voxel in question to the minimum value */ + set_volume_voxel_value(*vol, z, y, x, 0, 0, (Real) min_label); + counts[min_label]++; + break; + } /* end case */ + + } + } + } + update_progress_report(&progress, z + 1); + } + terminate_progress_report(&progress); + + /* reduce the equiv and counts array */ + num_groups = 0; + for(c = 0; c < group_idx; c++){ + + /* if this equivalence is not resolved yet */ + if(c != equiv[c]){ + + /* find the min label value */ + min_label = equiv[c]; + while(min_label != equiv[min_label]){ + min_label = equiv[min_label]; + } + + /* update the label and its counters */ + equiv[c] = min_label; + counts[min_label] += counts[c]; + counts[c] = 0; + } + else { + num_groups++; + } + } + + /* Allocate space for the array of groups */ + group_data = (Group_info *) malloc(num_groups * sizeof(Group_info)); + + num_groups = 0; + for(c = 0; c < group_idx; c++){ + if(counts[c] > 0){ + /* allocate space for this element */ + group_data[num_groups] = malloc(sizeof(group_info_struct)); + + group_data[num_groups]->orig_label = equiv[c]; + group_data[num_groups]->count = counts[c]; + num_groups++; + } + } + + /* sort the groups by the count size */ + if(verbose){ + fprintf(stdout, "Found %d unique groups from %d, sorting...\n", num_groups, + group_idx); + } + qsort(group_data, num_groups, sizeof(Group_info), &compare_groups); + + /* set up the transpose array */ + trans = (unsigned int *)malloc(sizeof(unsigned int) * group_idx); + for(c = 0; c < num_groups; c++){ + trans[group_data[c]->orig_label] = c + 1; /* +1 to bump past 0 */ + } + + /* pass 2 - resolve equivalences in the output data */ + if(verbose){ + fprintf(stdout, "Resolving equivalences...\n"); + } + for(z = sizes[0]; z--;){ + for(y = sizes[1]; y--;){ + for(x = sizes[2]; x--;){ + value = (unsigned int)get_volume_voxel_value(*vol, z, y, x, 0, 0); + if(value != 0){ + value = trans[equiv[value]]; + set_volume_voxel_value(*vol, z, y, x, 0, 0, (Real) value); + } + } + } + } + + /* tidy up */ + delete_volume(tmp_vol); + for(c = 0; c < num_groups; c++){ + free(group_data[c]); + } + free(group_data); + free(trans); + free(k1); + free(k2); + + return (vol); + } + +/* do local correlation to another volume */ +/* xcorr = sum((a*b)^2) / (sqrt(sum(a^2)) * sqrt(sum(b^2)) */ +VIO_Volume *lcorr_kernel(Kernel * K, VIO_Volume * vol, VIO_Volume *cmp) +{ + int x, y, z, c; + double value, v1, v2; + double ssum_v1, ssum_v2, sum_prd, denom; + int sizes[MAX_VAR_DIMS]; + progress_struct progress; + Volume tmp_vol; + + if(verbose){ + fprintf(stdout, "Local Correlation kernel\n"); + } + get_volume_sizes(*vol, sizes); + initialize_progress_report(&progress, FALSE, sizes[2], "Local Correlation"); + + /* copy the volume */ + tmp_vol = copy_volume(*vol); + + /* zero the output volume */ + for(z = sizes[0]; z--;){ + for(y = sizes[1]; y--;){ + for(x = sizes[2]; x--;){ + set_volume_voxel_value(*vol, z, y, x, 0, 0, 0); + } + } + } + + /* set output range */ + set_volume_real_range(*vol, 0.0, 1.0); + + for(z = -K->pre_pad[2]; z < sizes[0] - K->post_pad[2]; z++){ + for(y = -K->pre_pad[1]; y < sizes[1] - K->post_pad[1]; y++){ + for(x = -K->pre_pad[0]; x < sizes[2] - K->post_pad[0]; x++){ + + /* init counters */ + ssum_v1 = ssum_v2 = sum_prd = 0; + for(c = 0; c < K->nelems; c++){ + v1 = get_volume_real_value(tmp_vol, + z + K->K[c][2], + y + K->K[c][1], + x + K->K[c][0], 0 + K->K[c][3], + 0 + K->K[c][4]) * K->K[c][5]; + v2 = get_volume_real_value(*cmp, + z + K->K[c][2], + y + K->K[c][1], + x + K->K[c][0], 0 + K->K[c][3], + 0 + K->K[c][4]) * K->K[c][5]; + + /* increment counters */ + ssum_v1 += v1*v1; + ssum_v2 += v2*v2; + sum_prd += v1*v2; + } + + denom = sqrt(ssum_v1 * ssum_v2); + value = (denom == 0.0) ? 0.0 : sum_prd / denom; + + set_volume_real_value(*vol, z, y, x, 0, 0, value); + } + } + update_progress_report(&progress, z + 1); + } + terminate_progress_report(&progress); + + /* tidy up */ + delete_volume(tmp_vol); + + return (vol); + }
new file mode 100644 --- /dev/null +++ b/progs/mincmorph/kernel_ops.h @@ -0,0 +1,21 @@ +/* kernel_ops.h */ + +#ifndef KERNEL_OPS +#define KERNEL_OPS + +#include <volume_io.h> +#include "kernel_io.h" + +/* kernel functions */ +VIO_Volume *binarise(VIO_Volume * vol, double floor, double ceil, double fg, double bg); +VIO_Volume *clamp(VIO_Volume * vol, double floor, double ceil, double bg); +VIO_Volume *pad(Kernel * K, VIO_Volume * vol, double bg); +VIO_Volume *erosion_kernel(Kernel * K, VIO_Volume * vol); +VIO_Volume *dilation_kernel(Kernel * K, VIO_Volume * vol); +VIO_Volume *median_dilation_kernel(Kernel * K, VIO_Volume * vol); +VIO_Volume *convolve_kernel(Kernel * K, VIO_Volume * vol); +VIO_Volume *distance_kernel(Kernel * K, VIO_Volume * vol, double bg); +VIO_Volume *group_kernel(Kernel * K, VIO_Volume * vol, double bg); +VIO_Volume *lcorr_kernel(Kernel * K, VIO_Volume * vol, VIO_Volume *cmp); + +#endif
new file mode 100644 --- /dev/null +++ b/progs/mincmorph/kernels/3x3_4-conn.kern @@ -0,0 +1,20 @@ +MNI Morphology Kernel File +% +% 2D 4-connectivity ie: +% +% 1 +% 1 0 1 +% 1 +% +% 0 - center voxel +% Format: vector for voxel in 5 dimensions (x,y,z,t,v) then multiplier. + +Kernel_Type = Normal_Kernel; +Kernel = +% x y z t v coeff +% ----------------------------------- + 1.0 0.0 0.0 0.0 0.0 1.0 + -1.0 0.0 0.0 0.0 0.0 1.0 + 0.0 1.0 0.0 0.0 0.0 1.0 + 0.0 -1.0 0.0 0.0 0.0 1.0; +
new file mode 100644 --- /dev/null +++ b/progs/mincmorph/kernels/3x3_8-conn.kern @@ -0,0 +1,27 @@ +MNI Morphology Kernel File +% +% 2D 8-connectivity ie: +% +% 1 1 1 +% 1 0 1 +% 1 1 1 +% +% 0 - center voxel +% Format: vector for voxel in 5 dimensions (x,y,z,t,v) then multiplier. + +Kernel_Type = Normal_Kernel; +Kernel = +% x y z t v coeff +% ----------------------------------- + 1.0 1.0 0.0 0.0 0.0 1.0 + 1.0 0.0 0.0 0.0 0.0 1.0 + 1.0 -1.0 0.0 0.0 0.0 1.0 +% + 0.0 1.0 0.0 0.0 0.0 1.0 +%% 0.0 0.0 0.0 0.0 0.0 1.0 + 0.0 -1.0 0.0 0.0 0.0 1.0 +% + -1.0 1.0 0.0 0.0 0.0 1.0 + -1.0 0.0 0.0 0.0 0.0 1.0 + -1.0 -1.0 0.0 0.0 0.0 1.0; +
new file mode 100644 --- /dev/null +++ b/progs/mincmorph/kernels/3x3x3_26-conn.kern @@ -0,0 +1,53 @@ +MNI Morphology Kernel File +% +% 3D 26-connectivity ie: (in 2D) +% +% 1 1 1 +% 1 0 1 +% 1 1 1 +% +% 0 - center voxel +% Format: vector for voxel in 5 dimensions (x,y,z,t,v) then multiplier. + +Kernel_Type = Normal_Kernel; +Kernel = +% x y z t v coeff +% ----------------------------------- + 1.0 1.0 1.0 0.0 0.0 1.0 + 1.0 1.0 0.0 0.0 0.0 1.0 + 1.0 1.0 -1.0 0.0 0.0 1.0 +% + 1.0 0.0 1.0 0.0 0.0 1.0 + 1.0 0.0 0.0 0.0 0.0 1.0 + 1.0 0.0 -1.0 0.0 0.0 1.0 +% + 1.0 -1.0 1.0 0.0 0.0 1.0 + 1.0 -1.0 0.0 0.0 0.0 1.0 + 1.0 -1.0 -1.0 0.0 0.0 1.0 +% +% + 0.0 1.0 1.0 0.0 0.0 1.0 + 0.0 1.0 0.0 0.0 0.0 1.0 + 0.0 1.0 -1.0 0.0 0.0 1.0 +% + 0.0 0.0 1.0 0.0 0.0 1.0 +%% 0.0 0.0 0.0 0.0 0.0 1.0 + 0.0 0.0 -1.0 0.0 0.0 1.0 +% + 0.0 -1.0 1.0 0.0 0.0 1.0 + 0.0 -1.0 0.0 0.0 0.0 1.0 + 0.0 -1.0 -1.0 0.0 0.0 1.0 +% +% + -1.0 1.0 1.0 0.0 0.0 1.0 + -1.0 1.0 0.0 0.0 0.0 1.0 + -1.0 1.0 -1.0 0.0 0.0 1.0 +% + -1.0 0.0 1.0 0.0 0.0 1.0 + -1.0 0.0 0.0 0.0 0.0 1.0 + -1.0 0.0 -1.0 0.0 0.0 1.0 +% + -1.0 -1.0 1.0 0.0 0.0 1.0 + -1.0 -1.0 0.0 0.0 0.0 1.0 + -1.0 -1.0 -1.0 0.0 0.0 1.0; +
new file mode 100644 --- /dev/null +++ b/progs/mincmorph/kernels/3x3x3_6-conn.kern @@ -0,0 +1,22 @@ +MNI Morphology Kernel File +% +% 3D 6-connectivity ie: (in 2D) +% +% 1 +% 1 0 1 +% 1 +% +% 0 - center voxel +% Format: vector for voxel in 5 dimensions (x,y,z,t,v) then multiplier. + +Kernel_Type = Normal_Kernel; +Kernel = +% x y z t v coeff +% ----------------------------------- + 1.0 0.0 0.0 0.0 0.0 1.0 + -1.0 0.0 0.0 0.0 0.0 1.0 + 0.0 1.0 0.0 0.0 0.0 1.0 + 0.0 -1.0 0.0 0.0 0.0 1.0 + 0.0 0.0 1.0 0.0 0.0 1.0 + 0.0 0.0 -1.0 0.0 0.0 1.0; +
new file mode 100644 --- /dev/null +++ b/progs/mincmorph/mincmorph.c @@ -0,0 +1,748 @@ +/* mincmorph.c + + Copyright 2006-2012 + Andrew Janke - a.janke@gmail.com + + Permission to use, copy, modify, and distribute this software and its + documentation for any purpose and without fee is hereby granted, + provided that the above copyright notice appear in all copies. The + author and the University make no representations about the + suitability of this software for any purpose. It is provided "as is" + without express or implied warranty. +*/ + + +#include <config.h> +#include <stdlib.h> +#include <stdio.h> +#include <unistd.h> +#include <sys/param.h> +#include <float.h> +#include <string.h> + +#include <volume_io.h> + +#include <ParseArgv.h> +#include <time_stamp.h> +#include "kernel_io.h" +#include "kernel_ops.h" + +#define INTERNAL_PREC NC_FLOAT /* should be NC_FLOAT or NC_DOUBLE */ +#define DEF_DOUBLE -DBL_MAX + +/* function prototypes */ +char *get_real_from_string(char *string, double *value); +char *get_string_from_string(char *string, char **value); +void calc_volume_range(VIO_Volume * vol, double *min, double *max); +void print_version_info(void); + +/* kernel names for pretty output */ +char *KERN_names[] = { "NULL", "2D04", "2D08", "3D06", "3D26" }; + +/* typedefs */ +typedef enum { + UNDEF = 0, + BINARISE, CLAMP, PAD, ERODE, DILATE, MDILATE, + OPEN, CLOSE, LPASS, HPASS, CONVOLVE, DISTANCE, + GROUP, READ_KERNEL, WRITE, LCORR + } op_types; + +typedef struct { + op_types type; + char op_c; + char *kernel_fn; + kern_types kernel_id; + char *cmpfile; + char *outfile; + double range[2]; + double foreground; + double background; + } Operation; + +/* Argument variables */ +int verbose = FALSE; +int clobber = FALSE; +int is_signed = FALSE; +nc_type dtype = NC_SHORT; +double range[2] = { -DBL_MAX, DBL_MAX }; +double foreground = 1.0; +double background = 0.0; +kern_types kernel_id = K_NULL; +char *kernel_fn = NULL; +char *succ_txt = "B"; + +char successive_help[] = "Successive operations (Maximum: 100) \ +\n\tB[floor:ceil:fg:bg] - binarise in the range, using foreground and background \ +\n\tK[floor:ceil:bg] - clamp betwen the specified range. Set other voxels to 'bg' (default: 0) \ +\n\tP[bg] - pad volume with respect to the current kernel using 'bg' (default: 0)\ +\n\tE - erosion \ +\n\tD - dilation \ +\n\tM - median dilation \ +\n\tO - open \ +\n\tC - close \ +\n\tL - lowpass filter \ +\n\tH - highpass filter \ +\n\tX - convolve \ +\n\tF - distance transform (binary input only - not checked) \ +\n\tG - Label the groups in the volume in ascending order \ +\n\tR[TYPE|file.kern] - (2D04|2D08|3D06|3D26) or read in a kernel file \ +\n\tW[file.mnc] - write out current results \ +\n\tI[cmp.mnc] - local xcorr between current file and cmp.mnc \ +\n\tDefault: "; + +/* Argument table */ +ArgvInfo argTable[] = { + {NULL, ARGV_HELP, (char *)NULL, (char *)NULL, + "General options:"}, + {"-version", ARGV_FUNC, (char *)print_version_info, (char *)NULL, + "print version info and exit"}, + {"-verbose", ARGV_CONSTANT, (char *)TRUE, (char *)&verbose, + "be verbose"}, + {"-clobber", ARGV_CONSTANT, (char *)TRUE, (char *)&clobber, + "clobber existing files"}, + + {NULL, ARGV_HELP, NULL, NULL, + "\nOutfile Options"}, + {"-filetype", ARGV_CONSTANT, (char *)NC_UNSPECIFIED, (char *)&dtype, + "Use data type of the input file."}, + {"-byte", ARGV_CONSTANT, (char *)NC_BYTE, (char *)&dtype, + "Write out byte data."}, + {"-short", ARGV_CONSTANT, (char *)NC_SHORT, (char *)&dtype, + "Write out short integer data. (Default)"}, + {"-int", ARGV_CONSTANT, (char *)NC_INT, (char *)&dtype, + "Write out long integer data."}, + {"-float", ARGV_CONSTANT, (char *)NC_FLOAT, (char *)&dtype, + "Write out single-precision data."}, + {"-double", ARGV_CONSTANT, (char *)NC_DOUBLE, (char *)&dtype, + "Write out double-precision data."}, + {"-signed", ARGV_CONSTANT, (char *)TRUE, (char *)&is_signed, + "Write signed integer data."}, + {"-unsigned", ARGV_CONSTANT, (char *)FALSE, (char *)&is_signed, + "Write unsigned integer data."}, + + {NULL, ARGV_HELP, NULL, NULL, "\nKernel Options"}, + {"-2D04", ARGV_CONSTANT, (char *)K_2D04, (char *)&kernel_id, + "Use a 2D 4-connectivity kernel."}, + {"-2D08", ARGV_CONSTANT, (char *)K_2D08, (char *)&kernel_id, + "Use a 2D 8-connectivity kernel."}, + {"-3D06", ARGV_CONSTANT, (char *)K_3D06, (char *)&kernel_id, + "Use a 3D 6-connectivity kernel. (default)"}, + {"-3D26", ARGV_CONSTANT, (char *)K_3D26, (char *)&kernel_id, + "Use a 3D 26-connectivity kernel."}, + {"-kernel", ARGV_STRING, (char *)1, (char *)&kernel_fn, + "<kernel.kern> read in a custom kernel file"}, + + {NULL, ARGV_HELP, NULL, NULL, "\nMorphology Options"}, + {"-floor", ARGV_FLOAT, (char *)1, (char *)&range[0], + "lowwer value for binarising or clamping"}, + {"-ceil", ARGV_FLOAT, (char *)1, (char *)&range[1], + "upper value for binarising or clamping (incl)"}, + {"-range", ARGV_FLOAT, (char *)2, (char *)range, + "range for binarising or clamping (incl)"}, + {"-foreground", ARGV_FLOAT, (char *)1, (char *)&foreground, + "foreground value"}, + {"-background", ARGV_FLOAT, (char *)1, (char *)&background, + "background value"}, + + {NULL, ARGV_HELP, (char *)NULL, (char *)NULL, "\nSingle morphological operations:"}, + {"-binarise", ARGV_CONSTANT, (char *)"B", (char *)&succ_txt, + "binarise volume using the input range"}, + {"-clamp", ARGV_CONSTANT, (char *)"K", (char *)&succ_txt, + "clamp volume using the input range"}, + {"-pad", ARGV_CONSTANT, (char *)"P", (char *)&succ_txt, + "pad volume with respect to the current kernel (-background to specify value)"}, + {"-erosion", ARGV_CONSTANT, (char *)"E", (char *)&succ_txt, + "do a single erosion"}, + {"-dilation", ARGV_CONSTANT, (char *)"D", (char *)&succ_txt, + "do a single dilation"}, + {"-median_dilation", ARGV_CONSTANT, (char *)"M", (char *)&succ_txt, + "do a single median dilation (note: this is not a median filter!)"}, + {"-open", ARGV_CONSTANT, (char *)"O", (char *)&succ_txt, + "open: dilation(erosion(X))"}, + {"-close", ARGV_CONSTANT, (char *)"C", (char *)&succ_txt, + "close: erosion(dilation(X))"}, + {"-lowpass", ARGV_CONSTANT, (char *)"L", (char *)&succ_txt, + "lowpass filter: close(open(X))"}, + {"-highpass", ARGV_CONSTANT, (char *)"H", (char *)&succ_txt, + "highpass filter: X - lowpass(X)"}, + {"-convolve", ARGV_CONSTANT, (char *)"X", (char *)&succ_txt, + "convolve file with kernel"}, + {"-distance", ARGV_CONSTANT, (char *)"F", (char *)&succ_txt, + "distance transform"}, + {"-group", ARGV_CONSTANT, (char *)"G", (char *)&succ_txt, + "label groups in ascending order"}, + + {NULL, ARGV_HELP, (char *)NULL, (char *)NULL, + "\nSuccessive morphological operations:"}, + {"-successive", ARGV_STRING, (char *)1, (char *)&succ_txt, successive_help}, + + {NULL, ARGV_HELP, NULL, NULL, ""}, + {NULL, ARGV_END, NULL, NULL, NULL} + }; + +int main(int argc, char *argv[]) +{ + int c; + char *arg_string; + char *infile; + char *outfile; + + VIO_Volume *volume; + VIO_Volume *cmpvol; + Kernel *kernel; + int num_ops; + Operation operation[100]; + Operation *op; + char *tmp_str; + char ext_txt[256]; + char tmp_filename[MAXPATHLEN]; + double tmp_double[4]; + double min, max; + char *ptr; + + char *axis_order[3] = { MIzspace, MIyspace, MIxspace }; + + /* Save time stamp and args */ + arg_string = time_stamp(argc, argv); + + /* Get arguments */ + if(ParseArgv(&argc, argv, argTable, 0) || (argc < 2)){ + fprintf(stderr, "\nUsage: %s [options] <in.mnc> <out.mnc>\n", argv[0]); + fprintf(stderr, " %s -help\n\n", argv[0]); + exit(EXIT_FAILURE); + } + infile = argv[1]; + outfile = argv[2]; + + /* check for the infile */ + if(access(infile, F_OK) != 0){ + fprintf(stderr, "%s: Couldn't find %s\n\n", argv[0], infile); + exit(EXIT_FAILURE); + } + + /* check for the outfile */ + if(access(outfile, F_OK) == 0 && !clobber){ + fprintf(stderr, "%s: %s exists! (use -clobber to overwrite)\n\n", argv[0], outfile); + exit(EXIT_FAILURE); + } + + /* check kernel args */ + if(kernel_fn != NULL && kernel_id != K_NULL){ + fprintf(stderr, "%s: specify either a kernel file or a set kernel (not both)\n\n", + argv[0], outfile); + exit(EXIT_FAILURE); + } + + /* set the default kernel */ + if(kernel_fn == NULL && kernel_id == K_NULL){ + kernel_id = K_3D06; + } + + /* add the implicit read kernel operation */ + num_ops = 0; + op = &operation[num_ops++]; + + op->type = READ_KERNEL; + op->kernel_fn = kernel_fn; + op->kernel_id = kernel_id; + + /* setup operations and check them... */ + if(verbose){ + fprintf(stdout, "---Checking Operation(s): %s---\n", succ_txt); + } + + ptr = succ_txt; + while(ptr[0] != '\0'){ + + /* set up counters and extra text */ + strcpy(ext_txt, ""); + op = &operation[num_ops++]; + + /* get the operation type */ + op->op_c = ptr[0]; + ptr++; + + switch (op->op_c){ + case 'B': + op->type = BINARISE; + + /* get 4 possible values */ + ptr = get_real_from_string(ptr, &tmp_double[0]); + ptr = get_real_from_string(ptr, &tmp_double[1]); + ptr = get_real_from_string(ptr, &tmp_double[2]); + ptr = get_real_from_string(ptr, &tmp_double[3]); + op->range[0] = (tmp_double[0] == DEF_DOUBLE) ? range[0] : tmp_double[0]; + op->range[1] = (tmp_double[1] == DEF_DOUBLE) ? range[1] : tmp_double[1]; + op->foreground = (tmp_double[2] == DEF_DOUBLE) ? foreground : tmp_double[2]; + op->background = (tmp_double[3] == DEF_DOUBLE) ? background : tmp_double[3]; + + sprintf(ext_txt, "range: [%g:%g] fg/bg: [%g:%g]", op->range[0], + op->range[1], op->foreground, op->background); + break; + + case 'K': + op->type = CLAMP; + + /* get 3 possible values */ + ptr = get_real_from_string(ptr, &tmp_double[0]); + ptr = get_real_from_string(ptr, &tmp_double[1]); + ptr = get_real_from_string(ptr, &tmp_double[2]); + op->range[0] = (tmp_double[0] == DEF_DOUBLE) ? range[0] : tmp_double[0]; + op->range[1] = (tmp_double[1] == DEF_DOUBLE) ? range[1] : tmp_double[1]; + op->background = (tmp_double[2] == DEF_DOUBLE) ? background : tmp_double[2]; + + sprintf(ext_txt, "range: [%g:%g] fg/bg: [%g:%g]", op->range[0], + op->range[1], op->foreground, op->background); + break; + + case 'P': + op->type = PAD; + + /* get 1 possible value */ + ptr = get_real_from_string(ptr, &tmp_double[0]); + op->background = (tmp_double[0] == DEF_DOUBLE) ? background : tmp_double[0]; + + sprintf(ext_txt, "fill value: %g", op->background); + break; + + case 'E': + op->type = ERODE; + break; + + case 'D': + op->type = DILATE; + break; + + case 'M': + op->type = MDILATE; + break; + + case 'O': + op->type = OPEN; + break; + + case 'C': + op->type = CLOSE; + break; + + case 'L': + op->type = LPASS; + break; + + case 'H': + op->type = HPASS; + break; + + case 'X': + op->type = CONVOLVE; + break; + + case 'F': + op->type = DISTANCE; + break; + + case 'G': + op->type = GROUP; + break; + + case 'R': + op->type = READ_KERNEL; + + /* get the filename */ + ptr = get_string_from_string(ptr, &tmp_str); + if(tmp_str == NULL){ + fprintf(stderr, "%s: R[TYPE|file.kern] requires a file or kernel name\n\n", + argv[0]); + exit(EXIT_FAILURE); + } + + /* check if the input kernel is an inbuilt one */ + op->kernel_id = K_NULL; + for(c = 1; c <= n_inbuilt_kern; c++){ + if(strcmp(tmp_str, KERN_names[c]) == 0){ + op->kernel_id = (kern_types) c; + } + } + + /* if no inbuilt found, assume it's a file */ + if(op->kernel_id == K_NULL){ + + /* set up and check for the real filename */ + (void)realpath(tmp_str, tmp_filename); + if(access(tmp_filename, F_OK) != 0){ + fprintf(stderr, "%s: Couldn't find kernel file: %s\n\n", argv[0], + tmp_filename); + exit(EXIT_FAILURE); + } + + op->kernel_fn = strdup(tmp_filename); + sprintf(ext_txt, "kernel_fn: %s", op->kernel_fn); + } + else { + sprintf(ext_txt, "inbuilt_kernel[%d]: %s", op->kernel_id, + KERN_names[op->kernel_id]); + } + break; + + case 'W': + op->type = WRITE; + + /* get the filename */ + ptr = get_string_from_string(ptr, &op->outfile); + + if(op->outfile == NULL){ + fprintf(stderr, "%s: W[file.mnc] _requires_ a filename\n\n", argv[0]); + exit(EXIT_FAILURE); + } + + /* check for the outfile */ + if(access(op->outfile, F_OK) == 0 && !clobber){ + fprintf(stderr, "%s: %s exists! (use -clobber to overwrite)\n\n", argv[0], + op->outfile); + exit(EXIT_FAILURE); + } + + sprintf(ext_txt, "filename: %s", op->outfile); + break; + + case 'I': + op->type = LCORR; + + /* get the filenames */ + ptr = get_string_from_string(ptr, &op->cmpfile); + + if(op->cmpfile == NULL){ + fprintf(stderr, "%s: I[cmp.mnc] requires a filename\n\n", argv[0]); + exit(EXIT_FAILURE); + } + + /* check for cmpfile */ + if(access(op->cmpfile, F_OK) != 0){ + fprintf(stderr, "%s: Couldn't find compare file: %s\n\n", argv[0], op->cmpfile); + exit(EXIT_FAILURE); + } + + sprintf(ext_txt, "compare filename: %s", op->cmpfile); + break; + + default: + fprintf(stderr, "\nUnknown op: %c\n\n %s -help for operations\n\n", + op->op_c, argv[0]); + exit(EXIT_FAILURE); + } + + if(verbose){ + fprintf(stdout, " Op[%02d] %c = %d\t\t%s\n", num_ops, op->op_c, op->type, + ext_txt); + } + } + + /* add an implicit write statment to the end if needed */ + if(operation[num_ops - 1].type != WRITE){ + operation[num_ops].type = WRITE; + operation[num_ops].outfile = outfile; + num_ops++; + } + + /* malloc space for volume structure and read in infile */ + volume = (VIO_Volume *) malloc(sizeof(VIO_Volume)); + input_volume(infile, MAX_VAR_DIMS, axis_order, + INTERNAL_PREC, TRUE, 0.0, 0.0, TRUE, volume, NULL); + get_type_range(get_volume_data_type(*volume), &min, &max); + set_volume_real_range(*volume, min, max); + + /* init and then do some operations */ + kernel = new_kernel(0); + + if(verbose){ + fprintf(stdout, "\n---Doing %d Operation(s)---\n", num_ops); + } + for(c = 0; c < num_ops; c++){ + op = &operation[c]; + + switch (op->type){ + case BINARISE: + volume = binarise(volume, op->range[0], op->range[1], + op->foreground, op->background); + break; + + case CLAMP: + volume = clamp(volume, op->range[0], op->range[1], op->background); + break; + + case PAD: + volume = pad(kernel, volume, op->background); + break; + + case ERODE: + volume = erosion_kernel(kernel, volume); + break; + + case DILATE: + volume = dilation_kernel(kernel, volume); + break; + + case MDILATE: + volume = median_dilation_kernel(kernel, volume); + break; + + case OPEN: + volume = erosion_kernel(kernel, volume); + volume = dilation_kernel(kernel, volume); + break; + + case CLOSE: + volume = dilation_kernel(kernel, volume); + volume = erosion_kernel(kernel, volume); + break; + + case LPASS: + volume = erosion_kernel(kernel, volume); + volume = dilation_kernel(kernel, volume); + volume = dilation_kernel(kernel, volume); + volume = erosion_kernel(kernel, volume); + break; + + case HPASS: + fprintf(stderr, "%s: GNFARK! Highpass Not implemented yet..\n\n", argv[0]); + break; + + case CONVOLVE: + volume = convolve_kernel(kernel, volume); + break; + + case DISTANCE: + volume = distance_kernel(kernel, volume, background); + break; + + case GROUP: + volume = group_kernel(kernel, volume, background); + break; + + case READ_KERNEL: + /* free the existing kernel then set the pointer to the new one */ + free(kernel); + + /* read in the kernel or set the kernel to an inbuilt one */ + if(op->kernel_id == K_NULL){ + kernel = new_kernel(0); + if(input_kernel(op->kernel_fn, kernel) != OK){ + fprintf(stderr, "%s: Died reading in kernel file: %s\n\n", argv[0], + op->kernel_fn); + exit(EXIT_FAILURE); + } + } + else { + + switch (op->kernel_id){ + case K_2D04: + kernel = get_2D04_kernel(); + break; + + case K_2D08: + kernel = get_2D08_kernel(); + break; + + case K_3D06: + kernel = get_3D06_kernel(); + break; + + case K_3D26: + kernel = get_3D26_kernel(); + break; + + default: + fprintf(stderr, "%s: This shouldn't happen -- much bad\n\n", argv[0]); + exit(EXIT_FAILURE); + } + } + + setup_pad_values(kernel); + if(verbose){ + fprintf(stdout, "Input kernel:\n"); + print_kernel(kernel); + } + break; + + case WRITE: + if(op->outfile == NULL){ + fprintf(stdout, "%s: WRITE passed a NULL pointer! - this is bad\n\n", + argv[0]); + exit(EXIT_FAILURE); + } + + if(verbose){ + fprintf(stdout, "Outputting to %s\n", op->outfile); + } + + /* get the resulting range */ + calc_volume_range(volume, &min, &max); + + /* set the range to something sensible (if possible) */ + if(dtype == NC_BYTE && is_signed == FALSE){ + if(min >= 0 && max < 255){ + fprintf(stdout, "BYTE data, setting 1:1 mapping (0-256)\n"); + min = 0; + max = 255; + } + } + set_volume_real_range(*volume, min, max); + + output_modified_volume(op->outfile, + dtype, is_signed, + 0.0, 0.0, *volume, infile, arg_string, NULL); + break; + + case LCORR: + if(op->cmpfile == NULL){ + fprintf(stdout, "%s: LCORR passed a NULL pointer! - this is bad\n\n", + argv[0]); + exit(EXIT_FAILURE); + } + + if(verbose){ + fprintf(stdout, "Comparing to %s\n", op->cmpfile); + } + + /* malloc space for volume structure and read cmpfile */ + cmpvol = (VIO_Volume *) malloc(sizeof(VIO_Volume)); + input_volume(op->cmpfile, MAX_VAR_DIMS, axis_order, + INTERNAL_PREC, TRUE, 0.0, 0.0, TRUE, cmpvol, NULL); + + /* run the local correlation */ + volume = lcorr_kernel(kernel, volume, cmpvol); + + /* clean up */ + delete_volume(*cmpvol); + free(cmpvol); + + break; + + default: + fprintf(stderr, "\n%s: Unknown operation (This is very bad, call Houston)\n\n", argv[0]); + exit(EXIT_FAILURE); + } + } + + /* jump through operations freeing stuff */ + // free(op.kernel); + + delete_volume(*volume); + return (EXIT_SUCCESS); + } + +/* get a real from a char* stream */ +/* with possible trailing or leading square bracket */ +/* and possible leading ':' */ +/* return the string advanced to the next token or */ +/* as it was input if nothing found */ +char *get_real_from_string(char *string, double *value) +{ + char *ptr; + + /* skip a [ or : else we probably don't belong here */ + if(string[0] == '[' || string[0] == ':'){ + string++; + } + + /* get a double */ + *value = strtod(&string[0], &ptr); + + /* if nothing found return a default value */ + if(&string[0] == ptr){ + *value = DEF_DOUBLE; + } + + /* skip over a possible ']' */ + if(ptr[0] == ']'){ + ptr++; + } + + return ptr; + } + +/* get a copy of the string between [ and ] from a char* */ +/* return the string advanced to the next token or */ +/* as it was input if nothing found */ +char *get_string_from_string(char *string, char **value) +{ + int offset; + char *malloc_string; + + /* initialise the return values first */ + *value = NULL; + + /* get the string if there is one */ + if(string[0] == '['){ + string++; + + /* get the length of the string in question */ + offset = strcspn(string, "]"); + + /* alloc some space for the string */ + malloc_string = (char *)malloc((offset + 1) * sizeof(char)); + + /* get a copy of it and append the null charater */ + strncpy(malloc_string, string, offset); + malloc_string[offset] = '\0'; + + /* store the result */ + *value = malloc_string; + + /* increment the pointer and skip over a possible '[' */ + string += offset; + if(string[0] == ']'){ + string++; + } + } + + return string; + } + +void calc_volume_range(VIO_Volume * vol, double *min, double *max) +{ + + int x, y, z; + int sizes[MAX_VAR_DIMS]; + double value; + VIO_progress_struct progress; + + *min = DBL_MAX; + *max = -DBL_MIN; + + get_volume_sizes(*vol, sizes); + + initialize_progress_report(&progress, FALSE, sizes[2], "Finding Range"); + for(z = sizes[0]; z--;){ + for(y = sizes[1]; y--;){ + for(x = sizes[2]; x--;){ + + value = get_volume_voxel_value(*vol, z, y, x, 0, 0); + if(value < *min){ + *min = value; + } + else if(value > *max){ + *max = value; + } + } + } + update_progress_report(&progress, z + 1); + } + terminate_progress_report(&progress); + + if (*min == *max) { + *max = *min + 1.0; + } + + if(verbose){ + fprintf(stdout, "Found range of [%g:%g]\n", *min, *max); + } + } + +void print_version_info(void) +{ + fprintf(stdout, "%s version %s\n", PACKAGE, VERSION); + fprintf(stdout, "Comments to %s\n", PACKAGE_BUGREPORT); + fprintf(stdout, "\n"); + exit(EXIT_SUCCESS); + }