question about the dij structure

[chh105]

What is the question? Which file(s) does it concern?
I had a question about combining dij structures from two separate forward calculations. For example, say I compute dij for a one beam plan with (couch,gantry) = (0,90) and another dij for a one beam plan with (0,-90). If I wanted to save on compute for a two beam plan with couch angles [0,0] and gantry angles [90,-90], would it be possible to combine the dij structures?

If you can combine dij structures, would I just concatenate the physical dose matrices in the column-wise direction, concatenate the other vectors (i.e. bixelNum, rayNum, beamNum) in the row-wise direction, etc?

I have checked the following matRad Wiki pages and references:
documentation for the dij structure

[markbangert]

In principle you can combine DIJ structures by column-wise concatenation. Thereby you kind of increase the number of bixels. Be careful that you also modify all meta information (numberOfBixels etc.) accordingly if you want to use the matrix for optimization within matRad.

However, I do not really see how this maked sense in the context of a forward calculation. There you can simply add the dose cubes?!

[chh105]

Thanks for the quick response. Also sorry if I didnt describe my use case very well, but I'm searching through a bunch of candidate angles and iteratively adding the best angle to create a n-beam plan. Most of the computation time is spent computing dij for each beam, so I was hoping to save and reuse the information

[markbangert]

In case you want to use matRad for beam angle selection you certainly should precompute the Dij matrices and reuse them. Ideally you establish some grid before and then glue them together for the combinations that you want to investigate. For such a use case, however, it might even make sense to write your own backpropagation function because the sparse matrix concatenation will also take quite some time.

Also you may want to check out approaches like https://pubmed.ncbi.nlm.nih.gov/26936695/ to save something on the optimization side.

[chh105]

Very interesting. I'll probably start by just concatenating and see if it runs fast enough. For the optimization part, I've just been running IPOPT for ~15 iterations which finishes relatively quickly (I think this is similar to the early stopping method in your paper). I've also been using an objective function with a squared deviation for the ptv and squared overdosing to the body (with equal weights). I was wondering if you tried various objective functions besides squared over/under dosing and whether there was much difference between the ones you tried.

[markbangert]

No we have not experimented with different obj functions. I have investigated different methods on the beam angle selection problem itself, though. Iterative adding, however, was pretty much clinically equivalent.

https://aapm.onlinelibrary.wiley.com/doi/abs/10.1118/1.4771932

[chh105]

Ok I'll probably just stick with the current optimization approach then. Anyways thanks for your help!

[chh105]

Just had a follow up question about combining physicalDose matrices. Say we compute two physicalDose matrices for (couch=0, gantry=90) and (couch=0, gantry=-90) and combine them by concatenating column-wise. We can also compute a physicalDose matrix for (couch=[0,0], gantry=90,-90). I would expect these two matrices to be the same (isEqual returns 1), but it seems that they are very slightly different (isEqual returns 0 and the number of nonzero terms in the physicalDose matrices differs by around 100). Is this acceptable/expected behavior?

[markbangert]

that is interesting. are the steering setups the same? meaning isocenter, and position of beamlets, number of bixels, etc. I would not expect them to be binary identical but right now I do not see why they should not have the same non-zero pattern.

…

On Thu, Jan 14, 2021 at 10:11 PM chh105 ***@***.***> wrote: Just had a follow up question about combining physicalDose matrices. Say we compute two physicalDose matrices for (couch=0, gantry=90) and (couch=0, gantry=-90) and combine them by concatenating column-wise. We can also compute a physicalDose matrix for (couch=[0,0], gantry=90,-90). I would expect these two matrices to be the same (isEqual returns 1), but it seems that they are very slightly different (isEqual returns 0 and the number of nonzero terms in the physicalDose matrices differs by around 100). Is this acceptable/expected behavior? — You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub <#477 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACNPY7KRR45XFVX45FAQFBLSZ5MY7ANCNFSM4V6AZXLQ> .

[chh105]

I think the pln and stf structures should be unmodified. To reproduce this manually, I setup everything in the GUI, computed three dij structures, and then compared the two physicalDose matrices (one where things are combined and one copied from the dij structure). I'm wondering if maybe I overlooked something

[chh105]

Edit: Actually sorry after checking again, I think the pln and stf structures are being updated by the gui. I think this might be causing the isocenter to be different for the precomputed and regular dij calculations. Just so to clarify, should the couch/gantry angles be different for the precomputed/regular dij calculations while the isocenter remains the same?

[chh105]

I think I'm just making this more confusing with my comments, so I wrote a script to reproduce things.

gantry = [90 -90];
couch = [0 0];
pln.propStf.gantryAngles = gantry;
pln.propStf.couchAngles = couch;
pln.propStf.numOfBeams = numel(pln.propStf.gantryAngles);
pln.propStf.isoCenter = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0);
stf = matRad_generateStf(ct,cst,pln);

% compute dij for first angle
pln.propStf.gantryAngles = gantry(1);
pln.propStf.couchAngles = couch(1);
pln.propStf.numOfBeams = numel(pln.propStf.gantryAngles);
stf = matRad_generateStf(ct,cst,pln);
if strcmp(pln.radiationMode,'photons')
    dij_1 = matRad_calcPhotonDose(ct,stf,pln,cst);
    %dij = matRad_calcPhotonDoseVmc(ct,stf,pln,cst);
