1#! /usr/local/bin/perl
2#
3# Routines for converting mri files to minc format. Must be included
4# with routines for reading in data structures, specific to the input
5# file format.
6#
7
8# Global variable for indicating byte order
9$Little_endian = unpack("c",substr(pack("s",1),0,1));
10%Unpack_codes = ('a', 1, 'A', 1,
11                 's', 2, 'S', 2,
12                 'i', 4, 'I', 4,
13                 'f', 4, 'd', 8,
14                 );
15
16
17sub numeric_order { $a <=> $b;}
18
19# Routine to take absolute value
20sub abs {
21    local(@new, $val);
22    foreach $val (@_) {
23        push(@new, ($val<=>0) * $val);
24    }
25    return (scalar(@new) > 1) ? @new : $new[0];
26}
27
28# Subroutine to clean up files and exit
29sub cleanup_and_die {
30
31    # Get message to print and exit status
32    local($message,$status) = @_;
33    if (!defined($status)) {$status = 0;}
34    if (defined($message)) {
35        print STDERR $message;
36        if ($message !~ /\n$/) {print STDERR "\n";}
37    }
38
39    $SIG{'INT'}  = 'IGNORE';
40    $SIG{'TERM'} = 'IGNORE';
41    $SIG{'QUIT'} = 'IGNORE';
42    # Check for temp files
43    if (defined($tmpdir) && -e $tmpdir) {
44        print STDERR "Cleaning up temporary files.\n";
45        system "rm -rf $tmpdir";
46    }
47
48    exit($status);
49}
50
51# Subroutine to execute commands and check return status
52sub execute {
53    local($status);
54    if ($print_execution_statements) {
55        print join(' ',@_),"\n";
56    }
57    $status = system @_;
58    if ($status != 0) {
59        &cleanup_and_die("Error executing command \"".join(" ",@_)."\".\n",
60                         $status);
61    }
62}
63
64# Subroutine to remove a list of files
65sub remove_file {
66    unlink @_;
67}
68
69# Subroutine to create a temporary directory (global variable $tmpdir
70# is created)
71sub create_tmpdir {
72    local($subdir);
73    if (scalar(@_) >= 1) {
74        $subdir = $_[0];
75    }
76    else {
77        $subdir = $$;
78    }
79    if (defined($ENV{TMPDIR})) {
80        $tmpdir = $ENV{TMPDIR};
81    }
82    else {
83        $tmpdir = "/usr/tmp";
84    }
85    $tmpdir .= "/$subdir";
86    mkdir($tmpdir,0777) || die "Unable to create directory $tmpdir: $!";
87    $SIG{'INT'}  = 'cleanup_and_die';
88    $SIG{'TERM'} = 'cleanup_and_die';
89    $SIG{'QUIT'} = 'cleanup_and_die';
90}
91
92# Routine to execute external tape commands, returning the error
93sub catch_tape_command {
94    local($tapedrive) = shift(@_);
95    close(TAPEHANDLE);
96    local($status) = system(@_);
97    open(TAPEHANDLE, $tapedrive);
98    return $status;
99}
100
101# Routine to execute external tape commands, dying on error
102sub execute_tape_command {
103    local($status) = &catch_tape_command(@_);
104    if ($status != 0) {
105        shift(@_);
106        &cleanup_and_die("Error executing command \"".join(" ",@_)."\".\n",
107                         $status);
108    }
109}
110
111# Routine to get the current tape position
112sub get_tape_position {
113    if (defined($counter_for_read_next_file)) {
114        return $counter_for_read_next_file;
115    }
116    else {
117        return 0;
118    }
119}
120
121# Routine to set the current tape position
122sub set_tape_position {
123    local($tapedrive, $position) = @_;
124    if (length($tapedrive) <= 0) {return;}
125
126    # Loop to position tape
127    local($repos_sleep) = 1;
128    local($nreposition_retries) = 4;
129    local($tmp_status);
130    foreach $reposloop (0..$nreposition_retries-1) {
131        select(undef, undef, undef, $repos_sleep);
132        $tmp_status = &catch_tape_command($tapedrive,
133                                          "mt -t $tapedrive rewind");
134        select(undef, undef, undef, $repos_sleep);
135        if ($position > 0) {
136            $tmp_status = &catch_tape_command($tapedrive,
137                                              "mt -t $tapedrive fsf $position")
138                unless ($tmp_status != 0);
139        }
140        if ($tmp_status == 0) {last;}
141        print STDERR "Error repositioning tape - trying again.\n";
142    }
143    if ($tmp_status != 0) {
144        warn "\n\nWARNING!!!!! Unable to reposition the tape.\n\n";
145    }
146
147    # Set counter
148    $counter_for_read_next_file = $position;
149
150    return $tmp_status;
151}
152
153# Subroutine to read a file from tape
154sub read_next_file {
155    local($tapedrive, *input_list) = @_;
156    if (!defined($tapedrive) && !defined(@input_list)) {
157        $tapedrive = "/dev/nrtape";
158    }
159
160    # Constants
161    $tape_block_size = 8192;
162#    $tape_sleep = 1;
163    $tape_sleep = 0;
164    $retry_sleep = 1;
165    $nretries = 4;
166
167
168    # Get next value from list if no tape drive
169    if (length($tapedrive) == 0) {
170        local($filename) = shift(@input_list);
171        if (length($filename) > 0) {
172           print "Reading header for file $filename\n";
173        }
174        return $filename;
175    }
176
177    # Create file counting variable if it does not exist and rewind tape
178    # drive. Note that the rewind implicitly opens the TAPEHANDLE
179    if (!defined($counter_for_read_next_file)) {
180        print STDERR "Rewinding tape drive $tapedrive\n";
181        &execute_tape_command($tapedrive, "mt -t $tapedrive rewind");
182        $counter_for_read_next_file = 0;
183    }
184
185    # Create the filename
186    local($cur_file_number) = $counter_for_read_next_file;
187    $counter_for_read_next_file++;
188    local($filename) = "$tmpdir/datafile_".$cur_file_number;
189
190    # Try reading from the tape drive. We will try repeatedly if necessary
191    print STDERR "Retrieving file $filename from drive $tapedrive\n";
192    local($status, $oldsep, $tapemessage);
193    foreach $retryloop (0..$nretries-1) {
194
195        # Sleep for a moment, then read from tape drive (don't ask me why,
196        # it just works!)
197        select(undef, undef, undef, $tape_sleep);
198        $oldsep = $/; undef($/);
199        close(TAPEHANDLE);
200        $status = 1;
201        if (open(TP, "dd if=$tapedrive of=$filename bs=1000000 2>&1 |")) {
202           $tapemessage = <TP>;
203           close(TP);
204           $status = $?;
205           if (($status != 0) && (-z $filename) &&
206               ($tapemessage =~ /^Read error:\s*No space left on device/)) {
207              $status = 0;
208           }
209        }
210        open(TAPEHANDLE, $tapedrive);
211        $/ = $oldsep;
212
213        if ($status == 0) {last;}
214
215        # If we get to here then the read failed. Try to reposition the tape.
216        print STDERR "Error reading from tape - trying again.\n";
217        if (&set_tape_position($tapedrive, $cur_file_number) != 0) {last;}
218        $counter_for_read_next_file++;
219
220        # Sleep to let things settle after repositioning
221        select(undef, undef, undef, $retry_sleep);
222
223    }
224
225    # We've finished trying to read. Test to see if we failed or if
226    # we've reached the end of the tape.
227    if (($status!=0) || (-z $filename)) {
228        &remove_file($filename);
229        if ($status != 0) {
230            warn "\n\nWARNING!!!!! ".
231                "Error occurred while reading tape. Giving up.\n\n\n";
232            &cleanup_and_die("Not creating last minc file.\n", $status);
233        }
234        else {
235            print STDERR "End of tape.\n";
236        }
237        return "";
238    }
239    else {
240        return $filename;
241    }
242
243}
244
245# Subroutine to unpack a value from a string
246sub unpack_value {
247    local(*string, $offset, $type) = @_;
248
249    # Check for byte swapping
250    if ($Little_endian) {
251       if ($type !~ /^\s*([a-zA-Z])(\d+)?\s*$/) {
252          die "Unrecognized data type \"$type\" on little-endian machine.\n";
253       }
254       local($code) = $1;
255       local($number) = (defined($2) ? $2 : 1);
256       if (!defined($Unpack_codes{$code})) {
257          die "Unrecognized unpack code \"$code\" on little-endian machine.\n";
258       }
259       local($size) = $Unpack_codes{$code};
260       local($tempstring) = substr($string, $offset, $number * $size);
261       if ($size > 1) {
262          local($iloop);
263          local($max) = $number * $size;
264          for ($iloop=0; $iloop < $max; $iloop+=$size) {
265             substr($tempstring, $iloop, $size) =
266                reverse(substr($tempstring, $iloop, $size));
267          }
268       }
269       return unpack("$type", $tempstring);
270    }
271
272    # No byte swapping required
273    else {
274       return unpack("x$offset $type", $string);
275    }
276}
277
278# Subroutine to get a direction cosine from a vector, correcting for
279# magnitude and direction if needed (the direction cosine should point
280# along the positive direction of the nearest axis)
281sub get_dircos {
282    if (scalar(@_) != 3) {
283        die "Argument error in get_dircos\n";
284    }
285    local($xcos, $ycos, $zcos) = @_;
286
287    # Get magnitude
288    local($mag) = sqrt($xcos**2 + $ycos**2 + $zcos**2);
289    if ($mag <= 0) {$mag = 1};
290
291    # Make sure that direction cosine is pointing along positive axis
292    local($max) = $xcos;
293    if (&abs($ycos) > &abs($max)) {$max= $ycos;}
294    if (&abs($zcos) > &abs($max)) {$max= $zcos;}
295    if ($max < 0) {$mag *= -1;}
296
297    # Correct components
298    $xcos /= $mag;
299    $ycos /= $mag;
300    $zcos /= $mag;
301
302    return ($xcos, $ycos, $zcos);
303}
304
305# Subroutine to add an optional attribute to an array
306sub add_optional_attribute {
307    if (scalar(@_) != 4) {
308        die "Argument error in add_optional_attribute\n";
309    }
310    local(*array, $type, $name, $value) = @_;
311    local($key);
312
313    if ($type eq "string") {$key = "-sattribute";}
314    elsif ($type eq "double") {$key = "-dattribute";}
315    elsif ($type eq "none") {$key = "-attribute";}
316    else {die "Unknown type \"$type\" in add_optional_attribute\n";}
317
318    if (length($value) > 0) {
319        push(@array, $key, $name."=".$value);
320    }
321}
322
323# Subroutine to create a minc file
324sub create_mincfile {
325    if (scalar(@_) != 5) {
326        die "Argument error in create_mincfile";
327    }
328    local($mincfile, *file_list, *mincinfo, $image_list, $echo)
329        = @_;
330
331    # Sort images by z position
332    local($cur_image, %pos_to_image, @positions, $cur_slicepos);
333    foreach $cur_image (split($;, $image_list)) {
334        $cur_slicepos = $mincinfo{$cur_image, 'slicepos'};
335        if (!defined($pos_to_image{$cur_slicepos})) {
336            $pos_to_image{$cur_slicepos} = $cur_image;
337            push(@positions, $cur_slicepos);
338        }
339        else {
340           warn "Duplicate slice position: " .
341              "ignoring file $file_list{$cur_image}\n";
342        }
343    }
344    @positions = sort(numeric_order @positions);
345
346    # Get step size and get nslices and ordered list of images
347    local(@image_list, $slicestep, $nslices, $slicestart, $slicerange);
348    if (scalar(@positions) > 1) {
349
350        # Find smallest step size
351        $slicestep = $positions[1] - $positions[0];
352        local(@difflist, $diff);
353        foreach $i (2..$#positions) {
354            $diff = $positions[$i] - $positions[$i-1];
355            push(@difflist, $diff);
356            if (($diff < $slicestep) && ($diff > 0)) {
357                $slicestep = $diff;
358            }
359        }
360
361        # Find average step, accounting for stepping over multiple slices
362        # (based on rounded ratio to smallest step)
363        local($ndiffs) = 0;
364        local($totaldiff) = 0;
365        foreach $diff (@difflist) {
366           $totaldiff += $diff;
367           $ndiffs += int($diff/$slicestep + 0.5);
368        }
369        if ($ndiffs > 0) {
370           $slicestep = $totaldiff / $ndiffs;
371        }
372
373        # Work out number of slices and get the ordered list of images
374        if ($slicestep <= 0) {
375            $nslices = scalar(@positions);
376            $slicestep = 1;
377            @image_list = @pos_to_image{@positions};
378        }
379        else {
380            local($first_image, $last_image);
381            $first_image = $pos_to_image{$positions[0]};
382            $last_image = $pos_to_image{$positions[$#positions]};
383            $slicerange = $mincinfo{$last_image, 'slicepos'} -
384                $mincinfo{$first_image, 'slicepos'};
385            $nslices = int($slicerange / $slicestep + 1.5);
386            # Improve accuracy of $slicestep
387            $slicestep = $slicerange / ($nslices - 1);
388            foreach $i (0..$nslices) {
389                $image_list[$i] = '';
390            }
391            foreach $position (@positions) {
392                local($slice) =
393                    int(($position - $positions[0]) / $slicestep + 0.5);
394                @image_list[$slice] = $pos_to_image{$position};
395            }
396        }
397    }
398    else {
399        $slicestep = 1;
400        $nslices = 1;
401        @image_list = @pos_to_image{@positions};
402    }
403    $slicestart = $positions[0];
404
405    # Check for existence of file - create a new name if it exists
406    local($filebase);
407    if ($mincfile =~ /^(.*)\.mnc$/) {
408        $filebase = $1;
409    }
410    else {
411        $filebase = $mincfile;
412    }
413    local($version) = 1;
414    local(@glob);
415    while (scalar(@glob=<$mincfile*>) > 0) {
416        $version++;
417        $mincfile = $filebase."_version$version.mnc";
418        print STDERR "Minc file already exists. Trying name \"$mincfile\".\n";
419    }
420
421    # Set up rawtominc command
422
423    # Dimension info
424    local($nrows, $ncols, $orientation, $pixel_size);
425    local($xstep, $ystep, $zstep);
426    local($colstart, $rowstart);
427    $pixel_size = $mincinfo{'pixel_size'};
428    $nrows = $mincinfo{'height'};
429    $ncols = $mincinfo{'width'};
430    $xstep = $mincinfo{'colstep'};
431    $ystep = $mincinfo{'rowstep'};
432    $zstep = $slicestep;
433    $colstart = $mincinfo{'colstart'}+0;
434    $rowstart = $mincinfo{'rowstart'}+0;
435
436    # Orientation
437    local($xstart, $ystart, $zstart);
438    local(%world_axes);
439    local($orient_flag) = $mincinfo{'orientation'};
440    if ($orient_flag eq 'sagittal') {    # Sagittal
441        ($ystep, $zstep, $xstep) =
442            ($mincinfo{'colstep'}, $mincinfo{'rowstep'}, $slicestep);
443        ($ystart, $zstart, $xstart) = ($colstart, $rowstart, $slicestart);
444        $orientation = "-sagittal";
445        %world_axes = ('col', 'y', 'row', 'z', 'slc', 'x');
446    }
447    elsif ($orient_flag eq 'coronal') { # Coronal
448        ($xstep, $zstep, $ystep) =
449            ($mincinfo{'colstep'}, $mincinfo{'rowstep'}, $slicestep);
450        ($xstart, $zstart, $ystart) = ($colstart, $rowstart, $slicestart);
451        $orientation = "-coronal";
452        %world_axes = ('col', 'x', 'row', 'z', 'slc', 'y');
453    }
454    else {                      # Transverse
455        ($xstep, $ystep, $zstep) =
456            ($mincinfo{'colstep'}, $mincinfo{'rowstep'}, $slicestep);
457        ($xstart, $ystart, $zstart) = ($colstart, $rowstart, $slicestart);
458        $orientation = "-transverse";
459        %world_axes = ('col', 'x', 'row', 'y', 'slc', 'z');
460    }
461    local(@dircos_options) = ();
462    foreach $axis ('col', 'row', 'slc') {
463        if (defined($mincinfo{"${axis}_dircos_x"})) {
464            push(@dircos_options,"-${world_axes{$axis}}dircos");
465            foreach $component ('x', 'y', 'z') {
466                push(@dircos_options, $mincinfo{"${axis}_dircos_$component"});
467            }
468        }
469    }
470
471    # Optional attributes
472    local(@optional_attributes) = ();
473    &add_optional_attribute(*optional_attributes, "double",
474                            "acquisition:repetition_time", $mincinfo{'tr'});
475    &add_optional_attribute(*optional_attributes, "double",
476                            "acquisition:echo_time", $mincinfo{'te',$echo});
477    &add_optional_attribute(*optional_attributes, "double",
478                            "acquisition:inversion_time", $mincinfo{'ti'});
479    &add_optional_attribute(*optional_attributes, "double",
480                            "acquisition:flip_angle", $mincinfo{'mr_flip'});
481    &add_optional_attribute(*optional_attributes, "string",
482                            "acquisition:scanning_sequence",
483                            $mincinfo{'scan_seq'});
484    &add_optional_attribute(*optional_attributes, "string",
485                            "study:study_id", $mincinfo{'exam'});
486    &add_optional_attribute(*optional_attributes, "string",
487                            "study:acquisition_id", $mincinfo{'series'});
488    &add_optional_attribute(*optional_attributes, "string",
489                            "study:start_time", $mincinfo{'start_time'});
490    &add_optional_attribute(*optional_attributes, "string",
491                            "study:institution", $mincinfo{'institution'});
492    &add_optional_attribute(*optional_attributes, "string",
493                            "patient:birthdate",
494                            $mincinfo{'patient_birthdate'});
495    &add_optional_attribute(*optional_attributes, "double",
496                            "patient:age", $mincinfo{'patient_age'});
497    &add_optional_attribute(*optional_attributes, "string",
498                            "patient:sex", $mincinfo{'patient_sex'});
499    &add_optional_attribute(*optional_attributes, "string",
500                            "patient:identification", $mincinfo{'patient_id'});
501
502    # Check for history
503    if (length($history) > 0) {
504       &add_optional_attribute(*optional_attributes, "string",
505                               ":history", $mincinfo{'history'});
506    }
507
508    # GE specific stuff
509    local(@ge_specific_attributes);
510    @ge_specific_attributes = ();
511    &add_optional_attribute(*ge_specific_attributes, "string",
512                            "ge_mrimage:pseq", $mincinfo{'ge_pseq'});
513    &add_optional_attribute(*ge_specific_attributes, "string",
514                            "ge_mrimage:pseqmode", $mincinfo{'ge_pseqmode'});
515
516    # Dicom specific stuff
517    local(@dicom_specific_attributes);
518    @dicom_specific_attributes = ();
519    local(@elements) = sort(grep(/^dicom_/, keys(%mincinfo)));
520    local($element);
521    foreach $element (@elements) {
522       &add_optional_attribute(*dicom_specific_attributes, "string",
523                              $element, $mincinfo{$element});
524    }
525
526    # Run rawtominc with appropriate arguments
527    $| = 1;
528    local(@typeinfo);
529    if ($mincinfo{'pixel_size'} == 1) {
530       @typeinfo = ("-byte", "-unsigned", "-range", "0", "255");
531    }
532    else {
533       @typeinfo = ("-short", "-signed",  "-range", "0", "4095");
534    }
535    local(@minccommand) =
536        ("rawtominc", $mincfile, $nslices, $nrows, $ncols, "-noclobber",
537         "-scan_range", @typeinfo, $orientation,
538         "-xstep", $xstep, "-ystep", $ystep, "-zstep", $zstep,
539         "-xstart", $xstart, "-ystart", $ystart, "-zstart", $zstart,
540         @dircos_options,
541         "-mri",
542         "-attribute", "patient:full_name=".$mincinfo{'patient_name'},
543         @optional_attributes,
544         @ge_specific_attributes,
545         @dicom_specific_attributes
546         );
547    open(MINC, "|-") || exec(@minccommand);
548
549    # Get keys for machine specific file info
550    local(@specific_keys, %specific_keys, $key);
551    undef(%specific_keys);
552    foreach $key (keys(%mincinfo)) {
553        $key =~ /$;specific$;([^$;]*)$/;
554        $specific_keys{$1} = '';
555    }
556    @specific_keys = keys(%specific_keys);
557
558    # Loop through slices
559    local($cur_file, %specific_file_info, $cur_image);
560    foreach $slice (0..$nslices-1) {
561
562        # Get image number
563        $cur_image = $image_list[$slice];
564
565        # Get file name
566        $cur_file = $file_list{$cur_image};
567
568        # Print log message
569        if (length($cur_file) > 0) {
570            print STDERR
571                "Getting data from file \"$cur_file\" for slice $slice.\n";
572        }
573        else {
574            print STDERR "Using blank image for slice $slice.\n";
575        }
576
577        # Set up variables for getting the image
578        $image_data_len = $nrows * $ncols * $pixel_size;
579
580        # Read the image (if not defined, create a blank image)
581        if (length($cur_file) <= 0) {
582            $image_data = pack("x$image_data_len",());
583        }
584        else {
585            # Get machine specific file info
586            undef(%specific_file_info);
587            foreach $key (@specific_keys) {
588                $specific_file_info{$key} =
589                    $mincinfo{$cur_image, 'specific', $key};
590            }
591
592            local($image_cmd) =
593                &get_image_cmd($cur_file, *specific_file_info);
594            if ($Little_endian) {
595               $image_cmd .= " | byte_swap";
596            }
597            open(GEDAT, $image_cmd . " |");
598            read(GEDAT, $image_data, $image_data_len);
599            local($dummy, $nread);
600            do {
601               $nread = read(GEDAT, $dummy, 8192);
602            } while ($nread == 8192);
603            close(GEDAT);
604            if ($? != 0) {
605                warn("Error or signal while reading image.\n");
606                if ($ignore_image_errors) {
607                    warn "Using blank image instead.\n";
608                    $image_data = pack("x$image_data_len",());
609                }
610                else {
611                    &cleanup_and_die("Quitting.\n", $?);
612                }
613            }
614            if (length($image_data) < $image_data_len) {
615                warn "Error reading image from \"$cur_file\"\n";
616                if ($ignore_image_errors) {
617                    warn "Using blank image instead.\n";
618                    $image_data = pack("x$image_data_len",());
619                }
620                else {
621                    close(MINC);
622                    return;
623                }
624            }
625            elsif (length($image_data) > $image_data_len) {
626               warn "More image data than expected - truncating.\n";
627               $image_data = substr($image_data, 0, $image_data_len);
628            }
629        }
630
631        # Write out the image
632        print MINC $image_data;
633    }
634
635    # Close the minc file
636    close(MINC);
637    if ($? != 0) {
638        &cleanup_and_die("Error or signal while writing image.\n", $?);
639    }
640
641    return $mincfile;
642}
643
644# Routine to add to the list file
645sub save_list_info {
646    local($listfile, $mincfile, *mincinfo, $image_list, $echo) = @_;
647    local($pos) = $mincinfo{'tape_position'};
648    local($pat_name) = '"'.substr($mincinfo{'patient_name'},0,30).'"';
649    local($exam) = $mincinfo{'exam'};
650    local($acq) = $mincinfo{'series'};
651    local($nslc) = scalar(split($;, $image_list));
652    local($date) = '"'.substr($mincinfo{'start_time'},0,30).'"';
653    local($inst) = '"'.substr($mincinfo{'institution'},0,30).'"';
654    open(LIST, ">>$listfile") ||
655        warn "Error writing to list file \"$listfile\" ($!)\n";
656    write LIST;
657    close(LIST);
658    format LIST =
659@>>>> @<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< @>>>>>>> @>> @> @>> @<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< @<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
660$pos  $pat_name                        $exam    $acq $echo $nslc $date $inst
661.
662}
663
664# Routine to get argument values
665sub get_arguments {
666
667    # Usage line
668    $0 =~ /([^\/]+)$/;
669    local($prog) = $1;
670    local($usage) =
671        "Usage: $prog <outputdir> [<file1> <file2> ...] [<options>]\n";
672
673    # Set default values
674    local($outputdir) = undef;
675    local($tapedrive) = '';
676    local($listfile) = '';
677    local($nominc) = 0;
678    local($tape_position) = 0;
679    local($tape_end) = -1;
680    local(@input_list) = ();
681    local($do_compression) = 0;
682    local($need_diskfiles) = 0;
683    local($ignore_image_errors) = 0;
684    local($inputdir) = '';
685
686    # Loop through arguments
687    while (@_) {
688        $_ = shift;
689        if (/^-list$/) { $listfile = shift;}
690        elsif (/^-tape$/) { $tapedrive = shift;}
691        elsif (/^-notape$/) {$need_diskfiles = 1;}
692        elsif (/^-inputdir$/) {$inputdir = shift;}
693        elsif (/^-nominc$/) { $nominc = 1;}
694        elsif (/^-first$/) { $tape_position = shift;}
695        elsif (/^-last$/) {$tape_end = shift;}
696        elsif (/^-compress$/) {$do_compression = 1;}
697        elsif (/^-nocompress$/) {$do_compression = 0;}
698        elsif (/^-ignore_image_errors$/) {$ignore_image_errors = 1;}
699        elsif (/^-h(|elp)$/) {
700            die
701"Command-specific options:
702 -list:\t\t\tSpecify a file for listing contents of tape
703 -tape:\t\t\tSpecify an input tape drive (default=/dev/nrtape if no files given)
704 -notape:\t\tDo not try to use the tape drive
705 -inputdir:\t\tSpecify directory from which input file names should be read
706 -nominc:\t\tDo not produce output minc files
707 -first:\t\tSpecify a starting tape position (# files to skip)
708 -last:\t\t\tSpecify an ending tape position
709 -compress:\t\tCompress the output minc files
710 -nocompress:\t\tDo not compress the output minc files (default)
711 -ignore_image_errors:\tIgnore errors when reading images
712Generic options for all commands:
713 -help:\t\t\tPrint summary of comand-line options and abort
714
715$usage
716";
717        }
718        elsif (/^-./) {
719            die "Unknown option \"$_\"\n";
720        }
721        else {
722            if (!defined($outputdir)) {
723                $outputdir = $_;
724            }
725            else {
726                push(@input_list, $_);
727            }
728        }
729    }
730
731    # Get input file names from inputdir if needed
732    if ((length($inputdir) > 0) && (scalar(@input_list) > 0)) {
733       die "Do not use -inputdir option with input files\n";
734    }
735    if (length($inputdir) > 0) {
736       opendir(DIR, $inputdir) ||
737          die "Error opening input directory \"$inputdir\": $!\n";
738       @input_list = sort(grep(/^[^\.]/, readdir(DIR)));
739       closedir(DIR);
740       local($file);
741       foreach $file (@input_list) {
742          $file = "$inputdir/$file";
743       }
744    }
745
746    # Check expected values
747    if ($need_diskfiles && (scalar(@input_list) == 0)) {
748       die "Please specify disk files.\n";
749    }
750    if ((length($tapedrive) > 0) && (scalar(@input_list) > 0)) {
751        die "You cannot specify both a tape drive and a file list.\n";
752    }
753    elsif ((length($tapedrive) <= 0) && (scalar(@input_list) <= 0)) {
754        $tapedrive = "/dev/nrtape";
755    }
756    if (!defined($outputdir)) {
757        die $usage;
758    }
759    if ($tape_position < 0) {
760        die "Tape position must be >= 0\n";
761    }
762    if (!defined($tape_end)) {$tape_end = -1;}
763
764    # Return values
765    return($outputdir, $tapedrive, $listfile, $nominc,
766           $tape_position, $tape_end, $do_compression,
767           $ignore_image_errors, @input_list);
768}
769
770# Subroutine to do all the work - loops through files collecting info,
771# then calls routine to create minc file.
772# Because this was the main program, it uses global variables
773sub mri_to_minc {
774
775    $| = 1;
776
777    # Save arguments for history
778    @saved_args = @ARGV;
779
780    # Get arguments
781    ($outputdir, $tapedrive, $listfile, $nominc,
782     $tape_position, $tape_end, $do_compression,
783     $ignore_image_errors, @input_list) =
784         &get_arguments(@_);
785
786    # Save history
787    chop($history = `date`);
788    $0 =~ /([^\/]+)$/;
789    $history .= ">>>> $1";
790    foreach $file (@input_list[0..$#input_list]) {
791       $bad_args{$file} = 1;
792    }
793    foreach $arg (@saved_args) {
794       if (! $bad_args{$arg}) {
795          $history .= " $arg";
796       }
797    }
798    if (scalar(@input_list) > 0) {
799       $history .= " $input_list[0]";
800    }
801    if (scalar(@input_list) > 1) {
802       $history .= " ... $input_list[$#input_list]";
803    }
804    $history .= "\n";
805
806    # Check that we can write to the output directory
807    if (! -d $outputdir) {
808       die "Directory \"$outputdir\" does not exist.\n";
809    }
810    if (! -w $outputdir) {
811       die "Cannot write to \"$outputdir\".\n";
812    }
813
814    # Should we delete the files?
815    $delete_files = (length($tapedrive) > 0);
816
817    # Create a temporary directory
818    &create_tmpdir("mri_to_minc_$$");
819
820    # Check if list file already exists
821    if (length($listfile) > 0) {
822        if (-e $listfile) {
823            die "List file \"$listfile\" already exists.\n";
824        }
825        open(LIST_START, ">$listfile") ||
826            die "Unable to write to list file \"$listfile\" ($!)\n";
827        write LIST_START;
828        close(LIST_START);
829        format LIST_START =
830@>>>> @<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< @>>>>>>> @>> @> @>> @<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< @<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
831'File'  'Patient Name' 'Exam' 'Acq' 'E#' 'Slc' 'Date' 'Institution'
832
833.
834    }
835
836    # Rewind and initialize the tape
837    if (length($tapedrive) > 0) {
838        &initialize_tape_drive($tapedrive);
839        if ($tape_position > 0) {
840            &set_tape_position($tapedrive, $tape_position);
841        }
842    }
843
844    # Loop through files on tape
845    $keep_looping = 1;
846    while ($keep_looping) {
847
848        # Save current tape position
849        $tape_position = &get_tape_position;
850
851        # Check for reaching last tape file requested
852        if (($tape_end >= 0) && ($tape_position > $tape_end)) {
853            warn "Reached last file requested (file number $tape_end).\n";
854            $keep_looping = 0;
855        }
856
857        # Get next file
858        if ($keep_looping) {
859            $nextfile = &read_next_file($tapedrive, *input_list);
860            if (length($nextfile) <= 0) {
861                $keep_looping = 0;
862            }
863        }
864
865        # Read in headers
866        if ($keep_looping) {
867            undef(%file_info);
868            undef(%specific_file_info);
869            if (&read_file_info($nextfile, *file_info, *specific_file_info)) {
870                warn "Error reading file \"$nextfile\". Skipping to next.\n";
871                next;
872            }
873
874            # Get interesting values
875            $cur_numechos = $file_info{'numechos'};
876            $cur_exam = $file_info{'exam'};
877            $cur_series = $file_info{'series'};
878            $cur_echo = $file_info{'echo'};
879            $cur_image = $file_info{'image'};
880            if ($cur_numechos > 1) {
881               $cur_image = $cur_image * $cur_numechos +
882                  (($cur_echo <= $cur_numechos) ? $cur_echo : 1);
883            }
884            $cur_width = $file_info{'width'};
885            $cur_height = $file_info{'height'};
886            $cur_slicepos = $file_info{'slicepos'};
887        }
888
889        # Check if number exam or series has changed or if this is the
890        # last file
891        if ((scalar(keys(%file_list)) > 0) &&
892            (!$keep_looping ||
893             ($mincinfo{'exam'} ne $cur_exam) ||
894             ($mincinfo{'series'} != $cur_series))) {
895
896            # Loop through echos
897            @echos = sort(numeric_order keys(%echo_list));
898            if (scalar(@echos) > $mincinfo{'numechos'}) {
899                $mincinfo{'numechos'} = scalar(@echos);
900            }
901            foreach $echo (@echos) {
902
903                # Create minc file
904                ($patient_name = $mincinfo{'patient_name'}) =~
905                    tr/a-zA-Z0-9_/_/cs;             # Use only legal chars
906                $patient_name =~ tr/A-Z/a-z/;         # Lowercase
907                $patient_name =~ s/_*(.*[^_])_*$/$1/; # Remove _ on ends
908                $mincfile = "$outputdir/".$patient_name."_".
909                    $mincinfo{'exam'}."_".$mincinfo{'series'};
910                if ($mincinfo{'numechos'} > 1) {
911                    $mincfile .= "_e$echo";
912                }
913                $mincfile .= "_mri.mnc";
914                if (length($listfile) > 0) {
915                    &save_list_info($listfile, $mincfile, *mincinfo,
916                                    $image_list{$echo}, $echo);
917                }
918                if (!$nominc) {
919                    print STDERR "Creating minc file \"$mincfile\"\n";
920                    $mincfile =
921                       &create_mincfile($mincfile, *file_list, *mincinfo,
922                                        $image_list{$echo}, $echo);
923                    if ($do_compression) {
924                       print STDERR "Compressing file \"$mincfile\"\n";
925                       system("gzip $mincfile");
926                    }
927                }
928                else {
929                    print STDERR "Not creating minc file \"$mincfile\"\n";
930                }
931            }
932
933            # Delete files (if needed) and reset variables
934            if ($delete_files) {
935                &remove_file(values(%file_list));
936            }
937            undef(%file_list);
938            undef(%echo_list);
939            undef(%image_list);
940            undef(%mincinfo);
941        }
942
943        # Break out here if stopping loop
944        if (!$keep_looping) {next;}
945
946        # If first file, then save appropriate info
947        if (scalar(keys(%file_list)) <= 0) {
948            $mincinfo{'history'} = $history;
949            $mincinfo{'tape_position'} = $tape_position;
950            $mincinfo{'numechos'} = $cur_numechos;
951            $mincinfo{'exam'} = $cur_exam;
952            $mincinfo{'series'} = $cur_series;
953            $mincinfo{'width'} = $cur_width;
954            $mincinfo{'height'} = $cur_height;
955            $mincinfo{'pixel_size'} = $file_info{'pixel_size'};
956            $mincinfo{'patient_name'} = $file_info{'patient_name'};
957            $mincinfo{'orientation'} = $file_info{'orientation'};
958            $mincinfo{'tr'} = $file_info{'tr'};
959            $mincinfo{'ti'} = $file_info{'ti'};
960            $mincinfo{'mr_flip'} = $file_info{'mr_flip'};
961            $mincinfo{'scan_seq'} = $file_info{'scan_seq'};
962            $mincinfo{'ge_pseq'} = $file_info{'ge_pseq'};
963            $mincinfo{'ge_pseqmode'} = $file_info{'ge_pseqmode'};
964            $mincinfo{'start_time'} = $file_info{'start_time'};
965            $mincinfo{'institution'} = $file_info{'institution'};
966            $mincinfo{'patient_birthdate'} = $file_info{'patient_birthdate'};
967            $mincinfo{'patient_age'} = $file_info{'patient_age'};
968            $mincinfo{'patient_sex'} = $file_info{'patient_sex'};
969            $mincinfo{'patient_id'} = $file_info{'patient_id'};
970            $mincinfo{'colstep'} = $file_info{'colstep'};
971            $mincinfo{'rowstep'} = $file_info{'rowstep'};
972            $mincinfo{'colstart'} = $file_info{'colstart'};
973            $mincinfo{'rowstart'} = $file_info{'rowstart'};
974            $mincinfo{'col_dircos_x'} = $file_info{'col_dircos_x'};
975            $mincinfo{'col_dircos_y'} = $file_info{'col_dircos_y'};
976            $mincinfo{'col_dircos_z'} = $file_info{'col_dircos_z'};
977            $mincinfo{'row_dircos_x'} = $file_info{'row_dircos_x'};
978            $mincinfo{'row_dircos_y'} = $file_info{'row_dircos_y'};
979            $mincinfo{'row_dircos_z'} = $file_info{'row_dircos_z'};
980            $mincinfo{'slc_dircos_x'} = $file_info{'slc_dircos_x'};
981            $mincinfo{'slc_dircos_y'} = $file_info{'slc_dircos_y'};
982            $mincinfo{'slc_dircos_z'} = $file_info{'slc_dircos_z'};
983
984            # Save dicom element information
985            local(@elements) = sort(grep(/^dicom_/, keys(%file_info)));
986            local($element);
987            foreach $element (@elements) {
988               $mincinfo{$element} = $file_info{$element};
989            }
990
991        }
992        else {
993            if (($cur_width != $mincinfo{'width'}) || #
994                ($cur_height != $mincinfo{'height'})) {
995                warn "Width or height do not match first file ".
996                    "for file $nextfile. Skipping to next.\n";
997                next;
998            }
999        }
1000
1001        # Save echo number
1002        if (!defined($mincinfo{'te', $cur_echo})) {
1003            $mincinfo{'te', $cur_echo} = $file_info{'te'};
1004        }
1005
1006        # Save info for all files
1007        foreach $key (keys(%specific_file_info)) {
1008            $mincinfo{$cur_image, 'specific', $key} =
1009                $specific_file_info{$key};
1010        }
1011        $mincinfo{$cur_image, 'slicepos'} = $cur_slicepos;
1012
1013        # Keep track of files, image numbers and echo numbers
1014        $file_list{$cur_image} = $nextfile;
1015        $echo_list{$cur_echo} = '';
1016        if (defined($image_list{$cur_echo})) {$image_list{$cur_echo} .= $; ;}
1017        $image_list{$cur_echo} .= $cur_image;
1018
1019    }
1020
1021    # Rewind the tape
1022    if (length($tapedrive) > 0) {
1023        &execute_tape_command($tapedrive, "mt -t $tapedrive rewind");
1024    }
1025
1026    &cleanup_and_die();
1027}
1028
1029