1 /* ----------------------------- MNI Header -----------------------------------
2 @NAME : dicom_to_minc.c
3 @DESCRIPTION: Code to convert a list of DICOM files to minc
4 format.
5 @METHOD :
6 @GLOBALS :
7 @CALLS :
8 @CREATED : January 28, 1997 (Peter Neelin)
9 @MODIFIED :
10 * $Log: dicom_to_minc.c,v $
11 * Revision 1.30 2010-11-23 23:30:50 claude
12 * dcm2mnc: fixed seg fault bug (Claude) and added b-matrix (Ilana)
13 *
14 * Revision 1.29 2010-06-23 13:21:05 rotor
15 * * changed H5Acreate2 calls back to the H5Acreate macro
16 *
17 * Revision 1.28 2008/08/12 05:00:23 rotor
18 * * large number of changes from Claude (64 bit and updates)
19 *
20 * Revision 1.27 2008/01/17 06:20:54 stever
21 * * conversion/dcm2mnc/dicom_to_minc.c (copy_element_properties):
22 * Change return type from int to void; no callers require a return value.
23 *
24 * * conversion/micropet/upet2mnc.c (main): Return 0 at end of function.
25 *
26 * Revision 1.26 2008/01/17 02:33:01 rotor
27 * * removed all rcsids
28 * * removed a bunch of ^L's that somehow crept in
29 * * removed old (and outdated) BUGS file
30 *
31 * Revision 1.25 2008/01/12 19:08:14 stever
32 * Add __attribute__ ((unused)) to all rcsid variables.
33 *
34 * Revision 1.24 2008/01/11 07:17:07 stever
35 * Remove unused variables.
36 *
37 * Revision 1.23 2007/11/23 20:28:23 ilana
38 * condition on picking between DICOM slice spacing and coordinate spacing was wrong (extra ! in if statement)
39 *
40 * Revision 1.22 2007/08/13 16:34:52 ilana
41 * mods to handle naming scheme of more diffusion sequences
42 *
43 * Revision 1.21 2007/06/08 20:28:57 ilana
44 * added several fields to mincheader (dicom elements and found in ASCONV header)
45 *
46 * Revision 1.20 2007/05/30 15:17:34 ilana
47 * fix so that diffusion images all written into 1 4d volume, gradient directions and bvalues are written to mincheader, some fixes for TIM diffusion images
48 *
49 * Revision 1.19 2006/05/11 14:45:14 bert
50 * Fix endian-ness issues when parsing Siemens proprietary fields
51 *
52 * Revision 1.18 2006/04/09 15:39:04 bert
53 * Improve support for DTI
54 *
55 * Revision 1.17 2006/02/09 20:54:29 bert
56 * More changes to dcm2mnc
57 *
58 * Revision 1.16 2005/12/13 17:31:13 bert
59 * Ignore DICOM protocol errors. This change was necessitated by images from a Philips Intera scanner version 'NT 10.4.1\\PIIM V2.1.4.1 MIMIT MCS' that appears to set the DICOM length field incorrectly.
60 *
61 * Revision 1.15 2005/08/26 21:25:54 bert
62 * Latest changes ported from 1.0 branch
63 *
64 * Revision 1.13.2.5 2005/08/18 18:18:35 bert
65 * Implement the -usecoordinates option, also some minor cleanup and fixes for some warning messages.
66 *
67 * Revision 1.13.2.4 2005/08/18 16:38:44 bert
68 * Minor updates for dealing w/older numaris data
69 *
70 * Revision 1.13.2.3 2005/07/14 16:47:55 bert
71 * Handle additional classes of Numa3 mosaics by using the Siemens 'base raw matrix size' element
72 *
73 * Revision 1.13.2.2 2005/06/20 22:01:15 bert
74 * Add basic support for multiframe DICOM images
75 *
76 * Revision 1.13.2.1 2005/05/12 21:16:47 bert
77 * Initial checkin
78 *
79 * Revision 1.13 2005/05/09 15:33:20 bert
80 * Fix handing of descending mosaic sequences; don't print GEMS field missing warnings by default
81 *
82 * Revision 1.12 2005/04/28 17:38:06 bert
83 * Insert and retrieve width information when sorting a dimension
84 *
85 * Revision 1.11 2005/04/21 22:31:29 bert
86 * Copy SPI_Magnetic_field_strength to standard field in copy_spi_to_acr(); Relax some tests to avoid spurious warnings
87 *
88 * Revision 1.10 2005/04/20 23:25:50 bert
89 * Add copy_spi_to_acr() function, minor name and comment changes
90 *
91 * Revision 1.9 2005/04/20 17:48:16 bert
92 * Copy SPI_Number_of_slices_nominal to ACR_Image_in_acquisition for Siemens
93 *
94 * Revision 1.8 2005/04/18 21:04:25 bert
95 * Get rid of old is_reversed behavior, always use most natural sort of the image. Also tried to fix listing function somewhat.
96 *
97 * Revision 1.7 2005/04/05 21:55:33 bert
98 * Update handling of Siemens ASCCONV to reflect value of uc2DInterpolation field for mosaic images. Additional cleanup of mosaic code and minor tweak to suppress GEMS warning message for PET data.
99 *
100 * Revision 1.6 2005/03/29 20:20:42 bert
101 * Add checks for GE Medical Systems proprietary files
102 *
103 * Revision 1.5 2005/03/14 22:26:40 bert
104 * Dump the entire coordinate array for each MRI dimension after sorting. Also detect GE scans.
105 *
106 * Revision 1.4 2005/03/13 19:35:11 bert
107 * Lots of changes for dealing with some proprietary Philips stuff
108 *
109 * Revision 1.3 2005/03/03 18:59:15 bert
110 * Fix handling of image position so that we work with the older field (0020, 0030) as well as the new (0020, 0032)
111 *
112 * Revision 1.2 2005/03/02 20:16:24 bert
113 * Latest changes and cleanup
114 *
115 * Revision 1.1 2005/02/17 16:38:10 bert
116 * Initial checkin, revised DICOM to MINC converter
117 *
118 * Revision 1.1.1.1 2003/08/15 19:52:55 leili
119 * Leili's dicom server for sonata
120 *
121 * Revision 1.18 2002/09/26 15:24:33 rhoge
122 * Before was only skipping time sort for multi-slice N4 mosaics. Turns out
123 * this was also causing failure on single-slice scans.
124 * Changed slices>1 to slices>0 so that now N4 dicom scans never get sorted
125 * on their (apparently nonsensical) time value. The if statement should
126 * really be reworked, and should keep an eye on this. Seems like EPI
127 * time series sequences never sort properly on 'time'.
128 *
129 * Revision 1.17 2002/09/25 17:25:43 rhoge
130 * changed public void's to public int's
131 *
132 * Revision 1.16 2002/05/08 19:32:40 rhoge
133 * fixed handling of diffusion scans with separate series for each average
134 *
135 * Revision 1.15 2002/05/01 21:29:34 rhoge
136 * removed MrProt from minc header - encountered files with large strings,
137 * causing seg faults
138 *
139 * Revision 1.14 2002/04/30 12:36:35 rhoge
140 * fixes to handle current (and hopefully final) diffusion sequence
141 *
142 * Revision 1.13 2002/04/26 03:27:03 rhoge
143 * fixed MrProt problem - replaced fixed lenght char array with malloc
144 *
145 * Revision 1.12 2002/04/08 03:40:56 rhoge
146 * fixed mosaic extraction for non-square scans and 3D scans.
147 * added some new dicom elements
148 *
149 * Revision 1.11 2002/03/27 19:38:08 rhoge
150 * small comment change
151 *
152 * Revision 1.10 2002/03/27 18:57:50 rhoge
153 * added diffusion b value
154 *
155 * Revision 1.9 2002/03/23 13:17:53 rhoge
156 * added support for Bourget network pushed dicom files, cleaned up
157 * file check and read_numa4_dicom vr check/assignment
158 *
159 * Revision 1.8 2002/03/22 19:19:36 rhoge
160 * Numerous fixes -
161 * - handle Numaris 4 Dicom patient name
162 * - option to cleanup input files
163 * - command option
164 * - list-only option
165 * - debug mode
166 * - user supplied name, idstr
167 * - anonymization
168 *
169 * Revision 1.7 2002/03/21 13:31:56 rhoge
170 * updated comments
171 *
172 * Revision 1.6 2002/03/19 22:10:16 rhoge
173 * removed time sorting for N4DCM mosaics - time is random for mosaics
174 *
175 * Revision 1.5 2002/03/19 13:13:56 rhoge
176 * initial working mosaic support - I think time is scrambled though.
177 *
178 * Revision 1.4 2001/12/31 18:27:21 rhoge
179 * modifications for dicomreader processing of Numaris 4 dicom files - at
180 * this point code compiles without warning, but does not deal with
181 * mosaiced files. Also will probably not work at this time for Numaris
182 * 3 .ima files. dicomserver may also not be functional...
183 *
184 * Revision 1.3 2000/12/14 21:37:11 rhoge
185 * log message cleanup
186 *
187 * Revision 1.2 2000/12/14 21:36:22 rhoge
188 * changes to restore measurement loop support that was broken by changes
189 * to provide acquisition loop support
190 *
191 * Revision 1.1.1.1 2000/11/30 02:13:15 rhoge
192 * imported sources to CVS repository on amoeba
193 * -now support Siemens acquisition loop scans with and without correction
194 * on sending side
195 *
196 * Revision 6.1 1999/10/29 17:51:59 neelin
197 * Fixed Log keyword
198 *
199 * Revision 6.0 1997/09/12 13:24:27 neelin
200 * Release of minc version 0.6
201 *
202 * Revision 5.0 1997/08/21 13:25:26 neelin
203 * Release of minc version 0.5
204 *
205 * Revision 4.0 1997/05/07 20:06:20 neelin
206 * Release of minc version 0.4
207 *
208 * Revision 1.1 1997/03/04 20:56:47 neelin
209 * Initial revision
210 *
211 @COPYRIGHT : Copyright 1997 Peter Neelin, McConnell Brain Imaging
212 Centre, Montreal Neurological Institute, McGill University.
213 Permission to use, copy, modify, and distribute this software and
214 its documentation for any purpose and without fee is hereby
215 granted, provided that the above copyright notice appear in all
216 copies. The author and McGill University make no representations
217 about the suitability of this software for any purpose. It is
218 provided "as is" without express or implied warranty.
219 ---------------------------------------------------------------------------- */
220
221 #include "dcm2mnc.h"
222 const char *World_Names[WORLD_NDIMS] = { "X", "Y", "Z" };
223 const char *Volume_Names[VOL_NDIMS] = { "Slice", "Row", "Column" };
224 const char *Mri_Names[MRI_NDIMS] = {"Slice", "Echo", "Time", "Phase", "ChmSh"};
225
226 /* Private structure definitions. */
227
228 /* multi-image (mosaic) info */
229 typedef struct {
230 int packed;
231 mosaic_seq_t mosaic_seq;
232 int size[2];
233 int big[2];
234 int grid[2];
235 int pixel_size;
236 Acr_Element big_image;
237 Acr_Element small_image;
238 int sub_images;
239 int slice_count;
240 double normal[WORLD_NDIMS];
241 double step[WORLD_NDIMS];
242 double position[WORLD_NDIMS];
243 } Mosaic_Info;
244
245 /* DICOM Multiframe information. NOTE: This represents at best a very
246 * partial implementation of the DICOM multiframe specification, barely
247 * adequate for conversion of basic 3D files. More work is needed to
248 * support dynamic scans and other more complex objects.
249 */
250 typedef struct {
251 int frame_count;
252 int frame_size;
253 Acr_Element big_image;
254 Acr_Element sub_image;
255 double normal[WORLD_NDIMS];
256 double step[WORLD_NDIMS];
257 double position[WORLD_NDIMS];
258 } Multiframe_Info;
259
260 /* For dicom_to_minc(), these codes specify whether we are dealing
261 * with a "normal" DICOM sequence, a Siemens mosaic sequence, or a
262 * DICOM multiframe sequence.
263 */
264 #define SUBIMAGE_TYPE_NONE 0
265 #define SUBIMAGE_TYPE_MOSAIC 1
266 #define SUBIMAGE_TYPE_MULTIFRAME 2
267
268 /* Structure for sorting dimensions */
269 typedef struct {
270 int identifier;
271 int original_index;
272 double value;
273 double width;
274 } Sort_Element;
275
276 /* Private function definitions */
277 static int mosaic_init(Acr_Group, Mosaic_Info *, int);
278 static void mosaic_cleanup(Mosaic_Info *);
279 static int mosaic_insert_subframe(Acr_Group, Mosaic_Info *, int, int);
280
281 /* DICOM Multiframe conversion functions (see NOTE: above) */
282 static void multiframe_init(Acr_Group, Multiframe_Info *, int);
283 static void multiframe_cleanup(Multiframe_Info *);
284 static void multiframe_insert_subframe(Acr_Group, Multiframe_Info *, int, int);
285
286 static void free_info(General_Info *gi_ptr, File_Info *fi_ptr,
287 int num_files);
288 static int dimension_sort_function(const void *v1, const void *v2);
289 static void sort_dimensions(General_Info *gi_ptr);
290 static int prot_find_string(Acr_Element Protocol, const char *name,
291 char *value);
292 static char *dump_protocol_text(Acr_Element Protocol);
293
294 /* ----------------------------- MNI Header -----------------------------------
295 @NAME : dicom_to_minc
296 @INPUT : num_files - number of image files
297 file_list - list of file names
298 minc_file - name of minc file to create, or NULL to make one up.
299 clobber - if TRUE, then open the output with NC_CLOBBER
300 file_prefix - string providing any directory or prefix
301 for internally generated filename (if it is a directory,
302 then it must contain the last "/")
303 @OUTPUT : output_file_name - returns a pointer to an internal area
304 containing the file name of the created file if minc_file
305 is NULL, or simply a pointer to minc_file. If NULL, then
306 nothing is returned.
307 @RETURNS : EXIT_SUCCESS if no error, EXIT_FAILURE on error.
308 @DESCRIPTION: Routine to convert a list of Siemens dicom files to minc
309 format.
310 @METHOD :
311 @GLOBALS :
312 @CALLS :
313 @CREATED : November 25, 1993 (Peter Neelin)
314 @MODIFIED :
315 ---------------------------------------------------------------------------- */
316
317 int
dicom_to_minc(int num_files,const char * file_list[],const char * minc_file,int clobber,const char * file_prefix,char ** output_file_name)318 dicom_to_minc(int num_files,
319 const char *file_list[],
320 const char *minc_file,
321 int clobber,
322 const char *file_prefix,
323 char **output_file_name)
324 {
325 Acr_Group group_list; /* List of ACR/NEMA groups & elements */
326 File_Info *fi_ptr; /* Array of per-file information */
327 General_Info gi; /* General (common) DICOM file information */
328 int max_group; /* Maximum group number to read */
329 Image_Data image; /* Actual image data */
330 int icvid; /* MINC Image Conversion Variable */
331 int ifile; /* File index */
332 Mri_Index imri; /* MRI axis index */
333 char *out_file_name; /* Output MINC filename */
334 int isep; /* Loop counter */
335 const Loop_Type loop_type = NONE; /* MINC loop type always none for now */
336 int subimage; /* Loop counter for MOSAIC images per file */
337 int iimage; /* Loop counter for all files/images */
338 int num_images; /* Total number of slices (>= # files) */
339 Mosaic_Info mi; /* Mosaic (multi-image) information */
340 Multiframe_Info mfi; /* Multiframe information */
341 int n_slices_in_file; /* Number of slices in file */
342 int subimage_type; /* Subimage type */
343
344 subimage_type = SUBIMAGE_TYPE_NONE;
345
346 /* Allocate space for the file information */
347 fi_ptr = malloc(num_files * sizeof(*fi_ptr));
348 CHKMEM(fi_ptr);
349
350 num_images = num_files;
351
352 /* Last group needed for first pass
353 */
354 max_group = ACR_IMAGE_GID - 1;
355
356 /* Add all control characters as numeric array separators to handle
357 * odd behaviour with Siemens dicom files
358 */
359 for (isep = 0; isep < 31; isep++) {
360 acr_element_numeric_array_separator(isep);
361 }
362
363 /* Initialize some values for general info */
364 gi.initialized = FALSE;
365 gi.group_list = NULL;
366 for (imri = 0; imri < MRI_NDIMS; imri++) {
367 gi.cur_size[imri] = -1;
368 gi.indices[imri] = NULL;
369 gi.coordinates[imri] = NULL;
370 gi.widths[imri] = NULL;
371 }
372
373 /* Loop through file list getting information
374 * (note that we have to duplicate the handling
375 * of multiple images per file in this loop
376 * to accumulate dimension sizes correctly)
377
378 * need separate counter for images, since some files may
379 * contain more than one image!
380 */
381 iimage = 0;
382
383 for (ifile = 0; ifile < num_files; ifile++) {
384 if (G.Debug >= HI_LOGGING) {
385 printf("\nFile %s\n", file_list[ifile]);
386 }
387
388 if (!G.Debug) {
389 progress(ifile, num_files, "-Parsing series info");
390 }
391
392 /* Read the file
393 */
394 if (G.file_type == N4DCM) {
395 group_list = read_numa4_dicom(file_list[ifile], max_group);
396 }
397 else if (G.file_type == IMA) {
398 group_list = siemens_to_dicom(file_list[ifile], max_group);
399 }
400
401 if (group_list == NULL) {
402 fprintf(stderr, "Error parsing file '%s' on 1st pass.\n",
403 file_list[ifile]);
404 exit(-1);
405 }
406
407 if (G.opts & OPTS_NO_MOSAIC) {
408 n_slices_in_file = 1;
409 }
410 else {
411 n_slices_in_file = acr_find_int(group_list, EXT_Slices_in_file, 1);
412 }
413
414 /* initialize big and small images, if mosaic
415 */
416 if (n_slices_in_file > 1) {
417
418 subimage_type = SUBIMAGE_TYPE_MOSAIC;
419
420 mosaic_init(group_list, &mi, FALSE);
421
422 num_images += n_slices_in_file - 1;
423
424 fi_ptr = realloc(fi_ptr, num_images * sizeof(*fi_ptr));
425 CHKMEM(fi_ptr);
426 }
427 else {
428 /* See if we have an DICOM multiframe image, or at least a
429 * reasonable facsimile of one such as the files produced
430 * by the Shimadzu PET scanner.
431 *
432 * It is probably not correct to rely on the presence of the
433 * "Number of Frames" (0x0028, 0x0008) field, but the Shimadzu
434 * multiframe implementation is far from complete and does not
435 * seem to implement any of the other components of the
436 * specification.
437 */
438
439 n_slices_in_file = acr_find_int(group_list, ACR_Number_of_frames,
440 1);
441
442 if (n_slices_in_file > 1) {
443 subimage_type = SUBIMAGE_TYPE_MULTIFRAME;
444
445 multiframe_init(group_list, &mfi, FALSE);
446
447 num_images += n_slices_in_file - 1;
448
449 fi_ptr = realloc(fi_ptr, num_images * sizeof(*fi_ptr));
450 CHKMEM(fi_ptr);
451 }
452 }
453
454
455 /* loop over subimages in mosaic
456 */
457 for (subimage = 0; subimage < n_slices_in_file; subimage++) {
458
459 /* Modify the group list for this image if mosaic or multiframe
460 */
461 switch (subimage_type) {
462 case SUBIMAGE_TYPE_MOSAIC:
463 mosaic_insert_subframe(group_list, &mi, subimage, FALSE);
464 break;
465 case SUBIMAGE_TYPE_MULTIFRAME:
466 multiframe_insert_subframe(group_list, &mfi, subimage, FALSE);
467 break;
468 default:
469 break;
470 }
471
472 /* Get file-specific information
473 */
474 get_file_info(group_list, &fi_ptr[iimage], &gi);
475
476 //ilana debug
477 /*int acq=acr_find_int(group_list, ACR_Acquisition, 1);
478 printf("****WE HAVE ACQUISITION# %i",acq);*/
479
480 /* increment iimage here
481 */
482 iimage++;
483 }
484
485 /* Delete the group list
486 */
487 acr_delete_group_list(group_list);
488
489 /* cleanup mosaic struct if used
490 */
491 switch (subimage_type) {
492 case SUBIMAGE_TYPE_MOSAIC:
493 mosaic_cleanup(&mi);
494 break;
495
496 case SUBIMAGE_TYPE_MULTIFRAME:
497 multiframe_cleanup(&mfi);
498 break;
499
500 default:
501 break;
502 }
503 }
504
505 /* Sort the dimensions */
506 sort_dimensions(&gi);
507
508 /* Create the output file
509 */
510 if (gi.initialized) {
511 icvid = create_minc_file(minc_file,
512 clobber,
513 &gi,
514 file_prefix,
515 &out_file_name,
516 loop_type);
517 }
518 if (output_file_name != NULL) {
519 *output_file_name = out_file_name;
520 }
521
522 /* Check that we found the general info and that the minc file was
523 * created okay
524 */
525 if ((!gi.initialized) || (icvid == MI_ERROR)) {
526 if (gi.initialized) {
527 fprintf(stderr, "Error creating MINC file %s.\n", out_file_name);
528 }
529 free_info(&gi, fi_ptr, num_files);
530 free(fi_ptr);
531 return EXIT_FAILURE;
532 }
533
534 if (G.Debug) {
535 printf("Writing %d images to MINC file\n", num_files);
536 }
537
538 /* Last group needed for second pass
539 * we now have to read up to and including the image,
540 * since image pointers are needed in mosaic_init
541 */
542 max_group = ACR_IMAGE_GID;
543
544 /* Loop through the files again and put images into the minc file
545 */
546 iimage = 0;
547 for (ifile = 0; ifile < num_files; ifile++) {
548
549 if (!G.Debug) {
550 progress(ifile, num_files, "-Creating minc file");
551 }
552
553 /* Check that we have a valid file
554 */
555 if (!fi_ptr[ifile].valid) {
556 printf("WARNING: file %s was marked invalid\n",
557 file_list[ifile]);
558 continue;
559 }
560
561 /* Read the file
562 */
563 if (G.file_type == N4DCM) {
564 group_list = read_numa4_dicom(file_list[ifile], max_group);
565 }
566 else if (G.file_type == IMA) {
567 group_list = siemens_to_dicom(file_list[ifile], max_group);
568 }
569
570 if (group_list == NULL) {
571 fprintf(stderr, "Error parsing file '%s' during 2nd pass.\n",
572 file_list[ifile]);
573 exit(-1);
574 }
575
576 /* initialize big and small images, if mosaic
577 */
578 switch (subimage_type) {
579 case SUBIMAGE_TYPE_MOSAIC:
580 mosaic_init(group_list, &mi, TRUE);
581 break;
582
583 case SUBIMAGE_TYPE_MULTIFRAME:
584 multiframe_init(group_list, &mfi, TRUE);
585 break;
586
587 default:
588 break;
589 }
590
591 /* loop over subimages in mosaic
592 */
593 for (subimage = 0; subimage < n_slices_in_file; subimage++) {
594
595 /* Modify the group list for this image if mosaic
596 */
597 switch (subimage_type) {
598 case SUBIMAGE_TYPE_MOSAIC:
599 mosaic_insert_subframe(group_list, &mi, subimage, TRUE);
600 break;
601
602 case SUBIMAGE_TYPE_MULTIFRAME:
603 multiframe_insert_subframe(group_list, &mfi, subimage, TRUE);
604 break;
605
606 default:
607 break;
608 }
609
610 /* Get image
611 */
612 get_dicom_image_data(group_list, &image);
613
614 /* Save the image and any other information
615 */
616 save_minc_image(icvid, &gi, &fi_ptr[iimage], &image);
617
618 /* Free the image data */
619 if (image.data != NULL) {
620 free(image.data);
621 image.data = NULL;
622 }
623
624 /* increment image counter
625 */
626 iimage++;
627 }
628
629 /* Delete the group list
630 */
631 acr_delete_group_list(group_list);
632
633 /* cleanup mosaic struct if used
634 */
635 switch (subimage_type) {
636 case SUBIMAGE_TYPE_MOSAIC:
637 mosaic_cleanup(&mi);
638 break;
639
640 case SUBIMAGE_TYPE_MULTIFRAME:
641 multiframe_cleanup(&mfi);
642 break;
643
644 default:
645 break;
646 }
647 }
648
649 /* Close the output file */
650 close_minc_file(icvid);
651
652 /* Free the gi and fi_ptr stuff */
653 free_info(&gi, fi_ptr, num_files);
654 free(fi_ptr);
655
656 return EXIT_SUCCESS;
657 }
658
659 /* ----------------------------- MNI Header -----------------------------------
660 @NAME : read_std_dicom
661 @INPUT : filename - name of siemens Numaris 4 `dicom' file to read
662 max_group - maximum group number to read
663 @OUTPUT : (none)
664 @RETURNS : group list read in from file
665 @DESCRIPTION: Routine to read in a group list from a file.
666 @METHOD :
667 @GLOBALS :
668 @CALLS :
669 @CREATED : December 18, 2001 (Rick Hoge)
670 @MODIFIED :
671 ---------------------------------------------------------------------------- */
672
673 Acr_Group
read_std_dicom(const char * filename,int max_group)674 read_std_dicom(const char *filename, int max_group)
675 {
676 FILE *fp;
677 Acr_File *afp;
678 Acr_Group group_list;
679 int status;
680
681 /* Open the file
682 */
683 fp = fopen(filename, "rb");
684 if (fp == NULL) {
685 return NULL;
686 }
687
688 /* Connect to input stream
689 */
690 afp = acr_file_initialize(fp, 0, acr_stdio_read);
691 if (afp == NULL) {
692 return NULL;
693 }
694
695 acr_set_ignore_errors(afp, 1); /* ignore protocol errors */
696
697 if (acr_test_dicom_file(afp) != ACR_OK) {
698 return NULL;
699 }
700
701 // Read in group list
702 status = acr_input_group_list(afp, &group_list, max_group);
703 if (status != ACR_END_OF_INPUT && status != ACR_OK) {
704 return NULL;
705 }
706
707 // Close the file
708 acr_file_free(afp);
709 fclose(fp);
710
711 return (group_list);
712 }
713
714 /* Return non-zero if this image looks like it is a mosaic image.
715 */
716 int
is_siemens_mosaic(Acr_Group group_list)717 is_siemens_mosaic(Acr_Group group_list)
718 {
719 Acr_String str_ptr;
720
721 str_ptr = acr_find_string(group_list, ACR_Image_type, "");
722 if (strstr(str_ptr, "MOSAIC") != NULL)
723 return 1; /* slam dunk, this is mosaic data */
724
725 /* For some reason, some Localizer scans look like MOSAICs in that
726 * they appear to be upsampled in such a way that their
727 * rows/columns are larger than the acquisition matrix. The
728 * details of this are probably buried in the Siemens ASCCONV dump
729 * and therefore I have not been able to figure out how to tweak
730 * the MOSAIC stuff to get this right. Instead I am taking the
731 * cheesy way out and declaring that anything that is NOT
732 * explictly marked MOSAIC, and IS marked a localizer, should not
733 * be treated as a mosaic.
734 */
735 str_ptr = acr_find_string(group_list, ACR_Series_description, "");
736 if (strstr(str_ptr, "localizer") || strstr(str_ptr, "LOCALIZER")) {
737 return 0;
738 }
739
740 /* OK, we did not find the word "mosaic" in the image type. But this
741 * could still be a mosaic image. Let's look for some other clues.
742 */
743
744 /*the SPI_Number_of_slices_nominal and SPI_Number_of_raw_partitions_nomimal
745 are not always good fields with which to determine whether it's Mosaic
746 or not, these fields often get set whether we have a mosaic image or not*/
747
748
749 /*if (acr_find_int(group_list, SPI_Number_of_slices_nominal, 0) > 1)
750 return 1; */ /* probably mosaic. */
751
752 /*if (acr_find_int(group_list, SPI_Number_of_3D_raw_partitions_nominal, 0) > 1)
753 return 1; */ /* probably mosaic */
754
755 return 0; /* probably not mosaic */
756 }
757
758 #define MAXVM 32
759
760 Acr_Group
parse_siemens_proto2(Acr_Group group_list,Acr_Element element)761 parse_siemens_proto2(Acr_Group group_list, Acr_Element element)
762 {
763 int byte_cnt;
764 int byte_pos;
765 unsigned char *byte_ptr;
766 Acr_Long n_items;
767 Acr_Long n_values;
768 Acr_Long vm;
769 char vr[4];
770 Acr_Long syngodt;
771 int i;
772 char name[64];
773 char *value[MAXVM];
774 Acr_Long len;
775
776 byte_cnt = element->data_length;
777 byte_ptr = element->data_pointer;
778 byte_pos = 0;
779
780 /* See if the magic 8-byte header is present. If not, we start
781 * right off with the number of elements.
782 */
783 if (byte_ptr[0] == 'S' && byte_ptr[1] == 'V' &&
784 byte_ptr[2] == '1' && byte_ptr[3] == '0') {
785 byte_pos += 8; /* Skip header */
786 }
787
788 acr_get_long(ACR_LITTLE_ENDIAN, 1, byte_ptr + byte_pos, &n_items);
789 byte_pos += 8; /* Skip # items and 4 bytes of unknown junk */
790
791 while (n_items-- > 0) {
792 strncpy(name, (byte_ptr + byte_pos), 64);
793 byte_pos += 64;
794
795 acr_get_long(ACR_LITTLE_ENDIAN, 1, byte_ptr + byte_pos, &vm);
796 byte_pos += 4;
797
798 strncpy(vr, (byte_ptr + byte_pos), 4);
799 byte_pos += 4;
800
801 acr_get_long(ACR_LITTLE_ENDIAN, 1, byte_ptr + byte_pos, &syngodt);
802 byte_pos += 4;
803
804 acr_get_long(ACR_LITTLE_ENDIAN, 1, byte_ptr + byte_pos, &n_values);
805 byte_pos += 4;
806
807 byte_pos += 4; /* skip dummy */
808
809 if (n_values > 0) {
810 for (i = 0; i < (int)n_values; i++) {
811 byte_pos += 4; /* skip */
812
813 acr_get_long(ACR_LITTLE_ENDIAN, 1, byte_ptr + byte_pos, &len);
814
815 byte_pos += 4;
816
817 byte_pos += 8;
818
819 if (i < vm && i < MAXVM) {
820 value[i] = (byte_ptr + byte_pos);
821 }
822
823 byte_pos += ((len + 3) / 4) * 4;
824 }
825 }
826
827 if (!strcmp(name, "DiffusionDirectionality")) {
828 acr_insert_string(&group_list, ACR_Diffusion_directionality,
829 value[0]);
830 }
831 else if (!strcmp(name, "B_value")) {
832
833 /*double tmp = atof(value[0]); this atof makes the value null! ilana*/
834 Acr_Double tmp = atoi(value[0]);
835 /*need a hack for ICBM scan, see below ilana*/
836 acr_insert_double(&group_list, ACR_Diffusion_b_value, 1, &tmp);
837 }
838 else if (!strcmp(name, "DiffusionGradientDirection")) {
839 Acr_Double tmp[3];
840 if (vm == 3 && n_values >= vm) {
841 /*For the ICBM WIP scan, the b0 images do not have
842 B_value=0 or DiffusionGradientDirection=0 0 0, they
843 actually have DiffusionGradientDirection=-1.00010000 -1.00010000 -1.00010000
844 Use this to detect the b=0 images and hopefully this doesn't ever correspond
845 to a real gradient direction (have to change bvalues correspondingly)! ilana*/
846
847 if(!strcmp(value[0],"-1.00010000") && !strcmp(value[1],"-1.00010000") && !strcmp(value[2],"-1.00010000")){
848 value[0] = value[1] = value[2] = "0"; /*grad directions should be 0*/
849 Acr_Double tmp2 = atoi("0");
850 acr_insert_double(&group_list, ACR_Diffusion_b_value, 1, &tmp2); /*also have to modify B values ilana*/
851 }
852
853 tmp[0] = atof(value[0]);
854 tmp[1] = atof(value[1]);
855 tmp[2] = -1*atof(value[2]); /*dicom z = -minc z ilana*/
856
857 acr_insert_double(&group_list,
858 ACR_Diffusion_gradient_orientation,
859 3,
860 tmp);
861 }
862 }
863
864 #if 0
865 if (G.Debug >= HI_LOGGING) {
866 printf("%s VM=%d VR=%s %d ", name, vm, vr, n_values);
867
868 if (n_values != 0) {
869 for (i = 0; i < vm && i < MAXVM; i++) {
870 printf("%s ", value[i]);
871 }
872 }
873 printf("\n");
874 }
875 #endif
876 }
877 return (group_list);
878 }
879
880 static Acr_Group
add_siemens_info(Acr_Group group_list)881 add_siemens_info(Acr_Group group_list)
882 {
883 /* needed for group repair - some essential info
884 * only available in ascii dump of MrProt structure
885 */
886 Acr_Element protocol;
887 Acr_Element element;
888 int num_slices, num_partitions;
889 string_t str_buf;
890 char *str_ptr=NULL;
891 char *str_ptr2=NULL;
892 int interpolation_flag;
893 int enc_ix, num_encodings, num_b0;
894 int EXT=0; /*special handling when an external diffusion vector file is used*/
895
896 element = acr_find_group_element(group_list, SPI_Protocol2);
897 if (element != NULL) {
898 group_list = parse_siemens_proto2(group_list, element);
899 }
900
901 /* now fix the group list to provide essential info
902 * via standard dicom elements
903 * (this lets the rest of the code be more reusable)
904
905 * Note that these parameters are mostly dimension lengths,
906 * and are not usually supplied in standard DICOM implementations.
907 * We could do without them, except that the use of mosaics for
908 * multi-slice data makes at least the number of slices necessary.
909 * For such elements, we will spoof our own `private' entries
910 * using Numaris 3.5 coordinates. These should not be used elsewhere
911 * in the code unless there's no other way! Basically things
912 * should be done as follows:
913
914 * 1) use actual element in public dicom group, if possible
915 * 2) if not, then get info from MrProt and insert as
916 * correct public dicom element
917 * 3) if no public element exists, use an SPI element (careful!)
918 * 4) if no SPI element exists, use and EXT element (careful!)
919 */
920
921 /* read in Protocol group */
922
923 protocol = acr_find_group_element(group_list, SPI_Protocol);
924 if (protocol != NULL) {
925 /* parse_siemens_junk(protocol); */
926
927 if (G.Debug >= HI_LOGGING) {
928 printf("Incorporating Siemens protocol structure...\n");
929 }
930
931 /* Add number of dynamic scans:
932 */
933 prot_find_string(protocol, "lRepetitions", str_buf);
934
935 if(atoi(str_buf)==0){ /*lRepetitions not found in ASCONV header*/
936 int frames = acr_find_int(group_list,
937 ACR_Cardiac_number_of_images,0);
938 /* this seems to give number of dynamic
939 scans when lRepetitions not there*/
940 acr_insert_numeric(&group_list, ACR_Acquisitions_in_series, frames);
941 } else {
942 acr_insert_numeric(&group_list, ACR_Acquisitions_in_series,
943 atoi(str_buf) + 1);
944 }
945
946 /* add number of echoes:
947 */
948 prot_find_string(protocol, "lContrasts", str_buf);
949 acr_insert_numeric(&group_list, SPI_Number_of_echoes, (double)atoi(str_buf));
950
951 /* Add receiving coil (for some reason this isn't in generic groups)
952 */
953 prot_find_string(protocol,
954 "sCOIL_SELECT_MEAS.asList[0].sCoilElementID.tCoilID",
955 str_buf);
956
957 /*Adjust for change in naming convention (VB15 VA30)*/
958 if(atoi(str_buf)==0){
959 prot_find_string(protocol,
960 "sCoilSelectMeas[0].asList[0].sCoilElementID.tCoilID",
961 str_buf);
962 }
963 acr_insert_string(&group_list, ACR_Receive_coil_name, str_buf);
964
965 /* add MrProt dump
966 */
967 str_ptr = dump_protocol_text(protocol);
968 acr_insert_string(&group_list, EXT_MrProt_dump, (Acr_String)str_ptr);
969 free(str_ptr);
970
971 /* add number of slices: (called `Partitions' for 3D) */
972 prot_find_string(protocol, "sSliceArray.lSize", str_buf);
973 num_slices = atoi(str_buf);
974
975 prot_find_string(protocol, "sKSpace.lPartitions", str_buf);
976 num_partitions = atoi(str_buf);
977
978 /* This is a hack based upon the observation that for at least some
979 * conversions, this value seems to give the true number of slices
980 * rather than the sKSpace.lPartitions value (bert)
981 */
982 prot_find_string(protocol, "sKSpace.lImagesPerSlab", str_buf);
983 if (str_buf[0] != '\0') {
984 int num_images_per_slab = atoi(str_buf);
985 if (num_images_per_slab > num_partitions) {
986 num_partitions = num_images_per_slab;
987 }
988 }
989
990 /* NOTE: for some reason, lPartitions > 1 even for 2D scans
991 * (e.g. EPI, scouts)
992 */
993 if (!strncmp(acr_find_string(group_list, ACR_MR_acquisition_type,""),
994 "3D", 2)) {
995 /* Use partitions if 3D.
996 * (note that this gets more complicated if the 3D scan
997 * is mosaiced - see below)
998 */
999 acr_insert_numeric(&group_list, SPI_Number_of_slices_nominal, (double)1);
1000 acr_insert_numeric(&group_list,
1001 SPI_Number_of_3D_raw_partitions_nominal,
1002 (double)num_partitions);
1003 }
1004 else {
1005 /* use slices for 2D
1006 */
1007 acr_insert_numeric(&group_list, SPI_Number_of_slices_nominal,
1008 (double)num_slices);
1009 acr_insert_numeric(&group_list,
1010 SPI_Number_of_3D_raw_partitions_nominal, (double)1);
1011 }
1012
1013 /* See if the field of view is set, if not, see if we can find
1014 * it in the ASCCONV dump and copy it into the proprietary
1015 * fields.
1016 */
1017 if (acr_find_group_element(group_list, SPI_Field_of_view) == NULL) {
1018 /* TODO: This assumes a symmetric field of view. We need to
1019 * figure out what to do if this is not true!!
1020 */
1021 prot_find_string(protocol, "sSliceArray.asSlice[0].sReadoutFOV",
1022 str_buf);
1023 if (str_buf[0] != '\0') {
1024 int fov = atoi(str_buf);
1025
1026 sprintf(str_buf, "%d\\%d", fov, fov);
1027 acr_insert_string(&group_list, SPI_Field_of_view, str_buf);
1028 }
1029 }
1030
1031 prot_find_string(protocol, "sKSpace.uc2DInterpolation", str_buf);
1032 interpolation_flag = strtol(str_buf, NULL, 0);
1033
1034 /* find Delay time in TR in ASCONV header ilana */
1035 prot_find_string(protocol, "lDelayTimeInTR", str_buf);
1036 acr_insert_numeric(&group_list, EXT_Delay_in_TR,
1037 (double)atol(str_buf));
1038
1039 /* Find slice acquisition mode in ASCONV header,
1040 can't seem to find is in a standard dicom field
1041 sSliceArray.ucMode takes on 3 values:
1042 0x1 means ASCENDING
1043 0x2 means DESCENDING
1044 0x4 means INTERLEAVED */
1045 prot_find_string(protocol, "sSliceArray.ucMode", str_buf);
1046 acr_insert_string(&group_list,EXT_Slice_order,str_buf);
1047
1048 /*Modified by ilana to handle the common types of diffusion scans (ref: siemens_dicom_to_minc for dicomserver)
1049
1050 /* correct dynamic scan info if *MGH* diffusion scan:
1051 *
1052 * assumptions:
1053 *
1054 * - diffusion protocol indicated by sDiffusion.ulMode = 0x100
1055 * - bvalue is in sDiffusionalBValue[1]
1056 * - there are 10 b=0 scans and a user defined number of diffusion directions
1057 * - b=0 scans have sequence name "ep_b0#0, ep_b0#1... etc"
1058 * - encoded scans have seq names "ep_b1000#1, ep_b1000#2, ...,ep_b1300#1, ep_b1300#2, ..." etc.
1059 *
1060 * actions:
1061 *
1062 * - change number of dynamic scans to sDiffusion.lDiffDirections + num b0 images
1063 * - use sWiPMemBlock.alFree[8] for number of b=0 scans
1064 * - modify dynamic scan index to encoding index
1065 */
1066 /* correct dynamic scan info if standard *ep2d_diff* diffusion scan
1067 *
1068 * assumptions:
1069 *
1070 * - diffusion protocol indicated by sDiffusion.ulMode = 0x100
1071 * - bvalue is in sDiffusion.alBValue[0]
1072 * - there is 1 b=0 scans and 12 diffusion directions
1073 * - b=0 scan havs sequence name "ep_b0"
1074 * - encoded scans have seq names "ep_b1000#1, ep_b1000#2, ..." etc.
1075 *
1076 * actions:
1077 *
1078 * - change number of dynamic scans to sDiffusion.lDiffDirections + num b0 images
1079 * - modify dynamic scan index to encoding index
1080 */
1081 /* correct dynamic scan info if standard *ICBM_WIP* diffusion scan
1082 *
1083 * assumptions:
1084 *
1085 * - diffusion protocol indicated by sDiffusion.ulMode = 0x80
1086 * - bvalue is in sDiffusion.alBValue[1]
1087 * - there is 1 b=0 scans and user defined number of diffusion directions
1088 * - b=0 scan has sequence name "ep_b0"
1089 * - encoded scans have seq names "ep_b1000#1, ep_b1000#2, ..." etc.
1090 *
1091 * actions:
1092 *
1093 * - change number of dynamic scans to sDiffusion.lDiffDirections + num b0 images
1094 * - modify dynamic scan index to encoding index
1095 */
1096
1097 prot_find_string(protocol, "sDiffusion.ulMode", str_buf);
1098 if (!strcmp(str_buf, "0x100") | !strcmp(str_buf, "0x80")) {
1099 /*we have a diffusion scan; flag it*/
1100 acr_insert_string(&group_list, ACR_Acquisition_contrast, "DIFFUSION");
1101 /*----MGH-----*/
1102 prot_find_string(protocol,"sWiPMemBlock.alFree[8]", str_buf);
1103 if ((atoi ((char*)str_buf))!= 0 ) { /*num b0 images for MGH sequence*/
1104
1105 /* get number of b=0 images*/
1106 num_b0 = atoi ((char*)str_buf);
1107
1108 /* try to get b value */
1109 prot_find_string(protocol, "sDiffusion.alBValue[1]", str_buf);
1110
1111 acr_insert_numeric(&group_list, EXT_Diffusion_b_value,
1112 (double)atoi(str_buf));
1113
1114 }
1115 else
1116 {
1117 prot_find_string(protocol, "sDiffusion.ulMode", str_buf);
1118 /*-----ep2d_diff-----OR-----
1119 ep2d_diff_WIP_ICBM with "Diffusion mode"=MDDW & DiffusionWeightings=2-----------*/
1120 if (!strcmp(str_buf, "0x100")) {
1121 /* try to get b value */
1122 prot_find_string(protocol, "sDiffusion.alBValue[1]", str_buf);
1123 acr_insert_numeric(&group_list, EXT_Diffusion_b_value,
1124 (double)atoi(str_buf));
1125 num_b0=1;
1126
1127 }
1128 /*-----ICBM_WIP with "Diffusion mode"=Free (5 b=0 scans) or any time an external vectors file is used-----*/
1129 else if(!strcmp(str_buf, "0x80")) {
1130 EXT=1;
1131 /* try to get b value */
1132 prot_find_string(protocol, "sDiffusion.alBValue[0]", str_buf);
1133 acr_insert_numeric(&group_list, EXT_Diffusion_b_value,
1134 (double)atoi(str_buf));
1135 num_b0=0;
1136 /*For ICBM there are 5 b=0 scans but they are not identified
1137 any differently than the diffusion weighted images,
1138 sDiffusion.lDiffDirections includes the b=0 images. An external
1139 DiffusionVectors file can include other b=0 and this messes up
1140 the image count*/
1141 }
1142 }
1143
1144 /* if all averages in one series: */
1145 prot_find_string(protocol,"ucDixon",str_buf);
1146 if (!strcmp(str_buf,"0x1")) {
1147 prot_find_string(protocol, "sDiffusion.lDiffDirections", str_buf);
1148 num_encodings = atoi ((char*)str_buf) + num_b0;
1149
1150 /* number of 'time points' */
1151 acr_insert_numeric(&group_list, ACR_Acquisitions_in_series,
1152 num_encodings *
1153 acr_find_double(group_list, ACR_Nr_of_averages, 1));
1154
1155 /* time index of current scan: */
1156
1157 /* For multi-series scans, we DO USE THIS BECAUSE global
1158 * image number may be broken!!
1159 */
1160
1161 /*Have to also take care of numbered b=0 images (ep_b0#0, etc...) ilana*/
1162 /*the Siemems-based sequences, in MDDW mode with 2 diff weightings have b=0 images called ep_b0*/
1163
1164 str_ptr = strstr(acr_find_string(group_list, ACR_Sequence_name, ""), "b");
1165 str_ptr2 = strstr(acr_find_string(group_list, ACR_Sequence_name, ""), "#");
1166
1167 if (str_ptr == NULL || str_ptr2 == NULL) {
1168 enc_ix = 0;
1169 }
1170 else if(atoi(str_ptr +1) == 0){ /*a 0 after the b means b=0 image*/
1171 enc_ix = atoi(str_ptr2 + 1);
1172 }
1173 else{
1174 enc_ix = atoi(str_ptr2 + 1) + num_b0; /*should be in diffusion weighted images now*/
1175 }
1176 /* however with the current sequence, we get usable
1177 * time indices from floor(global_image_num/num_slices)*/ /*i'm not sure that works here*/
1178
1179 acr_insert_numeric(&group_list, ACR_Acquisition, (double)enc_ix);
1180 /*acr_insert_numeric(&group_list, ACR_Acquisition,
1181 (acr_find_int(group_list, ACR_Image, 1)-1) /
1182 num_slices);*/
1183
1184 }
1185 else { /* averages in different series - no special handling needed? */
1186
1187 prot_find_string(protocol, "sDiffusion.lDiffDirections", str_buf);
1188 num_encodings = atoi((char*)str_buf) + num_b0;
1189 //num_encodings = 7; /* for now assume 7 shots in diffusion scan */
1190
1191 /* number of 'time points' */
1192 acr_insert_numeric(&group_list, ACR_Acquisitions_in_series,
1193 (double)num_encodings);
1194
1195 /* For multi-series scans, we DO USE THIS BECAUSE global
1196 * image number may be broken!!
1197 */
1198 /*Have to also take care of numbered b=0 images (ep_b0#0, etc...) ilana*/
1199
1200 str_ptr = strstr(acr_find_string(group_list, ACR_Sequence_name, ""), "b");
1201 str_ptr2 = strstr(acr_find_string(group_list, ACR_Sequence_name, ""), "#");
1202
1203 if (str_ptr == NULL || str_ptr2 == NULL){
1204 enc_ix = 0;
1205 }
1206 else if(atoi(str_ptr +1) == 0){ /*a 0 after the b means b=0 image*/
1207 enc_ix = atoi(str_ptr2 + 1);
1208 }
1209 else{
1210 enc_ix = atoi(str_ptr2 + 1) + num_b0; /*should be in diffusion weighted images now*/
1211 }
1212
1213 acr_insert_numeric(&group_list, ACR_Acquisition, (double)enc_ix);
1214 }
1215 if (EXT==1) {
1216 /*if an external DiffusionVectors was used, the encoding index can be wrong
1217 because it does not take into account b=0 images within the acquisition.
1218 Have to rely on ACR_Image 0020x0013 (but this field won't work for MGH sequence)*/
1219 acr_insert_numeric(&group_list, ACR_Acquisition,acr_find_int(group_list, ACR_Image, 1) );
1220 }
1221
1222 /* BUG! TODO! FIXME! In dcm2mnc.c the sequence name is
1223 * used as one of the criteria for starting a new
1224 * file. For a DTI sequence, we don't want this to
1225 * happen. For now I replace the sequence name with a
1226 * bogus value in order to keep the DTI scan from being
1227 * split into pieces. This isn't right. We should deal
1228 * with this in a more graceful way.
1229 */
1230 acr_insert_string(&group_list, ACR_Sequence_name, "DTI");
1231
1232 /* Reset acquisition time to series time plus the series
1233 * index. Otherwise each slice may get its own
1234 * time. This is also probably wrong and in the ideal
1235 * world we should deal with this gracefully. But we have
1236 * no good way of storing per-slice timing information
1237 * right now.
1238 */
1239 acr_insert_numeric(&group_list, ACR_Acquisition_time, acr_find_double(group_list, ACR_Series_time, 0) + enc_ix);
1240 /* end of diffusion scan handling */
1241
1242
1243 /* There is another whole class of (probably newer) diffusion
1244 * weighted images that use a very different arrangement and
1245 * need better handling.
1246 */
1247
1248 }
1249 }
1250 else {
1251 /* If no protocol dump was found we have to assume no
1252 * interpolation since at the momemnt I have no idea where it
1253 * would live...
1254 */
1255 interpolation_flag = 0;
1256 num_partitions = acr_find_int(group_list,
1257 SPI_Number_of_3D_raw_partitions_nominal, 0);
1258 num_slices = acr_find_int(group_list, SPI_Number_of_slices_nominal, 0);
1259 }
1260
1261 /* Handle mosaic images if necessary.
1262 */
1263 if (is_siemens_mosaic(group_list)) {
1264
1265 int mosaic_rows, mosaic_cols;
1266 Acr_Short subimage_size[4];
1267 int subimage_rows, subimage_cols;
1268
1269 /* Now figure out mosaic rows and columns, and put in EXT
1270 * shadow group. Check for interpolation - will require 2x
1271 * scaling of rows and columns.
1272 */
1273
1274 /* Assign defaults in case something goes wrong below...
1275 */
1276 subimage_rows = subimage_cols = acr_find_int(group_list,
1277 SPI_Base_raw_matrix_size,
1278 0);
1279
1280 /* If we don't have an "SPI_Base_raw_matrix_size" field, try using
1281 * the acquisition matrix.
1282 */
1283 if (subimage_rows == 0) {
1284
1285 /* Compute mosaic rows and columns
1286 * Here is a hack to handle non-square mosaiced images
1287 *
1288 * WARNING: as far as I can tell, the phase-encoding dir (row/col)
1289 * is reversed for mosaic EPI scans (don't know if this is a
1290 * mosaic thing, an EPI thing, or whatever).
1291 * (Something is wrong here, it seems that it's not reversed, but
1292 * not sure obviously doesn't show up when the acquisition matrix is
1293 * square)
1294 * Get the array of sizes: freq row/freq col/phase row/phase col
1295 */
1296
1297 element = acr_find_group_element(group_list,
1298 ACR_Acquisition_matrix);
1299 if (element == NULL) {
1300 printf("WARNING: Can't find acquisition matrix\n");
1301 }
1302 else if (acr_get_element_short_array(element, 4,
1303 subimage_size) != 4) {
1304 printf("WARNING: Can't read acquisition matrix\n");
1305 }
1306 else {
1307 if (G.Debug >= HI_LOGGING) {
1308 printf(" * Acquisition matrix %d %d %d %d\n",
1309 subimage_size[0],
1310 subimage_size[1],
1311 subimage_size[2],
1312 subimage_size[3]);
1313 }
1314 /* Get subimage dimensions, assuming the OPPOSITE of the
1315 * reported phase-encode direction!!
1316 */
1317 str_ptr = (char *)acr_find_string(group_list,
1318 ACR_Phase_encoding_direction,
1319 "");
1320 if (!strncmp(str_ptr, "COL", 3)) { /*TODO test that this is right for non-square acquisition matrices*/
1321 subimage_rows = subimage_size[3];
1322 subimage_cols = subimage_size[0];
1323 }
1324 else if (!strncmp(str_ptr, "ROW", 3)) {
1325 /*not sure the subimage dimensions should be opposite of
1326 * the phase encoding direction, running into some problems,
1327 * so changed it to be the same but I also might be breaking
1328 * for other cases!*/
1329 /*subimage_rows = subimage_size[2];
1330 subimage_cols = subimage_size[1];*/
1331 subimage_rows = subimage_size[1];
1332 subimage_cols = subimage_size[2];
1333 }
1334 else {
1335 printf("WARNING: Unknown phase encoding direction '%s'\n",
1336 str_ptr);
1337 }
1338
1339 /* If interpolation, multiply rows and columns by 2 */
1340 if (interpolation_flag) {
1341 subimage_rows *= 2;
1342 subimage_cols *= 2;
1343 }
1344 }
1345 }
1346
1347 /* If these are not set or still zero, assume this is NOT a mosaic
1348 * format file.
1349 */
1350 if (subimage_rows == 0 || subimage_cols == 0) {
1351 subimage_rows = acr_find_int(group_list, ACR_Rows, 1);
1352 subimage_cols = acr_find_int(group_list, ACR_Columns, 1);
1353 mosaic_rows = 1;
1354 mosaic_cols = 1;
1355 }
1356 else {
1357 if (G.Debug >= HI_LOGGING) {
1358 printf("Assuming %dx%d mosaic\n", subimage_rows, subimage_cols);
1359 }
1360 mosaic_rows = acr_find_int(group_list, ACR_Rows, 1) / subimage_rows;
1361 mosaic_cols = acr_find_int(group_list, ACR_Columns, 1) / subimage_cols;
1362 }
1363
1364 acr_insert_numeric(&group_list, EXT_Mosaic_rows, (double)mosaic_rows);
1365 acr_insert_numeric(&group_list, EXT_Mosaic_columns, (double)mosaic_cols);
1366
1367 if (mosaic_rows * mosaic_cols > 1) {
1368
1369 str_ptr = (char *)acr_find_string(group_list, ACR_MR_acquisition_type, "");
1370
1371 /* assume any mosaiced file contains all slices
1372 * (we now support mosaics for 2D and 3D acquisitions,
1373 * so we may need to use partitions instead of slices)
1374 */
1375
1376 if (!strncmp(str_ptr, "2D", 2)) {
1377 acr_insert_numeric(&group_list, EXT_Slices_in_file, (double)num_slices);
1378 acr_insert_numeric(&group_list,
1379 SPI_Number_of_slices_nominal, (double)num_slices);
1380 }
1381
1382 /* if 3D mosaiced scan, write number of partitions to number of
1383 * slices in dicom group
1384 */
1385 else if (!strncmp(str_ptr, "3D", 2)) {
1386 acr_insert_numeric(&group_list, EXT_Slices_in_file,
1387 (double)num_partitions);
1388 acr_insert_numeric(&group_list, SPI_Number_of_slices_nominal,
1389 (double)num_partitions);
1390 /* also have to provide slice spacing - in case of 3D it's same
1391 * as slice thickness (and not provided in dicom header!)
1392 */
1393 acr_insert_numeric(&group_list, ACR_Spacing_between_slices,
1394 acr_find_double(group_list,
1395 ACR_Slice_thickness, 1.0));
1396 }
1397 }
1398 else {
1399 acr_insert_numeric(&group_list, EXT_Slices_in_file, (double)1);
1400 }
1401
1402 /* Correct the rows and columns values -
1403 * These will reflect those of the subimages in the mosaics
1404 * NOT the total image dimensions
1405 */
1406 acr_insert_short(&group_list, EXT_Sub_image_columns, subimage_cols);
1407 acr_insert_short(&group_list, EXT_Sub_image_rows, subimage_rows);
1408
1409 /* TODO: should also correct the image position here? */
1410 }
1411
1412 /* Hacks for specific numaris3 sequences.
1413 */
1414 if (is_numaris3(group_list)) {
1415 int tmp;
1416
1417 str_ptr = (char *)acr_find_string(group_list, ACR_Sequence_name, "");
1418
1419 /* This first hack is for Sebastian's Numaris-3 DTI scans.
1420 */
1421 if (!strcmp(str_ptr, "ep_d33a ")) {
1422 /* Numaris-3 DTI scan. Individual scans are identified by
1423 * echo number. We need to change echos to time.
1424 */
1425 tmp = acr_find_int(group_list, SPI_Number_of_echoes, 1);
1426 acr_insert_numeric(&group_list, SPI_Number_of_echoes, (double)1);
1427
1428 element = acr_find_group_element(group_list,
1429 ACR_Acquisitions_in_series);
1430 if (element != NULL) {
1431 int grp_id = acr_get_element_group(element);
1432 int elm_id = acr_get_element_element(element);
1433 acr_group_remove_element(acr_find_group(group_list, grp_id),
1434 elm_id);
1435 }
1436 acr_insert_short(&group_list, ACR_Number_of_time_slices, tmp);
1437
1438 tmp = acr_find_int(group_list, ACR_Echo_number, 0);
1439 acr_insert_numeric(&group_list, ACR_Echo_number, (double)1);
1440 acr_insert_numeric(&group_list, ACR_Frame_reference_time,
1441 (double)(tmp * 1000));
1442
1443 /* Just use defaults for these values. */
1444 acr_insert_numeric(&group_list, ACR_Actual_frame_duration, (double)1000);
1445 }
1446
1447 /* An additional hack required for Sebastian's Numaris-3 FMRI scans.
1448 */
1449 if (!strcmp(str_ptr, "ep_fid")) {
1450 /* Numaris-3 FMRI scan. Each time slice uses a different
1451 * series number!!
1452 */
1453 tmp = acr_find_int(group_list, ACR_Series_time, 0);
1454 acr_insert_numeric(&group_list, ACR_Series, (double)tmp);
1455
1456 tmp = acr_find_int(group_list, ACR_Acquisitions_in_series, 0);
1457 acr_insert_short(&group_list, ACR_Number_of_time_slices, tmp);
1458
1459 element = acr_find_group_element(group_list,
1460 ACR_Acquisitions_in_series);
1461 if (element != NULL) {
1462 int grp_id = acr_get_element_group(element);
1463 int elm_id = acr_get_element_element(element);
1464 acr_group_remove_element(acr_find_group(group_list, grp_id),
1465 elm_id);
1466 }
1467 }
1468 }
1469 return (group_list);
1470 }
1471
1472 #define PMS_SET_CREATOR(el, cr) \
1473 (el)->element_id = ((el)->element_id & 0xff) + ((cr)->element_id << 8)
1474
1475 static Acr_Group
add_philips_info(Acr_Group group_list)1476 add_philips_info(Acr_Group group_list)
1477 {
1478 struct Acr_Element_Id creator_id;
1479 Acr_Group pms_group;
1480 Acr_Element pms_element;
1481 Acr_Element pms_element_list;
1482 char *str_ptr;
1483 int slice_count;
1484 int slice_index;
1485 char str_buf[128];
1486
1487 /* To use the Philips proprietary group, we have to figure out the
1488 * DICOM private creator ID in use. The group ID is always 0x2001,
1489 * but the upper eight bits of the element ID are somewhat variable.
1490 * To tell what they are, we have to search through the element ID's
1491 * from 0x0001-0x00ff and find one that contains the text
1492 * "PHILIPS IMAGING DD 001" or "Philips Imaging DD 001". The value
1493 * of this element is then used for the upper eight bits of all
1494 * subsequent Philips private element ID's.
1495 */
1496 pms_group = acr_find_group(group_list, PMS_PRIVATE_GROUP_ID);
1497 if (pms_group != NULL) {
1498 pms_element_list = acr_get_group_element_list(pms_group);
1499 creator_id.group_id = PMS_PRIVATE_GROUP_ID;
1500 creator_id.vr_code = ACR_VR_LO;
1501 for (creator_id.element_id = 0x0001;
1502 creator_id.element_id <= 0x00ff;
1503 creator_id.element_id++) {
1504 pms_element = acr_find_element_id(pms_element_list, &creator_id);
1505 if (pms_element != NULL) {
1506 str_ptr = (char *)acr_get_element_string(pms_element);
1507 if (str_ptr != NULL &&
1508 (!strcmp(str_ptr, "Philips Imaging DD 001") ||
1509 !strcmp(str_ptr, "PHILIPS IMAGING DD 001"))) {
1510 /* Found it!!! */
1511 break;
1512 }
1513 }
1514 }
1515 }
1516 else {
1517 creator_id.element_id = 0;
1518 }
1519 if (creator_id.element_id > 0xff || creator_id.element_id < 0x01) {
1520 printf("WARNING: Can't find Philips private creator ID.\n");
1521
1522 /* OK, this may be an old Philips file with the SPI stuff in it.
1523 */
1524 str_ptr = (char *)acr_find_string(group_list, SPI_PMS_grp19_tag, "");
1525 if (!strncmp(str_ptr, "PHILIPS MR", 10)) {
1526 if (G.Debug >= HI_LOGGING) {
1527 printf("Found Philips SPI information\n");
1528 }
1529
1530 str_ptr = (char *)acr_find_string(group_list, ACR_Image_position_patient, "");
1531 if (*str_ptr == '\0') {
1532 double x, y, z;
1533
1534 x = (double)acr_find_double(group_list, SPI_PMS_lr_position2, 0);
1535 y = (double)acr_find_double(group_list, SPI_PMS_ap_position2, 0);
1536 z = (double)acr_find_double(group_list, SPI_PMS_cc_position2, 0);
1537 sprintf(str_buf, "%.15g\\%.15g\\%.15g", x, y, z);
1538 printf("New position: %s\n", str_buf);
1539 acr_insert_string(&group_list, ACR_Image_position_patient,
1540 (Acr_String)str_buf);
1541
1542 }
1543 str_ptr = (char *)acr_find_string(group_list, ACR_Number_of_slices, "");
1544 if (*str_ptr == '\0') {
1545 Acr_Short n;
1546 n = acr_find_short(group_list, SPI_PMS_slice_count, 0);
1547 acr_insert_short(&group_list, ACR_Number_of_slices, n);
1548 }
1549 str_ptr = (char*)acr_find_string(group_list,
1550 ACR_Image_orientation_patient, "");
1551 if (strchr(str_ptr, '\\') == NULL) {
1552 int orientation = acr_find_int(group_list,
1553 SPI_PMS_slice_orientation,
1554 1);
1555 switch (orientation) {
1556 case 1: /* transverse, slice=Z */
1557 strcpy(str_buf, "1\\0\\0\\0\\1\\0");
1558 break;
1559 case 2: /* sagittal, slice=X */
1560 strcpy(str_buf, "0\\1\\0\\0\\0\\1");
1561 break;
1562 case 3: /* coronal, slice=Y */
1563 strcpy(str_buf, "1\\0\\0\\0\\0\\1");
1564 break;
1565 }
1566
1567 acr_insert_string(&group_list, ACR_Image_orientation_patient,
1568 str_buf);
1569 }
1570 #if 0 /* not clear that this is needed */
1571 str_ptr = (char *)acr_find_string(group_list, ACR_Pixel_size, "");
1572 if (*str_ptr == '\0') {
1573 pms_element = acr_find_group_element(group_list,
1574 SPI_PMS_field_of_view);
1575 if (pms_element != NULL) {
1576 double fov[2];
1577 acr_get_element_numeric_array(pms_element, 2, fov);
1578 if (fov[0] != 0.0 && fov[1] != 0.0) {
1579 double derived_spacing[2];
1580 int rows = acr_find_int(group_list, ACR_Rows, 1);
1581 int cols = acr_find_int(group_list, ACR_Columns, 1);
1582
1583 derived_spacing[0] = fov[0] / cols;
1584 derived_spacing[1] = fov[1] / rows;
1585
1586 sprintf(str_buf, "%.15g\\%.15g",
1587 derived_spacing[0], derived_spacing[1]);
1588
1589 printf("PHILIPS: New pixel size %s\n", str_buf);
1590
1591 acr_insert_string(&group_list, ACR_Pixel_size,
1592 str_buf);
1593 }
1594 }
1595 }
1596 #endif
1597 }
1598 }
1599 else {
1600 if (G.Debug >= HI_LOGGING) {
1601 printf("Found Philips private creator ID at %#x\n",
1602 creator_id.element_id);
1603 }
1604
1605 PMS_SET_CREATOR(PMS_Number_of_Slices_MR, &creator_id);
1606 slice_count = acr_find_int(group_list, PMS_Number_of_Slices_MR, -1);
1607 if (slice_count < 0) {
1608 printf("WARNING: Can't find Philips slice count\n");
1609 }
1610 else {
1611 acr_insert_short(&group_list, ACR_Images_in_acquisition,
1612 slice_count);
1613 }
1614 PMS_SET_CREATOR(PMS_Slice_Number_MR, &creator_id);
1615 slice_index = acr_find_int(group_list, PMS_Slice_Number_MR, -1);
1616 if (slice_index < 0) {
1617 printf("WARNING: Can't find Philips slice index\n");
1618 }
1619 else {
1620 /* TODO: It is really quite gross that we have to resort to
1621 * using a Siemens-proprietary field to tell the rest of the
1622 * code which slice we are on. But such is life, for now...
1623 */
1624 acr_insert_numeric(&group_list, SPI_Current_slice_number,
1625 (double) slice_index);
1626 }
1627 }
1628 return (group_list);
1629 }
1630
1631 Acr_Group
add_gems_info(Acr_Group group_list)1632 add_gems_info(Acr_Group group_list)
1633 {
1634 Acr_String tmp_str;
1635 int image_count;
1636 int image_index;
1637 double tr;
1638 int ok;
1639
1640 ok = 1;
1641 tmp_str = acr_find_string(group_list, GEMS_Acqu_private_creator_id, "");
1642 if (strcmp(tmp_str, "GEMS_ACQU_01")) {
1643 ok = 0;
1644 }
1645
1646 tmp_str = acr_find_string(group_list, GEMS_Sers_private_creator_id, "");
1647 if (strcmp(tmp_str, "GEMS_SERS_01")) {
1648 ok = 0;
1649 }
1650
1651 if (!ok) {
1652 /* Warn only for non-PET things. We know the PET scanner doesn't
1653 * rely on this proprietary nonsense.
1654 */
1655 tmp_str = acr_find_string(group_list, ACR_Modality, "");
1656 if (strcmp(tmp_str, "PT") && G.Debug) {
1657 printf("WARNING: GEMS data not found\n");
1658 }
1659 }
1660
1661 tmp_str = acr_find_string(group_list, ACR_Image_type, "");
1662 image_index = acr_find_int(group_list, ACR_Image, -1);
1663 image_count = acr_find_int(group_list, GEMS_Images_in_series, -1);
1664 tr = (double)acr_find_double(group_list, ACR_Repetition_time, 0.0);
1665
1666 /* If we found a valid image count in the proprietary field, copy it
1667 * into the standard field if the standard field is not already set.
1668 */
1669 if (image_count > 0 &&
1670 acr_find_int(group_list, ACR_Images_in_acquisition, -1) < 0) {
1671 acr_insert_long(&group_list, ACR_Images_in_acquisition, image_count);
1672 }
1673
1674 /* Do this only for EPI images for now
1675 */
1676 if (strstr(tmp_str, "\\EPI") != NULL &&
1677 image_index >= 0 &&
1678 image_count > 0 &&
1679 tr != 0.0) {
1680 int frame_count = acr_find_int(group_list,
1681 GEMS_Fast_phases, -1);
1682 if (frame_count > 0) {
1683 if (acr_find_int(group_list, ACR_Acquisitions_in_series, -1) < 0) {
1684 acr_insert_long(&group_list, ACR_Acquisitions_in_series,
1685 (Acr_Long)frame_count);
1686 }
1687 }
1688
1689 if (acr_find_double(group_list, ACR_Frame_reference_time, -1.0) < 0) {
1690 acr_insert_numeric(&group_list, ACR_Frame_reference_time,
1691 (double)(((image_index - 1) / image_count) * tr) );
1692 }
1693 if (acr_find_double(group_list, ACR_Actual_frame_duration, -1.0) < 0) {
1694 acr_insert_numeric(&group_list, ACR_Actual_frame_duration,
1695 tr);
1696 }
1697
1698 if (G.Debug >= HI_LOGGING) {
1699 printf("EPI slice %d timing information: %f, %f\n",
1700 image_index, ((image_index - 1) / image_count) * tr,
1701 tr);
1702 }
1703 }
1704 return (group_list);
1705 }
1706
1707
1708 /* ----------------------------- MNI Header -----------------------------------
1709 @NAME : copy_spi_to_acr()
1710 @INPUT : group_list - read a standard DICOM file
1711 @OUTPUT : (none)
1712 @RETURNS : modified group list
1713 @DESCRIPTION: Routine to copy Siemens-specific (SPI) fields to generic
1714 ACR/NEMA fields.
1715 @METHOD :
1716 @GLOBALS :
1717 @CALLS :
1718 @CREATED : 2005 (Bert Vincent)
1719 @MODIFIED :
1720 ---------------------------------------------------------------------------- */
1721
1722 Acr_Group
copy_spi_to_acr(Acr_Group group_list)1723 copy_spi_to_acr(Acr_Group group_list)
1724 {
1725 int itemp;
1726 double dtemp;
1727
1728 itemp = acr_find_int(group_list, SPI_Number_of_slices_nominal, 0);
1729 if (itemp != 0) {
1730 acr_insert_numeric(&group_list, ACR_Images_in_acquisition, (double)itemp);
1731 }
1732
1733 itemp = acr_find_int(group_list, SPI_Number_of_echoes, 0);
1734 if (itemp != 0) {
1735 acr_insert_numeric(&group_list, ACR_Echo_train_length, (double)itemp);
1736 }
1737
1738 dtemp = (double)acr_find_double(group_list, SPI_Magnetic_field_strength, 0);
1739 if (dtemp != 0.0) {
1740 acr_insert_numeric(&group_list, ACR_Magnetic_field_strength, dtemp);
1741 }
1742 return (group_list);
1743 }
1744
1745 Acr_Group
add_shimadzu_info(Acr_Group group_list)1746 add_shimadzu_info(Acr_Group group_list)
1747 {
1748 Acr_String str_ptr;
1749
1750 /* TODO: Figure out what can be done here...
1751 */
1752 str_ptr = acr_find_string(group_list, ACR_Series_instance_UID, "");
1753 if (*str_ptr != '\0') {
1754 }
1755 return (group_list);
1756 }
1757
1758 /* ----------------------------- MNI Header -----------------------------------
1759 @NAME : read_numa4_dicom
1760 @INPUT : filename - read a standard DICOM file
1761 max_group - maximum group number to read
1762 @OUTPUT : (none)
1763 @RETURNS : group list read in from file
1764 @DESCRIPTION: Routine to read in a group list from a file.
1765 @METHOD :
1766 @GLOBALS :
1767 @CALLS :
1768 @CREATED : December 18, 2001 (Rick Hoge)
1769 @MODIFIED :
1770 ---------------------------------------------------------------------------- */
1771
1772 Acr_Group
read_numa4_dicom(const char * filename,int max_group)1773 read_numa4_dicom(const char *filename, int max_group)
1774 {
1775 Acr_Group group_list;
1776 Acr_String str_ptr;
1777
1778 group_list = read_std_dicom(filename, max_group);
1779 if (group_list == NULL) {
1780 return NULL;
1781 }
1782
1783 /* Check the manufacturer. If it is one we know, try to interpret
1784 * the private/proprietary data if present. This is converted to
1785 * standard DICOM fields whereever it makes sense to do so.
1786 */
1787 str_ptr = acr_find_string(group_list, ACR_Manufacturer, "");
1788 if (strstr(str_ptr, "SIEMENS") != NULL ||
1789 strstr(str_ptr, "Siemens") != NULL) {
1790
1791 group_list = add_siemens_info(group_list);
1792
1793 /* Now copy proprietary fields into the correct places in the
1794 * standard groups.
1795 */
1796 group_list = copy_spi_to_acr(group_list);
1797 }
1798 else if (strstr(str_ptr, "Philips") != NULL) {
1799 group_list = add_philips_info(group_list);
1800 }
1801 else if (strstr(str_ptr, "GEMS") != NULL ||
1802 strstr(str_ptr, "GE MEDICAL") != NULL) {
1803 group_list = add_gems_info(group_list);
1804 }
1805 else if (strstr(str_ptr, "shimadzu") != NULL) {
1806 group_list = add_shimadzu_info(group_list);
1807 }
1808 return (group_list);
1809 }
1810
1811 /* ----------------------------- MNI Header -----------------------------------
1812 @NAME : free_info
1813 @INPUT : gi_ptr
1814 fi_ptr
1815 num_files
1816 @OUTPUT : (none)
1817 @RETURNS : (nothing)
1818 @DESCRIPTION: Routine to free contents of general and file info structures.
1819 @METHOD :
1820 @GLOBALS :
1821 @CALLS :
1822 @CREATED : November 26, 1993 (Peter Neelin)
1823 @MODIFIED :
1824 ---------------------------------------------------------------------------- */
1825
1826 static void
free_info(General_Info * gi_ptr,File_Info * fi_ptr,int num_files)1827 free_info(General_Info *gi_ptr, File_Info *fi_ptr, int num_files)
1828 {
1829 Mri_Index imri;
1830
1831 /* Free the general info pointers */
1832 for (imri = 0; imri < MRI_NDIMS; imri++) {
1833 if (gi_ptr->indices[imri] != NULL) {
1834 free(gi_ptr->indices[imri]);
1835 }
1836 if (gi_ptr->coordinates[imri] != NULL) {
1837 free(gi_ptr->coordinates[imri]);
1838 }
1839 if (gi_ptr->widths[imri] != NULL) {
1840 free(gi_ptr->widths[imri]);
1841 }
1842 }
1843
1844 /* Free the group list */
1845 if (gi_ptr->group_list != NULL) {
1846 acr_delete_group_list(gi_ptr->group_list);
1847 }
1848
1849 return;
1850
1851 }
1852
1853 /* ----------------------------- MNI Header -----------------------------------
1854 @NAME : search_list
1855 @INPUT : value
1856 list_ptr
1857 list_length
1858 start_index - point from which search should start
1859 @OUTPUT : (none)
1860 @RETURNS : Index in list where value is found, or -1 is value not found.
1861 @DESCRIPTION: Routine to search a list for a value, returning the index
1862 into the list. If the value is not found, then -1 is returned.
1863 @METHOD :
1864 @GLOBALS :
1865 @CALLS :
1866 @CREATED : February 28, 1997 (Peter Neelin)
1867 @MODIFIED :
1868 ---------------------------------------------------------------------------- */
1869 int
search_list(int value,const int * list_ptr,int list_length,int start_index)1870 search_list(int value, const int *list_ptr, int list_length, int start_index)
1871 {
1872 int index;
1873
1874 /* If nothing on list, just return. */
1875 if (list_length <= 0) {
1876 return (-1);
1877 }
1878
1879 /* If starting point is invalid, start at zero.
1880 */
1881 if ((start_index >= list_length) || (start_index < 0)) {
1882 start_index = 0;
1883 }
1884
1885 /* Loop over indices, wrapping at the end of the list */
1886 index = start_index;
1887 do {
1888 if (list_ptr[index] == value) {
1889 return index; /* Found it. */
1890 }
1891 index++;
1892 if (index >= list_length) {
1893 index = 0;
1894 }
1895 } while (index != start_index);
1896
1897 return -1; /* Search failed. */
1898 }
1899
1900 /* ----------------------------- MNI Header -----------------------------------
1901 @NAME : sort_dimensions
1902 @INPUT : gi_ptr
1903 @OUTPUT : gi_ptr
1904 @RETURNS : (nothing)
1905 @DESCRIPTION: Routine to sort the MRI dimensions according to their
1906 coordinates. It also fills in the step and start values for
1907 the SLICE dimension.
1908 @METHOD :
1909 @GLOBALS :
1910 @CALLS :
1911 @CREATED : February 28, 1997 (Peter Neelin)
1912 @MODIFIED :
1913 ---------------------------------------------------------------------------- */
1914 static void
sort_dimensions(General_Info * gi_ptr)1915 sort_dimensions(General_Info *gi_ptr)
1916 {
1917 Mri_Index imri;
1918 Sort_Element *sort_array;
1919 int nvalues;
1920 int i,j;
1921 int reverse_array;
1922
1923 /* Sort the dimensions, if needed.
1924 */
1925 for (imri = 0; imri < MRI_NDIMS; imri++) {
1926 /* Sort each dimension iff the size is greater than 1 and it is not
1927 * a time dimension in an N4 mosaic.
1928 */
1929 if (gi_ptr->cur_size[imri] <= 1 ||
1930 /* Don't sort on time for N4 mosaics (TODO: Why not??)
1931 */
1932 ((G.file_type == N4DCM) &&
1933 (imri == TIME) &&
1934 (gi_ptr->num_slices_in_file > 1))) {
1935 if (G.Debug >= HI_LOGGING) {
1936 printf("Not sorting %s dimension\n", Mri_Names[imri]);
1937 }
1938 continue;
1939 }
1940
1941 if (G.Debug >= HI_LOGGING) {
1942 printf("Sorting %s dimension\n", Mri_Names[imri]);
1943 printf(" slice order is:%s\n",gi_ptr->acq.slice_order );
1944 }
1945
1946 /* Set up the array for sorting.
1947 */
1948 nvalues = gi_ptr->cur_size[imri];
1949 sort_array = malloc(nvalues * sizeof(*sort_array));
1950 CHKMEM(sort_array);
1951 for (i = 0; i < nvalues; i++) {
1952 sort_array[i].identifier = gi_ptr->indices[imri][i];
1953 sort_array[i].original_index = i;
1954 sort_array[i].value = gi_ptr->coordinates[imri][i];
1955 sort_array[i].width = gi_ptr->widths[imri][i];
1956 }
1957
1958 /* Sort the array.
1959 */
1960 qsort((void *) sort_array, (size_t) nvalues, sizeof(*sort_array),
1961 dimension_sort_function);
1962
1963 /* Added this from dicomserver conversion, it was probably removed.
1964 * Although it should not matter whether we start negative and
1965 * step positive or the other way around, we should try and stay
1966 * consistent with previous versions and with the dicom header.
1967 * It is especially important for slice times */
1968
1969 /* Figure out if we should reverse the array to keep something
1970 * similar to the original ordering.
1971 * Also need to check for slice ordering, MOSAIC images are
1972 * always sorted from bottom to top whether the sequence was
1973 * ascending or descending */
1974 reverse_array = (sort_array[0].original_index >
1975 sort_array[nvalues-1].original_index) ||
1976 (!strcmp(gi_ptr->acq.slice_order,"descending")&& !strcmp(Mri_Names[imri], "Slice"));
1977
1978 /* Copy the information back into the appropriate arrays */
1979 for (i=0; i < nvalues; i++) {
1980 j = (reverse_array ? nvalues - i - 1 : i);
1981 gi_ptr->indices[imri][i] = sort_array[j].identifier;
1982 gi_ptr->coordinates[imri][i] = sort_array[j].value;
1983 gi_ptr->widths[imri][i] = sort_array[j].width;
1984 }
1985
1986 if (G.Debug >= HI_LOGGING) {
1987 /* Print out all of the information about this dimension.
1988 */
1989 if(reverse_array){
1990 printf(" nvalues %d min %f max %f\n",
1991 nvalues,
1992 gi_ptr->coordinates[imri][nvalues - 1],
1993 gi_ptr->coordinates[imri][0]);
1994 }else{
1995 printf(" nvalues %d min %f max %f\n",
1996 nvalues,
1997 gi_ptr->coordinates[imri][0],
1998 gi_ptr->coordinates[imri][nvalues - 1]);
1999 }
2000 for (i = 0; i < nvalues; i++) {
2001 printf("%3d %6d %8.3f\n",
2002 i,
2003 gi_ptr->indices[imri][i],
2004 gi_ptr->coordinates[imri][i]);
2005 }
2006 }
2007
2008 /* Free the temporary array we used to sort the coordinate axis.
2009 */
2010 free(sort_array);
2011
2012 /* Now verify that the slice coordinate array looks sane. We
2013 * don't complain about things like missing time slices, since
2014 * many PET scanners produce irregular timings.
2015 */
2016 if (nvalues >= 2 && imri == SLICE) {
2017 /* Calculate the spacing between the first and second slice.
2018 */
2019 double delta = (gi_ptr->coordinates[imri][1] -
2020 gi_ptr->coordinates[imri][0]);
2021
2022 for (i = 1; i < nvalues; i++) {
2023 /* Check that each successive slice has roughly the same
2024 * spacing, to within 2 percent of the initial delta.
2025 */
2026
2027 /* Check against absolute value of epsilon, could be positive
2028 * or negative, because of reversal (fcmp does not handle this)
2029 */
2030 if (!fcmp(delta,
2031 (gi_ptr->coordinates[imri][i] -
2032 gi_ptr->coordinates[imri][i - 1]),
2033 fabs(2.0 * (delta / 100.0)))) {
2034 /* TODO: Perhaps this message could be improved?? */
2035 printf("WARNING: Missing %s data at index %d, %g %g\n",
2036 Mri_Names[imri],
2037 i,
2038 delta,
2039 (gi_ptr->coordinates[imri][i] -
2040 gi_ptr->coordinates[imri][i - 1]));
2041 }
2042 }
2043 }
2044
2045 /* Update slice step and start.
2046 */
2047 if (imri == SLICE) {
2048 if (gi_ptr->coordinates[imri][0] !=
2049 gi_ptr->coordinates[imri][nvalues-1]) {
2050 double dbl_tmp1;
2051 double dbl_tmp2;
2052
2053 dbl_tmp1 = gi_ptr->step[gi_ptr->slice_world];
2054 dbl_tmp2 = (gi_ptr->coordinates[imri][nvalues - 1] -
2055 gi_ptr->coordinates[imri][0]) /
2056 ((double) gi_ptr->cur_size[imri] - 1.0);
2057
2058 /* OK, so now we have to figure out who wins the great
2059 * battle to become the slice step value. If the value
2060 * we have been assuming up until now is the same as the
2061 * calculated value with only a sign change, or if the
2062 * assumed value is the default (1.0), we adopt the
2063 * calculated value.
2064 */
2065 /* Check against absolute value of epsilon, could be positive
2066 * or negative, because of reversal (fcmp does not handle this)
2067 */
2068 if (!fcmp(dbl_tmp1, dbl_tmp2, fabs((dbl_tmp1 / 1000.0)))) {
2069 printf("WARNING: Coordinate spacing (%g) differs from DICOM slice spacing (%g)\n", dbl_tmp2, dbl_tmp1);
2070 if (!G.prefer_coords) {
2071 printf(" (perhaps you should consider the -usecoordinates option)\n");
2072 }
2073 if (dbl_tmp1 == 1.0 || G.prefer_coords || fcmp(dbl_tmp1, -dbl_tmp2, fabs(dbl_tmp1 / 1000.0))) {
2074 /*Although in the comment above it says that we should use the
2075 calculated value if the assumed and calculated value only differ
2076 by a sign change, this was not included in this if statement, added
2077 by ilana because default was wrong (i.e. w/o using -usecoordinates option)*/
2078 gi_ptr->step[gi_ptr->slice_world] = dbl_tmp2;
2079 }
2080 printf(" Using %g for the slice spacing value.\n",
2081 gi_ptr->step[gi_ptr->slice_world]);
2082 }
2083
2084 }
2085 gi_ptr->start[gi_ptr->slice_world] = gi_ptr->coordinates[imri][0];
2086 if (G.Debug >= HI_LOGGING) {
2087 printf("Set slice %d step to %f, start to %f\n",
2088 gi_ptr->slice_world,
2089 gi_ptr->step[gi_ptr->slice_world],
2090 gi_ptr->start[gi_ptr->slice_world]);
2091 }
2092 } /* If slice dimension */
2093 } /* for all mri dimensions */
2094 }
2095
2096 /* ----------------------------- MNI Header -----------------------------------
2097 @NAME : dimension_sort_function
2098 @INPUT : v1, v2 - values to compare
2099 @OUTPUT : (none)
2100 @RETURNS : -1, 0 or 1 if v1 < v2, v1 == v2 or v1 > v2
2101 @DESCRIPTION: Function to compare to array elements for sorting. Elements are
2102 compared first on value, then on their original array index
2103 (this tries to preserve the original sequence).
2104 @METHOD :
2105 @GLOBALS :
2106 @CALLS :
2107 @CREATED : February 28, 1997 (Peter Neelin)
2108 @MODIFIED :
2109 ---------------------------------------------------------------------------- */
2110 static int
dimension_sort_function(const void * v1,const void * v2)2111 dimension_sort_function(const void *v1, const void *v2)
2112 {
2113 Sort_Element *value1, *value2;
2114
2115 value1 = (Sort_Element *) v1;
2116 value2 = (Sort_Element *) v2;
2117
2118 if (value1->value < value2->value)
2119 return -1;
2120 else if (value1->value > value2->value)
2121 return 1;
2122 else if (value1->original_index < value2->original_index)
2123 return -1;
2124 else if (value1->original_index > value2->original_index)
2125 return 1;
2126 else
2127 return 0;
2128 }
2129
2130 int
prot_find_string(Acr_Element elem_ptr,const char * name_str,char * value)2131 prot_find_string(Acr_Element elem_ptr, const char *name_str, char *value)
2132 {
2133 static const char prot_head[] = "### ASCCONV BEGIN ###";
2134 long cur_offset,tmp_offset;
2135 long max_offset;
2136 char *field_ptr;
2137 int ix1, ix2;
2138
2139 max_offset = elem_ptr->data_length - sizeof(prot_head);
2140 tmp_offset = max_offset;
2141
2142 // Scan through the element containing the protocol, to find the
2143 // ASCII dump of the MrProt structure.
2144 //
2145 /*For some reason, some dicom files have 2 "ASCCONV BEGIN" tags, this screws up the search.
2146 Keep the second one we find for now, don't know if this will always work for these dual-ASCCONV
2147 files. IRL*/
2148
2149 for (cur_offset = 0; cur_offset < max_offset; cur_offset++) {
2150 if (!memcmp(elem_ptr->data_pointer + cur_offset,
2151 prot_head, sizeof(prot_head) - 1)) {
2152 tmp_offset = cur_offset;
2153 }
2154 }
2155
2156 cur_offset = tmp_offset; /*set it to the last occurence of "ASCCONV BEGIN"*/
2157
2158 /* bail if we didn't find the protocol
2159 */
2160 if (cur_offset == max_offset) {
2161 return (0);
2162 }
2163
2164 field_ptr = strstr(elem_ptr->data_pointer + cur_offset, name_str);
2165 if (field_ptr != NULL) {
2166 sscanf(field_ptr, "%*s %*s %s", value);
2167 ix1 = 0;
2168 for (ix2 = 0; value[ix2] != '\0'; ix2++) {
2169 if (value[ix2] != '"') {
2170 value[ix1++] = value[ix2];
2171 }
2172 }
2173 }
2174 else {
2175 strcpy(value, "0");
2176 }
2177 return (1);
2178 }
2179
2180 static char *
dump_protocol_text(Acr_Element elem_ptr)2181 dump_protocol_text(Acr_Element elem_ptr)
2182 {
2183 const char prot_head[] = "### ASCCONV BEGIN ###";
2184 const char prot_tail[] = "### ASCCONV END ###";
2185 char *output = malloc(elem_ptr->data_length);
2186 int prot_found = FALSE;
2187 long cur_offset, tmp_offset;
2188 long max_offset;
2189
2190 CHKMEM(output);
2191
2192 // scan throught the group containing the protocol, to find the
2193 // ascii dump of the MrProt structure
2194 max_offset = elem_ptr->data_length - sizeof (prot_head);
2195
2196 /*For some reason, some dicom files have 2 "ASCCONV BEGIN" tags, this screws up the search.
2197 Keep the second one we find for now, don't know if this will always work for these dual-ASCCONV
2198 files. IRL*/
2199
2200 for (cur_offset = 0; cur_offset < max_offset; cur_offset++) {
2201 if (!memcmp(elem_ptr->data_pointer + cur_offset,
2202 prot_head, sizeof(prot_head) - 1)) {
2203 prot_found = TRUE;
2204 tmp_offset = cur_offset;
2205 }
2206 }
2207
2208 cur_offset = tmp_offset; /*set it to the last occurence of "ASCCONV BEGIN"*/
2209
2210 if (prot_found) {
2211 int ix1 = 0;
2212 char *tmp_ptr = elem_ptr->data_pointer + cur_offset;
2213
2214 for ( ; cur_offset < max_offset; cur_offset++) {
2215 if (!memcmp(tmp_ptr, prot_tail, sizeof(prot_tail) - 1)) {
2216 break;
2217 }
2218 if (*tmp_ptr != '"') {
2219 output[ix1++] = *tmp_ptr;
2220 }
2221 tmp_ptr++;
2222 }
2223 output[ix1] = '\0';
2224 }
2225 return (output);
2226 }
2227
2228 static void
copy_element_properties(Acr_Element new_element,Acr_Element old_element)2229 copy_element_properties(Acr_Element new_element, Acr_Element old_element)
2230 {
2231 acr_set_element_byte_order(new_element,
2232 acr_get_element_byte_order(old_element));
2233
2234 acr_set_element_vr_encoding(new_element,
2235 acr_get_element_vr_encoding(old_element));
2236 }
2237
2238 static int
mosaic_init(Acr_Group group_list,Mosaic_Info * mi_ptr,int load_image)2239 mosaic_init(Acr_Group group_list, Mosaic_Info *mi_ptr, int load_image)
2240 {
2241 int group_id, element_id;
2242 int grid_size;
2243 long new_image_size;
2244 void *data;
2245 Acr_Element element;
2246 int i;
2247 double pixel_spacing[2];
2248 Acr_Double separation;
2249 double RowColVec[6];
2250 double dircos[VOL_NDIMS][WORLD_NDIMS];
2251 Acr_String str_tmp,str_tmp2;
2252 int old = 1;
2253
2254 if (G.Debug >= HI_LOGGING) {
2255 printf("mosaic_init(%lx, %lx, %d)\n",
2256 (unsigned long) group_list, (unsigned long) mi_ptr,
2257 load_image);
2258 }
2259
2260 str_tmp = acr_find_string(group_list, SPI_Order_of_slices, "");
2261 /* Seems like this field is often not found,
2262 * the sSliceArray.ucMode field found in the ASCONV header
2263 * seems trustworthy. */
2264 /* 0x1 means ASCENDING
2265 * 0x2 means DESCENDING
2266 * 0x4 means INTERLEAVED*/
2267 str_tmp2 = acr_find_string(group_list, EXT_Slice_order, "");
2268
2269 if (!strncmp(str_tmp, "INTERLEAVED", 11) || !strncmp(str_tmp2, "0x4", 3)) {
2270 mi_ptr->mosaic_seq = MOSAIC_SEQ_INTERLEAVED;
2271 str_tmp="interleaved";
2272
2273 }
2274 else if (!strncmp(str_tmp, "ASCENDING", 9)|| !strncmp(str_tmp2, "0x1", 3)) {
2275 mi_ptr->mosaic_seq = MOSAIC_SEQ_ASCENDING;
2276 str_tmp="ascending";
2277 }
2278 else if (!strncmp(str_tmp, "DESCENDING", 10)|| !strncmp(str_tmp2, "0x2", 3)) {
2279 mi_ptr->mosaic_seq = MOSAIC_SEQ_DESCENDING;
2280 str_tmp="descending";
2281 }
2282 else {
2283 mi_ptr->mosaic_seq = G.mosaic_seq;
2284 str_tmp="default ascending";
2285 }
2286
2287 if (G.Debug >= HI_LOGGING) {
2288 printf(" ordering is %s\n", str_tmp);
2289 }
2290
2291 /* Get some basic image information.
2292 * big[0/1] is number of columns/rows in whole mosaic.
2293 */
2294 mi_ptr->big[0] = acr_find_int(group_list, ACR_Columns, 1);
2295 mi_ptr->big[1] = acr_find_int(group_list, ACR_Rows, 1);
2296 mi_ptr->pixel_size =
2297 (acr_find_int(group_list, ACR_Bits_allocated, 16) - 1) / 8 + 1;
2298
2299 /* Get the image size
2300 * (size[0/1] is number of columns/rows in a single slice)
2301 */
2302 mi_ptr->size[0] = (int)acr_find_short(group_list, EXT_Sub_image_columns, 1);
2303 mi_ptr->size[1] = (int)acr_find_short(group_list, EXT_Sub_image_rows, 1);
2304
2305 // Get the grid shape, checking that it is not too big if specified
2306 mi_ptr->grid[0] = mi_ptr->big[0] / mi_ptr->size[0];
2307 mi_ptr->grid[1] = mi_ptr->big[1] / mi_ptr->size[1];
2308 if ((mi_ptr->grid[0] < 1) || (mi_ptr->grid[0] < 1)) {
2309 fprintf(stderr, "Grid too small: %d x %d, size %d x %d, big %d x %d\n",
2310 mi_ptr->grid[0], mi_ptr->grid[1],
2311 mi_ptr->size[0], mi_ptr->size[1],
2312 mi_ptr->big[0], mi_ptr->big[1]);
2313 exit(EXIT_FAILURE);
2314 }
2315
2316 // Check whether we need to do anything (1x1 grid may be the whole image)
2317 mi_ptr->big_image = NULL;
2318 mi_ptr->small_image = NULL;
2319 mi_ptr->packed = FALSE;
2320 grid_size = mi_ptr->grid[0] * mi_ptr->grid[1];
2321 if ((grid_size == 1) &&
2322 (mi_ptr->size[0] == mi_ptr->big[0]) &&
2323 (mi_ptr->size[1] == mi_ptr->big[1])) {
2324 /* had to remove this as now ANY images acquired with
2325 the mosaic sequence need special treatment */
2326 mi_ptr->packed = FALSE;
2327 return 1;
2328 }
2329
2330 /* Update the number of image rows and columns
2331 */
2332 acr_insert_short(&group_list, ACR_Rows, mi_ptr->size[1]);
2333 acr_insert_short(&group_list, ACR_Columns, mi_ptr->size[0]);
2334
2335 /* Get image image index info (number of slices in file)
2336 */
2337 mi_ptr->slice_count = acr_find_int(group_list, EXT_Slices_in_file, 1);
2338
2339 /* sub_images is now just the number of mosaic elements, even if
2340 * they don't all contain slices
2341 */
2342 mi_ptr->sub_images = mi_ptr->grid[0] * mi_ptr->grid[1];
2343
2344
2345 /* Get the pixel size.
2346 */
2347 dicom_read_pixel_size(group_list, pixel_spacing);
2348
2349 /* Get step between slices
2350 */
2351 separation = acr_find_double(group_list, ACR_Slice_thickness, 0.0);
2352 if (separation == 0.0) {
2353 separation = acr_find_double(group_list, ACR_Spacing_between_slices,
2354 1.0);
2355 }
2356
2357 /* get image normal vector
2358 * (need to compute based on dicom field, which gives
2359 * unit vectors for row and column direction)
2360 * TODO: This code (and other code in this function) could probably
2361 * be broken out and made redundant with other code in the converter.
2362 */
2363 if (dicom_read_orientation(group_list, RowColVec)) {
2364 memcpy(dircos[VCOLUMN], RowColVec, sizeof(*RowColVec) * WORLD_NDIMS);
2365 memcpy(dircos[VROW], &RowColVec[3], sizeof(*RowColVec) * WORLD_NDIMS);
2366
2367 /* compute slice normal as cross product of row/column unit vectors
2368 * (should check for unit length?)
2369 */
2370 mi_ptr->normal[XCOORD] =
2371 dircos[VCOLUMN][YCOORD] * dircos[VROW][ZCOORD] -
2372 dircos[VCOLUMN][ZCOORD] * dircos[VROW][YCOORD];
2373
2374 mi_ptr->normal[YCOORD] =
2375 dircos[VCOLUMN][ZCOORD] * dircos[VROW][XCOORD] -
2376 dircos[VCOLUMN][XCOORD] * dircos[VROW][ZCOORD];
2377
2378 mi_ptr->normal[ZCOORD] =
2379 dircos[VCOLUMN][XCOORD] * dircos[VROW][YCOORD] -
2380 dircos[VCOLUMN][YCOORD] * dircos[VROW][XCOORD];
2381 }
2382
2383 /* compute slice-to-slice step vector
2384 */
2385 for (i = 0; i < WORLD_NDIMS; i++) {
2386 mi_ptr->step[i] = separation * mi_ptr->normal[i];
2387 }
2388
2389 /* Get position and correct to first slice
2390 */
2391 if (!dicom_read_position(group_list, 0, mi_ptr->position)) {
2392 if (G.Debug) {
2393 printf("WARNING: No image position found\n");
2394 }
2395 mi_ptr->position[XCOORD] = mi_ptr->position[YCOORD] =
2396 mi_ptr->position[ZCOORD] = 0.0;
2397 }
2398
2399 if (G.Debug >= HI_LOGGING) {
2400 printf(" step %.3f %.3f %.3f pos %.3f %.3f %.3f normal %.3f,%.3f,%.3f\n",
2401 mi_ptr->step[0],
2402 mi_ptr->step[1],
2403 mi_ptr->step[2],
2404 mi_ptr->position[0],
2405 mi_ptr->position[1],
2406 mi_ptr->position[2],
2407 mi_ptr->normal[0],
2408 mi_ptr->normal[1],
2409 mi_ptr->normal[2]
2410 );
2411 }
2412
2413 /* Treat slice ordering differently depending on software version
2414 - for newer software versions(>VA25 or >VB), the mosaic images seem to
2415 always be ordered ascending, regardless of acquisition order.*/
2416
2417 str_tmp=strstr(acr_find_string(group_list, ACR_Software_versions, ""),
2418 "syngo MR A");
2419 str_tmp2=strstr(acr_find_string(group_list, ACR_Software_versions, ""),
2420 "syngo MR B");
2421
2422 if(str_tmp !=NULL){
2423 if(atoi(str_tmp + 10) >= 25) old=0; /*software version >= VA25*/
2424 }
2425 else if(str_tmp2 !=NULL){
2426 if(atoi(str_tmp2 + 10) >= 11) old=0; /*software version >= VB11*/
2427 }
2428
2429 if(old==1){ /*old behavior*/
2430 if (mi_ptr->mosaic_seq != MOSAIC_SEQ_INTERLEAVED) {
2431 if (is_numaris3(group_list)) {
2432 for (i = 0; i < WORLD_NDIMS; i++) {
2433 mi_ptr->position[i] -=
2434 (double) (mi_ptr->sub_images-1) * mi_ptr->step[i];
2435 }
2436 }
2437 else {
2438 /* Numaris 4 mosaic correction:
2439 * - position given is edge of huge slice constructed as if
2440 * real slice was at center of mosaic
2441 * - mi_ptr->big[0,1] are number of columns and rows of mosaic
2442 * - mi_ptr->size[0,1] are number of columns and rows of sub-image
2443 */
2444
2445 for (i = 0; i < WORLD_NDIMS; i++) {
2446
2447 /* Correct offset from mosaic Center
2448 */
2449 mi_ptr->position[i] += (double)
2450 ((dircos[VCOLUMN][i] * mi_ptr->big[0] * pixel_spacing[0]/2.0) +
2451 (dircos[VROW][i] * mi_ptr->big[1] * pixel_spacing[1]/2));
2452
2453 /* Move from center to corner of slice
2454 */
2455 mi_ptr->position[i] -=
2456 ((dircos[VCOLUMN][i] * mi_ptr->size[0] * pixel_spacing[0]/2.0) +
2457 (dircos[VROW][i] * mi_ptr->size[1] * pixel_spacing[1]/2.0));
2458 }
2459
2460 }
2461 }
2462 }
2463 else{ /*new behavior*/
2464 for (i = 0; i < WORLD_NDIMS; i++) {
2465
2466 /* Correct offset from mosaic Center
2467 */
2468 mi_ptr->position[i] += (double)
2469 ((dircos[VCOLUMN][i] * mi_ptr->big[0] * pixel_spacing[0]/2.0) +
2470 (dircos[VROW][i] * mi_ptr->big[1] * pixel_spacing[1]/2));
2471
2472 /* Move from center to corner of slice
2473 */
2474 mi_ptr->position[i] -=
2475 ((dircos[VCOLUMN][i] * mi_ptr->size[0] * pixel_spacing[0]/2.0) +
2476 (dircos[VROW][i] * mi_ptr->size[1] * pixel_spacing[1]/2.0));
2477 }
2478 }
2479
2480 if (G.Debug >= HI_LOGGING) {
2481 printf(" corrected position %.3f %.3f %.3f\n",
2482 mi_ptr->position[0],
2483 mi_ptr->position[1],
2484 mi_ptr->position[2]);
2485 }
2486
2487 /* Now that we've corrected the _position_, we still might need to correct
2488 * the pixel spacing. For some gadawful reason, some (but not all) Siemens
2489 * mosaic files store the pixel spacing such that it is scaled by the ratio
2490 * between the mosaic image size and the subimage size. The best way I can
2491 * think of to detect this is by comparing the field of view to the product
2492 * of the pixel spacing and the actual subimage size, and if they are
2493 * wildly different, assume that we need to perform the scaling...
2494 */
2495 element = acr_find_group_element(group_list, SPI_Field_of_view);
2496 if (element != NULL) {
2497 double fov[2]; /* field of view (row/column) */
2498
2499 acr_get_element_numeric_array(element, 2, fov);
2500 if (fov[0] != 0.0 && fov[1] != 0.0) {
2501 double derived_spacing[2];
2502 int ratio[2];
2503 char str_buf[128];
2504
2505 /* Calculate the actual spacing required to cover the field of view.
2506 */
2507 derived_spacing[0] = fov[0] / mi_ptr->size[0];
2508 derived_spacing[1] = fov[1] / mi_ptr->size[1];
2509
2510 ratio[0] = (int) rint(derived_spacing[0] / pixel_spacing[0]);
2511 ratio[1] = (int) rint(derived_spacing[1] / pixel_spacing[1]);
2512
2513 /* If the ratio of the derived spacing to the pixel spacing is
2514 * the same as the ratio between the edge length of the large
2515 * and small images, adopt the derived spacing as the correct
2516 * value.
2517 */
2518 if (ratio[0] == (mi_ptr->big[0] / mi_ptr->size[0]) ||
2519 ratio[1] == (mi_ptr->big[1] / mi_ptr->size[1])) {
2520
2521 if (G.Debug) {
2522 printf("Updating mosaic pixel spacing from %f,%f to %f,%f\n",
2523 pixel_spacing[0],pixel_spacing[1],
2524 derived_spacing[0],derived_spacing[1]);
2525 }
2526
2527 /* Store the updated pixel spacing in the group list.
2528 */
2529 sprintf(str_buf, "%.15g\\%.15g",
2530 derived_spacing[0], derived_spacing[1]);
2531 acr_insert_string(&group_list, ACR_Pixel_size, str_buf);
2532 }
2533 }
2534 }
2535
2536 if (!load_image) {
2537 mi_ptr->big_image = NULL;
2538 mi_ptr->small_image = NULL;
2539 }
2540 else {
2541
2542 /* Steal the image element from the group list
2543 */
2544 mi_ptr->packed = TRUE;
2545 mi_ptr->big_image = acr_find_group_element(group_list, ACR_Pixel_data);
2546 if (mi_ptr->big_image == NULL) {
2547 fprintf(stderr, "Couldn't find an image\n");
2548 exit(EXIT_FAILURE);
2549 }
2550 group_id = acr_get_element_group(mi_ptr->big_image);
2551 element_id = acr_get_element_element(mi_ptr->big_image);
2552 acr_group_steal_element(acr_find_group(group_list, group_id),
2553 mi_ptr->big_image);
2554
2555 /* Add a small image
2556 */
2557 new_image_size = mi_ptr->size[0] * mi_ptr->size[1] * mi_ptr->pixel_size;
2558 data = malloc(new_image_size);
2559 CHKMEM(data);
2560 mi_ptr->small_image = acr_create_element(group_id, element_id,
2561 acr_get_element_vr(mi_ptr->big_image),
2562 new_image_size, data);
2563 acr_set_element_vr(mi_ptr->small_image,
2564 acr_get_element_vr(mi_ptr->big_image));
2565
2566 copy_element_properties(mi_ptr->small_image, mi_ptr->big_image);
2567
2568 acr_insert_element_into_group_list(&group_list, mi_ptr->small_image);
2569 }
2570
2571 /* Return number of sub-images in this image */
2572 return mi_ptr->sub_images;
2573 }
2574
2575 static int
mosaic_insert_subframe(Acr_Group group_list,Mosaic_Info * mi_ptr,int iimage,int load_image)2576 mosaic_insert_subframe(Acr_Group group_list, Mosaic_Info *mi_ptr,
2577 int iimage, int load_image)
2578 {
2579 int irow;
2580 int idim;
2581 int nbyte;
2582 int isub;
2583 int jsub;
2584 char *new;
2585 char *old;
2586 long old_offset;
2587 long new_offset;
2588 double position[WORLD_NDIMS];
2589 string_t string;
2590 int islice;
2591 char *str_tmp, *str_tmp2;
2592 int oldb=1;
2593
2594 if (G.Debug >= HI_LOGGING) {
2595 printf("mosaic_insert_subframe(%lx, %lx, %d, %d)\n",
2596 (unsigned long)group_list, (unsigned long)mi_ptr,
2597 iimage, load_image);
2598 }
2599
2600 /* Figure out what to do based upon the mosaic sequencing.
2601 */
2602 /* Since at least software version VA25 (and thus VA30, VB15),
2603 * the mosaic sequencing in the file is the same, regardless
2604 * of the acquisition (always ascending).
2605 * Also, the following code is based on a field that is
2606 * deprecated (0x0021 0x123f)... still keep old behavior in case
2607 * its an older image type*/
2608
2609 str_tmp=strstr(acr_find_string(group_list, ACR_Software_versions, ""),
2610 "syngo MR A");
2611 str_tmp2=strstr(acr_find_string(group_list, ACR_Software_versions, ""),
2612 "syngo MR B");
2613
2614 if(str_tmp !=NULL){
2615 if(atoi(str_tmp + 10) >= 25) oldb=0; /*software version >= VA25*/
2616 }
2617 else if(str_tmp2 !=NULL){
2618 if(atoi(str_tmp2 + 10) >= 11) oldb=0; /*software version >= VB11*/
2619 }
2620
2621 if(mi_ptr->mosaic_seq == MOSAIC_SEQ_INTERLEAVED && oldb==1){ //old behavior
2622 /*case MOSAIC_SEQ_INTERLEAVED:*/
2623 /* For interleaved sequences, we have to map the odd slices to
2624 * the range slice_count/2..slice_count-1 and the even slices
2625 * from zero to slice_count/2-1
2626 */
2627 if (iimage & 1) { /* Odd?? */
2628 islice = (mi_ptr->slice_count / 2) + (iimage / 2);
2629 }
2630 else {
2631 islice = iimage / 2;
2632 }
2633
2634 //break;
2635 }
2636 else{
2637 /*default:*/
2638
2639 /* Otherwise, just use the image number without modification for
2640 * ascending or unknown slice ordering.
2641 */
2642 islice = iimage;
2643 }
2644
2645 #if 0
2646 if (is_numaris3(group_list)) {
2647 islice = mi_ptr->slice_count - islice - 1;
2648 }
2649 #endif
2650
2651 /* Check the image number
2652 */
2653 if ((iimage < 0) || (iimage > mi_ptr->sub_images)) {
2654 fprintf(stderr, "Invalid image number to send: %d of %d\n",
2655 iimage, mi_ptr->sub_images);
2656 exit(EXIT_FAILURE);
2657 }
2658
2659 /* Update the index
2660 */
2661 acr_insert_numeric(&group_list, SPI_Current_slice_number, (double) iimage);
2662
2663 /* Update the position
2664 */
2665 for (idim = 0; idim < WORLD_NDIMS; idim++) {
2666 position[idim] = mi_ptr->position[idim] +
2667 (double) iimage * mi_ptr->step[idim];
2668 }
2669
2670 /* If the sequence is descending, invert the slice coordinate.
2671 * This involves subtracting the step * index from the mosaic
2672 * slice position.
2673 */
2674 /*This should already have been taken care of in sort_dimensions above*/
2675 /*if (mi_ptr->mosaic_seq == MOSAIC_SEQ_DESCENDING) {
2676 position[ZCOORD] = mi_ptr->position[ZCOORD] -
2677 (double) islice * mi_ptr->step[ZCOORD];
2678 }*/
2679
2680 sprintf(string, "%.15g\\%.15g\\%.15g",
2681 position[XCOORD], position[YCOORD], position[ZCOORD]);
2682
2683 acr_insert_string(&group_list, ACR_Image_position_patient, string);
2684 acr_insert_string(&group_list, SPI_Image_position, string);
2685
2686 /* HMM - is this necessary?? */
2687 if (is_numaris3(group_list)) {
2688 update_coordinate_info(group_list);
2689 }
2690
2691 if (G.Debug >= HI_LOGGING) {
2692 printf(" position %s\n", string);
2693 }
2694
2695 if (load_image) {
2696 /* Figure out the sub-image indices
2697 */
2698 isub = islice % mi_ptr->grid[0];
2699 jsub = islice / mi_ptr->grid[0];
2700
2701 /* Get pointers
2702 */
2703 old = acr_get_element_data(mi_ptr->big_image);
2704 new = acr_get_element_data(mi_ptr->small_image);
2705
2706 /* Copy the image
2707 */
2708 nbyte = mi_ptr->size[0] * mi_ptr->pixel_size;
2709 for (irow = 0; irow < mi_ptr->size[1]; irow++) {
2710 old_offset = isub * mi_ptr->size[0] +
2711 (jsub * mi_ptr->size[1] + irow) * mi_ptr->big[0];
2712 old_offset *= mi_ptr->pixel_size;
2713 new_offset = (irow * mi_ptr->size[0]) * mi_ptr->pixel_size;
2714 memcpy(&new[new_offset], &old[old_offset], nbyte);
2715 }
2716
2717 /* Reset the byte order and VR encoding. This will be modified on each
2718 * send according to what the connection needs.
2719 */
2720 copy_element_properties(mi_ptr->small_image, mi_ptr->big_image);
2721 }
2722 return 1;
2723
2724 }
2725
2726 static void
mosaic_cleanup(Mosaic_Info * mi_ptr)2727 mosaic_cleanup(Mosaic_Info *mi_ptr)
2728 {
2729 if (mi_ptr) {
2730 if (mi_ptr->packed && mi_ptr->big_image != NULL) {
2731 acr_delete_element(mi_ptr->big_image);
2732 mi_ptr->packed = FALSE;
2733 mi_ptr->big_image = NULL;
2734 }
2735 }
2736 }
2737
2738 /************************************************************************
2739 * Multiframe functions, which mimic the behavior of the older mosaic
2740 * functions. This is not yet a complete or correct implementation of
2741 * multiframe DICOM - see the NOTE: above!
2742 */
2743 /* ----------------------------- MNI Header -----------------------------------
2744 @NAME : multiframe_init()
2745 @INPUT : group_list - the list of DICOM groups/elements that make up
2746 this file.
2747 int load_image - a boolean value, non-zero if the function
2748 should actually load the data into memory.
2749 @OUTPUT : mfi_ptr - a pointer to a Multiframe_Info structure that will
2750 contain information used to expand this multiframe file into
2751 a series of slices.
2752 @RETURNS : void
2753 @DESCRIPTION: Initialize a multiframe conversion process.
2754 @METHOD :
2755 @GLOBALS :
2756 @CALLS :
2757 @CREATED : June 3, 2005 Bert Vincent
2758 @MODIFIED :
2759 ---------------------------------------------------------------------------- */
2760
2761 static void
multiframe_init(Acr_Group group_list,Multiframe_Info * mfi_ptr,int load_image)2762 multiframe_init(Acr_Group group_list, Multiframe_Info *mfi_ptr, int load_image)
2763 {
2764 int grp_id;
2765 int elm_id;
2766 void *data_ptr;
2767 Acr_Element element;
2768 int i;
2769 Acr_Double spacing;
2770 double RowColVec[6];
2771 double dircos[VOL_NDIMS][WORLD_NDIMS];
2772 int rows;
2773 int cols;
2774 int pixel_size;
2775
2776 if (G.Debug >= HI_LOGGING) {
2777 printf("multiframe_init(%lx, %lx, %d)\n",
2778 (unsigned long) group_list, (unsigned long) mfi_ptr,
2779 load_image);
2780 }
2781
2782 cols = acr_find_int(group_list, ACR_Columns, 1);
2783 rows = acr_find_int(group_list, ACR_Rows, 1);
2784 pixel_size =
2785 (acr_find_int(group_list, ACR_Bits_allocated, 16) - 1) / 8 + 1;
2786
2787 /* Get image image index info (number of slices in file)
2788 */
2789 mfi_ptr->frame_count = acr_find_int(group_list, ACR_Number_of_frames, 1);
2790
2791 /* Get spacing between slices
2792 */
2793 spacing = acr_find_double(group_list, ACR_Slice_thickness, 0.0);
2794 if (spacing == 0.0) {
2795 spacing = acr_find_double(group_list, ACR_Spacing_between_slices, 1.0);
2796 }
2797
2798 /* get image normal vector
2799 * (need to compute based on dicom field, which gives
2800 * unit vectors for row and column direction)
2801 * TODO: This code (and other code in this function) could probably
2802 * be broken out and made redundant with other code in the converter.
2803 */
2804 if (dicom_read_orientation(group_list, RowColVec)) {
2805 memcpy(dircos[VCOLUMN], RowColVec, sizeof(*RowColVec) * WORLD_NDIMS);
2806 memcpy(dircos[VROW], &RowColVec[3], sizeof(*RowColVec) * WORLD_NDIMS);
2807
2808 convert_dicom_coordinate(dircos[VROW]);
2809 convert_dicom_coordinate(dircos[VCOLUMN]);
2810
2811 /* compute slice normal as cross product of row/column unit vectors
2812 * (should check for unit length?)
2813 */
2814 mfi_ptr->normal[XCOORD] =
2815 dircos[VCOLUMN][YCOORD] * dircos[VROW][ZCOORD] -
2816 dircos[VCOLUMN][ZCOORD] * dircos[VROW][YCOORD];
2817
2818 mfi_ptr->normal[YCOORD] =
2819 dircos[VCOLUMN][ZCOORD] * dircos[VROW][XCOORD] -
2820 dircos[VCOLUMN][XCOORD] * dircos[VROW][ZCOORD];
2821
2822 mfi_ptr->normal[ZCOORD] =
2823 dircos[VCOLUMN][XCOORD] * dircos[VROW][YCOORD] -
2824 dircos[VCOLUMN][YCOORD] * dircos[VROW][XCOORD];
2825 }
2826
2827 /* If the normal is unreliable, use a default value.
2828 */
2829 if (mfi_ptr->normal[XCOORD] == mfi_ptr->normal[YCOORD] &&
2830 mfi_ptr->normal[XCOORD] == mfi_ptr->normal[ZCOORD]) {
2831 mfi_ptr->normal[ZCOORD] = -1.0;
2832 mfi_ptr->normal[YCOORD] = 0.0;
2833 mfi_ptr->normal[XCOORD] = 0.0;
2834 }
2835
2836 /* Compute slice-to-slice step vector
2837 */
2838 for (i = 0; i < WORLD_NDIMS; i++) {
2839 mfi_ptr->step[i] = spacing * mfi_ptr->normal[i];
2840 }
2841
2842 /* Get position and correct to first slice
2843 */
2844 if (!dicom_read_position(group_list, 0, mfi_ptr->position)) {
2845 if (G.Debug) {
2846 printf("WARNING: No image position found\n");
2847 }
2848 mfi_ptr->position[XCOORD] = mfi_ptr->position[YCOORD] =
2849 mfi_ptr->position[ZCOORD] = 0.0;
2850 }
2851 convert_dicom_coordinate(mfi_ptr->position);
2852
2853 if (G.Debug >= HI_LOGGING) {
2854 printf(" step %.3f %.3f %.3f position %.3f %.3f %.3f\n",
2855 mfi_ptr->step[0],
2856 mfi_ptr->step[1],
2857 mfi_ptr->step[2],
2858 mfi_ptr->position[0],
2859 mfi_ptr->position[1],
2860 mfi_ptr->position[2]);
2861 }
2862
2863 if (!load_image) {
2864 mfi_ptr->big_image = NULL;
2865 mfi_ptr->sub_image = NULL;
2866 }
2867 else {
2868 /* We need to load the image (we're probably on the second pass).
2869 */
2870
2871 /* Steal the image element from the group list
2872 */
2873 element = acr_find_group_element(group_list, ACR_Pixel_data);
2874 if (element == NULL) {
2875 fprintf(stderr, "Couldn't find an image\n");
2876 exit(EXIT_FAILURE);
2877 }
2878
2879 mfi_ptr->big_image = element; /* Save pointer to pixel data element. */
2880
2881 grp_id = acr_get_element_group(element);
2882 elm_id = acr_get_element_element(element);
2883 acr_group_steal_element(acr_find_group(group_list, grp_id), element);
2884
2885 /* Add a small image
2886 */
2887 mfi_ptr->frame_size = rows * cols * pixel_size;
2888 data_ptr = malloc(mfi_ptr->frame_size);
2889 CHKMEM(data_ptr);
2890
2891 mfi_ptr->sub_image = acr_create_element(grp_id,
2892 elm_id,
2893 acr_get_element_vr(element),
2894 mfi_ptr->frame_size,
2895 data_ptr);
2896
2897 copy_element_properties(mfi_ptr->sub_image, mfi_ptr->big_image);
2898
2899 acr_insert_element_into_group_list(&group_list, mfi_ptr->sub_image);
2900 }
2901 }
2902
2903 /* ----------------------------- MNI Header -----------------------------------
2904 @NAME : multiframe_insert_subframe()
2905 @INPUT : group_list - the list of DICOM groups/elements that make up
2906 this file.
2907 mfi_ptr - a pointer to a Multiframe_Info structure that will
2908 contain information used to expand this multiframe file into
2909 a series of slices.
2910 int iimage - the index of the subimage to parse and possibly
2911 to load.
2912 int load_image - a boolean value, non-zero if the function
2913 should actually load the data into memory.
2914 @OUTPUT : mfi_ptr - may be modified by the function.
2915 @RETURNS : void
2916 @DESCRIPTION: This function decomposes a multiframe DICOM image
2917 into a series of single-frame images. Modifies the
2918 group_list to include updated image and position
2919 information.
2920 @METHOD :
2921 @GLOBALS :
2922 @CALLS :
2923 @CREATED : June 3, 2005 Bert Vincent
2924 @MODIFIED :
2925 ---------------------------------------------------------------------------- */
2926 static void
multiframe_insert_subframe(Acr_Group group_list,Multiframe_Info * mfi_ptr,int iframe,int load_image)2927 multiframe_insert_subframe(Acr_Group group_list, Multiframe_Info *mfi_ptr,
2928 int iframe, int load_image)
2929 {
2930 int idim;
2931 char *new_ptr;
2932 char *old_ptr;
2933 double position[WORLD_NDIMS];
2934 string_t string;
2935 int result;
2936
2937 if (G.Debug >= HI_LOGGING) {
2938 printf("multiframe_insert_subframe(%lx, %lx, %d, %d)\n",
2939 (unsigned long)group_list, (unsigned long)mfi_ptr,
2940 iframe, load_image);
2941 }
2942
2943 /* Check the frame number
2944 */
2945 if ((iframe < 0) || (iframe > mfi_ptr->frame_count)) {
2946 fprintf(stderr, "Invalid image number to send: %d of %d\n",
2947 iframe, mfi_ptr->frame_count);
2948 exit(EXIT_FAILURE);
2949 }
2950
2951 /* Update the index in the file's group list.
2952 */
2953 acr_insert_numeric(&group_list, SPI_Current_slice_number, (double) iframe);
2954
2955 result = dicom_read_position(group_list, iframe, position);
2956 convert_dicom_coordinate(position);
2957
2958 if (result != DICOM_POSITION_LOCAL) {
2959 /* If either no position was found for this frame number, or if
2960 * only a global position was found, we need to update the
2961 * position for this particular frame number.
2962 *
2963 * If a local position is found, as in some multiframe files,
2964 * this step is unnecessary and possibly wrong.
2965 */
2966 for (idim = 0; idim < WORLD_NDIMS; idim++) {
2967 position[idim] = mfi_ptr->position[idim] +
2968 (double) iframe * mfi_ptr->step[idim];
2969 }
2970 }
2971
2972 sprintf(string, "%.15g\\%.15g\\%.15g",
2973 position[XCOORD], position[YCOORD], position[ZCOORD]);
2974
2975 acr_insert_string(&group_list, ACR_Image_position_patient, string);
2976
2977 if (G.Debug >= HI_LOGGING) {
2978 printf(" position %s\n", string);
2979 }
2980
2981 if (load_image) {
2982 /* Get pointers
2983 */
2984 old_ptr = acr_get_element_data(mfi_ptr->big_image);
2985 new_ptr = acr_get_element_data(mfi_ptr->sub_image);
2986
2987 /* Copy the image
2988 */
2989 memcpy(new_ptr, /* destination */
2990 old_ptr + (iframe * mfi_ptr->frame_size), /* source */
2991 mfi_ptr->frame_size); /* length */
2992
2993 /* Reset the byte order and VR encoding. This will be modified
2994 * on each send according to what the connection needs.
2995 */
2996 copy_element_properties(mfi_ptr->sub_image, mfi_ptr->big_image);
2997 }
2998 }
2999
3000 static void
multiframe_cleanup(Multiframe_Info * mfi_ptr)3001 multiframe_cleanup(Multiframe_Info *mfi_ptr)
3002 {
3003 if (mfi_ptr) {
3004 if (mfi_ptr->big_image != NULL) {
3005 acr_delete_element(mfi_ptr->big_image);
3006 mfi_ptr->big_image = NULL;
3007 }
3008 }
3009 }
3010
3011