Skip to content

handling DICOM-RT structure sets in oro.dicom

June 6, 2013

In my previous post I demonstrated how to read image data and meta-data from DICOM files, using the R package oro.dicom. In radiotherapy there is another kind of data that I work with, known as a DICOM-RT structure set. This is a special type of DICOM file that contains geometric contours of regions of interest, rather than image data.

library(oro.dicom)
dcmImages <- readDICOM("/dev/ITK4SampleDCMRT",
   recursive=FALSE, exclude="sql")
dcm.info <- dicomTable(dcmImages$hdr)
rtstruct.idx <-
   which(dcm.info["0008-0060-Modality"] == "RTSTRUCT")
rtstruct.files <- rownames(dcm.info[rtstruct.idx,])
print(rtstruct.files)

## [1] "/dev/ITK4SampleDCMRT/RS.1.2.246.352.71.4.886768594.5257.20090622110825.dcm"

As Dan Mueller has pointed out [1], RTSTRUCT files have been superceded by the Surface Segmentation Storage SOP Class (Supplement 132), which was finalised in 2008. However, RTSTRUCT is still widely used in clinical practice. Whether you prefer isometric contour lines or a polygonal 3D mesh largely depends on what application you are using the data for. For example, if you are developing a meshless algorithm for image segmentation, then the RTSTRUCT format is a better starting point. You can always convert to a 3D mesh using marching cubes [2], as in the hybrid algorithm of Chen & Metaxas [3].

The algorithm for reading a radiotherapy structure set from an RTSTRUCT file has been outlined by Jason Dowling et al. [4,5]. All of the contour information is stored under group 3006, amongst all of the other DICOM meta-data.

structureSetROISequence <- function(header, idx) {
  roiNum <- vector(mode="integer")
  roiName <- vector()
  sIdx <- 0
  while(header[idx,1] != "3006" || header[idx,2] != "0039")
  {
    if (header[idx,1] == "3006" && header[idx,2] == "0022") {
      sIdx <- sIdx + 1
      roiNum[sIdx] <- as.integer(header[idx,6])
    }
    else if (header[idx,1] == "3006" && header[idx,2] == "0026") {
      roiName[sIdx] <- header[idx,6]
    }
    idx <- idx+1
  }
  struct <- cbind(roiNum,roiName)
  colnames(struct) <- c("ROINumber","ROIName")
  roiDecl <- data.frame(struct, row.names=1,
                        stringsAsFactors=FALSE)  
  list(names=roiDecl,end=idx)
}

file <- rtstruct.files[1]
rtstruct <- readDICOMFile(file, skipSequence=FALSE,
                          pixelData=FALSE)
idx <- 1
while(rtstruct$hdr[idx,1] != "3006" || rtstruct$hdr[idx,2] != "0020")
   idx <- idx+1;
ssroi <- structureSetROISequence(rtstruct$hdr, idx)
print(ssroi$names)

## ROIName
## 1 GS1
## 6 MRI_RECTUM_RT
## 5 MRI_BONEOUTER_RT
## 2 GS2
## 3 GS3
## 4 MRI_BLADDER_RT
## 7 MRI_PROSTATE_RO

The RTSTRUCT file contains a sequence of structures (DICOM value representation SQ), which in turn contain a sequence of zero or more ROIs (one ROI for each image slice spanned by the structure), which themselves are defined in terms of a string of X,Y,Z coordinates (Decimal String, value representation DS). Further complicating all of this multi-level recursion is that all of these components can be stored out of order. The ROI number (Tag <3006,022>) needs to be matched up with the referenced ROI number (Tag <3006,0084>), as the following code illustrates:

roiContourSequence <- function(roiDecl, header, idx) {
  seqLen <- as.integer(header[idx,5])
  bytes <- 0
  items <- list()
  # until the end of the sequence of sequences...
  while(bytes < seqLen) {
    idx <- idx+1
    # new item
    if (header[idx,1] == "FFFE" && header[idx,2] ==  "E000" && header[idx,7] == "(3006,0039)") {
      item <- processItem(roiDecl, header, idx)
      bytes <- bytes + as.integer(header[idx,5]) + 8
      items <- list(unlist(items,recursive=FALSE),item)
      idx <- item$end # skip ahead
    }
  }
  list(items=unlist(items,recursive=FALSE),end=idx)
}