elseif strcmp(pln.radiationMode,'protons') || strcmp(pln.radiationMode,'carbon')
    dij_1 = matRad_calcParticleDose(ct,stf,pln,cst);
end

% compute dij for second angle
pln.propStf.gantryAngles = gantry(2);
pln.propStf.couchAngles = couch(2);
pln.propStf.numOfBeams = numel(pln.propStf.gantryAngles);
stf = matRad_generateStf(ct,cst,pln);
if strcmp(pln.radiationMode,'photons')
    dij_2 = matRad_calcPhotonDose(ct,stf,pln,cst);
    %dij = matRad_calcPhotonDoseVmc(ct,stf,pln,cst);
elseif strcmp(pln.radiationMode,'protons') || strcmp(pln.radiationMode,'carbon')
    dij_2 = matRad_calcParticleDose(ct,stf,pln,cst);
end

% compute dij for 2 beam config
pln.propStf.gantryAngles = gantry;
pln.propStf.couchAngles = couch;
pln.propStf.numOfBeams = numel(pln.propStf.gantryAngles);
stf = matRad_generateStf(ct,cst,pln);
if strcmp(pln.radiationMode,'photons')
    dij = matRad_calcPhotonDose(ct,stf,pln,cst);
    %dij = matRad_calcPhotonDoseVmc(ct,stf,pln,cst);
elseif strcmp(pln.radiationMode,'protons') || strcmp(pln.radiationMode,'carbon')
    dij = matRad_calcParticleDose(ct,stf,pln,cst);
end

% precomputed dose mat
dose_mat = [];
dose_mat = [dose_mat,dij_1.physicalDose{1}];
dose_mat = [dose_mat,dij_2.physicalDose{1}];
% regular dose mat
regular_dose_mat = dij.physicalDose{1};

%% comparisons
subplot(2,1,1)
spy(dose_mat)
subplot(2,1,2)
spy(regular_dose_mat)
sprintf('NNZ precomputed:%d',nnz(dose_mat))
sprintf('NNZ regular:%d',nnz(regular_dose_mat))
diff = full(sum(dose_mat(:)-regular_dose_mat(:))^2/sum(dose_mat(:))^2);
sprintf('Diff:%e',diff)

[markbangert]

Sorry for the late reply. Somehow your issue got covered in my mailbox. Are you still trying to figure this out?

…

On Thu, Jan 14, 2021 at 11:54 PM chh105 ***@***.***> wrote: I think I'm just making this more confusing with my comments, so I wrote a script to reproduce things. debug_precomputed_dij.zip <https://github.com/e0404/matRad/files/5817397/debug_precomputed_dij.zip> — You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub <#477 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACNPY7ICDLQEYOUC27S3CT3SZ5Y27ANCNFSM4V6AZXLQ> .

[chh105]

Yes I still haven't figured it out. There seems to be a small difference between the concatenated physical dose matrices and the regular physical dose matrix.

[markbangert]

Thank you for providing the script. The behaviour you observe is related to the dose influence matrix sampling. In order to save memory we perform a random sampling of the influence data in

matRad/matRad_calcPhotonDose.m

Line 279 in 6fad1b8

[ix,bixelDose] = matRad_DijSampling(ix,bixelDose,radDepthVdoseGrid{1}(ix),rad_distancesSq,Type,r0);

. As we are re-setting the state of the random engine only in the beginning of the entire dose calculation at

matRad/matRad_calcPhotonDose.m

Line 45 in 6fad1b8

rng('default');

we have the funny behaviour that the first beams are binary identical while the second beam differs. If you want to make this consistent you can reset the random engine before the calculation of every beam. @wahln what do you think about including this for the main release? I am somewhat hesitant because resetting the random engine is always fishy with regard to statistical independence... But we could even think about setting the random engine based on a hash computed on the index set of the beamlets within the dij sampling function itself. That way we would even have reproducability for individual bixels given they have the exact same geometry...

[chh105]

Thanks for figuring this out so quickly. Just to clarify, are you saying that the random engine needs to be reset for each beam when computing dij for a multiple beam plan? For example, would you paste rng('default'); after the for loop call in line 143 (for i = 1:dij.numOfBeams % loop over all beams)?

[markbangert]

yes! without having tried this out... :) if you reset the random engine before every beam, the dij sampling pattern should be the same for every beam as long as the stf structs are identical - which they are in your example. this I have checked with matlab's isequal function.

…

On Tue, Jan 26, 2021 at 11:49 PM chh105 ***@***.***> wrote: Thanks for figuring this out so quickly. Just to clarify, are you saying that the random engine needs to be reset for each beam when computing dij for a multiple beam plan? For example, would you paste rng('default'); after the for loop call in line 143 (for i = 1:dij.numOfBeams % loop over all beams)? — You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub <#477 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACNPY7N2LD6BS776CAAFSYLS35BGFANCNFSM4V6AZXLQ> .

[wahln]

Nice that you guys figured that out. I think we should allow some configuration of the sampling that enforces reproducible / consistent results. I will think about some configuration parameters in pln.propDoseCalc that allow people to enable/disable & steer the sampling process / set seeds & choose seed mechanic (per beam etc.).

[stale]

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. Thank you for your contributions.