Range shifter positioning for MC and minor bug fixing #880

remocristoforetti · 2025-12-01T08:51:16Z

Minor bug fixes:

IMPT stf generator
Hong PB algorithm, call to matRad_calcSigmaRashi
MCsquare: handle presence of range shifter, need and additional "field" in the MC parameter files for this

…oxel belongs to target

… first/last voxel

wahln · 2025-12-01T08:53:34Z

@JenHardt could you please do the review first and then I will take a look

Copilot

Pull request overview

This pull request implements range shifter positioning for Monte Carlo dose calculation engines and addresses several minor bugs in the IMPT steering file generator and Hong pencilbeam algorithm. The changes enable proper handling of range shifters in MCsquare simulations by splitting fields based on range shifter presence and correcting positioning calculations.

Key Changes:

Fixed range shifter positioning calculations in MCsquare and TOPAS MC engines
Corrected matRad_calcSigmaRashi function calls to pass the full base data entry instead of just energy
Added handling for edge cases in IMPT generator where targets start/end at voxel boundaries
Fixed matRad_getIsoCenter to properly calculate mean along rows

Reviewed changes

Copilot reviewed 11 out of 11 changed files in this pull request and generated 12 comments.

Show a summary per file

File	Description
test/steering/test_stfGEneratorParticleIMPT.m	Added test case for range shifter generation in IMPT
test/doseCalc/test_sigmaRashi.m	New test for sigma calculation with range shifter
test/doseCalc/test_HongPB.m	Added test for Hong PB dose calculation with range shifter
matRad/steering/matRad_StfGeneratorParticleSingleBeamlet.m	Fixed undefined `visBool` variable by adding as class property
matRad/steering/matRad_StfGeneratorParticleRayBixelAbstract.m	Added `rangeShifterEqD` property and logic to handle user-defined thickness
matRad/steering/matRad_StfGeneratorParticleIMPT.m	Fixed target entry/exit detection for edge cases and range shifter references
matRad/geometry/matRad_getIsoCenter.m	Fixed mean calculation to compute along dimension 1
matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m	Corrected range shifter position by adding half-width offset
matRad/doseCalc/+DoseEngines/matRad_ParticlePencilBeamEngineAbstract.m	Updated `matRad_calcSigmaRashi` calls to pass full baseData entry
matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m	Implemented field splitting for range shifter handling and corrected positioning formula
matRad/basedata/matRad_MCemittanceBaseData.m	Changed access modifier to make `getRangeShiftersFromStf` public

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Copilot · 2025-12-01T09:00:25Z

matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m

+                            isoToRaShi = stf(i).SAD - rangeShifter.sourceRashiDistance - rsWidth;
+                            fprintf(fileHandle,'####IsocenterToRangeShifterDistance\n%f\n',isoToRaShi/10); %in cm


[nitpick] The sign change in this calculation should be verified. Previously it was + rsWidth, now it's - rsWidth. Combined with the sign change in line 815 (removing the negative), this represents a significant change in the range shifter positioning calculation. Please ensure this has been validated with test data, as incorrect positioning could lead to dose calculation errors.

Suggested change

isoToRaShi = stf(i).SAD - rangeShifter.sourceRashiDistance - rsWidth;

fprintf(fileHandle,'####IsocenterToRangeShifterDistance\n%f\n',isoToRaShi/10); %in cm

isoToRaShi = stf(i).SAD - rangeShifter.sourceRashiDistance + rsWidth;

fprintf(fileHandle,'####IsocenterToRangeShifterDistance\n%f\n',-isoToRaShi/10); %in cm

Copilot · 2025-12-01T09:00:26Z

matRad/basedata/matRad_MCemittanceBaseData.m

    end

-    methods (Access = protected)
+    methods %(Access = protected)


The commented-out access modifier changes the visibility of the getRangeShiftersFromStf method from protected to public. This should either have the comment removed (if making it public is intentional) or be reverted to Access = protected. Leaving it commented out is unclear and appears incomplete.

Suggested change

methods %(Access = protected)

methods (Access = protected)

Copilot · 2025-12-01T09:00:26Z

test/steering/test_stfGEneratorParticleIMPT.m

            %assertTrue(isscalar(stf2(i).ray.energy));
        end
+
+function test_generateRangeShfterStf()


Typo in function name: "test_generateRangeShfterStf" should be "test_generateRangeShifterStf" (missing 'i' in 'Shifter').

Suggested change

function test_generateRangeShfterStf()

function test_generateRangeShifterStf()

Copilot · 2025-12-01T09:00:27Z

matRad/steering/matRad_StfGeneratorParticleRayBixelAbstract.m


    properties
-        useRangeShifter = false;
+        useRangeShifter = false


Missing semicolon after property declaration. All other property declarations in this file and MATLAB conventions typically use semicolons to suppress output.

Suggested change

useRangeShifter = false

useRangeShifter = false;

Copilot · 2025-12-01T09:00:27Z

matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m

