1 /* concorde microPET to minc */
2 #include "config.h"
3 
4 #define _XOPEN_SOURCE 1
5 #include <stdio.h>
6 #include <stdlib.h>
7 #include <stdarg.h>
8 #include <float.h>
9 #include <time.h>
10 #include <string.h>
11 #include <minc.h>
12 #include <time_stamp.h>
13 #include <ParseArgv.h>
14 
15 #define VERSIONSTR VERSION " built " __DATE__ " " __TIME__
16 
17 /*************************************************************************
18  * Concorde microPET definitions
19  */
20 /* file_type field */
21 #define UPET_FT_UNKNOWN 0
22 #define UPET_FT_LIST_MODE 1
23 #define UPET_FT_SINOGRAM 2
24 #define UPET_FT_NORMALIZATION 3
25 #define UPET_FT_ATTENUATION_CORRECTION 4
26 #define UPET_FT_IMAGE 5         /* Standard image data file (typical) */
27 #define UPET_FT_BLANK 6
28 #define UPET_FT_MU_MAP 8        /* Mu map data file */
29 #define UPET_FT_SCATTER_CORRECTION 9
30 
31 /* acquisition_mode field */
32 #define UPET_AM_UNKNOWN 0
33 #define UPET_AM_BLANK 1
34 #define UPET_AM_EMISSION 2
35 #define UPET_AM_DYNAMIC 3
36 #define UPET_AM_GATED 4
37 #define UPET_AM_CONTINUOUS_BED_MOTION 5
38 #define UPET_AM_SINGLES_TRANSMISSION 6
39 #define UPET_AM_WINDOWED_COINCIDENCE_TRANSMISSION 7
40 #define UPET_AM_NONWINDOWED_COINCIDENCE_TRANSMISSION 8
41 
42 #define UPET_DT_UNKNOWN 0
43 #define UPET_DT_BYTE 1
44 #define UPET_DT_II16 2           /* Intel 16-bit signed integer */
45 #define UPET_DT_II32 3           /* Intel 32-bit signed integer */
46 #define UPET_DT_IF32 4           /* Intel 32-bit float */
47 #define UPET_DT_MF32 5           /* Sun 32-bit float */
48 #define UPET_DT_MI16 6           /* Sun 16-bit signed integer */
49 #define UPET_DT_MI32 7           /* Sun 32-bit signed integer */
50 
51 #define DECLARE_FUNC(x) \
52 static int x(struct conversion_info *ci_ptr, char *val_str, char *new_var, char *new_att)
53 
54 struct conversion_info {
55     FILE *hdr_fp;
56     FILE *img_fp;
57     int mnc_fd;
58     int frame_index;
59     int frame_zero;
60     int data_type;
61     nc_type minc_type;
62     int dim_count;
63     int dim_lengths[5];
64     int dim_ids[5];
65     double dim_steps[5];
66     int frame_nbytes;
67     int frame_nvoxels;
68     void *frame_buffer;
69     double scale_factor;
70     double deadtime_correction;
71     double decay_correction;
72     double calibration_factor;
73     double isotope_branching_fraction;
74     int swap_size;              /* 0, 2, 4 */
75 };
76 
77 DECLARE_FUNC(upet_file_type);
78 DECLARE_FUNC(upet_acq_mode);
79 DECLARE_FUNC(upet_bed_motion);
80 DECLARE_FUNC(upet_data_type);
81 DECLARE_FUNC(upet_data_order);
82 DECLARE_FUNC(upet_ndims);
83 DECLARE_FUNC(upet_total_frames);
84 DECLARE_FUNC(upet_x_dim);
85 DECLARE_FUNC(upet_y_dim);
86 DECLARE_FUNC(upet_z_dim);
87 DECLARE_FUNC(upet_vector_dim);
88 DECLARE_FUNC(upet_injection_time);
89 DECLARE_FUNC(upet_scan_time);
90 DECLARE_FUNC(upet_axial_crystal_pitch);
91 DECLARE_FUNC(upet_pixel_size);
92 DECLARE_FUNC(upet_dose_units);
93 DECLARE_FUNC(upet_calibration_factor);
94 DECLARE_FUNC(upet_rotation);
95 DECLARE_FUNC(upet_isotope_branching_fraction);
96 
97 DECLARE_FUNC(upet_frame_no);
98 DECLARE_FUNC(upet_frame_start);
99 DECLARE_FUNC(upet_frame_duration);
100 DECLARE_FUNC(upet_frame_min);
101 DECLARE_FUNC(upet_frame_max);
102 DECLARE_FUNC(upet_frame_file_ptr);
103 DECLARE_FUNC(upet_frame_scale_factor);
104 DECLARE_FUNC(upet_frame_decay_correction);
105 DECLARE_FUNC(upet_frame_deadtime_correction);
106 
107 static void copy_init(struct conversion_info *ci_ptr);
108 static void copy_frame(struct conversion_info *ci_ptr);
109 
110 /* These values are used to represent the field types in the microPET
111  * header file.
112  */
113 #define UPET_TYPE_STR 1         /* String */
114 #define UPET_TYPE_INT 2         /* Integer */
115 #define UPET_TYPE_REAL 3        /* Floating-point */
116 #define UPET_TYPE_TIME 4        /* Timestamp */
117 #define UPET_TYPE_FILTER 5      /* Integer type followed by a float cutoff */
118 #define UPET_TYPE_3X64 6        /* 3 64-bit integers */
119 #define UPET_TYPE_FPTR 7        /* File pointer (2 32 bit integers) */
120 #define UPET_TYPE_SINGLE 8      /* Block #, singles/sec, raw singles/sec */
121 
122 /* concorde keywords */
123 struct keywd_entry {
124     char *upet_kwd;
125     short upet_type;
126     char *mnc_var;
127     char *mnc_att;
128     int (*func)(struct conversion_info *, char *, char *, char *);
129 };
130 
131 /* Per-frame attributes in the concorde microPET header */
132 struct keywd_entry frm_atts[] = {
133     { "frame",
134       UPET_TYPE_INT, NULL, NULL, upet_frame_no },
135     { "event_type",
136       UPET_TYPE_INT, NULL, NULL, NULL },
137     { "gate",
138       UPET_TYPE_INT, NULL, NULL, NULL },
139     { "bed",
140       UPET_TYPE_INT, NULL, NULL, NULL },
141     { "bed_offset",
142       UPET_TYPE_REAL, NULL, NULL, NULL },
143     { "ending_bed_offset",
144       UPET_TYPE_REAL, NULL, NULL, NULL },
145     { "vertical_bed_offset",
146       UPET_TYPE_REAL, NULL, NULL, NULL },
147     { "data_file_pointer",
148       UPET_TYPE_FPTR, NULL, NULL, upet_frame_file_ptr },
149     { "frame_start",
150       UPET_TYPE_REAL, NULL, NULL, upet_frame_start },
151     { "frame_duration",
152       UPET_TYPE_REAL, NULL, NULL, upet_frame_duration },
153     { "scale_factor",
154       UPET_TYPE_REAL, NULL, NULL, upet_frame_scale_factor },
155     { "minimum",
156       UPET_TYPE_REAL, NULL, NULL, upet_frame_min },
157     { "maximum",
158       UPET_TYPE_REAL, NULL, NULL, upet_frame_max },
159     { "deadtime_correction",
160       UPET_TYPE_REAL, NULL, NULL, upet_frame_deadtime_correction },
161     { "decay_correction",
162       UPET_TYPE_REAL, NULL, NULL, upet_frame_decay_correction },
163     { "prompts",
164       UPET_TYPE_3X64, NULL, NULL, NULL },
165     { "delays",
166       UPET_TYPE_3X64, NULL, NULL, NULL },
167     { "trues",
168       UPET_TYPE_3X64, NULL, NULL, NULL },
169     { "prompts_rate",
170       UPET_TYPE_INT, NULL, NULL, NULL },
171     { "delays_rate",
172       UPET_TYPE_INT, NULL, NULL, NULL },
173     { "singles",
174       UPET_TYPE_SINGLE, NULL, NULL, NULL },
175     { NULL, 0, NULL, NULL, NULL }
176 };
177 
178 /* Per-volume attributes in the concorde microPET header */
179 struct keywd_entry vol_atts[] = {
180     { "ROI_file:",
181       UPET_TYPE_STR, NULL, NULL, NULL },
182     { "version",
183       UPET_TYPE_STR, NULL, NULL, NULL },
184     { "model",
185       UPET_TYPE_INT, MIstudy, MIdevice_model, NULL },
186     { "institution",
187       UPET_TYPE_STR, MIstudy, MIinstitution, NULL },
188     { "study",
189       UPET_TYPE_STR, NULL, NULL, NULL },
190     { "file_name",
191       UPET_TYPE_STR, NULL, NULL, NULL }, /* actual file data */
192     { "file_type",
193       UPET_TYPE_INT, NULL, NULL, upet_file_type },
194     { "acquisition_mode",
195       UPET_TYPE_INT, NULL, NULL, upet_acq_mode},
196     { "bed_motion",
197       UPET_TYPE_INT, NULL, NULL, upet_bed_motion},
198     { "total_frames",
199       UPET_TYPE_INT, NULL, NULL, upet_total_frames },
200     { "isotope",
201       UPET_TYPE_STR, MIacquisition, MIradionuclide, NULL },
202     { "isotope_half_life",
203       UPET_TYPE_REAL, MIacquisition, MIradionuclide_halflife, NULL },
204     { "isotope_branching_fraction",
205       UPET_TYPE_REAL, NULL, NULL, upet_isotope_branching_fraction },
206     { "transaxial_crystals_per_block",
207       UPET_TYPE_INT, NULL, NULL, NULL },
208     { "axial_crystals_per_block",
209       UPET_TYPE_INT, NULL, NULL, NULL },
210     { "intrinsic_crystal_offset",
211       UPET_TYPE_INT, NULL, NULL, NULL },
212     { "transaxial_blocks",
213       UPET_TYPE_INT, NULL, NULL, NULL },
214     { "axial_blocks",
215       UPET_TYPE_INT, NULL, NULL, NULL },
216     { "transaxial_crystal_pitch",
217       UPET_TYPE_REAL, NULL, NULL, NULL },
218     { "axial_crystal_pitch",
219       UPET_TYPE_REAL, NULL, NULL, upet_axial_crystal_pitch },
220     { "radius",
221       UPET_TYPE_REAL, NULL, NULL, NULL },
222     { "radial_fov",
223       UPET_TYPE_REAL, NULL, NULL, NULL },
224     { "src_radius",
225       UPET_TYPE_REAL, NULL, NULL, NULL },
226     {"src_cm_per_rev",
227      UPET_TYPE_REAL, NULL, NULL, NULL },
228     {"src_steps_per_rev",
229      UPET_TYPE_INT, NULL, NULL, NULL },
230     {"tx_src_type",
231      UPET_TYPE_INT, NULL, NULL, NULL },
232     {"default_projections",
233      UPET_TYPE_INT, NULL, NULL, NULL },
234     {"default_transaxial_angles",
235      UPET_TYPE_INT, NULL, NULL, NULL },
236     {"crystal_thickness",
237      UPET_TYPE_REAL, NULL, NULL, NULL },
238     {"depth_of_interaction",
239      UPET_TYPE_REAL, NULL, NULL, NULL },
240     {"transaxial_bin_size",
241      UPET_TYPE_REAL, NULL, NULL, NULL },
242     {"axial_plane_size",
243      UPET_TYPE_REAL, NULL, NULL, NULL },
244     {"lld",
245      UPET_TYPE_REAL, NULL, NULL, NULL },
246     {"uld",
247      UPET_TYPE_REAL, NULL, NULL, NULL },
248     {"timing_window",
249      UPET_TYPE_INT, NULL, NULL, NULL },
250     {"data_type",
251      UPET_TYPE_INT, NULL, NULL, upet_data_type },
252     {"data_order",
253      UPET_TYPE_INT, NULL, NULL, upet_data_order },
254     {"span",
255      UPET_TYPE_INT, NULL, NULL, NULL },
256     {"ring_difference",
257      UPET_TYPE_INT,  NULL, NULL, NULL },
258     {"number_of_dimensions",
259      UPET_TYPE_INT, NULL, NULL, upet_ndims },
260     {"x_dimension",
261      UPET_TYPE_INT, NULL, NULL, upet_x_dim },
262     {"y_dimension",
263      UPET_TYPE_INT, NULL, NULL, upet_y_dim },
264     {"z_dimension",
265      UPET_TYPE_INT, NULL, NULL, upet_z_dim },
266     {"w_dimension",
267      UPET_TYPE_INT, NULL, NULL, upet_vector_dim },
268     {"x_filter",
269      UPET_TYPE_FILTER, NULL, NULL, NULL },
270     {"y_filter",
271      UPET_TYPE_FILTER, NULL, NULL, NULL },
272     {"z_filter",
273      UPET_TYPE_FILTER, NULL, NULL, NULL },
274     {"histogram_version",
275      UPET_TYPE_STR, NULL, NULL, NULL },
276     {"rebinning_type",
277      UPET_TYPE_INT, NULL, NULL, NULL },
278     {"rebinning_version",
279      UPET_TYPE_STR, NULL, NULL, NULL },
280     {"recon_algorithm",
281      UPET_TYPE_INT, NULL, NULL, NULL },
282     {"recon_version",
283      UPET_TYPE_STR, NULL, NULL, NULL },
284     {"map_subsets",
285      UPET_TYPE_INT, NULL, NULL, NULL },
286     {"map_osem3d_iterations",
287      UPET_TYPE_INT, NULL, NULL, NULL },
288     {"map_iterations",
289      UPET_TYPE_INT, NULL, NULL, NULL },
290     {"map_beta",
291      UPET_TYPE_REAL, NULL, NULL, NULL },
292     {"map_blur_type",
293      UPET_TYPE_INT, NULL, NULL, NULL },
294     {"map_prior_type",
295      UPET_TYPE_INT, NULL, NULL, NULL },
296     {"map_blur_file",
297      UPET_TYPE_STR, NULL, NULL, NULL },
298     {"map_pmatrix_file",
299      UPET_TYPE_STR, NULL, NULL, NULL },
300     {"deadtime_correction_applied",
301      UPET_TYPE_INT, NULL, NULL, NULL },
302     {"decay_correction_applied",
303      UPET_TYPE_INT, NULL, NULL, NULL },
304     {"normalization_applied",
305      UPET_TYPE_INT, NULL, NULL, NULL },
306     {"normalization_filename",
307      UPET_TYPE_STR, NULL, NULL, NULL },
308     {"attenuation_applied",
309      UPET_TYPE_INT, NULL, NULL, NULL },
310     {"attenuation_filename",
311      UPET_TYPE_STR, NULL, NULL, NULL },
312     {"scatter_correction",
313      UPET_TYPE_INT, NULL, NULL, NULL },
314     {"scatter_version",
315      UPET_TYPE_STR, NULL, NULL, NULL },
316     {"arc_correction_applied",
317      UPET_TYPE_INT, NULL, NULL, NULL },
318     {"rotation",
319      UPET_TYPE_REAL, NULL, NULL, upet_rotation },
320     {"x_offset",
321      UPET_TYPE_REAL, NULL, NULL, NULL },
322     {"y_offset",
323      UPET_TYPE_REAL, NULL, NULL, NULL },
324     {"z_offset",
325      UPET_TYPE_REAL, NULL, NULL, NULL },
326     {"zoom",
327      UPET_TYPE_REAL, NULL, NULL, NULL },
328     {"pixel_size",
329      UPET_TYPE_REAL, NULL, NULL, upet_pixel_size },
330     {"calibration_units",
331      UPET_TYPE_INT, NULL, NULL, NULL },
332     {"calibration_factor",
333      UPET_TYPE_REAL, NULL, NULL, upet_calibration_factor },
334     {"calibration_branching_fraction",
335      UPET_TYPE_REAL, NULL, NULL, NULL },
336     {"number_of_singles_rates",
337      UPET_TYPE_INT, NULL, NULL, NULL },
338     {"investigator",
339      UPET_TYPE_STR, MIstudy, "investigator", NULL },
340     {"operator",
341      UPET_TYPE_STR, MIstudy, MIoperator, NULL },
342     {"study_identifier",
343      UPET_TYPE_STR, MIstudy, MIstudy_id, NULL },
344     {"scan_time",
345      UPET_TYPE_TIME, NULL, NULL, upet_scan_time },
346     {"injected_compound",
347      UPET_TYPE_STR, NULL, NULL, NULL },
348     {"dose_units",
349      UPET_TYPE_INT, MIacquisition, MIdose_units, upet_dose_units },
350     {"dose",
351      UPET_TYPE_REAL, MIacquisition, MIinjection_dose, NULL },
352     {"injection_time",
353      UPET_TYPE_TIME, NULL, NULL, upet_injection_time },
354     {"injection_decay_correction",
355      UPET_TYPE_REAL, NULL, NULL, NULL },
356     {"subject_identifier",
357      UPET_TYPE_STR, NULL, NULL, NULL },
358     {"subject_genus",
359      UPET_TYPE_STR, NULL, NULL, NULL },
360     {"subject_orientation",
361      UPET_TYPE_INT, NULL, NULL, NULL },
362     {"subject_length_units",
363      UPET_TYPE_INT, NULL, NULL, NULL },
364     {"subject_length",
365      UPET_TYPE_REAL, NULL, NULL, NULL },
366     {"subject_weight_units",
367      UPET_TYPE_INT, NULL, NULL, NULL },
368     {"subject_weight",
369      UPET_TYPE_REAL, NULL, NULL, NULL },
370     {"subject_phenotype",
371      UPET_TYPE_STR, NULL, NULL, NULL },
372     {"study_model",
373      UPET_TYPE_STR, NULL, NULL, NULL },
374     {"anesthesia",
375      UPET_TYPE_STR, NULL, NULL, NULL },
376     {"analgesia",
377      UPET_TYPE_STR, NULL, NULL, NULL },
378     {"other_drugs",
379      UPET_TYPE_STR, NULL, NULL, NULL },
380     {"food_access",
381      UPET_TYPE_STR, NULL, NULL, NULL },
382     {"water_access",
383      UPET_TYPE_STR, NULL, NULL, NULL },
384     {NULL, 0, NULL, NULL, NULL }
385 };
386 
387 /* Reflects "normal" image data order */
388 
389 #define DIM_T 0
390 #define DIM_Z 1
391 #define DIM_Y 2
392 #define DIM_X 3
393 #define DIM_W 4
394 
395 static char *_dimnames[5];
396 
397 /* Calculate the overall scaling factor for the image data from the
398  * conversion information structure.
399  */
400 #define COMBINED_SCALE_FACTOR(ci_ptr) \
401      ((ci_ptr->scale_factor * ci_ptr->calibration_factor) / \
402       (ci_ptr->isotope_branching_fraction))
403 
404 
405 #define ORIENT_BODY 1
406 #define ORIENT_HEAD 2
407 
408 int _orient_flag = ORIENT_HEAD;
409 int _verbose_flag = 1;
410 
411 ArgvInfo argTable[] = {
412     {"-head", ARGV_CONSTANT, (char *) ORIENT_HEAD, (char *) &_orient_flag,
413      "Orient image for cerebral viewing (as with human brain)"},
414     {"-body", ARGV_CONSTANT, (char *) ORIENT_BODY, (char *) &_orient_flag,
415      "Orient image for whole-body viewing (Z along long axis)"},
416     {"-quiet", ARGV_CONSTANT, (char *) 0, (char *) &_verbose_flag,
417      "Turn off the various progress reporting messages."},
418      {NULL, ARGV_VERINFO, (char *) VERSIONSTR, (char *) NULL, NULL},
419      {NULL, ARGV_END, NULL, NULL, NULL}
420 };
421 
422 typedef enum {
423     MSG_INFO,
424     MSG_WARNING,
425     MSG_ERROR,
426     MSG_FATAL
427 } msg_level_t;
428 
429 static void
message(msg_level_t level,char * fmt,...)430 message(msg_level_t level, char *fmt, ...)
431 {
432     va_list ap;
433     const char *prefix_str;
434 
435     switch (level) {
436     case MSG_WARNING:
437         prefix_str = "WARNING: ";
438         break;
439     case MSG_ERROR:
440         prefix_str = "ERROR: ";
441         break;
442     case MSG_FATAL:
443         prefix_str = "FATAL: ";
444         break;
445     default:
446         prefix_str = NULL;
447         break;
448     }
449     if (_verbose_flag || level != MSG_INFO) {
450         if (level != MSG_INFO) {
451             if (prefix_str != NULL) {
452                 fprintf(stderr, "%s", prefix_str);
453             }
454             va_start(ap, fmt);
455             vfprintf(stderr, fmt, ap);
456             va_end(ap);
457         }
458         if (prefix_str != NULL) {
459             fprintf(stdout, "%s", prefix_str);
460         }
461         va_start(ap, fmt);
462         vfprintf(stdout, fmt, ap);
463         va_end(ap);
464     }
465 }
466 
467 static int
is_host_big_endian()468 is_host_big_endian()
469 {
470     long ltmp = 0x04030201;
471     char *ctmp = (char *) &ltmp;
472 
473     if (ctmp[0] == 0x01) {
474         return (0);
475     }
476     if (ctmp[0] == 0x04) {
477         return (1);
478     }
479     return (-1);
480 }
481 
482 static void
usage(const char * progname)483 usage(const char *progname)
484 {
485    fprintf(stderr, "\nUsage: %s [<options>] input.img[.hdr] [output.mnc]\n",
486            progname);
487    fprintf(stderr,   "       %s [-help]\n\n", progname);
488    exit(-1);
489 }
490 
491 int
upet_to_minc(char * hdr_fname,char * img_fname,char * out_fname,char * prog_name)492 upet_to_minc(char *hdr_fname, char *img_fname, char *out_fname,
493              char *prog_name)
494 {
495     char *line_ptr;
496     char line_buf[1024];
497     char *val_ptr;
498     int in_header;
499     double dbl_tmp;
500     int int_tmp;
501     struct conversion_info ci;
502     struct keywd_entry *ke_ptr;
503     int is_known;
504     char *argv_tmp[5];
505     char *out_history;
506 
507     ci.hdr_fp = fopen(hdr_fname, "r"); /* Text file */
508     if (ci.hdr_fp == NULL) {
509         perror(hdr_fname);
510         return (-1);
511     }
512 
513     ci.img_fp = fopen(img_fname, "rb"); /* Binary file */
514     if (ci.img_fp == NULL) {
515         perror(img_fname);
516         return (-1);
517     }
518 
519     ci.mnc_fd = micreate(out_fname, NC_NOCLOBBER);
520     if (ci.mnc_fd < 0) {
521         perror(out_fname);
522         return (-1);
523     }
524 
525     ci.frame_zero = -1;     /* Initial frame is -1 until set. */
526 
527     /* Define the basic MINC group variables.
528      */
529     micreate_group_variable(ci.mnc_fd, MIstudy);
530     micreate_group_variable(ci.mnc_fd, MIacquisition);
531     micreate_group_variable(ci.mnc_fd, MIpatient);
532     ncvardef(ci.mnc_fd, "micropet", NC_SHORT, 0, NULL);
533 
534     /* Fake the history here */
535     argv_tmp[0] = prog_name;
536     argv_tmp[1] = VERSIONSTR;
537     argv_tmp[2] = hdr_fname;
538     argv_tmp[3] = img_fname;
539     argv_tmp[4] = out_fname;
540 
541     out_history = time_stamp(5, argv_tmp);
542 
543     miattputstr(ci.mnc_fd, NC_GLOBAL, MIhistory, out_history);
544     free(out_history);
545 
546     in_header = 1;
547 
548     ci.frame_nbytes = 1;
549     ci.frame_nvoxels = 1;
550 
551     /* When we read voxels, we need COMBINED_SCALE_FACTOR() to have a sane
552      * value for all modalities. Set defaults for these in case the modality
553      * does not define one of these factors. For example, a CT (modality 2)
554      * will not define isotope_branching_fraction or calibration_factor.
555      */
556 
557     ci.scale_factor = 1.0;
558     ci.calibration_factor = 1.0;
559     ci.isotope_branching_fraction = 1.0;
560 
561     /* Collect the headers */
562     while (fgets(line_buf, sizeof(line_buf), ci.hdr_fp) != NULL) {
563         if (line_buf[0] == '#') /*  */
564             continue;
565         line_ptr = line_buf;
566         while (!isspace(*line_ptr)) {
567             line_ptr++;
568         }
569         *line_ptr++ = '\0';
570         val_ptr = line_ptr;
571         while (*line_ptr != '\n' && *line_ptr != '\r' && *line_ptr != '\0') {
572             line_ptr++;
573         }
574         *line_ptr = '\0';
575 
576         is_known = 0;
577 
578         if (in_header) {
579             if (*val_ptr != '\0') {
580                 /* Save the raw attribute into the file */
581                 ncattput(ci.mnc_fd, ncvarid(ci.mnc_fd, "micropet"),
582                          line_buf, NC_CHAR, strlen(val_ptr), val_ptr);
583             }
584 
585             for (ke_ptr = vol_atts; ke_ptr->upet_kwd != NULL; ke_ptr++) {
586                 if (!strcmp(ke_ptr->upet_kwd, line_buf)) {
587 
588                     is_known = 1;
589 
590                     if (ke_ptr->func != NULL) {
591                         (*ke_ptr->func)(&ci, val_ptr,
592                                         ke_ptr->mnc_var,
593                                         ke_ptr->mnc_att);
594                     }
595                     else if (ke_ptr->mnc_var != NULL &&
596                              ke_ptr->mnc_att != NULL) {
597 
598                         /* Interpret based upon type */
599                         switch (ke_ptr->upet_type) {
600                         case UPET_TYPE_INT:
601                             int_tmp = atoi(val_ptr);
602                             miattputint(ci.mnc_fd,
603                                         ncvarid(ci.mnc_fd, ke_ptr->mnc_var),
604                                         ke_ptr->mnc_att,
605                                         int_tmp);
606                             break;
607 
608                         case UPET_TYPE_REAL:
609                             dbl_tmp = atof(val_ptr);
610                             miattputdbl(ci.mnc_fd,
611                                         ncvarid(ci.mnc_fd, ke_ptr->mnc_var),
612                                         ke_ptr->mnc_att,
613                                         dbl_tmp);
614                             break;
615 
616                         case UPET_TYPE_STR:
617                             miattputstr(ci.mnc_fd,
618                                         ncvarid(ci.mnc_fd, ke_ptr->mnc_var),
619                                         ke_ptr->mnc_att,
620                                         val_ptr);
621                             break;
622 
623                         }
624 
625                     }
626                     break;
627                 }
628             }
629         }
630         else {
631             /* Not in the header any longer
632              */
633             for (ke_ptr = frm_atts; ke_ptr->upet_kwd != NULL; ke_ptr++) {
634                 if (!strcmp(ke_ptr->upet_kwd, line_buf)) {
635 
636                     is_known = 1;
637 
638                     if (ke_ptr->func != NULL) {
639                         (*ke_ptr->func)(&ci, val_ptr,
640                                         ke_ptr->mnc_var,
641                                         ke_ptr->mnc_att);
642                     }
643                     break;
644                 }
645             }
646         }
647 
648         if (!is_known) {
649             if (!strcmp(line_buf, "end_of_header")) {
650                 if (in_header) {
651                     in_header = 0;
652 
653                     copy_init(&ci);
654 
655                 }
656                 else {
657                     copy_frame(&ci);
658                 }
659             }
660             else {
661                 message(MSG_WARNING, "Unrecognized keyword %s\n", line_buf);
662             }
663         }
664     }
665 
666     fclose(ci.hdr_fp);
667     fclose(ci.img_fp);
668     miclose(ci.mnc_fd);
669     return (0);
670 }
671 
672 int
main(int argc,char ** argv)673 main(int argc, char **argv)
674 {
675     char *line_ptr;
676     int i;
677     char img_fname[1024];
678     char hdr_fname[1024];
679     char out_fname[1024];
680     int result;
681 
682     if (ParseArgv(&argc, argv, argTable, 0) || argc < 2) {
683         usage(argv[0]);
684         return (-1);
685     }
686 
687     ncopts = 0;
688 
689     /* Set the dimension names.  This is done here since the correct
690      * arrangement depends on the value of _orient_flag
691      */
692     _dimnames[DIM_T] = MItime;
693     _dimnames[DIM_X] = MIxspace;
694     _dimnames[DIM_Y] = MIyspace;
695     _dimnames[DIM_Z] = MIzspace;
696     _dimnames[DIM_W] = MIvector_dimension;
697 
698     if (_orient_flag == ORIENT_HEAD) {
699         /* If using head orientation, exchange Y and Z.
700          */
701         _dimnames[DIM_Y] = MIzspace;
702         _dimnames[DIM_Z] = MIyspace;
703     }
704 
705 
706     /* Open the header and the associated binary file. */
707 
708     for (i = 1; i < argc; i++) {
709         /* Here we try to be flexible about allowing the user to specify
710          * either the name of the .hdr file or the name of the .img file,
711          * or just the base name of the two files.  All three options
712          * should work.
713          */
714         strcpy(img_fname, argv[i]);
715         strcpy(hdr_fname, argv[i]);
716 
717         /* Find the last extension.
718          */
719         line_ptr = strrchr(argv[i], '.');
720 
721         /* Did the user specify the .hdr file??
722          */
723         if (line_ptr != NULL && !strcmp(line_ptr, ".hdr")) {
724             line_ptr = strrchr(img_fname, '.');
725             if (line_ptr != NULL) {
726                 *line_ptr = '\0';
727             }
728         }
729         /* Did the user specify the .img file??
730          */
731         else if (line_ptr != NULL && !strcmp(line_ptr, ".img")) {
732             strcat(hdr_fname, ".hdr");
733         }
734         /* Or perhaps just the base name??
735          */
736         else {
737             strcat(img_fname, ".img");
738             strcat(hdr_fname, ".img.hdr");
739         }
740 
741         /* See if there is a filename following this one, and if so, does it
742          * end with the ".mnc" extension.  If so, take that names as the
743          * output for this conversions.
744          */
745         if (i < argc - 1 &&
746             (line_ptr = strrchr(argv[i+1], '.')) != NULL &&
747             !strcmp(line_ptr, ".mnc")) {
748 
749             strcpy(out_fname, argv[i+1]);
750             i++;
751         }
752         else {
753             strcpy(out_fname, img_fname);
754             line_ptr = strrchr(out_fname, '.');
755             if (line_ptr != NULL) {
756                 strcpy(line_ptr, ".mnc");
757             }
758         }
759 
760         /* Perform the conversion.
761          */
762 
763         message(MSG_INFO, "Starting conversion\n");
764         message(MSG_INFO, "- Input header: %s\n", hdr_fname);
765         message(MSG_INFO, "- Input image:  %s\n", img_fname);
766         message(MSG_INFO, "- Output file:  %s\n", out_fname);
767 
768         result = upet_to_minc(hdr_fname, img_fname, out_fname, argv[0]);
769         if (result < 0) {
770             message(MSG_ERROR, "Error creating %s\n", out_fname);
771         }
772         else {
773             message(MSG_INFO, "Finished creating %s\n", out_fname);
774         }
775     }
776 
777     return 0;
778 }
779 
780 
DECLARE_FUNC(upet_file_type)781 DECLARE_FUNC(upet_file_type)
782 {
783     int file_type = atoi(val_str);
784     switch (file_type) {
785     case UPET_FT_IMAGE:         /* Image file */
786     case UPET_FT_MU_MAP:        /* Mu map file */
787         return (0);
788     default:
789         message(MSG_WARNING,
790                 "File type %d is not handled.  Conversion results may be problematic...\n", file_type);
791         break;
792     }
793     return (1);
794 }
795 
DECLARE_FUNC(upet_acq_mode)796 DECLARE_FUNC(upet_acq_mode)
797 {
798     int mode_int = atoi(val_str);
799     char *mode_str;
800 
801     switch (mode_int) {
802     case UPET_AM_UNKNOWN:
803         mode_str = "unknown";
804         break;
805     case UPET_AM_BLANK:
806         mode_str = "blank";
807         break;
808     case UPET_AM_EMISSION:
809         mode_str = "emission";
810         break;
811     case UPET_AM_DYNAMIC:
812         mode_str = "dynamic";
813         break;
814     case UPET_AM_GATED:
815         mode_str = "gated";
816         break;
817     case UPET_AM_CONTINUOUS_BED_MOTION:
818         mode_str = "continuous_bed_motion";
819         break;
820     case UPET_AM_SINGLES_TRANSMISSION:
821         mode_str = "singles_transmission";
822         break;
823     case UPET_AM_WINDOWED_COINCIDENCE_TRANSMISSION:
824         mode_str = "windowed_coincidence_transmission";
825         break;
826     case UPET_AM_NONWINDOWED_COINCIDENCE_TRANSMISSION:
827         mode_str = "non-windowed_coincidence_transmission";
828         break;
829     default:
830         message(MSG_WARNING, "Unknown acquisition mode %d\n", mode_int);
831         mode_str = NULL;
832         break;
833     }
834     if (mode_str != NULL) {
835         miattputstr(ci_ptr->mnc_fd, ncvarid(ci_ptr->mnc_fd, MIacquisition),
836                     "micropet_mode", mode_str);
837     }
838     return (0);
839 }
840 
DECLARE_FUNC(upet_bed_motion)841 DECLARE_FUNC(upet_bed_motion)
842 {
843     return (0);
844 }
845 
DECLARE_FUNC(upet_data_type)846 DECLARE_FUNC(upet_data_type)
847 {
848     ci_ptr->data_type = atoi(val_str);
849     switch (ci_ptr->data_type) {
850     case UPET_DT_BYTE:
851         ci_ptr->minc_type = NC_BYTE;
852         break;
853     case UPET_DT_II16:
854         ci_ptr->minc_type = NC_SHORT;
855         ci_ptr->frame_nbytes *= 2;
856         if (is_host_big_endian()) {
857             ci_ptr->swap_size = 2;
858         }
859         else {
860             ci_ptr->swap_size = 0;
861         }
862         break;
863     case UPET_DT_II32:
864         ci_ptr->minc_type = NC_INT;
865         ci_ptr->frame_nbytes *= 4;
866         if (is_host_big_endian()) {
867             ci_ptr->swap_size = 4;
868         }
869         else {
870             ci_ptr->swap_size = 0;
871         }
872         break;
873     case UPET_DT_IF32:
874         ci_ptr->minc_type = NC_FLOAT;
875         ci_ptr->frame_nbytes *= 4;
876         if (is_host_big_endian()) {
877             ci_ptr->swap_size = 4;
878         }
879         else {
880             ci_ptr->swap_size = 0;
881         }
882         break;
883     case UPET_DT_MF32:
884         ci_ptr->minc_type = NC_FLOAT;
885         ci_ptr->frame_nbytes *= 4;
886         if (!is_host_big_endian()) {
887             ci_ptr->swap_size = 4;
888         }
889         else {
890             ci_ptr->swap_size = 0;
891         }
892         break;
893     case UPET_DT_MI16:
894         ci_ptr->minc_type = NC_SHORT;
895         ci_ptr->frame_nbytes *= 2;
896         if (!is_host_big_endian()) {
897             ci_ptr->swap_size = 2;
898         }
899         else {
900             ci_ptr->swap_size = 0;
901         }
902         break;
903     case UPET_DT_MI32:
904         ci_ptr->minc_type = NC_INT;
905         ci_ptr->frame_nbytes *= 4;
906         if (!is_host_big_endian()) {
907             ci_ptr->swap_size = 4;
908         }
909         else {
910             ci_ptr->swap_size = 0;
911         }
912         break;
913     default:
914         message(MSG_ERROR, "Unknown data type %d\n", ci_ptr->data_type);
915         return (1);
916     }
917     if (ci_ptr->swap_size != 0) {
918         message(MSG_INFO, "Swapping groups of %d bytes.\n", ci_ptr->swap_size);
919     }
920     else {
921         message(MSG_INFO, "No byte-swapping required.\n");
922     }
923     return (0);
924 }
925 
DECLARE_FUNC(upet_data_order)926 DECLARE_FUNC(upet_data_order)
927 {
928     if (atoi(val_str) != 1) {
929         message(MSG_WARNING, "Unknown data order.\n");
930     }
931     return (0);
932 }
933 
DECLARE_FUNC(upet_ndims)934 DECLARE_FUNC(upet_ndims)
935 {
936     ci_ptr->dim_count = atoi(val_str);
937     return (0);
938 }
939 
940 static void
create_dimension(struct conversion_info * ci_ptr,int index,int length)941 create_dimension(struct conversion_info *ci_ptr, int index, int length)
942 {
943     ci_ptr->dim_lengths[index] = length;
944     if (index == DIM_W && length <= 1) {
945         return;
946     }
947 
948     ci_ptr->dim_ids[index] = ncdimdef(ci_ptr->mnc_fd,
949                                       _dimnames[index],
950                                       length);
951     if (index == DIM_T) {
952         micreate_std_variable(ci_ptr->mnc_fd, _dimnames[index],
953                               NC_DOUBLE, 1, &ci_ptr->dim_ids[index]);
954         micreate_std_variable(ci_ptr->mnc_fd, MItime_width,
955                               NC_DOUBLE, 1, &ci_ptr->dim_ids[index]);
956     }
957     else if (index != DIM_W) {
958         micreate_std_variable(ci_ptr->mnc_fd, _dimnames[index],
959                               NC_DOUBLE, 0, NULL);
960     }
961 }
962 
DECLARE_FUNC(upet_total_frames)963 DECLARE_FUNC(upet_total_frames)
964 {
965     create_dimension(ci_ptr, DIM_T, atoi(val_str));
966     return (0);
967 }
968 
DECLARE_FUNC(upet_x_dim)969 DECLARE_FUNC(upet_x_dim)
970 {
971     int x = atoi(val_str);
972     ci_ptr->frame_nbytes *= x;
973     ci_ptr->frame_nvoxels *= x;
974     create_dimension(ci_ptr, DIM_X, x);
975     return (0);
976 }
977 
DECLARE_FUNC(upet_y_dim)978 DECLARE_FUNC(upet_y_dim)
979 {
980     int y = atoi(val_str);
981     ci_ptr->frame_nbytes *= y;
982     ci_ptr->frame_nvoxels *= y;
983     create_dimension(ci_ptr, DIM_Y, y);
984     return (0);
985 }
986 
DECLARE_FUNC(upet_z_dim)987 DECLARE_FUNC(upet_z_dim)
988 {
989     int z = atoi(val_str);
990     ci_ptr->frame_nbytes *= z;
991     ci_ptr->frame_nvoxels *= z;
992     create_dimension(ci_ptr, DIM_Z, z);
993     return (0);
994 }
995 
DECLARE_FUNC(upet_vector_dim)996 DECLARE_FUNC(upet_vector_dim)
997 {
998     int w = atoi(val_str);
999     ci_ptr->frame_nbytes *= w;
1000     ci_ptr->frame_nvoxels *= w;
1001     create_dimension(ci_ptr, DIM_W, w);
1002     return (0);
1003 }
1004 
1005 /* Parse a micropet time string of the form: Ddd Mmm NN HH:MM:SS YYYY
1006  * e.g. Fri Jan 7 14:16:31 2005
1007  */
1008 static int
parse_time(char * str_ptr,struct tm * tm_ptr)1009 parse_time(char *str_ptr, struct tm *tm_ptr)
1010 {
1011     /* Just skip the first three characters. */
1012     while (*str_ptr != '\0' && *str_ptr != ' ') {
1013         str_ptr++;
1014     }
1015 
1016     while (*str_ptr == ' ') {
1017         str_ptr++;
1018     }
1019 
1020     /* Decode the month */
1021     if (str_ptr[0] == 'A') {
1022         if (str_ptr[1] == 'p') {
1023             tm_ptr->tm_mon = 4 - 1; /* April */
1024         }
1025         else {
1026             tm_ptr->tm_mon = 8 - 1; /* August */
1027         }
1028     }
1029     else if (str_ptr[0] == 'D') {
1030         tm_ptr->tm_mon = 12 - 1; /* December */
1031     }
1032     else if (str_ptr[0] == 'F') { /* February */
1033         tm_ptr->tm_mon = 2 - 1;
1034     }
1035     else if (str_ptr[0] == 'J') {
1036         if (str_ptr[1] == 'a') {
1037             tm_ptr->tm_mon = 1 - 1; /* January */
1038         }
1039         else if (str_ptr[2] == 'l') {
1040             tm_ptr->tm_mon = 7 - 1; /* July */
1041         }
1042         else {
1043             tm_ptr->tm_mon = 6 - 1; /* June */
1044         }
1045     }
1046     else if (str_ptr[0] == 'M') {
1047         if (str_ptr[2] == 'r') {
1048             tm_ptr->tm_mon = 3 - 1; /* March */
1049         }
1050         else {
1051             tm_ptr->tm_mon = 5 - 1; /* May */
1052         }
1053     }
1054     else if (str_ptr[0] == 'N') {
1055         tm_ptr->tm_mon = 11 - 1; /* November */
1056     }
1057     else if (str_ptr[0] == 'O') {
1058         tm_ptr->tm_mon = 10 - 1; /* October */
1059     }
1060     else if (str_ptr[0] == 'S') {
1061         tm_ptr->tm_mon = 9 - 1; /* September */
1062     }
1063     else {
1064         return 0;
1065     }
1066 
1067     /* Skip past the month */
1068     while (*str_ptr != ' ' && *str_ptr != '\0') {
1069         str_ptr++;
1070     }
1071 
1072     while (*str_ptr == ' ') {
1073         str_ptr++;
1074     }
1075 
1076     tm_ptr->tm_mday = 0;
1077     while (isdigit(*str_ptr)) {
1078         tm_ptr->tm_mday = (tm_ptr->tm_mday * 10) + (*str_ptr++ - '0');
1079     }
1080 
1081     while (*str_ptr == ' ') {
1082         str_ptr++;
1083     }
1084 
1085     tm_ptr->tm_hour = 0;
1086     while (isdigit(*str_ptr)) {
1087         tm_ptr->tm_hour = (tm_ptr->tm_hour * 10) + (*str_ptr++ - '0');
1088     }
1089 
1090     if (*str_ptr == ':') {
1091         str_ptr++;
1092     }
1093     else {
1094         return 0;
1095     }
1096 
1097     tm_ptr->tm_min = 0;
1098     while (isdigit(*str_ptr)) {
1099         tm_ptr->tm_min = (tm_ptr->tm_min * 10) + (*str_ptr++ - '0');
1100     }
1101 
1102     if (*str_ptr == ':') {
1103         str_ptr++;
1104     }
1105     else {
1106         return 0;
1107     }
1108 
1109     tm_ptr->tm_sec = 0;
1110     while (isdigit(*str_ptr)) {
1111         tm_ptr->tm_sec = (tm_ptr->tm_sec * 10) + (*str_ptr++ - '0');
1112     }
1113 
1114     while (*str_ptr == ' ') {
1115         str_ptr++;
1116     }
1117 
1118     tm_ptr->tm_year = 0;
1119     while (isdigit(*str_ptr)) {
1120         tm_ptr->tm_year = (tm_ptr->tm_year * 10) + (*str_ptr++ - '0');
1121     }
1122     tm_ptr->tm_year -= 1900;
1123 
1124     return 1;
1125 }
1126 
DECLARE_FUNC(upet_injection_time)1127 DECLARE_FUNC(upet_injection_time)
1128 {
1129     struct tm tmbuf;
1130     int id;
1131     char str_buf[128];
1132 
1133     id = ncvarid(ci_ptr->mnc_fd, MIacquisition);
1134 
1135     if (!parse_time(val_str, &tmbuf)) {
1136         strcpy(str_buf, "unknown");
1137         miattputstr(ci_ptr->mnc_fd, id, MIinjection_time, str_buf);
1138         miattputstr(ci_ptr->mnc_fd, id, MIinjection_year, str_buf);
1139         miattputstr(ci_ptr->mnc_fd, id, MIinjection_month, str_buf);
1140         miattputstr(ci_ptr->mnc_fd, id, MIinjection_day, str_buf);
1141     }
1142     else {
1143         sprintf(str_buf, "%02d%02d%02d",
1144                 tmbuf.tm_hour, tmbuf.tm_min, tmbuf.tm_sec);
1145         miattputstr(ci_ptr->mnc_fd, id, MIinjection_time, str_buf);
1146 
1147         sprintf(str_buf, "%d", tmbuf.tm_year + 1900);
1148         miattputstr(ci_ptr->mnc_fd, id, MIinjection_year, str_buf);
1149 
1150         sprintf(str_buf, "%d", tmbuf.tm_mon + 1);
1151         miattputstr(ci_ptr->mnc_fd, id, MIinjection_month, str_buf);
1152 
1153         sprintf(str_buf, "%d", tmbuf.tm_mday);
1154         miattputstr(ci_ptr->mnc_fd, id, MIinjection_day, str_buf);
1155     }
1156     return (0);
1157 }
1158 
DECLARE_FUNC(upet_scan_time)1159 DECLARE_FUNC(upet_scan_time)
1160 {
1161     struct tm tmbuf;
1162     int id;
1163     char str_buf[128];
1164 
1165     id = ncvarid(ci_ptr->mnc_fd, MIstudy);
1166 
1167     if (!parse_time(val_str, &tmbuf)) {
1168         strcpy(str_buf, "unknown");
1169         miattputstr(ci_ptr->mnc_fd, id, MIstart_time, str_buf);
1170         miattputstr(ci_ptr->mnc_fd, id, MIstart_year, str_buf);
1171         miattputstr(ci_ptr->mnc_fd, id, MIstart_month, str_buf);
1172         miattputstr(ci_ptr->mnc_fd, id, MIstart_day, str_buf);
1173     }
1174     else {
1175         sprintf(str_buf, "%02d%02d%02d",
1176                 tmbuf.tm_hour, tmbuf.tm_min, tmbuf.tm_sec);
1177         miattputstr(ci_ptr->mnc_fd, id, MIstart_time, str_buf);
1178 
1179         sprintf(str_buf, "%d", tmbuf.tm_year + 1900);
1180         miattputstr(ci_ptr->mnc_fd, id, MIstart_year, str_buf);
1181 
1182         sprintf(str_buf, "%d", tmbuf.tm_mon + 1);
1183         miattputstr(ci_ptr->mnc_fd, id, MIstart_month, str_buf);
1184 
1185         sprintf(str_buf, "%d", tmbuf.tm_mday);
1186         miattputstr(ci_ptr->mnc_fd, id, MIstart_day, str_buf);
1187     }
1188     return (0);
1189 }
1190 
DECLARE_FUNC(upet_axial_crystal_pitch)1191 DECLARE_FUNC(upet_axial_crystal_pitch)
1192 {
1193     double dbl_tmp = atof(val_str);
1194 
1195     /* dbl_tmp is in cm.  Convert to mm. */
1196     dbl_tmp *= 10.0;
1197 
1198     /* Now convert from crystal pitch to actual slice width */
1199     dbl_tmp /= 2.0;
1200 
1201     ci_ptr->dim_steps[DIM_Z] = dbl_tmp;
1202     return (0);
1203 }
1204 
DECLARE_FUNC(upet_pixel_size)1205 DECLARE_FUNC(upet_pixel_size)
1206 {
1207     double dbl_tmp = atof(val_str);
1208 
1209     /* dbl_tmp is in cm.  Convert to mm. */
1210     dbl_tmp *= 10.0;
1211 
1212     ci_ptr->dim_steps[DIM_X] = dbl_tmp;
1213     ci_ptr->dim_steps[DIM_Y] = dbl_tmp;
1214     return (0);
1215 }
1216 
DECLARE_FUNC(upet_dose_units)1217 DECLARE_FUNC(upet_dose_units)
1218 {
1219     int tmp = atoi(val_str);
1220     char *str_ptr;
1221 
1222     if (tmp == 0) {
1223         str_ptr = "unknown";
1224     }
1225     else if (tmp == 1) {
1226         str_ptr = "mCi";
1227     }
1228     else if (tmp == 2) {
1229         str_ptr = "MBq";
1230     }
1231     else {
1232         str_ptr = "???????";
1233         message(MSG_WARNING, "Unrecognized dose_units value %d\n", tmp);
1234     }
1235     miattputstr(ci_ptr->mnc_fd, ncvarid(ci_ptr->mnc_fd, new_var), new_att, str_ptr);
1236     return (0);
1237 }
1238 
DECLARE_FUNC(upet_calibration_factor)1239 DECLARE_FUNC(upet_calibration_factor)
1240 {
1241     double dbl_tmp = atof(val_str);
1242 
1243     ci_ptr->calibration_factor = dbl_tmp;
1244     return (0);
1245 }
1246 
DECLARE_FUNC(upet_isotope_branching_fraction)1247 DECLARE_FUNC(upet_isotope_branching_fraction)
1248 {
1249     double dbl_tmp = atof(val_str);
1250 
1251     ci_ptr->isotope_branching_fraction = dbl_tmp;
1252     return (0);
1253 }
1254 
1255 
DECLARE_FUNC(upet_rotation)1256 DECLARE_FUNC(upet_rotation)
1257 {
1258     double dbl_tmp = atof(val_str);
1259     if (dbl_tmp != 0.0) {
1260         message(MSG_WARNING, "Rotation is %f\n", dbl_tmp);
1261     }
1262     return (0);
1263 }
1264 
1265 /***********************/
1266 /* Per-frame functions */
DECLARE_FUNC(upet_frame_no)1267 DECLARE_FUNC(upet_frame_no)
1268 {
1269     ci_ptr->frame_index = atoi(val_str);
1270     /* Set index of "zeroth" frame if not already set.
1271      */
1272     if (ci_ptr->frame_zero < 0) {
1273         ci_ptr->frame_zero = ci_ptr->frame_index;
1274     }
1275     return (0);
1276 }
1277 
DECLARE_FUNC(upet_frame_start)1278 DECLARE_FUNC(upet_frame_start)
1279 {
1280     long index = ci_ptr->frame_index - ci_ptr->frame_zero;
1281     double dbl_tmp = atof(val_str);
1282     mivarput1(ci_ptr->mnc_fd, ncvarid(ci_ptr->mnc_fd, MItime), &index,
1283               NC_DOUBLE, MI_SIGNED, &dbl_tmp);
1284     return (0);
1285 }
1286 
DECLARE_FUNC(upet_frame_duration)1287 DECLARE_FUNC(upet_frame_duration)
1288 {
1289     long index = ci_ptr->frame_index - ci_ptr->frame_zero;
1290     double dbl_tmp = atof(val_str);
1291     mivarput1(ci_ptr->mnc_fd, ncvarid(ci_ptr->mnc_fd, MItime_width), &index,
1292               NC_DOUBLE, MI_SIGNED, &dbl_tmp);
1293     return (0);
1294 }
1295 
DECLARE_FUNC(upet_frame_min)1296 DECLARE_FUNC(upet_frame_min)
1297 {
1298     long index = ci_ptr->frame_index - ci_ptr->frame_zero;
1299     double dbl_tmp = atof(val_str);
1300     dbl_tmp *= COMBINED_SCALE_FACTOR(ci_ptr);
1301     mivarput1(ci_ptr->mnc_fd, ncvarid(ci_ptr->mnc_fd, MIimagemin), &index,
1302               NC_DOUBLE, MI_SIGNED, &dbl_tmp);
1303     return (0);
1304 }
1305 
DECLARE_FUNC(upet_frame_max)1306 DECLARE_FUNC(upet_frame_max)
1307 {
1308     long index = ci_ptr->frame_index - ci_ptr->frame_zero;
1309     double dbl_tmp = atof(val_str);
1310     dbl_tmp *= COMBINED_SCALE_FACTOR(ci_ptr);
1311     mivarput1(ci_ptr->mnc_fd, ncvarid(ci_ptr->mnc_fd, MIimagemax), &index,
1312               NC_DOUBLE, MI_SIGNED, &dbl_tmp);
1313     return (0);
1314 }
1315 
DECLARE_FUNC(upet_frame_file_ptr)1316 DECLARE_FUNC(upet_frame_file_ptr)
1317 {
1318     long hipart;
1319     long lopart;
1320     char *end_ptr;
1321 
1322     lopart = strtol(val_str, &end_ptr, 10);
1323     if (*end_ptr == ' ') {
1324         while (*end_ptr == ' ') {
1325             end_ptr++;
1326         }
1327         if (isdigit(*end_ptr)) {
1328             hipart = lopart;
1329             lopart = strtol(end_ptr, NULL, 10);
1330         }
1331     }
1332 
1333     /* Seek the image file to the data */
1334 
1335     fseek(ci_ptr->img_fp, lopart, SEEK_SET);
1336     return (0);
1337 }
1338 
DECLARE_FUNC(upet_frame_scale_factor)1339 DECLARE_FUNC(upet_frame_scale_factor)
1340 {
1341     ci_ptr->scale_factor = atof(val_str);
1342     return (0);
1343 }
1344 
1345 
DECLARE_FUNC(upet_frame_deadtime_correction)1346 DECLARE_FUNC(upet_frame_deadtime_correction)
1347 {
1348     ci_ptr->deadtime_correction = atof(val_str);
1349     return (0);
1350 }
1351 
DECLARE_FUNC(upet_frame_decay_correction)1352 DECLARE_FUNC(upet_frame_decay_correction)
1353 {
1354     ci_ptr->decay_correction = atof(val_str);
1355     return (0);
1356 }
1357 
1358 static void
swap_data(int swap_size,long nvox,unsigned char * data)1359 swap_data(int swap_size,
1360           long nvox,
1361           unsigned char *data)
1362 {
1363   unsigned char tmp;
1364 
1365   if (swap_size == 2) {
1366     while (nvox--) {
1367       tmp = data[0];
1368       data[0] = data[1];
1369       data[1] = tmp;
1370       data += 2;
1371     }
1372   }
1373   else if (swap_size == 4) {
1374     while (nvox--) {
1375       tmp = data[0];
1376       data[0] = data[3];
1377       data[3] = tmp;
1378       tmp = data[1];
1379       data[1] = data[2];
1380       data[2] = tmp;
1381       data += 4;
1382     }
1383   }
1384 }
1385 
1386 static void
scale_data(nc_type datatype,long nvox,void * data,double scale)1387 scale_data(nc_type datatype,
1388            long nvox,
1389            void *data,
1390            double scale)
1391 {
1392     long i;
1393     double tmp;
1394 
1395     switch (datatype) {
1396     case NC_BYTE:
1397         for (i = 0; i < nvox; i++) {
1398             tmp = (double) ((char *)data)[i];
1399             tmp *= scale;
1400             ((char *)data)[i] = tmp;
1401         }
1402         break;
1403     case NC_SHORT:
1404         for (i = 0; i < nvox; i++) {
1405             tmp = (double) ((short *)data)[i];
1406             tmp *= scale;
1407             ((short *)data)[i] = tmp;
1408         }
1409         break;
1410     case NC_INT:
1411         for (i = 0; i < nvox; i++) {
1412             tmp = (double) ((int *)data)[i];
1413             tmp *= scale;
1414             ((int *)data)[i] = tmp;
1415         }
1416         break;
1417     case NC_FLOAT:
1418         for (i = 0; i < nvox; i++) {
1419             tmp = (double) ((float *)data)[i];
1420             tmp *= scale;
1421             ((float *)data)[i] = tmp;
1422         }
1423         break;
1424     case NC_DOUBLE:
1425         for (i = 0; i < nvox; i++) {
1426             tmp = (double) ((double *)data)[i];
1427             tmp *= scale;
1428             ((double *)data)[i] = tmp;
1429         }
1430         break;
1431     default:
1432         message(MSG_ERROR, "Data type %d not handled\n", datatype);
1433         break;
1434     }
1435 }
1436 
1437 void
copy_init(struct conversion_info * ci_ptr)1438 copy_init(struct conversion_info *ci_ptr)
1439 {
1440     ci_ptr->frame_buffer = malloc(ci_ptr->frame_nbytes);
1441     if (ci_ptr->frame_buffer == NULL) {
1442         message(MSG_FATAL, "Out of memory\n");
1443         exit(-1);
1444     }
1445 
1446     /* Create the image, imagemax, and imagemin variables.
1447      */
1448     micreate_std_variable(ci_ptr->mnc_fd, MIimagemax, NC_DOUBLE,
1449                           1, ci_ptr->dim_ids);
1450     micreate_std_variable(ci_ptr->mnc_fd, MIimagemin, NC_DOUBLE,
1451                           1, ci_ptr->dim_ids);
1452     micreate_std_variable(ci_ptr->mnc_fd, MIimage, ci_ptr->minc_type,
1453                           ci_ptr->dim_count + 1, ci_ptr->dim_ids);
1454 
1455     /* Set up the dimension step and start values.  Because of the microPET
1456      * data orientation, we set Z and Y to be the inverse of the norm, to
1457      * put the animal's nose at the top of the display.
1458      * TODO: allow this behavior to be controlled on the command line.
1459      */
1460     miattputdbl(ci_ptr->mnc_fd, ncvarid(ci_ptr->mnc_fd, _dimnames[DIM_Z]),
1461                 MIstep, -ci_ptr->dim_steps[DIM_Z]);
1462     miattputdbl(ci_ptr->mnc_fd, ncvarid(ci_ptr->mnc_fd, _dimnames[DIM_X]),
1463                 MIstep, ci_ptr->dim_steps[DIM_X]);
1464     miattputdbl(ci_ptr->mnc_fd, ncvarid(ci_ptr->mnc_fd, _dimnames[DIM_Y]),
1465                 MIstep, -ci_ptr->dim_steps[DIM_Y]);
1466 
1467     miattputdbl(ci_ptr->mnc_fd, ncvarid(ci_ptr->mnc_fd, _dimnames[DIM_Z]),
1468                 MIstart, ci_ptr->dim_steps[DIM_Z] * ci_ptr->dim_lengths[DIM_Z]);
1469     miattputdbl(ci_ptr->mnc_fd, ncvarid(ci_ptr->mnc_fd, _dimnames[DIM_X]),
1470                 MIstart, 0.0);
1471     miattputdbl(ci_ptr->mnc_fd, ncvarid(ci_ptr->mnc_fd, _dimnames[DIM_Y]),
1472                 MIstart, ci_ptr->dim_steps[DIM_Y] * ci_ptr->dim_lengths[DIM_Y]);
1473 
1474     miattputstr(ci_ptr->mnc_fd, ncvarid(ci_ptr->mnc_fd, _dimnames[DIM_Z]),
1475                 MIunits, "mm");
1476     miattputstr(ci_ptr->mnc_fd, ncvarid(ci_ptr->mnc_fd, _dimnames[DIM_X]),
1477                 MIunits, "mm");
1478     miattputstr(ci_ptr->mnc_fd, ncvarid(ci_ptr->mnc_fd, _dimnames[DIM_Y]),
1479                 MIunits, "mm");
1480     miattputstr(ci_ptr->mnc_fd, ncvarid(ci_ptr->mnc_fd, _dimnames[DIM_T]),
1481                 MIunits, "s");
1482 
1483     ncendef(ci_ptr->mnc_fd);
1484 }
1485 
1486 static void
copy_frame(struct conversion_info * ci_ptr)1487 copy_frame(struct conversion_info *ci_ptr)
1488 {
1489     long start[5];
1490     long count[5];
1491     size_t nitems;
1492 
1493     message(MSG_INFO, "Inserting frame #%d\n", ci_ptr->frame_index);
1494 
1495     /* Actually read the data from the image file.
1496      */
1497     nitems = fread(ci_ptr->frame_buffer, ci_ptr->frame_nbytes, 1,
1498                    ci_ptr->img_fp);
1499 
1500     if (nitems != 1) {
1501         message(MSG_FATAL, "Read failed with error %d, return %d\n",
1502                 errno, nitems);
1503         exit(-1);
1504     }
1505 
1506     /* Setup the starts and counts for the data block.
1507      */
1508     start[DIM_T] = ci_ptr->frame_index - ci_ptr->frame_zero;
1509     start[DIM_X] = 0;
1510     start[DIM_Y] = 0;
1511     start[DIM_Z] = 0;
1512     start[DIM_W] = 0;
1513 
1514     count[DIM_T] = 1;
1515     count[DIM_X] = ci_ptr->dim_lengths[DIM_X];
1516     count[DIM_Y] = ci_ptr->dim_lengths[DIM_Y];
1517     count[DIM_Z] = ci_ptr->dim_lengths[DIM_Z];
1518     count[DIM_W] = ci_ptr->dim_lengths[DIM_W];
1519 
1520     /* Perform swapping if necessary.
1521      */
1522     if (ci_ptr->swap_size != 0) {
1523       swap_data(ci_ptr->swap_size, ci_ptr->frame_nvoxels,
1524                 ci_ptr->frame_buffer);
1525     }
1526 
1527     /* Scale the raw data into the final range.
1528      */
1529     scale_data(ci_ptr->minc_type, ci_ptr->frame_nvoxels, ci_ptr->frame_buffer,
1530                COMBINED_SCALE_FACTOR(ci_ptr));
1531 
1532     /* For now we perform no conversions on the data as it is stored.
1533      * This may be worth modifying in the future, to allow storage of
1534      * non-floating-point formats from a typical microPET file.
1535      */
1536     ncvarput(ci_ptr->mnc_fd, ncvarid(ci_ptr->mnc_fd, MIimage), start, count,
1537              ci_ptr->frame_buffer);
1538 }
1539 
1540