Merge 252a524 into 5f749b8

pylinac/__init__.py

@@ -13,6 +13,7 @@
 from pylinac.core import decorators, geometry, image, io, mask, profile, roi, utilities
 from pylinac.core.utilities import clear_data_files, assign2machine
 from pylinac.flatsym import FlatSym
+from pylinac.fieldparams import FieldParams
 from pylinac.planar_imaging import LeedsTOR, StandardImagingQC3, LasVegas, DoselabMC2kV, DoselabMC2MV
 from pylinac.log_analyzer import load_log, Dynalog, TrajectoryLog, MachineLogs
 from pylinac.picketfence import PicketFence  # must be after log analyzer

pylinac/core/hillreg.py

@@ -0,0 +1,46 @@
+"""Perform non-linear regression using a Hill function"""
+import numpy as np
+from scipy.optimize import curve_fit
+from scipy.optimize import differential_evolution
+import warnings
+
+
+def hill_reg(xData: np.ndarray, yData: np.ndarray):
+
+    def hill_func(x, a, b, c, d):  # Hill function
+        """
+            a: sigmoid low level
+            b: sigmoid high level
+            c: approximate inflection point
+            d: slope of the sigmoid
+        """
+        return a + (b-a)/(1.0 + (c/x)**d)
+
+    # function for genetic algorithm to minimize (sum of squared error)
+    def sumOfSquaredError(parameterTuple):
+        warnings.filterwarnings("ignore")  # do not print warnings by genetic algorithm
+        val = hill_func(xData, *parameterTuple)
+        return np.sum((yData - val) ** 2.0)
+
+    def generate_initial_parameters(xData, yData):
+        # min and max used for bounds
+        maxX = max(xData)
+        minX = min(xData)
+        maxY = max(yData)
+
+        parameterBounds = []
+        parameterBounds.append([0, maxY])  # search bounds for a
+        parameterBounds.append([0, maxY])  # search bounds for b
+        parameterBounds.append([minX, maxX])  # search bounds for c
+        parameterBounds.append([-100, 100])  # search bounds for d
+
+        # "seed" the numpy random number generator for repeatable results
+        result = differential_evolution(sumOfSquaredError, parameterBounds, seed=3)
+        return result.x
+
+    # generate initial parameter values
+    genetic_parameters = generate_initial_parameters(xData, yData)
+
+    # curve fit the data
+    fitted_parameters, pcov = curve_fit(hill_func, xData, yData, genetic_parameters)
+    return fitted_parameters

pylinac/core/profile.py

@@ -138,6 +138,7 @@ class SingleProfile(ProfileMixin):
     interpolation_factor: int = 100
     interpolation_type: str = 'linear'
     _values: np.ndarray
+    _dpmm: float
 
     def __init__(self, values: np.ndarray):
         """
@@ -159,6 +160,31 @@ def values(self, value):
             raise TypeError("Values must be a numpy array")
         self._values = value.astype(float)
 
+    @property
+    def dpmm(self) -> NumberLike:
+        """The Dots-per-mm of the profile, defined at isocenter. E.g. if an EPID image is taken at 150cm SID,
+        the dpmm will scale back to 100cm."""
+        return self._dpmm
+
+    @dpmm.setter
+    def dpmm(self, value):
+        """The Dots-per-mm of the profile, defined at isocenter. E.g. if an EPID image is taken at 150cm SID,
+        the dpmm will scale back to 100cm."""
+        if value > 0:
+            self._dpmm = value
+
+    def center(self) -> Tuple[NumberLike, NumberLike]:
+        """Returns the center index and value of the profile. If the profile has an even number of values the centre
+        lies between the two centre indices and the centre value is the average of the two centre values else the
+        centre index and value are returned. Added by ACC 3/12/2020"""
+        plen = self.values.shape[0]
+        if plen % 2 == 0:  # plen is even and central detectors straddle CAX
+            cax = (self.values[int(plen / 2)] + self.values[int(plen / 2) - 1]) / 2.0
+        else:  # plen is odd and we have a central detector
+            cax = self.values[int((plen - 1) / 2)]
+        plen = (plen - 1)/2.0
+        return plen, cax
+
     @property
     @lru_cache()
     def _values_interp(self) -> np.ndarray:
@@ -210,6 +236,64 @@ def fwxm_center(self, x: int=50, interpolate: bool=False) -> Tuple[NumberLike, N
             fwxm_center_idx = int(round(fwxm_center_idx))
             return fwxm_center_idx, self.values[fwxm_center_idx]
 
+    @argue.bounds(x=(0, 100), ifa=(0, 1.0))
+    @argue.options(norm=('max', 'max grounded', 'cax', 'cax grounded'), interpolate=(True, False))
+    def field_edges(self, ifa: float=1.0, x: int=50, norm='max grounded', interpolate=False) -> Tuple[float, float]:
+        """Return the width at X-Max, where X is the percentage height.
+
+        Parameters
+        ----------
+        x : int
+            The dose level of the profile. E.g. x = 50 is 50% height, 100 is the profile peak.
+            i.e. FWHM.
+        ifa : float
+            In Field Area: 1.0 is the entire field at the dose level.
+        norm : str
+            The type of normalisation to apply:
+                'max'          : normalises to the profile maximum with no grounding
+                'max grounded' : grounds each side to the profile minimum and normalises to the profile maximum
+                'cax'          : normalises to the profile centre value with no grounding
+                'cax grounded' : grounds each side to the profile minimum and normalises to the profile centre value
+            Default is 'max_grounded' as this is the intrinsic method of the scipy sgnal.find_peak function.
+        interpolate : bool
+            Whether to interpolate the profile array values to get subpixel precision.
+
+        Returns
+        -------
+        left index, right index
+            The left and right indices of the in field area at the dose level.
+        """
+
+        # prevent array from being grounded by setting ends to zero
+        if norm in ['max', 'cax']:
+            ydata = np.insert(self.values, 0, 0)
+            ydata = np.append(ydata, 0)
+        else:
+            ydata = self.values
+
+        # adjust x to give the equivalent level if peak val was at CAX
+        if norm == 'cax':
+            _, cax_val = self.center()
+            ymax = ydata.max()
+            x = x*cax_val/ymax
+
+        _, peak_props = find_peaks(ydata, fwxm_height=x/100, max_number=1)
+        left = peak_props['left_ips'][0]
+        right = peak_props['right_ips'][0]
+        if ifa < 1.0:
+            delta = (1.0 - ifa)*(right - left)/2.0
+            left = left + delta
+            right = right - delta
+
+        if not interpolate:
+            left = round(left)
+            right = round(right)
+
+        if norm in ['max', 'cax']:
+            left = left - 1
+            right = right - 1
+        return left, right
+
     @argue.bounds(lower=(0, 100), upper=(0, 100))
     def penumbra_width(self, lower: int=20, upper: int=80) -> Tuple[float, float]:
         """Return the penumbra width of the profile.
@@ -259,7 +343,7 @@ def field_values(self, field_width: float=0.8) -> np.ndarray:
         return field_values
 
     @argue.bounds(field_width=(0, 1))
-    def field_edges(self, field_width: float=0.8, interpolate: bool=False) -> Tuple[NumberLike, NumberLike]:
+    def field_edges_deprecated(self, field_width: float=0.8, interpolate: bool=False) -> Tuple[NumberLike, NumberLike]:
         """Return the indices of the field width edges, based on the FWHM.
 
         See Also