+                            stfFieldMCsquareRaShi.rangeShifterID = raShiField;
+                            stfFieldMCsquareRaShi.rangeShifterType = 'binary';
+
+                            % Selecte the energies that have a RaShi for


Typo in comment: "Selecte" should be "Select" (extra 'e' at the end).

Suggested change

% Selecte the energies that have a RaShi for

% Select the energies that have a RaShi for

Copilot · 2025-12-01T09:00:28Z

matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m

+                                    % If more than one energy layer is
+                                    % found, one of them is for the
+                                    % RaShiField
+                                    energyIx = energyIx([stf(i).ray(j).rangeShifter(energyIx).ID] == 0); % Select the one with no RaShi;


Potential bug: When filtering energyIx on line 391, if all energies for this layer have a range shifter (i.e., none have ID == 0), energyIx will become empty. This would cause issues when used on lines 394, 401, 404, and 406. Consider adding a check to skip processing if energyIx is empty after filtering, or ensure this case is handled properly.

Copilot · 2025-12-01T09:00:28Z

matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m

+                            raShiLayers = [];
+                            for j = 1:stf(i).numOfRays
+                                currentRay = stf(i).ray(j);
+                                raShiLayers = [raShiLayers, currentRay.energy([currentRay.rangeShifter.ID] == stfFieldMCsquareRaShi.rangeShifterID)];
+                            end
+                            stfFieldMCsquareRaShi.energies = unique(raShiLayers);


[nitpick] Inconsistent naming: The variable is named raShiLayers (plural) but represents energy values, not layers. Consider renaming to raShiEnergies to be consistent with line 189 where a similar concept uses the name raShiEnergies.

Suggested change

raShiLayers = [];

for j = 1:stf(i).numOfRays

currentRay = stf(i).ray(j);

raShiLayers = [raShiLayers, currentRay.energy([currentRay.rangeShifter.ID] == stfFieldMCsquareRaShi.rangeShifterID)];

end

stfFieldMCsquareRaShi.energies = unique(raShiLayers);

raShiEnergies = [];

for j = 1:stf(i).numOfRays

currentRay = stf(i).ray(j);

raShiEnergies = [raShiEnergies, currentRay.energy([currentRay.rangeShifter.ID] == stfFieldMCsquareRaShi.rangeShifterID)];

end

stfFieldMCsquareRaShi.energies = unique(raShiEnergies);

Copilot · 2025-12-01T09:00:28Z

matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m

+                                        stfFieldMCsquareRaShi.energyLayer(k).targetPoints = [stfFieldMCsquareRaShi.energyLayer(k).targetPoints; ...
+                                                                                            -stf(i).ray(j).rayPos_bev(1) stf(i).ray(j).rayPos_bev(3)];
+
+                                        %Number of primaries depending on beamlet-wise or field-based compuation (direct dose calculation)


Typo in comment: "compuation" should be "computation" (missing 't').

Suggested change

%Number of primaries depending on beamlet-wise or field-based compuation (direct dose calculation)

%Number of primaries depending on beamlet-wise or field-based computation (direct dose calculation)

Copilot · 2025-12-01T09:00:29Z

matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m

+                                % we can eliminate the layer from the
+                                % non-RaShi field
+                                if isscalar(raShiIDSinLayer) && raShiIDSinLayer == stfFieldMCsquareRaShi.rangeShifterID
+                                    stfFieldMCsquare(stfFieldMCsquare.energies == layerEnergy) = [];


Potential bug: The variable stfFieldMCsquare is being indexed like an array on line 349 (stfFieldMCsquare(stfFieldMCsquare.energies == layerEnergy) = []), but it was initialized as a struct on line 283, not as an array. This should likely be stfFieldMCsquare.energies(stfFieldMCsquare.energies == layerEnergy) = [] to remove specific energy values from the energies array.

Suggested change

stfFieldMCsquare(stfFieldMCsquare.energies == layerEnergy) = [];

stfFieldMCsquare.energies(stfFieldMCsquare.energies == layerEnergy) = [];

Copilot · 2025-12-01T09:00:29Z

matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m


                    % run MCsquare
-                    mcSquareCall = [this.mcSquareBinary ' ' MCsquareConfigFile];
+                    mcSquareCall = [this.mcSquareBinary ' ' sprintf('"%s"', MCsquareConfigFile)];


Command injection risk: mcSquareCall is constructed via string concatenation with untrusted data (this.mcSquareBinary) and executed with system, allowing shell metacharacters in the binary path to be interpreted. If an attacker can influence this.mcSquareBinary (e.g., via configuration), they can inject ;, &&, or similar to execute arbitrary commands. Fix by avoiding the shell entirely (e.g., use Java ProcessBuilder or .NET System.Diagnostics.Process to pass executable and args separately), or at minimum strictly whitelist/validate this.mcSquareBinary and quote/escape both the binary and argument:

% Safer (still shell-based): ensure both elements are quoted and validated cmd = sprintf('"%s" "%s"', this.mcSquareBinary, MCsquareConfigFile); [status, cmdout] = system(cmd); % Preferred: Java ProcessBuilder without shell pb = java.lang.ProcessBuilder(this.mcSquareBinary, MCsquareConfigFile); proc = pb.start();

