From 4d69d07fc933675ce247200969d1ff7d7a48a6f3 Mon Sep 17 00:00:00 2001 From: Southerby Date: Fri, 14 Feb 2025 14:09:32 +0000 Subject: [PATCH] Changed Gamma Analysis to include local Gamma calc. Uncertainty Stats are now also printed, and saved to uncertainty text file --- gate-pbt/analysis/analysis.py | 33 ++++- gate-pbt/analysis/gamma.py | 4 +- gate-pbt/analysis/gamma_index.py | 53 +++++-- gate-pbt/analysis/mergeresults.py | 2 +- gate-pbt/analysis/mhdtodicom.py | 6 +- gate-pbt/analysis/resampled_dosemask.mhd | 13 ++ gate-pbt/analysis/resampled_dosemask.raw | Bin 0 -> 495550 bytes gate-pbt/analysis/uncertainty_stats.py | 180 ++++++++++++++++++++++- 8 files changed, 262 insertions(+), 29 deletions(-) create mode 100644 gate-pbt/analysis/resampled_dosemask.mhd create mode 100644 gate-pbt/analysis/resampled_dosemask.raw diff --git a/gate-pbt/analysis/analysis.py b/gate-pbt/analysis/analysis.py index f916d89..7c9a4a1 100644 --- a/gate-pbt/analysis/analysis.py +++ b/gate-pbt/analysis/analysis.py @@ -25,6 +25,7 @@ import dicomtomhd import gamma import overrides +from uncertainty_stats import getUncertaintyPlots @@ -153,9 +154,10 @@ def full_analysis( outputdir ): print("Fields found: ", fieldnames) - for field in fieldnames: - + + + print("\nAnalyzing field: ", field) print(" Merging results...") @@ -181,7 +183,7 @@ def full_analysis( outputdir ): path_to_dcmdose = mhdtodicom.get_dcm_file_path( outputdir, beamref ) - tps_dose = dicomtomhd.dcm2mhd( path_to_dcmdose ) + tps_dose = dicomtomhd.dcm2mhd( path_to_dcmdose ) print(" Overriding TPS outside of zSurface to zero for gamma analysis") struct_file = mhdtodicom.get_struct_file_path( outputdir ) tps_dose = overrides.set_external_dose_zero( tps_dose, struct_file, "zSurface" ) ## hard-coded @@ -224,7 +226,8 @@ def full_analysis( outputdir ): #itk.imwrite(dose_none_ext, path_to_none_ext) #mhdtodicom.mhd2dcm(dose_none_ext, path_to_dcmdose, join(outputdir, field+"_DoseToWater_NoneExt.dcm") ) - + #''' + print(" Performing gamma analysis 3%/3mm: post-sim D2W vs Eclipse") gamma_img = gamma.gamma_image( d2wimg, tps_dose, 3, 3 ) itk.imwrite(gamma_img, join(outputdir, field+"_Gamma_33.mhd") ) @@ -242,7 +245,8 @@ def full_analysis( outputdir ): itk.imwrite(gamma_img_22, join(outputdir, field+"_Gamma_22.mhd") ) pass_rate = gamma.get_pass_rate( gamma_img_22 ) print(" *** Gamma pass rate @ 2%/2mm = {}%".format( round(pass_rate,2) )) - + + #''' dose2water = field+"_merged-DoseToWater.mhd" @@ -250,6 +254,7 @@ def full_analysis( outputdir ): print("\n") print(" Scaling merged-DoseToWater.mhd") doseimg_path = join(outputdir, dose2water) + # Read DoseMask to set dose outside zSurface to zero dosemask = itk.imread(join(parentdir,"data","DoseMask.mhd")) # TODO; read from config file @@ -272,6 +277,24 @@ def full_analysis( outputdir ): print(" *** Gamma GD2W pass rate @ 3%/3mm = {}%".format( round(pass_rate,2) )) print(" Converting gamma image to dicom") gamma_dcm = join(outputdir, field+"_Gamma_GD2W_33.dcm") + + + # Maybe we add uncertainy metrics here? + + dose = field+"_merged-Uncertainty.mhd" + if dose in [basename(f) for f in mergedfiles]: + print("\n") + print("Calculating uncertainty") + doseUncert_path = join(outputdir, dose) + dose2W = d2wimg + print("Number of Primaries Simulated", nsim) + results = getUncertaintyPlots(doseUncert_path, dose2W, nsim) + + + # Dose To Water MHD file + + + diff --git a/gate-pbt/analysis/gamma.py b/gate-pbt/analysis/gamma.py index 89f30cf..ba0a86b 100644 --- a/gate-pbt/analysis/gamma.py +++ b/gate-pbt/analysis/gamma.py @@ -20,7 +20,7 @@ import itk #from gatetools import gamma_index as gi -import gamma_index +import gamma_index as gamma_index import reorientate @@ -39,7 +39,7 @@ def gamma_image( ref_dose, target_dose, dta_val, dd_val ): Accepts ITK-like image, or path to image Returns image matching dimensions of target Set ref -> MC dose, target -> TPS dose - """ + """ ref, targ = None, None if type(ref_dose)==str: #Assume we have file path diff --git a/gate-pbt/analysis/gamma_index.py b/gate-pbt/analysis/gamma_index.py index 6c79327..8fd5925 100644 --- a/gate-pbt/analysis/gamma_index.py +++ b/gate-pbt/analysis/gamma_index.py @@ -17,7 +17,6 @@ import numpy as np import itk from scipy.spatial import cKDTree -from scipy.interpolate import RegularGridInterpolator def get_gamma_index(ref,target,**kwargs): """ @@ -99,7 +98,7 @@ def GetGamma(d0, d1, x0, x1, y0, y1, z0, z1, Max, dd, dta, gamma_method): else: print('Please specify gamma method correctly, local or global.') return False - + return np.sqrt( (d1 - d0) ** 2 / (0.01 * dd * norm_val) ** 2 + (x1 - x0) ** 2 / dta ** 2 + @@ -181,7 +180,7 @@ def gamma_index_3d(imgref, imgtarget, dta=3., dd=3., ddpercent=True, threshold=0 - gamma, total = 0,0 # for gamma calc in testing + gamma, gammaLocal, total = 0, 0, 0 # for gamma calc in testing # Loop over each target position for Xt in x_tgt: @@ -214,11 +213,13 @@ def gamma_index_3d(imgref, imgtarget, dta=3., dd=3., ddpercent=True, threshold=0 found = False # Flag to indicate when to exit all loops + foundLoc = False # iterate over 4 closest references voxels for idx in indices: - if found: - break # Exit the outer loop if found + if foundLoc: + if found: + break # Exit the outer loop if found closest_ref_point = ref_coords[idx] # (x, y, z) coordinates of closest reference point X, Y, Z = closest_ref_point[0], closest_ref_point[1], closest_ref_point[2] @@ -235,13 +236,21 @@ def gamma_index_3d(imgref, imgtarget, dta=3., dd=3., ddpercent=True, threshold=0 # strictly speaking, this is not a first filter, as we only loop over 4 voxels # in reality, we would need to loop over all voxels within dta and check for filter 1 Gamma = GetGamma(tgt_value, ref_value, Xt, X, Yt, Y, Zt, Z, max_targ, dd, dta, gamma_method) + GammaLocal = GetGamma(tgt_value, ref_value, Xt, X, Yt, Y, Zt, Z, max_targ, dd, dta, 'local') + if GammaLocal <=1: + gammaLocal +=1 + foundLoc = True + if Gamma <= 1: gamma += 1 - found = True # we have passed this target voxel + found = True # update gamma array values for each method gamma_array[tuple(tgt_index)] = Gamma - break # Exit the outer loop + #break # Exit the outer loop + if foundLoc: + if found: + break dirns = [+1, -1] @@ -251,6 +260,7 @@ def gamma_index_3d(imgref, imgtarget, dta=3., dd=3., ddpercent=True, threshold=0 # loop over both directions, i.e. +X or -X for dirn in dirns: x,y,z = X,Y,Z + xLoc,yLoc,zLoc = X,Y,Z # interpolating in X if D == 0: # check we are not at a boundary @@ -263,6 +273,7 @@ def gamma_index_3d(imgref, imgtarget, dta=3., dd=3., ddpercent=True, threshold=0 #calculate the optimum x position, and the corresponding dose value x, interpolated_value = min_gamma(X, X + dirn*referenceArraySpacing[D], ref_value, ref_value_next, Xt, tgt_value, dd, dta, max_targ, referenceArraySpacing[D], gamma_method) + xLoc, interpolated_valueLoc = min_gamma(X, X + dirn*referenceArraySpacing[D], ref_value, ref_value_next, Xt, tgt_value, dd, dta, max_targ, referenceArraySpacing[D], 'local') # simple check to see if x is within bounds of approximation if x == False: @@ -275,6 +286,7 @@ def gamma_index_3d(imgref, imgtarget, dta=3., dd=3., ddpercent=True, threshold=0 continue ref_value_next = referenceArray[i_ref, int(j_ref + dirn), k_ref] y, interpolated_value = min_gamma(Y, Y + dirn*referenceArraySpacing[D], ref_value, ref_value_next, Yt, tgt_value, dd, dta, max_targ, referenceArraySpacing[D], gamma_method) + yLoc, interpolated_valueLoc = min_gamma(Y, Y + dirn*referenceArraySpacing[D], ref_value, ref_value_next, Yt, tgt_value, dd, dta, max_targ, referenceArraySpacing[D], 'local') if y == False: continue # interpolating in Z @@ -287,25 +299,36 @@ def gamma_index_3d(imgref, imgtarget, dta=3., dd=3., ddpercent=True, threshold=0 if ref_value == ref_value_next: continue z, interpolated_value = min_gamma(Z, Z + dirn*referenceArraySpacing[D], ref_value, ref_value_next, Zt, tgt_value, dd, dta, max_targ, referenceArraySpacing[D], gamma_method) + zLoc, interpolated_valueLoc = min_gamma(Z, Z + dirn*referenceArraySpacing[D], ref_value, ref_value_next, Zt, tgt_value, dd, dta, max_targ, referenceArraySpacing[D], 'local') if z == False: continue # simple gamma calc function Gamma = GetGamma(tgt_value, interpolated_value, Xt, x, Yt, y, Zt, z, max_targ, dd, dta, gamma_method) - - # do they pass? + GammaLocal = GetGamma(tgt_value, interpolated_valueLoc, Xt, xLoc, Yt, yLoc, Zt, zLoc, max_targ, dd, dta, 'local') + + + # do they pass? + if GammaLocal <=1: + gammaLocal +=1 + foundLoc =True if Gamma <= 1: gamma += 1 found = True gamma_array[tuple(tgt_index)] = Gamma - break # Exit the inner loop if condition is met + if foundLoc: + if found: + break + if foundLoc: + if found: + break # Exit the outer loop if inner loop condition was met + if foundLoc: if found: - break # Exit the outer loop if inner loop condition was met - if found: - break + break - print('Linear optimum linear interpolation method = ',100.0 - 100*gamma_array[gamma_array>1].size / gamma_array[gamma_array>0].size, '%') + print('Global Gamma - linear interpolation method = ',100.0 - 100*gamma_array[gamma_array>1].size / gamma_array[gamma_array>0].size, '%') + print('Local Gamma - linear interpolation method = ', gammaLocal/total * 100, '%') # convert to image, note we copy information from targetimage gimg=itk.image_from_array(gamma_array.swapaxes(0,2).astype(np.float32).copy()) @@ -356,7 +379,7 @@ def min_gamma(x1, x2, y1, y2, tx, ty, dd, dta, Max, spacing, gamma_method): # optimal x position for local or global gamma if gamma_method == 'local': - x_opto = (grad*ty/(dd*0.01*ty)**2 + tx*1/(dta)**2 - grad*c/(dd*0.01*ty)**2)/(grad**2/(dd*0.01*ty)**2 + 1/dta**2) + x_opto= (grad*ty/(dd*0.01*ty)**2 + tx*1/(dta)**2 - grad*c/(dd*0.01*ty)**2)/(grad**2/(dd*0.01*ty)**2 + 1/dta**2) elif gamma_method == 'global': x_opto = (grad*ty/(dd*0.01*Max)**2 + tx*1/(dta)**2 - grad*c/(dd*0.01*Max)**2)/(grad**2/(dd*0.01*Max)**2 + 1/dta**2) else: diff --git a/gate-pbt/analysis/mergeresults.py b/gate-pbt/analysis/mergeresults.py index c173d67..7293c5f 100644 --- a/gate-pbt/analysis/mergeresults.py +++ b/gate-pbt/analysis/mergeresults.py @@ -84,7 +84,7 @@ def combine_uncertainty(dosefiles, dosesquaredfiles, statfiles, output): #treatment planning" and dividing by dose/N for relative uncertainty... if sumdosesq[i]!=0 and sumdose[i]!=0 and N>1: uncertainty[i] = math.sqrt( 1.0/(N-1) * (sumdosesq[i]/N - (sumdose[i]/N)**2) ) / (sumdose[i]/N) - + combinedUncert = uncertainty.reshape( shape ) uncertimg = itk.image_from_array( combinedUncert.astype(np.float32) ) ## ITK CANNOT WRITE DOUBLES, MUST CAST TO FLOAT uncertimg.CopyInformation( img1 ) diff --git a/gate-pbt/analysis/mhdtodicom.py b/gate-pbt/analysis/mhdtodicom.py index 050a1a8..9d57ab9 100644 --- a/gate-pbt/analysis/mhdtodicom.py +++ b/gate-pbt/analysis/mhdtodicom.py @@ -109,6 +109,7 @@ def mhd2dcm(mhdFile, dcmFile, output, dosescaling=None): print(" Dose scaling of {} used in mhd2dcm".format(dosescaling)) dcm = pydicom.dcmread(dcmFile) + mhd=None if type(mhdFile)==str: # Assume file path @@ -159,7 +160,7 @@ def mhd2dcm(mhdFile, dcmFile, output, dosescaling=None): # TODO: check this is safe # TODO: should this ever negative? - NO! dcm.GridFrameOffsetVector = [ x*mhd.GetSpacing()[2] for x in range(mhdpix.shape[0]) ] - + dose_abs = mhdpix * dosescaling scale_to_int = 1E4 # Dicom stores array of integers only @@ -174,5 +175,4 @@ def mhd2dcm(mhdFile, dcmFile, output, dosescaling=None): #dosescaling = dcm.DoseGridScaling #FrameIncrementPointer ? - dcm.save_as( output ) - + dcm.save_as( output) diff --git a/gate-pbt/analysis/resampled_dosemask.mhd b/gate-pbt/analysis/resampled_dosemask.mhd new file mode 100644 index 0000000..83a40d5 --- /dev/null +++ b/gate-pbt/analysis/resampled_dosemask.mhd @@ -0,0 +1,13 @@ +ObjectType = Image +NDims = 3 +BinaryData = True +BinaryDataByteOrderMSB = False +CompressedData = False +TransformMatrix = 1 0 0 0 1 0 0 0 1 +Offset = -102.29000000000001 33.5655 -44 +CenterOfRotation = 0 0 0 +AnatomicalOrientation = RAI +ElementSpacing = 1.96234 1.97088 2 +DimSize = 106 55 85 +ElementType = MET_UCHAR +ElementDataFile = resampled_dosemask.raw diff --git a/gate-pbt/analysis/resampled_dosemask.raw b/gate-pbt/analysis/resampled_dosemask.raw new file mode 100644 index 0000000000000000000000000000000000000000..9e6d9cd3ff9422519f819be8649f007338a84cf0 GIT binary patch literal 495550 zcmeI5?PBXZj)Q0K|EAaJX=3Y(1VBcwn{fWrqDVXjvQvHJ@*HoIPB88&7JZ96&ka#mZmA@35=Fs9II=u{ux5HEYv)DF=R;TFvDkRaDM2@yp;^{kn4*kT0)# z>uXv599)YxJLiCW)#xo3vixOmDc){eCe+JTaJf?DFLO)nZsRhdU3PlLQqw=D4C!v` zoYAg2Icuo_o>P``zxSNce)ekST0{JsGKGiDeWd^CsNt@?f{gNnhkcba2>AsO~Wz@qK14 z{qzn>5}gheX&MVytn|R3=56Nb)NB27g=6e4H_(v!b$B1ieOj)!6}$FK#ni9U zMY^WJURQc-V%Pp{>3JaQ4exrx{&=}_7MC|BcP4pWuPESDdz-|;-*kg+|?aP!TfMb(fAlr$P@+I zND81hOGmf-@^@pBYKpcb$0v|9Qwm6PwT+~K1XC)AN9iQ;sR)2LS0NCC3W9J{7<33g zFSZ#Vgq6g3$e6woIX_Y>ST^UM-avZ~=A!x_UY7U19O6+6u1!^VobyoQgjy+*n%rZq z(yR2E1iS~rejjS2?1c+*-tuhq>n=b_a)ppmP!K35zd<~hB)vtpH>DsSf z--@K$xK7^6Qp!t&`_}$r`I>=zjf9u-VuXP*P}@i-C~oSGq~Mc|TRJ-B$NMq3jYNS$ zrnZqtP{7oUN$TTATax1wNSY}HB)QT+iYbxAqf`)`pAKSC5fT+(2*kJwf{t&Pgswo~Mx>7)zDT$<@G?0>?N}^rKAeN+qj>+MT+bD!U3@Ql1QDG38UBN+DtI$i(l=x)n0D9%ptC6r*SK`Y20VZckENUtb_i%L%j?KM*Ay`tbRyzsoY z2bk!BvPeECi3EUhNC?P%AIA;~LSa%)*FF;TEw*#F;5K3kl=9jDlqZ4i`zTb&>noC% zkg@NFm-HgB?-vTOlw%4bWxSqmBY~iV{63HbTMEawjYNQgrnZqNP{`Ca5(x^JqDktj z9}n6{(olk>W5js#UT$e50VRTXS2Bn-rGt*i<$=sLLL>$i1L3ZsoX7kusG*htLRd+h zADO|#rI5+7VCkG+!Hbp0rH<)@0!!~60@B@Rk-9*rWu-^O)iqaoT_DWT3lDoQ0*Eds zi{yinNB}5@gn)8L*t-^R+CU593XlO>_K}eMfE$Kdd{>YR(6o=_=Lg&_)Z%)ARFEc@ zl=Q3XX9xaa;hXX;C7}>d4hkb>9Mv`w2ukpTlH#asBp4Ka{BB8#ofKgy7#~H7nBqVi zNji$O6gVlGq?+PEsw;&gxspH%NdqbQsU+H!3}RgcKpZLrI)?qvcE24_5R8h0Fk2~{ z9}zplJdO=Z=KPZzt;v{CO09Uepz4luMH(}~E0rFe)!13-jhVm|FFgEx&qID7$+eW0 z-knHUrhHFXB;X4^@99PNeN(8C*SC-4$LHNEyziNUh1|w%Bo73C-Va+!Bl(uHPzWjI zRRRKSBxRs5PYEOtw2>4?p`OAd7!>}Ux+N)=A7LpNA4Q6o;y@coI*PP(w7c(hi-v6^ z9;BL5NK$?hNO2{C1e8o-1JpsJs|bif#X!gO;y1J}hdB>{UUarj0RgNm&X3N_$Zsu` zW57g?ap&ZgcLSuH(Ox7Y9CxeIb0abuDm^0{cdOEKBU(3j;X!Y&Yx&ELOG-l(=@L7C z>Bh%PE++0nTHK&e_B;uai;??~KJqs=$&`0WKl_Ul=E~bh5GW-+jFf@eMgl?D-?iog z@_{7SQaHX_l41`!VJR3NMT(f>KoM8lNGwP<#goLN5Z6pb#m=ZxeN(Qm00O%Ms z-|^bb7>GbcL9nd^5X8#i{HU}?(ybCXCQRdWHyc=b{D5>fTBI%zYFX(~adpj=UKa?p ztn{e3y5{?&$9zKROms`nQ<1vBnx&Us(;Qkjq>y$y{Dq{vhkW0^WOt6~*lv$iU-_gg z)HV_TO3H5=2?C{<`alwBDIvZalA^!8u%&Q(1Sw{U0c|9OP?V<#i33F*3YH|hqCv8& zZ6pCCnvzJ$QE4E-l?tL!I*3F?K*ucp;W7wJwb6bvW?_`p!e~WiMIZQq&9%q%8oCl-A>pN;;_zBte!^;`=}nXelAS8PkJ|tE5k(@7=rMsUG?epZ2HvU4Tp*%~_2jlsre}W-Yv;l4TKTbBcE3f^KJ+hmCg^XHY)jyW z$VFUj*(0Up=Yl?xlx50KFNp+zHj;8sh^IeBf-I%HPTY`U2t8FA2B%2aQ+EFPW&6P$HT!|pwl?-B00njnpe%dxyF%W@@f>2ut z91O*S%!42{T{59%w|*kY{^A2Wdpot)t{2~{F?>;^NQ9r5p7WST;_7GZ`1vh`M+eRWm0aG+d z&5sAEt`w5wN&+daL=caXK`bf&;!q*bG3>k|oZC?ljEaLWTRi6sBX;GJWBTXxjY#KB z1+QFBLH*0Phq|-|FN+ij@DtNBzet2%Orccm_#g?e=_VxbZm~Ag6DEOnStwKk>oZ80 z0NiDe@DC>t?E^`crJb)4i6>=S$~1K!Qd)eTrM&p=M9PZK_w?sTfTg5o;%iA6&jb`+ z%0YpoxG5C0krYP3o}wfi6n!WdQpgnr+DHnZI7>%2^f7hTM3Yn$4^mwzB*~QoQb-!; zn3|oqS(8aDN(YgyA|M78<2AZ7kG9XGU% zWP!lXCSZPdAtj!aZ7DOp`;hYD^E{=IT+m07vP}8u{W%g~De0H@W~7Xhf-I%P2a*y@ zZ6l$ecz$;z1y2gM6pfD|g-lVPkSmg;yJA7Qt8FA6q?%Gl(osnu#gz!+Q8MTl1pkEj z*A@X$s2GSqML}qSG>1ADquB==!_?M|NYs5>&8{bFnA+P%GQVMY+y-^uX0vC?9Huru zk>ni|6_DL*$7j!SyZcCXeB50--KHsLmfOAuiMwa7+w$bia{C`hvMeROk{>ynq)e0z z+DJ-7Ii5a}lxNCIFOB4aHj=VXzNe2Q<(LA(OZvckGg3x;kfoIPKvIIKZ6p*F&u<$E z2SrUWq|ie_fkLjfkw{R$6irf(iU+B#6q4jh0x7OU5RZ~UEGhuvP$AGU4F7Q(1;Hp5 z1UqmVeXb(_Z+|5TPRDC>)AI!IeGi2zd7VY_5;FJg_{>pqb05k6GK{=!yWKQpkCMCl zNM?NGT|4cTDQA@1wm|qhIl|mVGC`@HK9FQvf(hNiRq@BW1(~SxSiyBqf;IMnXaH{O(8!o)m8B=v=4h%>k~=#OD-LeHDJbP4bpo>q9Z#8%<#}R@>9;ZDv}Y7 zw`GS1f8z2m$9|E3kd`g?NiXkr*YorTNv@@|->5Go2-9v@B$nj%0ONg)(t>F8a)Yx^4{Nk_3D-4#y~T?rtKq=2;iG?L&-1RZ1kGu?FP zAQBY;QK%S*K#`y(_Q1l&2E?1uuOu}IL>6v(5%G4`L$%$J)fbW?arLuye0_IhU8G8{ zO`x8LHV2Z0r-P_W$LBxl<_WmcUB27YpOZQ2S{ex+3(eso>azSa&`InM%P6`Mw z>7XD|im7cR43yysBqgA>kx)?l@xw{c{I-z@P%uA=6v>YRMOr&Wa<-1^)JvS2aT7U zM*AMvv=cE8jN%7r`loB80DT9Q8GEVA#q;P(yrFeWeDQJoS1zl03h${}X zk))$YOGmr=>1WZvlSGsN(p*U(1to$6luTmJN(YgyBA{b-@i0%y9B3k-lmMtAJDc%Q zvni$Kk`&SD%qNoSTkPpkqebVDKxL3V8-*yKJw-y3n)OwBvsBPprH3Xp>#Ou;si3v* zOON=#>V;^B$9y#YJEXk$y!8Gc$+eXBOx%i;bW(t)91;TBNJ>FLo-#-nDC0rDA1Rz4 zYAGHcPKuf$NHG)x+DM9^C`(7P_c>v`&yr*m4N_eRAPpsfB$PY(leTUlrcYW_shH-3w4oiVxRJQAdu*>0oEQF60LWMEBurPstk zE-F1Tu%^AzYvLdm-h zf1qHa0>Z;+k%mZ!!cEd6?$~cNJ^ABAN0ErIhE_8@<)h?4covHEOhWnjed#4We0F%$ zr*P-&JW`&i42DbbP$P(fG89Paz1UN+W-`m^7>qqCjjpk z365&k=*2gi(P|D!6P-?dBFXtj^V1x}q8*PtyUCk3RJLqb3~&&1tG z3Hf1`GU7u?aZ}q!FesegElIJHA}k%Pi+56}%#i{p5~RDLK{84pX#pxA$(05YP%4N= z=_K;30wB&61%h4#{3f62Yk(at3e@r2isWZxba?R@qvV!DBGs(+PbAg1*w%@=d*&nI zPBRTfV$$k+E4@AuX8m>Pk)Q7DEY0+aWQXUxKPARz9~1W}-1c-IQeJ$Xr!hDrq`%F=I9CK0Qh*(!n{}V~(V{WEFZ(8xuV&=dHk}OM! z?@mXoCMgwFB=aeDHtqB>4{9690VR6!NO`Cmk+Mv=mbkynjOYATq@?%&PdOw6lyg#d zBW2`=c}gIGpoEjUAt{<4Y$+TcL5i7TKrvU_NE|3)iY3WMMT1mV0!TzjBqcx%B)C#R z$DHw@WU>_iaVQEjxlyU)LM7GNP?2V-z%?&ESjcU3NbWq|P$WF7Nv9p(WI!%`9Z61V zHfW|d8<1r3jlJZcV zr;Vg6lxvCm%glH^x)muYKETsPQVI&Ol=DnrqzqHrNFXR7zZ;UGCk0yy$48KYrWjHP zMS(Vw0w~VX(QUl>m|HYRMhPU%l>(AnX&?cmf@qWuIwsBMtlMNj-Kf5Fp}KT8S)?uy zZt2Cxsaf3)$r>g1ibMw1w%hTwy|L9zNHG)LwCSnqjV+55$!}kJ^?iBdoc5)+?};!; z)Q3AS+3zD|#%HIONOC|MNqH#G(?(Jj%C*G(WoG>H-inkIAK+;tDFuaC%6TR*QiiE* zBoLI4-wjF8lY%XU<0D8xQw%ADqCgu-0TgHH=r-Pb%qWlvq)}4R)Za%HOcPXgv8#q*>8F>C)v#+ zv58&#E4^zP?Df~Br+q$N#A>Ejq(y4avYlRLeA`G4DAAKg%0qo1$+blOY4&`g$Dg+% zB|Q^>r;Vf(6k;jonZQUHrfx@yn*uE*oD@t7o7zUgLDBpeQYb$Pw2>4*ah8tm@jVBt z&>$HlkTh2cNOGlt1e6M*Q96i3F(B+g2D)4*pnR(fDi(>60b?UVSFQ~zDxi&((2eNKr7 z<$yMl@=%_qKS**dk$;-4AH!Rbl1>Wnw2_p8LM-Jx6BsGO)a^)dQ=p}UlY&WMQ`<;5 zD4O3k5(A2v;z)sqf&}TVXpoE&NLqjjNOGlt1e6M*Q96i3F`)IwDsIPTa(YI$i!x@p ztxqJG|Nq(#!(2pGw-cW=O70bj46JRp<7<0ktFI%);_9Z&^t#^IvPhBqD!odtO~8C2 z65n#a@muW7`0R5^Bsrjsq&$@8=_5&5rrhw-KA8igBvbbyrIJbff7%aj`) z`A1;?x^S)mQj)2Aky1<{mU2!CBW0Mn9Vu=Kw3KjCFez+m8wm$R^V>#ZKrvGsDezE` zAl(%Wl2HOl3s3<`t~8Kf?`@VdGW30 z__{;l1(y>=T1JCx*zq9>XwTP?&|siGGd)NFZ7LF))U2=4o27!*%=D;_TpqDVUP9Ks zonGcaZ6i6LL{Ape7*Z%d3bc_FKyj9iZsRS%+@e7;N+4;j6p-Xf0|_VeIzN%lp7xRH=yxr43Lsc-HViB z3bB-PQWz=2)a^)dQ=p}UlY&WMQ`<;5D4HKb3gt(EA|wtJ$&V$;u4s_zN&tx{iKGOm zfdp47=$JD;DonNlAPz-=CO4{-e5Rs0o7qQlz6m>SWOjd2%{3*d&W1jb#32eIbF zC#ufIiZn|Ft?f(CdHB#9M*2)shydJEBsi)`T>Li_c5%BS~4N-0;$#34oMj>RzN2Q;4OUlfp1e%v%3-g!kvLGq6ibrxqd}@G0VJX%k`kZ>5?raEV~+b>3V=8i1>!s`y8-5M zmfP@2oSxC`qKsK?t4Lg4UAGrs*B)3pq=0fiTBJS^X1y;x=L-quroCR)H;F`xJiF9k1_0%8SqQ^pT`2Q*L-^Pai-^GIcLf ziYdfW&PidU3{$rw#Z7^h5>5&xMNQ$LjieZguynMJZ=!gGBLz?-Xd_8Rv6kec5l-PAm5n6WZCz8B_ zvJ;_po%kp%ySYeqKKAY>rsq87b|BoY9UrG!1Kq66#*06jO+$oM!?fWtiGV0znD+-H;SLDcI7{XuYEO zd)r71C}xTy1@a?78%Z*XwRH6R&)BR(0BI-*B%wqS&yqnTD!_S+&W?Rv<3Mu_x(YyP z&H2EcNb@@ss;((Pb3Rg}c`jIWUwYnO)Y$z(U5+V818XW0meaiJ6Vr1Z^LY?j)s9co zoKF>Lo(op3^st=fUE9*LLB!|PxsOZLAtj=|k(3vom);*Fxt5qe%r1}gczG*Q(n$fH za!3d$=cMjN%E%A%lt2PO2`2@U!mhTFa8NWqh7`(=0&OG(P@JWsdweIrDl|w&Z6gUF z&6ET>hOm>lbg3X3rGrRR1Vkmsfi9T`E)ixx{yGYvuI+*APb4|t8XcAcfwiM0kAQSE(&R*tEVb~{|8HV2Z0r-IZRO2e(vqySGjBm|UmQgHB= zpm2T!DV84tin-cG;y@8oEJ;2p8nlrlqIgTkg!8sK-IhQKN(Aw?bP#DP0-|gMK{$&A zO`j;Th%rggc{_xfkCN4Uk;t20wF9PBqh$RPN#2*Emmih^L}+B(%o#GQx^!etn{e3x@IRm>CcM}@w+V2JQu9$rI&Y5+ej`b?bFehl9HYY zAiW$C0?K(N?nX+;53_W1VXrqv*ESLi3Y)qmDRxqXrK8n2Zg*=Ni30^pu^=7Alf8uMW_Au zB290RF9N1`(don|lDse7z*mE5UR}|Ruj|*BMbfi7CcXGV;?3}FNX=cv7LKQwcr#t3 zCJu6OoAj{HTaY`)Q<2OR?~`8ML3y6iNG@n2DGTL$+DJ-50iJS52q@>%<=sdL`C*na z;zLPsQ`<-|D4ZWbisi?EHj*MJ%F@wnzIB_MEJ;Supp7IE#aj}ON&+cQg$UwNGKftu z0-{(!5X_2m9t!4b!oN)c1gLr~Ak~ur_U}dlZrN{vO!+cE$Gu3v4I?dpDMSY7DH0IU zy2XueJxcBri43f1-58=!0Hx)=k?;1Sn{VB1Q7!KpROq zinMgJuOBUYhbM_B0i?N-K$0nqB)C#RG)e~@llFtw?1ngx;YQ6^gh3!HgY(cZto+-O zK(4CSfKojje}5mz{eI>5R{VQ*J1|qO4AKCh5FxkcTS;y(SOX9rCggT~A_@5Ab#nn5 z-1tBpzpY4qM(d7!>A4Tz`X<>O&lIb-?l00h8-M2;)AK*9`nmmdc;Nvbi;HB12&Ut; zS%}{~PXgQrNI6d*0@_GQK|!7}NEj&NxA%UeaDJ$zczif1Xo>&@U2P*#ppYq&q#qRv z(p_yM@gUKZ0@6?#Nq8zm5bsI`9fRPvpyah+Bu@p}hft7`*Yt@b+4sw*C}k))vnOUKCijc&Q4fdrHaqER}C zOfUpuSWytnN&rEu9L_^!`CR{OOl!P`IE+n{aB2a4)Xa^@ja#jalL9D9)=om$Q zg?@WtAeoC})U2lak*vD~_drvin$`A+BNTJ;W5QW_1E2 zYn<3C5)szW`labXpNxmFh60dYk%rWs1th(k_&$&XS;GDTe4fYR=j}-GX9Bd85Fbnm zo7zUgLDBrSkr+_S6h{g?6eMUPNk*}jj{e06|$%ifQhdTX1MjyY1dbg1*Rk z-A;R#-DXn1$iEkK^3PqA)`yKk$^vn<7aG*$fe~biK zvj2Vo3i!Q# z3ukbaK$+Y-SYqLz?5vGM%_P7YM*>BsZTm>j*DJ5vncgf~e5n8zsOF*%AUsHjyCEqYAM7bg!a>oK+D2kPF;g5Va#Y($ zBuLMX2Fda9BoQTm#QY?Z;z|Pvu2c|>(m^CD0-{hc&@t`2zrMD^pohRqZGDl#0a$6A zAE6C?&KsV=&fVmz7M3hNlLNm8n2mO|_SW%LCj%^T?XDBhZda5+=KR0KiZ0%ID zu2|Jszomcf%_iB}sb*cVs`i95Qln-qIN3?JCEo5$hN)7kX00gMvDuRD_QqkV(y5tCN_c2C zrMt~xm??Bx&XQ5?dmZU+yPsy(os_p^q}yIkdDw5q@zo&Y7I5vh(-j^YS^;!5>Hrl~ zyy*jUCoPQxwwS0OB~1L%bM)oe*n1 zPVKlXO3i6DV9unlJ%3IUtw z8(CBRJh3KGgmT}ykHno+f1lKvGDRu(&HG4HeEn^5Yr+&KJaq3P5%CRoNfvNZr0~!V zlZHvQu#W^E6l-DZM@@~BY^6wSV%Pp}P7nE~IM=MkYc#zgjoCr}!lN2(=Sc;t#M+%rCK^sXjinXMl6i*UO2_VtcHj)HVOo<>KC4*R00CbF6 z?+)Ew3`C%!AlOy{2x8@MepI?MX_tu{3zm)f)oZ;lpSP47deW5HZ<%>!ME#V)p$O0B zdy}4bN&XZ}k%VXaJ`(&z&g+(*_elSIDWDfL;WdB~o_()D`r@eo{-pKK0g&))I|uoz zrv~M|x=5+p&r8NpG_n6N045ZY^`(tAZATzq$W z%pbS2aPjm$k`y#Wq!%MGppB#uit-d8aUkMLoZ%Phz9cz6){-6{PZCWDAkCEol1ynN z0VRU?{A3V|3V=9N2y_e&uQsbf6a=HhEj-h%lRJ!$Wx zls1~cc^>c;D&6(L`e$ROfG$q~_zIP7`(XKNV=Do!t^)8isNJ@~>P=&-0j-*Od;w{< zO}Kd5*{UF`b{=0L+D#8G-E}r9$fBDkmWcYbhmh`CuPVrOH%=@O^-B+-JhWb9fXiW= zC;;-87D9Mvy~qHU!#Gm_^j|HE?%4a10WQXArUC$7EsXBi`;q}J#%ZPk0A4MO?xgpn z<6lgsnUW89buhA%-j|MlF_~scKH$~C$WD4+I>5zb94`U-XA3Sn={#kCtI0fG0`xB} zxbE0_kpV8pdAI_oUt3_^vGFQ{T#pK91yR5Dpt?ifRRp^pB)|fue(eHp?! z<)Y+VQnJ%#2k_jU&a%a)nDLX9T=pA`pRMJfvSGj!IeTUbvqu8EEiE8aqkk}g_VnGA z?$DL)8)Q^3bi*l?LxhcMBGHQ=stBwvWWeH{2{*piQB|L*IQ#!>z&v))Oo|H0~p@58rUBc!4y93-_HQ z$o(U**FV=RNkO2nIgi!%ht?J>Bs?3Q!)BccwH}rrX-GV2MoVo~ z0`UZ%51Rio$npWNj(sHdug><)UN?+fw4>w-0h~Rgc+-BWAXm*Ow?O2tJ)HWr`>KLm zUuU^%M7`1h_8JPXb_CeOskh>9x;9JrSFq2k$5@662%9plQwKsp%VC;37ikLd6i(!)W2+oQ5P-yavv%7b-D8r-kJ=W zYQ^kT3}tI~yAf0;$S6;E*jGt|kXzVCQlAOe!q_018im}7 zMsj5*4f{yiK}~MiFG;3+*-1x{eiN8v- zot*rodmqXEwrFx6m%Ah{Ct_BKRBzfX-EprOB^RTbQ6aUvc0+jB>PFe+ptmfr%0vHB zx!V{9`9i0-Sn|4)iG_Ako{UMQMs%`}L>EI7`Lnc`k;|uYdnF56gDv8*vbuwNS6EQ0 z#w&p@OyK$(w7AYmWg*^{Q&6qsYL$RmNq6;qr1(p1ckAnA*wk9O--)7nz4basS__Xm zG1REHuIwX)9>2AL*-pgNOmAJ=M+(HZzWw=*ZfdT#uI?l0@vZNEu_K#W=&j2tX|1^0 zDbkt1Dx-D%oaL_8imM$?|Lk7rt?TD3cePer?Qr^+?v>!OewpPimx`+mPX5xp5?t0V zv)tuUakat8U%OX|tNL}8yIv_S>YRAfy--}0n^|(P(wx^x@uq#QxTrRh^kN}8Ez{!7 zz*KTxZ03Z;Ty<6z#G8?s=CsyKNsFoMSd_H8p;2{G>1NdBM7b-e%EMfzJ1mERKw%KS z6aupI(TjF-Iqwcu8^XUuP|+A~U`Aik)`3e>Dw#^)>k?>wTPKUOKqlI4X&));1*5(H z=_<+8RJ$vtpgPUPG6S`c?#cRkU>RTS4VPu0oK^^mwTVLVErnlm=SoPo%81B(b@vu_Ix!&g@&PzhiO zG_MTZ0#g9HI_-7|YG7@xbicwhpsrrJTSFRHQ>Q#EF%_U|P#y|^2-GqN4+TiZZyIzb z1z!c}mDbgxu|g?4;~U*$orLRmm1|cUQ8rx+`Wj%oJBOQ^?(1 zsm|)InB6dwT-0nKcXy#WE4x~D!%T8fvxVHlMb% zR(~Z>2~+~{1g@__m+QFjpm0%dE{Dcwp;2!ZQRRNV5pNcKXK>wUcZ;BOJKkw`%eFVR z>ZH46OuL!xl!s;099j(0T_GfX9qx1|%WHS!dLrGGBI=jfPIbDxv4%*e-8GNtMxgdUage-x~hN9^@^+YIiX&yl>54> zf6n!atMxgdUage-x~hN9^@^+YIiX&yl>54>f0^qg7wgM}dbv>U%c}fku9sY_FB9tJ zLb)%i^4GaubH2V#sMmAlzN)HU=X%ZQ`Z}OrPo=w}u6~{DM5oJZf8ctm-4p}j&1fS! zT{P{H#Z>)LjA(bW7uCsS*Bo0;^k1cz@-Y3<9WI8}aA83FD};4t0u$D(Ev446QZ$z$ zCv4eL8Ot@N|8$xDPnEpI)#&I^c?W7vdIM?!$|tQr>G`ce z*{E$K3zTbW8_5P`n%YM4KzXJ^sp(rXjg_e4xT=ymw^-sDxMI$O%x z?chwhtF;6As*{~8b?tU?BHfnS_Q0x@ovZ@d-RMNTDRu3EWg|N-0?NbasD7;s?UBM& zxWA6*PG)=gb3JK}l}_pIDkeLd?f}2bS!cX9VZ}uk3O#GlmQud3QuJ(*2hHxN-AbSm zs01njC4pV9udjCvC?<{9r!a%q%`Z$IndJLF>YKGpZ%M(Nza}Uw{MXZKg)jyHMLm z6VRJH{kNp|n3{Pi((vlGD12givu`aI_M_!Vd{N28{65kc-{O|y<-t^3a;4CrR&5#e z29+UbX*^+8PwiC#l|Us>2~+}=KqXKKR05SiB~S@e0+qm!z|Wu5L-h@X-@o18pxjA5 zza)6i_o7$Wi}Y^P3HnIVn@wHPdo$_!ZQ}JodBgQ->J8$J (0.01 * max_dose) + mask_10pct = dose > (0.10 * max_dose) + mask_50pct = dose > (0.50 * max_dose) + + # Calculate statistics for different regions + stats = { + 'all_voxels': { + 'mean': np.mean(doseUncert), + 'median': np.median(doseUncert), + 'std': np.std(doseUncert) + }, + 'above_1pct': { + 'mean': np.mean(doseUncert[mask_1pct]), + 'median': np.median(doseUncert[mask_1pct]), + 'std': np.std(doseUncert[mask_1pct]) + }, + 'above_10pct': { + 'mean': np.mean(doseUncert[mask_10pct]), + 'median': np.median(doseUncert[mask_10pct]), + 'std': np.std(doseUncert[mask_10pct]) + }, + 'above_50pct': { + 'mean': np.mean(doseUncert[mask_50pct]), + 'median': np.median(doseUncert[mask_50pct]), + 'std': np.std(doseUncert[mask_50pct]) + } + } + + # Create figure with multiple plots + fig = plt.figure(figsize=(15, 15)) + + # Adjust the layout to prevent overlap + plt.subplots_adjust(hspace=0.4, wspace=0.3) + + # Plot 1: All voxels + ax1 = plt.subplot(221) + ax1.hist(doseUncert.flatten(), bins=100, density=True, alpha=0.75, color='skyblue') + ax1.set_title('All Voxels', pad=20, fontsize=12, fontweight='bold') + ax1.set_xlabel('Relative Uncertainty', fontsize=10) + ax1.set_ylabel('Density', fontsize=10) + mean_line = ax1.axvline(stats['all_voxels']['mean'], color='red', linestyle='--', + label=f"Mean: {stats['all_voxels']['mean']:.3f}") + ax1.legend(loc='upper right') + + # Plot 2: >1% max dose + ax2 = plt.subplot(222) + ax2.hist(doseUncert[mask_1pct].flatten(), bins=100, density=True, alpha=0.75, color='lightgreen') + ax2.set_title('Voxels >1% of Max Dose', pad=20, fontsize=12, fontweight='bold') + ax2.set_xlabel('Relative Uncertainty', fontsize=10) + ax2.set_ylabel('Density', fontsize=10) + mean_line = ax2.axvline(stats['above_1pct']['mean'], color='red', linestyle='--', + label=f"Mean: {stats['above_1pct']['mean']:.3f}") + ax2.legend(loc='upper right') + + # Plot 3: >10% max dose + ax3 = plt.subplot(223) + ax3.hist(doseUncert[mask_10pct].flatten(), bins=100, density=True, alpha=0.75, color='salmon') + ax3.set_title('Voxels >10% of Max Dose', pad=20, fontsize=12, fontweight='bold') + ax3.set_xlabel('Relative Uncertainty', fontsize=10) + ax3.set_ylabel('Density', fontsize=10) + mean_line = ax3.axvline(stats['above_10pct']['mean'], color='red', linestyle='--', + label=f"Mean: {stats['above_10pct']['mean']:.3f}") + ax3.legend(loc='upper right') + + # Plot 4: >50% max dose + ax4 = plt.subplot(224) + ax4.hist(doseUncert[mask_50pct].flatten(), bins=100, density=True, alpha=0.75, color='plum') + ax4.set_title('Voxels >50% of Max Dose', pad=20, fontsize=12, fontweight='bold') + ax4.set_xlabel('Relative Uncertainty', fontsize=10) + ax4.set_ylabel('Density', fontsize=10) + mean_line = ax4.axvline(stats['above_50pct']['mean'], color='red', linestyle='--', + label=f"Mean: {stats['above_50pct']['mean']:.3f}") + ax4.legend(loc='upper right') + + # Add grid to all plots + for ax in [ax1, ax2, ax3, ax4]: + ax.grid(True, alpha=0.3) + + # Add a main title to the figure + plt.suptitle('Dose Uncertainty Distribution Analysis', fontsize=14, fontweight='bold', y=0.95) + + plt.show() + ''' + + top_indices = np.argsort(dose.flatten())[-10:] + max_dose_median = np.median(dose.flatten()[top_indices]) + + fifty_mask = dose >= 0.5 * (max_dose_median) + fifty_uncert = doseUncert[fifty_mask] * 100 + + ninety_mask = dose >= 0.9 * (max_dose_median) + ninety_uncert = doseUncert[ninety_mask] * 100 + + + + results = { + 'top50%_mean': np.mean(fifty_uncert), + 'top50%_median': np.median(fifty_uncert), + 'top90%_mean': np.mean(ninety_uncert), + 'top90%_median': np.median(ninety_uncert), + 'voxel_size_x': doseSpacing[0], + 'voxel_size_y': doseSpacing[1], + 'voxel_size_z': doseSpacing[2], + 'n_primaries': nsim, + 'total_voxels': dose.size + } + + output_path = os.path.join(os.getcwd(), 'uncertainty_metrics.txt') + with open(output_path, 'w') as f: + headers = '\t'.join(results.keys()) + values = '\t'.join(str(x) for x in results.values()) + f.write(headers + '\n') + f.write(values + '\n') + + print(f"File saved to: {output_path}") + return results -if __name__=="__main__": - main() +#if __name__=="__main__": +# main()