ARRAY = 'Array'
DICOM = 'DICOM'
IMAGE = 'Image'
FILE_TYPE = 'file'
STREAM_TYPE = 'stream'
MM_PER_INCH = 25.4
def equate_images(image1: ImageLike, image2: ImageLike) -> Tuple[ImageLike, ImageLike]:
def is_image(path: Union[str, io.BytesIO, ImageLike, np.ndarray]) -> bool:
def retrieve_image_files(path: str) -> List[str]:
def load(path: Union[str, ImageLike, np.ndarray, BinaryIO], **kwargs) -> ImageLike:
def load_url(url: str, progress_bar: bool = True, **kwargs) -> ImageLike:
def load_multiples(image_file_list: Sequence, method: str = 'mean', stretch_each: bool = True, **kwargs) -> ImageLike:
def _is_dicom(path: Union[str, io.BytesIO, ImageLike, np.ndarray]) -> bool:
def _is_image_file(path: str) -> bool:
def _is_array(obj: Any) -> bool:
    def __init__(self, path: Union[str, BytesIO, ImageLike, np.ndarray, BufferedReader]):
    def truncated_path(self) -> str:
    def from_multiples(cls, filelist: List[str], method: str='mean', stretch: bool=True, **kwargs) -> ImageLike:
    def center(self) -> Point:
    def physical_shape(self) -> Tuple[float, float]:
    def date_created(self, format: str="%A, %B %d, %Y") -> str:
    def plot(self, ax: plt.Axes=None, show: bool=True, clear_fig: bool=False, **kwargs) -> plt.Axes:
    def filter(self, size: Union[float, int]=0.05, kind: str='median') -> None:
    def crop(self, pixels: int=15, edges: Tuple[str, ...]=('top', 'bottom', 'left', 'right')) -> None:
    def remove_edges(self, pixels: int=15, edges: Tuple[str, ...]=('top', 'bottom', 'left', 'right')) -> None:
    def flipud(self) -> None:
    def fliplr(self) -> None:
    def invert(self) -> None:
    def roll(self, direction: str='x', amount: int=1) -> None:
    def rot90(self, n: int=1) -> None:
    def threshold(self, threshold: int, kind: str='high') -> None:
    def as_binary(self, threshold: int) -> ImageLike:
    def dist2edge_min(self, point: Union[Point, Tuple]) -> float:
    def ground(self) -> float:
    def normalize(self, norm_val: Union[str, NumberLike]='max') -> None:
    def check_inversion(self, box_size: int=20, position: Sequence=(0.0, 0.0)) -> None:
    def check_inversion_by_histogram(self, percentiles=(5, 50, 95)) -> bool:
    def gamma(self, comparison_image: ImageLike, doseTA: NumberLike=1, distTA: NumberLike=1,
              threshold: NumberLike=0.1, ground: bool=True, normalize: bool=True) -> np.ndarray:
    def as_type(self, dtype) -> np.ndarray:
    def shape(self) -> Tuple[int, int]:
    def size(self) -> int:
    def ndim(self) -> int:
    def dtype(self) -> np.dtype:
    def sum(self) -> float:
    def ravel(self) -> np.ndarray:
    def flat(self) -> np.ndarray:
    def __len__(self):
    def __getitem__(self, item):
    def __init__(self, path: Union[str, BytesIO, BufferedReader], *, dtype=None, dpi: NumberLike=None, sid: NumberLike=None):
    def save(self, filename: str) -> str:
    def sid(self) -> NumberLike:
    def dpi(self) -> NumberLike:
    def dpmm(self) -> NumberLike:
    def cax(self) -> Point:
    def __init__(self, path: str, use_filenames: bool=False):
    def gantry_angle(self) -> float:
    def collimator_angle(self) -> float:
    def couch_angle(self) -> float:
    def _get_axis_value(self, axis_str: str, axis_dcm_attr: str) -> float:
    def __init__(self, path: str, *, dpi: NumberLike=None, sid: NumberLike=None, dtype=None):
    def dpi(self) -> float:
    def dpmm(self) -> Optional[float]:
    def __init__(self, array: np.array, *, dpi: NumberLike=None, sid: NumberLike=None, dtype=None):
    def dpmm(self) -> Optional[float]:
    def dpi(self) -> Optional[float]:
    def __sub__(self, other):
    def __init__(self, folder: str, dtype=None, min_number: int=39, check_uid: bool=True):
    def from_zip(cls, zip_path: str, dtype=None):
    def is_CT_slice(file: str) -> bool:
    def _check_number_and_get_common_uid_imgs(self, min_number: int) -> List:
    def plot(self, slice: int=0) -> None:
    def metadata(self) -> pydicom.FileDataset:
    def __getitem__(self, item) -> DicomImage:
    def __setitem__(self, key, value: DicomImage):
    def __len__(self):
