1 /* ----------------------------- MNI Header -----------------------------------
2 @NAME : scxtominc
3 @INPUT : argc, argv - command line arguments
4 @OUTPUT : (none)
5 @RETURNS : error status
6 @DESCRIPTION: Converts scanditronix format files to a minc format file.
7 @METHOD :
8 @GLOBALS :
9 @CALLS :
10 @CREATED : January 11, 1993 (Peter Neelin)
11 @MODIFIED :
12 * $Log: scxtominc.c,v $
13 * Revision 6.4 2008-01-17 02:33:02 rotor
14 * * removed all rcsids
15 * * removed a bunch of ^L's that somehow crept in
16 * * removed old (and outdated) BUGS file
17 *
18 * Revision 6.3 2008/01/12 19:08:14 stever
19 * Add __attribute__ ((unused)) to all rcsid variables.
20 *
21 * Revision 6.2 1999/11/09 13:34:48 neelin
22 * Year 2000 fix for date stored in minc file.
23 *
24 * Revision 6.1 1999/10/29 17:52:07 neelin
25 * Fixed Log keyword
26 *
27 * Revision 6.0 1997/09/12 13:23:31 neelin
28 * Release of minc version 0.6
29 *
30 * Revision 5.0 1997/08/21 13:24:32 neelin
31 * Release of minc version 0.5
32 *
33 * Revision 4.0 1997/05/07 20:00:13 neelin
34 * Release of minc version 0.4
35 *
36 * Revision 3.1 1996/01/04 13:41:48 neelin
37 * Added missing exit when user specifies a bad slice range.
38 *
39 * Revision 3.0 1995/05/15 19:31:05 neelin
40 * Release of minc version 0.3
41 *
42 * Revision 2.5 1995/02/09 14:11:43 neelin
43 * Mods to make irix 5 lint happy.
44 *
45 * Revision 2.4 1995/02/08 19:31:47 neelin
46 * Moved ARGSUSED statements for irix 5 lint.
47 *
48 * Revision 2.3 1995/01/23 09:21:13 neelin
49 * Changed ncclose to miclose
50 *
51 * Revision 2.2 95/01/23 08:57:37 neelin
52 * Changed nccreate to micreate.
53 *
54 * Revision 2.1 95/01/09 15:17:37 neelin
55 * Added code to check for generic reconstructions.
56 *
57 * Revision 2.0 94/09/28 10:33:50 neelin
58 * Release of minc version 0.2
59 *
60 * Revision 1.12 94/09/28 10:33:30 neelin
61 * Pre-release
62 *
63 * Revision 1.11 94/09/28 08:23:58 neelin
64 * Find max and min pixel values for both bytes and shorts.
65 * (Shorts are rescaled to full range).
66 *
67 * Revision 1.10 94/05/31 07:56:42 neelin
68 * Added insertblood.c to optionally insert blood data into minc file.
69 *
70 * Revision 1.9 93/11/17 12:21:55 neelin
71 * Changed default to -noclobber.
72 *
73 * Revision 1.8 93/08/11 15:27:34 neelin
74 * Added RCS logging in source.
75 *
76 @COPYRIGHT :
77 Copyright 1993 Peter Neelin, McConnell Brain Imaging Centre,
78 Montreal Neurological Institute, McGill University.
79 Permission to use, copy, modify, and distribute this
80 software and its documentation for any purpose and without
81 fee is hereby granted, provided that the above copyright
82 notice appear in all copies. The author and McGill University
83 make no representations about the suitability of this
84 software for any purpose. It is provided "as is" without
85 express or implied warranty.
86 ---------------------------------------------------------------------------- */
87
88 #include <stdlib.h>
89 #include <stdio.h>
90 #include <string.h>
91 #include <math.h>
92 #include <float.h>
93 #include <ParseArgv.h>
94 #include <time_stamp.h>
95 #include <scx_file.h>
96 #include <minc.h>
97 #include <minc_def.h>
98 #include "isotope_list.h"
99
100 /* Macro definitions */
101 #undef MALLOC
102 #undef REALLOC
103 #undef FREE
104 #define MALLOC( n_items, type ) \
105 ( (void *) malloc( (size_t) (n_items) * sizeof(type) ) )
106 #define REALLOC( ptr, n_items, type ) \
107 ( (void *) realloc( (void *) ptr, (size_t) (n_items) * sizeof(type) ) )
108 #define FREE( ptr ) \
109 free( (void *) ptr )
110
111 /* Type declarations */
112 typedef struct {
113 char name[8];
114 int mult;
115 nc_type type;
116 void *att_vector;
117 } scx_mnem_list_type;
118
119 typedef struct {
120 int nslices;
121 int low_slice;
122 int high_slice;
123 double scan_time;
124 double time_width;
125 double half_life;
126 double zstart;
127 double zstep;
128 char isotope[16];
129 int image_size;
130 int ordered_file;
131 int *ordered_slices;
132 char image_type[16];
133 } scx_file_info_type;
134
135 typedef struct {
136 int num_scx_slices;
137 int max_nslices;
138 int max_size;
139 float zwidth;
140 float xystep;
141 float xywidth;
142 long vmax;
143 char img_units[16];
144 char patient_name[32];
145 char patient_sex[8];
146 long patient_age;
147 char study_id[40];
148 char start_time[40];
149 long start_year;
150 long start_month;
151 long start_day;
152 long start_hour;
153 long start_minute;
154 float start_seconds;
155 char tracer[10];
156 char injection_time[40];
157 long injection_hour;
158 long injection_minute;
159 float injection_seconds;
160 float injection_dose;
161 scx_mnem_list_type *mnem_list;
162 int num_mnems;
163 int used_MNI_generic_reconstruction;
164 } scx_general_info_type;
165
166 typedef struct {
167 double sort_key;
168 void *sort_value;
169 } scx_sort_type;
170
171 /* Function declarations */
172 void usage_error(char *progname);
173 int get_scx_file_info(int num_scx_files, char **scx_files,
174 int slice_range[2],
175 scx_file_info_type *scx_file_info,
176 scx_general_info_type *scx_general_info);
177 void sort_scx_slices(int sort_over_time, int num_scx_files,
178 scx_file_info_type *scx_file_info,
179 scx_general_info_type *scx_general_info);
180 int sortcmp(const void *val1, const void *val2);
181 int setup_minc_file(int mincid, int write_byte_data, int copy_all_header,
182 int ndims, long count[], int num_scx_files,
183 scx_file_info_type *scx_file_info,
184 scx_general_info_type *scx_general_info,
185 char *blood_file);
186 int write_minc_slice(double scale, int write_byte_data,
187 int mincid, int icvid,
188 int ndims,long start[], long count[],
189 short *image, int image_size,
190 long pixel_max, float image_max,
191 double scan_time, double time_width, double zpos);
192 int get_scx_slice(scx_file *scx_fp, int slice_num,
193 long *pixel_max, float *image_max, short *image,
194 scx_file_info_type *scx_file_info,
195 scx_general_info_type *scx_general_info);
196 double decay_correction(double scan_time, double measure_time,
197 double start_time, double half_life);
198 void CreateBloodStructures (int mincHandle, int bloodHandle);
199 void FillBloodStructures (int mincHandle, int bloodHandle);
200
201 /* Constants */
202 #define TRUE 1
203 #define FALSE 0
204 #define MAX_DIMS 4
205 #define HOURS_PER_DAY 24
206 #define MIN_PER_HOUR 60
207 #define SEC_PER_MIN 60
208
209 /* Scanditronix mnemonics used */
210 #define SCX_NCS "NCS"
211 #define SCX_IMFM "IMFM"
212 #define SCX_DATY "DATY"
213 #define SCX_DATM "DATM"
214 #define SCX_DATD "DATD"
215 #define SCX_TIM "TIM"
216 #define SCX_TIMH "TIMH"
217 #define SCX_TIMM "TIMM"
218 #define SCX_TIMS "TIMS"
219 #define SCX_TIHU "TIHU"
220 #define SCX_ITM "ITM"
221 #define SCX_ITMH "ITMH"
222 #define SCX_ITMM "ITMM"
223 #define SCX_ITMS "ITMS"
224 #define SCX_MTM "MTM"
225 #define SCX_ISO "ISO"
226 #define SCX_CLV "CLV"
227 #define SCX_CDI "CDI"
228 #define SCX_CSZ "CSZ"
229 #define SCX_PXS "PXS"
230 #define SCX_MAX "MAX"
231 #define SCX_IMUN "IMUN"
232 #define SCX_MAG "MAG"
233 #define SCX_FWD "FWD"
234 #define SCX_PNM "PNM"
235 #define SCX_SEX "SEX"
236 #define SCX_AGE "AGE"
237 #define SCX_RIN "RIN"
238 #define SCX_CAR "CAR"
239 #define SCX_ACT "ACT"
240 #define SCX_IMTP "IMTP"
241 #define SCX_USC "USC"
242
243 /* Scanditronix constants */
244 #define SCX_ACTIVITY ""
245 #define SCX_MNI_GENERIC_RECONSTRUCTION_CODE "Generic 8"
246
247 /* Main program */
248
main(int argc,char * argv[])249 int main(int argc, char *argv[])
250 {
251 /* Variables for arguments */
252 static int write_byte_data=TRUE;
253 static int sort_over_time=TRUE;
254 static int clobber=FALSE;
255 static int verbose=TRUE;
256 static int decay_correct=TRUE;
257 static int slice_range[2]={0, 9999};
258 static int copy_all_header=FALSE;
259 static char *blood_file = NULL;
260
261 /* Argument option table */
262 static ArgvInfo argTable[] = {
263 {"-byte", ARGV_CONSTANT, (char *) TRUE, (char *) &write_byte_data,
264 "Write out data as bytes (default)."},
265 {"-short", ARGV_CONSTANT, (char *) FALSE, (char *) &write_byte_data,
266 "Write out data as short integers."},
267 {"-time", ARGV_CONSTANT, (char *) TRUE, (char *) &sort_over_time,
268 "Keep time ordering of data (default)."},
269 {"-zposition", ARGV_CONSTANT, (char *) FALSE, (char *) &sort_over_time,
270 "Sort data according to z position."},
271 {"-decay_correct", ARGV_CONSTANT, (char *) TRUE, (char *) &decay_correct,
272 "Do decay correction on images (default)."},
273 {"-nodecay_correct", ARGV_CONSTANT, (char *) FALSE,
274 (char *) &decay_correct, "Don't do decay correction."},
275 {"-clobber", ARGV_CONSTANT, (char *) TRUE, (char *) &clobber,
276 "Overwrite existing file."},
277 {"-noclobber", ARGV_CONSTANT, (char *) FALSE, (char *) &clobber,
278 "Don't overwrite existing file (default)."},
279 {"-verbose", ARGV_CONSTANT, (char *) TRUE, (char *) &verbose,
280 "List files as they are converted (default)"},
281 {"-quiet", ARGV_CONSTANT, (char *) FALSE, (char *) &verbose,
282 "Do not list files as they are converted."},
283 {"-small_header", ARGV_CONSTANT, (char *) FALSE,
284 (char *) ©_all_header,
285 "Copy only basic header information (default)."},
286 {"-all_header", ARGV_CONSTANT, (char *) TRUE, (char *) ©_all_header,
287 "Copy all scanditronix header information."},
288 {"-slices", ARGV_INT, (char *) 2, (char *) slice_range,
289 "Range of slices to copy."},
290 {"-bloodfile", ARGV_STRING, (char *) 1, (char *) &blood_file,
291 "Insert blood data from this file."},
292 {NULL, ARGV_END, NULL, NULL, NULL}
293 };
294
295 /* Other variables */
296 char *pname;
297 char *mincfile;
298 char **scx_files;
299 int num_scx_files;
300 scx_file_info_type *scx_file_info;
301 scx_general_info_type *scx_general_info;
302 long count[MAX_DIMS], start[MAX_DIMS];
303 int ndims;
304 int mincid, icvid, varid;
305 int islice, ifile, i, slice_num;
306 long pixel_max;
307 float image_max;
308 double scale;
309 short *image;
310 scx_file *scx_fp;
311 int status;
312 char *tm_stamp;
313 double first_z, last_z, zstep;
314
315 /* Get time stamp */
316 tm_stamp = time_stamp(argc, argv);
317
318 /* Check arguments */
319 pname = argv[0];
320 if (ParseArgv(&argc, argv, argTable, 0) || (argc < 3)) {
321 usage_error(pname);
322 }
323
324 /* Check the slice range */
325 if ((slice_range[0] < 0) || (slice_range[1] < 0) ||
326 (slice_range[1] < slice_range[0])) {
327 (void) fprintf(stderr, "%s: Error in slice range: %d to %d.\n",
328 pname, slice_range[0], slice_range[1]);
329 exit(EXIT_FAILURE);
330 }
331
332 /* Get file names */
333 mincfile = argv[argc-1];
334 argv[argc-1] = NULL; /* Null-terminate scx file list */
335 num_scx_files = argc - 2;
336 scx_files = &argv[1];
337
338 /* Print log message */
339 if (verbose) {
340 (void) fprintf(stderr, "Reading headers.\n");
341 }
342
343 /* Read the files to get basic information */
344 scx_general_info=MALLOC(1, *scx_general_info);
345 scx_file_info = MALLOC(num_scx_files, *scx_file_info);
346 status=get_scx_file_info(num_scx_files, scx_files, slice_range,
347 scx_file_info, scx_general_info);
348 if (status >= 0) {
349 (void) fprintf(stderr, "%s: Error reading scanditronix file %s.\n",
350 pname, scx_files[status]);
351 exit(EXIT_FAILURE);
352 }
353
354 /* Allocate space for ordered slice list if sorting over z position */
355 if (!sort_over_time) {
356 for (ifile=0; ifile<num_scx_files; ifile++) {
357 scx_file_info[ifile].ordered_slices =
358 MALLOC(scx_file_info[ifile].nslices, int);
359 }
360 }
361
362 /* Sort the slices */
363 sort_scx_slices(sort_over_time, num_scx_files,
364 scx_file_info, scx_general_info);
365
366 /* Setup the minc file using the first scanditronix file */
367 if (sort_over_time && (num_scx_files>1)) {
368 ndims = 4;
369 count[0]=num_scx_files;
370 count[1]=scx_general_info->max_nslices;
371 }
372 else {
373 ndims=3;
374 count[0]=scx_general_info->num_scx_slices;
375 }
376 count[ndims-1] = count[ndims-2] = scx_general_info->max_size;
377 mincid = micreate(mincfile, (clobber ? NC_CLOBBER : NC_NOCLOBBER));
378 (void) miattputstr(mincid, NC_GLOBAL, MIhistory, tm_stamp);
379 icvid=setup_minc_file(mincid, write_byte_data, copy_all_header,
380 ndims, count, num_scx_files,
381 scx_file_info, scx_general_info, blood_file);
382 if (icvid==MI_ERROR) {
383 (void) fprintf(stderr,
384 "%s: Error setting up minc file %s from scx file %s.\n",
385 pname, mincfile, scx_files[0]);
386 exit(EXIT_FAILURE);
387 }
388
389 /* Initialize minc start and count vectors */
390 for (i=0; i<ndims-2; i++) count[i]=1;
391 (void) miset_coords(ndims, (long) 0, start);
392
393 /* Set up values for decay correction */
394 scale = 1.0;
395
396 /* Allocate an image buffer */
397 image = MALLOC(scx_general_info->max_size * scx_general_info->max_size,
398 short);
399
400 /* Print log message */
401 if (verbose) {
402 (void) fprintf(stderr, "Copying files:\n");
403 }
404
405 /* Loop through files */
406 for (ifile=0; ifile<num_scx_files; ifile++) {
407
408 /* Open the scanditronix file */
409 if ((scx_fp=scx_open(scx_files[ifile]))==NULL) {
410 (void) fprintf(stderr, "%s: Error re-opening file %s.\n",
411 pname, scx_files[0]);
412 exit(EXIT_FAILURE);
413 }
414
415 /* Print log message */
416 if (verbose) {
417 (void) fprintf(stderr, " %s\n", scx_files[ifile]);
418 }
419
420 /* Get decay correction (scan time is already measured from injection
421 time) */
422 if (decay_correct &&
423 (strcmp(scx_file_info[ifile].image_type, SCX_ACTIVITY) == 0))
424 scale = decay_correction(scx_file_info[ifile].scan_time,
425 scx_file_info[ifile].time_width,
426 0.0,
427 scx_file_info[ifile].half_life);
428 else
429 scale = 1.0;
430
431 /* Loop through slices */
432 for (islice = 0; islice < scx_file_info[ifile].nslices; islice++) {
433
434 /* Set the minc coordinates */
435 if (sort_over_time && (num_scx_files>1)) {
436 start[0]=scx_file_info[ifile].ordered_file;
437 start[1]=islice;
438 }
439 else if (sort_over_time)
440 start[0]=islice;
441 else
442 start[0]=scx_file_info[ifile].ordered_slices[islice];
443 count[ndims-1] = count[ndims-2] = scx_file_info[ifile].image_size;
444
445 /* Copy the slice */
446 slice_num = islice + scx_file_info[ifile].low_slice;
447 if (get_scx_slice(scx_fp, slice_num, &pixel_max, &image_max,
448 image, &scx_file_info[ifile], scx_general_info) ||
449 write_minc_slice(scale, write_byte_data,
450 mincid, icvid, ndims, start, count,
451 image, scx_file_info[ifile].image_size,
452 pixel_max, image_max,
453 scx_file_info[ifile].scan_time,
454 scx_file_info[ifile].time_width,
455 scx_file_info[ifile].zstep * (double) slice_num +
456 scx_file_info[ifile].zstart)) {
457 (void) fprintf(stderr, "%s: Error copying slice from file %s.\n",
458 pname, scx_files[ifile]);
459 exit(EXIT_FAILURE);
460 }
461
462 } /* End slice loop */
463
464 /* Close the scanditronix file */
465 scx_close(scx_fp);
466
467 } /* End file loop */
468
469 /* Write out average z step and start for irregularly spaced slices */
470 if ((ndims!=MAX_DIMS) && (num_scx_files>1)) {
471 start[0] = 0;
472 varid = ncvarid(mincid, MIzspace);
473 (void) mivarget1(mincid, varid, start, NC_DOUBLE, NULL,
474 &first_z);
475 start[0] = scx_general_info->num_scx_slices - 1;
476 (void) mivarget1(mincid, varid, start, NC_DOUBLE, NULL,
477 &last_z);
478 if (start[0] > 0)
479 zstep = (last_z - first_z) / ((double) start[0]);
480 else
481 zstep = 1.0;
482 (void) miattputdbl(mincid, varid, MIstep, zstep);
483 (void) miattputdbl(mincid, varid, MIstart, first_z);
484 }
485
486 /* Close minc file */
487 (void) miattputstr(mincid, ncvarid(mincid, MIimage), MIcomplete, MI_TRUE);
488 (void) miclose(mincid);
489
490 FREE(image);
491 if (!sort_over_time) {
492 for (ifile=0; ifile<num_scx_files; ifile++) {
493 FREE(scx_file_info[ifile].ordered_slices);
494 }
495 }
496 FREE(scx_file_info);
497 FREE(scx_general_info);
498
499 exit(EXIT_SUCCESS);
500
501 }
502
503 /* ----------------------------- MNI Header -----------------------------------
504 @NAME : usage_error
505 @INPUT : progname - program name
506 @OUTPUT : (none)
507 @RETURNS : (nothing)
508 @DESCRIPTION: Prints a usage error message and exits.
509 @METHOD :
510 @GLOBALS :
511 @CALLS :
512 @CREATED : January 11, 1993 (Peter Neelin)
513 @MODIFIED :
514 ---------------------------------------------------------------------------- */
usage_error(char * progname)515 void usage_error(char *progname)
516 {
517 (void) fprintf(stderr,
518 "\nUsage: %s [<options>] <infile> [...] <outfile>\n", progname);
519 (void) fprintf(stderr,
520 " %s [-help]\n\n", progname);
521
522 exit(EXIT_FAILURE);
523 }
524
525 /* ----------------------------- MNI Header -----------------------------------
526 @NAME : get_scx_file_info
527 @INPUT : num_scx_files - number of scanditronix files.
528 scx_files - array of scanditronix file names.
529 slice_range - 2-component array giving range of slices
530 @OUTPUT : scx_file_info - array of structures containing information
531 about each file.
532 scx_general_info - general information about the scx files.
533 @RETURNS : (-1) if no error occurs, otherwise, the index (in scx_file)
534 of the first file that could not be read.
535 @DESCRIPTION: Reads information from each scanditronix file given by scx_files.
536 @METHOD :
537 @GLOBALS :
538 @CALLS :
539 @CREATED : January 12, 1993 (Peter Neelin)
540 @MODIFIED :
541 ---------------------------------------------------------------------------- */
get_scx_file_info(int num_scx_files,char ** scx_files,int slice_range[2],scx_file_info_type * scx_file_info,scx_general_info_type * scx_general_info)542 int get_scx_file_info(int num_scx_files, char **scx_files,
543 int slice_range[2],
544 scx_file_info_type *scx_file_info,
545 scx_general_info_type *scx_general_info)
546 {
547 static char *the_months[]=
548 {NULL, "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul",
549 "Aug", "Sep", "Oct", "Nov", "Dec"};
550 int *num_scx_slices, *max_nslices, *max_size;
551 scx_file *fp;
552 scx_file_info_type *sfip;
553 long lvalue, timh, timm, tims, tihu;
554 float fvalue;
555 char svalue[40];
556 double inj_time;
557 int ifile, i, imnem, imult, length;
558 scx_mnem_types mtype;
559 scx_mnem_list_type *mnem_list;
560 void *att_vector;
561
562 /* Initialize number of scanditronix slices */
563 num_scx_slices = &(scx_general_info->num_scx_slices);
564 max_nslices = &(scx_general_info->max_nslices);
565 max_size = &(scx_general_info->max_size);
566 *num_scx_slices = 0;
567 *max_nslices = 0;
568 *max_size = 0;
569
570 /* Loop through files, reading information */
571 for (ifile=0; ifile<num_scx_files; ifile++) {
572
573 /* Open file */
574 fp=scx_open(scx_files[ifile]);
575 if (fp==NULL) return ifile;
576 sfip = &scx_file_info[ifile];
577
578 /* Get number of slices */
579 if (scx_get_mnem(fp, SCX_NCS, 0, &lvalue, NULL, NULL)) return ifile;
580 sfip->low_slice = 0;
581 sfip->high_slice = lvalue - 1;
582 if (slice_range[0] > sfip->low_slice)
583 sfip->low_slice = slice_range[0];
584 if (slice_range[1] < sfip->high_slice)
585 sfip->high_slice = slice_range[1];
586 if (sfip->low_slice > sfip->high_slice)
587 sfip->low_slice = sfip->high_slice;
588 sfip->nslices = sfip->high_slice - sfip->low_slice + 1;
589 *num_scx_slices += sfip->nslices;
590 if (sfip->nslices > *max_nslices)
591 *max_nslices = sfip->nslices;
592
593 /* Get image width */
594 if (scx_get_mnem(fp, SCX_IMFM, 0, &lvalue, NULL, NULL)) return ifile;
595 sfip->image_size = lvalue;
596 if (lvalue > *max_size) *max_size = lvalue;
597
598 /* Get time (in seconds) */
599 if (scx_get_mnem(fp, SCX_TIMH, 0, &timh, NULL, NULL) ||
600 scx_get_mnem(fp, SCX_TIMM, 0, &timm, NULL, NULL) ||
601 scx_get_mnem(fp, SCX_TIMS, 0, &tims, NULL, NULL) ||
602 scx_get_mnem(fp, SCX_TIHU, 0, &tihu, NULL, NULL)) {
603 return ifile;
604 }
605 sfip->scan_time = (timh * MIN_PER_HOUR + timm) * SEC_PER_MIN
606 + tims + ((double) tihu)/100.0;
607
608 /* Get injection time (in seconds) (use date from before and check
609 for time before inj_time) */
610 if (scx_get_mnem(fp, SCX_ITMH, 0, &timh, NULL, NULL) ||
611 scx_get_mnem(fp, SCX_ITMM, 0, &timm, NULL, NULL) ||
612 scx_get_mnem(fp, SCX_ITMS, 0, &tims, NULL, NULL)) {
613 return ifile;
614 }
615 inj_time = (timh * MIN_PER_HOUR + timm) * SEC_PER_MIN + tims;
616
617 /* Set scan time to be from injection time. Check for scans over
618 midnight (cannot deal with scans longer than 1 day */
619 sfip->scan_time -= inj_time;
620 if (sfip->scan_time < 0.0)
621 sfip->scan_time += HOURS_PER_DAY * MIN_PER_HOUR * SEC_PER_MIN;
622
623 /* Get length of frame (in seconds) */
624 if (scx_get_mnem(fp, SCX_MTM, 0, NULL, &fvalue, NULL))
625 return ifile;
626 sfip->time_width = fvalue;
627
628 /* Get scan type */
629 if (scx_get_mnem(fp, SCX_IMTP, 0, NULL, NULL, sfip->image_type))
630 return ifile;
631
632 /* Get isotope and half-life */
633 if (scx_get_mnem(fp, SCX_ISO, 0, NULL, NULL,
634 sfip->isotope))
635 return ifile;
636 for (i=0; isotope_list[i].name !=NULL; i++) {
637 if (strncmp(isotope_list[i].name, sfip->isotope,
638 strlen(isotope_list[i].name))==0)
639 break;
640 }
641 sfip->half_life = isotope_list[i].half_life;
642
643 /* Get z start and step (correct start for non-zero first slice */
644 if (scx_get_mnem(fp, SCX_CDI, 0, NULL, &fvalue, NULL)) return ifile;
645 sfip->zstep = fvalue;
646 if (scx_get_mnem(fp, SCX_CLV, 0, NULL, &fvalue, NULL)) return ifile;
647 sfip->zstart = fvalue + sfip->low_slice * sfip->zstep;
648
649 /* Get general information from first scx file */
650 if (ifile==0) {
651 if (scx_get_mnem(fp, SCX_CSZ, 0,
652 NULL, &scx_general_info->zwidth, NULL))
653 return MI_ERROR;
654 if (scx_get_mnem(fp, SCX_PXS, 0,
655 NULL, &scx_general_info->xystep, NULL))
656 return MI_ERROR;
657 if (scx_get_mnem(fp, SCX_FWD, 0,
658 NULL, &scx_general_info->xywidth, NULL))
659 return MI_ERROR;
660 if (scx_get_mnem(fp, SCX_MAX, 0,
661 &scx_general_info->vmax, NULL, NULL))
662 return MI_ERROR;
663 if ((scx_general_info->vmax<=0) || (scx_general_info->vmax>32767))
664 scx_general_info->vmax=32000;
665 if (scx_get_mnem(fp, SCX_IMUN, 0,
666 NULL, NULL, scx_general_info->img_units))
667 return MI_ERROR;
668 if (scx_get_mnem(fp, SCX_PNM, 0,
669 NULL, NULL, scx_general_info->patient_name))
670 return MI_ERROR;
671 if (scx_get_mnem(fp, SCX_SEX, 0,
672 NULL, NULL, scx_general_info->patient_sex))
673 return MI_ERROR;
674 if (scx_general_info->patient_sex[0]=='M')
675 (void) strcpy(scx_general_info->patient_sex, MI_MALE);
676 else if (scx_general_info->patient_sex[0]=='F')
677 (void) strcpy(scx_general_info->patient_sex, MI_FEMALE);
678 else
679 (void) strcpy(scx_general_info->patient_sex, MI_OTHER);
680 if (scx_get_mnem(fp, SCX_AGE, 0,
681 &scx_general_info->patient_age, NULL, NULL))
682 return MI_ERROR;
683 if (scx_get_mnem(fp, SCX_RIN, 0,
684 NULL, NULL, scx_general_info->study_id))
685 return MI_ERROR;
686 if (scx_get_mnem(fp, SCX_DATY, 0,
687 &scx_general_info->start_year, NULL, NULL))
688 return MI_ERROR;
689 scx_general_info->start_year += 1900;
690 if (scx_general_info->start_year < 1950)
691 scx_general_info->start_year += 100;
692 if (scx_get_mnem(fp, SCX_DATM, 0,
693 &scx_general_info->start_month, NULL, NULL))
694 return MI_ERROR;
695 if (scx_get_mnem(fp, SCX_DATD, 0,
696 &scx_general_info->start_day, NULL, NULL))
697 return MI_ERROR;
698 if (scx_get_mnem(fp, SCX_TIMH, 0,
699 &scx_general_info->start_hour, NULL, NULL))
700 return MI_ERROR;
701 if (scx_get_mnem(fp, SCX_TIMM, 0,
702 &scx_general_info->start_minute, NULL, NULL))
703 return MI_ERROR;
704 if (scx_get_mnem(fp, SCX_TIMS, 0,
705 NULL, &scx_general_info->start_seconds, NULL))
706 return MI_ERROR;
707 if (scx_get_mnem(fp, SCX_TIHU, 0,
708 NULL, &fvalue, NULL))
709 return MI_ERROR;
710 scx_general_info->start_seconds += fvalue/100.0;
711 if (scx_get_mnem(fp, SCX_TIM, 0, NULL, NULL, svalue))
712 return MI_ERROR;
713 (void) sprintf(scx_general_info->start_time, "%d-%s-%d %s",
714 (int) scx_general_info->start_day,
715 the_months[scx_general_info->start_month],
716 (int) scx_general_info->start_year,
717 svalue);
718 if (scx_get_mnem(fp, SCX_CAR, 0,
719 NULL, NULL, scx_general_info->tracer))
720 return MI_ERROR;
721 if (scx_get_mnem(fp, SCX_ITM, 0,
722 NULL, NULL, scx_general_info->injection_time))
723 return MI_ERROR;
724 if (scx_get_mnem(fp, SCX_ITMH, 0,
725 &scx_general_info->injection_hour, NULL, NULL))
726 return MI_ERROR;
727 if (scx_get_mnem(fp, SCX_ITMM, 0,
728 &scx_general_info->injection_minute, NULL, NULL))
729 return MI_ERROR;
730 if (scx_get_mnem(fp, SCX_ITMS, 0,
731 NULL, &scx_general_info->injection_seconds, NULL))
732 return MI_ERROR;
733 if (scx_get_mnem(fp, SCX_ACT, 0,
734 NULL, &scx_general_info->injection_dose, NULL))
735 return MI_ERROR;
736
737 /* Find out if MNI generic reconstruction was used */
738 if (scx_get_mnem(fp, SCX_USC, 0, NULL, NULL, svalue))
739 return MI_ERROR;
740 scx_general_info->used_MNI_generic_reconstruction =
741 (strncmp(svalue, SCX_MNI_GENERIC_RECONSTRUCTION_CODE,
742 strlen(SCX_MNI_GENERIC_RECONSTRUCTION_CODE)) == 0);
743
744 /* Get list of header values */
745
746 /* Get space for first mnemonic */
747 mnem_list = MALLOC(1, *mnem_list);
748
749 /* Loop through mnemonics */
750 for (imnem=0;
751 scx_list_mnems(fp, imnem, mnem_list[imnem].name,
752 &mnem_list[imnem].mult, &mtype)!=NULL;
753 imnem++){
754
755 /* Get space for attributes (handle strings differently) */
756 switch (mtype) {
757 case scx_string:
758 mnem_list[imnem].type = NC_CHAR;
759 length = 2;
760 att_vector = MALLOC(length, char);
761 *((char *) att_vector) = '\0';
762 break;
763 case scx_long:
764 mnem_list[imnem].type = NC_LONG;
765 att_vector = MALLOC(mnem_list[imnem].mult, long);
766 break;
767 case scx_float:
768 mnem_list[imnem].type = NC_FLOAT;
769 att_vector = MALLOC(mnem_list[imnem].mult, float);
770 break;
771 }
772
773 /* Loop through multiplicity */
774 for (imult=0; imult < mnem_list[imnem].mult; imult++) {
775 if (scx_get_mnem(fp, mnem_list[imnem].name, imult,
776 &lvalue, &fvalue, svalue))
777 return MI_ERROR;
778 switch (mtype) {
779 case scx_string:
780 if (imult>0) {
781 *((char *) att_vector + length - 2) = '\n';
782 *((char *) att_vector + length - 1) = '\0';
783 length += 1;
784 }
785 length += strlen(svalue);
786 att_vector = REALLOC(att_vector, length, char);
787 att_vector = strcat((char *) att_vector, svalue);
788 break;
789 case scx_long:
790 *(((long *) att_vector) + imult) = lvalue;
791 break;
792 case scx_float:
793 *(((float *) att_vector) + imult) = fvalue;
794 break;
795 }
796
797 } /* Loop over multiplicity */
798
799 /* Save length of string */
800 if (mtype==scx_string) {
801 mnem_list[imnem].mult = length - 1;
802 }
803
804 /* Save pointer to attributes */
805 mnem_list[imnem].att_vector = att_vector;
806
807 /* Get space for next mnemonic */
808 mnem_list = REALLOC(mnem_list, imnem+2, *mnem_list);
809 }
810 scx_general_info->mnem_list=mnem_list;
811 scx_general_info->num_mnems=imnem;
812
813 } /* If first file */
814
815 /* Close file */
816 scx_close(fp);
817
818 } /* Loop over files */
819
820 return -1;
821 }
822
823 /* ----------------------------- MNI Header -----------------------------------
824 @NAME : sort_scx_slices
825 @INPUT : sort_over_time - boolean indicating whether sort should be
826 over time or z position.
827 num_scx_files - number of scanditronix files.
828 scx_file_info - array of scanditronix file information.
829 scx_general_info - general information about scx files.
830 @OUTPUT : scx_file_info - modified to give slice ordering information.
831 @RETURNS : (nothing)
832 @DESCRIPTION: Sorts the scanditronix slices for the output file.
833 @METHOD :
834 @GLOBALS :
835 @CALLS :
836 @CREATED : January 12, 1993 (Peter Neelin)
837 @MODIFIED :
838 ---------------------------------------------------------------------------- */
sort_scx_slices(int sort_over_time,int num_scx_files,scx_file_info_type * scx_file_info,scx_general_info_type * scx_general_info)839 void sort_scx_slices(int sort_over_time, int num_scx_files,
840 scx_file_info_type *scx_file_info,
841 scx_general_info_type *scx_general_info)
842 {
843 int ifile, islice, isort, num_sort, slice_num;
844 /* Variables for sorting */
845 scx_sort_type *scx_sort;
846 struct {
847 int file;
848 int slice;
849 } *slice_ptr;
850 scx_file_info_type *file_ptr;
851
852 /* Allocate array for sorting */
853 num_sort = (sort_over_time) ? num_scx_files :
854 scx_general_info->num_scx_slices;
855 scx_sort = MALLOC (num_sort, *scx_sort);
856
857 /* Are we sorting over time or z position */
858 if (sort_over_time) {
859 /* Go through the files */
860 for (ifile=0; ifile<num_scx_files; ifile++) {
861 scx_sort[ifile].sort_key = scx_file_info[ifile].scan_time;
862 scx_sort[ifile].sort_value = &scx_file_info[ifile];
863 }
864 }
865 /* Otherwise we are sorting over z position */
866 else {
867 /* Go through the files and slices */
868 isort=0;
869 for (ifile=0; ifile<num_scx_files; ifile++) {
870 for (islice = 0; islice < scx_file_info[ifile].nslices; islice++) {
871 slice_num = islice + scx_file_info[ifile].low_slice;
872 scx_sort[isort].sort_key = scx_file_info[ifile].zstart +
873 slice_num * scx_file_info[ifile].zstep;
874 slice_ptr = MALLOC(1, *slice_ptr);
875 slice_ptr->file = ifile;
876 slice_ptr->slice = islice;
877 scx_sort[isort].sort_value = slice_ptr;
878 isort++;
879 }
880 }
881 }
882
883 /* Sort the slices */
884 qsort(scx_sort, num_sort, sizeof(*scx_sort), sortcmp);
885
886 /* Loop through sorted list */
887 for (isort=0; isort<num_sort; isort++) {
888 if (sort_over_time) {
889 file_ptr=scx_sort[isort].sort_value;
890 file_ptr->ordered_file = isort;
891 }
892 else {
893 slice_ptr = scx_sort[isort].sort_value;
894 ifile = slice_ptr->file;
895 islice = slice_ptr->slice;
896 scx_file_info[ifile].ordered_slices[islice] = isort;
897 FREE(slice_ptr);
898 }
899 }
900
901 /* Free the sorting array */
902 FREE(scx_sort);
903
904 return;
905 }
906
907 /* ----------------------------- MNI Header -----------------------------------
908 @NAME : sortcmp
909 @INPUT : val1 - first value
910 val2 - second value
911 @OUTPUT : (none)
912 @RETURNS : 0 if values are the same, -1 if val1->sort_key < val2->sort_key
913 and +1 if val1->sort_key > val2->sort_key.
914 @DESCRIPTION: Compares two double precision values. If they are the same,
915 then return 0. If val1 < val2, return -1. If val1 > val2,
916 return +1.
917 @METHOD :
918 @GLOBALS :
919 @CALLS :
920 @CREATED : January 12, 1993 (Peter Neelin)
921 @MODIFIED :
922 ---------------------------------------------------------------------------- */
sortcmp(const void * val1,const void * val2)923 int sortcmp(const void *val1, const void *val2)
924 {
925
926 if (((scx_sort_type *)val1)->sort_key <
927 ((scx_sort_type *)val2)->sort_key)
928 return -1;
929 else if (((scx_sort_type *)val1)->sort_key >
930 ((scx_sort_type *)val2)->sort_key)
931 return 1;
932 else return 0;
933
934 }
935
936 /* ----------------------------- MNI Header -----------------------------------
937 @NAME : setup_minc_file
938 @INPUT : mincid - id of minc file
939 write_byte_data - boolean indicating whether data should be
940 written as bytes (TRUE) or shorts (FALSE).
941 copy_all_header - boolean indicating whether all of the
942 scanditronix header information should be copied or not.
943 ndims - number of dimensions for minc file
944 count - lengths of dimensions minc file
945 scx_filename - name of scanditronix file for header values
946 num_scx_files - number of scanditronix files.
947 scx_file_info - array of information on scanditronix files.
948 scx_general_info - general information about scx files
949 blood_file - name of blood file containing data to include
950 @OUTPUT : (nothing)
951 @RETURNS : Image conversion variable id or MI_ERROR if an error occurs.
952 @DESCRIPTION: Initializes the header of the minc file using information from
953 the scanditronix file and the other structures.
954 @METHOD :
955 @GLOBALS :
956 @CALLS :
957 @CREATED : January 12, 1993 (Peter Neelin)
958 @MODIFIED :
959 ---------------------------------------------------------------------------- */
setup_minc_file(int mincid,int write_byte_data,int copy_all_header,int ndims,long count[],int num_scx_files,scx_file_info_type * scx_file_info,scx_general_info_type * scx_general_info,char * blood_file)960 int setup_minc_file(int mincid, int write_byte_data, int copy_all_header,
961 int ndims, long count[], int num_scx_files,
962 scx_file_info_type *scx_file_info,
963 scx_general_info_type *scx_general_info,
964 char *blood_file)
965 {
966 static char *dim_names_array[]={MItime, MIzspace, MIyspace, MIxspace};
967 char **dim_names;
968 static char *dimwidth_names_array[]={
969 MItime_width, MIzspace_width, MIyspace_width, MIxspace_width};
970 char **dimwidth_names;
971 int dim[MAX_DIMS];
972 int img, imgmax, imgmin, dimvarid, widvarid, icv, varid, scx_var;
973 int idim, imnem;
974 int bloodid;
975 double vrange[2];
976
977 /* Create the dimensions */
978 dim_names = dim_names_array + MAX_DIMS - ndims;
979 dimwidth_names = dimwidth_names_array + MAX_DIMS - ndims;
980 for (idim=0; idim<ndims; idim++) {
981 dim[idim]=ncdimdef(mincid, dim_names[idim], count[idim]);
982 dimvarid=micreate_std_variable(mincid, dim_names[idim], NC_DOUBLE,
983 ((idim==0) && (num_scx_files>1)) ? 1 : 0,
984 &dim[idim]);
985 widvarid=micreate_std_variable(mincid, dimwidth_names[idim], NC_DOUBLE,
986 (strcmp(dim_names[idim], MItime)==0) ?
987 1 : 0,
988 &dim[idim]);
989
990 /* Add attributes to the dimension variables */
991 if (strcmp(dim_names[idim], MIzspace)==0) {
992 /* Write out step and start. We will rewrite this for irregularly
993 spaced files */
994 (void) miattputdbl(mincid, dimvarid, MIstep,
995 scx_file_info[0].zstep);
996 (void) miattputdbl(mincid, dimvarid, MIstart,
997 scx_file_info[0].zstart);
998 (void) miattputstr(mincid, dimvarid, MIunits, "mm");
999 (void) miattputstr(mincid, dimvarid, MIspacetype, MI_NATIVE);
1000 (void) miattputdbl(mincid, widvarid, MIwidth,
1001 (double) scx_general_info->zwidth);
1002 (void) miattputstr(mincid, widvarid, MIunits, "mm");
1003 (void) miattputstr(mincid, widvarid, MIfiltertype, MI_GAUSSIAN);
1004 }
1005 else if ((strcmp(dim_names[idim], MIyspace)==0) ||
1006 (strcmp(dim_names[idim], MIxspace)==0)) {
1007 (void) miattputstr(mincid, dimvarid, MIunits, "mm");
1008 (void) miattputdbl(mincid, dimvarid, MIstep,
1009 (double) scx_general_info->xystep);
1010 (void) miattputstr(mincid, dimvarid, MIspacetype, MI_NATIVE);
1011 (void) miattputdbl(mincid, widvarid, MIwidth,
1012 (double) scx_general_info->xywidth);
1013 (void) miattputstr(mincid, widvarid, MIfiltertype, MI_GAUSSIAN);
1014 }
1015 else if (strcmp(dim_names[idim], MItime)==0) {
1016 (void) miattputstr(mincid, dimvarid, MIunits, "seconds");
1017 (void) miattputstr(mincid, widvarid, MIunits, "seconds");
1018 }
1019 }
1020
1021 /* Create the image variable */
1022 if (write_byte_data) {
1023 img=micreate_std_variable(mincid, MIimage, NC_BYTE, ndims, dim);
1024 (void) miattputstr(mincid, img, MIsigntype, MI_UNSIGNED);
1025 vrange[0]=0; vrange[1]=255;
1026 }
1027 else {
1028 img=micreate_std_variable(mincid, MIimage, NC_SHORT, ndims, dim);
1029 (void) miattputstr(mincid, img, MIsigntype, MI_SIGNED);
1030 vrange[0] = -scx_general_info->vmax;
1031 vrange[1] = scx_general_info->vmax;
1032 }
1033 (void) ncattput(mincid, img, MIvalid_range, NC_DOUBLE, 2, vrange);
1034 (void) miattputstr(mincid, img, MIcomplete, MI_FALSE);
1035
1036 /* Create the image max and min variables */
1037 imgmax=micreate_std_variable(mincid, MIimagemax, NC_DOUBLE, ndims-2, dim);
1038 imgmin=micreate_std_variable(mincid, MIimagemin, NC_DOUBLE, ndims-2, dim);
1039 (void) miattputstr(mincid, imgmax, MIunits,
1040 scx_general_info->img_units);
1041 (void) miattputstr(mincid, imgmin, MIunits,
1042 scx_general_info->img_units);
1043
1044 /* Create the image conversion variable */
1045 icv=miicv_create();
1046 (void) miicv_setint(icv, MI_ICV_TYPE, NC_SHORT);
1047 (void) miicv_setdbl(icv, MI_ICV_VALID_MAX,
1048 (double) scx_general_info->vmax);
1049 (void) miicv_setdbl(icv, MI_ICV_VALID_MIN,
1050 (double) -scx_general_info->vmax);
1051
1052 /* Save patient info */
1053 varid = micreate_group_variable(mincid, MIpatient);
1054 (void) miattputstr(mincid, varid, MIfull_name,
1055 scx_general_info->patient_name);
1056 (void) miattputstr(mincid, varid, MIsex,
1057 scx_general_info->patient_sex);
1058 (void) ncattput(mincid, varid, MIage, NC_LONG, 1,
1059 &scx_general_info->patient_age);
1060
1061 /* Save study info */
1062 varid = micreate_group_variable(mincid, MIstudy);
1063 (void) miattputstr(mincid, varid, MImodality, MI_PET);
1064 (void) miattputstr(mincid, varid, MImanufacturer, "Scanditronix");
1065 (void) miattputstr(mincid, varid, MIstudy_id,
1066 scx_general_info->study_id);
1067 (void) miattputstr(mincid, varid, MIstart_time,
1068 scx_general_info->start_time);
1069 (void) ncattput(mincid, varid, MIstart_year, NC_LONG, 1,
1070 &scx_general_info->start_year);
1071 (void) ncattput(mincid, varid, MIstart_month, NC_LONG, 1,
1072 &scx_general_info->start_month);
1073 (void) ncattput(mincid, varid, MIstart_day, NC_LONG, 1,
1074 &scx_general_info->start_day);
1075 (void) ncattput(mincid, varid, MIstart_hour, NC_LONG, 1,
1076 &scx_general_info->start_hour);
1077 (void) ncattput(mincid, varid, MIstart_minute, NC_LONG, 1,
1078 &scx_general_info->start_minute);
1079 (void) ncattput(mincid, varid, MIstart_seconds, NC_FLOAT, 1,
1080 &scx_general_info->start_seconds);
1081
1082 /* Save acquisition info */
1083 varid = micreate_group_variable(mincid, MIacquisition);
1084 (void) miattputstr(mincid, varid, MIradionuclide,
1085 scx_file_info[0].isotope);
1086 if (scx_file_info[0].half_life > 0.0) {
1087 (void) miattputdbl(mincid, varid, MIradionuclide_halflife,
1088 (double) scx_file_info[0].half_life);
1089 }
1090 (void) miattputstr(mincid, varid, MItracer,
1091 scx_general_info->tracer);
1092 (void) miattputstr(mincid, varid, MIinjection_time,
1093 scx_general_info->injection_time);
1094 (void) ncattput(mincid, varid, MIinjection_hour, NC_LONG, 1,
1095 &scx_general_info->injection_hour);
1096 (void) ncattput(mincid, varid, MIinjection_minute, NC_LONG, 1,
1097 &scx_general_info->injection_minute);
1098 (void) ncattput(mincid, varid, MIinjection_seconds, NC_FLOAT, 1,
1099 &scx_general_info->injection_seconds);
1100 (void) ncattput(mincid, varid, MIinjection_dose, NC_FLOAT, 1,
1101 &scx_general_info->injection_dose);
1102 (void) miattputstr(mincid, varid, MIdose_units, "mCurie");
1103
1104 /* If the MNI generic reconstruction was used, indicate the fact */
1105 if (scx_general_info->used_MNI_generic_reconstruction) {
1106 varid = ncvardef(mincid, "MNI_PET_rec", NC_LONG, 0, NULL);
1107 (void) miattputstr(mincid, varid, MIvartype, MI_GROUP);
1108 (void) miattputstr(mincid, varid, MIvarid,
1109 "MNI PET reconstruction info");
1110 (void) miadd_child(mincid, ncvarid(mincid, MIrootvariable), varid);
1111 (void) miattputstr(mincid, varid, "generic_reconstruction", MI_TRUE);
1112 }
1113
1114 /* If we want all of the values from the scanditronix header, get them */
1115 if (copy_all_header) {
1116
1117 /* Create a variable for scx mnemonics */
1118 scx_var = ncvardef(mincid, "scanditronix", NC_LONG, 0, NULL);
1119 (void) miattputstr(mincid, scx_var, MIvartype, MI_GROUP);
1120 (void) miattputstr(mincid, scx_var, MIvarid, "MNI SCX variable");
1121 (void) miadd_child(mincid, ncvarid(mincid, MIrootvariable), scx_var);
1122
1123 /* Loop through mnemonics */
1124 for (imnem=0; imnem<scx_general_info->num_mnems; imnem++){
1125 ncattput(mincid, scx_var,
1126 scx_general_info->mnem_list[imnem].name,
1127 scx_general_info->mnem_list[imnem].type,
1128 scx_general_info->mnem_list[imnem].mult,
1129 scx_general_info->mnem_list[imnem].att_vector);
1130 } /* Loop through mnemonics */
1131 } /* If copy_all_header */
1132
1133 /* Open the blood file and create the variables if needed */
1134 if (blood_file != NULL) {
1135 bloodid = ncopen(blood_file, NC_NOWRITE);
1136 CreateBloodStructures(mincid, bloodid);
1137 }
1138
1139 /* Attach the icv */
1140 (void) ncendef(mincid);
1141 (void) miicv_attach(icv, mincid, img);
1142
1143 /* Copy the blood data */
1144 if (blood_file != NULL) {
1145 FillBloodStructures(mincid, bloodid);
1146 ncclose(bloodid);
1147 }
1148
1149 return icv;
1150 }
1151
1152 /* ----------------------------- MNI Header -----------------------------------
1153 @NAME : get_scx_slice
1154 @INPUT : scx_fp - file pointer for scanditronix file
1155 slice_num - slice to copy
1156 scx_file_info - information on scanditronix file.
1157 scx_general_info - general scx file information
1158 @OUTPUT : pixel_max - maximum pixel value
1159 image_max - real value to which pixel_max corresponds
1160 image - scanditronix image
1161 @RETURNS : Returns TRUE if an error occurs.
1162 @DESCRIPTION: Gets a scanditronix image.
1163 @METHOD :
1164 @GLOBALS :
1165 @CALLS :
1166 @CREATED : January 20, 1993 (Peter Neelin)
1167 @MODIFIED :
1168 ---------------------------------------------------------------------------- */
get_scx_slice(scx_file * scx_fp,int slice_num,long * pixel_max,float * image_max,short * image,scx_file_info_type * scx_file_info,scx_general_info_type * scx_general_info)1169 int get_scx_slice(scx_file *scx_fp, int slice_num,
1170 long *pixel_max, float *image_max, short *image,
1171 scx_file_info_type *scx_file_info,
1172 scx_general_info_type *scx_general_info)
1173 /* ARGSUSED */
1174 {
1175 long npix, ix, iy, y_offset, off1, off2;
1176 short temp;
1177
1178 /* Get the image from the scanditronix file */
1179 if (scx_get_image(scx_fp, slice_num, image)) return TRUE;
1180
1181 /* Flip the scanditronix image to give positive x & y axes */
1182 npix = scx_file_info->image_size * scx_file_info->image_size;
1183 for (iy=0; iy<scx_file_info->image_size/2; iy++) {
1184 y_offset = iy * scx_file_info->image_size;
1185 for (ix=0; ix<scx_file_info->image_size; ix++) {
1186 off1 = y_offset + ix;
1187 off2 = npix - off1 - 1;
1188 temp = image[off1];
1189 image[off1] = image[off2];
1190 image[off2] = temp;
1191 }
1192 }
1193
1194 /* Get image and pixel max */
1195 if (scx_get_mnem(scx_fp, SCX_MAG, slice_num, NULL, image_max, NULL))
1196 return TRUE;
1197 if (scx_get_mnem(scx_fp, SCX_MAX, 0, pixel_max, NULL, NULL)) return TRUE;
1198 if ((*pixel_max<=0) || (*pixel_max>32767)) *pixel_max = 32000;
1199
1200 return FALSE;
1201 }
1202
1203 /* ----------------------------- MNI Header -----------------------------------
1204 @NAME : write_minc_slice
1205 @INPUT : scale - scale for decay correcting image
1206 write_byte_data - boolean indicating whether data should be
1207 written as bytes (TRUE) or shorts (FALSE).
1208 mincid - id of minc file
1209 icvid - id of image conversion variable
1210 start - coordinate of slice in minc file
1211 count - edge lengths of image to write in minc file
1212 image - pointer to image buffer
1213 pixel_max - maximum pixel value
1214 image_max - real value to which pixel_max corresponds
1215 scan_time - time of slice
1216 time_width - time width of slice
1217 zpos - z position of slice
1218 scx_file_info - information on scanditronix file.
1219 scx_general_info - general scx file information
1220 @OUTPUT : (nothing)
1221 @RETURNS : Returns TRUE if an error occurs.
1222 @DESCRIPTION: Writes out the image to the minc file.
1223 @METHOD :
1224 @GLOBALS :
1225 @CALLS :
1226 @CREATED : January 12, 1993 (Peter Neelin)
1227 @MODIFIED :
1228 ---------------------------------------------------------------------------- */
write_minc_slice(double scale,int write_byte_data,int mincid,int icvid,int ndims,long start[],long count[],short * image,int image_size,long pixel_max,float image_max,double scan_time,double time_width,double zpos)1229 int write_minc_slice(double scale, int write_byte_data,
1230 int mincid, int icvid,
1231 int ndims,long start[], long count[],
1232 short *image, int image_size,
1233 long pixel_max, float image_max,
1234 double scan_time, double time_width, double zpos)
1235 /*ARGSUSED*/
1236 {
1237 double pixmin, pixmax;
1238 long ipix, npix;
1239 double maximum, minimum;
1240
1241 /* Search for pixel max and min */
1242 npix = image_size * image_size;
1243 pixmin = pixmax = image[0];
1244 for (ipix=1; ipix<npix; ipix++) {
1245 if (image[ipix]>pixmax) pixmax = image[ipix];
1246 if (image[ipix]<pixmin) pixmin = image[ipix];
1247 }
1248
1249 /* Get image max and min */
1250 maximum = image_max * pixmax / (double) pixel_max;
1251 minimum = image_max * pixmin / (double)pixel_max;
1252
1253 /* Change valid range on icv */
1254 (void) miicv_detach(icvid);
1255 (void) miicv_setdbl(icvid, MI_ICV_VALID_MAX, (double) pixmax);
1256 (void) miicv_setdbl(icvid, MI_ICV_VALID_MIN, (double) pixmin);
1257 (void) miicv_attach(icvid, mincid, ncvarid(mincid, MIimage));
1258
1259 /* Calculate real max and min */
1260 maximum = maximum * scale;
1261 minimum = minimum * scale;
1262
1263 /* Write out the image */
1264 (void) miicv_put(icvid, start, count, image);
1265
1266 /* Write out image max and min */
1267 (void) ncvarput1(mincid, ncvarid(mincid, MIimagemax), start, &maximum);
1268 (void) ncvarput1(mincid, ncvarid(mincid, MIimagemin), start, &minimum);
1269
1270 /* Write out time and z position */
1271 if (ndims==MAX_DIMS) {
1272 (void) mivarput1(mincid, ncvarid(mincid, MItime), start,
1273 NC_DOUBLE, NULL, &scan_time);
1274 (void) mivarput1(mincid, ncvarid(mincid, MItime_width), start,
1275 NC_DOUBLE, NULL, &time_width);
1276 }
1277 (void) mivarput1(mincid, ncvarid(mincid, MIzspace), &start[ndims-3],
1278 NC_DOUBLE, NULL, &zpos);
1279
1280 return FALSE;
1281 }
1282
1283 /* ----------------------------- MNI Header -----------------------------------
1284 @NAME : decay_correction
1285 @INPUT : scan_time - time of beginning of sample
1286 measure_time - length of sample
1287 start_time - time to which we should decay correct
1288 half_life - half life of isotope
1289 @OUTPUT : (nothing)
1290 @RETURNS : decay correction scaling factor
1291 @DESCRIPTION: Calculates the decay correction needed to get equivalent counts
1292 at start_time. Correction assumes constant activity (without
1293 decay) over measure_time.
1294 @METHOD :
1295 @GLOBALS :
1296 @CALLS :
1297 @CREATED : January 21, 1993 (Peter Neelin)
1298 @MODIFIED :
1299 ---------------------------------------------------------------------------- */
decay_correction(double scan_time,double measure_time,double start_time,double half_life)1300 double decay_correction(double scan_time, double measure_time,
1301 double start_time, double half_life)
1302 {
1303 double mean_life;
1304 double measure_correction;
1305 double decay;
1306
1307 /* Check for negative half_life and calculate mean life */
1308 if (half_life <= 0.0) return 1.0;
1309 mean_life = half_life/ log(2.0);
1310
1311 /* Normalize scan time and measure_time */
1312 scan_time = (scan_time - start_time) / mean_life;
1313 measure_time /= mean_life;
1314
1315 /* Calculate correction for decay over measuring time (assuming a
1316 constant activity). Check for possible rounding errors. */
1317 if ((measure_time*measure_time/2.0) < DBL_EPSILON) {
1318 measure_correction = 1.0 - measure_time/2.0;
1319 }
1320 else {
1321 measure_correction = (1.0 - exp(-measure_time)) / fabs(measure_time);
1322 }
1323
1324 /* Calculate decay */
1325 decay = exp(-scan_time) * measure_correction;
1326 if (decay<=0.0) decay = DBL_MAX;
1327 else decay = 1.0/decay;
1328
1329 return decay;
1330 }
1331
1332