DVH calculation in ICRP-145 mesh phantom

Dear FLUKA Experts,

I am working on ICRP-145 mesh phantom, I need to calculate the Dose-Volume Histogram (DVH) for several organs, including the Brain, Eye Lens, and Pituitary Gland.

Following previous discussions on the forum ( How to score DVH in mathematical phantom ), I understand the standard workflow is to score dose in a 3D grid using a Cart USRBIN. I have two specific questions :

1- Since ICRP-145 is a mesh phantom, how can I precisely identify which voxel in my Cart USRBIN belongs to which specific organ region?

2- When defining the X/Y/Z boundaries of the USRBIN, should the grid cover the entire volume containing all organs of interest, or should a separate USRBIN be defined for each organ? :

From my .log file, I extracted the following geometric information for some organs (region ID, density, volume, and bounding box):


> Organ	FLUKA_Region	AM_ID	Density(g/cmÂł)	Volume(cmÂł)	Mass(g)	Xmin	Xmax	Ymin	Ymax	Zmin	Zmax
> Brain	83	AM6100	1.041	1457.63	1517.3	-7.117	6.760	-8.742	8.890	73.183	86.548
> Pituitary_gland	158	AM11400	1.031	0.60312	0.622	-0.226	0.867	-2.469	-1.142	76.391	78.026
> Left_eye_components											
> Eye_lens_sensitive(L)	88	AM6600	1.060	0.03669	0.0389	3.048	4.048	-9.049	-8.897	77.186	78.186
> Eye_lens_insensitive(L)	89	AM6601	1.060	0.17856	0.1893	3.048	4.048	-9.049	-8.629	77.186	78.186
> Cornea(L)	90	AM6700	1.087	1.0121	1.100	2.338	4.758	-9.384	-6.845	76.476	78.896
> Aqueous_humour(L)	91	AM6701	1.014	0.30007	0.3043	2.759	4.337	-9.324	-8.897	76.896	78.475
> Vitreous_body(L)	92	AM6702	1.019	5.9379	6.052	2.395	4.701	-8.893	-6.899	76.533	78.839
> Right_eye_components											
> Eye_lens_sensitive(R)	93	AM6800	1.060	0.03669	0.0389	-3.277	-2.277	-9.243	-9.091	77.165	78.165
> Eye_lens_insensitive(R)	94	AM6801	1.060	0.17856	0.1893	-3.277	-2.277	-9.243	-8.823	77.165	78.165
> Cornea\_(R)	95	AM6900	1.087	1.0121	1.100	-3.987	-1.567	-9.578	-7.039	76.455	78.875
> Aqueous_humour(R)	96	AM6901	1.014	0.30007	0.3043	-3.566	-1.988	-9.518	-9.091	76.876	78.454
> Vitreous_body(R)	97	AM6902	1.019	5.9379	6.052	-3.930	-1.624	-9.088	-7.093	76.512	78.818
> Skull											
> Cranium_cortical	50	AM2600	1.904	298.57	568.6	-7.537	7.182	-10.916	9.623	72.426	87.432
> Cranium_spongiosa	51	AM2700	1.252	327.96	410.6	-7.277	6.891	-9.930	9.310	72.585	87.259
> Skull_skin											
> Skin_head_insensitive	166	AM12200	1.088	238.04	259.0	-8.895	9.578	-13.603	10.623	62.329	88.145
> Skin_head_sensitive	167	AM12201	1.088	7.7778	8.46	-8.890	9.577	-13.598	10.618	62.336	88.139

Dear @m.assalmi,

You could use the USRBIN scoring with

  • Type Region to score in each regions (Mesh or standard FLUKA regions) separately, or
  • Type UMesh to score in each tetrahedral element of the mesh. This is useful for the 2D color plots.

Cheers,
David

Dear David,

Thank you for your reply and suggestions.

I tried using USRBIN with Type REGION to score the dose in each organ region of the mesh phantom. However, this approach provides only a single average dose value per region. Since the Dose-Volume Histogram (DVH) requires plotting Volume (%) versus Dose (Gy), it needs a distribution of dose values inside the organ, not just one value.

Therefore, I also tested USRBIN with Type UMESH, which gives the dose scored in each tetrahedral element of the mesh.

The output is a matrix of dose values for the tetrahedral elements, for example in the text output format like:

root@DESKTOP-EAIEEGN:/mnt/c/Manal/S-Value/UterusAF# head -200 AF002_fort.25
 *****  AP geometry, 1 MeV                                                                *****

           DATE:  3/13/26,  TIME: 21:00:17

          Total number of particles followed          10000, for a total weight of  1.0000E+04


   Unstructured Mesh binning n.   3  "U_Wall_Ins" , generalized particle n.  228
      ***** bins corresponding to the element-id:
      from element      0 to element ****** in step of     1
      accurate deposition along the tracks requested
       0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00
       0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00
       0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00
       0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00
       0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00
       0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  4.8576E-09  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00
       0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00
       0.0000E+00  1.0207E-08  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00

This seems suitable for constructing a DVH, since each tetrahedron has its own dose value. My question is:

How can I identify which tetrahedral elements belong to a given organ (Uterus, urinary bladder, etc) in the ICRP-145 mesh phantom? is there a way in FLUKA to obtain the mapping between tetrahedral element IDs and the corresponding organ/region, so that I can group the elements belonging to each organ and compute the DVH?
AF.inp (138.1 KB)
source_newgen_with_regionsource.f (20.2 KB)

Best regards,
M. Assalmi

