1 /* ----------------------------- MNI Header -----------------------------------
2 @NAME : gyro_read.c
3 @DESCRIPTION: Code to read gyrocom files and get info from them.
4 @METHOD :
5 @GLOBALS :
6 @CALLS :
7 @CREATED : November 25, 1993 (Peter Neelin)
8 @MODIFIED :
9 * $Log: gyro_read.c,v $
10 * Revision 6.2 2001-04-09 23:02:49 neelin
11 * Modified copyright notice, removing permission statement since copying,
12 * etc. is probably not permitted by our non-disclosure agreement with
13 * Philips.
14 *
15 * Revision 6.1 1999/10/29 17:52:03 neelin
16 * Fixed Log keyword
17 *
18 * Revision 6.0 1997/09/12 13:23:50 neelin
19 * Release of minc version 0.6
20 *
21 * Revision 5.0 1997/08/21 13:24:50 neelin
22 * Release of minc version 0.5
23 *
24 * Revision 4.1 1997/06/13 22:08:13 neelin
25 * Modifications to get gcomserver working with modified Acr_nema library.
26 *
27 * Revision 4.0 1997/05/07 20:01:07 neelin
28 * Release of minc version 0.4
29 *
30 * Revision 3.1 1995/08/02 13:41:36 neelin
31 * Fixed bug in direction cosine inversion (in test cases, this code was never
32 * called, so it does not seem to be an important bug).
33 *
34 * Revision 3.0 1995/05/15 19:31:44 neelin
35 * Release of minc version 0.3
36 *
37 * Revision 2.7 1995/03/21 15:43:19 neelin
38 * Corrected setting of default image type.
39 *
40 * Revision 2.6 1995/02/14 18:12:26 neelin
41 * Added project names and defaults files (using volume name).
42 * Added process id to log file name.
43 * Moved temporary files to subdirectory.
44 *
45 * Revision 2.5 1995/02/08 19:31:47 neelin
46 * Moved ARGSUSED statements for irix 5 lint.
47 *
48 * Revision 2.4 1994/12/12 09:05:57 neelin
49 * Changed comment in calculate_slice_start (code is the same)
50 *
51 * Revision 2.3 94/12/07 08:20:22 neelin
52 * Fixed some lint messages.
53 *
54 * Revision 2.2 94/11/21 08:07:59 neelin
55 * Modified code to properly calculate start from centre locations, then
56 * changed calculation back to old way because it worked.
57 * Added a ncsetfill(mincid, NC_NOFILL).
58 *
59 * Revision 2.1 94/10/20 13:50:11 neelin
60 * Write out direction cosines to support rotated volumes.
61 * Store single slices as 1-slice volumes (3D instead of 2D).
62 * Changed storing of minc history (get args for gyrotominc).
63 *
64 * Revision 2.0 94/09/28 10:35:27 neelin
65 * Release of minc version 0.2
66 *
67 * Revision 1.10 94/09/28 10:34:59 neelin
68 * Pre-release
69 *
70 * Revision 1.9 94/05/24 15:06:35 neelin
71 * Break up multiple echoes or time frames into separate files for 2 echoes
72 * or 2 frames (put in 1 file for more).
73 * Changed units of repetition time, echo time, etc to seconds.
74 * Save echo times in dimension variable when appropriate.
75 * Changed to file names to end in _mri.mnc.
76 *
77 * Revision 1.8 94/03/15 14:24:01 neelin
78 * Changed image-max/min to use fp_scaled_max/min instead of ext_scale_max/min
79 * Added acquisition:comments attribute
80 * Changed reading of configuration file to allow execution of a command on
81 * the minc file.
82 *
83 * Revision 1.7 94/03/14 16:44:30 neelin
84 * Changed scale to be fp_scaled_min/max instead of ext_scale_min/max.
85 * Check units for millirad and change to radians.
86 *
87 * Revision 1.6 94/01/17 15:05:46 neelin
88 * Added some acquisition parameters (flip angle) and fixed error in writing
89 * of scanning sequence.
90 *
91 * Revision 1.5 94/01/14 11:37:18 neelin
92 * Fixed handling of multiple reconstructions and image types. Add spiinfo variable with extra info (including window min/max). Changed output
93 * file name to include reconstruction number and image type number.
94 *
95 * Revision 1.4 94/01/06 14:18:44 neelin
96 * Added image unpacking (3 bytes -> 2 pixels).
97 *
98 * Revision 1.3 93/12/14 16:36:49 neelin
99 * Fixed axis direction and start position.
100 *
101 * Revision 1.2 93/12/10 15:35:07 neelin
102 * Improved file name generation from patient name. No buffering on stderr.
103 * Added spi group list to minc header.
104 * Optionally read a defaults file to get output minc directory and owner.
105 *
106 * Revision 1.1 93/11/30 14:41:04 neelin
107 * Initial revision
108 *
109 @COPYRIGHT :
110 Copyright 1993 Peter Neelin, McConnell Brain Imaging Centre,
111 Montreal Neurological Institute, McGill University.
112 ---------------------------------------------------------------------------- */
113
114 #include <gcomserver.h>
115 #include <math.h>
116
117 /* Define strings for spi image types */
118 static struct {
119 int image_type;
120 char *name;
121 } Spi_image_type_strings[] = {
122 0, "SE_R",
123 1, "SE_I",
124 2, "SE_M",
125 3, "SE_P",
126 13, "SE/CR",
127 4, "SE_T2",
128 5, "SE_RHO",
129 6, "IR_R",
130 7, "IR_I",
131 8, "IR_M",
132 9, "IR_P",
133 14, "IR/CR",
134 10, "IR_T2",
135 11, "IR_RHO",
136 12, "T1",
137 15, "T2",
138 16, "RHO",
139 17, "FFE/CR",
140 18, "FFE/R",
141 19, "FFE/I",
142 20, "FFE/M",
143 21, "FFE/P",
144 22, "Spectra",
145 23, "PCA_R",
146 24, "PCA_I",
147 25, "PCA_M",
148 26, "PCA_P",
149 27, "Derived"
150 };
151
152 /* ----------------------------- MNI Header -----------------------------------
153 @NAME : read_gyro
154 @INPUT : filename - name of gyrocom file to read
155 max_group - maximum group number to read
156 @OUTPUT : (none)
157 @RETURNS : group list read in from file
158 @DESCRIPTION: Routine to read in a group list from a file.
159 @METHOD :
160 @GLOBALS :
161 @CALLS :
162 @CREATED : November 25, 1993 (Peter Neelin)
163 @MODIFIED :
164 ---------------------------------------------------------------------------- */
read_gyro(char * filename,int max_group)165 public Acr_Group read_gyro(char *filename, int max_group)
166 {
167 FILE *fp;
168 Acr_File *afp;
169 Acr_Group group_list;
170
171 /* Open the file */
172 fp = fopen(filename, "r");
173 if (fp == NULL) return NULL;
174
175 /* Connect to input stream */
176 afp=acr_file_initialize(fp, 0, acr_stdio_read);
177 (void) acr_test_byte_order(afp);
178
179 /* Read in group list */
180 (void) acr_input_group_list(afp, &group_list, max_group);
181
182 /* Close the file */
183 acr_file_free(afp);
184 (void) fclose(fp);
185
186 return group_list;
187
188 }
189
190 /* ----------------------------- MNI Header -----------------------------------
191 @NAME : get_file_info
192 @INPUT : group_list - input data
193 @OUTPUT : file_info - file-specific info
194 general_info - general information about files
195 @RETURNS : (nothing)
196 @DESCRIPTION: Routine to extract information from a group list
197 @METHOD :
198 @GLOBALS :
199 @CALLS :
200 @CREATED : November 25, 1993 (Peter Neelin)
201 @MODIFIED :
202 ---------------------------------------------------------------------------- */
get_file_info(Acr_Group group_list,File_Info * file_info,General_Info * general_info)203 public void get_file_info(Acr_Group group_list, File_Info *file_info,
204 General_Info *general_info)
205 {
206 Acr_Element element;
207 int data_set_type;
208 Mri_Index imri;
209 World_Index iworld, jworld;
210 char *string;
211 int nrows;
212 int ncolumns;
213 int study_id, acq_id, rec_num, image_type;
214 int image_group;
215 int length, maxlen;
216 int cur_index;
217 int index;
218 Orientation orientation;
219 double default_slice_pos;
220 double angulation_ap, angulation_lr, angulation_cc;
221 char *image_type_string;
222 char image_type_string_buf[20];
223 Acr_Element_Id mri_index_list[MRI_NDIMS];
224 Acr_Element_Id mri_total_list[MRI_NDIMS];
225 Acr_Element_Id off_center_elid[WORLD_NDIMS];
226
227 /* Directions of axes for different orientations */
228 static double axis_direction[NUM_ORIENTATIONS][WORLD_NDIMS] = {
229 -1.0, -1.0, +1.0,
230 +1.0, -1.0, -1.0,
231 -1.0, -1.0, -1.0};
232
233 /* Array of elements for mri dimensions */
234 mri_index_list[SLICE] = SPI_Slice_number;
235 mri_index_list[ECHO] = SPI_Echo_number;
236 mri_index_list[TIME] = SPI_Dynamic_scan_number;
237 mri_index_list[PHASE] = SPI_Phase_number;
238 mri_index_list[CHEM_SHIFT] = SPI_Chemical_shift_number;
239 mri_total_list[SLICE] = SPI_Number_of_slices;
240 mri_total_list[ECHO] = SPI_Number_of_echoes;
241 mri_total_list[TIME] = SPI_Number_of_dynamic_scans;
242 mri_total_list[PHASE] = SPI_Number_of_phases;
243 mri_total_list[CHEM_SHIFT] = SPI_Nr_of_chemical_shifts;
244
245 /* Array of elements giving off center for each axis */
246 off_center_elid[XCOORD] = SPI_Off_center_lr;
247 off_center_elid[YCOORD] = SPI_Off_center_ap;
248 off_center_elid[ZCOORD] = SPI_Off_center_cc;
249
250 /* Look for data set type */
251 element = acr_find_group_element(group_list, ACR_Data_set_type);
252 if (element != NULL)
253 data_set_type = acr_get_element_short(element);
254
255 if ((element == NULL) || (data_set_type != ACR_IMAGE_OBJECT)) {
256 file_info->valid = FALSE;
257 return;
258 }
259
260 /* Get stuff that must be in file */
261 file_info->bits_alloc = find_short(group_list, ACR_Bits_allocated, 0);
262 file_info->bits_stored = find_short(group_list, ACR_Bits_stored, 0);
263 nrows = find_short(group_list, ACR_Rows, 0);
264 ncolumns = find_short(group_list, ACR_Columns, 0);
265 if (nrows <= 0)
266 nrows = find_int(group_list, SPI_Recon_resolution, 0);
267 if (ncolumns <= 0)
268 ncolumns = find_int(group_list, SPI_Recon_resolution, 0);
269 image_group = find_short(group_list, ACR_Image_location, INT_MIN);
270
271 /* Check for necessary values not found */
272 if ((nrows <= 0) || (ncolumns <= 0) ||
273 (file_info->bits_stored <= 0) ||
274 (file_info->bits_alloc <= 0) ||
275 (image_group < 0)) {
276 file_info->valid = FALSE;
277 return;
278 }
279
280 /* Get study, acq, rec, image type id's */
281 element = acr_find_group_element(group_list, ACR_Study);
282 if (element != NULL)
283 study_id = (int) acr_get_element_numeric(element);
284 else
285 study_id = 0;
286 element = acr_find_group_element(group_list, ACR_Acquisition);
287 if (element != NULL)
288 acq_id = (int) acr_get_element_numeric(element);
289 else
290 acq_id = 0;
291 element = acr_find_group_element(group_list, SPI_Reconstruction_number);
292 if (element != NULL)
293 rec_num = (int) acr_get_element_numeric(element);
294 else
295 rec_num = 1;
296 element = acr_find_group_element(group_list, SPI_Image_type);
297 if (element != NULL)
298 image_type = (int) acr_get_element_numeric(element);
299 else
300 image_type = SPI_DEFAULT_IMAGE_TYPE;
301
302 /* Get pixel value information */
303 file_info->pixel_min = find_short(group_list, ACR_Smallest_pixel_value, 0);
304 file_info->pixel_max = find_short(group_list, ACR_Largest_pixel_value,
305 (1 << file_info->bits_stored) - 1);
306 file_info->slice_min = find_double(group_list, SPI_Fp_scaled_minimum,
307 (double) file_info->pixel_min);
308 file_info->slice_max = find_double(group_list, SPI_Fp_scaled_maximum,
309 (double) file_info->pixel_max);
310
311 /* Get window min and max */
312 file_info->window_min =
313 find_double(group_list, SPI_Fp_window_minimum, file_info->slice_min);
314 file_info->window_max =
315 find_double(group_list, SPI_Fp_window_maximum, file_info->slice_max);
316
317 /* Get image indices */
318 for (imri=0; imri < MRI_NDIMS; imri++) {
319 file_info->index[imri] =
320 find_int(group_list, mri_index_list[imri], 1) - 1;
321 if (file_info->index[imri] < 0)
322 file_info->index[imri] = 0;
323 }
324
325 /* Set up general info on first pass */
326 if (!general_info->initialized) {
327
328 /* Get row and columns sizes */
329 general_info->nrows = nrows;
330 general_info->ncolumns = ncolumns;
331
332 /* Save the study, acqusition, reconstruction and image type
333 identifiers */
334 general_info->study_id = study_id;
335 general_info->acq_id = acq_id;
336 general_info->rec_num = rec_num;
337 general_info->image_type = image_type;
338
339 /* Get a name for the image type (if not found, use the number) */
340 image_type_string = NULL;
341 for (index=0; index < ARRAY_SIZE(Spi_image_type_strings); index++) {
342 if (image_type == Spi_image_type_strings[index].image_type) {
343 image_type_string = Spi_image_type_strings[index].name;
344 break;
345 }
346 }
347 if (image_type_string == NULL) {
348 image_type_string = image_type_string_buf;
349 (void) sprintf(image_type_string, "%d", image_type);
350 }
351 (void) strncpy(general_info->image_type_string,
352 image_type_string, sizeof(Cstring) - 1);
353
354 /* Get dimension information */
355 for (imri=0; imri < MRI_NDIMS; imri++) {
356 general_info->size[imri] = 1;
357 general_info->first[imri] = file_info->index[imri];
358 general_info->total_size[imri] =
359 find_int(group_list, mri_total_list[imri], 1);
360 if (general_info->total_size[imri] < 1)
361 general_info->total_size[imri] = 1;
362 general_info->position[imri] = NULL;
363 general_info->image_index[imri] = -1;
364 }
365
366 /* Set direction cosines from rotation angles */
367 angulation_ap = find_double(group_list, SPI_Angulation_of_ap_axis, 0.0);
368 angulation_lr = find_double(group_list, SPI_Angulation_of_lr_axis, 0.0);
369 angulation_cc = find_double(group_list, SPI_Angulation_of_cc_axis, 0.0);
370 get_direction_cosines(angulation_ap, angulation_lr, angulation_cc,
371 general_info->dircos);
372
373 /* Get spatial orientation */
374 switch (find_int(group_list, SPI_Slice_orientation, 0)) {
375 case SPI_SAGITTAL_ORIENTATION:
376 orientation = SAGITTAL; break;
377 case SPI_CORONAL_ORIENTATION:
378 orientation = CORONAL; break;
379 case SPI_TRANSVERSE_ORIENTATION:
380 default:
381 orientation = TRANSVERSE; break;
382 }
383 get_orientation_info(orientation, general_info->dircos,
384 &general_info->slice_world,
385 &general_info->row_world,
386 &general_info->column_world);
387
388 /* Get step information. Use field-of-view since pixel size doesn't
389 seem to be correct */
390 #if 1
391 general_info->step[general_info->row_world] =
392 find_double(group_list, SPI_Field_of_view, (double) nrows)
393 / (double) nrows;
394 if (general_info->step[general_info->row_world] == 0.0) {
395 general_info->step[general_info->row_world] = 1.0;
396 }
397 general_info->step[general_info->column_world] =
398 general_info->step[general_info->row_world];
399 #else
400 string = find_string(group_list, ACR_Pixel_size, "");
401 general_info->step[general_info->column_world] = strtod(string, &end);
402 if (end == string) {
403 general_info->step[general_info->column_world] = 1.0;
404 }
405 else {
406 string = end;
407 if ((*string == '\\') || (*string == ',')) string++;
408 general_info->step[general_info->row_world] = strtod(string, &end);
409 if (end == string)
410 general_info->step[general_info->row_world] = 1.0;
411 }
412 #endif
413 general_info->step[general_info->slice_world] =
414 find_double(group_list, ACR_Slice_thickness, 1.0);
415 general_info->step[XCOORD] *= axis_direction[orientation][XCOORD];
416 general_info->step[YCOORD] *= axis_direction[orientation][YCOORD];
417 general_info->step[ZCOORD] *= axis_direction[orientation][ZCOORD];
418
419 /* Make sure that direction cosines point the right way (dot product
420 of direction cosine and axis is positive) and that step has proper
421 sign */
422 for (iworld = XCOORD; iworld < WORLD_NDIMS; iworld++) {
423 if (general_info->dircos[iworld][iworld] < 0.0) {
424 general_info->step[iworld] *= -1.0;
425 for (jworld = XCOORD; jworld < WORLD_NDIMS; jworld++) {
426 general_info->dircos[iworld][jworld] *= -1.0;
427 }
428 }
429 }
430
431 /* Get centre information (start info is calculated farther down) */
432 general_info->centre[XCOORD] =
433 find_double(group_list, SPI_Off_center_lr, 0.0);
434 general_info->centre[YCOORD] =
435 find_double(group_list, SPI_Off_center_ap, 0.0);
436 general_info->centre[ZCOORD] =
437 find_double(group_list, SPI_Off_center_cc, 0.0);
438 calculate_slice_start(general_info->slice_world,
439 general_info->row_world,
440 general_info->column_world,
441 general_info->nrows,
442 general_info->ncolumns,
443 general_info->centre,
444 general_info->step,
445 general_info->dircos,
446 general_info->start);
447
448 /* Keep track of range of slices */
449 general_info->slicepos_first =
450 general_info->centre[general_info->slice_world];
451 general_info->slicepos_last =
452 general_info->centre[general_info->slice_world];
453 general_info->sliceindex_first = file_info->index[SLICE];
454 general_info->sliceindex_last = file_info->index[SLICE];
455
456 /* Save slice step info */
457 iworld = general_info->slice_world;
458 general_info->slice_step = general_info->step[iworld];
459 general_info->slice_start = general_info->centre[iworld] -
460 file_info->index[SLICE] * general_info->slice_step;
461
462 /* Set data type and range */
463 if (file_info->bits_alloc <= 8)
464 general_info->datatype = NC_BYTE;
465 else
466 general_info->datatype = NC_SHORT;
467 general_info->is_signed = ((general_info->datatype == NC_SHORT) &&
468 (file_info->bits_stored < 16));
469 general_info->pixel_min = file_info->pixel_min;
470 general_info->pixel_max = file_info->pixel_max;
471
472 /* Save display window info */
473 general_info->window_min = file_info->window_min;
474 general_info->window_max = file_info->window_max;
475
476 /* Maximum length for strings */
477 maxlen = sizeof(Cstring) - 1;
478
479 /* Get intensity units */
480 (void) strncpy(general_info->units,
481 find_string(group_list, SPI_Ext_scale_units, ""), maxlen);
482 if (strcmp(general_info->units, "millirad") == 0) {
483 (void) strncpy(general_info->units, "radians", maxlen);
484 }
485
486 /* Get patient info */
487 (void) strncpy(general_info->patient.name,
488 find_string(group_list, ACR_Patient_name, ""), maxlen);
489 (void) strncpy(general_info->patient.identification,
490 find_string(group_list, ACR_Patient_identification, ""), maxlen);
491 (void) strncpy(general_info->patient.birth_date,
492 find_string(group_list, ACR_Patient_birth_date, ""), maxlen);
493 string = find_string(group_list, ACR_Patient_sex, "");
494 if (*string == 'M')
495 (void) strncpy(general_info->patient.sex, MI_MALE, maxlen);
496 else if (*string == 'F')
497 (void) strncpy(general_info->patient.sex, MI_FEMALE, maxlen);
498 else if (*string == 'O')
499 (void) strncpy(general_info->patient.sex, MI_OTHER, maxlen);
500 else
501 (void) strncpy(general_info->patient.sex, "", maxlen);
502 general_info->patient.weight =
503 find_double(group_list, ACR_Patient_weight, -DBL_MAX);
504
505 /* Get study info */
506 (void) strncpy(general_info->study.start_time,
507 find_string(group_list, ACR_Study_date, ""), maxlen - 1);
508 length = strlen(general_info->study.start_time);
509 general_info->study.start_time[length] = ' ';
510 length++;
511 (void) strncpy(&general_info->study.start_time[length],
512 find_string(group_list, ACR_Study_time, ""), maxlen - length);
513 string = find_string(group_list, ACR_Modality, "");
514 if (strcmp(string, ACR_MODALITY_MR) == 0)
515 (void) strncpy(general_info->study.modality, MI_MRI, maxlen);
516 (void) strncpy(general_info->study.manufacturer,
517 find_string(group_list, ACR_Manufacturer, ""), maxlen);
518 (void) strncpy(general_info->study.model,
519 find_string(group_list, ACR_Manufacturer_model, ""), maxlen);
520 (void) strncpy(general_info->study.institution,
521 find_string(group_list, ACR_Institution_id, ""), maxlen);
522 (void) strncpy(general_info->study.station_id,
523 find_string(group_list, ACR_Station_id, ""), maxlen);
524 (void) strncpy(general_info->study.ref_physician,
525 find_string(group_list, ACR_Referring_physician, ""), maxlen);
526 (void) strncpy(general_info->study.procedure,
527 find_string(group_list, ACR_Procedure_description, ""), maxlen);
528 (void) sprintf(general_info->study.study_id, "%d",
529 find_int(group_list, ACR_Study, 0));
530 (void) sprintf(general_info->study.acquisition_id, "%d",
531 find_int(group_list, ACR_Acquisition, 0));
532
533 /* Get acquisition information */
534 (void) strncpy(general_info->acq.scan_seq,
535 find_string(group_list, ACR_Scanning_sequence, ""), maxlen);
536 general_info->acq.rep_time =
537 find_double(group_list, ACR_Repetition_time, -DBL_MAX) / 1000.0;
538 general_info->acq.echo_time =
539 find_double(group_list, ACR_Echo_time, -DBL_MAX) / 1000.0;
540 general_info->acq.inv_time =
541 find_double(group_list, ACR_Inversion_time, -DBL_MAX) / 1000.0;
542 general_info->acq.flip_angle =
543 find_double(group_list, SPI_Flip_angle, -DBL_MAX);
544 general_info->acq.num_avg =
545 find_double(group_list, ACR_Nr_of_averages, -DBL_MAX);
546 general_info->acq.imaging_freq =
547 find_double(group_list, ACR_Imaging_frequency, -DBL_MAX) * 1000000.0;
548 (void) strncpy(general_info->acq.imaged_nucl,
549 find_string(group_list, ACR_Imaged_nucleus, ""), maxlen);
550 (void) strncpy(general_info->acq.comments,
551 find_string(group_list, ACR_Acq_comments, ""), maxlen);
552
553 /* Copy the group list */
554 general_info->group_list = acr_copy_group_list(group_list);
555
556 /* Set initialized flag */
557 general_info->initialized = TRUE;
558 } /* Set up file info */
559
560 /* Update general info and validate file on later passes */
561 else {
562
563 /* Check for consistent data type */
564 if (((general_info->datatype == NC_BYTE) &&
565 (file_info->bits_alloc > 8)) ||
566 ((general_info->datatype == NC_SHORT) &&
567 (file_info->bits_alloc <= 8))) {
568 file_info->valid = FALSE;
569 return;
570 }
571
572 /* Check row and columns sizes */
573 if ((nrows != general_info->nrows) &&
574 (ncolumns != general_info->ncolumns)) {
575 file_info->valid = FALSE;
576 return;
577 }
578
579 /* Check study, acq, reconstruction and image type id's */
580 if ((general_info->study_id != study_id) ||
581 (general_info->acq_id != acq_id) ||
582 (general_info->rec_num != rec_num) ||
583 (general_info->image_type != image_type)) {
584 file_info->valid = FALSE;
585 return;
586 }
587
588 /* Check indices for range */
589 for (imri=0; imri < MRI_NDIMS; imri++) {
590 cur_index = file_info->index[imri];
591 if ((cur_index < 0) ||
592 (cur_index >= general_info->total_size[imri])) {
593 file_info->valid = FALSE;
594 return;
595 }
596 }
597
598 /* Look to see if indices have changed */
599 for (imri=0; imri < MRI_NDIMS; imri++) {
600
601 /* Get current index */
602 cur_index = file_info->index[imri];
603
604 /* If we only have one index and this index is different, then
605 allocate a list */
606 if ((general_info->size[imri] == 1) &&
607 (cur_index != general_info->first[imri])) {
608 general_info->position[imri] =
609 MALLOC(general_info->total_size[imri] * sizeof(int));
610 for (index=0; index < general_info->total_size[imri]; index++)
611 general_info->position[imri][index] = -1;
612 general_info->position[imri][general_info->first[imri]] = 0;
613 general_info->position[imri][cur_index] = 0;
614 if (cur_index < general_info->first[imri])
615 general_info->first[imri] = cur_index;
616 general_info->size[imri]++;
617 }
618
619 /* If we have more than one index and this one has not been found
620 yet, then set it to true */
621 else if ((general_info->size[imri] > 1) &&
622 (general_info->position[imri][cur_index] < 0)) {
623 general_info->position[imri][cur_index] = 0;
624 if (cur_index < general_info->first[imri])
625 general_info->first[imri] = cur_index;
626 general_info->size[imri]++;
627 }
628
629 /* Update position list */
630 if (general_info->position[imri] != NULL) {
631 cur_index = 0;
632 for (index=0; index < general_info->total_size[imri]; index++) {
633 if (general_info->position[imri][index] >= 0) {
634 general_info->position[imri][index] = cur_index;
635 cur_index++;
636 }
637 }
638 }
639
640 } /* Loop over Mri_Index */
641
642 /* Update display window info */
643 if (general_info->window_min > file_info->window_min)
644 general_info->window_min = file_info->window_min;
645 if (general_info->window_max < file_info->window_max)
646 general_info->window_max = file_info->window_max;
647
648 } /* Update general info for this file */
649
650 /* Get other slice specific information */
651 default_slice_pos = general_info->slice_start +
652 file_info->index[SLICE] * general_info->slice_step;
653 file_info->slice_position =
654 find_double(group_list, off_center_elid[general_info->slice_world],
655 default_slice_pos);
656 file_info->dyn_begin_time =
657 find_double(group_list, SPI_dynamic_scan_begin_time, 0.0);
658 file_info->echo_time =
659 find_double(group_list, ACR_Echo_time, 0.0) / 1000.0;
660
661 /* Update slice position information for general_info. If we have a new
662 first slice, then update centre and start information */
663 if (file_info->index[SLICE] < general_info->sliceindex_first) {
664 general_info->sliceindex_first = file_info->index[SLICE];
665 general_info->slicepos_first = file_info->slice_position;
666
667 /* Update slice centre and start info */
668 general_info->centre[general_info->slice_world] =
669 general_info->slicepos_first;
670 calculate_slice_start(general_info->slice_world,
671 general_info->row_world,
672 general_info->column_world,
673 general_info->nrows,
674 general_info->ncolumns,
675 general_info->centre,
676 general_info->step,
677 general_info->dircos,
678 general_info->start);
679 }
680 if (file_info->index[SLICE] > general_info->sliceindex_last) {
681 general_info->sliceindex_last = file_info->index[SLICE];
682 general_info->slicepos_last = file_info->slice_position;
683 }
684 if (general_info->size[SLICE] > 1) {
685 general_info->step[general_info->slice_world] =
686 (general_info->slicepos_last - general_info->slicepos_first) /
687 (general_info->size[SLICE] - 1);
688 }
689
690 /* If we get to here, then we have a valid file */
691 file_info->valid = TRUE;
692
693 return;
694 }
695
696 /* ----------------------------- MNI Header -----------------------------------
697 @NAME : get_gyro_image
698 @INPUT : group_list - input data
699 @OUTPUT : image - image data structure (user must free data)
700 @RETURNS : (nothing)
701 @DESCRIPTION: Routine to get an image from a group list
702 @METHOD :
703 @GLOBALS :
704 @CALLS :
705 @CREATED : November 25, 1993 (Peter Neelin)
706 @MODIFIED :
707 ---------------------------------------------------------------------------- */
get_gyro_image(Acr_Group group_list,Image_Data * image)708 public void get_gyro_image(Acr_Group group_list, Image_Data *image)
709 {
710 /* Define some constants */
711 #define PACK_BITS 12
712 #define PACK_BYTES 3
713 #define PACK_MASK 0x0F
714 #define PACK_SHIFT 4
715
716 /* Variables */
717 Acr_Element element;
718 int nrows, ncolumns;
719 int bits_alloc;
720 int bits_stored;
721 int image_group;
722 void *data = NULL;
723 long imagepix, ipix;
724 struct Acr_Element_Id elid;
725 nc_type datatype;
726 unsigned char *packed;
727 unsigned char pixel[2][2];
728 int need_byte_flip;
729 unsigned char temp_byte;
730 Acr_byte_order byte_order;
731
732 /* Get the image information */
733 bits_alloc = find_short(group_list, ACR_Bits_allocated, 0);
734 bits_stored = find_short(group_list, ACR_Bits_stored, 0);
735 nrows = find_short(group_list, ACR_Rows, 0);
736 ncolumns = find_short(group_list, ACR_Columns, 0);
737 if (nrows <= 0)
738 nrows = find_int(group_list, SPI_Recon_resolution, 0);
739 if (ncolumns <= 0)
740 ncolumns = find_int(group_list, SPI_Recon_resolution, 0);
741 image_group = find_short(group_list, ACR_Image_location, INT_MIN);
742
743 /* Figure out type */
744 if (bits_alloc > CHAR_BIT)
745 datatype = NC_SHORT;
746 else
747 datatype = NC_BYTE;
748
749 /* Set image info */
750 image->nrows = nrows;
751 image->ncolumns = ncolumns;
752 imagepix = nrows * ncolumns;
753 image->data = (unsigned short *) MALLOC(imagepix * sizeof(short));
754 image->free = TRUE;
755
756 /* Get image pointer */
757 elid.group_id = image_group;
758 elid.element_id = SPI_IMAGE_ELEMENT;
759 element = acr_find_group_element(group_list, &elid);
760 if (element == NULL) {
761 (void) memset(image->data, 0, imagepix * sizeof(short));
762 return;
763 }
764 data = acr_get_element_data(element);
765
766 /* Convert the data according to type */
767
768 /* Look for byte data */
769 if (datatype == NC_BYTE) {
770 for (ipix=0; ipix < imagepix; ipix++) {
771 image->data[ipix] = *((unsigned char *) data + ipix);
772 }
773 }
774 else {
775
776 /* Get byte order */
777 byte_order = acr_get_element_byte_order(element);
778
779 /* Look for unpacked short data */
780 if (bits_alloc == nctypelen(datatype) * CHAR_BIT) {
781 acr_get_short(byte_order, nrows*ncolumns, data, image->data);
782 }
783
784 /* Get packed short data */
785 else if ((bits_alloc == PACK_BITS) && (bits_stored <= bits_alloc) &&
786 ((imagepix % 2) == 0)) {
787 need_byte_flip = acr_need_invert(byte_order);
788 packed = (unsigned char *) data;
789 for (ipix=0; ipix < imagepix; ipix+=2) {
790 pixel[0][0] = packed[0];
791 pixel[0][1] = packed[1] & PACK_MASK;
792 pixel[1][0] = (packed[1] >> PACK_SHIFT) |
793 ((packed[2] & PACK_MASK) << PACK_SHIFT);
794 pixel[1][1] = packed[2] >> PACK_SHIFT;
795 if (need_byte_flip) {
796 temp_byte = pixel[0][0];
797 pixel[0][0] = pixel[0][1];
798 pixel[0][1] = temp_byte;
799 temp_byte = pixel[1][0];
800 pixel[1][0] = pixel[1][1];
801 pixel[1][1] = temp_byte;
802 }
803 image->data[ipix] = *((unsigned short *) pixel[0]);
804 image->data[ipix+1] = *((unsigned short *) pixel[1]);
805 packed += PACK_BYTES;
806 }
807 }
808
809 /* Fill with zeros in any other case */
810 else {
811 (void) memset(image->data, 0, imagepix * sizeof(short));
812 }
813 }
814
815 return;
816 }
817
818 /* ----------------------------- MNI Header -----------------------------------
819 @NAME : get_direction_cosines
820 @INPUT : angulation_ap - Angle of rotation about ap (Y) axis
821 angulation_lr - Angle of rotation about lr (X) axis
822 angulation_cc - Angle of rotation about cc (Z) axis
823 @OUTPUT : dircos - array of direction cosines
824 @RETURNS : (nothing)
825 @DESCRIPTION: Routine to compute direction cosines from angles
826 @METHOD : The rotation matrices are designed to be pre-multiplied by row
827 vectors. The rows of the final rotation matrix are the
828 direction cosines.
829 @GLOBALS :
830 @CALLS :
831 @CREATED : October 18, 1994 (Peter Neelin)
832 @MODIFIED :
833 ---------------------------------------------------------------------------- */
get_direction_cosines(double angulation_ap,double angulation_lr,double angulation_cc,double dircos[WORLD_NDIMS][WORLD_NDIMS])834 public void get_direction_cosines(double angulation_ap, double angulation_lr,
835 double angulation_cc,
836 double dircos[WORLD_NDIMS][WORLD_NDIMS])
837 {
838 World_Index iloop, jloop, kloop, axis_loop, axis;
839 double rotation[WORLD_NDIMS][WORLD_NDIMS];
840 double temp[WORLD_NDIMS][WORLD_NDIMS];
841 double angle, sinangle, cosangle;
842 /* Array giving order of rotations - cc, lr, ap */
843 static World_Index rotation_axis[WORLD_NDIMS]={ZCOORD, XCOORD, YCOORD};
844
845 /* Start with identity matrix */
846 for (iloop=0; iloop < WORLD_NDIMS; iloop++) {
847 for (jloop=0; jloop < WORLD_NDIMS; jloop++) {
848 if (iloop == jloop)
849 dircos[iloop][jloop] = 1.0;
850 else
851 dircos[iloop][jloop] = 0.0;
852 }
853 }
854
855 /* Loop through angles */
856 for (axis_loop=0; axis_loop < WORLD_NDIMS; axis_loop++) {
857
858 /* Get axis */
859 axis = rotation_axis[axis_loop];
860
861 /* Get angle */
862 switch (axis) {
863 case XCOORD: angle = angulation_lr; break;
864 case YCOORD: angle = angulation_ap; break;
865 case ZCOORD: angle = angulation_cc; break;
866 }
867 angle = angle / 180.0 * M_PI;
868 if (angle == 0.0) {
869 cosangle = 1.0;
870 sinangle = 0.0;
871 }
872 else {
873 cosangle = cos(angle);
874 sinangle = sin(angle);
875 }
876
877 /* Build up rotation matrix (make identity, then stuff in sin and cos) */
878 for (iloop=0; iloop < WORLD_NDIMS; iloop++) {
879 for (jloop=0; jloop < WORLD_NDIMS; jloop++) {
880 if (iloop == jloop)
881 rotation[iloop][jloop] = 1.0;
882 else
883 rotation[iloop][jloop] = 0.0;
884 }
885 }
886 rotation[(YCOORD+axis)%WORLD_NDIMS][(YCOORD+axis)%WORLD_NDIMS]
887 = cosangle;
888 rotation[(YCOORD+axis)%WORLD_NDIMS][(ZCOORD+axis)%WORLD_NDIMS]
889 = sinangle;
890 rotation[(ZCOORD+axis)%WORLD_NDIMS][(YCOORD+axis)%WORLD_NDIMS]
891 = -sinangle;
892 rotation[(ZCOORD+axis)%WORLD_NDIMS][(ZCOORD+axis)%WORLD_NDIMS]
893 = cosangle;
894
895 /* Multiply this by the previous matrix and then save the result */
896 for (iloop=0; iloop < WORLD_NDIMS; iloop++) {
897 for (jloop=0; jloop < WORLD_NDIMS; jloop++) {
898 temp[iloop][jloop] = 0.0;
899 for (kloop=0; kloop < WORLD_NDIMS; kloop++)
900 temp[iloop][jloop] +=
901 (dircos[iloop][kloop] * rotation[kloop][jloop]);
902 }
903 }
904 for (iloop=0; iloop < WORLD_NDIMS; iloop++)
905 for (jloop=0; jloop < WORLD_NDIMS; jloop++)
906 dircos[iloop][jloop] = temp[iloop][jloop];
907
908 }
909
910 }
911
912 /* ----------------------------- MNI Header -----------------------------------
913 @NAME : get_orientation_info
914 @INPUT : orientation - orientation of acquisition
915 dircos - direction cosines for each world axis
916 @OUTPUT : slice_world - world axis for slices
917 row_world - world axis for rows
918 column_world - world axis for columns
919 dircos - reordered direction cosines
920 @RETURNS : (nothing)
921 @DESCRIPTION: Routine to get the orientation information for a volume,
922 taking into account rotation of the axes. We need to
923 get the mapping from slice, row or column to world axis,
924 as well as re-ordering the direction cosines so that
925 they correspond to the appropriate axes.
926 @METHOD :
927 @GLOBALS :
928 @CALLS :
929 @CREATED : October 19, 1994 (Peter Neelin)
930 @MODIFIED :
931 ---------------------------------------------------------------------------- */
get_orientation_info(Orientation orientation,double dircos[WORLD_NDIMS][WORLD_NDIMS],World_Index * slice_world,World_Index * row_world,World_Index * column_world)932 public void get_orientation_info(Orientation orientation,
933 double dircos[WORLD_NDIMS][WORLD_NDIMS],
934 World_Index *slice_world,
935 World_Index *row_world,
936 World_Index *column_world)
937 {
938 World_Index original_to_rotated[WORLD_NDIMS];
939 World_Index original_axis;
940 static World_Index orientation_axes[NUM_ORIENTATIONS][VOL_NDIMS] = {
941 ZCOORD, YCOORD, XCOORD, /* TRANSVERSE */
942 XCOORD, ZCOORD, YCOORD, /* SAGITTAL */
943 YCOORD, ZCOORD, XCOORD /* CORONAL */
944 };
945
946 /* Get mapping from original world axis to rotated world axis */
947 for (original_axis = XCOORD; original_axis < WORLD_NDIMS; original_axis++) {
948 original_to_rotated[original_axis] =
949 get_nearest_world_axis(dircos[original_axis]);
950 }
951
952 /* Get mapping from volume axis to world axis (volume axes were defined
953 for original orientation, so we have to map those axes through the
954 rotation) */
955 *slice_world =
956 original_to_rotated[orientation_axes[orientation][VSLICE]];
957 *row_world =
958 original_to_rotated[orientation_axes[orientation][VROW]];
959 *column_world =
960 original_to_rotated[orientation_axes[orientation][VCOLUMN]];
961
962 /* Re-order direction cosines if needed */
963 if (XCOORD == original_to_rotated[YCOORD])
964 swap_dircos(dircos[XCOORD], dircos[YCOORD]);
965 else if (XCOORD == original_to_rotated[ZCOORD])
966 swap_dircos(dircos[XCOORD], dircos[ZCOORD]);
967 if (YCOORD == original_to_rotated[ZCOORD])
968 swap_dircos(dircos[YCOORD], dircos[ZCOORD]);
969
970 return;
971 }
972
973 /* ----------------------------- MNI Header -----------------------------------
974 @NAME : swap_dircos
975 @INPUT : dircos1 - first direction cosine
976 dircos2 - second direction cosine
977 @OUTPUT : dircos1 - first direction cosine
978 dircos2 - second direction cosine
979 @RETURNS : (nothing)
980 @DESCRIPTION: Routine to swap direction cosines in place
981 @METHOD :
982 @GLOBALS :
983 @CALLS :
984 @CREATED : October 19, 1994 (Peter Neelin)
985 @MODIFIED :
986 ---------------------------------------------------------------------------- */
swap_dircos(double dircos1[WORLD_NDIMS],double dircos2[WORLD_NDIMS])987 public void swap_dircos(double dircos1[WORLD_NDIMS],
988 double dircos2[WORLD_NDIMS])
989 {
990 World_Index iaxis;
991 double temp;
992
993 for (iaxis = XCOORD; iaxis < WORLD_NDIMS; iaxis++) {
994 temp = dircos1[iaxis];
995 dircos1[iaxis] = dircos2[iaxis];
996 dircos2[iaxis] = temp;
997 }
998
999 return;
1000 }
1001
1002 /* ----------------------------- MNI Header -----------------------------------
1003 @NAME : get_nearest_world_axis
1004 @INPUT : dircos - direction cosines
1005 @OUTPUT : (none)
1006 @RETURNS : Index label for nearest axis (XCOORD, YCOORD or ZCOORD).
1007 @DESCRIPTION: Routine to find the nearest axis to a direction cosine.
1008 @METHOD :
1009 @GLOBALS :
1010 @CALLS :
1011 @CREATED : October 19, 1994 (Peter Neelin)
1012 @MODIFIED :
1013 ---------------------------------------------------------------------------- */
get_nearest_world_axis(double dircos[WORLD_NDIMS])1014 public World_Index get_nearest_world_axis(double dircos[WORLD_NDIMS])
1015 {
1016 World_Index nearest_axis, iaxis;
1017
1018 nearest_axis = XCOORD;
1019 for (iaxis = YCOORD; iaxis < WORLD_NDIMS; iaxis++) {
1020 if (fabs(dircos[iaxis]) > fabs(dircos[nearest_axis]))
1021 nearest_axis = iaxis;
1022 }
1023
1024 return nearest_axis;
1025 }
1026
1027 /* ----------------------------- MNI Header -----------------------------------
1028 @NAME : calculate_start_from_centre
1029 @INPUT : slice_world - world index for slices
1030 row_world - world index for rows
1031 column_world - world index for columns
1032 nrows - number of rows in slice
1033 ncolumns - number of columns in slice
1034 centre - coordinate of centre of slice
1035 step - vector of steps (pixel separations)
1036 dircos - direction cosines for axes
1037 @OUTPUT : start - calculated start coordinate for slice
1038 @RETURNS : (nothing)
1039 @DESCRIPTION: Routine to calculate the start coordinate for a slice given
1040 its centre, step and direction cosines.
1041 @METHOD :
1042 @GLOBALS :
1043 @CALLS :
1044 @CREATED : November 18, 1994 (Peter Neelin)
1045 @MODIFIED :
1046 ---------------------------------------------------------------------------- */
1047 /* ARGSUSED */
calculate_slice_start(World_Index slice_world,World_Index row_world,World_Index column_world,int nrows,int ncolumns,double centre[WORLD_NDIMS],double step[WORLD_NDIMS],double dircos[WORLD_NDIMS][WORLD_NDIMS],double start[WORLD_NDIMS])1048 public void calculate_slice_start(World_Index slice_world,
1049 World_Index row_world,
1050 World_Index column_world,
1051 int nrows, int ncolumns,
1052 double centre[WORLD_NDIMS],
1053 double step[WORLD_NDIMS],
1054 double dircos[WORLD_NDIMS][WORLD_NDIMS],
1055 double start[WORLD_NDIMS])
1056 {
1057 World_Index iworld;
1058 double offset[WORLD_NDIMS];
1059
1060 /* Get offsets along each axis */
1061 offset[slice_world] = 0.0;
1062 offset[row_world] = (-step[row_world]) * (nrows - 1.0) / 2.0;
1063 offset[column_world] = (-step[column_world]) * (ncolumns - 1.0) / 2.0;
1064
1065 /* Transform offsets but don't use direction cosines since these files
1066 seem to give centres along the rotated axes, not along real axes. */
1067 for (iworld=0; iworld < WORLD_NDIMS; iworld++) {
1068 start[iworld] = centre[iworld] + offset[iworld];
1069 }
1070
1071 }
1072
find_short(Acr_Group group_list,Acr_Element_Id elid,int default_value)1073 public int find_short(Acr_Group group_list, Acr_Element_Id elid,
1074 int default_value)
1075 {
1076 return acr_find_short(group_list, elid, default_value);
1077 }
1078
find_int(Acr_Group group_list,Acr_Element_Id elid,int default_value)1079 public int find_int(Acr_Group group_list, Acr_Element_Id elid,
1080 int default_value)
1081 {
1082 return acr_find_int(group_list, elid, default_value);
1083 }
1084
find_double(Acr_Group group_list,Acr_Element_Id elid,double default_value)1085 public double find_double(Acr_Group group_list, Acr_Element_Id elid,
1086 double default_value)
1087 {
1088 return acr_find_double(group_list, elid, default_value);
1089 }
1090
find_string(Acr_Group group_list,Acr_Element_Id elid,char * default_value)1091 public char *find_string(Acr_Group group_list, Acr_Element_Id elid,
1092 char *default_value)
1093 {
1094 return acr_find_string(group_list, elid, default_value);
1095 }
1096
1097