class BaseImage:
class DicomImage(BaseImage):
class LinacDicomImage(DicomImage):
class FileImage(BaseImage):
class ArrayImage(BaseImage):
class DicomImageStack: class PylinacCanvas:
    def __init__(self, filename: str, page_title: str, font: str="Helvetica", metadata: dict=None, metadata_location: tuple=(2, 25.5)):
    def add_new_page(self) -> None:
    def _add_metadata(self) -> None:
    def _generate_pylinac_template_theme(self) -> None:
    def add_text(self, text: Union[str, List[str]], location: Tuple[NumberLike, NumberLike], font_size: int=10) -> None:
    def add_image(self, image_data: io.BytesIO, location: Sequence, dimensions: Sequence, preserve_aspect_ratio: bool=True) -> None:
    def finish(self) -> None: class Point:
    """A geometric point with x, y, and z coordinates/attributes."""
    z: Union[int, float]
    y: Union[int, float]
    x: Union[int, float]
    _attr_list: List[str] = ['x', 'y', 'z', 'idx', 'value']
    _coord_list: List[str] = ['x', 'y', 'z']

    def __init__(self, x: Union[NumberLike, tuple, 'Point'] = 0, y: NumberLike = 0, z: NumberLike = 0, idx: Optional[int]=None,
                 value: Optional[NumberLike]=None, as_int: bool=False):
        """
        Parameters
        ----------
        x : number-like, Point, iterable
            x-coordinate or iterable type containing all coordinates. If iterable, values are assumed to be in order: (x,y,z).
        y : number-like, optional
            y-coordinate
        idx : int, optional
            Index of point. Useful for sequential coordinates; e.g. a point on a circle profile is sometimes easier to describe
            in terms of its index rather than x,y coords.
        value : number-like, optional
            value at point location (e.g. pixel value of an image)
        as_int : boolean
            If True, coordinates are converted to integers.
        """
        if isinstance(x, Point):
            for attr in self._attr_list:
                item = getattr(x, attr, None)
                setattr(self, attr, item)
        elif is_iterable(x):
            for attr, item in zip_longest(self._attr_list, x, fillvalue=0):
                setattr(self, attr, item)
        else:
            self.x = x
            self.y = y
            self.z = z
            self.idx = idx
            self.value = value

        if as_int:
            self.x = int(round(self.x))
            self.y = int(round(self.y))
            self.z = int(round(self.z))

    def distance_to(self, thing: Union['Point', 'Circle']) -> float:
        """Calculate the distance to the given point.

        Parameters
        ----------
        thing : Circle, Point, 2 element iterable
            The other thing to calculate distance to.
        """
        if isinstance(thing, Circle):
            return abs(np.sqrt((self.x - thing.center.x)**2 + (self.y - thing.center.y)**2) - thing.radius)
        p = Point(thing)
        return math.sqrt((self.x - p.x)**2 + (self.y - p.y)**2 + (self.z - p.z)**2)

    def as_array(self, only_coords: bool=True) -> np.array:
        """Return the point as a numpy array."""
        if only_coords:
            return np.array([getattr(self, item) for item in self._coord_list])
        else:
            return np.array([getattr(self, item) for item in self._attr_list if (getattr(self, item) is not None)])

    def __repr__(self):
        return f"Point(x={self.x:3.2f}, y={self.y:3.2f}, z={self.z:3.2f})"

    def __eq__(self, other):
        # if all attrs equal, points considered equal
        return all(getattr(self, attr) == getattr(other, attr) for attr in self._attr_list)

    def __sub__(self, other):
        p = Point()
        for attr in self._attr_list:
            try:
                diff = getattr(self, attr) - getattr(other, attr)
            except TypeError:
                diff = None
            setattr(p, attr, diff)
        return p

    def __mul__(self, other):
        for attr in self._attr_list:
            try:
                self.__dict__[attr] *= other
            except TypeError:
                pass

    def __truediv__(self, other):
        for attr in self._attr_list:
            val = getattr(self, attr)
            try:
                setattr(self, attr, val/other)
            except TypeError:
                pass
        return self class Line:
    """A line that is represented by two points or by m*x+b.

    Notes
    -----
    Calculations of slope, etc are from here:
    http://en.wikipedia.org/wiki/Linear_equation
    and here:
    http://www.mathsisfun.com/algebra/line-equation-2points.html
    """
    point1: Point
    point2: Point

    def __init__(self, point1: Union[Point, Tuple[float, float]], point2: Union[Point, Tuple[float, float]]):
        """
        Parameters
        ----------
        point1 : Point
            One point of the line
        point2 : Point
            Second point along the line.
        """
        self.point1 = Point(point1)
        self.point2 = Point(point2)

    def __repr__(self):
        return f'Line: p1:(x={self.point1.x:.1f}, y={self.point1.y:.1f}, z={self.point1.z:.1f}), ' \
                   f'p2:(x={self.point2.x:.1f}, y={self.point2.y:.1f}, z={self.point2.z:.1f})'

    @property
    def m(self) -> float:
        """Return the slope of the line.

        m = (y1 - y2)/(x1 - x2)

        From: http://www.purplemath.com/modules/slope.htm
        """
        return (self.point1.y - self.point2.y) / (self.point1.x - self.point2.x)

    @property
    def b(self) -> float:
        """Return the y-intercept of the line.

        b = y - m*x
        """
        return self.point1.y - (self.m * self.point1.x)

    def y(self, x) -> float:
        """Return y-value along line at position x."""
        return self.m * x + self.b

    def x(self, y) -> float:
        """Return x-value along line at position y."""
        return (y - self.b)/self.m

    @property
    def center(self) -> Point:
        """Return the center of the line as a Point."""
        mid_x = np.abs((self.point2.x - self.point1.x)/2 + self.point1.x)
        mid_y = (self.point2.y - self.point1.y) / 2 + self.point1.y
        return Point(mid_x, mid_y)

    @property
    def length(self) -> float:
        """Return length of the line, if finite."""
        return self.point1.distance_to(self.point2)

    def distance_to(self, point) -> float:
        """Calculate the minimum distance from the line to a point.

        Equations are from here: http://mathworld.wolfram.com/Point-LineDistance2-Dimensional.html #14

        Parameters
        ----------
        point : Point, iterable
            The point to calculate distance to.
        """
        point = Point(point).as_array()
        lp1 = self.point1.as_array()
        lp2 = self.point2.as_array()
        numerator = np.sqrt(np.sum(np.power(np.cross((lp2 - lp1), (lp1 - point)), 2)))
        denominator = np.sqrt(np.sum(np.power(lp2 - lp1, 2)))
        return numerator/denominator

    def plot2axes(self, axes: plt.Axes, width: NumberLike=1, color: str='w') -> None:
        """Plot the line to an axes.

        Parameters
        ----------
        axes : matplotlib.axes.Axes
            An MPL axes to plot to.
        color : str
            The color of the line.
        """
        axes.plot((self.point1.x, self.point2.x), (self.point1.y, self.point2.y), linewidth=width, color=color) class Circle:
    """A geometric circle with center Point, radius, and diameter."""
    center: Point
    radius: float

    def __init__(self, center_point: Union[Point, Iterable]=(0, 0), radius: float = 0):
        """
        Parameters
        ----------
        center_point : Point, optional
            Center point of the wobble circle.
        radius : float, optional
            Radius of the wobble circle.
        """
        if center_point is None:
            center_point = Point()
        elif isinstance(center_point, Point) or is_iterable(center_point):
            center_point = Point(center_point)
        else:
            raise TypeError("Circle center must be of type Point or iterable")

        self.center = center_point
        self.radius = radius

    @property
    def diameter(self) -> float:
        """Get the diameter of the circle."""
        return self.radius*2

    def plot2axes(self, axes, edgecolor: str='black', fill: bool=False) -> None:
        """Plot the Circle on the axes.

        Parameters
        ----------
        axes : matplotlib.axes.Axes
            An MPL axes to plot to.
        edgecolor : str
            The color of the circle.
        fill : bool
            Whether to fill the circle with color or leave hollow.
        """
        axes.add_patch(mpl_Circle((self.center.x, self.center.y), edgecolor=edgecolor, radius=self.radius, fill=fill)) def get_url(url: str, destination: str=None, progress_bar: bool=True) -> str:
    """Download a URL to a local file.

    Parameters
    ----------
    url : str
        The URL to download.
    destination : str, None
        The destination of the file. If None is given the file is saved to a temporary directory.
    progress_bar : bool
        Whether to show a command-line progress bar while downloading.

    Returns
    -------
    filename : str
        The location of the downloaded file.

    Notes
    -----
    Progress bar use/example adapted from tqdm documentation: https://github.com/tqdm/tqdm
    """

    def my_hook(t):
        last_b = [0]

        def inner(b=1, bsize=1, tsize=None):
            if tsize is not None:
                t.total = tsize
            if b > 0:
                t.update((b - last_b[0]) * bsize)
            last_b[0] = b
        return inner

    try:
        if progress_bar:
            with tqdm(unit='B', unit_scale=True, miniters=1, desc=url.split('/')[-1]) as t:
                filename, _ = urlretrieve(url, filename=destination, reporthook=my_hook(t))
        else:
            filename, _ = urlretrieve(url, filename=destination)
    except (HTTPError, URLError, ValueError) as e:
        raise e
    return filename class TemporaryZipDirectory(TemporaryDirectory):
    """Creates a temporary directory that unpacks a ZIP archive."""
    def __init__(self, zfile):
        """
        Parameters
        ----------
        zfile : str
            String that points to a ZIP archive.
        """
        super().__init__()
        zfiles = zipfile.ZipFile(zfile)
        zfiles.extractall(path=self.name) def retrieve_demo_file(url: str, force: bool = False) -> str:
    """Retrieve the demo file either by getting it from file or from a URL.

    If the file is already on disk it returns the file name. If the file isn't
    on disk, get the file from the URL and put it at the expected demo file location
    on disk for lazy loading next time.

    Parameters
    ----------
    url : str
        The suffix to the url (location within the S3 bucket) pointing to the demo file.
    """
    true_url = r'https://storage.googleapis.com/pylinac_demo_files/' + url
    demo_file = osp.join(osp.dirname(osp.dirname(__file__)), 'demo_files', url)
    demo_dir = osp.dirname(demo_file)
    if not osp.exists(demo_dir):
        os.makedirs(demo_dir)
    if force or not osp.isfile(demo_file):
        get_url(true_url, destination=demo_file)
    return demo_file class SingleProfile(ProfileMixin):
    """A profile that has one large signal, e.g. a radiation beam profile.
    Signal analysis methods are given, mostly based on FWXM and on Hill function calculations.
    Profiles with multiple peaks are better suited by the MultiProfile class.
    """

    def __init__(self, values: np.ndarray, dpmm: float = None,
                 interpolation: Union[Interpolation, str, None] = Interpolation.LINEAR,
                 ground: bool = True,
                 interpolation_resolution_mm: float = 0.1,
                 interpolation_factor: float = 10,
                 normalization_method: Union[Normalization, str] = Normalization.BEAM_CENTER,
                 edge_detection_method: Union[Edge, str] = Edge.FWHM,
                 edge_smoothing_ratio: float = 0.003,
                 hill_window_ratio: float = 0.1):
        """
        Parameters
        ----------
        values
            The profile numpy array. Must be 1D.
        dpmm
            The dots (pixels) per mm. Pass to get info like beam width in distance units in addition to pixels
        interpolation
            Interpolation technique.
        ground
            Whether to ground the profile (set min value to 0). Helpful most of the time.
        interpolation_resolution_mm
            The resolution that the interpolation will scale to. **Only used if dpmm is passed and interpolation is set**.
            E.g. if the dpmm is 0.5 and the resolution is set to 0.1mm the data will be interpolated to have a new dpmm of 10 (1/0.1).
        interpolation_factor
            The factor to multiply the data by. **Only used if interpolation is used and dpmm is NOT passed**. E.g. 10
            will perfectly decimate the existing data according to the interpolation method passed.
        normalization_method
            How to pick the point to normalize the data to.
        edge_detection_method
            The method by which to detect the field edge. FWHM is reasonable most of the time except for FFF beams.
            Inflection-derivative will use the max gradient to determine the field edge. Note that this may not be the
            50% height. In fact, for FFF beams it shouldn't be. Inflection methods are better for FFF and other unusual
            beam shapes.
        edge_smoothing_ratio
            **Only applies to INFLECTION_DERIVATIVE and INFLECTION_HILL.**

            The ratio of the length of the values to use as the sigma for a Gaussian filter applied before searching for
            the inflection. E.g. 0.005 with a profile of 1000 points will result in a sigma of 5.
            This helps make the inflection point detection more robust to noise. Increase for noisy data.
        hill_window_ratio
            The ratio of the field size to use as the window to fit the Hill function. E.g. 0.2 will using a window
            centered about each edge with a width of 20% the size of the field width. **Only applies when the edge
            detection is INFLECTION_HILL**.
        """
        self._interp_method = convert_to_enum(interpolation, Interpolation)
        self._interpolation_res = interpolation_resolution_mm
        self._interpolation_factor = interpolation_factor
        self._norm_method = convert_to_enum(normalization_method, Normalization)
        self._edge_method = convert_to_enum(edge_detection_method, Edge)
        self._edge_smoothing_ratio = edge_smoothing_ratio
        self._hill_window_ratio = hill_window_ratio
        self.values = values  # set initial data so we can do things like find beam center
        self.dpmm = dpmm
        fitted_values, new_dpmm, x_indices = self._interpolate(values, dpmm, interpolation_resolution_mm,
                                                               interpolation_factor, self._interp_method)
        self.dpmm = new_dpmm  # update as needed
        self.values = fitted_values
        self.x_indices = x_indices
        if ground:
            fitted_values -= fitted_values.min()
        norm_values = self._normalize(fitted_values, self._norm_method)
        self.values = norm_values  # update values

    @staticmethod
    def _interpolate(values, dpmm, interpolation_resolution, interpolation_factor, interp_method: Interpolation) -> (
            np.ndarray, float, float, float):
        """Fit the data to the passed interpolation method. Will also calculate the new values to correct the measurements such as dpmm"""
        x_indices = list(range(len(values)))
        if interp_method == Interpolation.NONE:
            return values, dpmm, x_indices  # do nothing
        elif interp_method == Interpolation.LINEAR:
            if dpmm is not None:
                samples = int(round((len(x_indices)-1)/(dpmm*interpolation_resolution)))
                new_dpmm = 1/interpolation_resolution
            else:
                samples = int(round((len(x_indices)-1)*interpolation_factor))
                new_dpmm = None
            f = interp1d(x_indices, values, kind='linear', bounds_error=True)
            new_x = np.linspace(0, len(x_indices)-1, num=samples)
            return f(new_x), new_dpmm, new_x
        elif interp_method == Interpolation.SPLINE:
            if dpmm is not None:
                samples = int(round((len(x_indices)-1)/(dpmm*interpolation_resolution)))
                new_dpmm = 1 / interpolation_resolution
            else:
                samples = int(round((len(x_indices)-1)*interpolation_factor))
                new_dpmm = None
            f = interp1d(x_indices, values, kind='cubic')
            new_x = np.linspace(0, len(x_indices)-1, num=samples)
            return f(new_x), new_dpmm, new_x

    def _normalize(self, values, method: Normalization) -> np.ndarray:
        """Normalize the data given a method."""
        if method == Normalization.NONE:
            return values
        elif method == Normalization.MAX:
            return values / values.max()
        elif method == Normalization.GEOMETRIC_CENTER:
            return values / self._geometric_center(values)['value (exact)']
        elif method == Normalization.BEAM_CENTER:
            return values / self.beam_center()['value (@rounded)']

    def _geometric_center(self, values) -> dict:
        """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."""
        plen = values.shape[0]
        # buffer overflow can cause the below addition to give strange results
        values = values.astype(np.float64)
        if plen % 2 == 0:  # plen is even and central detectors straddle CAX
            cax = (values[int(plen / 2)] + values[int(plen / 2) - 1]) / 2.0
        else:  # plen is odd and we have a central detector
            cax = values[int((plen - 1) / 2)]
        plen = (plen - 1)/2.0
        return {'index (exact)': plen, 'value (exact)': cax}

    def geometric_center(self) -> dict:
        """The geometric center (i.e. the device center)"""
        return self._geometric_center(self.values)

    def beam_center(self) -> dict:
        """The center of the detected beam. This can account for asymmetries in the beam position (e.g. offset jaws)"""
        if self._edge_method == Edge.FWHM:
            data = self.fwxm_data(x=50)
            return {'index (rounded)': data['center index (rounded)'],
                    'index (exact)': data['center index (exact)'],
                    'value (@rounded)': data['center value (@rounded)']}
        elif self._edge_method in (Edge.INFLECTION_DERIVATIVE, Edge.INFLECTION_HILL):
            infl = self.inflection_data()
            mid_point = infl['left index (exact)'] + (infl['right index (exact)'] - infl['left index (exact)']) / 2
            return {'index (rounded)': int(round(mid_point)),
                    'index (exact)': mid_point,
                    'value (@rounded)': self.values[int(round(mid_point))]}

    @argue.bounds(x=(0, 100))
    def fwxm_data(self, x: int = 50) -> dict:
        """Return the width at X-Max, where X is the percentage height.

        Parameters
        ----------
        x
            The percent height of the profile. E.g. x = 50 is 50% height,
            i.e. FWHM.
        """
        _, peak_props = find_peaks(self.values, fwxm_height=x/100, max_number=1)
        left_idx = peak_props['left_ips'][0]
        right_idx = peak_props['right_ips'][0]
        fwxm_center_idx = ((peak_props['right_ips'][0] - peak_props['left_ips'][0]) / 2 + peak_props['left_ips'][0])

        data = {'width (exact)': peak_props['widths'][0],
                'width (rounded)': int(round(right_idx)) - int(round(left_idx)),
                'center index (rounded)': int(round(fwxm_center_idx)),
                'center index (exact)': fwxm_center_idx,
                'center value (@rounded)': self.values[int(round(fwxm_center_idx))],
                'left index (exact)': left_idx,
                'left index (rounded)': int(round(left_idx)),
                'left value (@rounded)': self.values[int(round(left_idx))],
                'right index (exact)': right_idx,
                'right index (rounded)': int(round(right_idx)),
                'right value (@rounded)': self.values[int(round(right_idx))],
                'field values': self.values[int(round(left_idx)):
                                            int(round(right_idx))],
                'peak_props': peak_props}
        if self.dpmm:
            data['width (exact) mm'] = data['width (exact)'] / self.dpmm
            data['left distance (exact) mm'] = abs(data['center index (exact)'] - data['left index (exact)']) / self.dpmm
            data['right distance (exact) mm'] = abs(data['right index (exact)'] - data['center index (exact)']) / self.dpmm

        return data

    @argue.bounds(in_field_ratio=(0, 1.0), slope_exclusion_ratio=(0, 1.0))
    def field_data(self, in_field_ratio: float = 0.8, slope_exclusion_ratio=0.2) -> dict:
        """Return the width at X-Max, where X is the percentage height.

        Parameters
        ----------
        in_field_ratio
            In Field Ratio: 1.0 is the entire detected field; 0.8 would be the central 80%, etc.
        slope_exclusion_ratio
            Ratio of the field width to use as the cutoff between "top" calculation and "slope" calculation. Useful for FFF beams.
            This area is centrally located in the field. E.g. 0.2 will use the central 20% of the field to calculate
            the "top" value. To calculate the slope of each side, the field width between the edges of the in_field_ratio
            and the slope exclusion region are used.

            .. warning:: The "top" value is always calculated. For FFF beams this should be reasonable, but for flat beams
                         this value may end up being non-sensible.
        """
        if slope_exclusion_ratio >= in_field_ratio:
            raise ValueError("The exclusion region must be smaller than the field ratio")
        if self._edge_method == Edge.FWHM:
            data = self.fwxm_data(x=50)
            beam_center_idx = data['center index (exact)']
            full_width = data['width (exact)']

        elif self._edge_method in (Edge.INFLECTION_DERIVATIVE, Edge.INFLECTION_HILL):
            infl_data = self.inflection_data()
            beam_center_idx = self.beam_center()['index (exact)']
            full_width = infl_data['right index (exact)'] - infl_data['left index (exact)']
        beam_center_idx_r = int(round(beam_center_idx))

        cax_idx = self.geometric_center()['index (exact)']
        cax_idx_r = int(round(cax_idx))

        field_left_idx = beam_center_idx - in_field_ratio * full_width / 2
        field_left_idx_r = int(round(field_left_idx))
        field_right_idx = beam_center_idx + in_field_ratio * full_width / 2
        field_right_idx_r = int(round(field_right_idx))
        field_width = field_right_idx - field_left_idx

        # slope calcs
        inner_left_idx = beam_center_idx - slope_exclusion_ratio*field_width/2
        inner_left_idx_r = int(round(inner_left_idx))
        inner_right_idx = beam_center_idx + slope_exclusion_ratio*field_width/2
        inner_right_idx_r = int(round(inner_right_idx))
        left_fit = linregress(range(field_left_idx_r, inner_left_idx_r),
                              self.values[field_left_idx_r:inner_left_idx_r])
        right_fit = linregress(range(inner_right_idx_r, field_right_idx_r),
                               self.values[inner_right_idx_r:field_right_idx_r])

        # top calc
        fit_params = np.polyfit(range(inner_left_idx_r, inner_right_idx_r),
                                self.values[inner_left_idx_r:inner_right_idx_r], deg=2)
        width = abs(inner_right_idx_r - inner_left_idx_r)

        def poly_func(x):
            # return the negative since we're MINIMIZING and want the top value
            return -(fit_params[0] * (x ** 2) + fit_params[1] * x + fit_params[2])

        # minimize the polynomial function
        min_f = minimize(poly_func, x0=(inner_left_idx_r+width/2,), bounds=((inner_left_idx_r, inner_right_idx_r),))
        top_idx = min_f.x[0]
        top_val = -min_f.fun

        data = {'width (exact)': field_width,
                'beam center index (exact)': beam_center_idx,
                'beam center index (rounded)': beam_center_idx_r,
                'beam center value (@rounded)': self.values[int(round(beam_center_idx))],
                'cax index (exact)': cax_idx,
                'cax index (rounded)': cax_idx_r,
                'cax value (@rounded)': self.values[int(round(cax_idx))],
                'left index (exact)': field_left_idx,
                'left index (rounded)': field_left_idx_r,
                'left value (@rounded)': self.values[int(round(field_left_idx))],
                'left slope': left_fit.slope,
                'left intercept': left_fit.intercept,
                'right slope': right_fit.slope,
                'right intercept': right_fit.intercept,
                'left inner index (exact)': inner_left_idx,
                'left inner index (rounded)': inner_left_idx_r,
                'right inner index (exact)': inner_right_idx,
                'right inner index (rounded)': inner_right_idx_r,
                '"top" index (exact)': top_idx,
                '"top" index (rounded)': int(round(top_idx)),
                '"top" value (@exact)': top_val,
                'top params': fit_params,
                'right index (exact)': field_right_idx,
                'right index (rounded)': field_right_idx_r,
                'right value (@rounded)': self.values[int(round(field_right_idx))],
                'field values': self.values[int(round(field_left_idx)):
                                            int(round(field_right_idx))]}
        if self.dpmm:
            data['width (exact) mm'] = data['width (exact)'] / self.dpmm
            data['left slope (%/mm)'] = data['left slope'] * self.dpmm * 100
            data['right slope (%/mm)'] = data['right slope'] * self.dpmm * 100
            data['left distance->beam center (exact) mm'] = abs(data['beam center index (exact)'] - data['left index (exact)']) / self.dpmm
            data['right distance->beam center (exact) mm'] = abs(data['right index (exact)'] - data['beam center index (exact)']) / self.dpmm
            data['left distance->CAX (exact) mm'] = abs(data['cax index (exact)'] - data['left index (exact)']) / self.dpmm
            data['right distance->CAX (exact) mm'] = abs(data['cax index (exact)'] - data['right index (exact)']) / self.dpmm

            data['left distance->top (exact) mm'] = abs(data['"top" index (exact)'] - data['left index (exact)']) / self.dpmm
            data['right distance->top (exact) mm'] = abs(data['"top" index (exact)'] - data['right index (exact)']) / self.dpmm
            data['"top"->beam center (exact) mm'] = (data['"top" index (exact)'] - data['beam center index (exact)']) / self.dpmm
            data['"top"->CAX (exact) mm'] = abs(data['"top" index (exact)'] - data['cax index (exact)']) / self.dpmm
        return data

    def inflection_data(self) -> dict:
        """Calculate the profile inflection values using either the 2nd derivative or a fitted Hill function.

        .. note::
            This only applies if the edge detection method is `INFLECTION_...`.

        Parameters
        ----------

        """
        # get max/min of the gradient, which is basically the same as the 2nd deriv 0-crossing
        if self._edge_method == Edge.FWHM:
            raise ValueError("FWHM edge method does not have inflection points. Use a different edge detection method")
        d1 = np.gradient(gaussian_filter1d(self.values, sigma=self._edge_smoothing_ratio * len(self.values)))
        (peak_idxs, _) = MultiProfile(d1).find_peaks(threshold=0.8)
        (valley_idxs, _) = MultiProfile(d1).find_valleys(threshold=0.8)
        left_idx = peak_idxs[0]  # left-most index
        right_idx = valley_idxs[-1]  # right-most index
        if self._edge_method == Edge.INFLECTION_DERIVATIVE:
            data = {'left index (rounded)': left_idx,
                    'left index (exact)': left_idx,
                    'right index (rounded)': right_idx,
                    'right index (exact)': right_idx,
                    'left value (@rounded)': self.values[int(round(left_idx))],
                    'left value (@exact)': self.values[int(round(left_idx))],
                    'right value (@rounded)': self.values[int(round(right_idx))],
                    'right value (@exact)': self.values[int(round(right_idx))]
                    }
            return data
        else:  # Hill
            # the 2nd deriv is a good approximation for the inflection point. Start there and fit Hill about it
            # penum_half_window = self.field_data()['width (exact)'] * self._hill_window_ratio / 2
            penum_half_window = int(round(self._hill_window_ratio * abs(right_idx - left_idx) / 2))

            # left side
            x_data = np.array([x for x in np.arange(left_idx - penum_half_window, left_idx + penum_half_window) if x >=0])
            y_data = self.values[x_data]
            # y_data = self.values[left_idx - penum_half_window: left_idx + penum_half_window]
            left_hill = Hill.fit(x_data, y_data)
            left_infl = left_hill.inflection_idx()

            # right side
            x_data = np.array([x for x in np.arange(right_idx - penum_half_window, right_idx + penum_half_window) if x < len(d1)])
            y_data = self.values[x_data]
            right_hill = Hill.fit(x_data, y_data)
            right_infl = right_hill.inflection_idx()

            data = {'left index (rounded)': left_infl['index (rounded)'],
                    'left index (exact)': left_infl['index (exact)'],
                    'right index (rounded)': right_infl['index (rounded)'],
                    'right index (exact)': right_infl['index (exact)'],
                    'left value (@exact)': left_hill.y(left_infl['index (exact)']),
                    'right value (@exact)': right_hill.y(right_infl['index (exact)']),
                    'left Hill params': left_hill.params, 'right Hill params': right_hill.params}
            return data

    def penumbra(self, lower: int = 20, upper: int = 80):
        """Calculate the penumbra of the field. Dependent on the edge detection method.

        Parameters
        ----------
        lower
            The lower % of the beam to use. If the edge method is FWHM, this is the typical % penumbra you're thinking.
            If the inflection method is used it will be the value/50 of the inflection point value. E.g. if the inflection
            point is perfectly at 50% with a ``lower`` of 20, then the penumbra value here will be 20% of the maximum.
            If the inflection point is at 30% of the max value (say for a FFF beam) then the lower penumbra will be ``lower/50``
            of the inflection point or ``0.3*lower/50``.
        upper
            Upper % of the beam to use. See lower for details.
        """
        if lower > upper:
            raise ValueError("Upper penumbra value must be larger than the lower penumbra value")
        if self._edge_method == Edge.FWHM:
            upper_data = self.fwxm_data(x=upper)
            lower_data = self.fwxm_data(x=lower)
            data = {f'left {lower}% index (exact)': lower_data['left index (exact)'],
                    f'left {lower}% value (@rounded)': lower_data['left value (@rounded)'],
                    f'left {upper}% index (exact)': upper_data['left index (exact)'],
                    f'left {upper}% value (@rounded)': upper_data['left value (@rounded)'],

                    f'right {lower}% index (exact)': lower_data['right index (exact)'],
                    f'right {lower}% value (@rounded)': lower_data['right value (@rounded)'],
                    f'right {upper}% index (exact)': upper_data['right index (exact)'],
                    f'right {upper}% value (@rounded)': upper_data['right value (@rounded)'],

                    'left values': self.values[lower_data['left index (rounded)']:upper_data['left index (rounded)']],
                    'right values': self.values[upper_data['right index (rounded)']:lower_data['right index (rounded)']],

                    f'left penumbra width (exact)': abs(upper_data['left index (exact)'] - lower_data['left index (exact)']),
                    f'right penumbra width (exact)': abs(upper_data['right index (exact)'] - lower_data['right index (exact)']),
                    }
            if self.dpmm:
                data['left penumbra width (exact) mm'] = data['left penumbra width (exact)'] / self.dpmm
                data['right penumbra width (exact) mm'] = data['right penumbra width (exact)'] / self.dpmm
            return data
        elif self._edge_method == Edge.INFLECTION_DERIVATIVE:
            infl_data = self.inflection_data()
            lower_left_value = infl_data['left value (@rounded)']*lower/50*100
            upper_left_value = infl_data['left value (@rounded)']*upper/50*100
            upper_left_data = self.fwxm_data(x=upper_left_value)
            lower_left_data = self.fwxm_data(x=lower_left_value)

            lower_right_value = infl_data['right value (@exact)']*lower/50*100
            upper_right_value = infl_data['right value (@exact)']*upper/50*100
            upper_right_data = self.fwxm_data(x=upper_right_value)
            lower_right_data = self.fwxm_data(x=lower_right_value)

            data = {f'left {lower}% index (exact)': lower_left_data['left index (exact)'],
                    f'left {upper}% index (exact)': upper_left_data['left index (exact)'],

                    f'right {lower}% index (exact)': lower_right_data['right index (exact)'],
                    f'right {upper}% index (exact)': upper_right_data['right index (exact)'],

                    'left values': self.values[lower_left_data['left index (rounded)']:upper_left_data['left index (rounded)']],
                    'right values': self.values[upper_right_data['right index (rounded)']:lower_right_data['right index (rounded)']],

                    f'left penumbra width (exact)': abs(upper_left_data['left index (exact)'] - lower_left_data['left index (exact)']),
                    f'right penumbra width (exact)': abs(upper_right_data['right index (exact)'] - lower_right_data['right index (exact)']),
                    }
            if self.dpmm:
                data['left penumbra width (exact) mm'] = data['left penumbra width (exact)'] / self.dpmm
                data['right penumbra width (exact) mm'] = data['right penumbra width (exact)'] / self.dpmm
            return data
        elif self._edge_method == Edge.INFLECTION_HILL:
            infl_data = self.inflection_data()
            left_hill = Hill.from_params(infl_data['left Hill params'])
            right_hill = Hill.from_params(infl_data['right Hill params'])

            lower_left_value = infl_data['left value (@exact)']*lower/50
            lower_left_index = left_hill.x(lower_left_value)
            upper_left_value = infl_data['left value (@exact)']*upper/50
            upper_left_index = left_hill.x(upper_left_value)

            lower_right_value = infl_data['right value (@exact)']*lower/50
            lower_right_index = right_hill.x(lower_right_value)
            upper_right_value = infl_data['right value (@exact)']*upper/50
            upper_right_index = right_hill.x(upper_right_value)

            data = {f'left {lower}% index (exact)': lower_left_index,
                    f'left {lower}% value (exact)': lower_left_value,
                    f'left {upper}% index (exact)': upper_left_index,
                    f'left {upper}% value (exact)': upper_left_value,

                    f'right {lower}% index (exact)': lower_right_index,
                    f'right {lower}% value (exact)': lower_right_value,
                    f'right {upper}% index (exact)': upper_right_index,
                    f'right {upper}% value (exact)': upper_right_value,

                    'left values': self.values[int(round(lower_left_index)):int(round(upper_left_index))],
                    'right values': self.values[int(round(upper_right_index)):int(round(lower_right_index))],

                    f'left penumbra width (exact)': abs(upper_left_index - lower_left_index),
                    f'right penumbra width (exact)': abs(upper_right_index - lower_right_index),

                    f'left gradient (exact)': left_hill.gradient_at(infl_data['left index (exact)']),
                    r'right gradient (exact)': right_hill.gradient_at(infl_data['right index (exact)']),
                    }
            if self.dpmm:
                data['left penumbra width (exact) mm'] = data['left penumbra width (exact)'] / self.dpmm
                data['left gradient (exact) %/mm'] = data['left gradient (exact)'] * self.dpmm * 100  # 100 to convert to %
                data['right penumbra width (exact) mm'] = data['right penumbra width (exact)'] / self.dpmm
                data['right gradient (exact) %/mm'] = data['right gradient (exact)'] * self.dpmm * 100
            return data

    @argue.options(calculation=('mean', 'median', 'max', 'min', 'area'))
    def field_calculation(self, in_field_ratio: float=0.8, calculation: str='mean') -> Union[float, Tuple[float, float]]:
        """Perform an operation on the field values of the profile.
        This function is useful for determining field symmetry and flatness.

        Parameters
        ----------
        in_field_ratio
            Ratio of the field width to use in the calculation.
        calculation : {'mean', 'median', 'max', 'min', 'area'}
            Calculation to perform on the field values.
        """
        field_values = self.field_data(in_field_ratio)

        if calculation == 'mean':
            return field_values['field values'].mean()
        elif calculation == 'median':
            return float(np.median(field_values['field values']))
        elif calculation == 'max':
            return field_values['field values'].max()
        elif calculation == 'min':
            return field_values['field values'].min()

    def plot(self) -> None:
        """Plot the profile."""
        plt.plot(self.values)
        plt.show() class CollapsedCircleProfile(CircleProfile):
    """A circular profile that samples a thick band around the nominal circle, rather than just a 1-pixel-wide profile
    to give a mean value.
    """
    width_ratio: float
    num_profiles: int

    @argue.bounds(width_ratio=(0, 1))
    def __init__(self, center: Point, radius: NumberLike, image_array: Union[np.ndarray, 'ArrayImage'], start_angle: int=0,
                 ccw: bool=True, sampling_ratio: float=1.0, width_ratio: float=0.1, num_profiles: int=20):
        """
        Parameters
        ----------
        width_ratio : float
            The "thickness" of the band to sample. The ratio is relative to the radius. E.g. if the radius is 20
            and the width_ratio is 0.2, the "thickness" will be 4 pixels.
        num_profiles : int
            The number of profiles to sample in the band. Profiles are distributed evenly within the band.

        See Also
        --------
        :class:`~pylinac.core.profile.CircleProfile` : Further parameter info.
        """
        self.width_ratio = width_ratio
        self.num_profiles = num_profiles
        super().__init__(center, radius, image_array, start_angle, ccw, sampling_ratio)

    @property
    def _radii(self) -> np.ndarray:
        return np.linspace(start=self.radius * (1 - self.width_ratio), stop=self.radius * (1 + self.width_ratio),
                  num=self.num_profiles)

    @property
    def size(self) -> float:
        return np.pi * max(self._radii) * 2 * self.sampling_ratio

    @property
    def _multi_x_locations(self) -> List:
        """List of x-locations of the sampling profiles"""
        x = []
        cos = np.cos(self._radians)
        # extract profile for each circle radii
        for radius in self._radii:
            x.append(cos * radius + self.center.x)
        return x

    @property
    def _multi_y_locations(self) -> List:
        """List of x-locations of the sampling profiles"""
        y = []
        sin = np.sin(self._radians)
        # extract profile for each circle radii
        for radius in self._radii:
            y.append(sin * radius + self.center.y)
        return y

    @property
    def _profile(self) -> np.ndarray:
        """The actual profile array; private attr that is passed to MultiProfile."""
        profile = np.zeros(len(self._multi_x_locations[0]))
        for radius, x, y in zip(self._radii, self._multi_x_locations, self._multi_y_locations):
            profile += ndimage.map_coordinates(self.image_array, [y, x], order=0)
        profile /= self.num_profiles
        return profile

    def plot2axes(self, axes: plt.Axes=None, edgecolor: str='black', fill: bool=False, plot_peaks: bool=True) -> None:
        """Add 2 circles to the axes: one at the maximum and minimum radius of the ROI.

        See Also
        --------
        :meth:`~pylinac.core.profile.CircleProfile.plot2axes` : Further parameter info.
        """
        if axes is None:
            fig, axes = plt.subplots()
            axes.imshow(self.image_array)
        axes.add_patch(mpl_Circle((self.center.x, self.center.y), edgecolor=edgecolor, radius=self.radius*(1+self.width_ratio),
                                  fill=fill))
        axes.add_patch(mpl_Circle((self.center.x, self.center.y), edgecolor=edgecolor, radius=self.radius*(1-self.width_ratio),
                                  fill=fill))
        if plot_peaks:
            x_locs = [peak.x for peak in self.peaks]
            y_locs = [peak.y for peak in self.peaks]
            axes.autoscale(enable=False)
            axes.scatter(x_locs, y_locs, s=20, marker='x', c=edgecolor) class Interpolation(enum.Enum):
    """Interpolation Enum"""
    NONE = None  #:
    LINEAR = 'Linear'  #:
    SPLINE = 'Spline'  #: def open_path(path: str) -> None:
    """Open the specified path in the system default viewer."""

    if os.name == 'darwin':
        launcher = "open"
    elif os.name == 'posix':
        launcher = "xdg-open"
    elif os.name == 'nt':
        launcher = "explorer"
    subprocess.call([launcher, path]) class ResultBase:
    pylinac_version: str = field(init=False)  #:
    date_of_analysis: datetime = field(init=False)  #:

    def __post_init__(self):
        self.pylinac_version = __version__
        self.date_of_analysis = datetime.today() def get_dicom_cmap():
    """Return the DICOM colormap. Passed to cmap parameter in matplotlib calls."""
    return DICOM_COLORMAP
