[Gate-users] Problems/inconsistencies calculating dose uncertainty

[Matthew Strugari]

Hi again,

I'm following up with some further calculations/derivations of the relative statistical uncertainty. According to the docs, GATE outputs the relative statistical uncertainty. Using the usual statistical formulae, this would be the standard deviation divided by the dose. Given that the standard deviation is the square root of the sample variance, the final formula should be
        sampleVariance = (1/(N-1)) * (sq_dose - pow(dose,2)/N)
        standardDeviation = np.sqrt(sampleVariance)
        std_dose = standardDeviation/dose
        # or inline:
        std_dose = np.sqrt( (1/(N-1))*(sq_dose - pow(dose,2)/N) )/dose

When comparing my calculated uncertainty to the reference uncertainty, I am now 58% different at the origin (opposed to >8000%) but the average percent difference is 77% (opposed to 1% using the GATE formula). I am overall quite confused how the uncertainty is calculated by GATE, especially considering that the docs (https://opengate.readthedocs.io/en/latest/voxelized_source_and_phantom.html#dose-collection) outline the dose calculations using the sample variance with a formula that represents the population variance.

If anyone would be so kind to instruct me how GATE calculates the uncertainty, and/or how I can reproduce the uncertainty using the Dose, Dose-Squared, and NbOfHits files, it would be greatly appreciated!

Positively,
Matthew



________________________________
From: Gate-users <gate-users-bounces at lists.opengatecollaboration.org> on behalf of Matthew Strugari <matthew.strugari at dal.ca>
Sent: February 6, 2020 2:59 PM
To: gate-users at lists.opengatecollaboration.org <gate-users at lists.opengatecollaboration.org>
Subject: [Gate-users] Problems/inconsistencies calculating dose uncertainty

Hi all,

I am in need of assistance when calculating the dose uncertainty using output from the dose actor. I would like to merge multiple dose point kernel outputs from the dose actor to speed up my overall simulation time. Using the output from a single simulation to verify the calculation of the uncertainty, my results are ~1% different from the reference GATE output on average with a >8000% difference at the origin. I am using the formula found in GateDoseActor.cc:
      double std_dose = sqrt( (1.0/(N-1))*(sq_dose/N - pow(dose/N, 2)) )/(dose/N);
      if( dose == 0.0 || N == 1 || sq_dose == 0 )
        std_dose = 1.0;

Here, the corresponding image data (attached in doseActorOutput.zip) are assigned as dose = Dose.mhd, sq_dose = Dose-Squared.mhd, and N = NbOfHits.mhd and then element-wise operations are performed. My approach is slightly different as I mask the matrices to avoid nan/inf results:
        N = np.ma.masked_equal(N, 0)
        N = np.ma.masked_equal(N, 1)
        dose = np.ma.masked_equal(dose, 0.0)
        sq_dose = np.ma.masked_equal(sq_dose,0)

        # calculate uncertainty using formula from GateDoseActor.cc
        std_dose = np.sqrt( (1.0/(N-1))*(sq_dose/N - pow(dose/N, 2)) )/(dose/N)
        std_dose = np.ma.filled(std_dose, 1.0).astype(np.float32)

Some/very few of the calculated std_dose values are greater than one (1.0000001 so I am not sure if this is indicative of an error or if I should just normalize std_dose before filling the masked values with 1. Also, the Dose-Uncertainty.mhd images appear identical in ImageJ but the actual values appear to be scaled differently (see attached doseKernel.zip).

Please correct me in my calculation or provide any information on the correct way to calculate the uncertainty using the Dose, Dose-Squared, and NbOfHits files.

Positively,
Matthew


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.opengatecollaboration.org/pipermail/gate-users/attachments/20200206/a243e97a/attachment.html>