Question about the dose

[justingm]

Hi.

I was looking at the code on how you calculate the dose and I am wondering what is the 0.997 constant value. I am guessing this is the x in the formula but I don't understand where you got that from. Thanks!

Best regards, Justin

[jerichooconnell]

Before you try to make it make sense too much just know that the dose ends up having an imperical correction factor at the end to make sure that it agrees with Topas, so if something looks wrong my hope is that it will still be corrected by the final correction to Topas if that makes sense.

The 0.997 is a density correction for water as the density of water isn't actually 1. Perhaps a little pedantic but it's a fun fact I like to throw in some places to see if things line up better.

[justingm]

Okay, now I understand. However, it is not clear to me whether the dose calculation in fastCAT is still applicable if I have a different geometry. I read the paper and you have M in formula and I'm assuming that that mass is for a 16 cm diameter water or am I confusing things?

[jerichooconnell]

Yes you are right. I don't believe the absolute dose to line up in other geometries. I do think that the noise statistics will be correct if the same dose is used for two different setups. So what you can get is the dose equivalent for a certain number of photons in a CTDI-ish phantom using the verified setup and then using the same number of photons on a different phantom. Not the best but dose is tricky to do without particle transport so I've had trouble implementing a good method.

[justingm]

I actually set a certain dose for my mouse phantom but the resulting CTs were too good in terms of noise so I figured the dose calculation in fastCAT is not applicable when the object and the distances are changed. I am thinking of just doing some simulations to quantify the number of photons that would correspond to the dose that I want. When you did your simulations in TOPAS for the dose, was it just simply scoring the total dose in the water phantom?

Also, I did have some follow up question on an issue that you already closed. Can you also maybe clarify that or should I open another issue? Thanks!

[jerichooconnell]

Hey just comment on the other issue and I'll take a look! I would do just that. If you had a number of photons to mean dose conversion factor I think things would work out! Sorry for having such buggy code

[jerichooconnell]

I also have an example of the type of Topas file I used if you like. I'm quite disorganized on this project at the moment but maybe some lines would be helpful.

Catphan_515.txt

[jerichooconnell]

We are actually working on a dose calculation tool that uses kernel superposition right now, so hopefully we can implement that some time in the future and save you some grief!

[justingm]

I commented on it just now. Don't worry about the buggy code. It actually forced me to really understand it which has been quite educational. Great tool, by the way.

I also do a lot of Monte Carlo simulations so calculating the dose shouldn't be too hard. And it is already a relief that I don't have to do the CT simulation using Monte Carlo because that would take a lot of time.

[jerichooconnell]

Thank you! Yes that has been my experience as well; that it makes you really think deeply about the relationship between the primary and scatter and how energy effects integrating detectors and such. It has been a lot of fun!

[justingm]

Hi. I am a little bit lost again. It says here (see below) that 2e7 is the number in the reference MC simulations. Is that the number of photon counts in the detector that results to the no noise configuration? What is the dose in the water phantom for that case?

https://github.com/jerichooconnell/fastcat/blob/03d32e6aa9fcdff30bb405eae1c18f84bbcaebd3/fastcat/simulate.py#L821-L824

[jerichooconnell]

Sorry for the late reply, AAPM preparation and then AAPM kept me a little busy there.

2e7 is the reference number of initial photons in the MC simulation. From that I get the average number of photons incident per pixel of the detector. Knowing the number of photons that hit the detector I use Poisson statistics on the number of counts. I'm not sure how this articulates with the detector size change and changing the source distance. As I said in another comment we are working on a better dose engine.

I think to fix this I should make a test of a different geometry in the automated tests to make sure that the dose changes correctly with the geometry as well as changing the detector pixel size, which after the stuff you said about re-scaling the fluence with geometry I'm not sure it will.

[justingm]

Hi! Additional question. Can you point me to the part of the code where you then estimate the average number of photons incident per pixel of the detector? Thanks!

[jerichooconnell]

