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

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)