Skip to content

Advanced Tutorial

Available modules

Making SR (Structured Report)

Utility Functions

Making a SR

from pygrpm.dicom.sr import make_sr

sr_dicom = make_sr('<text_to_add>', '<optional_path_to_reference_dicom_file>', '<optional_purpose_of_reference_code>')

Storing a text into a DICOM SR

from pygrpm.dicom.sr import make_sr_from_text

# Make SR DICOM dataset from text
sr_with_text = make_sr_from_text("My Text", ['./CT1.dcm', './CT2.dcm', './CT3.dcm', ...])

# Afterward, can be written to a file like this:
sr_with_text.save_as('./my_sr.dcm')

SR Builder

Generate DICOM SR File

from pygrpm.dicom.sr import SRBuilder

sr_builder = SRBuilder('<sr_target_file_path>', '<purpose_of_reference_code>')

sr_builder.build()

sr_builder.dicom_sr #How to access the Dicom SR object
  1. A DicomSRBuilder object is created with parameters:
  2. sr_target_file_path: Path of the reference DICOM file to use to populate the DICOM SR. Dicom files are also accepted directly

  3. purpose_of_reference_code: Dictionary with information tags on WHY the target file is used to generate the DICOM SR For example: python { "CodeValue": "SR", "CodeMeaning": "Structured Report Document", "CodingSchemeDesignator": "DCM", }

    For more information on what this dictionary might contain: https://dicom.nema.org/medical/dicom/2018b/output/chtml/part03/sect_8.8.html#table_8.8-1a and https://dicom.nema.org/medical/dicom/2018b/output/chtml/part03/chapter_8.html

    The most common CodingSchemeDesignator is "DCM", here is a reference to the values accepted for this coding scheme: https://dicom.nema.org/medical/dicom/current/output/chtml/part16/chapter_D.html

  4. Calling the build() method then builds the DICOM SR, assembling MetaData with available data

Adding a content sequence

sr_dicom.add_content_sequence(<content_sequence>)
  • Adds the content sequence given as parameter to the DICOMSR object. A content sequence is a structured content (a list generally) giving relevant information on the associated document. Includes, minimally:
    • The Value Type (TEXT, NUM, CODE, DATETIME, etc.)
    • The Value (must match in type with the Value Type)
    • The ConceptNameCodeSequence in the same format as the purpose_of_reference_code used above

