From da9831346a167dfbb07495f93772e64e13f7c518 Mon Sep 17 00:00:00 2001 From: bert Date: Tue, 1 Mar 2005 17:21:53 +0000 Subject: [PATCH] Added normalize_pet to the conglomerate --- ChangeLog | 3 + Makefile.am | 2 + normalize_pet.c | 341 ++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 346 insertions(+) create mode 100644 normalize_pet.c diff --git a/ChangeLog b/ChangeLog index f61ae64..2be5b98 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,6 @@ +2005-03-01 Bert Vincent + * Added normalize_pet + 2004-04-22 Bert Vincent * Added find_peaks * Update version to 1.4 diff --git a/Makefile.am b/Makefile.am index b31378d..c35136e 100644 --- a/Makefile.am +++ b/Makefile.am @@ -117,6 +117,7 @@ bin_PROGRAMS = \ merge_polygons \ minctotag \ minc_to_rgb \ + normalize_pet \ perturb_surface \ place_images \ plane_polygon_intersect \ @@ -411,6 +412,7 @@ merge_polygons_SOURCES = merge_polygons.c minc_to_rgb_SOURCES = minc_to_rgb.c minctotag_SOURCES = minctotag.c #normalize_intensity_SOURCES = normalize_intensity.c +normalize_pet_SOURCES = normalize_pet.c #obj_to_rgb_SOURCES = obj_to_rgb.c #offset_volume2_SOURCES = offset_volume2.c #offset_volume_SOURCES = offset_volume.c diff --git a/normalize_pet.c b/normalize_pet.c new file mode 100644 index 0000000..652f076 --- /dev/null +++ b/normalize_pet.c @@ -0,0 +1,341 @@ +/* ----------------------------- MNI Header ----------------------------------- +@NAME : normalize_pet.c +@INPUT : argc, argv - command line arguments +@OUTPUT : (nothing) +@RETURNS : status +@DESCRIPTION: Program to normalize a volume to a standard intensity. +@CREATED : October 27, 1993 (Peter Neelin) +@MODIFIED : $Log: normalize_pet.c,v $ +@MODIFIED : Revision 1.1 2005-03-01 17:21:54 bert +@MODIFIED : Added normalize_pet to the conglomerate +@MODIFIED : + * Revision 1.9 1998/06/25 16:23:57 neelin + * Added check for error on writing output volume. + * + * Revision 1.8 1996/05/06 11:57:59 neelin + * Changed name to normalize_pet/ + * + * Revision 1.7 1995/04/21 13:23:12 neelin + * Added options to use mask file and normalize to a value other than 100. + * + * Revision 1.6 1995/04/11 12:36:50 neelin + * *** empty log message *** + * + * Revision 1.5 93/10/29 15:31:40 neelin + * Delete norm volume before loading input volume (so that type and range + * don't get changed). + * + * Revision 1.4 93/10/29 10:44:24 neelin + * Pass options to input_volume. + * + * Revision 1.3 93/10/28 16:53:00 neelin + * Removed extra GET_VOXEL_3D call. + * + * Revision 1.2 93/10/28 16:12:14 neelin + * Modified usage line. + * + * Revision 1.1 93/10/28 15:21:14 neelin + * Initial revision + * +---------------------------------------------------------------------------- */ + +#include +#include +#include +#include +#include +#include + +/* Constants */ +#define RANGE_EPSILON (10.0 * FLT_EPSILON) + +/* Function prototypes */ +int match_volume_dimensions(Volume vol1, Volume vol2); +int match_reals(Real value1, Real value2); + +/* Variables for arguments */ +char *maskfile = NULL; +double normalized_mean = 100.0; + +/* Argument table */ +ArgvInfo argTable[] = { + {"-mask", ARGV_STRING, (char *) 1, (char *) &maskfile, + "Use the specified file as a mask for normalization."}, + {"-normalized_mean", ARGV_FLOAT, (char *) 1, (char *) &normalized_mean, + "Specify a target value for the normalized mean."}, + {NULL, ARGV_END, NULL, NULL, NULL} +}; + +int main(int argc, char *argv[]) +{ + char *normfile, *infile, *outfile; + Volume volume; + Volume mask; + int sizes[3]; + int ix, iy, iz; + Real sum1, sumv, value, mean, threshold, norm_factor, mask_value; + Real maximum, minimum, mask_min, mask_max; + char *history, *pname; + minc_input_options options; + int need_threshold; + + /* Check arguments */ + history = time_stamp(argc, argv); + pname=argv[0]; + if (ParseArgv(&argc, argv, argTable, 0)) { + (void) fprintf(stderr, + "\nUsage: %s [] [] \n", + pname); + return EXIT_FAILURE; + } + if (argc == 3) { + normfile = argv[1]; + infile = NULL; + outfile = argv[2]; + } + else if (argc == 4) { + normfile = argv[1]; + infile = argv[2]; + outfile = argv[3]; + } + else { + (void) fprintf(stderr, + "\nUsage: %s [] [] \n", + pname); + return EXIT_FAILURE; + } + need_threshold = (maskfile == NULL); + + /* Set input options */ + set_default_minc_input_options(&options); + set_minc_input_promote_invalid_to_min_flag(&options, FALSE); + set_minc_input_vector_to_scalar_flag(&options, FALSE); + + /* Read in the normalization file */ + if (input_volume(normfile, 3, NULL, NC_UNSPECIFIED, FALSE, 0.0, 0.0, + TRUE, &volume, &options) != OK) { + (void) fprintf(stderr, "Error loading volume %s\n", normfile); + return EXIT_FAILURE; + } + get_volume_real_range(volume, &minimum, &maximum); + minimum -= ABS(minimum) * RANGE_EPSILON; + maximum += ABS(maximum) * RANGE_EPSILON; + get_volume_sizes(volume, sizes); + + /* Read in the mask file and check that the dimensions match */ + if (maskfile != NULL) { + if (input_volume(maskfile, 3, NULL, NC_UNSPECIFIED, FALSE, 0.0, 0.0, + TRUE, &mask, &options) != OK) { + (void) fprintf(stderr, "Error loading volume %s\n", maskfile); + return EXIT_FAILURE; + } + if (!match_volume_dimensions(volume, mask)) { + (void) fprintf(stderr, + "Mask volume dimensions do not match normalization volume.\n"); + return EXIT_FAILURE; + } + get_volume_real_range(mask, &mask_min, &mask_max); + mask_min -= ABS(mask_min) * RANGE_EPSILON; + mask_max += ABS(mask_max) * RANGE_EPSILON; + } + + /* Calculate the threshold if needed */ + if (need_threshold) { + sum1 = 0.0; + sumv = 0.0; + for (iz=0; iz= minimum) && (value <= maximum)) { + sum1 += 1.0; + sumv += value; + } + } + } + } + if (sum1 > 0.0) + mean = sumv / sum1; + else + mean = 0.0; + threshold = 1.5 * mean; + } + else { + threshold = 0.0; + } + + /* Calculate the normalization factor */ + sum1 = 0.0; + sumv = 0.0; + mask_value = 1.0; + for (iz=0; iz mask_max)) + mask_value = 0.0; + } + if ((!need_threshold || (value > threshold)) && + (value >= minimum) && (value <= maximum)) { + sum1 += mask_value; + sumv += value * mask_value; + } + } + } + } + if (sum1 > 0.0) + mean = sumv / sum1; + else + mean = 0.0; + if (mean <= 0.0) { + (void) fprintf(stderr, "Thresholded mean <= 0.0\n"); + return EXIT_FAILURE; + } + norm_factor = normalized_mean / mean; + + /* Get rid of the mask volume */ + if (maskfile != NULL) { + delete_volume(mask); + } + + /* Write out the numbers */ + (void) printf("File = %s\n", normfile); + (void) printf("Threshold = %.20g\n", (double) threshold); + (void) printf("Mean = %.20g\n", (double) mean); + (void) printf("Norm factor = %.20g\n", (double) norm_factor); + + /* Load file to normalize if needed */ + if (infile != NULL) { + delete_volume(volume); + if (input_volume(infile, 3, NULL, NC_UNSPECIFIED, FALSE, 0.0, 0.0, + TRUE, &volume, &options) != OK) { + (void) fprintf(stderr, "Error loading volume %s\n", infile); + return EXIT_FAILURE; + } + get_volume_sizes(volume, sizes); + } + else { + infile = normfile; + } + + /* Normalize volume */ + get_volume_real_range(volume, &minimum, &maximum); + minimum *= norm_factor; + maximum *= norm_factor; + set_volume_real_range(volume, minimum, maximum); + if ((volume->nc_data_type == NC_FLOAT) || + (volume->nc_data_type == NC_DOUBLE)) { + for (iz=0; iz= 0) + voxel[idim] = sizes1[idim] - 1; + convert_voxel_to_world(vol1, voxel, &x1, &y1, &z1); + convert_voxel_to_world(vol2, voxel, &x2, &y2, &z2); + if (!(match_reals(x1, x2) || match_reals(y1, y2) || + match_reals(z1, z2))) { + return FALSE; + } + if (idim >= 0) + voxel[idim] = 0.0; + } + + return TRUE; +} + +/* ----------------------------- MNI Header ----------------------------------- +@NAME : match_reals +@INPUT : value1 - real value 1 + value2 - real value 2 +@OUTPUT : (nothing) +@RETURNS : TRUE if values match, FALSE if not +@DESCRIPTION: Routine to compare two Real values, allowing a tolerance. +@METHOD : +@GLOBALS : +@CALLS : +@CREATED : April 21, 1995 (Peter Neelin) +@MODIFIED : +---------------------------------------------------------------------------- */ +int match_reals(Real value1, Real value2) +{ +#define TOLERANCE ( 100.0 * FLT_EPSILON ) + + Real diff, denom; + + diff = value1 - value2; + denom = (value1 + value2) / 2.0; + if (denom != 0.0) diff /= denom; + diff = ABS(diff); + if (diff < TOLERANCE) + return TRUE; + else + return FALSE; +}