processItem <- function(roiDecl, header, idx) {
  len <- as.integer(header[idx,5])
  bytes <- 0
  color <- NA
  refID <- NA
  pts <- vector()
  while(bytes < len) {
    idx <- idx + 1
    bytes <- bytes + as.integer(header[idx,5]) + 8
    # ROIDisplayColor
    if (header[idx,1] == "3006" && header[idx,2] ==  "002A") {
      color <- header[idx,6]
    }
    # ContourSequence
    else if (header[idx,1] == "3006" && header[idx,2] ==  "0040") {
      pts <- contourSQ(header, idx)
      idx <- pts$end
    }
    # ReferencedROINumber
    else if (header[idx,1] == "3006" && header[idx,2] ==  "0084") {
      refID <- header[idx,6]
    }
  }
  list(col=color,ref=refID,pts=pts,end=idx)
}

contourSQ <- function(header, idx) {
  seqLen <- as.integer(header[idx,5])
  seqBytes <- 0
  pts <- vector()
  numPts <- vector()
  types <- vector()
  # until the end of the sequence of sequences...
  while(seqBytes < seqLen) {
    idx <- idx+1
    # new item
    if (header[idx,1] == "FFFE" && header[idx,2] ==  "E000") {
      len <- as.integer(header[idx,5])
      seqBytes <- seqBytes + len + 8
      bytes <- 0
      while (bytes < len) {
        idx <- idx+1
        # ContourImageSequence
        if (header[idx,1] == "3006" && header[idx,2] ==  "0016") { }
        # ContourGeometricType
        else if (header[idx,1] == "3006" && header[idx,2] ==  "0042") {
          types <- c(types, header[idx,6])
        }
        # NumberOfContourPoints
        else if (header[idx,1] == "3006" && header[idx,2] ==  "0046") {
          numPts <- c(numPts, header[idx,6])
        }
        # ContourData
        else if (header[idx,1] == "3006" && header[idx,2] ==  "0050") {
          pts <- c(pts, header[idx,6])
        }
        else next; # ignore this line
        bytes <- bytes + as.integer(header[idx,5]) + 8
      }
    }
  }
  list(num=numPts,pts=pts,types=types,end=idx)
}

contours <- roiContourSequence(ssroi$names, rtstruct$hdr, ssroi$end)
print(length(contours$items))

## [1] 28

print(ssroi$names[(contours$items[[26]]),])

## [1] "MRI_PROSTATE_RO"

print(contours$items[c(25,28)])

## $col
## [1] "255 255 128"
##
## $end
## [1] 2282

print(contours$items[[27]]$num)

## [1] "52" "66" "84" "92" "98" "104" "108" "106" "98" "82" "64"
## [12] "54"

print(contours$items[[27]]$types)

## [1] "CLOSED_PLANAR" "CLOSED_PLANAR" "CLOSED_PLANAR" "CLOSED_PLANAR"
## [5] "CLOSED_PLANAR" "CLOSED_PLANAR" "CLOSED_PLANAR" "CLOSED_PLANAR"
## [9] "CLOSED_PLANAR" "CLOSED_PLANAR" "CLOSED_PLANAR" "CLOSED_PLANAR"

print(contours$items[[27]]$pts[1])

## [1] "-4.25 -38.72 -17.2 -2.77 -38.57 -17.2 -2.72 -38.53 -17.2 -1.28 -37.79 -17.2 2.0e-1 -37.15 -17.2 3.2e-1 -37.05 -17.2 1.68 -36.34 -17.2 2.81 -35.56 -17.2 3.17 -35.17 -17.2 4.06 -34.08 -17.2 4.65 -32.86 -17.2 4.81 -32.59 -17.2 5 -31.11 -17.2 4.95 -29.62 -17.2 4.84 -28.14 -17.2 4.65 -27.5 -17.2 4.41 -26.65 -17.2 3.94 -25.17 -17.2 3.17 -24.27 -17.2 2.39 -23.69 -17.2 1.68 -23.19 -17.2 2.0e-1 -22.49 -17.2 -9.9e-1 -22.2 -17.2 -1.28 -22.12 -17.2 -2.77 -21.68 -17.2 -4.25 -21.49 -17.2 -5.74 -21.38 -17.2 -7.22 -21.27 -17.2 -8.71 -21.72 -17.2 -10.19 -22.1 -17.2 -10.32 -22.2 -17.2 -11.67 -22.94 -17.2 -12.82 -23.69 -17.2 -13.16 -23.97 -17.2 -14.31 -25.17 -17.2 -14.64 -25.69 -17.2 -15.17 -26.65 -17.2 -15.31 -28.14 -17.2 -15.28 -29.62 -17.2 -14.87 -31.11 -17.2 -14.64 -31.54 -17.2 -14.16 -32.59 -17.2 -13.52 -34.08 -17.2 -13.16 -34.53 -17.2 -12.17 -35.56 -17.2 -11.67 -35.96 -17.2 -10.19 -36.52 -17.2 -8.76 -37.05 -17.2 -8.71 -37.07 -17.2 -7.22 -37.62 -17.2 -5.74 -38.15 -17.2 -4.7 -38.53 -17.2"

