Mercurial > hg > minc-tools
changeset 994:875912d658f7
Various changes including changed coordinates in x and y and computing
FWHM from cutoff frequency.
author | neelin <neelin> |
---|---|
date | Tue, 26 Mar 1996 15:58:18 +0000 |
parents | c8d2b3e04b6c |
children | 26eb4481c907 |
files | conversion/ecattominc/ecattominc.c |
diffstat | 1 files changed, 176 insertions(+), 28 deletions(-) [+] |
line wrap: on
line diff
--- a/conversion/ecattominc/ecattominc.c +++ b/conversion/ecattominc/ecattominc.c @@ -9,9 +9,13 @@ @CALLS : @CREATED : January 3, 1996 (Peter Neelin) @MODIFIED : $Log: ecattominc.c,v $ -@MODIFIED : Revision 1.1 1996-01-18 14:52:14 neelin -@MODIFIED : Initial revision +@MODIFIED : Revision 1.2 1996-03-26 15:58:18 neelin +@MODIFIED : Various changes including changed coordinates in x and y and computing +@MODIFIED : FWHM from cutoff frequency. @MODIFIED : + * Revision 1.1 1996/01/18 14:52:14 neelin + * Initial revision + * @COPYRIGHT : Copyright 1996 Peter Neelin, McConnell Brain Imaging Centre, Montreal Neurological Institute, McGill University. @@ -25,12 +29,13 @@ ---------------------------------------------------------------------------- */ #ifndef lint -static char rcsid[]="$Header: /private-cvsroot/minc/conversion/ecattominc/ecattominc.c,v 1.1 1996-01-18 14:52:14 neelin Exp $"; +static char rcsid[]="$Header: /private-cvsroot/minc/conversion/ecattominc/ecattominc.c,v 1.2 1996-03-26 15:58:18 neelin Exp $"; #endif #include <stdlib.h> #include <stdio.h> #include <string.h> +#include <ctype.h> #include <math.h> #include <limits.h> #include <float.h> @@ -56,6 +61,7 @@ double half_life; double zstart; double zstep; + double decay_correction; char isotope[16]; int image_xsize; int image_ysize; @@ -71,11 +77,13 @@ int max_nslices; int max_xsize; int max_ysize; - double zwidth; + double xstart; + double ystart; double xstep; double ystep; double xwidth; double ywidth; + double zwidth; int decay_corrected; char img_units[16]; char patient_name[40]; @@ -96,6 +104,7 @@ long injection_minute; double injection_seconds; double injection_dose; + int septa_retracted; ecat_header_data_type *main_field_list; int num_main_fields; ecat_header_data_type *subhdr_field_list; @@ -141,11 +150,17 @@ #define FALSE 0 #define MAX_DIMS 4 #define ECAT_ACTIVITY "ACTIVITY" +#define ECAT_CALIB_UNITS_UNKNOWN 0 +#define ECAT_CALIB_UNITS_BECQUEREL 1 #define ECAT_CALIB_UNITS_CPS 3 #define MM_PER_CM 10.0 -#define BECQUEREL_PER_MCURIE 3.7e7 +#define BECQUEREL_PER_NCURIE 37 +#define BECQUEREL_PER_MCURIE (BECQUEREL_PER_NCURIE * 1e6) +#define NCURIE_PER_CC_STRING "nCi/cc" #define SECONDS_PER_HOUR 3600 #define DEFAULT_RANGE INT_MIN +#define MINIMUM_HALFLIFE 0.1 +#define FWHM_SCALE_FOR_HANN 1.082 /* Main program */ @@ -157,7 +172,7 @@ static int verbose = TRUE; static int decay_correct = TRUE; static int slice_range[2] = {DEFAULT_RANGE, DEFAULT_RANGE}; - static int copy_all_header = FALSE; + static int copy_all_header = TRUE; static char *blood_file = NULL; static int frame_range[2] = {DEFAULT_RANGE, DEFAULT_RANGE}; @@ -187,9 +202,9 @@ "Single frame to copy (counting from 0)."}, {"-small_header", ARGV_CONSTANT, (char *) FALSE, (char *) ©_all_header, - "Copy only basic header information (default)."}, + "Copy only basic header information."}, {"-all_header", ARGV_CONSTANT, (char *) TRUE, (char *) ©_all_header, - "Copy all header information."}, + "Copy all header information (default)."}, {"-bloodfile", ARGV_STRING, (char *) 1, (char *) &blood_file, "Insert blood data from this file."}, {NULL, ARGV_END, NULL, NULL, NULL} @@ -383,6 +398,9 @@ 0.0, frame_info[iframe].half_life); } + else if (!decay_correct && general_info->decay_corrected) { + scale = 1.0 / frame_info[iframe].decay_correction; + } else { scale = 1.0; } @@ -541,6 +559,10 @@ Ecat_field_name field; char *description; time_t the_time; + char *ptr; + int isotope_name_okay; + int septa_state; + double cutoff, binsize; /* Initialize number of slices */ num_slices = &(general_info->num_slices); @@ -596,7 +618,9 @@ /* Get scan type */ if (ecat_get_main_value(ecat_fp, ECAT_Calibration_Units, 0, &lvalue, NULL, NULL)) return curframe; - if (lvalue == ECAT_CALIB_UNITS_CPS) { + if ((lvalue == ECAT_CALIB_UNITS_BECQUEREL) || + (lvalue == ECAT_CALIB_UNITS_UNKNOWN) || + (lvalue == ECAT_CALIB_UNITS_CPS)) { (void) strcpy(fip->image_type, ECAT_ACTIVITY); } else { @@ -611,6 +635,23 @@ NULL, &fip->half_life, NULL)) return curframe; + /* Check that they are reasonable */ + isotope_name_okay = TRUE; + for (ptr=fip->isotope; + (*ptr != '\0') && (ptr < &fip->isotope[sizeof(fip->isotope)]); + ptr++) { + if (!isprint((int) *ptr)) { + isotope_name_okay = FALSE; + } + } + if (ptr == fip->isotope) + isotope_name_okay = FALSE; + if (!isotope_name_okay || (fip->half_life < MINIMUM_HALFLIFE)) { + (void) fprintf(stderr, "Ignoring bad isotope name or half-life.\n"); + fip->isotope[0] = '\0'; + fip->half_life = 0.0; + } + /* Get z start and step (correct start for non-zero first slice */ if (ecat_get_main_value(ecat_fp, ECAT_Plane_Separation, 0, NULL, &fip->zstep, NULL)) return curframe; @@ -627,8 +668,13 @@ fip->zstart += fip->low_slice * fip->zstep; /* Check to see if file has been decay corrected */ - if (!ecat_get_subhdr_value(ecat_fp, curframe, 0, ECAT_Decay_Corr_Fctr, 0, - NULL, &fvalue, NULL) && fvalue != 1.0) { + fvalue = 1.0; + (void) ecat_get_subhdr_value(ecat_fp, curframe, 0, + ECAT_Decay_Corr_Fctr, 0, + NULL, &fvalue, NULL); + if (fvalue <= 0.0) fvalue = 1.0; + fip->decay_correction = fvalue; + if (fip->decay_correction != 1.0) { general_info->decay_corrected = TRUE; } @@ -645,6 +691,20 @@ general_info->xstep *= MM_PER_CM; general_info->ystep *= MM_PER_CM; + /* Get location of first voxel */ + if (ecat_get_subhdr_value(ecat_fp, curframe, 0, ECAT_X_Offset, 0, + NULL, &general_info->xstart, NULL)) + return curframe; + if (ecat_get_subhdr_value(ecat_fp, curframe, 0, ECAT_Y_Offset, 0, + NULL, &general_info->ystart, NULL)) + return curframe; + general_info->xstart *= -MM_PER_CM; + general_info->ystart *= MM_PER_CM; + general_info->xstart -= + general_info->xstep * ((double) fip->image_xsize - 1.0) / 2.0; + general_info->ystart -= + general_info->ystep * ((double) fip->image_ysize - 1.0) / 2.0; + /* Get resolution in each direction (or zero if not found) */ general_info->xwidth = general_info->ywidth = general_info->zwidth = -1.0; @@ -658,11 +718,44 @@ general_info->ywidth *= MM_PER_CM; general_info->zwidth *= MM_PER_CM; + /* If resolution is not found, then use cutoff frequency and + assume FWHM for Hann filter */ + cutoff = -1.0; + binsize = -1.0; + (void) ecat_get_subhdr_value(ecat_fp, curframe, 0, + ECAT_Rfilter_Cutoff, + 0, NULL, &cutoff, NULL); + if (cutoff <= 0) { + (void) ecat_get_subhdr_value(ecat_fp, curframe, 0, + ECAT_Filter_Cutoff_Frequency, + 0, NULL, &cutoff, NULL); + } + (void) ecat_get_main_value(ecat_fp, ECAT_Bin_Size, 0, + NULL, &binsize, NULL); + binsize *= MM_PER_CM; + if ((general_info->xwidth <= 0.0) && + (cutoff > 0.0) && (binsize > 0.0)) { + general_info->xwidth = FWHM_SCALE_FOR_HANN * binsize / cutoff; + } + if ((general_info->ywidth <= 0.0) && + (cutoff > 0.0) && (binsize > 0.0)) { + general_info->ywidth = FWHM_SCALE_FOR_HANN * binsize / cutoff; + } + if ((general_info->zwidth <= 0.0) && + !ecat_get_subhdr_value(ecat_fp, curframe, 0, + ECAT_Zfilter_Cutoff, + 0, NULL, &cutoff, NULL) && + !ecat_get_main_value(ecat_fp, ECAT_Plane_Separation, 0, + NULL, &binsize, NULL)) { + binsize *= MM_PER_CM; + general_info->zwidth = FWHM_SCALE_FOR_HANN * binsize / cutoff; + } + /* Get image range and units */ if (ecat_get_main_value(ecat_fp, ECAT_Calibration_Units, 0, &lvalue, NULL, NULL)) return curframe; - if (lvalue == ECAT_CALIB_UNITS_CPS) { - (void) strcpy(general_info->img_units, "nCi/cc"); + if (lvalue == ECAT_CALIB_UNITS_BECQUEREL) { + (void) strcpy(general_info->img_units, NCURIE_PER_CC_STRING); } else { (void) strcpy(general_info->img_units, ""); @@ -754,9 +847,14 @@ (int) general_info->start_hour, (int) general_info->start_minute, (int) general_info->start_seconds); - if (ecat_get_main_value(ecat_fp, ECAT_Radiopharmaceutical, 0, - NULL, NULL, general_info->tracer)) - return curframe; + if ((int) strlen(fip->isotope) > 0) { + if (ecat_get_main_value(ecat_fp, ECAT_Radiopharmaceutical, 0, + NULL, NULL, general_info->tracer)) + return curframe; + } + else { + general_info->tracer[0] = '\0'; + } if (!ecat_get_main_value(ecat_fp, ECAT_Dose_Start_Time, 0, &lvalue, NULL, NULL) && (lvalue != 0)) { the_time = (time_t) lvalue; @@ -778,6 +876,13 @@ general_info->injection_dose = -1.0; } + /* Get septa state */ + septa_state = 0; + if (!ecat_get_main_value(ecat_fp, ECAT_Septa_State, 0, + &septa_state, NULL, NULL)) { + general_info->septa_retracted = (septa_state == 1); + } + /* Get list of header values */ for (iheader=0; iheader < 2; iheader++) { @@ -1031,11 +1136,13 @@ int dim[MAX_DIMS]; int img, imgmax, imgmin, dimvarid, widvarid, icv, varid, ecat_var; int bloodid; - int idim, ifield, num_fields, iheader; + int idim, ifield, num_fields, iheader, iframe; double vrange[2]; char varname[MAX_NC_NAME]; ecat_header_data_type *field_list; double dimwidths[MAX_DIMS]; + double *frame_times; + double *frame_lengths; /* Set up dimension arrays for looping */ dim_names = dim_names_array + MAX_DIMS - ndims; @@ -1053,7 +1160,7 @@ for (idim=0; idim<ndims; idim++) { dim[idim]=ncdimdef(mincid, dim_names[idim], count[idim]); dimvarid=micreate_std_variable(mincid, dim_names[idim], NC_DOUBLE, - ((idim==0) && (num_frames>1)) ? 1 : 0, + (((idim==0) && (num_frames>1)) ? 1 : 0), &dim[idim]); if (dimwidths[idim] > 0) { widvarid=micreate_std_variable(mincid, dimwidth_names[idim], @@ -1084,6 +1191,8 @@ } else if (strcmp(dim_names[idim], MIyspace)==0) { (void) miattputstr(mincid, dimvarid, MIunits, "mm"); + (void) miattputdbl(mincid, dimvarid, MIstart, + (double) general_info->ystart); (void) miattputdbl(mincid, dimvarid, MIstep, (double) general_info->ystep); (void) miattputstr(mincid, dimvarid, MIspacetype, MI_NATIVE); @@ -1095,6 +1204,8 @@ } else if (strcmp(dim_names[idim], MIxspace)==0) { (void) miattputstr(mincid, dimvarid, MIunits, "mm"); + (void) miattputdbl(mincid, dimvarid, MIstart, + (double) general_info->xstart); (void) miattputdbl(mincid, dimvarid, MIstep, (double) general_info->xstep); (void) miattputstr(mincid, dimvarid, MIspacetype, MI_NATIVE); @@ -1173,14 +1284,18 @@ /* Save acquisition info */ varid = micreate_group_variable(mincid, MIacquisition); - (void) miattputstr(mincid, varid, MIradionuclide, - frame_info[0].isotope); + if ((int) strlen(frame_info[0].isotope) > 0) { + (void) miattputstr(mincid, varid, MIradionuclide, + frame_info[0].isotope); + } if (frame_info[0].half_life > 0.0) { (void) miattputdbl(mincid, varid, MIradionuclide_halflife, (double) frame_info[0].half_life); } - (void) miattputstr(mincid, varid, MItracer, - general_info->tracer); + if ((int) strlen(general_info->tracer) > 0) { + (void) miattputstr(mincid, varid, MItracer, + general_info->tracer); + } if ((int) strlen(general_info->injection_time) > 0) { (void) miattputstr(mincid, varid, MIinjection_time, general_info->injection_time); @@ -1197,6 +1312,28 @@ (void) miattputstr(mincid, varid, MIdose_units, "mCurie"); } + /* Save the septa state in a special ECAT variable, along with frame + starts and lengths if we are not creating a time dimension. */ + varid = ncvardef(mincid, "ecat_acquisition", NC_LONG, 0, NULL); + (void) miattputstr(mincid, varid, MIvartype, MI_GROUP); + (void) miattputstr(mincid, varid, MIvarid, + "ECAT-specific acquisition information"); + (void) miadd_child(mincid, ncvarid(mincid, MIrootvariable), varid); + (void) miattputstr(mincid, varid, "septa_retracted", + (general_info->septa_retracted ? MI_TRUE : MI_FALSE)); + if (ndims < MAX_DIMS) { + frame_times = MALLOC(sizeof(*frame_times) * num_frames); + frame_lengths = MALLOC(sizeof(*frame_lengths) * num_frames); + for (iframe=0; iframe < num_frames; iframe++) { + frame_times[iframe] = frame_info[iframe].scan_time; + frame_lengths[iframe] = frame_info[iframe].time_width; + } + (void) ncattput(mincid, varid, "frame_times", NC_DOUBLE, num_frames, + frame_times); + (void) ncattput(mincid, varid, "frame_lengths", NC_DOUBLE, num_frames, + frame_lengths); + } + /* If we want all of the values from the header, get them */ if (copy_all_header) { @@ -1277,6 +1414,7 @@ { long npix, ix, iy, y_offset, off1, off2; int pmax; + int lvalue; short temp; double scale, global_scale; @@ -1302,10 +1440,20 @@ ecat_get_subhdr_value(ecat_fp, frame_num, slice_num, ECAT_Scale_Factor, 0, NULL, &scale, NULL)) return TRUE; - global_scale = 1.0; - (void) ecat_get_main_value(ecat_fp, ECAT_Calibration_Factor, 0, - NULL, &global_scale, NULL); - if (global_scale <= 0.0) global_scale = 1.0; + if (!ecat_get_main_value(ecat_fp, ECAT_Calibration_Factor, 0, + NULL, &global_scale, NULL) && + (global_scale > 0.0)) { + if (!ecat_get_main_value(ecat_fp, ECAT_Calibration_Units, 0, + &lvalue, NULL, NULL) && + (lvalue == ECAT_CALIB_UNITS_BECQUEREL)) { + global_scale /= BECQUEREL_PER_NCURIE; + } + } + else if (ecat_get_subhdr_value(ecat_fp, frame_num, slice_num, + ECAT_Calibration_Factor, 0, + NULL, &global_scale, NULL)) { + global_scale = 1.0; + } *pixel_max = pmax; *image_max = (double) pmax * scale * global_scale; @@ -1342,7 +1490,7 @@ ---------------------------------------------------------------------------- */ int write_minc_slice(double scale, int write_byte_data, int mincid, int icvid, - int ndims,long start[], long count[], + int ndims, long start[], long count[], short *image, int image_xsize, int image_ysize, long pixel_max, double image_max, double scan_time, double time_width, double zpos) @@ -1382,7 +1530,7 @@ (void) ncvarput1(mincid, ncvarid(mincid, MIimagemin), start, &minimum); /* Write out time and z position */ - if (ndims==MAX_DIMS) { + if (ndims == MAX_DIMS) { (void) mivarput1(mincid, ncvarid(mincid, MItime), start, NC_DOUBLE, NULL, &scan_time); (void) mivarput1(mincid, ncvarid(mincid, MItime_width), start,