Simple display of DICOM format file image

Question

I am trying to read a simple dicom file named "CR-MONO1-10-chest" from http://www.barre.nom.fr/medical/samples/ which is a 440x440 sized image.

As expained on http://www.dclunie.com/medical-image-faq/html/part1.html , the image data is at the end of the file:

In other words, if an image is 256 by 256, uncompressed, and 12-16 bits deep (and hence usually, though not always, stored as two bytes per pixel), then we all know that the file is going to contain 256*256*2=131072 bytes of pixel data at the end of the file. If the file is say 145408 bytes long, as all GE Signa 3X/4X files are for example, then you need to skip 14336 bytes of header before you get to the data. Presume row by row starting with top left hand corner raster order, try both alternatives for byte order, deal with the 16 to 8 bit windowing problem, and very soon you have your image on the screen of your workstation.

A figure on http://people.cas.sc.edu/rorden/dicom/index.html also shows that the image data is located at the end of the file.

I am using following code to read and display image in this file:

(define in (open-input-file "CR-MONO1-10-chest" #:mode 'binary))
(define size (* 2 440 440))                        ; width and ht are 440
(define ignore (read-bytes (- 387976 size) in))    ; size of file is 387976
(define imgdata (read-bytes size in))
(close-input-port in)

(define img (make-object bitmap% imgdata 440 440))
img

However, it only shows a random mix of black and white pixels:

Using 440*440 instead of 2*440*440 also does not work.

Following code also does not read the image:

(define img (make-object bitmap% in 'unknown))

This does not display any image at all.

Where is the problem and how can I solve it?

Answer 1

You are calculating the image data offset correctly and this data appears to be raw and uncompressed. The problem appears to be just that Racket does not support this kind of image data. This is 16-bit, single channel, intensity data. Each two bytes represents a pixel in grayscale (and actually for this image only 10 bits is used, the other 6 bits should be ignored).

The make-object make-bitmap functions appear to support only color (24 or 32) or monochrome (1) bit depth. This is what jumps out at you in your example: nowhere when you create the bitmap do you state the pixels are 16-bit. Also absent is anything about byte order. And nowhere in the Racket documentation does it appear that it allows you to specify either of these.

Lack of support for 16-bit grayscale data seems evident in the get-depth function documentation:

(send a-bitmap get-depth) → exact-nonnegative-integer? Gets the color depth of the bitmap, which is 1 for a monochrome bitmap and 32 for a color bitmap. See also is-color?.

Here's one solution, looping through and converting each 16-bit pixel into an ARGB pixel.

#lang racket/gui
(define in (open-input-file "CR-MONO1-10-chest" #:mode 'binary))
(define size (* 2 440 440))                        ; width and ht are 440
(define ignore (read-bytes (- 387976 size) in))    ; size of file is 387976
(define imgdata (read-bytes size in))
(close-input-port in)

(define rgbdata (make-bytes (* 4 440 440)))
(define img (make-object bitmap% 440 440))

(define max 1024.0) ; 10 bits valid

(for ([y 440])
  (for ([x 440])
      (define index (+ (* y 440) x))
      (define b1 (bytes-ref imgdata (+ (* index 2) 0))) ; first byte
      (define b2 (arithmetic-shift (bytes-ref imgdata (+ (* index 2) 1)) 8)) ; second byte
      (define val (bitwise-xor b1 b2)) ; combine bytes
      (define screenval (exact-floor (* 255 (/ val max)))) ; convert to 8-bit screen value

      ; create ARGB pixel
      (define b (bytes (bytes-ref (make-bytes 1 255) 0) (bytes-ref (make-bytes 1 screenval) 0) (bytes-ref (make-bytes 1 screenval) 0) (bytes-ref (make-bytes 1 screenval) 0)))    
      (send img set-argb-pixels x y 1 1 b)
  ))

img

Producing:

A note of caution though: be aware that there are many, many ways image data can be stored in DICOM files (some good examples on the page you linked). And there are many parts of the DICOM header you would need to pay attention to decode the image data correctly.

Answer 2

The problem is that the Racket drawing library does not recognize that image encoding. I suggest converting it to 32-bit ARGB:

(require racket/draw)

(define (dicom->bitmap path x y)
  (let* ([bmp            (make-object bitmap% x y)]
         [dc             (send bmp make-dc)]
         [dicom          (file->bytes path #:mode 'binary)]
         [header-size    (- (file-size path) (* 2 x y))]
         [dicom-img/raw  (subbytes dicom header-size)]
         [dicom-img/argb (dicom-img->argb dicom-img/raw)])
    (send dc set-argb-pixels 0 0 x y dicom-img/argb)
    (send dc get-bitmap)))

(define (dicom-img->argb bytes)
  (let* ([len         (bytes-length bytes)]
         [pixel-count (/ len 2)]
         [argb        (make-bytes (* 4 pixel-count))])
    (define (set-pixel! value ix)
      (let ([offset (* 4 ix)])
        (bytes-set! argb offset 0)
        (bytes-set! argb (+ 1 offset) value)
        (bytes-set! argb (+ 2 offset) value)
        (bytes-set! argb (+ 3 offset) value)))
    (for ([ix (in-range pixel-count)])
      (let* ([offset       (* 2 ix)]
             [pixel-value  (+ (bytes-ref bytes offset)
                              (arithmetic-shift (bytes-ref bytes (+ 1 offset)) 8))]
             [scaled-value (arithmetic-shift pixel-value -2)])
        (set-pixel! scaled-value ix)))
    argb))

Then you can call it like so:

(dicom->bitmap "CR-MONO1-10-chest" 440 440)

This particular program works only for 10 bits/pixel stored in 2 bytes each in little-endian order, but with modest effort you could parameterize it for other encodings.

If you have a file with multiple images in that format, all at the end, this program should be able to extract them as a list of bitmaps.

(require racket/draw)

(define (dicom->bitmap* path x y z)
  (let* ([dicom           (file->bytes path #:mode 'binary)]
         [img-size        (* 2 x y z)]
         [header-size     (- (file-size path) img-size)]
         [dicom-img/raw*  (for/list ([z^ (in-range z)])
                            (let* ([offset    (+ header-size (* z^ img-size))]
                                   [img-bytes (subbytes dicom offset (+ offset img-size))])
                              img-bytes))])
    (map (λ (raw) (raw->bitmap raw x y)) dicom-img/raw*)))

(define (raw->bitmap bytes x y)
  (let* ([bmp             (make-object bitmap% x y)]
         [drawing-context (send bmp make-dc)]
         [dicom-img/argb  (raw->argb bytes)])
    (send drawing-context set-argb-pixels 0 0 x y dicom-img/argb)
    (send drawing-context get-bitmap)))

(define (raw->argb bytes)
  (let* ([len         (bytes-length bytes)]
         [pixel-count (/ len 2)]
         [argb        (make-bytes (* 4 pixel-count))])
    (define (set-pixel! value ix)
      (let ([offset (* 4 ix)])
        (bytes-set! argb offset 0)
        (bytes-set! argb (+ 1 offset) value)
        (bytes-set! argb (+ 2 offset) value)
        (bytes-set! argb (+ 3 offset) value)))
    (for ([ix (in-range pixel-count)])
      (let* ([offset       (* 2 ix)]
             [pixel-value  (+ (bytes-ref bytes offset)
                              (arithmetic-shift (bytes-ref bytes (+ 1 offset)) 8))]
             [scaled-value (arithmetic-shift pixel-value -2)])
        (set-pixel! scaled-value ix)))
    argb))

Where z is the number of images. I've only been able to test it with z=1 though: (dicom->bitmap* "CR-MONO1-10-chest" 440 440 1)