1 /* ----------------------------- MNI Header -----------------------------------
2 @NAME : fix_dicom_coords.c
3 @DESCRIPTION: Routine to corretly set the dicom image coordinate information
4 based on the Philips MR shadow fields.
5 @GLOBALS : (none)
6 @CREATED : April 9, 2001 (Peter Neelin)
7 @MODIFIED :
8 ---------------------------------------------------------------------------- */
9
10 #include <stdio.h>
11 #include <stdlib.h>
12 #include <string.h>
13 #include <math.h>
14
15 #ifndef TRUE
16 # define TRUE 1
17 #endif
18 #ifndef FALSE
19 # define FALSE 0
20 #endif
21
22 /* World dimensions */
23 typedef enum { XCOORD = 0, YCOORD, ZCOORD, WORLD_NDIMS } World_Index;
24
25 /* Define spi constants */
26 #define SPI_TRANSVERSE_ORIENTATION 1
27 #define SPI_SAGITTAL_ORIENTATION 2
28 #define SPI_CORONAL_ORIENTATION 3
29
30 /* Define other constants */
31 #define private static
32
33 /* Public Function prototypes */
34 void calculate_dicom_coords(int orientation,
35 double angulation_lr, double angulation_ap,
36 double angulation_cc,
37 double off_centre_lr, double off_centre_ap,
38 double off_centre_cc,
39 double field_of_view,
40 double position[],
41 double row_dircos[], double col_dircos[],
42 double *location);
43 int convert_modified_imagenum(char *string);
44
45 /* Private function prototypes */
46 private void convert_coordinate(double coord[WORLD_NDIMS]);
47 private void get_direction_cosines(int orientation,
48 double angulation_ap, double angulation_lr,
49 double angulation_cc,
50 int use_safe_orientations,
51 double dircos[WORLD_NDIMS][WORLD_NDIMS]);
52 private void calculate_image_position(int orientation,
53 double row_fov,
54 double column_fov,
55 double centre[WORLD_NDIMS],
56 double dircos[WORLD_NDIMS][WORLD_NDIMS],
57 double position[WORLD_NDIMS]);
58 private void calculate_slice_location(int orientation,
59 double position[WORLD_NDIMS],
60 double dircos[WORLD_NDIMS][WORLD_NDIMS],
61 double *location);
62
63 /* ----------------------------- MNI Header -----------------------------------
64 @NAME : get_direction_cosines
65 @INPUT : orientation - shadow field (0x19, 0x100a)
66 angulation_lr - Angle of rotation about lr (X) axis
67 shadow field (0x19, 0x1007)
68 angulation_ap - Angle of rotation about ap (Y) axis
69 shadow field (0x19, 0x1006)
70 angulation_cc - Angle of rotation about cc (Z) axis
71 shadow field (0x19, 0x1005)
72 off_centre_lr - Centre of image along lr (X)
73 shadow field (0x19, 0x100b)
74 off_centre_ap - Centre of image along ap (Y)
75 shadow field (0x19, 0x100d)
76 off_centre_cc - Centre of image along cc (Z)
77 shadow field (0x19, 0x100c)
78 field_of_view - Width of field of view
79 shadow field (0x19, 0x1000)
80 @OUTPUT : position - dicom field (0x20, 0x32)
81 row_dircos, col_dircos - dicom field (0x20, 0x37)
82 location - retired dicom field (0x20, 0x1041)
83 @RETURNS : (nothing)
84 @DESCRIPTION: Routine to calculate dicom image position, orientation (row
85 and column direction cosines) and slice location from
86 Philips Gyroscan MR shadow fields.
87 @METHOD :
88 @GLOBALS :
89 @CALLS : get_direction_cosines
90 calculate_image_position
91 calculate_slice_location
92 @CREATED : April 9, 2001 (Peter Neelin)
93 @MODIFIED :
94 ---------------------------------------------------------------------------- */
calculate_dicom_coords(int orientation,double angulation_lr,double angulation_ap,double angulation_cc,double off_centre_lr,double off_centre_ap,double off_centre_cc,double field_of_view,double position[],double row_dircos[],double col_dircos[],double * location)95 void calculate_dicom_coords(int orientation,
96 double angulation_lr, double angulation_ap,
97 double angulation_cc,
98 double off_centre_lr, double off_centre_ap,
99 double off_centre_cc,
100 double field_of_view,
101 double position[],
102 double row_dircos[], double col_dircos[],
103 double *location)
104 {
105 double dircos[WORLD_NDIMS][WORLD_NDIMS];
106 World_Index row_world, column_world, iworld;
107 double centre[3];
108
109 /*
110 * Work out the axis direction cosines
111 */
112 get_direction_cosines(orientation,
113 angulation_ap, angulation_lr, angulation_cc,
114 FALSE, dircos);
115
116 /*
117 * Map the world axes to the image axes
118 */
119 switch (orientation) {
120 case SPI_SAGITTAL_ORIENTATION:
121 row_world = YCOORD;
122 column_world = ZCOORD;
123 break;
124 case SPI_CORONAL_ORIENTATION:
125 row_world = XCOORD;
126 column_world = ZCOORD;
127 break;
128 case SPI_TRANSVERSE_ORIENTATION:
129 default:
130 row_world = XCOORD;
131 column_world = YCOORD;
132 break;
133 }
134 for (iworld=0; iworld < WORLD_NDIMS; iworld++) {
135 row_dircos[iworld] = dircos[row_world][iworld];
136 col_dircos[iworld] = dircos[column_world][iworld];
137 }
138
139 /*
140 * Calculate the dicom image position
141 */
142 centre[XCOORD] = off_centre_lr;
143 centre[YCOORD] = off_centre_ap;
144 centre[ZCOORD] = off_centre_cc;
145 calculate_image_position(orientation, field_of_view, field_of_view,
146 centre, dircos, position);
147
148 /*
149 * Calculate the slice location.
150 */
151 calculate_slice_location(orientation, position, dircos, location);
152 }
153
154 /* ----------------------------- MNI Header -----------------------------------
155 @NAME : convert_coordinate
156 @INPUT : coord - coordinate according to old convention
157 @OUTPUT : coord - coordinate converted to DICOM convention
158 @RETURNS : (nothing)
159 @DESCRIPTION: Routine to convert coordinates to and from DICOM.
160 @METHOD : Since this operation is its own inverse, we do not provide a
161 flag to indicate which direction we are going. If this
162 ever changes, then the interface will have to be changed.
163 (It would have been a good idea to begin with, but we all
164 know how code evolves...)
165 @GLOBALS :
166 @CALLS :
167 @CREATED : October 18, 1994 (Peter Neelin)
168 @MODIFIED :
169 ---------------------------------------------------------------------------- */
convert_coordinate(double coord[WORLD_NDIMS])170 private void convert_coordinate(double coord[WORLD_NDIMS])
171 {
172 coord[XCOORD] *= -1.0;
173 coord[YCOORD] *= -1.0;
174 }
175
176 /* ----------------------------- MNI Header -----------------------------------
177 @NAME : get_direction_cosines
178 @INPUT : orientation
179 angulation_ap - Angle of rotation about ap (Y) axis
180 angulation_lr - Angle of rotation about lr (X) axis
181 angulation_cc - Angle of rotation about cc (Z) axis
182 use_safe_orientations - TRUE if dir cos should be realigned
183 to safe values for machines that do not handle angled axes
184 @OUTPUT : dircos - array of direction cosines
185 @RETURNS : (nothing)
186 @DESCRIPTION: Routine to compute direction cosines from angles
187 @METHOD : The rotation matrices are designed to be pre-multiplied by row
188 vectors. The rows of the final rotation matrix are the
189 direction cosines.
190 @GLOBALS :
191 @CALLS : sin, cos, fabs, convert_coordinate
192 @CREATED : October 18, 1994 (Peter Neelin)
193 @MODIFIED :
194 ---------------------------------------------------------------------------- */
get_direction_cosines(int orientation,double angulation_ap,double angulation_lr,double angulation_cc,int use_safe_orientations,double dircos[WORLD_NDIMS][WORLD_NDIMS])195 private void get_direction_cosines(int orientation,
196 double angulation_ap, double angulation_lr,
197 double angulation_cc,
198 int use_safe_orientations,
199 double dircos[WORLD_NDIMS][WORLD_NDIMS])
200 {
201 World_Index iloop, jloop, kloop, axis_loop, axis;
202 double rotation[WORLD_NDIMS][WORLD_NDIMS];
203 double temp[WORLD_NDIMS][WORLD_NDIMS];
204 double angle, sinangle, cosangle;
205 /* Array giving order of rotations - cc, lr, ap */
206 static World_Index rotation_axis[WORLD_NDIMS]={ZCOORD, XCOORD, YCOORD};
207
208 /* Start with identity matrix */
209 for (iloop=0; iloop < WORLD_NDIMS; iloop++) {
210 for (jloop=0; jloop < WORLD_NDIMS; jloop++) {
211 if (iloop == jloop)
212 dircos[iloop][jloop] = 1.0;
213 else
214 dircos[iloop][jloop] = 0.0;
215 }
216 }
217
218 /* Loop through angles */
219 for (axis_loop=0; axis_loop < WORLD_NDIMS; axis_loop++) {
220
221 /* Get axis */
222 axis = rotation_axis[axis_loop];
223
224 /* Get angle */
225 switch (axis) {
226 case XCOORD: angle = angulation_lr; break;
227 case YCOORD: angle = angulation_ap; break;
228 case ZCOORD: angle = angulation_cc; break;
229 }
230 angle = angle / 180.0 * M_PI;
231 if (angle == 0.0) {
232 cosangle = 1.0;
233 sinangle = 0.0;
234 }
235 else {
236 cosangle = cos(angle);
237 sinangle = sin(angle);
238 }
239
240 /* Build up rotation matrix (make identity, then stuff in sin and cos) */
241 for (iloop=0; iloop < WORLD_NDIMS; iloop++) {
242 for (jloop=0; jloop < WORLD_NDIMS; jloop++) {
243 if (iloop == jloop)
244 rotation[iloop][jloop] = 1.0;
245 else
246 rotation[iloop][jloop] = 0.0;
247 }
248 }
249 rotation[(YCOORD+axis)%WORLD_NDIMS][(YCOORD+axis)%WORLD_NDIMS]
250 = cosangle;
251 rotation[(YCOORD+axis)%WORLD_NDIMS][(ZCOORD+axis)%WORLD_NDIMS]
252 = sinangle;
253 rotation[(ZCOORD+axis)%WORLD_NDIMS][(YCOORD+axis)%WORLD_NDIMS]
254 = -sinangle;
255 rotation[(ZCOORD+axis)%WORLD_NDIMS][(ZCOORD+axis)%WORLD_NDIMS]
256 = cosangle;
257
258 /* Multiply this by the previous matrix and then save the result */
259 for (iloop=0; iloop < WORLD_NDIMS; iloop++) {
260 for (jloop=0; jloop < WORLD_NDIMS; jloop++) {
261 temp[iloop][jloop] = 0.0;
262 for (kloop=0; kloop < WORLD_NDIMS; kloop++)
263 temp[iloop][jloop] +=
264 (dircos[iloop][kloop] * rotation[kloop][jloop]);
265 }
266 }
267 for (iloop=0; iloop < WORLD_NDIMS; iloop++)
268 for (jloop=0; jloop < WORLD_NDIMS; jloop++)
269 dircos[iloop][jloop] = temp[iloop][jloop];
270
271 }
272
273 /* Convert the coordinates to the DICOM standard */
274 for (iloop=0; iloop < WORLD_NDIMS; iloop++) {
275 convert_coordinate(dircos[iloop]);
276 }
277
278 /* Invert direction cosines to account for flipped dimensions */
279 for (iloop=0; iloop < WORLD_NDIMS; iloop++) {
280 if ((iloop == XCOORD) && (orientation == SPI_SAGITTAL_ORIENTATION))
281 continue;
282 if ((iloop == ZCOORD) && (orientation == SPI_TRANSVERSE_ORIENTATION))
283 continue;
284 for (jloop=0; jloop < WORLD_NDIMS; jloop++) {
285 dircos[iloop][jloop] *= -1.0;
286 }
287 }
288
289 /* Kludge to handle the fact the the Picker software does not seem
290 to handle rotated volumes properly */
291 if (use_safe_orientations) {
292 double biggest;
293 int dimused[WORLD_NDIMS];
294
295 /* Force all direction cosines to be along major axes */
296
297 /* Keep track of axes that have been used (in case of 45 degree
298 rotations) */
299 for (iloop=0; iloop < WORLD_NDIMS; iloop++)
300 dimused[iloop] = FALSE;
301
302 /* Loop over all direction cosines */
303 for (iloop=0; iloop < WORLD_NDIMS; iloop++) {
304
305 /* Start with the first unused dimension */
306 for (kloop=0; dimused[kloop]; kloop++) {}
307 biggest = dircos[iloop][kloop];
308
309 /* Loop over all dimensions, looking for the biggest unused one */
310 for (jloop=0; jloop < WORLD_NDIMS; jloop++) {
311 if ((fabs(biggest) < fabs(dircos[iloop][jloop])) &&
312 !dimused[jloop]) {
313 kloop = jloop;
314 biggest = dircos[iloop][kloop];
315 }
316 dircos[iloop][jloop] = 0.0;
317 }
318
319 /* Set the biggest coordinate to +/- 1 */
320 if (biggest > 0.0)
321 dircos[iloop][kloop] = 1.0;
322 else
323 dircos[iloop][kloop] = -1.0;
324 dimused[kloop] = TRUE;
325
326 } /* End of loop over dimension cosines */
327
328 } /* End if use_safe_orientations */
329
330 }
331
332 /* ----------------------------- MNI Header -----------------------------------
333 @NAME : calculate_image_position
334 @INPUT : orientation - indicates orientation of slice
335 row_fov - field-of-view for rows
336 column_fov - field-of-view for columns
337 centre - coordinate of centre of slice as a rotated Philips
338 coordinate
339 dircos - direction cosines for axes
340 @OUTPUT : position - calculated position coordinate for slice
341 @RETURNS : (nothing)
342 @DESCRIPTION: Routine to calculate the position (top left) coordinate for
343 a slice given its centre, field-of-view and direction cosines.
344 @METHOD : Direction cosines are converted back to Philips coordinate
345 convention before being used to rotate the centre. The final
346 position is converted to DICOM convention at the end.
347 @GLOBALS :
348 @CALLS : fabs, convert_coordinate
349 @CREATED : September 2, 1997 (Peter Neelin)
350 @MODIFIED :
351 ---------------------------------------------------------------------------- */
calculate_image_position(int orientation,double row_fov,double column_fov,double centre[WORLD_NDIMS],double dircos[WORLD_NDIMS][WORLD_NDIMS],double position[WORLD_NDIMS])352 private void calculate_image_position(int orientation,
353 double row_fov,
354 double column_fov,
355 double centre[WORLD_NDIMS],
356 double dircos[WORLD_NDIMS][WORLD_NDIMS],
357 double position[WORLD_NDIMS])
358 {
359 double coord[WORLD_NDIMS], pos_dircos[WORLD_NDIMS][WORLD_NDIMS];
360 double row_offset, column_offset;
361 double biggest, flip;
362 int slice_world, row_world, column_world;
363 World_Index iloop, jloop;
364
365 /* Figure out indices */
366 switch (orientation) {
367 case SPI_SAGITTAL_ORIENTATION:
368 slice_world = XCOORD;
369 row_world = ZCOORD;
370 column_world = YCOORD;
371 break;
372 case SPI_CORONAL_ORIENTATION:
373 slice_world = YCOORD;
374 row_world = ZCOORD;
375 column_world = XCOORD;
376 break;
377 case SPI_TRANSVERSE_ORIENTATION:
378 default:
379 slice_world = ZCOORD;
380 row_world = YCOORD;
381 column_world = XCOORD;
382 break;
383 }
384
385 /* Work out positive direction cosines and direction of row and column
386 offsets. */
387 for (iloop=0; iloop < WORLD_NDIMS; iloop++) {
388
389 /* Switch the direction cosine back to Philips orientation */
390 for (jloop=0; jloop < WORLD_NDIMS; jloop++) {
391 pos_dircos[iloop][jloop] = dircos[iloop][jloop];
392 }
393 convert_coordinate(pos_dircos[iloop]);
394
395 /* Find biggest component and figure out if we need to flip the
396 direction cosine */
397 biggest = pos_dircos[iloop][0];
398 for (jloop=1; jloop < WORLD_NDIMS; jloop++) {
399 if (fabs(biggest) < fabs(pos_dircos[iloop][jloop])) {
400 biggest = pos_dircos[iloop][jloop];
401 }
402 }
403 flip = ((biggest >= 0.0) ? +1.0 : -1.0);
404
405 /* Make the direction cosine positive */
406 for (jloop=0; jloop < WORLD_NDIMS; jloop++) {
407 pos_dircos[iloop][jloop] *= flip;
408 }
409
410 /* Work out row and column offsets */
411 if (iloop == row_world) {
412 row_offset = flip * row_fov / 2.0;
413 }
414 if (iloop == column_world) {
415 column_offset = flip * column_fov / 2.0;
416 }
417
418 }
419
420 /* Calculate position - note that centres are given in rotated frame */
421 coord[slice_world] = centre[slice_world];
422 coord[row_world] = centre[row_world] - row_offset;
423 coord[column_world] = centre[column_world] - column_offset;
424
425 /* Rotate the coordinate according to the direction cosines */
426 for (iloop=0; iloop < WORLD_NDIMS; iloop++) {
427 position[iloop] = 0.0;
428 for (jloop=0; jloop < WORLD_NDIMS; jloop++) {
429 position[iloop] += pos_dircos[jloop][iloop] * coord[jloop];
430 }
431 }
432
433 /* Convert the coordinate back to DICOM */
434 convert_coordinate(position);
435
436 }
437
438 /* ----------------------------- MNI Header -----------------------------------
439 @NAME : calculate_slice_location
440 @INPUT : orientation - indicates orientation of slice
441 position - position of slice
442 dircos - direction cosines for axes
443 @OUTPUT : location - offset perpendicular to slice
444 @RETURNS : (nothing)
445 @DESCRIPTION: Routine to calculate the location of the slice along an
446 axis perpendicular to it.
447 @METHOD :
448 @GLOBALS :
449 @CALLS :
450 @CREATED : July 11, 2000 (Peter Neelin)
451 @MODIFIED :
452 ---------------------------------------------------------------------------- */
calculate_slice_location(int orientation,double position[WORLD_NDIMS],double dircos[WORLD_NDIMS][WORLD_NDIMS],double * location)453 private void calculate_slice_location(int orientation,
454 double position[WORLD_NDIMS],
455 double dircos[WORLD_NDIMS][WORLD_NDIMS],
456 double *location)
457 {
458 double distance;
459 int index, slice_world;
460
461 /* Figure out which direction cosine to use */
462 switch (orientation) {
463 case SPI_SAGITTAL_ORIENTATION:
464 slice_world = XCOORD;
465 break;
466 case SPI_CORONAL_ORIENTATION:
467 slice_world = YCOORD;
468 break;
469 case SPI_TRANSVERSE_ORIENTATION:
470 default:
471 slice_world = ZCOORD;
472 break;
473 }
474
475 /* Project the position onto the appropriate direction cosine */
476 distance = 0.0;
477 for (index=0; index < WORLD_NDIMS; index++) {
478 distance += dircos[slice_world][index] * position[index];
479 }
480
481 /* Check for a negative direction cosine */
482 if (dircos[slice_world][slice_world] < 0.0) {
483 distance *= -1.0;
484 }
485
486 /* Save the result */
487 *location = distance;
488
489 }
490
491 /* ----------------------------- MNI Header -----------------------------------
492 @NAME : convert_modified_imagenum
493 @INPUT : string - string to be converted in place
494 @OUTPUT : string - modified date string
495 @RETURNS : TRUE if string was modified, FALSE otherwise.
496 @DESCRIPTION: Routine to convert an image number from gyroscan format to
497 a valid dicom integer by making it smaller, if necessary.
498 This function differs from the original in that it assumes that
499 the string has already had two characters chopped off the end.
500 @METHOD :
501 @GLOBALS :
502 @CALLS :
503 @CREATED : March 5, 2001 (Peter Neelin)
504 @MODIFIED :
505 ---------------------------------------------------------------------------- */
convert_modified_imagenum(char * string)506 int convert_modified_imagenum(char *string)
507 {
508 int istart, iend, iold, inew, oldlen, ipos;
509 static Acr_max_IS_len = 8;
510 /* static int chars_to_remove[] = {1,2,5,6,13}; */
511 static int chars_to_remove[] = {3,4,11};
512 static int nchars_to_remove =
513 sizeof(chars_to_remove) / sizeof(chars_to_remove[0]);
514 int ichar;
515
516 /* Check for NULL string */
517 if (string == NULL) return FALSE;
518
519 /* Get the original length of the string */
520 oldlen = strlen(string);
521
522 /* Check for zero-length string */
523 if (oldlen == 0) return FALSE;
524
525 /* Find the start of the number */
526 for (istart=0; (string[istart] != '\0') &&
527 (string[istart] < '0' || string[istart] > '9'); istart++) {}
528
529 /* Find the end of the string */
530 for (iend=istart; (string[iend] >= '0') && (string[iend] <= '9'); iend++) {}
531
532 /* Does anything need to be fixed? If not, then return. */
533 if (iend-istart <= Acr_max_IS_len) {
534 return FALSE;
535 }
536
537 /* Copy the string */
538 inew = 0;
539 iold = istart;
540 ichar = nchars_to_remove - 1;
541 ipos = iend-iold;
542 while ((chars_to_remove[ichar] > ipos) && (ichar > 0)) {
543 ichar--;
544 }
545 while (string[iold] >= '0' && string[iold] <= '9') {
546
547 /* Get position from end of string */
548 ipos = iend-iold;
549
550 /* Skip any chars to the left of the highest position, and skip
551 the specified character positions and skip any leading zeros */
552 if ((ipos > chars_to_remove[ichar]) && (ichar == nchars_to_remove-1)) {}
553 else if (ipos == chars_to_remove[ichar]) {
554 if (ichar > 0) ichar--;
555 }
556 else if ((inew == 0) && (string[iold] == '0')) {}
557 else {
558
559 if (inew != iold)
560 string[inew] = string[iold];
561 inew++;
562
563 }
564
565 iold++;
566 }
567
568 /* Check for an empty string */
569 if (inew == 0) {
570 string[inew] = '0';
571 inew++;
572 }
573
574 /* Pad the string with NULs */
575 do {
576 string[inew] = '\0';
577 inew++;
578 } while (inew < oldlen);
579
580 return TRUE;
581 }
582
583