Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

vol.load_dicom() implementation #30

Merged
merged 8 commits into from
Apr 5, 2021
Merged

vol.load_dicom() implementation #30

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
8f45483
Remove parentheses from Macro
Mar 26, 2021
36a2025
Started implementing from scratch
Mar 30, 2021
2036aef
Outsourcing lps_from_ijk
Mar 30, 2021
e7f87e3
Reconsidered. Back to Manual
Mar 31, 2021
814d6a6
Add description and remove parameter
Mar 31, 2021
06d355e
Change data layout in memory
Apr 3, 2021
63163ad
Removed unnecessary import and minor refactoring
maxrohleder Apr 3, 2021
8425b6e
Remove CUDA commit from this pull-request
maxrohleder Apr 3, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
62 changes: 31 additions & 31 deletions deepdrr/projector/project_kernel.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,69 +52,69 @@ texture<float, 3, cudaReadModeElementType> seg(12);
texture<float, 3, cudaReadModeElementType> seg(13);
#endif

#define UPDATE(multiplier, n) ({\
#define UPDATE(multiplier, n) {\
output[idx + (n)] += (multiplier) * tex3D(volume, px, py, pz) * round(cubicTex3D(seg(n), px, py, pz));\
})
}

#if NUM_MATERIALS == 1
#define INTERPOLATE(multiplier) ({\
#define INTERPOLATE(multiplier) {\
UPDATE(multiplier, 0);\
})
}
#elif NUM_MATERIALS == 2
#define INTERPOLATE(multiplier) ({\
#define INTERPOLATE(multiplier) {\
UPDATE(multiplier, 0);\
UPDATE(multiplier, 1);\
})
}
#elif NUM_MATERIALS == 3
#define INTERPOLATE(multiplier) ({\
#define INTERPOLATE(multiplier) {\
UPDATE(multiplier, 0);\
UPDATE(multiplier, 1);\
UPDATE(multiplier, 2);\
})
}
#elif NUM_MATERIALS == 4
#define INTERPOLATE(multiplier) ({\
#define INTERPOLATE(multiplier) {\
UPDATE(multiplier, 0);\
UPDATE(multiplier, 1);\
UPDATE(multiplier, 2);\
UPDATE(multiplier, 3);\
})
}
#elif NUM_MATERIALS == 5
#define INTERPOLATE(multiplier) ({\
#define INTERPOLATE(multiplier) {\
UPDATE(multiplier, 0);\
UPDATE(multiplier, 1);\
UPDATE(multiplier, 2);\
UPDATE(multiplier, 3);\
UPDATE(multiplier, 4);\
})
}
#elif NUM_MATERIALS == 6
#define INTERPOLATE(multiplier) ({\
#define INTERPOLATE(multiplier) {\
UPDATE(multiplier, 0);\
UPDATE(multiplier, 1);\
UPDATE(multiplier, 2);\
UPDATE(multiplier, 4);\
UPDATE(multiplier, 5);\
})
}
#elif NUM_MATERIALS == 7
#define INTERPOLATE(multiplier) ({\
#define INTERPOLATE(multiplier) {\
UPDATE(multiplier, 0);\
UPDATE(multiplier, 1);\
UPDATE(multiplier, 2);\
UPDATE(multiplier, 4);\
UPDATE(multiplier, 5);\
UPDATE(multiplier, 6);\
})
}
#elif NUM_MATERIALS == 8
#define INTERPOLATE(multiplier) ({\
#define INTERPOLATE(multiplier) {\
UPDATE(multiplier, 0);\
UPDATE(multiplier, 1);\
UPDATE(multiplier, 2);\
UPDATE(multiplier, 4);\
UPDATE(multiplier, 5);\
UPDATE(multiplier, 6);\
UPDATE(multiplier, 7);\
})
}
#elif NUM_MATERIALS == 9
#define INTERPOLATE(multiplier) ({\
#define INTERPOLATE(multiplier) {\
UPDATE(multiplier, 0);\
UPDATE(multiplier, 1);\
UPDATE(multiplier, 2);\
Expand All @@ -123,9 +123,9 @@ texture<float, 3, cudaReadModeElementType> seg(13);
UPDATE(multiplier, 6);\
UPDATE(multiplier, 7);\
UPDATE(multiplier, 8);\
})
}
#elif NUM_MATERIALS == 10
#define INTERPOLATE(multiplier) ({\
#define INTERPOLATE(multiplier) {\
UPDATE(multiplier, 0);\
UPDATE(multiplier, 1);\
UPDATE(multiplier, 2);\
Expand All @@ -135,9 +135,9 @@ texture<float, 3, cudaReadModeElementType> seg(13);
UPDATE(multiplier, 7);\
UPDATE(multiplier, 8);\
UPDATE(multiplier, 9);\
})
}
#elif NUM_MATERIALS == 11
#define INTERPOLATE(multiplier) ({\
#define INTERPOLATE(multiplier) {\
UPDATE(multiplier, 0);\
UPDATE(multiplier, 1);\
UPDATE(multiplier, 2);\
Expand All @@ -148,9 +148,9 @@ texture<float, 3, cudaReadModeElementType> seg(13);
UPDATE(multiplier, 8);\
UPDATE(multiplier, 9);\
UPDATE(multiplierl, 10);\
})
}
#elif NUM_MATERIALS == 12
#define INTERPOLATE(multiplier) ({\
#define INTERPOLATE(multiplier) {\
UPDATE(multiplier, 0);\
UPDATE(multiplier, 1);\
UPDATE(multiplier, 2);\
Expand All @@ -162,9 +162,9 @@ texture<float, 3, cudaReadModeElementType> seg(13);
UPDATE(multiplier, 9);\
UPDATE(multiplier, 10);\
UPDATE(multiplier, 11);\
})
}
#elif NUM_MATERIALS == 13
#define INTERPOLATE(multiplier) ({\
#define INTERPOLATE(multiplier) {\
UPDATE(multiplier, 0);\
UPDATE(multiplier, 1);\
UPDATE(multiplier, 2);\
Expand All @@ -177,9 +177,9 @@ texture<float, 3, cudaReadModeElementType> seg(13);
UPDATE(multiplier, 10);\
UPDATE(multiplier, 11);\
UPDATE(multiplier, 12);\
})
}
#elif NUM_MATERIALS == 14
#define INTERPOLATE(multiplier) ({\
#define INTERPOLATE(multiplier) {\
UPDATE(multiplier, 0);\
UPDATE(multiplier, 1);\
UPDATE(multiplier, 2);\
Expand All @@ -193,11 +193,11 @@ texture<float, 3, cudaReadModeElementType> seg(13);
UPDATE(multiplier, 11);\
UPDATE(multiplier, 12);\
UPDATE(multiplier, 13);\
})
}
#else
#define INTERPOLATE(multiplier) {\
fprintf(stderr, "NUM_MATERIALS not in [1, 14]");\
)
}
#endif

// the CT volume (used to be tex_density)
Expand Down
143 changes: 136 additions & 7 deletions deepdrr/vol.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
import numpy as np
from pathlib import Path
import nibabel as nib
from pydicom.filereader import dcmread, InvalidDicomError

from . import load_dicom
from . import geo
Expand Down Expand Up @@ -56,7 +57,7 @@ def from_parameters(
):
"""Create a volume object with a segmentation of the materials, with its own anatomical coordinate space, from parameters.

Note that the anatomical coordinate system is not the world coordinate system (which is cartesion).
Note that the anatomical coordinate system is not the world coordinate system (which is cartesian).

Suggested anatomical coordinate space units is milimeters.
A helpful introduction to the geometry is can be found [here](https://www.slicer.org/wiki/Coordinate_systems).
Expand All @@ -75,7 +76,7 @@ def from_parameters(

assert spacing.dim == 3

# define anatomical_from_indices FrameTransform
# define anatomical_from_ijk FrameTransform
if anatomical_coordinate_system is None or anatomical_coordinate_system == 'none':
anatomical_from_ijk = geo.FrameTransform.from_scaling(scaling=spacing, translation=origin)
elif anatomical_coordinate_system == 'LPS':
Expand Down Expand Up @@ -159,11 +160,139 @@ def from_nifti(

@classmethod
def from_dicom(
cls,
path: Union[str, Path],
) -> Volume:
"""Create the volume from a DICOM file."""
raise NotImplementedError('load a volume from a dicom file')
cls,
path: Path,
use_thresholding: bool = True,
world_from_anatomical: Optional[geo.FrameTransform] = None,
use_cached: bool = True,
cache_dir: Optional[Path] = None
):
"""
load a volume from a dicom file and compute the anatomical_from_ijk transform from metadata
https://www.slicer.org/wiki/Coordinate_systems
Args:
path: path-like to a multi-frame dicom file. (Currently only Multi-Frame from Siemens supported)
use_thresholding (bool, optional): segment the materials using thresholding (faster but less accurate). Defaults to True.
world_from_anatomical (Optional[geo.FrameTransform], optional): position the volume in world space. If None, uses identity. Defaults to None.
use_cached (bool, optional): [description]. Use a cached segmentation if available. Defaults to True.
cache_dir (Optional[Path], optional): Where to load/save the cached segmentation. If None, use the parent dir of `path`. Defaults to None.

Returns:
Volume: an instance of a deepdrr volume
"""
path = Path(path)
stem = path.name.split('.')[0]

if cache_dir is None:
cache_dir = path.parent

# Multi-frame dicoms store all slices of a volume in one file.
# they must specify the necessary dicom tags under
# https://dicom.innolitics.com/ciods/enhanced-ct-image/enhanced-ct-image-multi-frame-functional-groups
assert path.is_file(), 'Currently only multi-frame dicoms are supported. Path must refer to a file.'
logger.info(f'loading Dicom volume from {path}')

# reading the dicom dataset object
ds = dcmread(path)

# extracting all needed tags TODO add try exepts
frames = ds.PerFrameFunctionalGroupsSequence
shared = ds.SharedFunctionalGroupsSequence[0]
num_slices = len(frames)
first_slice_position = np.array(frames[0].PlanePositionSequence[0].ImagePositionPatient)
last_slice_position = np.array(frames[-1].PlanePositionSequence[0].ImagePositionPatient)
maxrohleder marked this conversation as resolved.
Show resolved Hide resolved
RC = np.array(shared.PlaneOrientationSequence[0].ImageOrientationPatient).reshape(2, 3).T
PixelSpacing = np.array(ds.SharedFunctionalGroupsSequence[0].PixelMeasuresSequence[0].PixelSpacing)
SliceThickness = np.array(ds.SharedFunctionalGroupsSequence[0].PixelMeasuresSequence[0].SliceThickness)
offset = ds.SharedFunctionalGroupsSequence[0].PixelValueTransformationSequence[0].RescaleIntercept
scale = ds.SharedFunctionalGroupsSequence[0].PixelValueTransformationSequence[0].RescaleSlope

# make user aware that this is only tested on windows
if ds.Manufacturer != "SIEMENS":
logger.warning("Multi-frame loading has only been tested on Siemens Enhanced CT DICOMs."
"Please verify everything works as expected.")

# read the 'raw' data array
raw_data = ds.pixel_array.astype(np.float32)
hu_values = raw_data * scale + offset

'''
EXPLANATION - indexing conventions

According to dicom (C.7.6.3.1.4 - Pixel Data) slices are of shape (Rows, Columns)
=> must be (j, i) indexed if we define i == horizontal and j == vertical.
=> we want to conform to the (i, j, k) layout and therefore move the axis of the data array
'''

# convert data to our indexing convention (k, j, i) -> (j, i, k)
hu_values = hu_values.transpose((2, 1, 0)).copy()
benjamindkilleen marked this conversation as resolved.
Show resolved Hide resolved

# transform the volume in HU to densities
data = load_dicom.conv_hu_to_density(hu_values)

# obtain materials analogous to nifti
if use_thresholding:
materials_path = cache_dir / f'{stem}_materials_thresholding.npz'
if use_cached and materials_path.exists():
logger.info(f'found materials segmentation at {materials_path}.')
materials = dict(np.load(materials_path))
else:
logger.info(f'segmenting materials in volume')
materials = load_dicom.conv_hu_to_materials_thresholding(hu_values)
np.savez(materials_path, **materials)
else:
materials_path = cache_dir / f'{stem}_materials.npz'
if use_cached and materials_path.exists():
logger.info(f'found materials segmentation at {materials_path}.')
materials = dict(np.load(materials_path))
else:
logger.info(f'segmenting materials in volume')
materials = load_dicom.conv_hu_to_materials(hu_values)
np.savez(materials_path, **materials)

'''
EXPLANATION - 3d affine transform

DICOM does not offer a 3d transform to locate the voxel data in world space for historic reasons.
However we can construct it from some related DICOM tags. See this resource for more information:
https://nipy.org/nibabel/dicom/dicom_orientation.html

Note, that we do not modify the affine transform to account for the differences in indexing, but
instead modified the data in memory to be in (i, j, k) order.
'''
# construct column for index k
k = np.array((last_slice_position - first_slice_position) / (num_slices - 1)).reshape(3, 1)

# check if the calculated increment matches the SliceThickness (allow .1 millimeters deviations)
assert np.allclose(np.abs(k[2]), SliceThickness, atol=0.1, rtol=0)

# apply scaling to mm
RC_scaled = RC * PixelSpacing

maxrohleder marked this conversation as resolved.
Show resolved Hide resolved
# construct rotation matrix from three columns for (i, j, k)
rot = np.hstack((RC_scaled, k))

# construct affine matrix
affine = np.zeros((4, 4))
affine[:3, :3] = rot # rotation and scaling
affine[:3, 3] = first_slice_position # translation
affine[3, 3] = 1 # homogenous component

# log affine matrix in debug mode
logger.debug(f"manually constructed affine matrix: \n{affine}")
logger.debug(
f"volume_center_xyz : {np.mean([affine @ np.array([*data.shape, 1]), affine @ [0, 0, 0, 1]], axis=0)}")

# cast to FrameTransform
lps_from_ijk = geo.FrameTransform(affine)

# constructing the volume
return cls(
data,
materials,
lps_from_ijk,
world_from_anatomical,
)

def to_dicom(self, path: Union[str, Path]):
"""Write the volume to a DICOM file.
Expand Down