The next step will be to translate these coordinates from physical to image space and use R plotting commands to create a mask for each contour.

References

[1]: S Gorthi, M Bach Cuadra, J Thiran (2009) Exporting Contours to DICOM-RT Structure Sets The Insight Journal

[2]: W E Lorensen, H E Cline (1987) Marching Cubes: A high resolution 3D surface construction algorithm In Proc. 14th SIGGRAPH

[3]: T Chen, D Metaxas (2005) A hybrid framework for 3D medical image segmentation Medical Image Analysis 9(6)

[4]: J Dowling, M Malaterre, P B Greer, O Salvado (2009) Importing Contours from DICOM-RT Structure Sets The Insight Journal

[5]: J Dowling (2013) Importing Contours from DICOM-RT Structure Sets with ITK4 The Insight Journal

Advertisements

From → Imaging, R

7 Comments
  1. Kang xiaochun permalink

    need help.

    RStudio Console:
    Error in if (header[idx, 1] == “FFFE” && header[idx, 2] == “E000” && header[idx, :
    missing value where TRUE/FALSE needed

  2. Hi Kang, unfortunately the above code no longer works with the current version of oro.dicom

    I would recommend trying the RadOnc R package instead, which is available on CRAN:

    https://cran.r-project.org/web/packages/RadOnc/index.html

    • Kang xiaochun permalink

      Could you please send me a copy of the older version oro.dixcom,thanks

  3. Kang xiaochun permalink

    I installed the ” oro.dicom_0.3.5.tar.gz 06-May-2012″,
    when test the code, I got the same error,what’s wrong?

  4. Have you tried out this code with the example DICOM-RT file from Jason Dowling’s paper?

    http://www.insight-journal.org/browse/publication/887

    The following is the setup that I used (apologies, very outdated):

    R version 2.15.3 (2013-03-01)
    Platform: x86_64-w64-mingw32/x64 (64-bit)

    oro.dicom_0.3.7
    oro.nifti_0.3.9
    bitops_1.0-6

    It is possible that the DICOM-RT file that you are trying to read is in a slightly different format. The code that I wrote was for treatment plans that were created in Pinnacle and exported from Elekta MOSAIQ. Are you able to view the contours in other open source applications, such as PyDICOM or RadOnc?

Trackbacks & Pingbacks

  1. Interactive image stacks in R | Matt Moores
  2. Reading DICOM-RT in R using RadOnc | Matt Moores

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Richard Everitt's blog

Computational Bayesian statistics

Let's Look at the Figures

David Firth's blog

Nicholas Tierney

Computational Bayesian statistics

One weiRd tip

Computational Bayesian statistics

Series B'log

discussion blog for JRSS Series B papers

Mad (Data) Scientist

Musings, useful code etc. on R and data science

R-bloggers

R news and tutorials contributed by (750) R bloggers

Another Astrostatistics Blog

The random musings of a reformed astronomer ...

Darren Wilkinson's research blog

Statistics, computing, data science, Bayes, stochastic modelling, systems biology and bioinformatics

CHANCE

Computational Bayesian statistics

StatsLife - Significance magazine

Computational Bayesian statistics

(badness 10000)

Computational Bayesian statistics

Igor Kromin

Computational Bayesian statistics

Statisfaction

I can't get no

Xi'an's Og

an attempt at bloggin, nothing more...

Sam Clifford

Postdoctoral Fellow, Bayesian Statistics, Aerosol Science

Bayesian Research & Applications Group

Frontier Research in Bayesian Methodology & Computation

%d bloggers like this: