Dicom.pm

=pod

=head1 NAME

PDL::IO::Dicom - a module for reading DICOM images.

=head1 DESCRIPTION

The PDL::IO::Dicom module enables reading 16-bit gray level Dicom
images into PDL. As Dicom is an extremely complex format, this module
can unfortunately not handle all different image types included in the
DICOM standard. One common format that is currently not supported is
the Papyrus format.

=head1 USAGE

    use PDL;
    use PDL::IO::Dicom;

    $img = rdcm("image.dcm");

=head1 SEE ALSO

L<PDL::IO::Dcm>

=head1 AUTHOR

Copyright (C) Dov Grobgeld <dov@imagic.weizmann.ac.il> 2002.
All rights reserved. There is no warranty. You are allowed
to redistribute this software / documentation under certain
conditions. For details, see the file COPYING in the PDL
distribution. If this file is separated from the PDL distribution,
the copyright notice should be included in the file.

=cut

package PDL::IO::Dicom;

use strict;
use warnings;
use PDL;
use PDL::Core;
use PDL::IO::Misc;

use Exporter;
our @ISA = qw( Exporter ); 
our @EXPORT = qw( rdcm );
our $VERSION = '2.098';

my $debug = 0;

my %element_table =
    (
     "0002,0000" => ["Group length","UL"],
     "0002,0001" => ["File Meta Information Version","OB"],
     "0002,0002" => ["Media Storage SOP Class UID","UI"],
     "0002,0003" => ["Media Storage SOP Instance UID","UI"],
     "0002,0010" => ["Transfer Syntax UID","UI"],
     "0002,0012" => ["Implementation Class UID","UI"],
     "0002,0013" => ["Implementation Version Name","SH"],
     "0002,0016" => ["Source Application Entity Title","AE"],
     "0008,0000" => ["Identifying Group Length", "UL"],
     "0008,0001" => ["Length to End"],
     "0008,0005" => ["Specific Character Set","CS"],
     "0008,0008" => ["Image Type","CS"],
     "0008,0010" => ["Recognition Code"],
     "0008,0012" => ["Instance Creation Date"],
     "0008,0013" => ["Instance Creation Time"],
     "0008,0014" => ["Instance Creator UID"],
     "0008,0016" => ["SOP Class UID","UI"],
     "0008,0018" => ["SOP Instance UID","UI"],
     "0008,0020" => ["Study Date","DA"],
     "0008,0021" => ["Series date","DA"],
     "0008,0022" => ["Acquisition Date","DA"],
     "0008,0023" => ["Image Date","DA"],
     "0008,0030" => ["Study Time","TM"],
     "0008,0031" => ["Series Time","TM"],
     "0008,0032" => ["Acquisition Time","TM"],
     "0008,0033" => ["Image Time","TM"],
     "0008,103e" => ["Series Description"],
     "0008,1030" => ["Study Description"],
     "0008,0050" => ["Accession Number","SH"],
     "0008,0060" => ["Modality","CS"],
     "0008,0070" => ["Manufacturer","LO"],
     "0008,0080" => ["Institution Name","LO"],
     "0008,0090" => ["Referring Physician's Name","PN"],
     "0008,1010" => ["Station Name","SH"],
     "0008,1030" => ["Study Description","LO"],
     "0008,103e" => ["Series Description","LO"],
     "0008,1060" => ["Name of Physician(s)","PN"],
     "0008,1090" => ["Manufacturers Model Name",],
     "0010,0000" => ["Patient Information Group"],
     "0010,0010" => ["Patient Name","PN"],
     "0010,0020" => ["Patient ID"],
     "0010,0030" => ["Patient Birth Date"],
     "0010,0040" => ["Patient's Sex"],
     "0010,1010" => ["Patient Age"],
     "0010,1030" => ["Patient Weight"],
     "0018,0000" => ["Acquisition Gnformation Group"],
     "0018,0010" => ["Contrast/Bolus Agent"],
     "0018,0020" => ["Scanning Sequence"],
     "0018,0022" => ["Scan Options"],
     "0018,0050" => ["Slice thickness"],
     "0018,1050" => ["Spatial Resolution", "lO"],
     "0008,1060" => ["Name of Physician(s) Reading Study"],
     "0008,1070" => ["Operator's Name"],
     "0018,0080" => ["Repetition Time"],
     "0018,0081" => ["Echo Time"],
     "0018,0082" => ["Inversion Time"],
     "0018,0083" => ["Number of Averages"],
     "0018,0084" => ["Imaging Frequency"],
     "0018,0085" => ["Imaged Nucleus"],
     "0018,0086" => ["Echo Number(s)"],
     "0018,0087" => ["Magnetic Field Strength"],
     "0018,0088" => ["Spacing Between Slices"],
     "0018,0089" => ["Number of Phase Encoding Steps"],
     "0018,0090" => ["Data Collection Diameter"],
     "0018,0091" => ["Echo Train Length"],
     "0018,0093" => ["Percent Sampling"],
     "0018,0094" => ["Percent Phase Field of View"],
     "0018,0095" => ["Pixel Bandwidth"],
     "0018,1088" => ["Heart Rate","US"],
     "0018,1090" => ["Cardiac Number of Images","US"],
     "0018,1094" => ["Trigger Window","US"],
     "0018,1100" => ["Reconstruction Diameter"],
     "0018,1314" => ["Flip Angle"],
     "0018,1315" => ["Variable Flip Angle Flag"],
     "0018,1316" => ["SAR"],
     "0020,0000" => ["Relationship information group"],
     "0020,000d" => ["Study Instance UID"],
     "0020,000e" => ["Series Instance UID"],
     "0020,0010" => ["Study ID"],
     "0020,0011" => ["Series Number"],
     "0020,0012" => ["Acquisition Number"],
     "0020,0013" => ["Image Number"],
     "0020,0014" => ["Isotope Number"],
     "0020,0015" => ["Phase Number"],
     "0020,0016" => ["Interval Number"],
     "0020,0017" => ["Time Slot Number"],
     "0020,0018" => ["Angle Number"],
     "0020,0020" => ["Patient Orientation"],
     "0020,0022" => ["Overlay Number"],
     "0020,0024" => ["Curve Number"],
     "0020,0026" => ["LUT Number"],
     "0020,0030" => ["Image Position"],
     "0020,0032" => ["Image Position (Patient)"],
     "0020,0035" => ["Image Orientation"],
     "0020,0037" => ["Image Orientation (Patient)"],
     "0020,0050" => ["Location"],
     "0020,0052" => ["Frame of Reference UID"],
     "0020,0060" => ["Laterality"],
     "0020,1002" => ["Images in acqusition"],
     "0020,1040" => ["Position Reference Indicator"],
     "0020,1041" => ["Slice Location"],
     "0028,0000" => ["Relationship information group"],
     "0028,0002" => ["Samples per Pixel"],
     "0028,0004" => ["Photometric Interpretation"],
     "0028,0005" => ["Image Dimensions"],
     "0028,0010" => ["Rows", "US"],
     "0028,0011" => ["Columns", "US"],
     "0028,0030" => ["Pixel Spacing"],
     "0028,0100" => ["Bits Allocated","US"],
     "0028,0101" => ["Bits Stored", "US"],
     "0028,0102" => ["High Bit","US"],
     "0028,0103" => ["Pixel Representation","US"],
     "0028,1052" => ["Rescale Intercept"],
     "0028,1053" => ["Rescale Slope"],
     "0033,1002" => ["IMGF"],
     "7fe0,0000" => ["Pixel Data Information Group", "UL"],
     "7fe0,0010" => ["Pixel Data","OW"]
     );

