ImagingDataCommons / slim

Interoperable web-based slide microscopy viewer and annotation tool
https://imagingdatacommons.github.io/slim/
Apache License 2.0
118 stars 36 forks source link

Parametric maps of exemplary datasets #154

Closed maxfscher closed 1 year ago

maxfscher commented 1 year ago

Hey, first of all thank you for this very nice viewer and the very good documentation. I have a question regarding the example datasets that you provide here. Can you share some information how the parametric maps from the Study ID C3N-01016 were generated? Is the process somewhere described or can you provide further hints on which tool did you use? I am especially interested in how you defined the real_world_value_mappings in the parametric map dicom file.

Thanks in advance Max

fedorov commented 1 year ago

We discussed this today with @cgorman - he is the one who generated the samples you mentioned using highdicom. Chris can share the part of the code that generated those.

cgorman commented 1 year ago

@maxfscher Hi Max,

Here is the code I used to convert the pixel array into a parametric map (in Python). The conversion to unsigned int is not strictly necessary since PMs support floating point pixel data, but we found that converting to int makes things faster and less likely to cause problems with other tools.

This snippet requires the packages highdicom, opencv-python, pydicom, and numpy.

from random import randint
from socket import gethostname

import cv2
import highdicom as hd
import numpy as np
import pydicom
from pydicom.sr.codedict import codes

def parametric_map(
    source_image: pydicom.Dataset,
    input_map: np.ndarray,
    class_label: str,
    manufacturer: str,
    manufacturer_model_name: str,
):
    # Encode pixel data as unsigned integers for improved interoperability
    pixel_data_bits = 8
    dtype = np.dtype(f"uint{pixel_data_bits}")
    slope = 2.0 / (2**pixel_data_bits - 1)
    intercept = -1.0
    rescaled_map = ((input_map - intercept) / slope).astype(dtype)

    rows = source_image.TotalPixelMatrixRows
    cols = source_image.TotalPixelMatrixColumns

    pixel_array = np.zeros((1, rows, cols, 1), dtype=dtype)
    resized_plane = cv2.resize(
        rescaled_map,
        # OpenCV arrays are in column-major order: (width, height)
        dsize=(cols, rows),
        interpolation=cv2.INTER_CUBIC,
    )
    pixel_array = resized_plane

    value_range = [pixel_array.min(), pixel_array.max()]
    window_width = value_range[1] - value_range[0]
    # Make sure this is an integer
    window_center = int(int(value_range[1]) + int(value_range[0]) // 2)

    real_world_value_mappings = [
        hd.pm.RealWorldValueMapping(
            lut_label="0",
            lut_explanation=class_label,
            unit=codes.UCUM.NoUnits,
            value_range=value_range,
            intercept=intercept,
            slope=slope,
        )
    ]

    parametric_map = hd.pm.ParametricMap(
        [source_image],
        pixel_array=pixel_array,
        series_instance_uid=hd.UID(),
        series_number=randint(1, 100),
        sop_instance_uid=hd.UID(),
        instance_number=1,
        manufacturer=manufacturer,
        manufacturer_model_name=manufacturer_model_name,
        software_versions="0.0.1",
        device_serial_number=f"{gethostname()}",
        contains_recognizable_visual_features=False,
        real_world_value_mappings=real_world_value_mappings,
        window_center=window_center,
        window_width=window_width,
        content_label="HEATMAP",
    )

    return parametric_map

Let me know if you have any further questions!