Improve RTSTRUCT contour slice matching

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

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

matRad/dicom/matRad_findRtssContourSlicesInCt.m

@@ -0,0 +1,84 @@
+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.
+%
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+matRad_cfg = MatRad_Config.instance();
+
+if ~isnumeric(contourZ) || ~isscalar(contourZ) || ~isfinite(contourZ)
+    matRad_cfg.dispError('contourZ must be a finite numeric scalar.');
+end
+
+sliceAxis = matRad_getCtSliceAxis(ct);
+sliceThickness = matRad_getCtSliceThickness(ct);
+
+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)
+
+matRad_cfg = MatRad_Config.instance();
+
+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
+    matRad_cfg.dispError('ct must contain a numeric z-axis or dicomInfo.SlicePositions.');
+end
+
+if any(~isfinite(sliceAxis))
+    matRad_cfg.dispError('CT slice positions must be finite.');
+end
+end
+
+function sliceThickness = matRad_getCtSliceThickness(ct)
+
+matRad_cfg = MatRad_Config.instance();
+
+if ~isstruct(ct) || ~isfield(ct, 'dicomInfo') || ~isstruct(ct.dicomInfo) || ...
+        ~isfield(ct.dicomInfo, 'SliceThickness')
+    matRad_cfg.dispError('ct.dicomInfo.SliceThickness is required.');
+end
+
+sliceThickness = ct.dicomInfo.SliceThickness;
+if ~isnumeric(sliceThickness) || ~isscalar(sliceThickness) || ...
+        ~isfinite(sliceThickness) || sliceThickness <= 0
+    matRad_cfg.dispError('ct.dicomInfo.SliceThickness must be a positive numeric scalar.');
+end
+end

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];