On Feb 25, 2007, at 2:52 PM, Scott White wrote:

> I am writing an application that reads medical imaging files (CT,  
> MRI, etc in DICOM format) and simply displays them. However, I am  
> having trouble using Gtk2::Gdk::Pixmap. I must use a pixmap since  
> DICOM images are usually 16-bit, and almost always grayscale.
[...]
> Unfortunately, I cannot make use of the much simpler pixbuf, since  
> it only supports 8-bit pixels.

Stop right there.  GdkPixmap's option for 16 on the bits-per-pixel  
has nothing to do with what you're doing here.

GdkPixmap is a wrapper for XPixmap, which is a server-side resource.  
+ 
The pixel data is at the server's resolution and bitdepth.  When you  
think of a 16-bit pixmap, you are usually talking about packed 16-bit 
+ 
RGB data, or possibly 16-bit indexed color.  This is very, very  
different from a 16-bit grayscale image.

A 16-bit grayscale image contains grayscale data with 64k possible  
shades -- 256 times as much information as 8-bit data.  You monitor  
can't show you this much detail, and your video card deals only with  
8 bits per channel (8 x 3 = 24, sometimes with an unused padding byte 
+ 
for 32), so in order to display this image, you *must* mangle it to 8-
+ 
bits-per-channel somehow.  The art it is how you choose to compress  
the dynamic range to make it visible.

The very simplest thing you can do is just show the most significant  
byte of each pixel.  That is:

        foreach pixel
                output = input >> 8

or

        $data = pack "C*", map { $_ >> 8 } unpack "S*", $raw;

(To display this in a native GdkPixbuf, you'd have to explode that to 
+ 
24-bit by replicating each pixel three times, which rather sucks,  
memory-wise.  The method draw_gray_image() on Gtk2::Gdk::Drawable  
would allow you not to have to expand your memory usage like this,  
but would require you to write your own image display widget.)

If the image uses its full dynamic range, you'll be able to see  
*something*.

However, it's often the case that 16-bit images don't use the full  
spectrum.  If you get a histogram, you may see hotspots in different  
parts of the spectrum.  You might also be losing the detail that you  
really need to see.  So, to be really fancy, you'd use different  
algorithms to map the 16-bit data to 8-bit.  That's where the fun  
comes in... gamma curves, sliders for the top and bottom of the  
selected portion of the dynamic range, etc, etc.


Here's a quick, very crappy example of a few different ways you can  
mangle the data.  Anything truly good will require more finesse than  
i was willing to put into 130 lines of code.  :-)


##################################
#!/usr/bin/perl -w
use strict;
use Gtk2 -init;
use Glib qw(:constants);

# input file is raw 16-bit grayscale image data at 512x512, from
# http://zentara.net/perlplay/DICOM/CT-MONO2-16-brain.raw
my $file = "CT-MONO2-16-brain.raw";
my $bits = 16;
my $width = 512;
my $height = 512;
open IN, $file or die "$file: $!\n";
my $expect_size = $width * $height * ($bits/8);
my $data;
my $n_read = sysread IN, $data, $expect_size;
die "File is the wrong size -- got $n_read, expected $expect_size\n"
    unless $n_read == $expect_size;


my $window = Gtk2::Window->new;
$window->signal_connect (destroy => sub { Gtk2->main_quit });

my $hbox = Gtk2::HBox->new;
$window->add ($hbox);

my $image = Gtk2::Image->new;
$image->set_size_request ($width, $height);
$hbox->add ($image);

my $vbox = Gtk2::VBox->new;
$hbox->add ($vbox);

foreach (
    { label => "out = low byte", func => \&low_byte },
    { label => "out = in >>  1", func => sub { shift_by (1) } },
    { label => "out = in >>  2", func => sub { shift_by (2) } },
    { label => "out = in >>  3", func => sub { shift_by (3) } },
    { label => "out = in >>  4", func => sub { shift_by (4) } },
    { label => "out = in >>  5", func => sub { shift_by (5) } },
    { label => "out = in >>  6", func => sub { shift_by (6) } },
    { label => "out = in >>  7", func => sub { shift_by (7) } },
    { label => "out = in >>  8", func => sub { shift_by (8) } },
    { label => "top 12 bits", func => sub { stretch ((1<<4)-1, (1<<16)
+-1) } },
    { label => "bottom 12 bits", func => sub { stretch (0, (1<<12)-1) 
+} },
    { label => "bottom 10 bits", func => sub { stretch (0, (1<<10)-1) 
+} },
    { label => "bottom 6 bits", func => sub { stretch (0, (1<<6)-1) } 
+},
    { label => "sqrt curve", func => sub { curve () } },
) {
    my $button = Gtk2::Button->new ($_->{label});
    $button->signal_connect (clicked => $_->{func});
    $vbox->pack_start ($button, FALSE, FALSE, 0);
}

$window->show_all;

Gtk2->main;

sub set_from_rgb {
    my $rgb = shift;
    my $pixbuf = Gtk2::Gdk::Pixbuf->new_from_data ($rgb, 'rgb', FALSE,
+ 8,
                                                   $width, $height,
                                                   $width * 3);
    $image->set_from_pixbuf ($pixbuf);
}


#
# Dead simple -- just mask off bits 9-16 and show only the low byte.
# This will result in interesting wrapping effects for any pixel >255.
#
sub low_byte {
    my $rgb = pack "C*",
                  map { ($_ & 0x00ff) x 3 }
                      unpack "S*", $data;

    set_from_rgb ($rgb);
}

#
# Like low_byte(), but shift right by n first.
#
sub shift_by {
    my ($n) = @_;
    my $rgb = pack "C*",
                  map { (($_ >> $n) & 0x00ff) x 3 }
                      unpack "S*", $data;

    set_from_rgb ($rgb);
}

#
# Attempt to do a linear scaling of the input pixels from the range
# [$min, $max] to [0,255].  This is implemented in a very slow manner.
+..
#
sub stretch {
    my ($min, $max) = @_;

    sub scale_one {
        my ($min, $max, $this) = @_;
        my $val = int (($this - $min) / ($max - $min) * 255);
        $val = 255 if $val > 255;
        $val = 0 if $val < 0;
        return 0x00ff & $val;
    }

    my $rgb = pack "C*",
                  map { (scale_one ($min, $max, $_)) x 3 }
                      unpack "S*", $data;

    set_from_rgb ($rgb);
}

#
# Apply a square root transformation to each pixel.
# The result is like looking at the high byte, but with more detail
# in the dark areas.
#
sub curve {

    my $rgb = pack "C*",
                  map { (int (255 * sqrt ($_ / 65535.0))) x 3 }
                      unpack "S*", $data;

    my $pixbuf = Gtk2::Gdk::Pixbuf->new_from_data ($rgb, 'rgb', FALSE,
+ 8,
                                                   $width, $height,
                                                   $width * 3);
    $image->set_from_pixbuf ($pixbuf);
}
[download]