Changed Gamma Analysis to include local Gamma calc. Uncertainty Stats…

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,14 +245,16 @@ 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"
         if dose2water in [basename(f) for f in mergedfiles]:
             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
+
+
+
 
 
 

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

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:

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 )

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)

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