For example: python { "ValueType": "TEXT", "Value": "This is some text", "ConceptNameCodeSequence": { "CodeValue": "1233", "CodeMeaning": "Test", "CodingSchemeDesignator": "DCM", }, "SomeOtherTag" : "Some additionnal information" }

  • The content sequence is validated to appropriate standards before adding it (see https://dicom.nema.org/medical/dicom/2020b/output/chtml/part03/sect_C.17.3.2.html for more information)

Containers

Container objects are common content sequences that must respect additionnal format restrictions. Here is an example of a valid container object for the method add_content_sequence:

{
  "ValueType": "CONTAINER",
  "ConceptNameCodeSequence": {
      "CodeValue": "DOC",
      "CodeMeaning": "Document",
      "CodingSchemeDesignator": "DCM",
  },
  "ContinuityOfContent": "SEPERATE",
  "Value": [
      {
          "RelationshipType": "HAS PROPERTIES",
          "ValueType": "TEXT",
          "ConceptNameCodeSequence": {
              "CodeValue": "113012",
              "CodeMeaning": "Key Object Description",
              "CodingSchemeDesignator": "DCM",
          },
          "Value": "This is a text",
      }
  ],
}
  • Refer to https://dicom.nema.org/medical/dicom/current/output/chtml/part03/sect_C.17.3.2.4.html for more information

Additionnal Information

Users who wish to have more information on the creation of SR are invited to read documentation concerning SR dose content sequences and DICOM part 03 section 17.

Acknowledgements

  • This work is based on a previous project from Samuel Ouellet (https://gitlab.chudequebec.ca/sam23/dicom-sr-builder)
  • DICOM SR files assembled by this module are in accordance with PS3.6 2023b standards. PS3.6 standard can (and most certainly will) change over time. See https://dicom.nema.org/medical/dicom/current/output/html/part06.html for current standards.
  • Thanks to Gabriel Couture for some much needed revision work

Contact

Pierre-Luc Asselin, B. Ing., plass2@ulaval.ca

TG43

Python package to calculate the TG43 formalism based on xarray and the data from ESTRO.

Available seeds

  • SelectSeed
  • MicroSelectronV2
  • MBDCA
  • Flexisource

Usage

import matplotlib.pyplot as plt
from pygrpm.tg43 import Seed

seed = Seed("MBDCA")
print(seed.gL)
print(seed.F)

# Display dosimetric grid
grid = seed.grid()
# Plot arbitrary 100th slice
plt.imshow(grid.data[:, :, 100])
plt.show()

# Plot mean dose vs radius with std
mean = seed.get_mean(grid)
std = seed.get_std(grid)
plt.errorbar(mean["r"].data, mean.data, std.data, fmt='o')
plt.show()

Index tracker

Python package to allow scrolling through slices of 3-D images using the matplotlib backend.

Usage

Here is shortly described the main usage of this module. See the docstring of each function, class, method and attribute for more details.

In the first example shown below, a 3-D array is populated with random values, and the object tracker is initialized with a few arguments. The first one represents the unpacked tuple of the objects (fig,ax), followed by the 3-D array that represents the image, and other parameters to customize how the image is shown by the matplotlib backend. It is important to provide slice_axis, in order to define where the slices are located in the image. The object contains as members ax, im, and cb (colorbar), which allows further customization:

# 3-D array with shape (512,512,10)
random_array = np.random.rand(512,512,10)
tracker = IndexTracker(*plt.subplots() ,     # plt objects
                        random_array , # 3-D array of image
                        vmin = 0 ,   # Minimum value (colormap)
                        vmax = 0.8 , # Maximum value (colormap)
                        slice_axis = 2 , # Axis containing slices
                        cmap='gray') # Colormap
tracker.ax.set_title('My 3-D random image')
tracker.show()

In the following example, it is shown how easily one can format the image style, either before or after the object tracker is initialized: 

# Formatting the plot
random_array = np.random.rand(10,512,512)
fig, ax = plt.subplots()
fig.set_tight_layout(True)
ax.set_title('My first 3-D random image')
tracker = IndexTracker(fig, ax, random_array, slice_axis=0)
# Change min and maximum value of colormap
minv = 0.0
maxv = 0.03
tracker.im.set_clim(minv, maxv)
# Change colorbar configuration
tracker.colorbar(extend='both')
tracker.show()

# Passing a color normalization
from matplotlib.colors import LogNorm

random_array = np.random.rand(10,512,512)
tracker = IndexTracker(*plt.subplots(),
                       random_array,
                       slice_axis=0,
                       cmap='viridis',
                       norm=LogNorm(vmin=0.01, vmax=1))
tracker.show()

Jupyter environment

In order to use IndexTracker inside the jupyter environment (e.g. Jupyter lab, and Jupyter Notebook*), it is necessary to add the following statement prior to your code:

import matplotlib
matplotlib.use('QtAgg')
tracker = IndexTracker(...)

This statement changes the matplotlib's backend, allowing to show images in a new window.

See this for more information

NISTParser

A simple class to extract information from certain NIST webpages. At the time of writing this covers atomic and electronic cross-sections, material attenuation, as well as material composition.

get_cross_sections

This method retrieves the desired cross-sections of an element at given energies on the NIST website in (barns/electron), barn=10^-24cm^2.

Simple use example:

import numpy as np
from pygrpm.material.nistparser import get_cross_sections

# Define the energies in KeV
# Numpy array is not mandatory, can be any sequence
energies = np.linspace(30, 200, 200)

# Prints the returned list
print(get_cross_sections("H", energies))

get_electronic_cross_section

See get_cross_section(), method use is identical minus the options argument

get_attenuations

Method to retrieve the attenuation of a material in cm^2/g at given energies on the NIST website. * Note the option parameter which can specify the attenuation physics types * Note the outopt parameter which can alter the returned information

Example is similar to get_cross_sections()

get_composition

This method is used to get and parse material composition from https://physics.nist.gov/cgi-bin/Star/compos.pl

Simple use example:

from pygrpm.material.nistparser import get_composition, NISTMaterials

# Prints the returned dictionary
print(get_composition(NISTMaterials.M3_WAX))

Note that get_composition expects the material to be of instance NISTMaterials

Acknowledgements

  • This submodule makes use of the HTMLTablePasrer built by Josua Schmid, further information can be found in the pygrpm.nistparser.nistparser.py file header.
  • This submodule is also dependent on the data provided by https://www.nist.gov/pml

Hounsfield conversion

A helper class meant to offer a quick and simple means to convert an HU array to density values based on a provided conversion table. Note that the conversion factors for HU to density are generally provided by the CT manufacturer. This class is currently only able to be read under csv format type.

General notes

The conversion factors for HU to density are generally provided by the CT manufacturer. This file is currently only able to be read under csv format type.

Usage

Using the following sample data as ./curve.csv file,

HU,Mass density [g/cm^3]
-1000,0.00121
-804,0.217
-505,0.508
-72,0.967
-32,0.99
7,1.018
44,1.061
52,1.071
254,1.159
4000,3.21

A call through the class can be made to rescale an arbitrary array of Hounsfield unit to density values.

import numpy as np
from pygrpm.material.hounsfield_conversion import ConvertHUToDensity

fithu = ConvertHUToDensity()
my_curve = fithu.load_curve_csv("./curve.csv")

# Note that setting plot to True generates a matplotlib plot of the curve fitting
data, fit_params, labels = fithu.fit_curve(my_curve, plot=True)
# Fit returns unused in this example

my_real_image = np.array([-980, -1000., -823., 1, 20, 700, 2900])
densities = fithu.apply_fit(my_real_image, my_curve)

print(densities)  # [0.02702014 0.00652871 0.18787786 1.0291972 1.03955733 1.41034076, 2.60993423]

Note that the class will automatically "wrap" any HU values beyond its set extrema to their nearest allowed extrema value. These values can be altered using the .set_extrema(min_val, max_val) method.

TG263

Validating the structure name

from pygrpm.tg263.nomenclature import is_structure_name_allowed

result = is_structure_name_allowed('Prostate')
# Result is `True`

This will simply specify if the parameter string is an allowed structure Primary Name in the TG263 structure database.

Finding a structure

from pygrpm.tg263.nomenclature import find_structures

result = find_structures('SpinalCord')

This will print all information on every structure that as a TG263 primary name matching the parameter string. Structure information is return in dictionary format, and all structures are returned into a list.

Printing structure info

from pygrpm.tg263.nomenclature import print_structure_info

print_structure_info('<structure_to_print>')

This will print relevant structure information on the structure given as a parameter.

Validating a structure's format

from pygrpm.tg263.nomenclature import is_structure_valid

validation, log = is_structure_valid('<structure_to_validate>')

This will return a boolean (validation) wheiter the structure is valid. If not, the log will explain more precisely the nature of the problem.

Criterias are as follows:

  • Structures are python libraries
  • Structures must contain keys* for:
    • tg263PrimaryName
    • tg263ReverseOrderName
    • targetType
    • majorCategory
    • minorCategory
    • anatomicGroup
    • description
    • fmaid
    • nCharacters
  • Every key as a string value, except for fmaid and nCharacters, which must be integers.

  • String values cannot contain any of these characters: "<>:"'/|?*()#$"

  • The containt of requisite keys is detailled in the RPT_263 report (link below).

DicomReader

A helper class meant to ingest a given folder of Dicom files and sort them into studies and series. This is aided with the DicomSeries and DicomStudy dataclasses.

Basic Usage

The class constructor is to be fed a system path for a folder containing, under any folder structure, the desired Dicom files.

from pygrpm.dicom import DicomReader

my_path = "/path/to/dicoms/folder/"
reader = DicomReader(my_path)

reader.studies # list of DicomStudy classes
reader.studies[0].series # list of DicomSeries classes

Dicoms can also be fed through unsorted through a generator with the following

from pygrpm.dicom import DicomReader

my_path = "/path/to/dicoms/folder/"
reader = DicomReader(my_path)

for study in reader.yield_dicoms():
    do_something(study)

DicomReader

There are a few helper methods available. Notably <reader>.get_study(<study_UID>) and <reader>.get_series(<study_UID>, <series_UID>) for easy access to given series or studies.

DicomStudy

It is also possible to filter a study by present modalities, <study>.filter_by_modality("CT") which would return a list of DicomSeries dataclasses of modality "CT".

DicomSeries

The DicomSeries dataclass offers easy access to the first instance of the series through <series>.first_instance as sorted by Z-slice position. Custom sorting is also available through <series>.sort_instances_by_tag(<dicom_tag>).

<series>.numpy returns a constructed numpy array from multiple 2D PixelData components as in the case of CT images, or fast-forwards pydicom's pixel_array in the case of 3D arrays such as in RTDose or RTStruct.

The pydicom FileDatasets are available within their respective DicomSeries dataclass within the <series>.instances list.

Egsphant Generator

A class designed to transform data parsed from a CT volume in numpy format, and provide an egsphant type text file.

Usage

The following example will assume that the user has already determined volume center, pixel spacing, and stacked the CT volume in a 3D numpy array. Note that the volume is expected to be in HU values for proper material assignment.

import numpy as np
from pygrpm.geometry.egsphant import make_egsphant_from_numpy

spacing = np.array(...)
center = np.array(...)
volume = np.array(...)  # Assuming slices in last position

# Select built-in materials list CT01 and specify slices in last axis
myegs = make_egsphant_from_numpy(volume, spacings=spacing, center=center,
                                 materials="CT01", slice_axis=-1)
# Can see the full string contents directly if desired
contents = myegs.content

myegs.write_to_file("./path/to/myvolume.egsphant")

Individual sections of the egsphant can also be accessed via the header, voxel_position_string, material_string, and density_string properties of the provided Egsphant class.

Acknowledgements

This work is a basic implementation of the TG263 (https://www.aapm.org/pubs/reports/RPT_263.pdf)

The allowed structure names (and corresponding information) were taken from https://www.aapm.org/pubs/reports/RPT_263_Supplemental/ .

The initial allowed structure names were taken from the ESAPIX project (https://github.com/rexcardan/ESAPIX), made by Rex Cardan. The ESAPIX license is included in the LICENSE file.