Mercurial > hg > minc-tools
changeset 733:61750edba670
Changed calculation of direction cosines (use angles instead).
author | neelin <neelin> |
---|---|
date | Thu, 02 Feb 1995 21:14:07 +0000 |
parents | f434f9a1262f |
children | f166724c3e1c |
files | conversion/mri_to_minc/siemens_to_minc.pl |
diffstat | 1 files changed, 163 insertions(+), 19 deletions(-) [+] |
line wrap: on
line diff
--- a/conversion/mri_to_minc/siemens_to_minc.pl +++ b/conversion/mri_to_minc/siemens_to_minc.pl @@ -150,6 +150,43 @@ return (1,0,0, 0,1,0, 0,0,1); } +sub matrix_xrot { + local($angle) = @_; + local(@result); + local($sin, $cos); + + $sin = sin($angle); + $cos = cos($angle); + return (1,0,0, 0,$cos,$sin, 0,-$sin,$cos); +} + +sub matrix_yrot { + local($angle) = @_; + local(@result); + local($sin, $cos); + + $sin = sin($angle); + $cos = cos($angle); + return ($cos,0,-$sin, 0,1,0, $sin,0,$cos); +} + +sub matrix_zrot { + local($angle) = @_; + local(@result); + local($sin, $cos); + + $sin = sin($angle); + $cos = cos($angle); + return ($cos,$sin,0, -$sin,$cos,0, 0,0,1); +} + +sub matrix_extract_row { + local(*matrix, $row) = @_; + + local($offset) = $row*3; + return @matrix[$offset..$offset+2]; +} + sub get_transform_from_normal { local(*major, *normal) = @_; local(@result); @@ -173,23 +210,9 @@ return @result; } -sub get_direction_cosines { +sub old_get_direction_cosines { local($orientation, *normal, *col_dircos, *row_dircos, *slc_dircos) = @_; - # Check for normal that is not in one of the major planes - local($danger_normal) = 1; - foreach $coord (@normal) { - if ($coord == 0.0) { - $danger_normal = 0; - } - } - if ($danger_normal) { - warn "\n\nWARNING!!!!! ". - "Found old Siemens file with dangerous image normal (2 rotations).\n". - "Getting out of here fast! Time to fix this program!\n\n\n"; - &cleanup_and_die("Sorry!\n", -1); - } - # Get transformation from col,row,slice (axis) to x,y,z (world) local(@axis_to_world, @world_to_axis); if ($orientation eq "transverse") { @@ -221,6 +244,126 @@ @col_dircos = &matrix_vector_multiply(*axis_to_world, *col_dircos); @row_dircos = &matrix_vector_multiply(*axis_to_world, *row_dircos); @slc_dircos = &matrix_vector_multiply(*axis_to_world, *slc_dircos); + @col_dircos = &get_dircos(@col_dircos); + @row_dircos = &get_dircos(@row_dircos); + @slc_dircos = &get_dircos(@slc_dircos); + +} + +sub get_transform_from_axes { + local($first, $second, $angle) = @_; + + $angle *= 3.14159265358979323846/180.0; + + if (($first eq 'tra') && ($second eq 'sag')) { + return &matrix_yrot($angle); + } + elsif (($first eq 'sag') && ($second eq 'tra')) { + return &matrix_yrot(-$angle); + } + elsif (($first eq 'sag') && ($second eq 'cor')) { + return &matrix_zrot($angle); + } + elsif (($first eq 'cor') && ($second eq 'sag')) { + return &matrix_zrot(-$angle); + } + elsif (($first eq 'cor') && ($second eq 'tra')) { + return &matrix_xrot($angle); + } + elsif (($first eq 'tra') && ($second eq 'cor')) { + return &matrix_xrot(-$angle); + } + else { + return &matrix_identity(); + } +} + +sub get_transform_from_orient_strings { + local($orient_str1, $orient_str2) = @_; + local(@result); + + # Parse the orientation strings + local($axis_list) = 'tra|sag|cor'; + $orient_str1 =~ tr/A-Z/a-z/; + $orient_str2 =~ tr/A-Z/a-z/; + if ($orient_str1 =~ /^\s*($axis_list)>($axis_list)\s*(\S*)/) { + $major_axis = $1; + $first_axis = $2; + $first_angle = $3; + if ($orient_str2 =~ /^\s*>($axis_list)\s*(\S*)/) { + $second_axis = $1; + $second_angle = $2; + } + else { + $second_axis = ""; + $second_angle = 0; + } + } + else { + $major_axis = ""; + $first_axis = ""; + $first_angle = 0; + $second_axis = ""; + $second_angle = 0; + } + + # Create rotation matrix + local(@mat1, @mat2); + @mat1 = &get_transform_from_axes($major_axis, $first_axis, $first_angle); + @mat2 = &get_transform_from_axes($major_axis, $second_axis, $second_angle); + @result = &matrix_multiply(*mat1, *mat2); + + return @result; +} + +sub get_direction_cosines { + local($orientation, $orient_str1, $orient_str2, *normal, + *col_dircos, *row_dircos, *slc_dircos) = @_; + + # Get transformation from col,row,slice (axis) to x,y,z (world) + local(@axis_row); + if ($orientation eq "transverse") { + @axis_row = (0,1,2); + } + elsif ($orientation eq "sagittal") { + @axis_row = (1,2,0); + } + elsif ($orientation eq "coronal") { + @axis_row = (0,2,1); + } + else {die "Unknown orientation \"$orientation\"";} + + # Extract rotation angles from strings + local(@transform) = + &get_transform_from_orient_strings($orient_str1, $orient_str2); + + # Extract direction cosines + local(@transpose) = &matrix_transpose(*transform); + @col_dircos = &matrix_extract_row(*transpose, $axis_row[0]); + @row_dircos = &matrix_extract_row(*transpose, $axis_row[1]); + @slc_dircos = &matrix_extract_row(*transpose, $axis_row[2]); + + # Fix up the coordinates + @col_dircos = &convert_coordinates(@col_dircos); + @row_dircos = &convert_coordinates(@row_dircos); + @slc_dircos = &convert_coordinates(@slc_dircos); + @col_dircos = &get_dircos(@col_dircos); + @row_dircos = &get_dircos(@row_dircos); + @slc_dircos = &get_dircos(@slc_dircos); + + # Compare to normal + @normal = &get_dircos(@normal); + local($diff) = 0; + foreach $i (0..2) { + $diff += ($normal[$i] - $slc_dircos[$i]) ** 2; + } + $diff = sqrt($diff); + if ($diff > 1.0e-6) { + warn "\n\nWARNING!!!!! ". + "Calculated image normal doesn't match file value.\n". + "Something went wrong in dircos calculation for Numaris 2!\n\n\n"; + &cleanup_and_die("Sorry!\n", -1); + } } @@ -294,6 +437,7 @@ local($patient_hdr) = substr($header, 6*$blksize, 512); local($meas_hdr) = substr($header, 7*$blksize, 3072); local($image_hdr) = substr($header, 13*$blksize, 512); + local($text_hdr) = substr($header, 14*$blksize, 512); # Get interesting values $file_info{'numechos'} = &unpack_int(*meas_hdr, 512); @@ -335,11 +479,11 @@ if ($fovy == 0) {$fovy = $fovx}; $file_info{'colstep'} = -$fovx / $file_info{'width'}; $file_info{'rowstep'} = -$fovy / $file_info{'height'}; - &get_direction_cosines($file_info{'orientation'}, *normal, + local($orient_str1) = &unpack_value(*text_hdr, 332, "A11"); + local($orient_str2) = &unpack_value(*text_hdr, 343, "A11"); + &get_direction_cosines($file_info{'orientation'}, + $orient_str1, $orient_str2, *normal, *col_dircos, *row_dircos, *slc_dircos); - @col_dircos = &get_dircos(@col_dircos); - @row_dircos = &get_dircos(@row_dircos); - @slc_dircos = &get_dircos(@slc_dircos); local($xcentre, $ycentre, $zcentre) = &convert_coordinates(&unpack_float(*meas_hdr, 368), &unpack_float(*meas_hdr, 372),