=head1 FUNCTIONS

=head2 rdcm

=for ref

Read a dicom image.

=for usage

 $dcm = rdcm("filename")

=cut
  
sub rdcm {
    my $file = shift;
    my $options = shift;
    my $do_print_info = 0;
    my (%info, %bin);
    my $do_raw = 0;   # Only for debugging
    my ($rescale_intercept, $rescale_slope) = (0, 1);
    my $do_explicit = 0;
    my $do_guess_endian = 1;

    # options
    if ($options) {
	if (defined $options->{do_print_info}) {
	    $do_print_info = $options->{do_print_info};
	}
    }

    open(IN, $file) || die "Failed opening image $file!\n";
    binmode IN;
    # read the whole image
    my $header;
    read(IN, $header, -s $file);
    
    # File preamble - a fixed 128 byte field
    my $hpos = 0x80;
    
    # Next four bytes should be DICM
    if (substr($header, $hpos, 4) ne 'DICM') {
	die "This is not a DICM file!\n";
    }
    $hpos+= 4;

    # Precheck if the first entry has explicit vr. Unfortunately this
    # is not enough to determine if the file always has explicit value
    # representation.
    if (substr($header, $hpos+4, 2)=~ /[A-Z]{2}/) {
	$do_explicit++;
    }

    while($hpos < length($header)) {
	my $is_binary = 0;
	my $groupword = unpack("v", substr($header, $hpos, 2));     $hpos+=2;
	my $elementword = unpack("v", substr($header, $hpos, 2));   $hpos+=2;
	my $value_rep = substr($header, $hpos, 2);
	my $key = sprintf("%04x,%04x", $groupword, $elementword);
	my ($lookup) = $element_table{$key};
	my $override_vr = 0;
	
	# Check for explicit value representation. There must be a different
	# way to figure this out, but I still haven't figured out how!
	if ($value_rep =~ /[A-Z][A-Z]/) {
	    $hpos+=2;
	}
	else {
	    $override_vr++;
	    $value_rep = $lookup->[1] || "UN";
	}

	my $elementlength;
	if (substr($header, $hpos, 4) eq "IMGF") {
	    die "No support for IMGF files at the moment!\n";
	}

	# The following calculation agrees with dicom3tools
	if ($override_vr) {
	    $elementlength = unpack("V", substr($header, $hpos, 4));
	    $hpos+= 4;
	}
	elsif (grep($value_rep eq $_, qw(OB OW SQ UN))) {
	    # Long length
	    $is_binary = 1;
	    $hpos+= 2;         # Always zero
	    $elementlength = unpack("V", substr($header, $hpos, 4));
	    $hpos+= 4;
	}
	else {
	    # Short length
	    $elementlength = unpack("v", substr($header, $hpos, 2));
	    $hpos+= 2;
	}

	my ($descr) = "";
	my $contents = substr($header, $hpos, $elementlength) . "";
	
	($descr) = @$lookup if $lookup;
	
	# recode contents
	if    ($value_rep eq "UL") {   $contents = unpack("V", $contents); }
	elsif ($value_rep eq "US") {   $contents = unpack("v", $contents); }
	elsif ($value_rep eq "TM") {   $contents = clean_time($contents); }
	elsif ($value_rep eq "DA") {   $contents = clean_date($contents); }
	elsif (!$is_binary) {
	    $contents=~ s/\s+$//;
	}
	
	# store the contents
	if ($key eq "7fe0,0010") {  # Pixel data
	    $bin{$descr} = $contents;
	}
	else {
	    $info{$descr} = $contents if $descr;
	}
	
	if ($do_print_info) {
	    $contents = "" unless $contents;
	    $contents = "<<Binary>>" if $is_binary;
	    $contents=~ tr/\0-\037\200-\377//d;
	    $contents= substr($contents, 0, 40) unless $do_raw;
	    $contents=~ s/[\0-\037\200]/?/g;
	    $value_rep=~ s/[\0-\037]/?/g;
	    printf STDERR "%04x> %04x,%04x (%04x,$value_rep) %-30s : %s\n",
	    $hpos, $groupword, $elementword, $elementlength, $descr, $contents;
	}
	
	$hpos+= $elementlength;
    }
    
    my($width, $height) = ($info{Columns}, $info{Rows});
    my($bs) = $info{"Bits Allocated"}/8;
    my $img = $bin{"Pixel Data"};
    # The following logics works on Intel endian
    my $do_toggle_endian = $info{"Pixel Representation"};
    $rescale_intercept = $info{"Rescale Intercept"} if defined $info{"Rescale Intercept"};
    $rescale_slope = $info{"Rescale Slope"} if defined $info{"Rescale Slope"};

    my $pdl;
    
    # Create a pdl from the raw data
    if ($bs == 2) {
	$pdl = zeroes(ushort, $width,$height);
	my $hdr = $pdl->gethdr;
	$pdl->make_physical();

	# Store the pixel data
	${$pdl->get_dataref()} = $img;

        # Guess endian
        if ($do_guess_endian) {
	    # Compare spread of high byte with low byte
	    my $high_byte = ($pdl>>8)->byte;
	    my $low_byte = ($pdl & 0xff)->byte;
	    my $max_high = $high_byte->max;
	    my $max_low = $low_byte->max;

	    # print STDERR "max_high max-low = $max_high $max_low\n";

	    # The following might need to be adjusted on different
	    # architectures.
	    $do_toggle_endian = $max_high > $max_low;
        }
    
        # Endian swap
        if ($do_toggle_endian) {
            bswap2($pdl);
	}

        # Rescale and convert to double
        if ($rescale_intercept != 0
	    || $rescale_slope != 1) {
	    print STDERR "scaling with $rescale_slope and $rescale_intercept\n";
	    $pdl = 1.0*($pdl * $rescale_slope) + $rescale_intercept;
	}
    } else {
	die "Sorry! PDL::IO::Dicom currently only supported DICOM with bs=2. bs = $bs\n";
    }

    # Store the info in the pdl header
    $pdl->sethdr(\%info);
    
    return $pdl;
}

sub clean_time {
    my $time = shift;
    my ($hour, $min, $sec, $sec_frac);
    if ($time=~ /(\d+):(\d+):(\d+)\.?(\d*)/) {
	($hour,$min,$sec,$sec_frac) = ($1,$2,$3,$4);
    }
    elsif ($time=~ /(\d\d)(\d\d)(\d\d)\.?(\d+)/) {
	($hour,$min,$sec,$sec_frac) = ($1,$2,$3,$4);
    }

    if (defined $hour) {
	$time = sprintf("%02d:%02d:%02d", $hour, $min, $sec);
	$time .= ".$sec_frac" if $sec_frac;
    }
    return $time;
}

sub clean_date {
    my $date = shift;
    $date=~ s/(\d\d\d\d)(\d\d)(\d\d)/$1-$2-$3/;
    return $date;
}

=head1 AUTHOR

Copyright (C) Dov Grobgeld <dov@imagic.weizmann.ac.il> 1997.
All rights reserved. There is no warranty. You are allowed
to redistribute this software / documentation under certain
conditions. For details, see the file COPYING in the PDL
distribution. If this file is separated from the PDL distribution,
the copyright notice should be included in the file.

=cut

1;