prof_funct.py

# Copyright (C) 2019 Paul King

# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version (the "AGPL-3.0+").

# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License and the additional terms for more
# details.

# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.

# ADDITIONAL TERMS are also included as allowed by Section 7 of the GNU
# Affero General Public License. These additional terms are Sections 1, 5,
# 6, 7, 8, and 9 from the Apache License, Version 2.0 (the "Apache-2.0")
# where all references to the definition "License" are instead defined to
# mean the AGPL-3.0+.

# You should have received a copy of the Apache-2.0 along with this
# program. If not, see <http://www.apache.org/licenses/LICENSE-2.0>.

""" For importing, analyzing, and comparing dose or intensity profiles
    from different sources."""

import os
import copy
import sys

from typing import Callable
from scipy import interpolate

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg

import tkinter as tk
from tkinter.filedialog import askopenfilename
from matplotlib.backends.backend_tkagg import (
    FigureCanvasTkAgg, NavigationToolbar2Tk)
from matplotlib.backend_bases import key_press_handler
from matplotlib.figure import Figure

from functools import partial
import PIL
import csv
import re
import time
import pwlf

# NumpyFunction = Callable[[np.ndarray], np.ndarray]

# pylint: disable = C0103, C0121, W0102

class Profile():
    """  One-dimensional distribution of intensity vs position.

    Attributes
    ----------
    x : np.array
        position, +/- in cm
    y : np.array
        intensity in unspecified units
    meta : dict, optional
        metadata

    Notes
    -----
    Requires PIL.  https://pypi.org/project/PIL/

    """

    def __init__(self, x=np.array([]),
                 y=np.array([]), meta={}):
        """ create profile

        Parameters
        ----------
        x : np.array, optional
        y : np.array, optional
        meta : dict, optional

        Notes
        -----
        Normally created empty, then filled using a method, which returns
        a new Profile.

        """
        self.x = np.array(x)
        self.y = np.array(y)
        self.meta = meta
        if len(self.x) < 2:
            self.interp = None
        else:
            self.interp = interpolate.interp1d(self.x, self.y,
                                               bounds_error=False, fill_value=0.0)

    def __len__(self):
        """ # data points  """
        return len(self.x)

    def __eq__(self, other):  # SAME DATA POINTS
        """ same data points """
        if np.array_equal(self.x, other.x) and \
           np.array_equal(self.y, other.y) and \
           self.meta == other.meta:
            return True
        else:
            return False

    def __copy__(self):
        """ deep copy """
        return copy.deepcopy(self)

    def __str__(self):
        """
        Examples
        --------
        ``Profile object: 83 pts | x (-16.4 cm -> 16.4 cm) | y (0.22 -> 45.54)``

        """
        try:
            fmt_str = 'Profile object: '
            fmt_str += '{} pts | x ({} cm -> {} cm) | y ({} -> {})'
            return fmt_str.format(len(self.x),
                                  min(self.x), max(self.x),
                                  min(self.y), max(self.y))
        except ValueError:
            return ''  # EMPTY PROFILE

    def __add__(self, other):
        """ shift right """
        new_x = self.x + other
        return Profile(x=new_x, y=self.y, meta=self.meta)
    __radd__ = __add__
    __iadd__ = __add__

    def __sub__(self, other):
        """ shift left """
        self.x -= other
        return Profile(x=self.x, y=self.y, meta=self.meta)
    __rsub__ = __sub__
    __isub__ = __sub__

    def __mul__(self, other):
        """ scale y """
        self.y *= other
        return self
    __rmul__ = __mul__
    __imul__ = __mul__


    def get_y(self, x):
        """ y-value at distance x

        Return a y value based on interpolation of source data for a
        supplied distance.

        Parameters
        ----------
        x : float

        Returns
        -------
        float

        """
        try:
            return self.interp(x)
        except ValueError:
            return np.nan

    def get_x(self, y):
        """ tuple of x-values at intensity y

        Return distance values based on interpolation of source data for a
        supplied y value.

        Parameters
        ----------
        y : float

        Returns
        -------
        tuple : (x1, x2, ...)

         """

        dose_step = (max(self.y)-min(self.y)) / 100
        x_ = self.resample_y(dose_step).x
        y_ = self.resample_y(dose_step).y
        dists = []
        for i in range(1, len(x_)):
            val = None
            if (y_[i]-y)*(y_[i-1]-y) < 0:
                val = (x_[i]-((y_[i]-y)/(y_[i]-y_[i-1]))*(x_[i]-x_[i-1]))
            elif np.isclose(y_[i], y):
                val = x_[i]
            if val and (val not in dists):
                dists.append(val)
        return tuple(dists)

    def get_increment(self):
        """ minimum step-size increment

        Returns
        -------
        increment : float

        """
        steps = np.diff(self.x)
        if np.isclose(steps.min(), steps.mean()):
            return steps.mean()
        else:
            return steps.min()

    def plot(self, marker='o-'):
        """ profile plot

        Parameters
        ----------
        marker : string, optional

        Returns
        -------
        None

        """
        plt.plot(self.x, self.y, marker)
        plt.show()
        return

    def slice_segment(self, start=-np.inf, stop=np.inf):
        """ slice between given end-points

        Resulting profile is comprised of those points in the source
        profile whose distance values are not-less-than start and
        not-greater-than stop.

        Parameters
        ----------
        start : float, optional
        stop : float, optional

        Returns
        -------
        Profile

        """
        try:
            start = max(start, min(self.x))  # default & limit to curve ends
            stop = min(stop, max(self.x))
            new_x = self.x[np.logical_and(self.x >= start, self.x <= stop)]
            new_y = self.interp(new_x)
        except ValueError:
            new_x = []
            new_y = []

        return Profile(new_x, new_y)


    def resample_x(self, step=None, begin=None, end=None, num_points=None):
        """ resampled x-values at a given increment

        Resulting profile has stepsize of the indicated step based on
        linear interpolation over the points of the source profile.

        Parameters
        ----------
        step : float, optional
        begin : float, optional
        end : float, optional
        num_points : int, optional

        Returns
        -------
        Profile

        """

        if not begin:
            begin = min(self.x)
        if not end:
            end = max(self.x)
        if not step:
            try: 
                step = (end-begin)/num_points
            except TypeError:
                step = np.average(np.diff(self.x))

        new_x = np.arange(begin, end, step)
        new_y = self.get_y(new_x)
        return Profile(new_x, new_y, self.meta)


    def resample_y(self, step):
        """ resampled y-values at a given increment

        Resulting profile has nonuniform step-size, but each step
        represents and approximately equal step in dose.

        Parameters
        ----------
        step : float
            sampling increment

        Returns
        -------
        Profile

        """

        temp_x = np.arange(min(self.x), max(self.x),
                           0.01*self.get_increment())
        temp_y = self.interp(temp_x)

        resamp_x = [temp_x[0]]
        resamp_y = [temp_y[0]]

        last_y = temp_y[0]

        for i, _ in enumerate(temp_x):
            if np.abs(temp_y[i] - last_y) >= step:
                resamp_x.append(temp_x[i])
                resamp_y.append(temp_y[i])
                last_y = temp_y[i]

        if temp_x[-1] not in resamp_x:
            resamp_x.append(temp_x[-1])
            resamp_y.append(temp_y[-1])

        return Profile(x=np.array(resamp_x), y=np.array(resamp_y), meta=self.meta)

    def make_normal_y(self, x=0.0, y=1.0):
        """ normalised to dose at distance

        Source profile values multiplied by scaling factor to yield the specified dose at
        the specified distance. If distance is not specified, the central axis value is
        used. If dose is not specified, then normalization is to unity. With neither
        specified, resulting curve is the conventional off-center-ratio.

        Parameters
        ----------
        x : float, optional
        y : float, optional

        Returns
        -------
        Profile

        """
        norm_factor = y / self.get_y(x)
        new_x = self.x
        new_y = norm_factor * self.y
        return Profile(new_x, new_y, meta=self.meta)

    def get_edges(self):
        """ x-values of profile edges (left, right)

        Notes
        -----
        Points of greatest positive and greatest negative gradient.

        Returns
        -------
        tuple

        """
        dydx = list(np.gradient(self.y, self.x))
        lt_edge = self.x[dydx.index(max(dydx))]
        rt_edge = self.x[dydx.index(min(dydx))]
        return (lt_edge, rt_edge)

    def make_normal_x(self):
        """ normalised to distance at edges

        Source profile distances multiplied by scaling factor to yield unit distance
        at beam edges. [1]_ [2]_

        Returns
        -------
        Profile


        References
        ----------
        .. [1] Milan & Bentley, BJR Feb-74, The Storage and manipulationof radiation dose data
           in a small digital computer
        .. [2] Heintz, King, & Childs, May-95, User Manual, Prowess 3000 CT Treatment Planning


        """

        lt_edge, rt_edge = self.get_edges()
        cax = 0.5*(lt_edge + rt_edge)

        new_x = []
        for _, dist in enumerate(self.x):
            if dist < cax:
                new_x.append(-dist/lt_edge)
            elif dist > cax:
                new_x.append(dist/rt_edge)
            else:
                new_x.append(0.0)

        return Profile(new_x, self.y, meta=self.meta)

    def slice_umbra(self):
        """ umbra central 80%

        Source dose profile sliced to include only the central region between beam edges.

        Returns
        -------
        Profile

        """
        lt, rt = self.get_edges()
        idx = [i for i, d in enumerate(
            self.x) if d >= 0.8 * lt and d <= 0.8 * rt]
        new_x = self.x[idx[0]:idx[-1]+1]
        new_y = self.y[idx[0]:idx[-1]+1]

        return Profile(x=new_x, y=new_y, meta=self.meta)

    def slice_penumbra(self):
        """ penumbra (20 -> 80%, 80 -> 20%)

        Source dose profile sliced to include only the penumbral edges, where the dose
        transitions from 20% - 80% of the umbra dose, as precent at the umbra edge,
        to support wedged profiles.

        Returns
        -------
        tuple
            (left penumbra Profile, right penumbra Profile)

        """

        not_umbra = {'lt': self.slice_segment(stop=self.slice_umbra().x[0]),
                     'rt': self.slice_segment(start=self.slice_umbra().x[-1])}
        result = []
        for side in not_umbra:
            min_val = min(not_umbra[side].y)
            max_val = max(not_umbra[side].y)
            incr_val = 0.2 * (max_val - min_val)
            lo_x = not_umbra[side].get_x(min_val + incr_val)
            hi_x = not_umbra[side].get_x(max_val - incr_val)
            coords = [lo_x, hi_x]
            coords.sort()
            penum = not_umbra[side].slice_segment(start=coords[0], stop=coords[1])
            result.append(penum)
        return tuple(result)

    def slice_shoulders(self):
        """ shoulders (penumbra -> umbra, umbra -> penumbra)

        Source dose profile sliced to include only the profile shoulders,
        outside the central 80% of of the profile but inside the region bounded
        by the 20-80% transition.

        Returns
        -------
        tuple
            (left shoulder Profile, right shoulder Profile)

        """
        try:
            lt_start = self.slice_penumbra()[0].x[-1]
        except IndexError:
            lt_start = self.slice_umbra().x[-1]
        lt_stop = self.slice_umbra().x[0]

        rt_start = self.slice_umbra().x[-1]
        try:
            rt_stop = self.slice_penumbra()[-1].x[0]
        except IndexError:
            rt_stop = self.slice_umbra().x[0]

        lt_should = self.slice_segment(start=lt_start, stop=lt_stop)
        rt_should = self.slice_segment(start=rt_start, stop=rt_stop)

        return (lt_should, rt_should)

    def slice_tails(self):
        """ tails (-> penumbra, penumbra ->)

        Source dose profile sliced to include only the profile tail,
        outside the beam penumbra.

        Returns
        -------
        tuple
            (left tail Profile, right tail Profile)

        """
        lt_start = self.x[0]
        try:
            lt_stop = self.slice_penumbra()[0].x[0]
        except IndexError:
            lt_stop = self.slice_shoulders()[0].x[0]

        try:
            rt_start = self.slice_penumbra()[-1].x[-1]
        except IndexError:
            rt_start = self.slice_shoulders()[-1].x[-1]
        rt_stop = self.x[-1]

        lt_tail = self.slice_segment(start=lt_start, stop=lt_stop)
        rt_tail = self.slice_segment(start=rt_start, stop=rt_stop)

        return (lt_tail, rt_tail)

    def get_flatness(self):
        """ dose range relative to mean

        Calculated as the dose range normalized to mean dose.

        Returns
        -------
        float

        """
        dose = self.slice_umbra().y
        return (max(dose)-min(dose))/np.average(dose)

    def get_symmetry(self):
        """ max point diff relative to mean

        Calculated as the maximum difference between corresponding points
        on opposite sides of the profile center, relative to mean dose.

        Returns
        -------
        float

        """
        dose = self.slice_umbra().y
        return max(np.abs(np.subtract(dose, dose[::-1])/np.average(dose)))

    def make_symmetric(self):
        """ avg of corresponding points

        Created by averaging over corresponding +/- distances,
        except at the endpoints.

        Returns
        -------
        Profile

        """

        reflected = Profile(x=-self.x[::-1], y=self.y[::-1])

        step = self.get_increment()
        new_x = np.arange(min(self.x), max(self.x), step)
        new_y = [self.y[0]]
        for n in new_x[1:-1]:  # AVOID EXTRAPOLATION
            new_y.append(0.5*self.interp(n) + 0.5*reflected.interp(n))
        new_y.append(reflected.y[0])

        return Profile(x=new_x, y=new_y, meta=self.meta)

    def make_centered(self):
        """ shift to align edges

        Created by shifting the profile based on edge locations.

        Returns
        -------
        Profile

        """

        return self - np.average(self.get_edges())

    def make_flipped(self):
        """ flip L -> R

        Created by reversing the sequence of y values.

        Returns
        -------
        Profile

        """

        return Profile(x=self.x, y=self.y[::-1], meta=self.meta)

    def align_to(self, other):
        """ shift self to align to other

        Calculated using shift that produces greatest peak correlation between
        the curves. Flips the curve left-to-right, if this creates a better fit.

        Parameters
        ----------
        other : Profile
            profile to be be shifted to

        Returns
        -------
        Profile

        """

        dist_step = min(self.get_increment(), other.get_increment())

        dist_vals_fixed = np.arange(
            -3*abs(min(list(self.x) + list(other.x))),
            3*abs(max(list(self.x) + list(other.x))),
            dist_step)
        dose_vals_fixed = other.interp(dist_vals_fixed)
        fixed = Profile(x=dist_vals_fixed, y=dose_vals_fixed)

        possible_offsets = np.arange(
            max(min(self.x), min(other.x)),
            min(max(other.x), max(self.x)) + dist_step,
            dist_step)

        best_fit_qual, best_offset, flipped = 0, -np.inf, False
        for offset in possible_offsets:
            fit_qual_norm = max(np.correlate(
                dose_vals_fixed,
                (self + offset).interp(fixed.x)))
            fit_qual_flip = max(np.correlate(
                dose_vals_fixed,
                (self.make_flipped() + offset).interp(fixed.x)))

            if fit_qual_norm > best_fit_qual:
                best_fit_qual = fit_qual_norm
                best_offset = offset
                flipped = False
            if fit_qual_flip > best_fit_qual:
                best_fit_qual = fit_qual_flip
                best_offset = offset
                flipped = True

        if flipped:
            return self.make_flipped() + best_offset
        else:
            return self + best_offset

if __name__ == "__main__":
    import prof_gui
    root = tk.Tk()
    prof_gui.GUI(root).pack(side="top", fill="both", expand=True)
    root.mainloop()