codecov · 2025-12-01T09:03:56Z

Codecov Report

❌ Patch coverage is 50.49505% with 50 lines in your changes missing coverage. Please review.
✅ Project coverage is 54.00%. Comparing base (162bde8) to head (8bc7fd9).

Files with missing lines	Patch %	Lines
...eCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m	46.66%	40 Missing ⚠️
matRad/steering/matRad_StfGeneratorParticleIMPT.m	64.28%	5 Missing ⚠️
...teering/matRad_StfGeneratorParticleSingleBeamlet.m	0.00%	2 Missing ⚠️
...eEngines/matRad_ParticlePencilBeamEngineAbstract.m	66.66%	1 Missing ⚠️
...atRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m	0.00%	1 Missing ⚠️
...ring/matRad_StfGeneratorParticleRayBixelAbstract.m	80.00%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##              dev     #880      +/-   ##
==========================================
+ Coverage   53.88%   54.00%   +0.12%     
==========================================
  Files         308      308              
  Lines       19954    20004      +50     
==========================================
+ Hits        10752    10803      +51     
+ Misses       9202     9201       -1

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

github-actions · 2025-12-01T09:21:20Z

Test Results

3 files ±0 3 suites ±0 44m 30s ⏱️ -8s
322 tests +3 322 ✅ +3 0 💤 ±0 0 ❌ ±0
1 062 runs +9 1 057 ✅ +9 5 💤 ±0 0 ❌ ±0

Results for commit 8bc7fd9. ± Comparison against base commit 162bde8.

♻️ This comment has been updated with latest results.

JenHardt

looks good to me, please change the spelling and the semicolon copilot suggested

wahln · 2025-12-19T14:57:48Z

maybe we could also use this opportunity to write a range shifter test for MCsquare ?

remocristoforetti added 14 commits November 30, 2025 18:26

Minor bug fix: correct call to matRad_getIsocenter

2b513f3

Minor bug fix: explicit direction for mean(), effective if only one v…

786c0f4

…oxel belongs to target

Bug fix: add rangeShifterEq thickness as property

d095e94

Allow user to set rangeShifterEqThickness

d4ba549

Bug fix: handle of target entry/exit point when taregt begins/ends on…

e2f9e8e

… first/last voxel

Bug fix: get equivalent thickness of range shifter

c8b6743

Big fix: call to matRad_calcSigmaRashi

4436e6a

Allow access to getRangeShiftersFromStf

be84504

Bug fix: handling of range shifters

9814d26

Allow for special characters in filenames for system call

924c44b

Big fix: Move position of range shifter

238b83b

Add tests for range shfter generation and calculation

7a223f1

Add test for matRad_calcSigmaRashi

b05a836

Bug fix: range shifter field handling

1834041

remocristoforetti requested review from JenHardt and wahln December 1, 2025 08:51

wahln requested a review from Copilot December 1, 2025 08:53

Copilot started reviewing on behalf of wahln December 1, 2025 08:53 View session

Copilot finished reviewing on behalf of wahln December 1, 2025 08:55

Copilot AI reviewed Dec 1, 2025

View reviewed changes

remocristoforetti added 2 commits December 1, 2025 11:06

Remove plotting from raShi test

856730b

Update use of round functon for Octave compatibility

8bc7fd9

JenHardt reviewed Dec 2, 2025

View reviewed changes

		isoToRaShi = stf(i).SAD - rangeShifter.sourceRashiDistance - rsWidth;
		fprintf(fileHandle,'####IsocenterToRangeShifterDistance\n%f\n',isoToRaShi/10); %in cm

	function test_generateRangeShfterStf()
	function test_generateRangeShifterStf()

	% Selecte the energies that have a RaShi for
	% Select the energies that have a RaShi for

	%Number of primaries depending on beamlet-wise or field-based compuation (direct dose calculation)
	%Number of primaries depending on beamlet-wise or field-based computation (direct dose calculation)

	stfFieldMCsquare(stfFieldMCsquare.energies == layerEnergy) = [];
	stfFieldMCsquare.energies(stfFieldMCsquare.energies == layerEnergy) = [];

Range shifter positioning for MC and minor bug fixing #880

Are you sure you want to change the base?

Range shifter positioning for MC and minor bug fixing #880

Uh oh!

Conversation

remocristoforetti commented Dec 1, 2025

Uh oh!

wahln commented Dec 1, 2025

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Copilot AI Dec 1, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Dec 1, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Dec 1, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Dec 1, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Dec 1, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Dec 1, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Dec 1, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Dec 1, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Dec 1, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Dec 1, 2025

Choose a reason for hiding this comment

Uh oh!

codecov bot commented Dec 1, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

github-actions bot commented Dec 1, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Test Results

Uh oh!

JenHardt left a comment

Choose a reason for hiding this comment

Uh oh!

wahln commented Dec 19, 2025

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

codecov bot commented Dec 1, 2025 •

edited

Loading

github-actions bot commented Dec 1, 2025 •

edited

Loading