The final photon count is scaled so that the noise is Poisson distributed through the counts, this happens around line 817 of simulate.py

        if test_intensity:
            return intensity

        # Sum over the image dimesions to get the energy
        # intensity and multiply by fluence, 2e7 was number
        # in reference MC simulations
        def get_dose_nphoton(nphot):
            return nphot / 2e7

        def get_dose_mgy(mgy, doses, w_fluence):
            nphoton = mgy / (
                get_dose_per_photon(doses, w_fluence)
                * (1.6021766e-13)
                * 1000
            )
            return get_dose_nphoton(nphoton)

        def get_dose_per_photon(doses, w_fluence):
            # linear fit of the data mod Mar 2021
            pp = np.array([0.87810143, 0.01136471])
            return ((np.array(doses) / 1000) @ (w_fluence)) * pp[0] + pp[1]

        ratio = None

        self.doses = doses
        # Dose in micro grays
        if mgy != 0.0:
            ratio = get_dose_mgy(mgy, np.array(doses), w_fluence)
        elif nphoton is not None:
            ratio = get_dose_nphoton(nphoton)

        # --- Noise and Scatter Calculation ---
        # Now I interpolate deposition
        # and get the average photons reaching the detector
        p_detected_times_energy_long = np.interp(
            spectra.x,
            MC_energies_keV,
            deposition[0] / (MC_energies_keV / 1000) / 1000000,
        )

        # Use the long normalized fluence
        fluence_norm_long = spectra.y / np.sum(spectra.y)
        nphotons_at_energy = fluence_norm_long * p_detected_times_energy_long

        nphotons_av = np.sum(nphotons_at_energy)
        self.nphoton_av = nphotons_av
        self.nphotons_at = np.array(doses)@w_fluence
        self.ratio = ratio

        if return_dose:
            pp = np.array([0.87810143, 0.01136471])
            return (
                np.array(doses),
                spectra.y,
                ((np.array(doses) / 1000) @ (w_fluence)),
                ((np.array(doses) / 1000) @ (w_fluence)) * pp[0] + pp[1],
            )

        # ----------------------------------------------
        # ----------- Add Noise ------------------------
        # ----------------------------------------------

        if ratio is not None:

            # The noise scales in quadrature while the intensity
            # scales linearly to give the right noise level

            # The real noise level is the ratio of the calculated dose
            # to the reference dose data from MC times nphotons_av which
            # is a factor accounting for the detector efficiency.
            adjusted_ratio = ratio * nphotons_av
            logging.info(f"    Added noise {adjusted_ratio} times reference")
            intensity = (
                intensity * adjusted_ratio
                + noise * (adjusted_ratio) ** (1 / 2)
            ) / adjusted_ratio

[justingm]

Hi Jericho,

Thanks! Probably my stupidity but I am a bit confused so I'll just lay out here how I understood what is happening in the code. When the fastCAT generates a no noise and no scatter image, the intensity is taken from here

https://github.com/jerichooconnell/fastcat/blob/03d32e6aa9fcdff30bb405eae1c18f84bbcaebd3/fastcat/simulate.py#L699-L706

flood_photon_counts is the primary profile (curvature * 660) that was generated for 2e7 initial photons from a Monte Carlo simulation.

Now let's say that I want to use a different geometry and I have also determined my primary field for the 2e7 initial photons. I also know how many initial photons correspond to a certain level of dose. Can I then just directly specify in the "nphoton" parameter the initial photons that correspond to the dose that I want?

Let's say that now I want to use the "mgy" option, I then have to change the pp and pp just contains the fitting coefficients for the dose calibration. It says in this part that you are getting the dose per photon. What does this (np.array(doses) / 1000) @ (w_fluence) give you then? It seems like the "doses" refer to the energy deposited by the beam and w_fluence seems like the normalized x-ray spectra. How do you do the calibration to get pp then?

https://github.com/jerichooconnell/fastcat/blob/03d32e6aa9fcdff30bb405eae1c18f84bbcaebd3/fastcat/simulate.py#L834-L837

[jerichooconnell]

Hi Justin, Thanks for your continued interest! And apologies for the late replies. nphoton should give you back the dose that you calculated. I have a test that checks this for the catphan geometry with a few beams:

https://github.com/jerichooconnell/fastcat/blob/03d32e6aa9fcdff30bb405eae1c18f84bbcaebd3/fastcat/tests/test_fastcat.py#L9-L49

pp is a linear fit that I made using the doses I calculated for 2e7 photons for a number of beams spectra in monte carlo as compared to the dose calculated from the total attenuation in the phantom multiplied by the mass energy absorption coefficient in water or something like that. I think I describe the method in one of the fastcat papers.

So if you calculate the dose in MC for a few of different beams and then get change pp to [1,0] and return the doses. Calculate a linear fit for the two arrays, you should get a new pp.

And then you can change the test referenced above to fit your new code if you like!

[justingm]

I actually did figure it out when I read the paper again. I was just unsure what beams you actually simulated. I initially wanted to just use the option of providing the number of photons that corresponds to the dose I had but that is a bit problematic. For instance, I was comparing 50kV and 90 kV filtered with Al and according to my simulations, the dose for the 90kV is two times higher than 50 kV for the same number of initial particles. That means that I would need less starting particles for 90 kV to achieve my target dose. But when you look at the number of incident photons on the detector at the same dose, 90 kV has more. That is not reflected in the intensity scaling. It appeared that my 90kV scan is noisier than the 50 kV when it should’ve been the opposite.