# -*- coding: utf-8 -*-

"""

Created on Tue Nov 08 17:45:23 2016

@author: Shreyas_V

"""

# Copyright (C) 2013 Oskar Maier

# 

# This program is free software: you can redistribute it and/or modify

# it under the terms of the GNU General Public License as published by

# the Free Software Foundation, either version 3 of the License, or

# (at your option) any later version.

# 

# 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 General Public License for more details.

# 

# You should have received a copy of the GNU General Public License

# along with this program.  If not, see <http://www.gnu.org/licenses/>.

#

# author Oskar Maier

# version r0.1.2

# since 2013-09-04

# status Release

# build-in modules

# third-party modules

import numpy

from scipy.interpolate.interpolate import interp1d

# path changes

# own modules

# code

class IntensityRangeStandardization (object):

    r"""

    Class to standardize intensity ranges between a number of images.

    **Short description:**

    Often images containing similar objects or scenes have different intensity ranges

    that make it difficult to compare them manually as well as to process them

    further.

    IntensityRangeStandardization offers a way to transform a number of such images

    intensity ranges to a common standard intensity space without any loss of

    information using a multi-segment linear transformation model.

    Once learned, this model can be applied to other, formerly unseen images to map

    them to the same standard intensity space.    

    **Concept of similar images:**

    IntensityRangeStandardization is limited to similar images. Images containing

    different object or different compositions of objects are not suitable to be

    transformed to a common intensity space (and it would furthermore not make much

    sense).

    A typical application of IntensityRangeStandardization are MRI images showing the

    same body region. These often have different intensity ranges, even when acquired

    from the same patient and using the same scanner. For further processing, e.g.

    for training a classifier, they have to be mapped to a common intensity space.

    **Failure of the transformation:**

    The method implemented in IntensityRangeStandardization ensures that no

    information is lost i.e. a lossless transformation is performed. This can be

    assured when there exists a one-to-one mapping between the images original

    intensity values and their values mapped to the standard intensity space.

    But since the transformation model is trained on, and the standard intensity

    space range selected over the training images, this can not be guaranteed for all

    formerly unseen image. If they differ greatly from the training set images, a

    lossless transformation can not be assured anymore. In this case the transform()

    method will throw an InformationLossException.

    Should this happen, the model needs to be re-trained with the original training

    images and additionally the images which caused the failure. Since this will lead

    to a new intensity standard space, all already transformed images have to be

    processed again.

    **Setting the training parameters:**

    The method comes with a set of default parameters, that are suitable for most

    cases. But for some special cases, it might be better to set them on your own. Ti

    understand the working of the parameters, it is recommended to read the detailed

    method description first.

    **The method depends on three parameters:**

    cutoffp, i.e. the cut-off percentiles

        These are used to the define the intensity outliers, both during training and

        image transformation. The default values are usualy a good choice.

        (in [1]_ these are called the minimum and maximum percentile values pc1 and pc2 respectively)

    landmarkp, i.e. the landmark percentiles

        These percentiles define the landmark positions. The more supplied, the more

        exact but less general becomes the model. It is common to supply equally

        spaced percentiles between 0 and 100.

        (in [1]_ these are called the landmark locations mu_1, .., mu_l)

    strange, i.e. the standard intensity space range

        These two intensity values define roughly the standard intensity space (or

        common intensity space of the images; or even target intensity space) to

        which each images intensities are mapped. This space can be supplied, but it

        is usually recommended to let the method select it automatically during the

        training process. It is additionally possible to supply only the lower or

        upper range border and set the other to ''auto'', in which case the method

        chooses the range automatically, but not the position. 

        (in [1]_ these are called the minimum and maximum intensities on the standard scale of the IOI s1 resp. s2)

    **Details of the method:**

    In the following the method is described in some more detail. For even more

    information see [1]_.

    Essentially the method is based on a multi-segment linear transformation model. A

    standard intensity space (or common intensity space) is defined by an intensity

    value range ''stdrange''.

    During the training phase, the intensity values at certain cut-off percentiles of

    each image are computed and a single-segment linear mapping from them to the

    standard intensity space range limits created. Then the images intensity values

    at a number of landmark percentiles are extracted and passed to the linear

    mapping to be transfered roughly to the standard intensity space. The mean of all

    these mapped landmark intensities form the model learned.

    When presented with an image to transform, these images intensity values are

    extracted at the cut-off percentile as well as at the landmark percentile

    positions. This results in a number of segments. Using these and the

    corresponding standard intensity space range values and learned mean landmark

    values, a multi-segment linear transformation model is created for the image.

    This is then applied to the images intensity values to map them to the standard

    intensity space.

    Outliers, i.e. the images intensity values that lie outside of the cut-off

    percentiles, are treated separately. They are transformed like the first resp.

    last segmented of the transformation model. Not that this means the transformed

    images intensity values do not always lie inside the standard intensity space

    range, but are fitted as best as possible inside.

    Parameters

    ----------

    cutoffp : (float, float)

        Lower and upper cut-off percentiles to exclude outliers.

    landmarkp : sequence of floats

        List of percentiles serving as model landmarks, must lie

        between the cutoffp values.

    stdrange : string or (float, float)

        The range of the standard intensity space for which a

        transformation is learned; when set to 'auto, automatically

        determined from the training image upon training; it is also

        possible to fix either the upper or the lower border value and

        setting the other to 'auto'.

    Examples

    --------

    We have a number of similar images with varying intensity ranges. To make them

    comparable, we would like to transform them to a common intensity space. Thus we

    run:

        >>> from medpy.filter import IntensityRangeStandardization

        >>> irs = IntensityRangeStandardization()

        >>> trained_model, transformed_images = irs.train_transform(images)

    Let us assume we now obtain another, new image, that we would like to make

    comparable to the others. As long as it does not differ to much from these, we

    can simply call:

        >>> transformed_image = irs.transform(new_image)

    For many application, not all images are already available at the time of

    execution. It would therefore be good to be able to preserve a once trained

    model. The solution is to just pickle the once trained model:

        >>> import pickle

        >>> with open('my_trained_model.pkl', 'wb') as f:

        >>>     pickle.dump(irs, f)

    And load it again when required with:

        >>> with open('my_trained_model.pkl', 'r') as f:

        >>>     irs = pickle.load(f)

    References

    ----------

    .. [1] Nyul, L.G.; Udupa, J.K.; Xuan Zhang, "New variants of a method of MRI scale

       standardization," Medical Imaging, IEEE Transactions on , vol.19, no.2, pp.143-150,

       Feb. 2000

    """    

    # static member variables

    L2 = [50]

    """1-value landmark points model."""

    L3 = [25, 50, 75]

    """3-value landmark points model."""

    L4 = [10, 20, 30, 40, 50, 60, 70, 80, 90]

    """9-value landmark points model."""

    def __init__(self, cutoffp = (1, 99), landmarkp = L4, stdrange = 'auto'):

        # check parameters

        if not IntensityRangeStandardization.is_sequence(cutoffp):

            raise ValueError('cutoffp must be a sequence')

        if not 2 == len(cutoffp):

            raise ValueError('cutoffp must be of length 2, not {}'.format(len(cutoffp)))

        if not IntensityRangeStandardization.are_numbers(cutoffp):

            raise ValueError('cutoffp elements must be numbers')

        if not IntensityRangeStandardization.are_in_interval(cutoffp, 0, 100, 'included'):

            raise ValueError('cutoffp elements must be in [0, 100]')

        if not cutoffp[1] > cutoffp[0]:

            raise ValueError('the second element of cutoffp must be larger than the first')

        if not IntensityRangeStandardization.is_sequence(landmarkp):

            raise ValueError('landmarkp must be a sequence')

        if not 1 <= len(landmarkp):

            raise ValueError('landmarkp must be of length >= 1, not {}'.format(len(landmarkp)))

        if not IntensityRangeStandardization.are_numbers(landmarkp):

            raise ValueError('landmarkp elements must be numbers')

        if not IntensityRangeStandardization.are_in_interval(landmarkp, 0, 100, 'included'):

            raise ValueError('landmarkp elements must be in [0, 100]')

        if not IntensityRangeStandardization.are_in_interval(landmarkp, cutoffp[0], cutoffp[1], 'excluded'):

            raise ValueError('landmarkp elements must be in between the elements of cutoffp')

        if not len(landmarkp) == len(numpy.unique(landmarkp)):

            raise ValueError('landmarkp elements must be unique')

        if 'auto' == stdrange:

            stdrange = ('auto', 'auto')

        else:

            if not IntensityRangeStandardization.is_sequence(stdrange):

                raise ValueError('stdrange must be a sequence or \'auto\'')

            if not 2 == len(stdrange):

                raise ValueError('stdrange must be of length 2, not {}'.format(len(stdrange)))

            if not 'auto' in stdrange:

                if not IntensityRangeStandardization.are_numbers(stdrange):

                    raise ValueError('stdrange elements must be numbers or \'auto\'')

                if not stdrange[1] > stdrange[0]:

                    raise ValueError('the second element of stdrange must be larger than the first')

            elif 'auto' == stdrange[0] and not IntensityRangeStandardization.is_number(stdrange[1]):

                raise ValueError('stdrange elements must be numbers or \'auto\'')

            elif 'auto' == stdrange[1] and not IntensityRangeStandardization.is_number(stdrange[0]):

                raise ValueError('stdrange elements must be numbers or \'auto\'')

        # process parameters

        self.__cutoffp = IntensityRangeStandardization.to_float(cutoffp)

        self.__landmarkp = IntensityRangeStandardization.to_float(sorted(landmarkp))

        self.__stdrange = ['auto' if 'auto' == x else float(x) for x in stdrange]

        # initialize remaining instance parameters

        self.__model = None

        self.__sc_umins = None

        self.__sc_umaxs = None

    def train(self, images):

        r"""

        Train a standard intensity space and an associated transformation model.

        Note that the passed images should be masked to contain only the foreground.

        Parameters

        ----------

        images : sequence of array_likes

            A number of images.

        Returns

        -------

        IntensityRangeStandardization : IntensityRangeStandardization

            This instance of IntensityRangeStandardization

        """

        self.__stdrange = self.__compute_stdrange(images)

        lim = []

        for idx, i in enumerate(images):

            ci = numpy.array(numpy.percentile(i, self.__cutoffp))

            li = numpy.array(numpy.percentile(i, self.__landmarkp))

            ipf = interp1d(ci, self.__stdrange)

            lim.append(ipf(li))

            # treat single intensity accumulation error            

            if not len(numpy.unique(numpy.concatenate((ci, li)))) == len(ci) + len(li):

                raise SingleIntensityAccumulationError('Image no.{} shows an unusual single-intensity accumulation that leads to a situation where two percentile values are equal. This situation is usually caused, when the background has not been removed from the image. Another possibility would be to reduce the number of landmark percentiles landmarkp or to change their distribution.'.format(idx))

        self.__model = [self.__stdrange[0]] + list(numpy.mean(lim, 0)) + [self.__stdrange[1]]

        self.__sc_umins = [self.__stdrange[0]] + list(numpy.min(lim, 0)) + [self.__stdrange[1]]

        self.__sc_umaxs = [self.__stdrange[0]] + list(numpy.max(lim, 0)) + [self.__stdrange[1]]

        return self

    def transform(self, image, surpress_mapping_check = False):

        r"""

        Transform an images intensity values to the learned standard intensity space.

        Note that the passed image should be masked to contain only the foreground.

        The transformation is guaranteed to be lossless i.e. a one-to-one mapping between

        old and new intensity values exists. In cases where this does not hold, an error

        is thrown. This can be suppressed by setting ``surpress_mapping_check`` to 'True'.

        Do this only if you know what you are doing.

        Parameters

        ----------

        image : array_like

            The image to transform.

        surpress_mapping_check : bool

            Whether to ensure a lossless transformation or not.

        Returns

        -------

        image : ndarray

            The transformed image

        Raises

        -------

        InformationLossException 

            If a lossless transformation can not be ensured

        Exception

            If no model has been trained before

        """

        if None == self.__model:

            raise UntrainedException('Model not trained. Call train() first.')

        image = numpy.asarray(image)

        # determine image intensity values at cut-off percentiles & landmark percentiles

        li = numpy.percentile(image, [self.__cutoffp[0]] + self.__landmarkp + [self.__cutoffp[1]])

        # treat single intensity accumulation error            

        if not len(numpy.unique(li)) == len(li):

            raise SingleIntensityAccumulationError('The image shows an unusual single-intensity accumulation that leads to a situation where two percentile values are equal. This situation is usually caused, when the background has not been removed from the image. The only other possibility would be to re-train the model with a reduced number of landmark percentiles landmarkp or a changed distribution.')

        # create linear mapping models for the percentile segments to the learned standard intensity space  

        ipf = interp1d(li, self.__model, bounds_error = False)

        # transform the input image intensity values

        output = ipf(image)

        # treat image intensity values outside of the cut-off percentiles range separately

        llm = IntensityRangeStandardization.linear_model(li[:2], self.__model[:2])

        rlm = IntensityRangeStandardization.linear_model(li[-2:], self.__model[-2:])

        output[image < li[0]] = llm(image[image < li[0]])

        output[image > li[-1]] = rlm(image[image > li[-1]])

        if not surpress_mapping_check and not self.__check_mapping(li):

            raise InformationLossException('Image can not be transformed to the learned standard intensity space without loss of information. Please re-train.')

        return output

    def train_transform(self, images, surpress_mapping_check = False):

        r"""

        See also

        --------

        train, transform

        """

        ret = self.train(images)

        outputs = [self.transform(i, surpress_mapping_check) for i in images]

        return ret, outputs

    @property

    def stdrange(self):

        """Get the set resp. learned standard intensity range."""

        return self.__stdrange

    @property

    def cutoffp(self):

        """Get the cut-off percentiles."""

        return self.__cutoffp

    @property

    def landmarkp(self):

        """Get the landmark percentiles."""

        return self.__landmarkp

    @property

    def model(self):

        """Get the model (the learned percentile values)."""

        return self.__model

    def __compute_stdrange(self, images):

        r"""

        Computes a common standard intensity range over a number of images.

        Depending on the settings of the internal self.__stdrange variable,

        either (1) the already fixed values are returned, (2) a complete standard

        intensity range is computed from the supplied images, (3) an intensity range

        fixed at the lower end or (4) an intensity range fixed at the upper end is

        returned.

        Takes into account the maximum length of each percentile segment over all

        images, then adds a security margin defined by the highest variability among

        all segments over all images.

        Be

        .. math::

            L = (cop_l, lp_1, lp_2, ..., lp_n, cop_u)

        the set formed by the two cut-off percentiles :math:`cop_l` and :math:`cop_u` and the

        landmark percentiles :math:`lp_1, ..., lp_n`. The corresponding intensity values of

        an image :math:`i\in I` are then

        .. math::

            V_i = (v_{i,1}, v_{i,2}, ..., v_{i,n+2})

        The distance between each of these intensity values forms a segment along the

        images :math:`i` intensity range denoted as

        ..math ::

            S_i = (s_{i,1}, s_{i,2}, ..., s_{i, n+1})

        The common standard intensity range :math:`sir` over the set of images :math:`I` is

        then defined as

        ..math ::

            sir = \sum_{l=1}^{n+1}\max_{i=1}^I s_{i,l} * \max_{l=1}^{n+1} \left(\frac{\max_{i=1}^I s_{i,l}}{\min_{i=1}^I s_{i,l}}\right)

        Parameters

        ----------

        images : sequence of array_like

            A number of images.

        Returns

        -------

        stdrange : (float, float)

            The borders of the computed standard intensity range.

        """

        if not 'auto' in self.__stdrange:

            return self.__stdrange

        copl, copu = self.__cutoffp

        # collect cutoff + landmark percentile segments and image mean intensity values 

        s = []

        m = []

        for idx, i in enumerate(images):

            li = numpy.percentile(i, [copl] + self.__landmarkp + [copu])

            s.append(numpy.asarray(li)[1:] - numpy.asarray(li)[:-1])

            m.append(i.mean())

            # treat single intensity accumulation error

            if 0 in s[-1]:

                raise SingleIntensityAccumulationError('Image no.{} shows an unusual single-intensity accumulation that leads to a situation where two percentile values are equal. This situation is usually caused, when the background has not been removed from the image. Another possibility would be to reduce the number of landmark percentiles landmarkp or to change their distribution.'.format(idx))

        # select the maximum and minimum of each percentile segment over all images

        maxs = numpy.max(s, 0)

        mins = numpy.min(s, 0)

        # divide them pairwise

        divs = numpy.divide(numpy.asarray(maxs, dtype=numpy.float), mins) 

        # compute interval range according to generalized theorem 2 of [1]

        intv = numpy.sum(maxs) + numpy.max(divs)

        # compute mean intensity value over all images (assuming equal size)

        im = numpy.mean(m)

        # return interval with borders according to settings

        if 'auto' == self.__stdrange[0] and 'auto' == self.__stdrange[1]:

            return im - intv / 2, im + intv / 2

        elif 'auto' == self.__stdrange[0]:

            return self.__stdrange[1] - intv, self.__stdrange[1]

        else:

            return self.__stdrange[0], self.__stdrange[0] + intv

    def __check_mapping(self, landmarks):

        """

        Checks whether the image, from which the supplied landmarks were extracted, can

        be transformed to the learned standard intensity space without loss of

        information.

        """

        sc_udiff = numpy.asarray(self.__sc_umaxs)[1:] - numpy.asarray(self.__sc_umins)[:-1]

        l_diff = numpy.asarray(landmarks)[1:] - numpy.asarray(landmarks)[:-1]

        return numpy.all(sc_udiff > numpy.asarray(l_diff))

    @staticmethod

    def is_sequence(arg):

        """

        Checks via its hidden attribute whether the passed argument is a sequence (but

        excluding strings).

        Credits to Steve R. Hastings a.k.a steveha @ http://stackoverflow.com

        """

        return (not hasattr(arg, "strip") and

            hasattr(arg, "__getitem__") or

            hasattr(arg, "__iter__"))

    @staticmethod

    def is_number(arg):

        """

        Checks whether the passed argument is a valid number or not.

        """

        import numbers

        return isinstance(arg, numbers.Number)

    @staticmethod

    def are_numbers(arg):

        """

        Checks whether all elements in a sequence are valid numbers.

        """

        return numpy.all([IntensityRangeStandardization.is_number(x) for x in arg])

    @staticmethod

    def is_in_interval(n, l, r, border = 'included'):

        """

        Checks whether a number is inside the interval l, r.

        """

        if 'included' == border:

            return (n >= l) and (n <= r)

        elif 'excluded' == border:

            return (n > l) and (n < r)

        else:

            raise ValueError('borders must be either \'included\' or \'excluded\'')

    @staticmethod

    def are_in_interval(s, l, r, border = 'included'):

        """

        Checks whether all number in the sequence s lie inside the interval formed by

        l and r.

        """

        return numpy.all([IntensityRangeStandardization.is_in_interval(x, l, r, border) for x in s])

    @staticmethod

    def to_float(s):

        """

        Cast a sequences elements to float numbers.

        """

        return [float(x) for x in s]

    @staticmethod

    def linear_model((x1, x2), (y1, y2)):

        """

        Returns a linear model transformation function fitted on the two supplied points.

        y = m*x + b

        Note: Assumes that slope > 0, otherwise division through zero might occur.

        """

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

        b = y1 - (m * x1)

        return lambda x: m * x + b

class SingleIntensityAccumulationError(Exception):

    """

    Thrown when an image shows an unusual single-intensity peaks which would obstruct

    both, training and transformation.

    """

class InformationLossException(Exception):

    """

    Thrown when a transformation can not be guaranteed to be lossless.

    """

    pass

class UntrainedException(Exception):

    """

    Thrown when a transformation is attempted before training.

    """

    pass