Mercurial > hg > minc-tools
changeset 730:067a424554b0
Added code to automatically distinguish between old and new file types
(Numaris 2 and 3). Added code to do world coordinates for Numaris 2 (old).
| author | neelin <neelin> |
|---|---|
| date | Tue, 31 Jan 1995 13:36:08 +0000 |
| parents | b277e0f8384a |
| children | d253c4719c4c |
| files | conversion/mri_to_minc/siemens_to_minc.pl |
| diffstat | 1 files changed, 326 insertions(+), 51 deletions(-) [+] |
line wrap: on
line diff
--- a/conversion/mri_to_minc/siemens_to_minc.pl +++ b/conversion/mri_to_minc/siemens_to_minc.pl @@ -41,33 +41,214 @@ scalar(reverse(substr($string, $offset+2, 2)))) / 4.0; } -# Routine to initialize tape drive -sub siemens_initialize_tape_drive { +# Routines to do vector and matrix operations + +sub vector_normalize { + local(*vector) = @_; + local(@result); + $#result = 2; + local($scale); + local($mag) = 0; + for $i (0..2) { + $mag += ($vector[$i] ** 2); + } + $mag = sqrt($mag); + if ($mag > 0) {$scale = 1 / $mag;} + else {$scale = 1;} + for $i (0..2) { + $result[$i] = $scale * $vector[$i]; + } + return @result; +} + +sub vector_equal { + local(*vec1, *vec2) = @_; + local($result) =1; + for $i (0..2) { + $result = ($result && ($vec1[$i] == $vec2[$i])); + } + return $result; +} + +sub vector_dot_product { + local(*vec1, *vec2) = @_; + local(@result); + $#result = 2; + for $i (0..2) { + $result[$i] = $vec1[$i] * $vec2[$i]; + } + return @result; +} + +sub vector_cross_product { + local(*vec1, *vec2) = @_; + local(@result); + $#result = 2; + $result[0] = $vec1[1] * $vec2[2] - $vec1[2] * $vec2[1]; + $result[1] = $vec1[2] * $vec2[0] - $vec1[0] * $vec2[2]; + $result[2] = $vec1[0] * $vec2[1] - $vec1[1] * $vec2[0]; + return @result; +} + +sub matrix_multiply { + local(*mat1, *mat2) = @_; + local(@result); + $#result = 8; + local($temp); + for $i (0..2) { + for $j (0..2) { + $temp = 0; + for $k (0..2) { + $temp += $mat1[$i*3 + $k] * $mat2[$k*3 + $j]; + } + $result[$i*3+$j] = $temp; + } + } + return @result; +} + +sub matrix_vector_multiply { + local(*matrix, *vector) = @_; + local(@result); + $#result = 2; + local($temp); + for $i (0..2) { + $temp = 0; + for $j (0..2) { + $temp += $matrix[$i*3 + $j] * $vector[$j]; + } + $result[$i] = $temp; + } + return @result; +} + +sub matrix_transpose { + local(*matrix) = @_; + local(@result); + $#result = 8; + for $i (0..2) { + for $j (0..2) { + $result[$i*3 + $j] = $matrix[$j*3 + $i]; + } + } + return @result; } -# Subroutine to read the siemens file headers -sub siemens_read_headers { +sub matrix_create_transform { + local(*vec1, *vec2, *vec3) = @_; + local(@result); + $#result = 8; + for $i (0..2) { + $result[$i*3+0] = $vec1[$i]; + $result[$i*3+1] = $vec2[$i]; + $result[$i*3+2] = $vec3[$i]; + } + return @result; +} + +sub matrix_identity { + return (1,0,0, 0,1,0, 0,0,1); +} + +sub get_transform_from_normal { + local(*major, *normal) = @_; + local(@result); + local(@a, @b, @bprime); + local(@z, @zprime); + local(@first, @second); + local(@zero_vector) = (0,0,0); + @z = &vector_normalize(*major); + @zprime = &vector_normalize(*normal); + @a = &vector_cross_product(*z, *zprime); + @a = &vector_normalize(*a); + if (&vector_equal(*a, *zero_vector)) { + return &matrix_identity(); + } + @b = &vector_cross_product(*z, *a); + @bprime = &vector_cross_product(*zprime, *a); + @first = &matrix_create_transform(*a, *b, *z); + @first = &matrix_transpose(*first); + @second = &matrix_create_transform(*a, *bprime, *zprime); + @result = &matrix_multiply(*second, *first); + return @result; +} + +sub 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") { + @world_to_axis = (1,0,0, 0,1,0, 0,0,1); + } + elsif ($orientation eq "sagittal") { + @world_to_axis = (0,1,0, 0,0,1, 1,0,0); + } + elsif ($orientation eq "coronal") { + @world_to_axis = (1,0,0, 0,0,1, 0,1,0); + } + else {die "Unknown orientation \"$orientation\"";} + @axis_to_world = &matrix_transpose(*world_to_axis); + + # Transform normal and set closest axis + local(@new_normal) = &matrix_vector_multiply(*world_to_axis, *normal); + local(@major) = (0,0,1); + + # Get transform + local(@transform) = &get_transform_from_normal(*major, *new_normal); + + # Extract direction cosines + local(@transpose) = &matrix_transpose(*transform); + @col_dircos = @transpose[0..2]; + @row_dircos = @transpose[3..5]; + @slc_dircos = @transpose[6..8]; + + # Transform the direction cosines back to x,y,z (world) coordinates + @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); + +} + +sub convert_coordinates { + local(@coords) = @_; + $coords[0] *= -1; + $coords[2] *= -1; + return @coords; +} + +# NUMARIS 2 ROUTINES + +# Subroutine to read the siemens Numaris 2 file headers +sub numaris2_read_headers { # Set constants for reading file - local($blksize) = 512; - local($ident_hdr_off) = 0; - local($ident_hdr_len) = 244; + local($header_off) = 0; + local($header_len) = 8192; local($magic_off) = 244+42; local($magic_len) = 12; - local($patient_hdr_off) = 6*$blksize; - local($patient_hdr_len) = 512; - local($meas_hdr_off) = 7*$blksize; - local($meas_hdr_len) = 3072; - local($image_hdr_off) = 13*$blksize; - local($image_hdr_len) = 512; local($magic_value) = "SIEMENS MED "; + # Check arguements - if (scalar(@_) != 5) { - &cleanup_and_die("Argument error in siemens_read_headers",1); + if (scalar(@_) != 2) { + &cleanup_and_die("Argument error in numaris2_read_headers",1); } - local($filename, *ident_hdr, *patient_hdr, *meas_hdr, *image_hdr) = @_; - local($magic); + local($filename, *header) = @_; # Open the file if (!open(SMNF, "<".$filename)) { @@ -75,27 +256,17 @@ return 1; } - # Check the image file magic number - if (!seek(SMNF, $magic_off, 0) || - (read(SMNF, $magic, $magic_len) != $magic_len) || - ($magic ne $magic_value)) { - warn "Bad image file magic number in \"$filename\""; + # Read in the header + if (!seek(SMNF, $header_off, 0) || + (read(SMNF, $header, $header_len) != $header_len)) { return 1; } - # Read in file headers - return 1 - if (!seek(SMNF, $ident_hdr_off, 0) || - (read(SMNF, $ident_hdr, $ident_hdr_len) != $ident_hdr_len)); - return 1 - if (!seek(SMNF, $patient_hdr_off, 0) || - (read(SMNF, $patient_hdr, $patient_hdr_len) != $patient_hdr_len)); - return 1 - if (!seek(SMNF, $meas_hdr_off, 0) || - (read(SMNF, $meas_hdr, $meas_hdr_len) != $meas_hdr_len)); - return 1 - if (!seek(SMNF, $image_hdr_off, 0) || - (read(SMNF, $image_hdr, $image_hdr_len) != $image_hdr_len)); + # Check the image file magic number + if (substr($header, $magic_off, $magic_len) ne $magic_value) { + warn "Bad image file magic number in \"$filename\""; + return 1; + } # Close input file close(SMNF); @@ -103,38 +274,47 @@ return 0; } -# Routine to get file info -sub siemens_read_file_info { +# Routine to get Numaris 2 file info +sub numaris2_read_file_info { if (scalar(@_) != 3) { &cleanup_and_die("Argument error in read_file_info",1); } local($filename, *file_info, *specific_file_info) = @_; # Get headers - local($ident_hdr, $patient_hdr, $meas_hdr, $image_hdr); - if (&siemens_read_headers($filename, *ident_hdr, *patient_hdr, - *meas_hdr, *image_hdr)) { + local($header); + if (&numaris2_read_headers($filename, *header)) { return 1; } + # Get the file headers + local($blksize) = 512; + local($ident_hdr) = substr($header, 0, 244); + local($relat_hdr) = substr($header, 1*$blksize+162, 166); + local($patient_hdr) = substr($header, 6*$blksize, 512); + local($meas_hdr) = substr($header, 7*$blksize, 3072); + local($image_hdr) = substr($header, 13*$blksize, 512); + # Get interesting values - local($fovx, $fovy); $file_info{'numechos'} = &unpack_int(*meas_hdr, 512); if ($file_info{'numechos'} <= 0) {$file_info{'numechos'} = 1;} # We cheat and use patient id for exam, getting rid of weird characters ($file_info{'exam'} = &unpack_value(*patient_hdr, 28, 'A12')) =~ s/^\W//g; $file_info{'series'} = &unpack_int(*patient_hdr, 112); - $file_info{'image'} = &unpack_int(*meas_hdr, 364); + $file_info{'image'} = &unpack_value(*relat_hdr, 56, 'A6') + 0; + $file_info{'echo'} = &unpack_int(*meas_hdr, 516); $file_info{'width'} = &unpack_short(*image_hdr, 138); $file_info{'height'} = &unpack_short(*image_hdr, 140); $file_info{'pixel_size'} = 2; - $file_info{'slicepos'} = &unpack_float(*meas_hdr, 404); $file_info{'patient_name'} = &unpack_value(*patient_hdr, 0, "A26"); - local($norm_r) = &unpack_float(*meas_hdr, 380); - local($norm_a) = &unpack_float(*meas_hdr, 384); - local($norm_s) = &unpack_float(*meas_hdr, 388); + local(@normal) = &convert_coordinates(&unpack_float(*meas_hdr, 380), + &unpack_float(*meas_hdr, 384), + &unpack_float(*meas_hdr, 388)); + local($norm_r) = &abs($normal[0]); + local($norm_a) = &abs($normal[1]); + local($norm_s) = &abs($normal[2]); local($plane) = 'transverse'; local($max) = $norm_s; if ($norm_r > $max) { @@ -146,11 +326,47 @@ $max = $norm_a; } $file_info{'orientation'} = $plane; + + # Get coordinate information + local($fovx, $fovy); + local(@col_dircos, @row_dircos, @slc_dircos); $fovx = &unpack_float(*meas_hdr, 400); $fovy = &unpack_float(*meas_hdr, 400); + if ($fovy == 0) {$fovy = $fovx}; $file_info{'colstep'} = -$fovx / $file_info{'width'}; - $file_info{'rowstep'} = - -(($fovy != 0) ? $fovy : $fovx) / $file_info{'height'}; + $file_info{'rowstep'} = -$fovy / $file_info{'height'}; + &get_direction_cosines($file_info{'orientation'}, *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), + &unpack_float(*meas_hdr, 376)); + $file_info{'slicepos'} = $xcentre * $slc_dircos[0] + + $ycentre * $slc_dircos[1] + $zcentre * $slc_dircos[2]; + $file_info{'colstart'} = + ($xcentre * $col_dircos[0] + + $ycentre * $col_dircos[1] + + $zcentre * $col_dircos[2]) - + $file_info{'colstep'} * ($file_info{'width'} - 1) / 2; + $file_info{'rowstart'} = + ($xcentre * $row_dircos[0] + + $ycentre * $row_dircos[1] + + $zcentre * $row_dircos[2]) - + $file_info{'rowstep'} * ($file_info{'height'} - 1) / 2; + $file_info{'col_dircos_x'} = $col_dircos[0]; + $file_info{'col_dircos_y'} = $col_dircos[1]; + $file_info{'col_dircos_z'} = $col_dircos[2]; + $file_info{'row_dircos_x'} = $row_dircos[0]; + $file_info{'row_dircos_y'} = $row_dircos[1]; + $file_info{'row_dircos_z'} = $row_dircos[2]; + $file_info{'slc_dircos_x'} = $slc_dircos[0]; + $file_info{'slc_dircos_y'} = $slc_dircos[1]; + $file_info{'slc_dircos_z'} = $slc_dircos[2]; + + # Get other info $file_info{'tr'} = &unpack_float(*meas_hdr, 268)/1000; $file_info{'te'} = &unpack_float(*meas_hdr, 520)/1000; $file_info{'ti'} = &unpack_float(*meas_hdr, 264)/1000; @@ -158,7 +374,6 @@ $file_info{'scan_seq'} = &unpack_value(*meas_hdr, 172, 'a8'); $file_info{'scan_seq'} =~ s/\0.*$//; $file_info{'scan_seq'} =~ s/\n//; - ($file_info{'patient_birthdate'} = &unpack_value(*patient_hdr, 40, 'a10')) =~ s/\./-/g; local($sex_flag) = &unpack_value(*patient_hdr, 52, 'a2'); @@ -183,9 +398,9 @@ return 0; } -sub siemens_get_image_cmd { +sub numaris2_get_image_cmd { if (scalar(@_) != 2) { - &cleanup_and_die("Argument error in siemens_get_image_cmd",1); + &cleanup_and_die("Argument error in numaris2_get_image_cmd",1); } local($cur_file, *specific_file_info) = @_; @@ -198,6 +413,66 @@ return $cmd; } +# NUMARIS 3 ROUTINES + +# GENERAL SIEMENS ROUTINES + +# Routine to initialize tape drive +sub siemens_initialize_tape_drive { +} + +# Routine to read the file info +sub siemens_read_file_info { + if (scalar(@_) != 3) { + &cleanup_and_die("Argument error in read_file_info",1); + } + local($filename, *file_info, *specific_file_info) = @_; + + local($numaris2_magic) = "SIEMENS MED "; + local($numaris2_magic_id) = "0x0009 0x0011"; + local($numaris3_magic) = "SIEMENS CM VA0 CMS "; + local($numaris3_magic_id) = "0x0009 0x0012"; + + local($magic); + + # Check for a numaris 2 file + $magic = `extract_acr_nema $filename $numaris2_magic_id 2>/dev/null`; + if ($magic eq $numaris2_magic) { + $specific_file_info{'siemens_file_type'} = "Numaris2"; + return &numaris2_read_file_info(@_); + } + + # Check for a numaris 3 file + $magic = `extract_acr_nema $filename $numaris3_magic_id 2>/dev/null`; + if ($magic eq $numaris3_magic) { + $specific_file_info{'siemens_file_type'} = "Numaris3"; + return &numaris3_read_file_info(@_); + } + + # Can't figure out file type + warn "Unable to identify siemens file type in \"$filename\""; + return 1; + +} + +# Routine to get the image +sub siemens_get_image_cmd { + if (scalar(@_) != 2) { + &cleanup_and_die("Argument error in siemens_get_image_cmd",1); + } + local($cur_file, *specific_file_info) = @_; + + if ($specific_file_info{'siemens_file_type'} eq "Numaris2") { + return &numaris2_get_image_cmd(@_); + } + elsif ($specific_file_info{'siemens_file_type'} eq "Numaris3") { + return &numaris2_get_image_cmd(@_); + } + else { + die "Unknown siemens file type"; + } +} + # MAIN PROGRAM *initialize_tape_drive = *siemens_initialize_tape_drive;