From ‘OsiriX’ ROIs to PET voxels and radiomics features

Overview

Regions of interest (ROIs) drawn slice by slice in the ‘OsiriX’ DICOM viewer can be exported as an XML file. For each ROI on each slice, the file stores the polygon corner points, the value of every voxel inside the ROI and summary statistics, but not the voxel positions. osirad finds the voxel positions:

1. read_osirix_rois() reads the XML file.
2. suv_factor_from_header() and read_pet_slice() read the PET DICOM slices and convert them to standardised uptake values (SUV).
3. match_roi_voxels() finds, near the ROI corner points, the pixels whose SUV equals the voxel values in the XML file.
4. extract_roi_voxels() runs steps 1 to 3 for all ROIs and returns one table with the position (pixels and mm) and the SUV of every ROI voxel.
5. pyradiomics_features() sends the PET series and the ROI voxels to the ‘Python’ library ‘PyRadiomics’ to compute radiomics features.

Export the ROIs from ‘OsiriX’

1. Draw the ROIs slice by slice on the original PET series (closed polygon or another ROI tool).
2. Export the ROIs: Plugins > ROI Tools > Export ROIs, XML format.
3. Export the PET series as DICOM files (File > Export > Current image only, not the fused PET/CT). Keep the files of this series, and only these, in one folder.

Example data

The package contains a small synthetic PET phantom (32 x 32 pixels, 8 slices) with two ROIs, “Tumour” and “Node”, exported in the ‘OsiriX’ XML format.

library(osirad)
xml_file = system.file("extdata", "phantom", "ROI.xml", package = "osirad")
pet_dir = system.file("extdata", "phantom", "pet", package = "osirad")

ROI voxels

x = extract_roi_voxels(xml_file, pet_dir)
x
#> 'OsiriX' ROI voxels: 8 ROI(s) on 6 slice(s), 222 voxels
#> ROI names: Tumour, Node 
#> SUV factor: 0.000276308 
#> XML values: 222 | clear match: 221 | listed twice: 0 | tie by neighbours: 1 | tie by distance: 0 | not found: 0
head(x$voxels)
#>   roi_id roi_name image_index x_px y_px x_mm y_mm  z_mm      suv
#> 1      1   Tumour           1   15   12   -2  -14 56.73 2.343108
#> 2      1   Tumour           1   16   12    2  -14 56.73 2.226960
#> 3      1   Tumour           1   14   13   -6  -10 56.73 1.887495
#> 4      1   Tumour           1   15   13   -2  -10 56.73 2.404036
#> 5      1   Tumour           1   16   13    2  -10 56.73 2.685511
#> 6      1   Tumour           1   17   13    6  -10 56.73 2.356956

The printed report shows how each XML value was matched to a pixel:

• clear match: the value is found at exactly as many pixels as ‘OsiriX’ lists it.
• listed twice: ‘OsiriX’ lists the same pixel more than once.
• tie by neighbours / by distance: the same SUV is found at more pixels than listed; the pixels touching the most other ROI voxels (then the ones nearest the ROI centre) are kept.
• not found: no pixel has this value. This happens when pet_dir is not the series the ROIs were drawn on; extract_roi_voxels() then gives a warning.

The match per ROI is in x$match:

x$match
#>   roi_id roi_name image_index xml_values clear listed_twice tie_by_neighbours
#> 1      1   Tumour           1         14    14            0                 0
#> 2      2   Tumour           2         33    32            0                 1
#> 3      3   Tumour           3         54    54            0                 0
#> 4      4   Tumour           4         54    54            0                 0
#> 5      5     Node           4         10    10            0                 0
#> 6      6   Tumour           5         33    33            0                 0
#> 7      7     Node           5         10    10            0                 0
#> 8      8   Tumour           6         14    14            0                 0
#>   tie_by_distance not_found voxels
#> 1               0         0     14
#> 2               0         0     33
#> 3               0         0     54
#> 4               0         0     54
#> 5               0         0     10
#> 6               0         0     33
#> 7               0         0     10
#> 8               0         0     14

Check the voxels on the images:

plot_roi_slices(x, image_index = 2:5)

Save the voxel table for other software:

write.csv(x$voxels, "roi_voxels.csv", row.names = FALSE)

Radiomics features with ‘PyRadiomics’

pyradiomics_features() needs ‘Python’ with the modules pyradiomics, SimpleITK and numpy, and the R package ‘reticulate’. One way to set up ‘Python’ with ‘conda’ (in a terminal):

conda create -n radiomics python=3.9 "numpy<2"
conda activate radiomics
pip install pyradiomics

On computers where pip install pyradiomics fails (for example Apple silicon), install it from the source code:

pip install versioneer "setuptools<70"
pip install --no-build-isolation git+https://github.com/AIM-Harvard/pyradiomics.git@v3.1.0

Then give the path of this ‘Python’ to pyradiomics_features() (conda activate radiomics; which python shows it). The default settings are a fixed bin width of 0.25 SUV, no intensity normalisation and no resampling; all ROIs form one mask.

rad = pyradiomics_features(x, pet_dir, python = Sys.getenv("OSIRAD_PYTHON"))
head(rad$features)
rad$check

rad$check compares the ‘Python’ mask with the R voxels: the number of voxels, the largest distance between a voxel and its slice, and the largest SUV difference between ‘Python’ and R.