import copy
import dataclasses
import io
import argue
import matplotlib.pyplot as plt
import numpy as np
from dataclasses import dataclass
from typing import Union, List, Optional, Tuple, BinaryIO
from scipy import optimize
from typing.io import IO
from .core import image, pdf
from .core.geometry import Point, Line, Circle
from .core.io import get_url, TemporaryZipDirectory, retrieve_demo_file
from .core.profile import SingleProfile, CollapsedCircleProfile, Interpolation
from .core.utilities import open_path, ResultBase
from .settings import get_dicom_cmap
* **Accurate, FWHM star line detection** - Pylinac uses not simply the maximum value to find the center of a star line,
  but analyzes the entire star profile to determine the center of the FWHM, ensuring small noise or maximum value bias is avoided.
* **Adaptive searching** - If you passed pylinac a set of parameters and a good result wasn't found, pylinac can recover and
  do an adaptive search by adjusting parameters to find a "reasonable" wobble.
"""




@dataclass
class StarshotResults(ResultBase):
    """This class should not be called directly. It is returned by the ``results_data()`` method.
    It is a dataclass under the hood and thus comes with all the dunder magic.

    Use the following attributes as normal class attributes.
    """
    tolerance_mm: float  #:
    circle_diameter_mm: float  #:
    circle_radius_mm: float  #:
    passed: bool  #:
    circle_center_x_y: Tuple[float, float]  #:


class Starshot:
    """Class that can determine the wobble in a "starshot" image, be it gantry, collimator,
    couch or MLC. The image can be a scanned film (TIF, JPG, etc) or a sequence of EPID DICOM images.

    Attributes
    ----------
    image : :class:`~pylinac.core.image.Image`
    circle_profile : :class:`~pylinac.starshot.StarProfile`
    lines : :class:`~pylinac.starshot.LineManager`
    wobble : :class:`~pylinac.starshot.Wobble`
    tolerance : :class:`~pylinac.starshot.Tolerance`

    Examples
    --------
    Run the demo:
        >>> Starshot.run_demo()

    Typical session:
        >>> img_path = r"C:/QA/Starshots/Coll.jpeg"
        >>> mystar = Starshot(img_path, dpi=105, sid=1000)
        >>> mystar.analyze()
        >>> print(mystar.results())
        >>> mystar.plot_analyzed_image()
    """
    def __init__(self, filepath: Union[str, BinaryIO], **kwargs):
        """
        Parameters
        ----------
        filepath
            The path to the image file.
        kwargs
            Passed to :func:`~pylinac.core.image.load`.
        """
        self.image = image.load(filepath, **kwargs)
        self.wobble = Wobble()
        self.tolerance = 1
        if self.image.dpmm is None:
            raise ValueError("DPI was not a tag in the image nor was it passed in. Please pass a DPI value")
        if self.image.sid is None:
            raise ValueError("Source-to-Image distance was not an image tag and was not passed in. Please pass an SID value.")

    @classmethod
    def from_url(cls, url: str, **kwargs):
        """Instantiate from a URL.

        Parameters
        ----------
        url : str
            URL of the raw file.
        kwargs
            Passed to :func:`~pylinac.core.image.load`.
        """
        filename = get_url(url)
        return cls(filename, **kwargs)

    @classmethod
    def from_demo_image(cls):
        """Construct a Starshot instance and load the demo image."""
        demo_file = retrieve_demo_file(url='starshot.tif')
        return cls(demo_file, sid=1000)

    @classmethod
    def from_multiple_images(cls, filepath_list: list, **kwargs):
        """Construct a Starshot instance and load in and combine multiple images.

        Parameters
        ----------
        filepath_list : iterable
            An iterable of file paths to starshot images that are to be superimposed.
        kwargs
            Passed to :func:`~pylinac.core.image.load_multiples`.
        """
        obj = cls.from_demo_image()
        obj.image = image.load_multiples(filepath_list, **kwargs)
        return obj

    @classmethod
    def from_zip(cls, zip_file: str, **kwargs):
        """Construct a Starshot instance from a ZIP archive.

        Parameters
        ----------
        zip_file : str
            Points to the ZIP archive. Can contain a single or multiple images. If multiple images
            the images are combined and thus should be from the same test sequence.
        kwargs
            Passed to :func:`~pylinac.core.image.load_multiples`.
        """
        with TemporaryZipDirectory(zip_file) as tmpdir:
            image_files = image.retrieve_image_files(tmpdir)
            if not image_files:
                raise IndexError(f"No valid starshot images were found in {zip_file}")
            if len(image_files) > 1:
                return cls.from_multiple_images(image_files, **kwargs)
            else:
                return cls(image_files[0], **kwargs)


    def _get_reasonable_start_point(self) -> Point: