From b448fcbf754017fca50a2bdf948a64861b9311ce Mon Sep 17 00:00:00 2001 From: Camilo Sevilla Date: Tue, 19 May 2026 12:35:56 -0500 Subject: [PATCH] Improve RTSTRUCT contour slice matching --- .../matRad_importDicomRtss.m | 7 +- .../dicom/matRad_convRtssContours2Indices.m | 32 ++++---- .../dicom/matRad_findRtssContourSlicesInCt.m | 78 +++++++++++++++++++ test/dicom/test_rtssContourSliceMatching.m | 71 +++++++++++++++++ 4 files changed, 167 insertions(+), 21 deletions(-) create mode 100644 matRad/dicom/matRad_findRtssContourSlicesInCt.m create mode 100644 test/dicom/test_rtssContourSliceMatching.m diff --git a/matRad/dicom/@matRad_DicomImporter/matRad_importDicomRtss.m b/matRad/dicom/@matRad_DicomImporter/matRad_importDicomRtss.m index 63c6b0aed..ba93b885e 100644 --- a/matRad/dicom/@matRad_DicomImporter/matRad_importDicomRtss.m +++ b/matRad/dicom/@matRad_DicomImporter/matRad_importDicomRtss.m @@ -126,8 +126,11 @@ end % sanity check 2 - if unique(structZ) > max(obj.ct.dicomInfo.SlicePositions) || unique(structZ) < min(obj.ct.dicomInfo.SlicePositions) - matRad_cfg.dispWarning(['Omitting contour data for ' obj.importRtss.structures(i).structName ' at slice position ' num2str(unique(structZ)) 'mm - no ct data available.\n']); + contourZ = unique(structZ); + if isempty(matRad_findRtssContourSlicesInCt(contourZ, obj.ct)) + matRad_cfg.dispWarning(['Omitting contour data for ' obj.importRtss.structures(i).structName ... + ' at slice position ' num2str(contourZ) ... + 'mm - no ct data available.\n']); else obj.importRtss.structures(i).item(j).points = [structX, structY, structZ]; end diff --git a/matRad/dicom/matRad_convRtssContours2Indices.m b/matRad/dicom/matRad_convRtssContours2Indices.m index 8f436e19f..9d9d2aade 100644 --- a/matRad/dicom/matRad_convRtssContours2Indices.m +++ b/matRad/dicom/matRad_convRtssContours2Indices.m @@ -43,17 +43,15 @@ matRad_cfg.dispError('Contour defined over multiple planes!'); end - round2 = @(a,b) round(a*10^b)/10^b; - dicomCtSliceThickness = ct.dicomInfo.SliceThickness; - - %Sanity check - msg = checkSliceThickness(dicomCtSliceThickness); - if ~isempty(msg) - matRad_cfg.dispError('Slice Thickness of slice at %f could not be identified: %s',dicomCtSlicePos,msg); + slicesInMatradCt = matRad_findRtssContourSlicesInCt(dicomCtSlicePos, ct); + if isempty(slicesInMatradCt) + structureName = matRad_getStructureName(structure); + matRad_cfg.dispWarning(['Omitting contour data for ' structureName ... + ' at slice position ' num2str(dicomCtSlicePos) ... + 'mm - no ct data available.\n']); + continue end - slicesInMatradCt = find(dicomCtSlicePos+dicomCtSliceThickness/2 > ct.z & dicomCtSlicePos-dicomCtSliceThickness/2 <= ct.z); - coords1 = interp1(ct.x,1:ct.cubeDim(2),structure.item(i).points(:,1),'linear','extrap'); coords2 = interp1(ct.y,1:ct.cubeDim(1),structure.item(i).points(:,2),'linear','extrap'); @@ -72,14 +70,10 @@ end -function msg = checkSliceThickness(dicomCtSliceThickness) - if isempty(dicomCtSliceThickness) - msg = 'Slice could not be identified (empty)'; - elseif ~isscalar(dicomCtSliceThickness) - msg = 'Slice thickness not unique'; - elseif ~isnumeric(dicomCtSliceThickness) - msg = 'unexpected value'; - else - msg = ''; - end +function structureName = matRad_getStructureName(structure) +if isfield(structure, 'structName') && ~isempty(structure.structName) + structureName = structure.structName; +else + structureName = 'structure'; +end end diff --git a/matRad/dicom/matRad_findRtssContourSlicesInCt.m b/matRad/dicom/matRad_findRtssContourSlicesInCt.m new file mode 100644 index 000000000..b5f50ac54 --- /dev/null +++ b/matRad/dicom/matRad_findRtssContourSlicesInCt.m @@ -0,0 +1,78 @@ +function sliceIndices = matRad_findRtssContourSlicesInCt(contourZ, ct) +% matRad function to find CT slices compatible with an RTSTRUCT contour plane +% +% call +% sliceIndices = matRad_findRtssContourSlicesInCt(contourZ,ct) +% +% input +% contourZ: z-position of one RTSTRUCT contour plane in mm +% ct: matRad ct struct with z-axis and DICOM slice thickness +% +% output +% sliceIndices: row vector with CT slice indices intersected by the +% physical slice slab of the contour plane +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2026 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSE.md. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +matRadCfg = MatRad_Config.instance(); + +if ~isnumeric(contourZ) || ~isscalar(contourZ) || ~isfinite(contourZ) + matRadCfg.dispError('contourZ must be a finite numeric scalar.'); +end + +sliceAxis = matRad_getCtSliceAxis(ct, matRadCfg); +sliceThickness = matRad_getCtSliceThickness(ct, matRadCfg); + +halfSliceThickness = sliceThickness / 2; +tol = max(1e-6, 1e-6 * sliceThickness); + +sliceIndices = find(contourZ + halfSliceThickness + tol > sliceAxis & ... + contourZ - halfSliceThickness - tol <= sliceAxis); +sliceIndices = sliceIndices(:)'; + +end + +function sliceAxis = matRad_getCtSliceAxis(ct, matRadCfg) +if isstruct(ct) && isfield(ct, 'z') && isnumeric(ct.z) && isvector(ct.z) && ... + ~isempty(ct.z) + sliceAxis = ct.z(:)'; +elseif isstruct(ct) && isfield(ct, 'dicomInfo') && isstruct(ct.dicomInfo) && ... + isfield(ct.dicomInfo, 'SlicePositions') && ... + isnumeric(ct.dicomInfo.SlicePositions) && isvector(ct.dicomInfo.SlicePositions) && ... + ~isempty(ct.dicomInfo.SlicePositions) + sliceAxis = ct.dicomInfo.SlicePositions(:)'; +else + matRadCfg.dispError('ct must contain a numeric z-axis or dicomInfo.SlicePositions.'); +end + +if any(~isfinite(sliceAxis)) + matRadCfg.dispError('CT slice positions must be finite.'); +end +end + +function sliceThickness = matRad_getCtSliceThickness(ct, matRadCfg) +if ~isstruct(ct) || ~isfield(ct, 'dicomInfo') || ~isstruct(ct.dicomInfo) || ... + ~isfield(ct.dicomInfo, 'SliceThickness') + matRadCfg.dispError('ct.dicomInfo.SliceThickness is required.'); +end + +sliceThickness = ct.dicomInfo.SliceThickness; +if ~isnumeric(sliceThickness) || ~isscalar(sliceThickness) || ... + ~isfinite(sliceThickness) || sliceThickness <= 0 + matRadCfg.dispError('ct.dicomInfo.SliceThickness must be a positive numeric scalar.'); +end +end diff --git a/test/dicom/test_rtssContourSliceMatching.m b/test/dicom/test_rtssContourSliceMatching.m new file mode 100644 index 000000000..0dfeaaa5b --- /dev/null +++ b/test/dicom/test_rtssContourSliceMatching.m @@ -0,0 +1,71 @@ +function test_suite = test_rtssContourSliceMatching + +test_functions = localfunctions(); + +initTestSuite; + +function test_acceptsBoundaryContourInsideSliceSlab +ct = helper_ctFixture(); + +sliceIndices = matRad_findRtssContourSlicesInCt(91.530000000000, ct); + +assertEqual(sliceIndices, 3); + +function test_rejectsContourOutsideSliceSlab +ct = helper_ctFixture(); +outsideZ = ct.z(end) + ct.dicomInfo.SliceThickness / 2 + 0.01; + +sliceIndices = matRad_findRtssContourSlicesInCt(outsideZ, ct); + +assertTrue(isempty(sliceIndices)); + +function test_invalidSliceThicknessThrowsError +ct = helper_ctFixture(); + +ct.dicomInfo.SliceThickness = []; +assertExceptionThrown(@() matRad_findRtssContourSlicesInCt(0, ct), 'matRad:Error'); + +ct.dicomInfo.SliceThickness = [3 3]; +assertExceptionThrown(@() matRad_findRtssContourSlicesInCt(0, ct), 'matRad:Error'); + +ct.dicomInfo.SliceThickness = '3'; +assertExceptionThrown(@() matRad_findRtssContourSlicesInCt(0, ct), 'matRad:Error'); + +ct.dicomInfo.SliceThickness = 0; +assertExceptionThrown(@() matRad_findRtssContourSlicesInCt(0, ct), 'matRad:Error'); + +function test_voxelizesBoundaryContourOnExtremeSlice +ct = helper_ctFixture(); +structure = helper_structureFixture(91.530000000000); + +indices = matRad_convRtssContours2Indices(structure, ct); +[~, ~, iz] = ind2sub(ct.cubeDim, indices); + +assertTrue(~isempty(indices)); +assertEqual(unique(iz), 3); + +function test_voxelizationOmitsContourOutsideCt +ct = helper_ctFixture(); +outsideZ = ct.z(end) + ct.dicomInfo.SliceThickness / 2 + 0.01; +structure = helper_structureFixture(outsideZ); + +indices = matRad_convRtssContours2Indices(structure, ct); + +assertTrue(isempty(indices)); + +function ct = helper_ctFixture() +ct.x = 1:5; +ct.y = 1:5; +ct.z = [-91.527328693131 0 91.527328693131]; +ct.cubeDim = [5 5 3]; +ct.cubeHU = {zeros(ct.cubeDim)}; +ct.dicomInfo.SlicePositions = ct.z; +ct.dicomInfo.SliceThickness = 3; + +function structure = helper_structureFixture(contourZ) +structure.structName = 'BODY'; +structure.item(1).points = [2 2 contourZ; ... + 4 2 contourZ; ... + 4 4 contourZ; ... + 2 4 contourZ; ... + 2 2 contourZ];