Hi @m.assalmi

I will look into your problem on the week just before Easter. Sorry for the delay!

Cheers,

Jerzy

Hi Jerzy,

Thank you for your reply. I will wait for your update.

Best regards,
m.assalmi

Hi @m.assalmi ,

First of all, I advise you to process your output data using for example a python script. As it was already mentioned in the forum thread that you referred to, it is not possible to exctract the information that you are asking for simply using native FLUKA tools or Flair.
What can be quite easily extracted from the UMESH type USRBIN is the dose value in each tetrahedron present in the MCRP_AF.ele file. I think that the order of bins should be the same as the order of elements in the MCRP_AF.ele file. If this is true, then it’s easy to group the results by organ just by using the last column of the aforementioned phantom file. I am not sure if this would be enough for you, because I am not sure if all the polygons have the same volume (I assume not). However, the ids of the vertices present in each row of MCRP_AF.ele allow you extract the vertex coordinates from MCRP_AF.node and then it should be possible to calculate the volume of each element of a given organ.
Can you tell me a bit more in detail, what is exactly the data that you need for the Dose-Volume histogram? You wrote

Since the Dose-Volume Histogram (DVH) requires plotting Volume (%) versus Dose (Gy), it needs a distribution of dose values inside the organ, not just one value.

but what is this volume in percentage? Is it a fraction fo the organ volume?

Cheers,
Jerzy

Hi Jerzy,

Thank you for your explanation.

For the Dose-Volume Histogram (DVH), the “volume (%)” corresponds to the fraction of the total organ volume receiving at least a given dose level. In other words, it is the cumulative percentage of the organ volume as a function of dose.

I understand that a Python script is needed, but I still need help on how to identify which tetrahedral elements belong to each organ or region. If I can map each element to its corresponding organ, then I can compute the DVH using volume-weighted dose values by selecting the tetrahedral elements associated with each region and processing their dose contributions.

Thank you again for your help.

Best regards,
M. Assalmi

Hi @m.assalmi

What you are asking is in priciple outside of the scope of this forum since it is more of a topic of data processing. There a few details about the output that you should keep in mind while creating your histogram:

  • The bin order in the output of UMESH type USRBIN is the same as the element order in your .ele file (I verified this)
  • The number that you are getting in a given bin of a scored quantity will be per primary and per cm^3, but the volume of each voxel might be different!

Now about the identification of tetrahedral elements, here are a few hints. First of all, it is useful to transform your USRBIN output into a single vector of rows. If you have a file like this

   Unstructured Mesh binning n.   1  "Dose      " , generalized particle n.  228
      ***** bins corresponding to the element-id:
      from element      0 to element ****** in step of     1
      accurate deposition along the tracks requested
      this is a track-length binning
       0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00
       0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00
       0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00
       0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  4.8981E-09  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00
       0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00
       0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  7.0474E-08  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00
       0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00
       0.0000E+00  8.4074E-09  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00
       0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00
       0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00

first you need to get rid of the error part. For example (if you are using python)

input_file = "UMESH_output.bnn.lis"
output_file = "UMESH_output_NoErrors.lis"
marker = "errors"

with open(input_file) as fin, open(output_file, "w") as fout:
    for line in fin:
        if marker in line:
            break
        fout.write(line) 

then you load the file without the error information to pandas dataframe:

df = pd.read_csv(input_file, skiprows=6, header=None, sep='\s+')

to flatten all the values into a single row.

vec = df.values.flatten()

You might need to drop a few last elements if they are NaNs (this can happen if the number of voxels does not match exactly the shape of the output matrix shape)
then you can load your label file with the phantom elements

phantom = pd.read_csv("MRCP_AF.ele", header=None, sep="\s+", skiprows=1, names=['index', 'v1', 'v2', 'v3', 'v4', 'organ'])

then you can add your vector with dose values

phantom.insert(loc=len(phantom.columns), column="dose", value=vec)

and finally you can create for example groups by organ

groups = {
    k: g.copy()
    for k, g in phantom.groupby("organ")
}

Now if you define a dictionary with your organs, like for example here:

organs = {
    "Brain" : 6100,
    "Pituitary_gland" : 11400,
    "Aqueous_humour" : 6901,
    "Cranium_spongiosa" : 2700,
    "Scapulae_cortical" : 4500,
}

you can get the voxels corresponding only to the Brain for example by taking

brain = groups[organs["Brain"]]

and you should end up with something like this

brain.head
<bound method NDFrame.head of            index      v1      v2      v3      v4  organ  dose
982          982  613780  613256  613249  612673   6100   0.0
1306        1306  615388  615079  611752  611751   6100   0.0
2669        2669  612920  611487  611170  613321   6100   0.0
4660        4660  614729  614718  614725  614726   6100   0.0
4682        4682  611554  612274  612273  613224   6100   0.0
...          ...     ...     ...     ...     ...    ...   ...
8572690  8572690  610590  611953  612938  610261   6100   0.0
8574033  8574033  610962  610472  614065  611014   6100   0.0
8574849  8574849  610290  611610  612100  613149   6100   0.0
8578402  8578402  610463  609945  610227  610064   6100   0.0
8578404  8578404  613030  610463  610064  612413   6100   0.0

[19398 rows x 7 columns]>

Now in the columns v1, v2, v3, v4 you have the indices of the vertices that you can find in the .node file. From these you can calculate the volume of each bin and add it to the dataframe.
The rest I will leave to you to figure out.

Cheers,
Jerzy

2 Likes