1 /* ----------------------------- MNI Header -----------------------------------
2 @NAME : dicom_read.c
3 @DESCRIPTION: Code to read siemens dicom files and get info from them.
4 @METHOD :
5 @GLOBALS :
6 @CALLS :
7 @CREATED : January 28, 1997 (Peter Neelin)
8 @MODIFIED :
9 * $Log: dicom_read.c,v $
10 * Revision 1.29 2010-11-23 23:30:50 claude
11 * dcm2mnc: fixed seg fault bug (Claude) and added b-matrix (Ilana)
12 *
13 * Revision 1.28 2008-08-12 05:00:23 rotor
14 * * large number of changes from Claude (64 bit and updates)
15 *
16 * Revision 1.27 2008/01/13 09:38:54 stever
17 * Avoid compiler warnings about functions and variables that are defined
18 * but not used. Remove some such functions and variables,
19 * conditionalize some, and move static declarations out of header files
20 * into C files.
21 *
22 * Revision 1.26 2008/01/11 07:17:07 stever
23 * Remove unused variables.
24 *
25 * Revision 1.25 2007/12/09 19:52:15 alex
26 * Fixed error in echo_time extraction
27 *
28 * Revision 1.24 2007/06/08 20:28:57 ilana
29 * added several fields to mincheader (dicom elements and found in ASCONV header)
30 *
31 * Revision 1.23 2007/05/30 15:17:34 ilana
32 * 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
33 *
34 * Revision 1.22 2006/04/09 15:38:02 bert
35 * Add support for standard DTI fields
36 *
37 * Revision 1.21 2005/12/05 16:50:08 bert
38 * Deal with weird XMedCon images
39 *
40 * Revision 1.20 2005/10/26 23:43:35 bert
41 * Latest attempt at getting slice spacing consistent
42 *
43 * Revision 1.19 2005/09/16 22:13:33 bert
44 * Fix slice spacing handling
45 *
46 * Revision 1.18 2005/08/26 21:25:54 bert
47 * Latest changes ported from 1.0 branch
48 *
49 * Revision 1.16.2.5 2005/08/26 03:50:16 bert
50 * Use ACR_Number_of_temporal_positions for number of time slices
51 *
52 * Revision 1.16.2.4 2005/08/18 16:38:43 bert
53 * Minor updates for dealing w/older numaris data
54 *
55 * Revision 1.16.2.3 2005/06/20 22:03:01 bert
56 * Add functions for traversing DICOM sequences and recursively hunting for needed fields.
57 *
58 * Revision 1.16.2.2 2005/06/02 17:04:26 bert
59 * 1. Fix handling of signed values for pixel minimum and maximum, 2. Set found_dircos[] when we adopt the default direction cosines for files with null direction cosines
60 *
61 * Revision 1.16.2.1 2005/05/12 21:16:47 bert
62 * Initial checkin
63 *
64 * Revision 1.16 2005/05/09 15:30:32 bert
65 * Don't allow a rescale slope value of zero
66 *
67 * Revision 1.15 2005/04/28 17:17:57 bert
68 * Set and update new width information fields in a manner analogous to the coordinate fields in the General_Info and File_Info structures
69 *
70 * Revision 1.14 2005/04/20 23:14:04 bert
71 * Remove most SPI_ references
72 *
73 * Revision 1.13 2005/04/20 17:47:38 bert
74 * Fairly major restructuring, added init_general_info() function
75 *
76 * Revision 1.12 2005/04/18 21:43:04 bert
77 * Properly set default minimum and maximum values based on the pixel representation
78 *
79 * Revision 1.11 2005/04/18 21:01:51 bert
80 * Set signed/unsigned flag correctly
81 *
82 * Revision 1.10 2005/04/18 20:43:25 bert
83 * Added some additional debugging information for image position and orientation
84 *
85 * Revision 1.9 2005/04/18 16:22:13 bert
86 * Don't allow non-slice MRI dimensions to grow arbitrarily
87 *
88 * Revision 1.8 2005/04/05 21:49:52 bert
89 * Use rescale slope and intercept to determine the proper slice minimum and maximum
90 *
91 * Revision 1.7 2005/03/29 20:21:44 bert
92 * Fix use of slice spacing; fully check for position information if possible, otherwise create a reasonable position from the slice index
93 *
94 * Revision 1.6 2005/03/14 23:29:35 bert
95 * Support basic dynamic PET fields. Also allocate indices and coordinates arrays for all dimensions, even those we won't use.
96 *
97 * Revision 1.5 2005/03/13 19:37:42 bert
98 * Try to use slice location for coordinate when all else fails, also added one debugging message and a check for PET modality
99 *
100 * Revision 1.4 2005/03/03 20:11:00 bert
101 * Consider patient_id and patient_name when sorting into series. Fix handling of missing direction cosines
102 *
103 * Revision 1.3 2005/03/03 18:59:15 bert
104 * Fix handling of image position so that we work with the older field (0020, 0030) as well as the new (0020, 0032)
105 *
106 * Revision 1.2 2005/03/02 20:18:09 bert
107 * Latest fixes and tweaks
108 *
109 * Revision 1.1 2005/02/17 16:38:10 bert
110 * Initial checkin, revised DICOM to MINC converter
111 *
112 * Revision 1.1.1.1 2003/08/15 19:52:55 leili
113 * Leili's dicom server for sonata
114 *
115 * Revision 1.12 2002/05/01 21:29:34 rhoge
116 * removed MrProt from minc header - encountered files with large strings,
117 * causing seg faults
118 *
119 * Revision 1.11 2002/04/26 03:27:03 rhoge
120 * fixed MrProt problem - replaced fixed lenght char array with malloc
121 *
122 * Revision 1.10 2002/04/08 17:26:34 rhoge
123 * added additional sequence info to minc header
124 *
125 * Revision 1.9 2002/03/27 18:57:50 rhoge
126 * added diffusion b value
127 *
128 * Revision 1.8 2002/03/19 13:13:56 rhoge
129 * initial working mosaic support - I think time is scrambled though.
130 *
131 * Revision 1.7 2001/12/31 18:27:21 rhoge
132 * modifications for dicomreader processing of Numaris 4 dicom files - at
133 * this point code compiles without warning, but does not deal with
134 * mosaiced files. Also will probably not work at this time for Numaris
135 * 3 .ima files. dicomserver may also not be functional...
136 *
137 * Revision 1.6 2000/12/14 21:33:13 rhoge
138 * code modifications to restore measurement loop support that was broken
139 * by changes to support acqusition loop scanning
140 *
141 * Revision 1.5 2000/12/13 13:25:36 rhoge
142 * removed dummy printf from convert_time_to_seconds - turns out that
143 * buggy behaviour was an optimization problem
144 *
145 * Revision 1.4 2000/12/12 19:27:52 rhoge
146 * changed atof(acr_find_string) back to acr_find_double after realizing
147 * that these functions assume string representation
148 *
149 * Revision 1.3 2000/12/12 14:43:22 rhoge
150 * fixed syntax error (dangling comment symbol) from removal of debug
151 * printfs - compiles now
152 *
153 * Revision 1.2 2000/12/11 20:01:44 rhoge
154 * fixed code for frame time computation - ACR vals are strings. Also
155 * inserted dummy fprintf statement in convert_time_to_seconds as this
156 * seems to salvage Linux problems
157 *
158 * Revision 1.1.1.1 2000/11/30 02:13:15 rhoge
159 * imported sources to CVS repository on amoeba
160 * -now always use ACR_Series for run number (seemed to be
161 * problems with test introduced in 6.1)
162 * -added code to detect use of acquisition loop and also to handle
163 * `corrected' siemens acq loop scans
164 * -changed code to use registration time instead of scan time for
165 * session id
166 * -got rid of extraneous acquisition id digit
167 * -added technical information about data acquisition
168 *
169 * Revision 6.2 1999/10/29 17:51:58 neelin
170 * Fixed Log keyword
171 *
172 * Revision 6.1 1999/08/05 20:00:34 neelin
173 * Get acquisition id from series or study element, depending on the
174 * version of the Siemens software.
175 *
176 * Revision 6.0 1997/09/12 13:24:27 neelin
177 * Release of minc version 0.6
178 *
179 * Revision 5.1 1997/09/10 19:36:13 neelin
180 * Small fix to set default direction cosines when they are absent from the
181 * dicom data.
182 *
183 * Revision 5.0 1997/08/21 13:25:26 neelin
184 * Release of minc version 0.5
185 *
186 * Revision 4.1 1997/06/13 12:51:21 neelin
187 * Changed definition of time index and acquisition id to match change
188 * in Siemens dicom software.
189 *
190 * Revision 4.0 1997/05/07 20:06:20 neelin
191 * Release of minc version 0.4
192 *
193 * Revision 1.2 1997/03/11 13:10:48 neelin
194 * Working version of dicomserver.
195 *
196 * Revision 1.1 1997/03/04 20:56:47 neelin
197 * Initial revision
198 *
199 @COPYRIGHT :
200 Copyright 1997 Peter Neelin, McConnell Brain Imaging Centre,
201 Montreal Neurological Institute, McGill University.
202 Permission to use, copy, modify, and distribute this
203 software and its documentation for any purpose and without
204 fee is hereby granted, provided that the above copyright
205 notice appear in all copies. The author and McGill University
206 make no representations about the suitability of this
207 software for any purpose. It is provided "as is" without
208 express or implied warranty.
209 ---------------------------------------------------------------------------- */
210
211 #include "dcm2mnc.h"
212
213 #define IMAGE_NDIMS 2
214
215 static double convert_time_to_seconds(double dicom_time);
216 static void get_intensity_info(Acr_Group group_list, File_Info *file_info);
217 static void get_coordinate_info(Acr_Group group_list, File_Info *file_info,
218 Orientation *orientation,
219 World_Index volume_to_world[VOL_NDIMS],
220 const int sizes[VOL_NDIMS],
221 double dircos[VOL_NDIMS][WORLD_NDIMS],
222 double steps[VOL_NDIMS],
223 double starts[VOL_NDIMS],
224 double coordinate[WORLD_NDIMS]);
225 static void get_general_header_info(Acr_Group group_list,
226 General_Info *general_info);
227 //static void convert_numa3_coordinate(double coordinate[WORLD_NDIMS]);
228 static void get_identification_info(Acr_Group group_list,
229 double *study_id, int *acq_id,
230 int *rec_num, int *image_type);
231
irnd(double x)232 static int irnd(double x)
233 {
234 if (x > 0.0) {
235 x += 0.5;
236 }
237 else {
238 x -= 0.5;
239 }
240 return (int) floor(x);
241 }
242
243 int
is_numaris3(Acr_Group group_list)244 is_numaris3(Acr_Group group_list)
245 {
246 return (strstr(acr_find_string(group_list, ACR_Manufacturer, ""),
247 "SIEMENS") != NULL &&
248 strstr(acr_find_string(group_list, ACR_Software_versions, ""),
249 "VB33") != NULL);
250 }
251
252 /* ----------------------------- MNI Header -----------------------------------
253 @NAME : init_general_info
254 @INPUT : fi_ptr - file-specific info
255 group_list - input data
256 volume_to_world - correspondence of volume to world dimensions
257 spatial_sizes - 3D voxel counts
258 dircos - direction cosines
259 steps - width of each voxel for each spatial dimension
260 starts - starting position for each spatial dimension
261 coordinate -
262 @OUTPUT : gi_ptr - general information about files in this series
263
264
265 @RETURNS : (nothing)
266 @DESCRIPTION: Initializes a "General_Info" structure based upon several
267 bits of information, including a previously initialized
268 File_Info structure (fi_ptr) and the DICOM group_list.
269 Broken out from the get_file_info() function to help
270 simplify that function.
271 @METHOD :
272 @GLOBALS :
273 @CALLS :
274 @CREATED : April 19, 2005 (Bert Vincent)
275 @MODIFIED :
276 ---------------------------------------------------------------------------- */
277 static void
init_general_info(General_Info * gi_ptr,const File_Info * fi_ptr,const Acr_Group group_list,const World_Index volume_to_world[VOL_NDIMS],const int spatial_sizes[VOL_NDIMS],double dircos[VOL_NDIMS][WORLD_NDIMS],const double steps[VOL_NDIMS],const double starts[VOL_NDIMS],double study_id,int acq_id,int rec_num)278 init_general_info(General_Info *gi_ptr, /* OUT */
279 const File_Info *fi_ptr, /* IN */
280 const Acr_Group group_list, /* IN */
281 const World_Index volume_to_world[VOL_NDIMS], /* IN */
282 const int spatial_sizes[VOL_NDIMS], /* IN */
283 double dircos[VOL_NDIMS][WORLD_NDIMS], /* IN */
284 const double steps[VOL_NDIMS], /* IN */
285 const double starts[VOL_NDIMS], /* IN */
286 double study_id, /* IN */
287 int acq_id, /* IN */
288 int rec_num) /* IN */
289 {
290 Acr_Element_Id mri_total_list[MRI_NDIMS];
291 int ivalue; /* For pixel representation value. */
292 World_Index iworld; /* World coordinate index (XCOORD, YCOORD...) */
293 World_Index jworld; /* World coordinate index */
294 Volume_Index ivolume; /* Voxel coordinate index (VROW, VCOLUMN...) */
295 Mri_Index imri;
296 int index;
297
298 // Initialize array for MRI dimension lengths
299 //
300 mri_total_list[SLICE] = ACR_Images_in_acquisition;
301 mri_total_list[ECHO] = ACR_Echo_train_length;
302 mri_total_list[TIME] = ACR_Acquisitions_in_series;
303 mri_total_list[PHASE] = NULL;
304 mri_total_list[CHEM_SHIFT] = NULL;
305
306 // Get row and columns sizes
307 gi_ptr->nrows = spatial_sizes[VROW];
308 gi_ptr->ncolumns = spatial_sizes[VCOLUMN];
309
310 // Save the study, acquisition, reconstruction and image type
311 // identifiers
312 gi_ptr->study_id = study_id;
313 gi_ptr->acq_id = acq_id;
314 gi_ptr->rec_num = rec_num;
315
316 strcpy(gi_ptr->image_type_string, acr_find_string(group_list,
317 ACR_Image_type,
318 ""));
319
320 /* Get dimension information
321 */
322 for (imri = 0; imri < MRI_NDIMS; imri++) {
323
324 /* Get sizes along "MRI" dimensions...
325 */
326 gi_ptr->cur_size[imri] = 1;
327
328 if (mri_total_list[imri] != NULL) {
329 int def_val = 1;
330
331 /* Special case for slice index - need to look at the
332 * new standard element 0x0020:0x1002 "Images in
333 * Acquisition". We use this as a default, but override
334 * it with the Siemens-specific value if present.
335 */
336 if (imri == SLICE) {
337 /* Look for the standard slice count fields first. We
338 * start with the (0054, 0081) first, and if that fails
339 * we retry with (0020, 1002).
340 */
341 def_val = acr_find_int(group_list, ACR_Number_of_slices, 0);
342 if (def_val == 0) {
343 def_val = acr_find_int(group_list,
344 ACR_Images_in_acquisition,
345 1);
346 }
347 }
348
349 if (imri == TIME) {
350 /* Look for the official time slice count field first.
351 */
352 def_val = acr_find_int(group_list,
353 ACR_Number_of_temporal_positions,
354 0);
355
356 def_val = acr_find_int(group_list,
357 ACR_Number_of_time_slices,
358 def_val);
359
360 /*Sometimes need to look further for time count (ex:segmented FLASH)*/
361 def_val = acr_find_int(group_list,
362 ACR_Cardiac_number_of_images,
363 def_val);
364
365 }
366 gi_ptr->max_size[imri] = acr_find_int(group_list,
367 mri_total_list[imri],
368 def_val);
369 }
370 else {
371 gi_ptr->max_size[imri] = 1;
372 }
373
374 if (gi_ptr->max_size[imri] < 1) {
375 gi_ptr->max_size[imri] = 1;
376 }
377
378 /* Check for 3D partitions for slice dimensions */
379 if (imri == SLICE) {
380 /* Get number of 3D partitions for working out number of
381 * slices
382 */
383 int number_of_3D_partitions =
384 acr_find_int(group_list, SPI_Number_of_3D_raw_partitions_nominal, 1);
385 if (number_of_3D_partitions < 1) {
386 number_of_3D_partitions = 1;
387 }
388
389 gi_ptr->max_size[imri] *= number_of_3D_partitions;
390 }
391
392 gi_ptr->default_index[imri] = fi_ptr->index[imri];
393 gi_ptr->image_index[imri] = -1;
394
395 /* Allocate space for index and coordinate arrays.
396 * Set the first values.
397 */
398
399 gi_ptr->indices[imri] = malloc(gi_ptr->max_size[imri] * sizeof(int));
400
401 gi_ptr->coordinates[imri] =
402 malloc(gi_ptr->max_size[imri] * sizeof(double));
403
404 gi_ptr->widths[imri] =
405 malloc(gi_ptr->max_size[imri] * sizeof(double));
406
407 for (index = 0; index < gi_ptr->max_size[imri]; index++) {
408 gi_ptr->indices[imri][index] = -1;
409 gi_ptr->coordinates[imri][index] = 0;
410 gi_ptr->widths[imri][index] = 0; /* default */
411 }
412 gi_ptr->search_start[imri] = 0;
413 gi_ptr->indices[imri][0] = fi_ptr->index[imri];
414 gi_ptr->coordinates[imri][0] = fi_ptr->coordinate[imri];
415 gi_ptr->widths[imri][0] = fi_ptr->width[imri];
416
417 if (G.Debug) {
418 printf("%2d. %s axis length %d\n",
419 imri, Mri_Names[imri], gi_ptr->max_size[imri]);
420 }
421 } /* Loop over dimensions */
422
423 /* Get spatial coordinate information */
424 gi_ptr->slice_world = volume_to_world[VSLICE];
425 gi_ptr->row_world = volume_to_world[VROW];
426 gi_ptr->column_world = volume_to_world[VCOLUMN];
427 for (ivolume = 0; ivolume < VOL_NDIMS; ivolume++) {
428 iworld = volume_to_world[ivolume];
429 gi_ptr->step[iworld] = steps[ivolume];
430 gi_ptr->start[iworld] = starts[ivolume];
431 for (jworld = 0; jworld < WORLD_NDIMS; jworld++) {
432 gi_ptr->dircos[iworld][jworld] = dircos[ivolume][jworld];
433 }
434
435 if (G.Debug) {
436 printf("%2d. %s axis length %d step %.3f, start %.3f, cosines %.3f,%.3f,%.3f\n",
437 ivolume,
438 World_Names[iworld],
439 spatial_sizes[ivolume],
440 gi_ptr->step[iworld],
441 gi_ptr->start[iworld],
442 gi_ptr->dircos[iworld][XCOORD],
443 gi_ptr->dircos[iworld][YCOORD],
444 gi_ptr->dircos[iworld][ZCOORD]
445 );
446 }
447 }
448
449 /* Set data type and range */
450 if (fi_ptr->bits_alloc <= 8) {
451 gi_ptr->datatype = NC_BYTE;
452 }
453 else {
454 gi_ptr->datatype = NC_SHORT;
455 }
456
457 /* bert- modify code to correctly read the pixel
458 * representation if available and use that to set the
459 * signed/unsigned flag.
460 */
461
462 ivalue = acr_find_int(group_list, ACR_Pixel_representation, -1);
463 if (ivalue == ACR_PIXEL_REP_UNSIGNED) {
464 gi_ptr->is_signed = 0;
465 }
466 else if (ivalue == ACR_PIXEL_REP_SIGNED) {
467 gi_ptr->is_signed = 1;
468 }
469 else {
470 if (ivalue != -1) {
471 printf("WARNING: Unknown pixel representation value %d\n",
472 ivalue);
473 }
474 gi_ptr->is_signed = ((gi_ptr->datatype == NC_SHORT) &&
475 (fi_ptr->bits_stored < 16));
476 }
477
478 gi_ptr->pixel_min = fi_ptr->pixel_min;
479 gi_ptr->pixel_max = fi_ptr->pixel_max;
480
481 /* Save display window info */
482 gi_ptr->window_min = fi_ptr->window_min;
483 gi_ptr->window_max = fi_ptr->window_max;
484
485 /* Get the rest of the header information */
486 get_general_header_info(group_list, gi_ptr);
487
488 /* Copy the group list */
489 gi_ptr->group_list = acr_copy_group_list(group_list);
490
491 // note that number of slices will be wrong here for mosaics
492 // we add some other mosaic-relevant fields...
493
494 gi_ptr->num_mosaic_rows =
495 acr_find_int(group_list, EXT_Mosaic_rows, 1);
496 gi_ptr->num_mosaic_cols =
497 acr_find_int(group_list, EXT_Mosaic_columns, 1);
498 gi_ptr->num_slices_in_file =
499 acr_find_int(group_list, EXT_Slices_in_file, 1);
500 gi_ptr->sub_image_rows =
501 (int)acr_find_short(group_list, EXT_Sub_image_rows,
502 acr_find_short(group_list, ACR_Rows, 0));
503 gi_ptr->sub_image_columns =
504 (int)acr_find_short(group_list, EXT_Sub_image_columns,
505 acr_find_short(group_list, ACR_Rows, 0));
506
507 /* Set initialized flag */
508 gi_ptr->initialized = TRUE;
509
510 if (G.Debug) {
511 printf("Pixel minimum %.10f maximum %.10f\n",
512 gi_ptr->pixel_min, gi_ptr->pixel_max);
513 printf("Window minimum %.10f maximum %.10f\n",
514 gi_ptr->window_min, gi_ptr->window_max);
515 }
516 }
517
518 /* ----------------------------- MNI Header -----------------------------------
519 @NAME : get_file_info
520 @INPUT : group_list - input data
521 @OUTPUT : file_inf_p - file-specific info
522 general_info - general information about files
523 @RETURNS : (nothing)
524 @DESCRIPTION: Routine to extract information from a group list
525 @METHOD :
526 @GLOBALS :
527 @CALLS :
528 @CREATED : November 25, 1993 (Peter Neelin)
529 @MODIFIED : Modified Feb. 2000 by Rick Hoge to handle Acquisition loop mode
530 : if assume_acq_loop is FALSE the other added variables don't
531 : matter
532 ---------------------------------------------------------------------------- */
533 void
get_file_info(Acr_Group group_list,File_Info * fi_ptr,General_Info * gi_ptr)534 get_file_info(Acr_Group group_list, File_Info *fi_ptr, General_Info *gi_ptr)
535 {
536 Mri_Index imri; /* MRI index (SLICE, ECHO, TIME, PHASE...) */
537 int nrows; /* Row count in this file */
538 int ncolumns; /* Column count in this file */
539 int spatial_sizes[VOL_NDIMS]; /* Voxel coordinate extents */
540 double study_id; /* Study identifier */
541 int acq_id; /* Acquisition identifier */
542 int rec_num; /* ? Seems to be a dummy */
543 int cur_index; /* Index of slice(s) in current file */
544 int index; /* General index value */
545 Orientation orientation; /* TRANSVERSE, SAGITTAL, or CORONAL */
546 World_Index volume_to_world[VOL_NDIMS]; /* Maps voxel to world indices */
547 double coordinate[WORLD_NDIMS]; /* Slice coordinates */
548 double dircos[VOL_NDIMS][WORLD_NDIMS]; /* Direction cosines */
549 double steps[VOL_NDIMS]; /* Step (spacing) coordinates */
550 double starts[VOL_NDIMS]; /* Start (origin) coordinates */
551 Acr_Element_Id mri_index_list[MRI_NDIMS];
552 Acr_Element_Id mri_total_list[MRI_NDIMS]; /*added by ilana*/
553 Acr_Element element;
554
555 // Array of elements for mri dimensions
556 mri_index_list[SLICE] = SPI_Current_slice_number;
557 mri_index_list[ECHO] = ACR_Echo_number;
558 // dicom time element may be different on old systems
559 mri_index_list[TIME] = ACR_Acquisition;
560 mri_index_list[PHASE] = NULL;
561 mri_index_list[CHEM_SHIFT] = NULL;
562 mri_total_list[SLICE] = SPI_Number_of_slices_nominal;
563 mri_total_list[ECHO] = SPI_Number_of_echoes;
564 mri_total_list[TIME] = ACR_Acquisitions_in_series;
565 mri_total_list[PHASE] = NULL;
566 mri_total_list[CHEM_SHIFT] = NULL;
567
568 /* Get image dimensions
569 */
570 nrows = (int)acr_find_short(group_list, ACR_Rows, 0);
571 ncolumns = (int)acr_find_short(group_list, ACR_Columns, 0);
572
573 spatial_sizes[VROW] = nrows;
574 spatial_sizes[VCOLUMN] = ncolumns;
575
576 spatial_sizes[VSLICE] = acr_find_int(group_list, ACR_Images_in_acquisition,
577 1);
578
579 /* Get intensity information
580 */
581 get_intensity_info(group_list, fi_ptr);
582
583 /* Check for necessary values not found
584 */
585 if ((nrows <= 0) || (ncolumns <= 0) ||
586 (fi_ptr->bits_stored <= 0) ||
587 (fi_ptr->bits_alloc <= 0)) {
588 if (G.Debug) {
589 printf("ERROR: Needed values missing, marking invalid\n");
590 }
591 fi_ptr->valid = FALSE;
592 return;
593 }
594
595 /* Get study, acq, rec, image type id's
596 */
597 get_identification_info(group_list, &study_id, &acq_id, &rec_num, NULL);
598
599 /* Get indices for image in current file
600 */
601 for (imri=0; imri < MRI_NDIMS; imri++) {
602
603 if (mri_index_list[imri] != NULL) {
604 fi_ptr->index[imri] =
605 /* note that for TIME this will use ACR_Acqusition,
606 which does not work for measurement loop scans */
607 acr_find_int(group_list, mri_index_list[imri], 1);
608 }
609 else {
610 fi_ptr->index[imri] = 1;
611 }
612 }
613
614 /* Get coordinate information
615 */
616 get_coordinate_info(group_list, fi_ptr, &orientation, volume_to_world,
617 spatial_sizes, dircos, steps, starts, coordinate);
618
619 /*
620 * Use the coordinate information rather than the slice or time
621 * position derived above. This seems to be much more reliable.
622 */
623
624 fi_ptr->index[SLICE] = irnd(fi_ptr->coordinate[SLICE] * 100.0);
625
626 /* For non-IMA files, use the time coordinate to order the time.
627 * index. IMA files do not seem to have a reliable slice time
628 * indicator.
629 */
630 /*if (G.file_type != IMA) {
631 fi_ptr->index[TIME] = irnd(fi_ptr->coordinate[TIME] * 100.0);
632 }*/
633
634 //ilana debug
635 //int test = fi_ptr->coordinate[TIME];
636 //print ("*******Time coordinate IS: %i\n\n",test);
637
638 /* Set up general info on first pass
639 */
640 if (!gi_ptr->initialized) {
641 init_general_info(gi_ptr,
642 fi_ptr,
643 group_list,
644 volume_to_world,
645 spatial_sizes,
646 dircos,
647 steps,
648 starts,
649 study_id,
650 acq_id,
651 rec_num);
652 }
653
654 /* Set up file info */
655
656 /* Update general info and validate file on later passes
657 */
658 else {
659
660 /* Check for consistent pixel minimum and maximum. */
661 if ((gi_ptr->pixel_max != fi_ptr->pixel_max) ||
662 (gi_ptr->pixel_min != fi_ptr->pixel_min)) {
663 printf("WARNING: Inconsistent pixel minimum and maximum\n");
664 printf(" %f %f, %f %f\n",
665 gi_ptr->pixel_min, gi_ptr->pixel_max,
666 fi_ptr->pixel_min, fi_ptr->pixel_max);
667 #if 0
668 if (gi_ptr->pixel_max < fi_ptr->pixel_max) {
669 gi_ptr->pixel_max = fi_ptr->pixel_max;
670 }
671 if (gi_ptr->pixel_min > fi_ptr->pixel_min) {
672 gi_ptr->pixel_min = fi_ptr->pixel_min;
673 }
674 #endif
675 }
676
677 /* Check for consistent data type */
678 if (((gi_ptr->datatype == NC_BYTE) && (fi_ptr->bits_alloc > 8)) ||
679 ((gi_ptr->datatype == NC_SHORT) && (fi_ptr->bits_alloc <= 8))) {
680 printf("Inconsistent datatype, marking invalid\n");
681 fi_ptr->valid = FALSE;
682 return;
683 }
684
685 /* Check row and columns sizes */
686 if ((nrows != gi_ptr->nrows) && (ncolumns != gi_ptr->ncolumns)) {
687 printf("Mismatched rows/columns, marking invalid\n");
688 fi_ptr->valid = FALSE;
689 return;
690 }
691
692 /* Check study and acquisition id's */
693 if ((gi_ptr->study_id != study_id) || (gi_ptr->acq_id != acq_id)) {
694 printf("Mismatched acquisition/study, marking invalid\n");
695 fi_ptr->valid = FALSE;
696 return;
697 }
698
699 /* Look to see if indices have changed */
700 for (imri = 0; imri < MRI_NDIMS; imri++) {
701 /* If a dimension is known to have a maximum size of one
702 * or less, we do NOT allow it to grow in any way. An
703 * exception is made for the slice dimension, however,
704 * since it appears that it is common for it to be
705 * unspecified and can be guessed only by the number of
706 * distinct locations discovered.
707 */
708 if (imri != SLICE && gi_ptr->max_size[imri] <= 1) {
709 continue;
710 }
711
712 /* Get current index */
713 cur_index = fi_ptr->index[imri];
714
715 //ilana debug
716 //print ("************CURRENT INDEX IS: %i\n\n",cur_index);
717
718 /* Check whether this index is in the list.
719 */
720 if (gi_ptr->cur_size[imri] == 1) {
721 index = ((cur_index == gi_ptr->default_index[imri]) ? 0 : -1);
722 }
723 else {
724 /* Search list of indices for 'cur_index'. Search is
725 started at search_start[] and has maximum length of
726 size[imri].
727 */
728 index = search_list(cur_index,
729 gi_ptr->indices[imri],
730 gi_ptr->cur_size[imri],
731 gi_ptr->search_start[imri]);
732 }
733
734 /* If it is not, then add it */
735 if (index < 0) {
736 if (G.Debug >= HI_LOGGING) {
737 printf("Need to add index %d to %s list, %d/%d\n",
738 cur_index, Mri_Names[imri],
739 gi_ptr->cur_size[imri],
740 gi_ptr->max_size[imri]);
741 }
742
743 /* Check whether we can add a new index */
744 if (gi_ptr->cur_size[imri] >= gi_ptr->max_size[imri]) {
745 gi_ptr->max_size[imri]++;
746 gi_ptr->indices[imri] =
747 realloc(gi_ptr->indices[imri],
748 gi_ptr->max_size[imri] * sizeof(int));
749
750 gi_ptr->coordinates[imri] =
751 realloc(gi_ptr->coordinates[imri],
752 gi_ptr->max_size[imri] * sizeof(double));
753
754 gi_ptr->widths[imri] =
755 realloc(gi_ptr->widths[imri],
756 gi_ptr->max_size[imri] * sizeof(double));
757 }
758
759
760 /* Add the index and coordinate to the lists */
761 index = gi_ptr->cur_size[imri];
762 gi_ptr->search_start[imri] = index;
763 gi_ptr->indices[imri][index] = cur_index;
764 gi_ptr->coordinates[imri][index] = fi_ptr->coordinate[imri];
765 gi_ptr->widths[imri][index] = fi_ptr->width[imri];
766 gi_ptr->cur_size[imri]++;
767
768 }
769 } /* Loop over Mri_Index */
770
771 // Update display window info
772 if (gi_ptr->window_min > fi_ptr->window_min)
773 gi_ptr->window_min = fi_ptr->window_min;
774 if (gi_ptr->window_max < fi_ptr->window_max)
775 gi_ptr->window_max = fi_ptr->window_max;
776
777 } // Update general info for this file
778
779 /* Get DTI information if available.
780 */
781 fi_ptr->b_value = (double)acr_find_double(group_list, ACR_Diffusion_b_value, -1);
782
783 element = acr_find_group_element(group_list,
784 ACR_Diffusion_gradient_orientation);
785
786 Acr_Double grad_direction[WORLD_NDIMS];
787 if (element == NULL ||
788 acr_get_element_double_array(element, WORLD_NDIMS, grad_direction) != WORLD_NDIMS) {
789 grad_direction[XCOORD] =
790 grad_direction[YCOORD] =
791 grad_direction[ZCOORD] = 0; /*this should be 0 for the b=0 images*/
792 }
793 fi_ptr->grad_direction[XCOORD] = (double)grad_direction[XCOORD];
794 fi_ptr->grad_direction[YCOORD] = (double)grad_direction[YCOORD];
795 fi_ptr->grad_direction[ZCOORD] = (double)grad_direction[ZCOORD];
796
797 /*Want to add B-matrix if it exists, the problem is that
798 the b=0 image does not even have the field (0019x1027) while
799 the other directions might (AND b=0 is ofen the first image!).
800 Initialize to all 0s and decide later whether to include in header? */
801 if (gi_ptr->acq.dti) {
802 element = acr_find_group_element(group_list, ACR_B_matrix);
803 int num_bmatrix_elements=6;/*bmatrix should have 6 values*/
804 if (element == NULL ||
805 acr_get_element_double_array(element, num_bmatrix_elements,
806 fi_ptr->b_matrix) != num_bmatrix_elements) {
807 int i=0;
808 for(i;i<num_bmatrix_elements;i++){
809 fi_ptr->b_matrix[i]=0; /*bmatrix for b=0 should be 0*/
810 }
811 }
812 }
813
814 // If we get to here, then we have a valid file
815 fi_ptr->valid = TRUE;
816 return;
817 }
818
819 /* ----------------------------- MNI Header -----------------------------------
820 @NAME : get_identification_info
821 @INPUT : group_list - input data
822 @OUTPUT : study_id
823 acq_id
824 rec_num
825 image_type
826 @RETURNS : (nothing)
827 @DESCRIPTION: Routine to get image identification information.
828 @METHOD :
829 @GLOBALS :
830 @CALLS :
831 @CREATED : February 28, 1997 (Peter Neelin)
832 @MODIFIED :
833 ---------------------------------------------------------------------------- */
834 static void
get_identification_info(Acr_Group group_list,double * study_id,int * acq_id,int * rec_num,int * image_type)835 get_identification_info(Acr_Group group_list,
836 double *study_id, int *acq_id,
837 int *rec_num, int *image_type)
838 {
839 int number_of_frames;
840 int number_of_averages;
841
842 if (study_id != NULL) {
843 // generate a study ID number from date & time: yyyymmdd.hhmmss
844 // (should be unique enough for our application)
845 *study_id = (((float)acr_find_int(group_list, ACR_Study_date, 0)) +
846 ((float)acr_find_int(group_list, ACR_Study_time, 0))/1e6);
847 }
848 if (acq_id != NULL) {
849
850 *acq_id = acr_find_int(group_list, ACR_Series, 0);
851
852 number_of_frames =
853 acr_find_int(group_list, ACR_Acquisitions_in_series, 1);
854
855 number_of_averages =
856 acr_find_int(group_list, ACR_Nr_of_averages, 1);
857
858 /* Determine if measurement loop was used (rhoge) -
859
860 if so, we replace the different series numbers with
861 ACR_Study_time, which is the same for all frames in a run.
862 This will aid in grouping the files later on.
863
864 Criteria used for identification of meast loop:
865
866 1) more than one dynamic scan
867
868 2) number of dynamic scans NOT equal to nominal number of signal
869 averages (if they're equal, we assume acquisition loop scan)
870
871 WARNING: it is possible that someone might use the
872 measurement loop do serial scans which each have multiple signal
873 averages. If NSA = the number of measts. in this case, then
874 the scan will not be recognized as a Meast loop scan and the
875 different frames will be placed in different series. To fix
876 this, we'd really need to look at the series numbers of
877 future and past files. It's also unlikely to happen but I'm
878 sure it will...
879
880 This also means we should NOT correct the number of signal
881 averages on the sending side */
882
883
884 /* if ((number_of_frames > 1) || (*acq_id == 0)) { (orig test) */
885
886 if ( (G.file_type == N3DCM) &&
887 (number_of_frames > 1) &&
888 (number_of_frames != number_of_averages)) {
889
890 *acq_id = acr_find_int(group_list, ACR_Study_time, 0);
891
892 }
893 }
894 if (rec_num != NULL)
895 *rec_num = 0;
896 if (image_type != NULL)
897 *image_type = 0;
898 }
899
900 /* ----------------------------- MNI Header -----------------------------------
901 @NAME : get_intensity_info
902 @INPUT : group_list - input data
903 @OUTPUT : fi_ptr - file-specific info
904 @RETURNS : (nothing)
905 @DESCRIPTION: Routine to get intensity information from a group list
906 @METHOD :
907 @GLOBALS :
908 @CALLS :
909 @CREATED : February 28, 1997 (Peter Neelin)
910 @MODIFIED :
911 ---------------------------------------------------------------------------- */
912 static void
get_intensity_info(Acr_Group group_list,File_Info * fi_ptr)913 get_intensity_info(Acr_Group group_list, File_Info *fi_ptr)
914 {
915 double window_centre, window_width;
916 double rescale_intercept, rescale_slope;
917 int ivalue; /* 0000 for unsigned, 0001 for signed */
918 int imin, imax; /* Default minimum and maximum values */
919
920 /* Get pixel storage information */
921 fi_ptr->bits_alloc = (int)acr_find_short(group_list, ACR_Bits_allocated, 0);
922 fi_ptr->bits_stored = (int)acr_find_short(group_list, ACR_Bits_stored, 0);
923
924 /* bert- properly set the minimum and maximum pixel values depending
925 * on whether or not this file specifies signed pixel values.
926 */
927 ivalue = acr_find_int(group_list, ACR_Pixel_representation, -1);
928
929 if (ivalue == ACR_PIXEL_REP_SIGNED) {
930 imin = -(1 << (fi_ptr->bits_stored - 1));
931 imax = (1 << (fi_ptr->bits_stored - 1)) - 1;
932 }
933 else {
934 imin = 0;
935 imax = (1 << fi_ptr->bits_stored) - 1;
936 }
937
938 if (G.useMinMax) {
939 int pmin;
940 int pmax;
941
942 /* Get pixel value information
943 * I think this might wrongly assume that the ACR values min/max
944 * apply to the whole volume - it actually appears that they apply
945 * to the current slice.
946 */
947
948 pmin = (int)acr_find_short(group_list, ACR_Smallest_pixel_value, imin);
949 pmax = (int)acr_find_short(group_list, ACR_Largest_pixel_value, imax);
950
951 /* Hack to convert to signed representation if indicated - if
952 * bit 15 is set, we have to "promote" to 2's complement by
953 * sign-extending the result before converting it to a double.
954 * Perhaps this sort of thing needs to be pushed down into the
955 * ACR-NEMA library somehow, so that the representations of
956 * short values are automatically converted properly?
957 */
958 if (ivalue == ACR_PIXEL_REP_SIGNED) {
959 if (pmin & 0x8000) {
960 pmin -= 0x10000;
961 }
962 if (pmax & 0x8000) {
963 pmax -= 0x10000;
964 }
965 }
966
967 fi_ptr->pixel_min = (double) pmin;
968 fi_ptr->pixel_max = (double) pmax;
969 }
970 else {
971 /* for now, use bits_stored to determine dynamic range
972 * DICOM info on largest pixel applies to first slice,
973 * not whole volume - this caused problems (roundoff?)
974 * in Siemens Numaris 4 scans
975 */
976 fi_ptr->pixel_min = imin;
977 fi_ptr->pixel_max = imax;
978 }
979
980 /* Get the rescale intercept and slope. If they are not present,
981 * we use the default values of 0.0 for the intercept and 1.0 for
982 * the slope.
983 */
984 rescale_intercept = (double)acr_find_double(group_list, ACR_Rescale_intercept, 0);
985 rescale_slope = (double)acr_find_double(group_list, ACR_Rescale_slope, 1);
986
987 /* If the rescale slope is set to zero, force the default value of
988 * one and issue a warning.
989 */
990 if (rescale_slope == 0.0) {
991 printf("WARNING: File contains a rescale slope value of zero.\n");
992 rescale_slope = 1.0;
993 }
994
995 fi_ptr->slice_min = fi_ptr->pixel_min * rescale_slope + rescale_intercept;
996 fi_ptr->slice_max = fi_ptr->pixel_max * rescale_slope + rescale_intercept;
997
998 /* Get window min and max */
999 window_centre = (fi_ptr->slice_max + fi_ptr->slice_min) / 2.0;
1000 window_width = fi_ptr->slice_max - fi_ptr->slice_min;
1001 window_centre =
1002 (double)acr_find_double(group_list, ACR_Window_centre, window_centre);
1003 window_width =
1004 (double)acr_find_double(group_list, ACR_Window_width, window_width);
1005 fi_ptr->window_min = window_centre - window_width / 2.0;
1006 fi_ptr->window_max = window_centre + window_width / 2.0;
1007
1008 }
1009
1010 /* Function to recursively search an element list for a specific
1011 * element, skipping a specified number of occurrences before
1012 * returning. This is only called by acr_recurse_for_element().
1013 */
1014 static Acr_Element
acr_recursive_search(Acr_Element el_lst,int * nskip,Acr_Element_Id srch_id)1015 acr_recursive_search(Acr_Element el_lst, int *nskip, Acr_Element_Id srch_id)
1016 {
1017 Acr_Element el_ret = NULL;
1018 Acr_Element el_tmp;
1019
1020 for (el_tmp = el_lst; el_tmp != NULL;
1021 el_tmp = acr_get_element_next(el_tmp)) {
1022
1023 /* If we find what we're looking for, return it.
1024 */
1025 if (acr_get_element_group(el_tmp) == srch_id->group_id &&
1026 acr_get_element_element(el_tmp) == srch_id->element_id) {
1027 if (*nskip <= 0) {
1028 el_ret = el_tmp;
1029 break;
1030 }
1031 else {
1032 --(*nskip);
1033 }
1034 }
1035 /* See if we need to recurse.
1036 */
1037 if (acr_element_is_sequence(el_tmp)) {
1038 el_lst = (Acr_Element) acr_get_element_data(el_tmp);
1039 el_ret = acr_recursive_search(el_lst, nskip, srch_id);
1040 if (el_ret != NULL) {
1041 break;
1042 }
1043 }
1044 }
1045 return (el_ret);
1046 }
1047
1048 /* acr_recurse_for_element()
1049 *
1050 * Function to search a group list for a particular element. Unlike other
1051 * functions along these lines, this function will recursively descend into
1052 * compound datatypes (DICOM sequences) to hunt for instances of a particular
1053 * element.
1054 *
1055 * The search proceeds in two stages: The first is to search for
1056 * a particular sequence object. This must be found or else the search is
1057 * called off. Once the expected sequences is found, the function will
1058 * recursively search all of the substructure of that sequence for the
1059 * requested subelement. The "nskip" parameter tells the function to ignore
1060 * the first "nskip" matches that it locates.
1061 */
1062 Acr_Element
acr_recurse_for_element(Acr_Group group_list,int nskip,Acr_Element_Id seq_id,Acr_Element_Id srch_id)1063 acr_recurse_for_element(Acr_Group group_list,
1064 int nskip,
1065 Acr_Element_Id seq_id,
1066 Acr_Element_Id srch_id)
1067 {
1068 Acr_Element el_seq;
1069
1070 /* Hunt for the necessary sequence object.
1071 */
1072 el_seq = acr_find_group_element(group_list, seq_id);
1073 if (el_seq == NULL || !acr_element_is_sequence(el_seq)) {
1074 /* If not found, or not a sequence, abort the search.
1075 */
1076 return (NULL);
1077 }
1078
1079 /* Otherwise proceed to "stage 2" and hunt for the requested subelement.
1080 */
1081 return acr_recursive_search((Acr_Element) acr_get_element_data(el_seq),
1082 &nskip, srch_id);
1083 }
1084
1085 int
dicom_read_position(Acr_Group group_list,int n,double coordinate[3])1086 dicom_read_position(Acr_Group group_list, int n, double coordinate[3])
1087 {
1088 Acr_Element element;
1089 int result;
1090
1091 /* Try to read a unique element from the sequences. If this
1092 * succeeds, we need to flag this fact so that the higher-level
1093 * processing can adapt accordingly.
1094 */
1095 element = acr_recurse_for_element(group_list, n,
1096 ACR_Perframe_func_groups_seq,
1097 ACR_Image_position_patient);
1098 if (element != NULL) {
1099 result = DICOM_POSITION_LOCAL; /* Found a slice-specific position */
1100 }
1101 else {
1102 result = DICOM_POSITION_GLOBAL; /* Found a global position */
1103
1104 /* bert-look for field in weird XMedCon location
1105 */
1106 element = acr_recurse_for_element(group_list, 0,
1107 ACR_Detector_information_seq,
1108 ACR_Image_position_patient);
1109
1110 if (element == NULL) {
1111 element = acr_find_group_element(group_list,
1112 ACR_Image_position_patient);
1113 }
1114
1115 if (element == NULL) {
1116 element = acr_find_group_element(group_list,
1117 ACR_Image_position_patient_old);
1118 }
1119 }
1120
1121 if (element == NULL) {
1122 printf("WARNING: Failed to find image position\n");
1123 }
1124 else {
1125 if (acr_get_element_numeric_array(element,
1126 WORLD_NDIMS,
1127 coordinate) == WORLD_NDIMS) {
1128 return (result);
1129 }
1130
1131 if (G.Debug) {
1132 printf("WARNING: Failed to read image position ('%s')\n",
1133 (char *)acr_get_element_string(element));
1134 }
1135 }
1136 return DICOM_POSITION_NONE;
1137 }
1138
1139 int
dicom_read_orientation(Acr_Group group_list,double orientation[6])1140 dicom_read_orientation(Acr_Group group_list, double orientation[6])
1141 {
1142 Acr_Element element;
1143 int result;
1144
1145 /* read in row/col vectors:
1146 */
1147
1148 /* Look for single global header with slices all in same file. */
1149
1150 element = acr_recurse_for_element(group_list, 0,
1151 ACR_Shared_func_groups_seq,
1152 ACR_Image_orientation_patient);
1153
1154 /* Look for single local header with slices all in same file,
1155 with one header per slice (Philips Intera). */
1156
1157 if (element == NULL) {
1158 element = acr_recurse_for_element(group_list, 0,
1159 ACR_Perframe_func_groups_seq,
1160 ACR_Image_orientation_patient);
1161 }
1162
1163 /* bert - deal with weird XMedCon images...
1164 */
1165 if (element == NULL) {
1166 element = acr_recurse_for_element(group_list, 0,
1167 ACR_Detector_information_seq,
1168 ACR_Image_orientation_patient);
1169 }
1170
1171 /* Look for header with slices in separate files. */
1172 if (element == NULL) {
1173 element = acr_find_group_element(group_list,
1174 ACR_Image_orientation_patient);
1175 }
1176 if (element == NULL) {
1177 /* If we failed to find the newer, better patient orientation
1178 * information, try to use the obsolete information if present.
1179 */
1180 element = acr_find_group_element(group_list,
1181 ACR_Image_orientation_patient_old);
1182 }
1183 if (element == NULL) {
1184 printf("WARNING: Failed to find image orientation!\n");
1185 return (0);
1186 }
1187 else if ((result = acr_get_element_numeric_array(element, 6,
1188 orientation)) != 6) {
1189 printf("WARNING: Failed to read image orientation! (%d, '%s')\n",
1190 result, (char *)acr_get_element_string(element));
1191 return (0);
1192 }
1193 return (1);
1194 }
1195
1196 /*
1197 * Read the pixel size, an array of 2 floating point numbers, from the
1198 * DICOM group list.
1199 */
1200 int
dicom_read_pixel_size(Acr_Group group_list,double pixel_size[IMAGE_NDIMS])1201 dicom_read_pixel_size(Acr_Group group_list, double pixel_size[IMAGE_NDIMS])
1202 {
1203 Acr_Element element;
1204 int result = 0;
1205 int i;
1206
1207 for (i = 0; i < IMAGE_NDIMS; i++) {
1208 pixel_size[i] = -DBL_MAX;
1209 }
1210
1211 /* Look for single global header with slices all in same file. */
1212
1213 element = acr_recurse_for_element(group_list, 0,
1214 ACR_Shared_func_groups_seq,
1215 ACR_Pixel_size);
1216
1217 /* Look for single local header with slices all in same file,
1218 with one header per slice (Philips Intera). */
1219
1220 if (element == NULL) {
1221 element = acr_recurse_for_element(group_list, 0,
1222 ACR_Perframe_func_groups_seq,
1223 ACR_Pixel_size);
1224 }
1225
1226 /* Look for header with slices in separate files. */
1227
1228 if (element == NULL) {
1229 element = acr_find_group_element(group_list, ACR_Pixel_size);
1230 }
1231
1232 if (element == NULL) {
1233 printf("WARNING: Can't find pixel size element\n");
1234 } else {
1235 /* Extract sizes from header string if found. */
1236 if (acr_get_element_numeric_array(element, IMAGE_NDIMS,
1237
1238 pixel_size) != IMAGE_NDIMS) {
1239 printf("WARNING: Can't read pixel size element\n");
1240 }
1241 else {
1242 result = 1;
1243 }
1244 }
1245
1246 /* If the values are still uninitialized, set them to some reasonable
1247 * defaults.
1248 */
1249 if (pixel_size[0] == -DBL_MAX) {
1250 pixel_size[0] = 1.0; /* Assume 1mm spacing */
1251 }
1252
1253 if (pixel_size[1] == -DBL_MAX) {
1254 pixel_size[1] = pixel_size[0]; /* Assume uniform sample grid. */
1255 }
1256 return (result);
1257 }
1258
1259 /* ----------------------------- MNI Header -----------------------------------
1260 @NAME : get_coordinate_info
1261 @INPUT : group_list - input data
1262 sizes - size of each spatial dimension
1263 @OUTPUT : fi_ptr - file-specific info
1264 volume_to_world - volume index to world coordinate index mapping
1265 dircos - direction cosines for spatial dimensions
1266 steps - step sizes for spatial dimensions
1267 starts - start positions for spatial dimensions (for a slice)
1268 coordinate - coordinate of centre of slice
1269 @RETURNS : (nothing)
1270 @DESCRIPTION: Routine to get coordinate information for a slice from
1271 a group list
1272 @METHOD :
1273 @GLOBALS :
1274 @CALLS :
1275 @CREATED : February 28, 1997 (Peter Neelin)
1276 @MODIFIED :
1277 ---------------------------------------------------------------------------- */
1278
1279 static void
get_coordinate_info(Acr_Group group_list,File_Info * fi_ptr,Orientation * orientation,World_Index volume_to_world[VOL_NDIMS],const int sizes[VOL_NDIMS],double dircos[VOL_NDIMS][WORLD_NDIMS],double steps[VOL_NDIMS],double starts[VOL_NDIMS],double coordinate[WORLD_NDIMS])1280 get_coordinate_info(Acr_Group group_list,
1281 File_Info *fi_ptr,
1282 Orientation *orientation,
1283 World_Index volume_to_world[VOL_NDIMS],
1284 const int sizes[VOL_NDIMS],
1285 double dircos[VOL_NDIMS][WORLD_NDIMS],
1286 double steps[VOL_NDIMS],
1287 double starts[VOL_NDIMS],
1288 double coordinate[WORLD_NDIMS])
1289 {
1290 Volume_Index ivolume;
1291 World_Index iworld;
1292 int found_dircos[VOL_NDIMS];
1293 int found_coordinate;
1294 double frame_time;
1295 double start_time;
1296 double magnitude;
1297 double largest;
1298 double psize[IMAGE_NDIMS];
1299 double slice_thickness, slice_spacing;
1300
1301 double RowColVec[6]; /* row/column unit vectors in public dicom element */
1302
1303 const Orientation orientation_list[WORLD_NDIMS] = {
1304 SAGITTAL, CORONAL, TRANSVERSE
1305 };
1306
1307 if (G.Debug >= HI_LOGGING) {
1308 printf("get_coordinate_info(%lx, ...)\n", (unsigned long) group_list);
1309 }
1310
1311 /* Initialize a few things... */
1312 for (ivolume = 0; ivolume < VOL_NDIMS; ivolume++) {
1313 found_dircos[ivolume] = FALSE;
1314 }
1315 found_coordinate = FALSE;
1316
1317 #if 0
1318 /* TODO: For now this appears to be necessary. In cases I don't fully
1319 * understand, the Siemens Numaris 3 DICOM image orientation does not
1320 * give the correct direction cosines, so we use the nonstandard Siemens
1321 * fields instead. Someday I should figure out the relation (if any)
1322 * between the standard fields and these fields, and try to normalize
1323 * this mess.
1324 *
1325 * We only attempt this for files that are clearly marked as SIEMENS
1326 * files, with a version string that looks like VB33 (VB33D, VB33G, etc.)
1327 * Later versions do not seem to use these fields.
1328 */
1329 if (is_numaris3(group_list)) {
1330 Acr_Element_Id dircos_elid[VOL_NDIMS];
1331
1332 /* Set direction cosine element ids. Note that the reversal of
1333 * rows and columns is intentional - their idea of the meaning
1334 * of theses labels is different from ours. (Their row vector
1335 * points along the row and not along the row dimension.)
1336 */
1337
1338 dircos_elid[VSLICE] = SPI_Image_normal;
1339 dircos_elid[VROW] = SPI_Image_column;
1340 dircos_elid[VCOLUMN] = SPI_Image_row;
1341
1342 /* Get direction cosines
1343 */
1344 for (ivolume = 0; ivolume < VOL_NDIMS; ivolume++) {
1345 element = acr_find_group_element(group_list, dircos_elid[ivolume]);
1346 if (element == NULL) {
1347 continue;
1348 }
1349 if (acr_get_element_numeric_array(element, WORLD_NDIMS,
1350 dircos[ivolume]) != WORLD_NDIMS) {
1351 continue;
1352 }
1353 /* negate the X and Z coordinates
1354 */
1355 convert_numa3_coordinate(dircos[ivolume]);
1356 found_dircos[ivolume] = TRUE;
1357 }
1358 }
1359 #endif
1360
1361 /* If we did not find the Siemens proprietary image vectors, try
1362 * the DICOM standard image position.
1363 */
1364 if (!found_dircos[VCOLUMN] || !found_dircos[VROW] || !found_dircos[VSLICE]) {
1365 if (dicom_read_orientation(group_list, RowColVec)) {
1366 dircos[VCOLUMN][XCOORD] = RowColVec[0];
1367 dircos[VCOLUMN][YCOORD] = RowColVec[1];
1368 dircos[VCOLUMN][ZCOORD] = RowColVec[2];
1369
1370 dircos[VROW][XCOORD] = RowColVec[3];
1371 dircos[VROW][YCOORD] = RowColVec[4];
1372 dircos[VROW][ZCOORD] = RowColVec[5];
1373
1374 found_dircos[VCOLUMN] = TRUE;
1375 found_dircos[VROW] = TRUE;
1376
1377 convert_dicom_coordinate(dircos[VROW]);
1378 convert_dicom_coordinate(dircos[VCOLUMN]);
1379
1380 /* slice direction unit vector is cross product of row,
1381 col vectors:
1382 */
1383 dircos[VSLICE][XCOORD] =
1384 dircos[VCOLUMN][YCOORD] * dircos[VROW][ZCOORD] -
1385 dircos[VCOLUMN][ZCOORD] * dircos[VROW][YCOORD];
1386
1387 dircos[VSLICE][YCOORD] =
1388 dircos[VCOLUMN][ZCOORD] * dircos[VROW][XCOORD] -
1389 dircos[VCOLUMN][XCOORD] * dircos[VROW][ZCOORD];
1390
1391 dircos[VSLICE][ZCOORD] =
1392 dircos[VCOLUMN][XCOORD] * dircos[VROW][YCOORD] -
1393 dircos[VCOLUMN][YCOORD] * dircos[VROW][XCOORD];
1394 found_dircos[VSLICE] = TRUE;
1395 }
1396 }
1397
1398 if (G.Debug >= HI_LOGGING) {
1399 printf("dircos %f %f %f %f %f %f %f %f %f\n",
1400 dircos[VSLICE][XCOORD],
1401 dircos[VSLICE][YCOORD],
1402 dircos[VSLICE][ZCOORD],
1403 dircos[VROW][XCOORD],
1404 dircos[VROW][YCOORD],
1405 dircos[VROW][ZCOORD],
1406 dircos[VCOLUMN][XCOORD],
1407 dircos[VCOLUMN][YCOORD],
1408 dircos[VCOLUMN][ZCOORD]);
1409 }
1410
1411 /* Normalize the direction cosines
1412 */
1413 for (ivolume = 0; ivolume < VOL_NDIMS; ivolume++) {
1414 magnitude = 0.0;
1415 for (iworld=0; iworld < WORLD_NDIMS; iworld++) {
1416 magnitude += dircos[ivolume][iworld] * dircos[ivolume][iworld];
1417 }
1418 if (magnitude <= 0) {
1419 found_dircos[ivolume] = FALSE;
1420 continue;
1421 }
1422 magnitude = sqrt(magnitude);
1423 for (iworld=0; iworld < WORLD_NDIMS; iworld++) {
1424 dircos[ivolume][iworld] /= magnitude;
1425 }
1426 }
1427
1428 /* If we don't find direction cosines, then assume transverse volume
1429 */
1430 if (!found_dircos[VSLICE] ||
1431 !found_dircos[VROW] ||
1432 !found_dircos[VCOLUMN]) {
1433
1434 if (G.Debug) {
1435 printf("Using default direction cosines\n");
1436 }
1437
1438 for (ivolume = 0; ivolume < VOL_NDIMS; ivolume++) {
1439 for (iworld = 0; iworld < WORLD_NDIMS; iworld++) {
1440 dircos[ivolume][iworld] =
1441 ((ivolume == (WORLD_NDIMS-iworld-1)) ? -1.0 : 0.0);
1442 found_dircos[iworld] = TRUE;
1443 }
1444 }
1445 }
1446
1447 /* Figure out volume index to world index mapping and sign of direction
1448 * cosines - the code below figures out the primary direction in x,y,z
1449 * of each volume coordinate (row,col,slice)
1450 */
1451 for (ivolume = 0; ivolume < VOL_NDIMS; ivolume++) {
1452 largest = -1.0;
1453 for (iworld = 0; iworld < WORLD_NDIMS; iworld++) {
1454 magnitude = dircos[ivolume][iworld];
1455 if (magnitude < 0.0) magnitude = -magnitude;
1456 if (magnitude > largest) {
1457 largest = magnitude;
1458 volume_to_world[ivolume] = iworld;
1459 }
1460 }
1461 }
1462
1463 if (G.Debug >= HI_LOGGING) {
1464 printf(" Volume_to_world slice=%s row=%s column=%s\n",
1465 World_Names[volume_to_world[VSLICE]],
1466 World_Names[volume_to_world[VROW]],
1467 World_Names[volume_to_world[VCOLUMN]]);
1468 }
1469
1470 /* Get orientation (depends on primary direction of slice normal)
1471 */
1472 *orientation = orientation_list[volume_to_world[VSLICE]];
1473 if (G.Debug >= HI_LOGGING) {
1474 printf(" Orientation is %s\n",
1475 (*orientation == SAGITTAL) ? "SAGITTAL" :
1476 (*orientation == CORONAL) ? "CORONAL" : "TRANSVERSE");
1477 }
1478
1479 /* Get step information
1480 */
1481
1482 dicom_read_pixel_size(group_list, psize);
1483
1484 steps[VCOLUMN] = psize[0];
1485 steps[VROW] = psize[1]; /* anisotropic resolution? */
1486
1487 /* Figure out the slice thickness. It could be from either one of
1488 * two possible places in the file.
1489 *
1490 * This code has changed several times, and there may be no single
1491 * correct way of deriving the true slice spacing from the official
1492 * DICOM slice thickness and slice spacing fields. My best guess is
1493 * to look for both fields, and to adopt the
1494 */
1495 slice_thickness = (double)acr_find_double(group_list, ACR_Slice_thickness, 0);
1496 slice_spacing = (double)acr_find_double(group_list, ACR_Spacing_between_slices, 0);
1497
1498 if (slice_thickness == 0.0) {
1499 /* No slice thickness value found. */
1500 if (slice_spacing == 0.0) {
1501 if (G.Debug >= HI_LOGGING) {
1502 printf("Using default slice thickness of 1.0\n");
1503 }
1504 steps[VSLICE] = 1.0;
1505 }
1506 else {
1507 if (G.Debug >= HI_LOGGING) {
1508 printf("Using (0018,0088) for slice thickness\n");
1509 }
1510 steps[VSLICE] = slice_spacing;
1511 }
1512 }
1513 else if (slice_spacing == 0.0) {
1514 /* No slice spacing value found. */
1515 if (G.Debug >= HI_LOGGING) {
1516 printf("Using (0018,0050) for slice thickness\n");
1517 }
1518 steps[VSLICE] = slice_thickness;
1519 }
1520 else {
1521 /* Both fields are set. I choose the slice spacing rather
1522 * than the slice thickness in this case. However, I believe
1523 * there is some evidence that this can cause problems in rare
1524 * cases.
1525 */
1526 if (G.Debug && !NEARLY_EQUAL(slice_thickness, slice_spacing)) {
1527 printf("WARNING: slice thickness conflict: ");
1528 printf("old = %.10f, new = %.10f\n",
1529 slice_thickness, slice_spacing);
1530 }
1531 steps[VSLICE] = slice_spacing;
1532 }
1533
1534 /* Make sure that direction cosines point the right way (dot
1535 * product of direction cosine and axis is positive) and that step
1536 * has proper sign.
1537 */
1538 for (ivolume = 0; ivolume < VOL_NDIMS; ivolume++) {
1539 iworld = volume_to_world[ivolume];
1540 if (dircos[ivolume][iworld] < 0.0) {
1541 if (G.Debug >= HI_LOGGING) {
1542 printf("Swapping direction of %s %s\n",
1543 Volume_Names[ivolume],
1544 World_Names[iworld]);
1545 }
1546 steps[ivolume] *= -1.0;
1547 for (iworld = 0; iworld < WORLD_NDIMS; iworld++) {
1548 dircos[ivolume][iworld] *= -1.0;
1549 }
1550 }
1551 }
1552
1553 /* Find 3D coordinate of slice - ACR_Image_position_patient gives
1554 * the *corner* of the slice!
1555 *
1556 * Start by assuming that we didn't find it.
1557 */
1558 found_coordinate = FALSE;
1559 for (iworld = 0; iworld < WORLD_NDIMS; iworld++) {
1560 coordinate[iworld] = 0.0;
1561 }
1562
1563 if (G.opts & OPTS_NO_LOCATION) {
1564 /* If the coordinates are untrustworthy, just generate something
1565 * reasonable for the slice coordinate. Ignore the rest.
1566 */
1567 coordinate[volume_to_world[VSLICE]] =
1568 (steps[VSLICE] * fi_ptr->index[SLICE]);
1569 found_coordinate = TRUE;
1570 }
1571 else {
1572 found_coordinate = dicom_read_position(group_list,
1573 fi_ptr->index[SLICE] - 1,
1574 coordinate);
1575 if (!found_coordinate) {
1576 /* Last gasp - try to interpret the slice location as our slice
1577 * position. It might work.
1578 */
1579 if (!found_coordinate) {
1580 coordinate[volume_to_world[VSLICE]] =
1581 (double)acr_find_double(group_list, ACR_Slice_location, 1.0);
1582 }
1583
1584 found_coordinate = TRUE;
1585 }
1586 }
1587
1588 convert_dicom_coordinate(coordinate);
1589
1590 /* Work out start positions in volume coordinates
1591 */
1592 for (ivolume=0; ivolume < VOL_NDIMS; ivolume++) {
1593
1594 if (found_coordinate &&
1595 found_dircos[VSLICE] &&
1596 found_dircos[VROW] &&
1597 found_dircos[VCOLUMN]) {
1598 starts[ivolume] =
1599 coordinate[XCOORD] * dircos[ivolume][XCOORD] +
1600 coordinate[YCOORD] * dircos[ivolume][YCOORD] +
1601 coordinate[ZCOORD] * dircos[ivolume][ZCOORD];
1602 }
1603 else {
1604 starts[ivolume] = 0.0;
1605 }
1606 }
1607
1608 if (G.Debug >= HI_LOGGING) {
1609 printf(" coordinate %f %f %f, start %f %f %f\n",
1610 coordinate[XCOORD], coordinate[YCOORD], coordinate[ZCOORD],
1611 starts[VROW], starts[VCOLUMN], starts[VSLICE]);
1612 }
1613
1614 /* Find position along each dimension
1615 */
1616 fi_ptr->coordinate[SLICE] = starts[VSLICE];
1617 fi_ptr->coordinate[ECHO] =
1618 (double)acr_find_double(group_list, ACR_Echo_time, 0.0) / MS_PER_SECOND;
1619
1620
1621 /* Get the dimension width for time, if available. The units are in
1622 * milliseconds in DICOM, whereas we use seconds in MINC.
1623 */
1624 fi_ptr->width[TIME] = (double)acr_find_double(group_list,
1625 ACR_Actual_frame_duration,
1626 0.0) / MS_PER_SECOND;
1627
1628 /* PET scan times (bert)
1629 */
1630 start_time = (double)acr_find_double(group_list, ACR_Frame_reference_time, -1.0);
1631 frame_time = (double)acr_find_double(group_list, ACR_Actual_frame_duration, -1.0);
1632 if (start_time > 0.0 && frame_time > 0.0) {
1633 frame_time = start_time / 1000.0; /* Convert msec to seconds. */
1634 }
1635 else {
1636 /* time section (rhoge)
1637 * now assume that time has been fixed when file was read
1638 */
1639 start_time = (double)acr_find_double(group_list, ACR_Series_time, 0.0);
1640 frame_time = (double)acr_find_double(group_list, ACR_Acquisition_time, 0.0);
1641 start_time = convert_time_to_seconds(start_time);
1642 frame_time = convert_time_to_seconds(frame_time) - start_time;
1643
1644 /* check for case where scan starts right before midnight,
1645 * but frame is after midnight
1646 */
1647 if (frame_time < 0.0) {
1648 frame_time += SECONDS_PER_DAY;
1649 }
1650 }
1651 fi_ptr->coordinate[TIME] = frame_time;
1652
1653 /* end of time section */
1654
1655 fi_ptr->coordinate[PHASE] = 0.0;
1656 fi_ptr->coordinate[CHEM_SHIFT] = 0.0;
1657
1658 }
1659
1660 #if 0
1661 /* ----------------------------- MNI Header -----------------------------------
1662 @NAME : convert_numa3_coordinate
1663 @INPUT : coordinate
1664 @OUTPUT : coordinate
1665 @RETURNS : (nothing)
1666 @DESCRIPTION: Routine to convert a coordinate to the correct orientation
1667 @METHOD :
1668 @GLOBALS :
1669 @CALLS :
1670 @CREATED : February 28, 1997 (Peter Neelin)
1671 @MODIFIED : made version specific to Numaris 3 SPI data (rhoge)
1672 ---------------------------------------------------------------------------- */
1673 static void
1674 convert_numa3_coordinate(double coordinate[WORLD_NDIMS])
1675 {
1676 coordinate[XCOORD] = -coordinate[XCOORD];
1677 coordinate[ZCOORD] = -coordinate[ZCOORD];
1678 }
1679 #endif
1680
1681 /* ----------------------------- MNI Header -----------------------------------
1682 @NAME : convert_dicom_coordinate
1683 @INPUT : coordinate
1684 @OUTPUT : coordinate
1685 @RETURNS : (nothing)
1686 @DESCRIPTION: Routine to convert a coordinate to the correct orientation
1687 @METHOD :
1688 @GLOBALS :
1689 @CALLS :
1690 @CREATED : February 28, 1997 (Peter Neelin)
1691 @MODIFIED : made new dicom version (rhoge)
1692 ---------------------------------------------------------------------------- */
1693 void
convert_dicom_coordinate(double coordinate[WORLD_NDIMS])1694 convert_dicom_coordinate(double coordinate[WORLD_NDIMS])
1695 {
1696 /* Allow the user to override this, if only for debugging purposes...
1697 */
1698 if (G.opts & OPTS_KEEP_COORD) {
1699 return;
1700 }
1701
1702 coordinate[XCOORD] = -coordinate[XCOORD];
1703 coordinate[YCOORD] = -coordinate[YCOORD];
1704 }
1705
1706 /* ----------------------------- MNI Header -----------------------------------
1707 @NAME : get_general_header_info
1708 @INPUT : group_list - input data
1709 @OUTPUT : gi_ptr - general information about files
1710 @RETURNS : (nothing)
1711 @DESCRIPTION: Routine to extract general header information from a group list
1712 @METHOD :
1713 @GLOBALS :
1714 @CALLS :
1715 @CREATED : February 28, 1997 (Peter Neelin)
1716 @MODIFIED :
1717 ---------------------------------------------------------------------------- */
1718
1719 static void
get_string_field(char * out_str,Acr_Group group_list,Acr_Element_Id element_id)1720 get_string_field(char *out_str, Acr_Group group_list,
1721 Acr_Element_Id element_id)
1722 {
1723 strncpy(out_str, acr_find_string(group_list, element_id, ""),
1724 STRING_T_LEN);
1725 }
1726
1727 void
get_general_header_info(Acr_Group group_list,General_Info * gi_ptr)1728 get_general_header_info(Acr_Group group_list, General_Info *gi_ptr)
1729 {
1730 int length;
1731 char *string;
1732
1733 if (G.Debug) {
1734 printf("SOP Class UID: %s\n",
1735 acr_find_string(group_list, ACR_SOP_Class_UID, ""));
1736 printf("Images in acquisition: %d\n",
1737 acr_find_int(group_list, ACR_Images_in_acquisition, -1));
1738 printf("Acquisitions in series: %d\n",
1739 acr_find_int(group_list, ACR_Acquisitions_in_series, -1));
1740 printf("3D raw partitions: %d\n",
1741 acr_find_int(group_list, SPI_Number_of_3D_raw_partitions_nominal, -1));
1742 }
1743 /* Get intensity units */
1744 strncpy(gi_ptr->units, "", STRING_T_LEN);
1745
1746 /* Get patient info */
1747 get_string_field(gi_ptr->patient.name, group_list, ACR_Patient_name);
1748
1749 get_string_field(gi_ptr->patient.identification,
1750 group_list, ACR_Patient_identification);
1751 get_string_field(gi_ptr->patient.birth_date,
1752 group_list, ACR_Patient_birth_date);
1753
1754 get_string_field(gi_ptr->patient.age,
1755 group_list, ACR_Patient_age);
1756
1757 string = (char*)acr_find_string(group_list, ACR_Patient_sex, "");
1758 if (*string == 'M')
1759 strncpy(gi_ptr->patient.sex, MI_MALE, STRING_T_LEN);
1760 else if (*string == 'F')
1761 strncpy(gi_ptr->patient.sex, MI_FEMALE, STRING_T_LEN);
1762 else if (*string == 'O')
1763 strncpy(gi_ptr->patient.sex, MI_OTHER, STRING_T_LEN);
1764 else
1765 strncpy(gi_ptr->patient.sex, "", STRING_T_LEN);
1766
1767 gi_ptr->patient.weight =
1768 (double)acr_find_double(group_list, ACR_Patient_weight, -DBL_MAX);
1769
1770 /* added by rhoge - registration timing info */
1771 get_string_field(gi_ptr->patient.reg_date,
1772 group_list, ACR_Study_date);
1773
1774 get_string_field(gi_ptr->patient.reg_time,
1775 group_list, ACR_Study_time);
1776
1777 get_string_field(gi_ptr->patient.position, /*position of patient added by ilana*/
1778 group_list, ACR_Patient_position);
1779
1780 /* Get study info */
1781
1782 /*some more timing info added by ilana*/
1783 get_string_field(gi_ptr->acq.series_time, group_list, ACR_Series_time);
1784
1785 get_string_field(gi_ptr->study.start_time,
1786 group_list, ACR_Study_date);
1787
1788 length = strlen(gi_ptr->study.start_time);
1789 gi_ptr->study.start_time[length] = ' ';
1790 length++;
1791 strncpy(&gi_ptr->study.start_time[length],
1792 acr_find_string(group_list, ACR_Study_time, ""), STRING_T_LEN - length);
1793 string = (char*)acr_find_string(group_list, ACR_Modality, "");
1794 if (strcmp(string, ACR_MODALITY_MR) == 0)
1795 strncpy(gi_ptr->study.modality, MI_MRI, STRING_T_LEN);
1796 else if (strcmp(string, ACR_MODALITY_PT) == 0)
1797 strncpy(gi_ptr->study.modality, MI_PET, STRING_T_LEN);
1798 get_string_field(gi_ptr->study.manufacturer,
1799 group_list, ACR_Manufacturer);
1800 get_string_field(gi_ptr->study.model,
1801 group_list, ACR_Manufacturer_model);
1802 gi_ptr->study.field_value =
1803 (double)acr_find_double(group_list, ACR_Magnetic_field_strength, -DBL_MAX);
1804 get_string_field(gi_ptr->study.software_version,
1805 group_list, ACR_Software_versions);
1806 get_string_field(gi_ptr->study.serial_no,
1807 group_list, ACR_Device_serial_number);
1808 get_string_field(gi_ptr->study.calibration_date,
1809 group_list, ACR_Calibration_date);
1810 get_string_field(gi_ptr->study.calibration_time, /*add time as well ilana*/
1811 group_list, ACR_Calibration_time);
1812 get_string_field(gi_ptr->study.institution,
1813 group_list, ACR_Institution_id);
1814 get_string_field(gi_ptr->study.station_id,
1815 group_list, ACR_Station_id);
1816 get_string_field(gi_ptr->study.referring_physician,
1817 group_list, ACR_Referring_physician);
1818 get_string_field(gi_ptr->study.performing_physician,
1819 group_list, ACR_Performing_physician);
1820 get_string_field(gi_ptr->study.operator,
1821 group_list, ACR_Operators_name);
1822 get_string_field(gi_ptr->study.procedure,
1823 group_list, ACR_Procedure_description);
1824 sprintf(gi_ptr->study.study_id, "%.6f",gi_ptr->study_id);
1825
1826 /* Acquisition id modified by rhoge to get rid of first digit that
1827 is not required for identification of run */
1828 /* sprintf(gi_ptr->study.acquisition_id, "%d_%d",
1829 acr_find_int(group_list, ACR_Series, 0), gi_ptr->acq_id); */
1830 sprintf(gi_ptr->study.acquisition_id, "%d", gi_ptr->acq_id);
1831
1832
1833
1834 /* Get acquisition information */
1835
1836 get_string_field(gi_ptr->acq.acquisition_time, group_list, ACR_Acquisition_time); /*add acquisition start time ilana*/
1837 get_string_field(gi_ptr->acq.image_time, group_list, ACR_Image_time);
1838 get_string_field(gi_ptr->acq.scan_seq, group_list, ACR_Sequence_name);
1839 get_string_field(gi_ptr->acq.protocol_name, group_list, ACR_Protocol_name);
1840 get_string_field(gi_ptr->acq.series_description, /*add series description ilana*/
1841 group_list, ACR_Series_description);
1842 get_string_field(gi_ptr->acq.receive_coil, group_list,
1843 ACR_Receive_coil_name);
1844 get_string_field(gi_ptr->acq.transmit_coil, group_list,
1845 ACR_Transmit_coil_name);
1846 get_string_field(gi_ptr->acq.slice_order, group_list,
1847 EXT_Slice_order);
1848 /*0x1 means ASCENDING
1849 0x2 means DESCENDING
1850 0x4 means INTERLEAVED*/
1851 if(!strcmp(gi_ptr->acq.slice_order,"0x1 "))
1852 strncpy(gi_ptr->acq.slice_order, "ascending", STRING_T_LEN);
1853 else if(!strcmp(gi_ptr->acq.slice_order,"0x2 "))
1854 strncpy(gi_ptr->acq.slice_order, "descending", STRING_T_LEN);
1855 else if(!strcmp(gi_ptr->acq.slice_order,"0x4 "))
1856 strncpy(gi_ptr->acq.slice_order, "interleaved", STRING_T_LEN);
1857
1858 gi_ptr->acq.rep_time =
1859 (double)acr_find_double(group_list, ACR_Repetition_time, -DBL_MAX);
1860 if (gi_ptr->acq.rep_time != -DBL_MAX)
1861 gi_ptr->acq.rep_time /= 1000.0;
1862
1863 gi_ptr->acq.echo_time =
1864 (double)acr_find_double(group_list, ACR_Echo_time, -DBL_MAX);
1865 if (gi_ptr->acq.echo_time != -DBL_MAX)
1866 gi_ptr->acq.echo_time /= 1000.0;
1867
1868 gi_ptr->acq.echo_train_length =
1869 (double)acr_find_double(group_list, ACR_Echo_train_length, -DBL_MAX); /*added echo train length ilana*/
1870
1871 gi_ptr->acq.echo_number =
1872 (double)acr_find_double(group_list, ACR_Echo_number, -DBL_MAX);
1873
1874 gi_ptr->acq.inv_time =
1875 (double)acr_find_double(group_list, ACR_Inversion_time, -DBL_MAX);
1876 if (gi_ptr->acq.inv_time != -DBL_MAX)
1877 gi_ptr->acq.inv_time /= 1000.0;
1878 gi_ptr->acq.delay_in_TR =
1879 (double)acr_find_double(group_list, EXT_Delay_in_TR, -DBL_MAX); /*added delay in TR ilana*/
1880 if (gi_ptr->acq.delay_in_TR != -DBL_MAX)
1881 gi_ptr->acq.delay_in_TR /= 1000000.0; /*write in seconds*/
1882 gi_ptr->acq.b_value =
1883 (double)acr_find_double(group_list, EXT_Diffusion_b_value, -DBL_MAX);
1884 gi_ptr->acq.flip_angle =
1885 (double)acr_find_double(group_list, ACR_Flip_angle, -DBL_MAX);
1886 gi_ptr->acq.slice_thickness =
1887 (double)acr_find_double(group_list, ACR_Slice_thickness, -DBL_MAX);
1888 gi_ptr->acq.num_slices =
1889 (double)acr_find_double(group_list, ACR_Images_in_acquisition, -DBL_MAX);
1890 gi_ptr->acq.num_dyn_scans =
1891 (double)acr_find_double(group_list, ACR_Acquisitions_in_series, -DBL_MAX);
1892 gi_ptr->acq.num_avg =
1893 (double)acr_find_double(group_list, ACR_Nr_of_averages, -DBL_MAX);
1894 gi_ptr->acq.imaging_freq =
1895 (double)acr_find_double(group_list, ACR_Imaging_frequency, -DBL_MAX);
1896 if (gi_ptr->acq.imaging_freq != -DBL_MAX)
1897 gi_ptr->acq.imaging_freq *= 1e6;
1898 get_string_field(gi_ptr->acq.imaged_nucl,
1899 group_list, ACR_Imaged_nucleus);
1900 gi_ptr->acq.win_center =
1901 (double)acr_find_double(group_list, ACR_Window_centre, -DBL_MAX);
1902 gi_ptr->acq.win_width =
1903 (double)acr_find_double(group_list, ACR_Window_width, -DBL_MAX);
1904
1905 gi_ptr->acq.num_phase_enc_steps =
1906 (double)acr_find_double(group_list, ACR_Number_of_phase_encoding_steps, -DBL_MAX);
1907 gi_ptr->acq.percent_sampling =
1908 (double)acr_find_double(group_list, ACR_Percent_sampling, -DBL_MAX); /*don't need to multiply by 100 ilana*/
1909
1910 gi_ptr->acq.percent_phase_fov =
1911 (double)acr_find_double(group_list, ACR_Percent_phase_field_of_view, -DBL_MAX); /*don't need to multiply by 100 ilana*/
1912
1913 gi_ptr->acq.pixel_bandwidth =
1914 (double)acr_find_double(group_list, ACR_Pixel_bandwidth, -DBL_MAX);
1915
1916 gi_ptr->acq.sar = (double)acr_find_double(group_list, ACR_SAR, -DBL_MAX);
1917
1918 get_string_field(gi_ptr->acq.mr_acq_type,
1919 group_list, ACR_MR_acquisition_type);
1920
1921 get_string_field(gi_ptr->acq.image_type, group_list, ACR_Image_type);
1922 if (G.Debug) {
1923 if (strstr(gi_ptr->acq.image_type, "MOSAIC") != NULL) {
1924 printf("This appears to be a Mosaic image\n");
1925 }
1926 }
1927
1928 get_string_field(gi_ptr->acq.phase_enc_dir,
1929 group_list, ACR_Phase_encoding_direction);
1930
1931 /*Add image comments*/
1932 /*strncpy(gi_ptr->acq.comments, "", STRING_T_LEN);*/
1933 get_string_field(gi_ptr->acq.comments,
1934 group_list, ACR_Image_comments);
1935
1936 /* Siemens Numaris 4 specific!
1937 */
1938
1939 #if 0
1940 gi_ptr->acq.MrProt = strdup(acr_find_string(group_list, EXT_MrProt_dump,
1941 ""));
1942 #else
1943 gi_ptr->acq.MrProt = strdup("");
1944 #endif
1945
1946 string = (char*)acr_find_string(group_list, ACR_Acquisition_contrast, "");
1947 gi_ptr->acq.dti = (strstr(string, "DIFFUSION") != NULL);
1948 }
1949
1950 /* ----------------------------- MNI Header -----------------------------------
1951 @NAME : convert_time_to_seconds
1952 @INPUT : dicom_time
1953 @OUTPUT : (none)
1954 @RETURNS : real time in seconds from beginning of day
1955 @DESCRIPTION: Routine to convert dicom seconds (decimal hhmmss.xxxxx) to
1956 real seconds since the start of day.
1957 @METHOD :
1958 @GLOBALS :
1959 @CALLS :
1960 @CREATED : February 28, 1997 (Peter Neelin)
1961 @MODIFIED :
1962 ---------------------------------------------------------------------------- */
1963 static double
convert_time_to_seconds(double dicom_time)1964 convert_time_to_seconds(double dicom_time)
1965 {
1966 /* Constants */
1967 #define DICOM_SECONDS_PER_HOUR 10000.0
1968 #define DICOM_SECONDS_PER_MINUTE 100.0
1969
1970 /* Variables */
1971 double hh, mm, ss;
1972
1973 /* Get the components of the time */
1974
1975 hh = (int) (dicom_time / DICOM_SECONDS_PER_HOUR);
1976 dicom_time -= hh * DICOM_SECONDS_PER_HOUR;
1977 mm = (int) (dicom_time / DICOM_SECONDS_PER_MINUTE);
1978 dicom_time -= mm * DICOM_SECONDS_PER_MINUTE;
1979 ss = dicom_time;
1980
1981 /* Work out the number of seconds */
1982
1983 return (hh * 3600.0) + (mm * 60.0) + ss;
1984 }
1985
1986 /* ----------------------------- MNI Header -----------------------------------
1987 @NAME : get_dicom_image_data
1988 @INPUT : group_list - input data
1989 @OUTPUT : image - image data structure (user must free data)
1990 @RETURNS : (nothing)
1991 @DESCRIPTION: Routine to get an image from a group list
1992 @METHOD :
1993 @GLOBALS :
1994 @CALLS :
1995 @CREATED : November 25, 1993 (Peter Neelin)
1996 @MODIFIED :
1997 ---------------------------------------------------------------------------- */
1998 void
get_dicom_image_data(Acr_Group group_list,Image_Data * image)1999 get_dicom_image_data(Acr_Group group_list, Image_Data *image)
2000 {
2001
2002 /* Variables */
2003 Acr_Element element;
2004 int nrows, ncolumns;
2005 int bits_alloc;
2006 int image_group;
2007 void *data = NULL;
2008 long imagepix, ipix;
2009 struct Acr_Element_Id elid;
2010 nc_type datatype;
2011
2012 /* Get the image information */
2013 bits_alloc = (int)acr_find_short(group_list, ACR_Bits_allocated, 0);
2014 nrows = (int)acr_find_short(group_list, ACR_Rows, 0);
2015 ncolumns = (int)acr_find_short(group_list, ACR_Columns, 0);
2016 image_group = (int)acr_find_short(group_list, ACR_Image_location, ACR_IMAGE_GID);
2017
2018 /* Figure out type */
2019 if (bits_alloc > CHAR_BIT)
2020 datatype = NC_SHORT;
2021 else
2022 datatype = NC_BYTE;
2023
2024 /* Set image info */
2025 imagepix = nrows * ncolumns;
2026 image->data = (unsigned short *) malloc(imagepix * sizeof(short));
2027 CHKMEM(image->data);
2028
2029 /* Get image pointer */
2030 elid.group_id = image_group;
2031 elid.element_id = ACR_IMAGE_EID;
2032 element = acr_find_group_element(group_list, &elid);
2033 if (element == NULL) {
2034 memset(image->data, 0, imagepix * sizeof(short));
2035 return;
2036 }
2037 data = acr_get_element_data(element);
2038
2039 /* Convert the data according to type */
2040
2041 /* Look for byte data */
2042 if (datatype == NC_BYTE) {
2043 for (ipix=0; ipix < imagepix; ipix++) {
2044 image->data[ipix] = *((unsigned char *) data + ipix);
2045 }
2046 }
2047 else {
2048
2049 /* Look for unpacked short data */
2050 if (bits_alloc == nctypelen(datatype) * CHAR_BIT) {
2051 acr_get_short(acr_get_element_byte_order(element),
2052 nrows*ncolumns, data, image->data);
2053 }
2054
2055 /* Fill with zeros in any other case */
2056 else {
2057 memset(image->data, 0, imagepix * sizeof(short));
2058 }
2059 }
2060
2061 return;
2062 }
2063
2064
2065 /* ----------------------------- MNI Header -----------------------------------
2066 @NAME : parse_dicom_groups
2067 @INPUT : group_list - list of acr-nema groups that make up object
2068 @OUTPUT : data_info - information about data object
2069 @RETURNS : (nothing)
2070 @DESCRIPTION: Routine to parse dicom object
2071 @METHOD :
2072 @GLOBALS :
2073 @CALLS :
2074 @CREATED : June 2001 (Rick Hoge)
2075 @MODIFIED :
2076 ---------------------------------------------------------------------------- */
2077
2078 #define IDEFAULT (-1)
2079
2080 void
parse_dicom_groups(Acr_Group group_list,Data_Object_Info * di_ptr)2081 parse_dicom_groups(Acr_Group group_list, Data_Object_Info *di_ptr)
2082 {
2083 Acr_Element element;
2084 Acr_Short AcqMat[4];
2085 Acr_Short freq_rows;
2086 Acr_Short freq_cols;
2087 Acr_Short phase_rows;
2088 Acr_Short phase_cols;
2089 double slice_coord[WORLD_NDIMS];
2090
2091 /* Get info to construct unique identifiers for study, series/acq
2092 * for file processing
2093 */
2094 get_identification_info(group_list,
2095 &(di_ptr->study_id), &(di_ptr->acq_id),
2096 &(di_ptr->rec_num), &(di_ptr->image_type));
2097
2098 /* Get number of echos, echo number, number of dynamic scans and
2099 * dynamic_scan_number
2100 */
2101
2102 di_ptr->num_echoes = acr_find_int(group_list,
2103 ACR_Echo_train_length,
2104 IDEFAULT);
2105
2106 di_ptr->echo_number = acr_find_int(group_list,
2107 ACR_Echo_number,
2108 IDEFAULT);
2109
2110 di_ptr->num_dyn_scans = acr_find_int(group_list,
2111 ACR_Acquisitions_in_series,
2112 IDEFAULT);
2113
2114 di_ptr->dyn_scan_number = acr_find_int(group_list,
2115 ACR_Acquisition,
2116 IDEFAULT);
2117
2118 di_ptr->global_image_number = acr_find_int(group_list,
2119 ACR_Image,
2120 IDEFAULT);
2121
2122 /* rhoge:
2123 new info added to di_ptr by rhoge: nominal number of slices;
2124 this is used in detection of a stream of files with the same
2125 acquisition ID number in which there are more files than
2126 slices. If the number of signal averages is greater than one,
2127 we will assume that this means the acquisition loop was used for
2128 dynamic scanning.
2129
2130 WARNINGS: the same thing may need to be done with `number of
2131 partitions' for it to work with 3D scans */
2132
2133 di_ptr->num_slices_nominal = acr_find_int(group_list,
2134 ACR_Images_in_acquisition,
2135 IDEFAULT);
2136
2137 di_ptr->slice_number = acr_find_int(group_list,
2138 SPI_Current_slice_number,
2139 IDEFAULT);
2140
2141 di_ptr->slice_location = (double)acr_find_double(group_list,
2142 ACR_Slice_location,
2143 0.0);
2144
2145 di_ptr->coord_found = dicom_read_position(group_list, 0, slice_coord);
2146
2147 /* identification info needed to generate unique session id
2148 * for file names
2149 */
2150 di_ptr->study_date = acr_find_int(group_list, ACR_Study_date,
2151 IDEFAULT);
2152
2153 di_ptr->study_time = acr_find_int(group_list, ACR_Study_time,
2154 IDEFAULT);
2155
2156 di_ptr->scanner_serialno = acr_find_int(group_list,
2157 ACR_Device_serial_number,
2158 IDEFAULT);
2159
2160 /* identification info needed to determine if mosaics used
2161 */
2162
2163 element = acr_find_group_element(group_list, ACR_Acquisition_matrix);
2164
2165 if (element != NULL) {
2166 acr_get_element_short_array(element, 4, AcqMat);
2167
2168 freq_rows = AcqMat[0];
2169 freq_cols = AcqMat[1];
2170
2171 phase_rows = AcqMat[2];
2172 phase_cols = AcqMat[3];
2173
2174 /* rows in acq matrix is larger of freq rows and freq columns:
2175 */
2176 di_ptr->acq_rows = ( freq_rows > freq_cols ? freq_rows : freq_cols );
2177
2178 /* all images are square, at this time
2179 */
2180 di_ptr->acq_cols = di_ptr->acq_rows;
2181 }
2182 else {
2183 di_ptr->acq_rows = IDEFAULT;
2184 di_ptr->acq_cols = IDEFAULT;
2185 }
2186
2187 di_ptr->rec_rows = acr_find_int(group_list, ACR_Rows, IDEFAULT);
2188 di_ptr->rec_cols = acr_find_int(group_list, ACR_Columns, IDEFAULT);
2189
2190 di_ptr->num_mosaic_rows = acr_find_int(group_list, EXT_Mosaic_rows,
2191 IDEFAULT);
2192 di_ptr->num_mosaic_cols = acr_find_int(group_list, EXT_Mosaic_columns,
2193 IDEFAULT);
2194 di_ptr->num_slices_in_file = acr_find_int(group_list, EXT_Slices_in_file,
2195 IDEFAULT);
2196
2197 /* sequence, protocol names (useful for debugging):
2198 */
2199
2200 get_string_field(di_ptr->sequence_name, group_list, ACR_Sequence_name);
2201 get_string_field(di_ptr->protocol_name, group_list, ACR_Protocol_name);
2202 get_string_field(di_ptr->patient_name, group_list, ACR_Patient_name);
2203 get_string_field(di_ptr->patient_id, group_list, ACR_Patient_identification);
2204 }
2205
2206
2207