<kaggle_start><data_title>OCTant project<data_name>octant-project
<code># # Enhancing an open-source OCT segmentation tool with manual segmentation capabilities
# ## Introduction
# Optical coherence tomography (OCT) is a way of capturing images for the diagnosis and monitoring of the retina. The retina is complex tissue composed of ten major layers of distinct cells. An important part of the routine use of OCT retinal images is the segmentation of the tissue layers. The segmentation consists of labelling the image pixels with their corresponding layer and remains a complex problem despite the existence of some approximate solutions.
# Around the central task of segmentation, additional functions need to be performed on the image, for instance to quantify the thickness of the tissue layers, or to stitch several images together to see a wider region of the eye. For ophthalmologists to be able to use the algorithms needed to process the OCT image, these functions are often encapsulated in user friendly applications. But alas, software to segment retinal OCT images is expensive and not always available in developing countries and/or economically disfavored areas.
# OCTant (https://github.com/forihuelaespina/OCTant/) is a python based free open-source software for the analysis of OCT images originally proposed by the National Institute of Astrophysics, Optics and Electronics (INAOE in Spanish) in Mexico but which is still under development. In its current version, OCTant can flatten or mosaic the images, and measure the thickness of each layer.
# Because the performance of even state-of-the-art segmentation algorithms is always suboptimal, ophthalmologists often benefit from being able to make manual adjustments. This project will improve the OCTant tool with the capacity to manually adjust a given segmentation
# The project will teach you:
# * Basic medical image processing.
# * To develop code valid within a collaborative project (i.e. several programming rules ought to be respected so that other collaborators can use and understand you code off-the-shelf. For instance, you will familiarize yourself with GitHub versioning platform.
# * To interact with software in development, where bugs are likely to exist, and errors are likely to appear. For instance, you will learn to report bugs.
# * To produce technical documentation in Sphynx.
# To start off with, let's take a look at an example image:
# import os
# import glob
# import numpy as np
from PIL import Image

# import matplotlib.pyplot as plt
img = Image.open("../input/octant-project/_1323904_.jpeg")
img
# what can we do to the image?
dir(img)
print(img.width, img.height)
spinPicture = img.rotate(30)
spinPicture
from PIL import Image
import matplotlib.pyplot as plt


def getRed(redVal):
    return "#%02x%02x%02x" % (redVal, 0, 0)


def getGreen(greenVal):
    return "#%02x%02x%02x" % (0, greenVal, 0)


def getBlue(blueVal):
    return "#%02x%02x%02x" % (0, 0, blueVal)


# Create an Image with specific RGB value
image = img
# Display the image
image.show()
# Get the color histogram of the image
histogram = image.histogram()
# Take only the Red counts
l1 = histogram[0:256]
# Take only the Blue counts
l2 = histogram[256:512]
# Take only the Green counts
l3 = histogram[512:768]
plt.figure(0)
# R histogram
for i in range(0, 256):
    plt.bar(i, l1[i], color=getRed(i), edgecolor=getRed(i), alpha=0.3)
# G histogram
plt.figure(1)
for i in range(0, 256):
    plt.bar(i, l2[i], color=getGreen(i), edgecolor=getGreen(i), alpha=0.3)
# B histogram
plt.figure(2)
for i in range(0, 256):
    plt.bar(i, l3[i], color=getBlue(i), edgecolor=getBlue(i), alpha=0.3)
plt.show()
# %matplotlib inline
# from matplotlib import pyplot as plt
# from gluoncv import model_zoo, data, utils
# net = model_zoo.get_model('mask_rcnn_fpn_resnet101_v1d_coco', pretrained=True)
# from PIL import Image
# import pandas as pd
# import os
# import numpy as np
# pth = '../input/octant-project/_1323904_.jpeg'
# x, orig_img = data.transforms.presets.rcnn.load_test(pth)
# ids, scores, bboxes, masks = [xx[0].asnumpy() for xx in net(x)]
# width, height = orig_img.shape[1], orig_img.shape[0]
# masks = utils.viz.expand_mask(masks, bboxes, (width, height), scores)
# orig_img = utils.viz.plot_mask(orig_img, masks)
# # identical to Faster RCNN object detection
# fig = plt.figure(figsize=(20, 20))
# ax = fig.add_subplot(1, 1, 1)
# ax = utils.viz.plot_bbox(orig_img, bboxes, scores, ids,
#                          class_names=net.classes, ax=ax)
# plt.show()
from skimage import io
import copy
import octant as oc

fileName = "../input/octant-project/_1323904_.jpeg"
img = io.imread(fileName)
scan = OCTscan(img)
scanSegmentation = OCTscanSegmentation(scan)
print("Creating document.")
study = OCTvolume()  # Initialize a document
study.addScans(scan)
segmentationVol = OCTvolumeSegmentation()  # Initialize a document
segmentationVol.addScanSegmentations(scanSegmentation)
doc = Document()  # Initialize a document
doc.name = fileName
tmp, _ = os.path.split(fileName)
doc.folderName = tmp
doc.fileName = fileName
doc.study = study
doc.segmentation = segmentationVol
# Keep reference image.
print("Replicating image.")
doc2 = copy.deepcopy(doc)
# Flattening
print("-- Flattening.")
flt = OpScanFlatten()
flt.addOperand(doc2.getCurrentScan())
imFlattened = flt.execute()
doc2.setCurrentScan(imFlattened)
# Segmentation
print("-- Segmenting.")
doc3 = copy.deepcopy(doc2)
seg = OpScanSegment()
seg.addOperand(doc3.getCurrentScan())
imSegmented = seg.execute()
doc3.setCurrentScanSegmentation(imSegmented)
# Load colormap
print("-- Plotting.")
appsettings = Settings()
appsettingsfile = "..\\resources\\OCTantApp.config"
appsettings.read(appsettingsfile)
cmap = appsettings.retinallayerscolormap
# hFig = myInitFigure()
# myPaint(hFig,doc) #Plot raw with dummy segmentation
# hFig = myPaint(myInitFigure(),doc) #Plot raw
hFig = myPaint(myInitFigure(), doc2)  # Plot flattened
hFig = myPaint(myInitFigure(), doc3)  # Plot segmented


# imgin = OCTscan(img)
class OCTscan(object):
    # Sphinx documentation
    """A single OCT scan.
    A single OCT scan. A scan is a grayscale image.
    Check .scantype property for the scan type (A, B or C)
    .. seealso:: :class:`data.OCTvolume`
    .. note:: None
    .. todo:: Upgrade to color scans.
    """

    # Class constructor
    def __init__(self, *args):
        """The class constructor.
        The class constructor.
        tmp = OCTscan() - Creates an black scan sized 480x640.
        tmp = OCTScan(img) - Creates scan from the (grayscale) image. Assumed to be an A scan.
        tmp = OCTScan(img,type) - Creates scan from the (grayscale) image.
        :param img: The scan image
        :type img: numpy.ndarray
        :param type: The scan type ('A' -default-, 'B' or 'C')
        :type type: char
        """
        if len(args) > 1:
            warnMsg = (
                self.getClassName() + ":__init__: Unexpected number of input arguments."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        # Initialize attributes (without decorator @property)
        # Initialize properties (with decorator @property)
        self.data = np.zeros(shape=(480, 640), dtype=np.uint8)  # The default scan image
        self.scantype = "A"
        if len(args) > 0:
            self.data = args[0]
        if len(args) > 1:
            self.scantype = args[1]
        return

    # Properties getters/setters
    #
    # Remember: Sphinx ignores docstrings on property setters so all
    # documentation for a property must be on the @property method
    @property
    def data(self):  # data getter
        """
        The (grayscale) scan image. The image is expected to be
        a grayscale image. Colour images will be converted to grayscale.
        :getter: Gets the OCT scan image
        :setter: Sets the OCT scan image.
        :type: numpy.ndarray shaped [width,height]
        """
        return self.__data

    @data.setter
    def data(self, img):  # data setter
        if img is not None:
            # Check whether the image is in RGB (ndim=3) or in grayscale (ndim=2)
            # and convert to grayscale if necessary
            if img.ndim == 2:
                # Dimensions are only width and height. The image is already in grayscale.
                pass
            elif img.ndim == 3:
                # Image is in RGB. Convert.
                img = color.rgb2gray(img)
            else:  # Unexpected case. Return warning
                warnMsg = self.getClassName() + ":data: Unexpected image shape."
                warnings.warn(warnMsg, SyntaxWarning)
            self.__data = img
        return None

    @property
    def shape(self):  # shape getter
        """
        The scan shape [width,height].
        :getter: Gets the scan shape
        :setter: None. This is a read-only property.
        :type: Tuple [width,height]
        """
        return self.__data.shape

    @shape.setter
    def shape(self, *args):  # shape setter
        # Catching attempts to set the shape of the scan
        warnMsg = self.getClassName() + ":shape: shape is a read-only property."
        warnings.warn(warnMsg, UserWarning)
        return

    @property
    def scantype(self):  # scantype getter
        """
        The scan type; 'A', 'B' or 'C'.
        :getter: Gets the scan type.
        :setter: Sets the scan type.
        :type: char 'A', 'B' or 'C'.
        """
        return self.__scantype

    @scantype.setter
    def scantype(self, *args):  # shape setter
        stype = args[0].upper()  # Uppercase
        if not stype in ("A", "B", "C"):
            warnMsg = (
                self.getClassName() + ":scantype: Scan type can only be "
                "A"
                ", "
                "B"
                " or "
                "C"
                "."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        self.__scantype = stype
        return None

    # Private methods
    def __str__(self):
        s = (
            "<"
            + self.getClassName()
            + "(["
            + "  data: "
            + format(self.data)
            + ","
            + "  shape: "
            + format(self.shape)
            + ","
            + "  scantype: "
            + self.scantype
            + "])>"
        )
        return s

    # Public methods
    def getClassName(self):
        """Get the class name as a string.
        Get the class name as a string.
        :returns: The class name.
        :rtype: string
        """
        return type(self).__name__


class OCTscanSegmentation(object):
    # Sphinx documentation
    """A retinal layer segmentation over a :class:`data.OCTscan`
    A retinal layer segmentation over a :class:`data.OCTscan`. A segmentation
    assigns every pixel of the scan a class label.
    Please note that this is a data model class; it keeps the segmentation
    but it is NOT capable of "computing" such segmentation. To compute a
    segmentation please refer to :class:`op.OpScanSegment`.
    The segmentation is sized and shaped equal to its base
    :class:`data.OCTscan`.
    A default segmentation sets the whole segmentation to BACKGROUND.
    .. seealso:: :class:`data.OCTscan`, :class:`op.OpScanSegment`
    .. note:: None
    .. todo:: None
    """
    _BACKGROUND = 0  # The background label identifier

    # Class constructor
    def __init__(self, *args):
        """The class constructor.
        The class constructor.
        tmp = OCTscanSegmentation(theOCTScan) - Creates a default
            segmentation for the given :class:`data.OCTscan`
        :param theOCTScan: The OCT scan to be segmented
        :type img: :class:`data.OCTscan`
        """
        refImage = OCTscan()
        # Dummy reference
        if len(args) == 0:
            warnMsg = (
                self.getClassName()
                + ":__init__: Unexpected number of input arguments. Generating a dummy reference scan."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        else:
            refImage = args[0]
            # if type(refImage) is not OCTscan:
            #     raise ErrorValue #Throw error
        # Initialize attributes (without decorator @property)
        # Initialize properties (with decorator @property)
        self.scan = refImage  # The OCT scan over which the segmentation is made
        self.data = self._BACKGROUND * np.ones(
            refImage.shape
        )  # The segmentation itself
        self.classMap = RetinalLayers().layers  # The map of class labels
        return

    # Properties getters/setters
    #
    # Remember: Sphinx ignores docstrings on property setters so all
    # documentation for a property must be on the @property method
    @property
    def data(self):  # data getter
        """
        The segmentation labels map. Please refer to :py:attr:`classMap` for
        classes.
        ..note: WARNING! This method is not currently checking whether the
            data is sized equal to the scan. This may become a problem later.
            The problem is that trying to check scan.shape will raise an
            error during object creation, when attemting to set the data
            but because the object has not been created yet, it still lacks
            the scan property even if declared in advance.
        :getter: Gets the segmentation map
        :setter: Sets the segmentation map
        :type: numpy.ndarray shaped [width,height]
        """
        return self.__data

    @data.setter
    def data(self, segmentedImg):  # data setter
        self.__data = segmentedImg
        # if segmentedImg is not None:
        #     #Check whether the image is in RGB (ndim=3) or in grayscale (ndim=2)
        #     #and convert to grayscale if necessary
        #     if ((segmentedImg.ndim == 2) & (segmentedImg.shape == self.scan.shape)):
        #         #Dimensions are only width and height, and matches that of
        #         #the scan.
        #         self.__data = segmentedImg;
        #     else: #Unexpected case. Return warning
        #         warnMsg = self.getClassName() + ':data: Unexpected segmentation shape.'
        #         warnings.warn(warnMsg,SyntaxWarning)
        return None

    @property
    def scan(self):  # scan getter
        """
        The base OCT scan. Please refer to :py:attr:`data` for
        the segmentation map.
        :getter: Gets the base OCT scan
        :setter: Sets the base OCT scan
        :type: :class:`data.OCTscan`
        """
        return self.__scan

    @scan.setter
    def scan(self, octScan):  # scan setter
        if octScan is not None:
            # Check whether the image is in RGB (ndim=3) or in grayscale (ndim=2)
            # and convert to grayscale if necessary
            if type(octScan) is OCTscan:
                # Dimensions are only width and height, and matches that of
                # the scan.
                self.__scan = octScan
                self.clear()
            else:  # Unexpected case. Return warning
                warnMsg = self.getClassName() + ":data: Unexpected type for OCT scan."
                warnings.warn(warnMsg, SyntaxWarning)
        return None

    @property
    def shape(self):  # shape getter
        """
        The scan segmentation shape [width,height].
        :getter: Gets the scan segmentation shape
        :setter: None. This is a read-only property.
        :type: Tuple [width,height]
        """
        return self.__data.shape

    @shape.setter
    def shape(self, *args):  # shape setter
        # Catching attempts to set the shape of the scan
        warnMsg = self.getClassName() + ":shape: shape is a read-only property."
        warnings.warn(warnMsg, UserWarning)
        return

    @property
    def classMap(self):  # classMap getter
        """
        The map of classes.
        The map of classes; the list of class names associated to each
        value in the segmentation map.
        ..note: This list does NOT include the BACKGROUND class.
        :getter: Gets the base OCT scan
        :setter: Sets the base OCT scan
        :type: :class:`data.OCTscan`
        """
        return self.__classMap

    @classMap.setter
    def classMap(self, cm):  # classMap setter
        if cm is not None:
            # Check that we are receiving the correct type
            if type(cm) is dict:
                self.__classMap = cm
            else:  # Unexpected case. Return warning
                warnMsg = (
                    self.getClassName() + ":classMap: Unexpected type for classMap."
                )
                warnings.warn(warnMsg, SyntaxWarning)
        return None

    # Private methods
    def __str__(self):
        s = (
            "<"
            + self.getClassName()
            + "(["
            + "  scan: "
            + format(self.scan)
            + ","
            + "  data: "
            + format(self.data)
            + ","
            + "  classMap: "
            + format(self.classMap)
            + "])>"
        )
        return s

    # Public methods
    def getClassName(self):
        """Get the class name as a string.
        Get the class name as a string.
        :returns: The class name.
        :rtype: string
        """
        return type(self).__name__

    def clear(self):
        """Clears/Resets the segmentation map to _BACKGROUND.
        Clears/Resets the segmentation map to _BACKGROUND. All pixels are
        assigned the background label.
        """
        self.data = self._BACKGROUND * np.ones(self.scan.shape)
        return None


class OCTscan(object):
    # Sphinx documentation
    """A single OCT scan.
    A single OCT scan. A scan is a grayscale image.
    Check .scantype property for the scan type (A, B or C)
    .. seealso:: :class:`data.OCTvolume`
    .. note:: None
    .. todo:: Upgrade to color scans.
    """

    # Class constructor
    def __init__(self, *args):
        """The class constructor.
        The class constructor.
        tmp = OCTscan() - Creates an black scan sized 480x640.
        tmp = OCTScan(img) - Creates scan from the (grayscale) image. Assumed to be an A scan.
        tmp = OCTScan(img,type) - Creates scan from the (grayscale) image.
        :param img: The scan image
        :type img: numpy.ndarray
        :param type: The scan type ('A' -default-, 'B' or 'C')
        :type type: char
        """
        if len(args) > 1:
            warnMsg = (
                self.getClassName() + ":__init__: Unexpected number of input arguments."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        # Initialize attributes (without decorator @property)
        # Initialize properties (with decorator @property)
        self.data = np.zeros(shape=(480, 640), dtype=np.uint8)  # The default scan image
        self.scantype = "A"
        if len(args) > 0:
            self.data = args[0]
        if len(args) > 1:
            self.scantype = args[1]
        return

    # Properties getters/setters
    #
    # Remember: Sphinx ignores docstrings on property setters so all
    # documentation for a property must be on the @property method
    @property
    def data(self):  # data getter
        """
        The (grayscale) scan image. The image is expected to be
        a grayscale image. Colour images will be converted to grayscale.
        :getter: Gets the OCT scan image
        :setter: Sets the OCT scan image.
        :type: numpy.ndarray shaped [width,height]
        """
        return self.__data

    @data.setter
    def data(self, img):  # data setter
        if img is not None:
            # Check whether the image is in RGB (ndim=3) or in grayscale (ndim=2)
            # and convert to grayscale if necessary
            if img.ndim == 2:
                # Dimensions are only width and height. The image is already in grayscale.
                pass
            elif img.ndim == 3:
                # Image is in RGB. Convert.
                img = color.rgb2gray(img)
            else:  # Unexpected case. Return warning
                warnMsg = self.getClassName() + ":data: Unexpected image shape."
                warnings.warn(warnMsg, SyntaxWarning)
            self.__data = img
        return None

    @property
    def shape(self):  # shape getter
        """
        The scan shape [width,height].
        :getter: Gets the scan shape
        :setter: None. This is a read-only property.
        :type: Tuple [width,height]
        """
        return self.__data.shape

    @shape.setter
    def shape(self, *args):  # shape setter
        # Catching attempts to set the shape of the scan
        warnMsg = self.getClassName() + ":shape: shape is a read-only property."
        warnings.warn(warnMsg, UserWarning)
        return

    @property
    def scantype(self):  # scantype getter
        """
        The scan type; 'A', 'B' or 'C'.
        :getter: Gets the scan type.
        :setter: Sets the scan type.
        :type: char 'A', 'B' or 'C'.
        """
        return self.__scantype

    @scantype.setter
    def scantype(self, *args):  # shape setter
        stype = args[0].upper()  # Uppercase
        if not stype in ("A", "B", "C"):
            warnMsg = (
                self.getClassName() + ":scantype: Scan type can only be "
                "A"
                ", "
                "B"
                " or "
                "C"
                "."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        self.__scantype = stype
        return None

    # Private methods
    def __str__(self):
        s = (
            "<"
            + self.getClassName()
            + "(["
            + "  data: "
            + format(self.data)
            + ","
            + "  shape: "
            + format(self.shape)
            + ","
            + "  scantype: "
            + self.scantype
            + "])>"
        )
        return s

    # Public methods
    def getClassName(self):
        """Get the class name as a string.
        Get the class name as a string.
        :returns: The class name.
        :rtype: string
        """
        return type(self).__name__


class RetinalLayers(object):
    # Sphinx documentation
    """A collection of constants for identifying retinal layers.
    A collection of constants for identifying retinal layers.
    It is basically a dictionary of pairs key:value
    The following retinal layers are considered:
    * Inner Limiting Membrane (ILM)
    * Nerve Fiber Layer (NFL)
    * Ganglion Cell Layer (GCL)
    * Inner Plexiform Layer (IPL)
    * Inner Nuclear Layer (INL)
    * Outer Plexiform Layer (OPL)
    * Outner Nuclear Layer (ONL)
    * External Limiting Membrane (ELM)
    * Rods and Cones layers (RCL)
    * RetinalPigmentEpithelium (RPE)
    Also the non-retinal layer is indexed:
    * Choroid (CHR)
    For pathological conditions, the following is also encoded:
    * Fluid (FLD)
    .. seealso:: None
    .. note:: None
    .. todo:: None
    """

    # Private class attributes shared by all instances
    # Class constructor
    def __init__(self):
        """The class constructor.
        The class constructor.
        tmp = RetinalLayers() - Creates an RetinalLayers object.
        """
        # Initialize attributes
        # Layer constants
        self.layers = {
            "Inner Limiting Membrane": 1,
            "Nerve Fiber Layer": 2,
            "Ganglion Cell Layer": 3,
            "Inner Plexiform Layer": 4,
            "Inner Nuclear Layer": 5,
            "Outer Plexiform Layer": 6,
            "Outer Nuclear Layer": 7,
            "External Limiting Membrane": 8,
            "Rods and Cones Layers": 9,
            "Retinal Pigment Epithelium": 10,
            "Choroid": 20,
            "Fluid": 1001,
        }
        # Currently:
        # RPE includes Bruch's membrane even though most books will
        # 	consider BM as part of the choroid already.
        # ONL includes Henle’s layer
        # ELM is also known as Inner Segment Layer
        # The RCL includes connecting cilia (CL), outer segment layer (OSL) and Verhoeff membrane (VM)
        # ...and so that one can have different names for each layer
        self.layerNames = {
            "ilm": self.layers["Inner Limiting Membrane"],
            "innerlimitingmembrane": self.layers["Inner Limiting Membrane"],
            "nfl": self.layers["Nerve Fiber Layer"],
            "nervefiberlayer": self.layers["Nerve Fiber Layer"],
            "gcl": self.layers["Ganglion Cell Layer"],
            "ganglioncelllayer": self.layers["Ganglion Cell Layer"],
            "ipl": self.layers["Inner Plexiform Layer"],
            "innerplexiformlayer": self.layers["Inner Plexiform Layer"],
            "inl": self.layers["Inner Nuclear Layer"],
            "innernuclearlayer": self.layers["Inner Nuclear Layer"],
            "opl": self.layers["Outer Plexiform Layer"],
            "outerplexiformlayer": self.layers["Outer Plexiform Layer"],
            "onl": self.layers["Outer Nuclear Layer"],
            "outernuclearlayer": self.layers["Outer Nuclear Layer"],
            "elm": self.layers["External Limiting Membrane"],
            "externallimitingmembrane": self.layers["External Limiting Membrane"],
            "rcl": self.layers["Rods and Cones Layers"],
            "rodsandconeslayers": self.layers["Rods and Cones Layers"],
            "rpe": self.layers["Retinal Pigment Epithelium"],
            "retinalpigmentepithelium": self.layers["Retinal Pigment Epithelium"],
            "chr": self.layers["Choroid"],
            "fld": self.layers["Fluid"],
        }
        # Layer acronyms
        self.layerAcronyms = {
            self.layers["Inner Limiting Membrane"]: "ILM",
            self.layers["Nerve Fiber Layer"]: "NFL",
            self.layers["Ganglion Cell Layer"]: "GCL",
            self.layers["Inner Plexiform Layer"]: "IPL",
            self.layers["Inner Nuclear Layer"]: "INL",
            self.layers["Outer Plexiform Layer"]: "OPL",
            self.layers["Outer Nuclear Layer"]: "ONL",
            self.layers["External Limiting Membrane"]: "ELM",
            self.layers["Rods and Cones Layers"]: "RCL",
            self.layers["Retinal Pigment Epithelium"]: "RPE",
            self.layers["Choroid"]: "Choroid",
            self.layers["Fluid"]: "Fluid",
        }
        return

    # Private methods
    # Public methods
    def getClassName(self):
        """Gets the class name
        return: The class name
        rtype: string
        """
        return type(self).__name__

    def getAllLayersIndexes(self):
        """Gets the list of layer values
        return: The list of layers values
        rtype: list
        """
        return list(self.layers.values())

    def getAllLayersNames(self):
        """Gets the list of layer keys
        return: The list of layers keys
        rtype: list
        """
        # Retrieves a list of layer keys
        return list(self.layers.keys())

    def getLayerAcronym(self, idx):
        """Gets the acronym of the i-th layer
        return: The layer acronym e.g. NFL
        rtype: string
        """
        lacronym = "NaN"
        try:
            lacronym = self.layerAcronyms[idx]
        except:
            lacronym = "Unknown"
            print(
                self.getClassName(),
                ':getLayerAcronym: Unexpected layer index. Returning name "',
                lacronym,
                '"',
            )
        return lacronym

    def getLayerIndex(self, layerName):
        """Retrieve the index of a given layer
        return: The index of the layer
        rtype: int
        """
        r = -1
        try:
            layerName = layerName.replace(" ", "")  # Remove whitespaces
            r = self.layerNames[layerName.lower()]  # Ignore case
        except:
            print(
                self.getClassName(),
                ":getLayerIndex: Unknown layer name. Returning index ",
                r,
            )
        return r

    def getLayerName(self, idx):
        """Retrieve the i-th layer name
        return: The name of the i-th layer
        rtype: string
        """
        lname = "Default"
        try:
            # There is no 'direct' method to access the keys given the value.
            lname = list(self.layers.keys())[list(self.layers.values()).index(idx)]
        except:
            lname = "Unknown"
            print(
                self.getClassName(),
                ':getLayerName: Unexpected layer index. Returning name "',
                lname,
                '"',
            )
        return lname

    def getNumLayers(self):
        """Return the number of known layers.
        Return the number of known layers. Please note that this also
        include known non-retinal layers like the choroid.
        return: Length of property map :func:`data.OCTvolume.layers`
        rtype: int
        """
        return len(self.layers)


## Import
import os
import warnings

# from deprecated import deprecated
import deprecation
import math
import numpy as np
from scipy import signal, ndimage  # Used for 2D convolution
from skimage import feature, color
from functools import reduce
import cv2  # That's OpenCV

# import matlab.engine
import matplotlib.pyplot as plt
import _tkinter
from PyQt5.QtCore import Qt  # Imports constants
from PyQt5.QtWidgets import QProgressDialog

# from matplotlib.backends.qt_compat import QtCore, QtWidgets, is_pyqt5
# from version import __version__
import octant

# from data import OCTscan, OCTscanSegmentation, RetinalLayers
# from util import segmentationUtils
# from .Operation import Operation
from scipy.optimize import curve_fit


class OpScanSegment(Operation):
    # Private class attributes shared by all instances
    # Class constructor
    def __init__(self):
        # Call superclass constructor
        super().__init__()
        # Set the operation name
        self.name = "Segmentation"
        # Initialize private attributes unique to this instance
        # self._imgin = np.zeros(shape = (0,0,0), dtype = np.uint8 ); #Input image
        # self._imgout = np.zeros(shape = (0,0,0), dtype = np.uint8 ); #The segmented image

    # Private methods
    @staticmethod
    def anisotropicDiffusionFilter(img):
        """
        Non-linear anisotropic diffusion filtering
        #Original algorithm in [Perona and Malik, (1990) TPAMI 12(7):629-639]
        #Parameters according to [WangRK2005, Proc. of SPIE 5690:380-385]
        :param img: ndarray. The image to be filtered
        :return: the image filtered
        :rtype: ndarray
        """
        sigma = 1
        GaussMask = generateGaussianMask2D(shape=(3, 3), sigma=sigma)
        # print(GaussMask)
        m = 8
        l = 10  # contrast parameter; structures with s > λ are regarded as edges,
        # while with s < λ are assumed to belong to the interior of a region
        # [Wang2005] uses l=5 for the porcine trachea, but for us that value
        # does not produce as good results as l=10.
        timestep = 0.24
        niterations = 10  # diffusion time (iteration)
        Cm = 3.31488
        img2 = img  # Initialize image
        # Get rid of the "large" noise with morphological closing after opening
        morphkernel = np.ones((5, 5), np.uint8)
        # img2 = cv2.morphologyEx(img2, cv2.MORPH_OPEN, morphkernel)
        # img2 = cv2.morphologyEx(img2, cv2.MORPH_CLOSE, morphkernel)
        if np.max(img2) <= 1:
            # Scale parameters
            # m=m/255
            l = l / 255
            # Cm = Cm/255
        for tau in range(1, niterations):
            # 			#Progress bar
            # 			if progress.wasCanceled():
            # 				break
            # 			progress.setValue(round(100*tau/niterations))
            # Morphological removal of noise
            # Note that this is NOT part of the original diffusion filter
            # but the results are clearly enhanced!.
            img2 = cv2.morphologyEx(img2, cv2.MORPH_OPEN, morphkernel)
            img2 = cv2.morphologyEx(img2, cv2.MORPH_CLOSE, morphkernel)
            # Estimate gradient vector at sigma scale
            gradI = np.gradient(img2)
            tmpGradI = np.add.reduce(gradI)  # For illustrations purposes only
            # Regularize gradient
            # gradientScaleSigma=math.sqrt(2*(tau*timestep))
            gradientScaleSigma = math.sqrt(2 * (tau))
            GaussMask = generateGaussianMask2D(shape=(7, 7), sigma=gradientScaleSigma)
            # Individually convolve the Gaussian filter with each gradient component
            # s = signal.convolve2d(gradI, GaussMask, boundary='symm', mode='same')
            s = [None] * len(gradI)  # Preallocate list
            for dim in range(0, len(gradI)):
                s[dim] = signal.convolve2d(
                    gradI[dim], GaussMask, boundary="symm", mode="same"
                )
            s = np.linalg.norm(s, ord=2, axis=0)  # Norm of the gradient.
            # Calculate diffusivity
            tmp = (s / l) ** m
            tmp = np.divide(
                -Cm, tmp, out=np.zeros(tmp.shape), where=tmp != 0
            )  # Avoid division by 0 when estimating diffusivity
            # D = 1-np.exp(-Cm/((s/l)**m)) #diffusivity or conduction coefficient
            D = 1 - np.exp(tmp)  # diffusivity or conduction coefficient
            # Update image
            img2 = img2 + divergence(np.multiply(D, gradI))  # Update the image
            # Reminder> The divergence of gradient is the Laplacian operator
            # See: https://math.stackexchange.com/questions/690493/what-is-divergence-in-image-processing
        return img2

    @staticmethod
    def divergence(f):
        # """Compute the divergence of n-D SCALAR field `f`."""
        # See: https://stackoverflow.com/questions/11435809/compute-divergence-of-vector-field-using-python
        # return reduce(np.add,np.gradient(f))
        """
        Computes the divergence of the VECTOR field f, corresponding to dFx/dx + dFy/dy + ...
        :param f: List of ndarrays, where every item of the list is one dimension of the vector field
        :return: Single ndarray of the same shape as each of the items in f, which corresponds to a scalar field
        """
        num_dims = len(f)
        return np.ufunc.reduce(
            np.add, [np.gradient(f[i], axis=i) for i in range(num_dims)]
        )

    @staticmethod
    def generateGaussianMask2D(shape=(3, 3), sigma=0.5):
        """
        Generates a 2D gaussian mask
        It should give the same result as MATLAB's
        fspecial('gaussian',[shape],[sigma])
        See: https://stackoverflow.com/questions/17190649/how-to-obtain-a-gaussian-filter-in-python
        """
        m, n = [(ss - 1.0) / 2.0 for ss in shape]
        y, x = np.ogrid[-m : m + 1, -n : n + 1]
        h = np.exp(-(x * x + y * y) / (2.0 * sigma * sigma))
        h[h < np.finfo(h.dtype).eps * h.max()] = 0
        sumh = h.sum()
        if sumh != 0:
            h /= sumh
        return h

    # Public methods
    def execute(self, *args, **kwargs):
        """Executes the operation on the :py:attr:`operands`.
        Executes the operation on the :py:attr:`operands` and stores the outcome
        in :py:attr:`result`. Preload operands using method
        :func:`addOperand()`.
        :returns: Result of executing the operation.
        :rtype: :class:`data.OCTscanSegmentation`
        """
        # Ensure the operand has been set.
        if len(self.operands) < 1:
            warnMsg = self.getClassName() + ":execute: Operand not set."
            warnings.warn(warnMsg, SyntaxWarning)
            return None
        # Establish mode of operation
        # NOT WORKING YET... :(
        # Now, it always detect:
        # 	Terminal for the stdin
        # 	GUI for the stdout
        mode = "terminal"
        if os.isatty(0):
            mode = "gui"
        print("Executing segmentation in mode " + mode + ".")
        if mode == "gui":
            MainWindow = QtWidgets.QWidget()
            progress = QProgressDialog("Segmentation...", "Cancel", 0, 100, MainWindow)
            progress.setWindowModality(Qt.WindowModal)
            progress.setAutoReset(True)
            progress.setAutoClose(True)
            progress.setMinimum(0)
            progress.setMaximum(100)
            progress.setWindowTitle("Automatic segmentation")
            progress.setLabelText("Progress:")
            progress.setMinimumDuration(0)
            progress.resize(500, 100)
            progress.forceShow()
            progress.setValue(0)
        imgin = self.operands[0]
        if type(imgin) is OCTscan:
            imgin = imgin.data
        # Define a default output
        segmentedImage = np.zeros(shape=(0, 0, 0), dtype=np.uint8)
        # Check whether the image is in RGB (ndim=3) or in grayscale (ndim=2)
        # and convert to grayscale if necessary
        # img = cv2.cvtColor(self._imgin, cv2.COLOR_BGR2GRAY)
        if imgin.ndim == 2:
            # Dimensions are only width and height. The image is already in grayscale.
            img = imgin
        elif imgin.ndim == 3:
            # Image is in RGB. Convert.
            img = color.rgb2gray(imgin)
        else:  # Unexpected case. Return warning
            print(self.getClassName(), ": Unexpected image shape.")
            return None
        ## Algoritmo 1: Basado en el código de Matlab CASeReL
        # https://pangyuteng.github.io/caserel/
        #
        # Llama externamente a retSegment.exe
        #
        # Requiere de Matlab Runtime environment
        #  #Remove external noise from image
        #  img = segmentationUtils.ejecuta_close(img,4,4) #Clausura
        #  #segmentationUtils.mostrar_imagen(img)
        #
        #  #Generate temporal intermediate image file to be processed externally
        #  tmpFilename = "paso2.tiff"
        #  cv2.imwrite(tmpFilename, img)
        #
        #  #Rely externally on the matlab algorithm for segmentation
        # # eng = matlab.engine.start_matlab()
        # # eng.retSegment(img)
        #  #segmentationUtils.mostrar_imagen(img)
        #
        #  #Delete the intermediate image file
        #  segmentationUtils.elimina_imagen(tmpFilename)
        ## Algoritmo 3 Felipe
        # Step 1) Non-linear anisotropic diffusion filtering
        # Original algorithm in [Perona and Malik, (1990) TPAMI 12(7):629-639]
        # Parameters according to [WangRK2005, Proc. of SPIE 5690:380-385]
        if mode == "gui":
            progress.setLabelText("Applying diffusion filter:")
        img2 = anisotropicDiffusionFilter(img)
        # The output of the anisotropic filter is a float64 2D array
        # with values between [0 1] (although note that IT MAY not
        # include the segment boundaries 0 and/or 1, so for instance
        # while testing I was getting [0,0.85]). Anyway, in this
        # range, anything below 1% (~0.099) is background, and the
        # rest is tissue.
        if mode == "gui":
            progress.setValue(100)
            if progress.wasCanceled():
                return img  # Return the original image
        # print('Diffusion Filter output range: ['+str(np.amin(img2))+','+str(np.amax(img2))+']')
        # Detect background mask
        BGmask = np.where(img2 <= 0.099, True, False)
        # See above for threshold 0.099 on anisotropic diffusion filter
        # output.
        # Note that there might be "background" pixels
        # found somewhere within the tissue. This is VERY LIKELY fluid!
        # Note as well that the lower part also include a bit of the
        # choroid with this threshold.
        # Finally, the top part should be the vitreous.
        # print(BGmask)
        # Just to see the mask
        # img2 = np.where(BGmask,10,OCTscanSegmentation._BACKGROUND)
        # 		#Plot histogram of image
        # 		binSize = 0.01
        # 		bins=np.arange(0,np.amax(img2)+binSize,binSize)
        # 		hist = np.histogram(img2, bins=bins)
        # 		hfig, ax = plt.subplots(1, 1)
        # 		ax.hist(img2.flatten(),bins)
        # 		ax.set_yscale('log')
        # 		hfig.show()
        # PENDING
        # Assign intratissue BG pixels to fluid.
        # Normalize
        img2 = np.floor(255 * (img2 / np.amax(img2))).astype(int)
        img2 = np.where(BGmask == True, OCTscanSegmentation._BACKGROUND, img2)
        # Assign integers to nRetinalLayers
        r = RetinalLayers()
        nRetinalLayers = r.getLayerIndex("RPE")
        # Note that I'm not counting for choroid or fluid.
        # I can't use pixel intensity alone here for segmenting as several
        # (separated) layers exhibit similar reflectivity, so segmentation
        # criteria has to be a mix of location and intensity.
        # Estimate upper and lower boundary by checking first and last
        # non-BG pixels per column
        # UNFINISHED.
        # upperBoundary =
        # [uVals,uIdx,uInv,uCounts] = np.unique(img2, \
        # 						return_index=True, \
        # 						return_inverse=True, \
        # 						return_counts=True)
        # print(len(uVals))
        # 		## Algoritmo 2: Arlem
        # 		img2 = img
        # 		print(np.amax(img2))
        # 		tmp=np.reshape(img2,(np.size(img2),1))
        # 		print(type(tmp))
        # 		print(np.size(tmp))
        # 		print(np.amax(tmp))
        #
        #
        # 		hfig, ax = plt.subplots(1, 1)
        # 		ax.hist(tmp,255)
        # 		hfig.show()
        #
        # 		#Elimina ruido
        # 		img = segmentationUtils.ejecuta_close(img,4,4)
        #
        # 		#Amplifica capas
        # 		img = segmentationUtils.ejecuta_dilate(img,5,20,1)
        #
        # 		#Tensor
        # 		Axx, Axy, Ayy = feature.structure_tensor(img)
        #
        # 		#Elimina mas ruido
        # 		Ayy = segmentationUtils.ejecuta_close(Ayy,6,1)
        #
        # 		#Resalta las capas que sean mayores a la media
        # 		Ayy = segmentationUtils.resalta_bordes(Ayy,True,0)
        #
        # 		#Elimina aun mas ruido
        # 		Ayy = segmentationUtils.ejecuta_open(Ayy,1,1)
        #
        # 		#Binarizacion
        # 		binary = segmentationUtils.ejecuta_OTSU(Ayy)
        #
        # 		#elimina ruido del posible borde superior
        # 		binary = segmentationUtils.ejecuta_elimina_ruido_extremos(True,0,0,binary)
        #
        # 		#elimina ruido del posible borde inferior
        # 		binary = segmentationUtils.ejecuta_elimina_ruido_extremos(False,0,0,binary)
        #
        # 		#obtiene bordes exteriores
        # 		arraySuperior, arrayInferior = segmentationUtils.obten_bordes_externos(binary)
        #
        # 		#elimina ruido a la imagen original
        # 		img2 = segmentationUtils.elimina_desde_arreglos(img2, arraySuperior, arrayInferior)
        # 		img2 = segmentationUtils.ejecuta_close(img2,2,1)
        # 		img2 = feature.canny(img2,sigma = 2.5)
        # 		img2 = segmentationUtils.elimina_ruido_canny(img2,1)
        # Hide and close progress bar.
        if mode == "gui":
            progress.hide()
        # if isinstance(imgin,(OCTscan,)):
        if type(imgin) is OCTscan:
            self.result = OCTscanSegmentation(imgin)
        else:
            self.result = OCTscanSegmentation(OCTscan(imgin))
        self.result.data = img2
        return self.result


# 	#@deprecated(version='0.2', reason="Deprecated. Use method execute() instead.")
# 	@deprecation.deprecated(deprecated_in="0.2", removed_in="0.3",
# 						current_version=__version__,
# 						details="Use method execute() instead.")
# 	def segmentar(self,image):
# 		#Encapsulate the image as an OCTscan
# 		tmp=OCTscan(image)
# 		self.clear()
# 		self.addOperand(tmp)
# 		#Execute
# 		self.execute()
# 		return None
def anisotropicDiffusionFilter(img):
    """
    Non-linear anisotropic diffusion filtering
    #Original algorithm in [Perona and Malik, (1990) TPAMI 12(7):629-639]
    #Parameters according to [WangRK2005, Proc. of SPIE 5690:380-385]
    :param img: ndarray. The image to be filtered
    :return: the image filtered
    :rtype: ndarray
    """
    sigma = 1
    GaussMask = generateGaussianMask2D(shape=(3, 3), sigma=sigma)
    # print(GaussMask)
    m = 8
    l = 10  # contrast parameter; structures with s > λ are regarded as edges,
    # while with s < λ are assumed to belong to the interior of a region
    # [Wang2005] uses l=5 for the porcine trachea, but for us that value
    # does not produce as good results as l=10.
    timestep = 0.24
    niterations = 10  # diffusion time (iteration)
    Cm = 3.31488
    img2 = img  # Initialize image
    # Get rid of the "large" noise with morphological closing after opening
    morphkernel = np.ones((5, 5), np.uint8)
    # img2 = cv2.morphologyEx(img2, cv2.MORPH_OPEN, morphkernel)
    # img2 = cv2.morphologyEx(img2, cv2.MORPH_CLOSE, morphkernel)
    if np.max(img2) <= 1:
        # Scale parameters
        # m=m/255
        l = l / 255
        # Cm = Cm/255
    for tau in range(1, niterations):
        # 			#Progress bar
        # 			if progress.wasCanceled():
        # 				break
        # 			progress.setValue(round(100*tau/niterations))
        # Morphological removal of noise
        # Note that this is NOT part of the original diffusion filter
        # but the results are clearly enhanced!.
        img2 = cv2.morphologyEx(img2, cv2.MORPH_OPEN, morphkernel)
        img2 = cv2.morphologyEx(img2, cv2.MORPH_CLOSE, morphkernel)
        # Estimate gradient vector at sigma scale
        gradI = np.gradient(img2)
        tmpGradI = np.add.reduce(gradI)  # For illustrations purposes only
        # Regularize gradient
        # gradientScaleSigma=math.sqrt(2*(tau*timestep))
        gradientScaleSigma = math.sqrt(2 * (tau))
        GaussMask = generateGaussianMask2D(shape=(7, 7), sigma=gradientScaleSigma)
        # Individually convolve the Gaussian filter with each gradient component
        # s = signal.convolve2d(gradI, GaussMask, boundary='symm', mode='same')
        s = [None] * len(gradI)  # Preallocate list
        for dim in range(0, len(gradI)):
            s[dim] = signal.convolve2d(
                gradI[dim], GaussMask, boundary="symm", mode="same"
            )
        s = np.linalg.norm(s, ord=2, axis=0)  # Norm of the gradient.
        # Calculate diffusivity
        tmp = (s / l) ** m
        tmp = np.divide(
            -Cm, tmp, out=np.zeros(tmp.shape), where=tmp != 0
        )  # Avoid division by 0 when estimating diffusivity
        # D = 1-np.exp(-Cm/((s/l)**m)) #diffusivity or conduction coefficient
        D = 1 - np.exp(tmp)  # diffusivity or conduction coefficient
        # Update image
        img2 = img2 + divergence(np.multiply(D, gradI))  # Update the image
        # Reminder> The divergence of gradient is the Laplacian operator
        # See: https://math.stackexchange.com/questions/690493/what-is-divergence-in-image-processing
    return img2


segmentedImage = np.zeros(shape=(0, 0, 0), dtype=np.uint8)
img2 = anisotropicDiffusionFilter(img)
BGmask = np.where(img2 <= 0.099, True, False)
img2 = np.floor(255 * (img2 / np.amax(img2))).astype(int)
img2 = np.where(BGmask == True, OCTscanSegmentation._BACKGROUND, img2)
# Assign integers to nRetinalLayers
r = RetinalLayers()
nRetinalLayers = r.getLayerIndex("RPE")
OCTscanSegmentation(imgin)


def divergence(f):
    # """Compute the divergence of n-D SCALAR field `f`."""
    # See: https://stackoverflow.com/questions/11435809/compute-divergence-of-vector-field-using-python
    # return reduce(np.add,np.gradient(f))
    """
    Computes the divergence of the VECTOR field f, corresponding to dFx/dx + dFy/dy + ...
    :param f: List of ndarrays, where every item of the list is one dimension of the vector field
    :return: Single ndarray of the same shape as each of the items in f, which corresponds to a scalar field
    """
    num_dims = len(f)
    return np.ufunc.reduce(np.add, [np.gradient(f[i], axis=i) for i in range(num_dims)])


def generateGaussianMask2D(shape=(3, 3), sigma=0.5):
    """
    Generates a 2D gaussian mask
    It should give the same result as MATLAB's
    fspecial('gaussian',[shape],[sigma])
    See: https://stackoverflow.com/questions/17190649/how-to-obtain-a-gaussian-filter-in-python
    """
    m, n = [(ss - 1.0) / 2.0 for ss in shape]
    y, x = np.ogrid[-m : m + 1, -n : n + 1]
    h = np.exp(-(x * x + y * y) / (2.0 * sigma * sigma))
    h[h < np.finfo(h.dtype).eps * h.max()] = 0
    sumh = h.sum()
    if sumh != 0:
        h /= sumh
    return h


generateGaussianMask2D()


def execute(self, *args, **kwargs):
    """Executes the operation on the :py:attr:`operands`.
    Executes the operation on the :py:attr:`operands` and stores the outcome
    in :py:attr:`result`. Preload operands using method
    :func:`addOperand()`.
    :returns: Result of executing the operation.
    :rtype: :class:`data.OCTscanSegmentation`
    """
    # Ensure the operand has been set.
    if len(self.operands) < 1:
        warnMsg = self.getClassName() + ":execute: Operand not set."
        warnings.warn(warnMsg, SyntaxWarning)
        return None
    # Establish mode of operation
    # NOT WORKING YET... :(
    # Now, it always detect:
    # 	Terminal for the stdin
    # 	GUI for the stdout
    mode = "terminal"
    if os.isatty(0):
        mode = "gui"
    print("Executing segmentation in mode " + mode + ".")
    if mode == "gui":
        MainWindow = QtWidgets.QWidget()
        progress = QProgressDialog("Segmentation...", "Cancel", 0, 100, MainWindow)
        progress.setWindowModality(Qt.WindowModal)
        progress.setAutoReset(True)
        progress.setAutoClose(True)
        progress.setMinimum(0)
        progress.setMaximum(100)
        progress.setWindowTitle("Automatic segmentation")
        progress.setLabelText("Progress:")
        progress.setMinimumDuration(0)
        progress.resize(500, 100)
        progress.forceShow()
        progress.setValue(0)
    imgin = self.operands[0]
    if type(imgin) is OCTscan:
        imgin = imgin.data
    # Define a default output
    segmentedImage = np.zeros(shape=(0, 0, 0), dtype=np.uint8)
    # Check whether the image is in RGB (ndim=3) or in grayscale (ndim=2)
    # and convert to grayscale if necessary
    # img = cv2.cvtColor(self._imgin, cv2.COLOR_BGR2GRAY)
    if imgin.ndim == 2:
        # Dimensions are only width and height. The image is already in grayscale.
        img = imgin
    elif imgin.ndim == 3:
        # Image is in RGB. Convert.
        img = color.rgb2gray(imgin)
    else:  # Unexpected case. Return warning
        print(self.getClassName(), ": Unexpected image shape.")
        return None
    ## Algoritmo 1: Basado en el código de Matlab CASeReL
    # https://pangyuteng.github.io/caserel/
    #
    # Llama externamente a retSegment.exe
    #
    # Requiere de Matlab Runtime environment
    #  #Remove external noise from image
    #  img = segmentationUtils.ejecuta_close(img,4,4) #Clausura
    #  #segmentationUtils.mostrar_imagen(img)
    #
    #  #Generate temporal intermediate image file to be processed externally
    #  tmpFilename = "paso2.tiff"
    #  cv2.imwrite(tmpFilename, img)
    #
    #  #Rely externally on the matlab algorithm for segmentation
    # # eng = matlab.engine.start_matlab()
    # # eng.retSegment(img)
    #  #segmentationUtils.mostrar_imagen(img)
    #
    #  #Delete the intermediate image file
    #  segmentationUtils.elimina_imagen(tmpFilename)
    ## Algoritmo 3 Felipe
    # Step 1) Non-linear anisotropic diffusion filtering
    # Original algorithm in [Perona and Malik, (1990) TPAMI 12(7):629-639]
    # Parameters according to [WangRK2005, Proc. of SPIE 5690:380-385]
    if mode == "gui":
        progress.setLabelText("Applying diffusion filter:")
    img2 = anisotropicDiffusionFilter(img)
    # The output of the anisotropic filter is a float64 2D array
    # with values between [0 1] (although note that IT MAY not
    # include the segment boundaries 0 and/or 1, so for instance
    # while testing I was getting [0,0.85]). Anyway, in this
    # range, anything below 1% (~0.099) is background, and the
    # rest is tissue.
    if mode == "gui":
        progress.setValue(100)
        if progress.wasCanceled():
            return img  # Return the original image
    # print('Diffusion Filter output range: ['+str(np.amin(img2))+','+str(np.amax(img2))+']')
    # Detect background mask
    BGmask = np.where(img2 <= 0.099, True, False)
    # See above for threshold 0.099 on anisotropic diffusion filter
    # output.
    # Note that there might be "background" pixels
    # found somewhere within the tissue. This is VERY LIKELY fluid!
    # Note as well that the lower part also include a bit of the
    # choroid with this threshold.
    # Finally, the top part should be the vitreous.
    # print(BGmask)
    # Just to see the mask
    # img2 = np.where(BGmask,10,OCTscanSegmentation._BACKGROUND)
    # 		#Plot histogram of image
    # 		binSize = 0.01
    # 		bins=np.arange(0,np.amax(img2)+binSize,binSize)
    # 		hist = np.histogram(img2, bins=bins)
    # 		hfig, ax = plt.subplots(1, 1)
    # 		ax.hist(img2.flatten(),bins)
    # 		ax.set_yscale('log')
    # 		hfig.show()
    # PENDING
    # Assign intratissue BG pixels to fluid.
    # Normalize
    img2 = np.floor(255 * (img2 / np.amax(img2))).astype(int)
    img2 = np.where(BGmask == True, OCTscanSegmentation._BACKGROUND, img2)
    # Assign integers to nRetinalLayers
    r = RetinalLayers()
    nRetinalLayers = r.getLayerIndex("RPE")
    # Note that I'm not counting for choroid or fluid.
    # I can't use pixel intensity alone here for segmenting as several
    # (separated) layers exhibit similar reflectivity, so segmentation
    # criteria has to be a mix of location and intensity.
    # Estimate upper and lower boundary by checking first and last
    # non-BG pixels per column
    # UNFINISHED.
    # upperBoundary =
    # [uVals,uIdx,uInv,uCounts] = np.unique(img2, \
    # 						return_index=True, \
    # 						return_inverse=True, \
    # 						return_counts=True)
    # print(len(uVals))
    # 		## Algoritmo 2: Arlem
    # 		img2 = img
    # 		print(np.amax(img2))
    # 		tmp=np.reshape(img2,(np.size(img2),1))
    # 		print(type(tmp))
    # 		print(np.size(tmp))
    # 		print(np.amax(tmp))
    #
    #
    # 		hfig, ax = plt.subplots(1, 1)
    # 		ax.hist(tmp,255)
    # 		hfig.show()
    #
    # 		#Elimina ruido
    # 		img = segmentationUtils.ejecuta_close(img,4,4)
    #
    # 		#Amplifica capas
    # 		img = segmentationUtils.ejecuta_dilate(img,5,20,1)
    #
    # 		#Tensor
    # 		Axx, Axy, Ayy = feature.structure_tensor(img)
    #
    # 		#Elimina mas ruido
    # 		Ayy = segmentationUtils.ejecuta_close(Ayy,6,1)
    #
    # 		#Resalta las capas que sean mayores a la media
    # 		Ayy = segmentationUtils.resalta_bordes(Ayy,True,0)
    #
    # 		#Elimina aun mas ruido
    # 		Ayy = segmentationUtils.ejecuta_open(Ayy,1,1)
    #
    # 		#Binarizacion
    # 		binary = segmentationUtils.ejecuta_OTSU(Ayy)
    #
    # 		#elimina ruido del posible borde superior
    # 		binary = segmentationUtils.ejecuta_elimina_ruido_extremos(True,0,0,binary)
    #
    # 		#elimina ruido del posible borde inferior
    # 		binary = segmentationUtils.ejecuta_elimina_ruido_extremos(False,0,0,binary)
    #
    # 		#obtiene bordes exteriores
    # 		arraySuperior, arrayInferior = segmentationUtils.obten_bordes_externos(binary)
    #
    # 		#elimina ruido a la imagen original
    # 		img2 = segmentationUtils.elimina_desde_arreglos(img2, arraySuperior, arrayInferior)
    # 		img2 = segmentationUtils.ejecuta_close(img2,2,1)
    # 		img2 = feature.canny(img2,sigma = 2.5)
    # 		img2 = segmentationUtils.elimina_ruido_canny(img2,1)
    # Hide and close progress bar.
    if mode == "gui":
        progress.hide()
    # if isinstance(imgin,(OCTscan,)):
    if type(imgin) is OCTscan:
        self.result = OCTscanSegmentation(imgin)
    else:
        self.result = OCTscanSegmentation(OCTscan(imgin))
    self.result.data = img2
    return self.result


class OCTvolume(object):
    # Sphinx documentation
    """A set of :class:`data.OCTscan`
    A set of :class:`data.OCTscan`
    .. seealso:: :class:`data.OCTscan`
    .. note:: None
    .. todo:: None
    """

    # Class constructor
    def __init__(self):
        """The class constructor.
        The class constructor.
        tmp = OCTvolume() - Creates an empty volume with no scans.
        """
        self.scans = list()
        return

    # Properties getters/setters
    @property
    def scans(self):  # scans getter
        """
        The set of OCT scans.
        :getter: Gets the set of scans
        :setter: Sets the set of scans
        :type: list. All scans are of :class:`data.OCTscan`
        .. seealso:: :func:`data.OCTvolume.addScan` , :func:`data.OCTvolume.removeScan` , :func:`data.OCTvolume.clear`
        """
        return self.__scans

    @scans.setter
    def scans(self, *args):  # scans setter
        tmpScanSet = args[0]
        if type(tmpScanSet) is not list:
            warnMsg = (
                self.getClassName() + ":scans: Unexpected type. "
                "Please provide a list of data.OCTscan."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        else:
            self.__scans = list()
            for x in tmpScanSet:
                if type(x) is not OCTscan:
                    warnMsg = (
                        self.getClassName() + ":scans: Unexpected scan type "
                        "for object " + x + ". Skipping object."
                    )
                    warnings.warn(warnMsg, SyntaxWarning)
                else:
                    self.__scans.append(x)
        return None

    # Public methods
    def getClassName(self):
        return type(self).__name__

    def clear(self):
        """
        Clears the OCT volume; Removes all scans.
        :return: None
        """
        self.__scans = list()
        return None

    @deprecation.deprecated(
        deprecated_in="0.3",
        removed_in="1.0",
        #                         current_version=__version__,
        details="Use method addScans() instead.",
    )
    def addScan(self, theScan):
        """
        Add an OCT scan to the volume.
        :param theScan: The OCT scan.
        :type theScan: :class:`data.OCTscan`
        :return: None
        """
        if type(theScan) is not OCTscan:
            warnMsg = (
                self.getClassName() + ":addScan: Unexpected scan type. "
                "Nothing will be added."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        else:
            self.__scans.append(theScan)
        return None

    def addScans(self, theScans):
        """
        Add one or more OCT scans to the volume at once.
        :param theScans: The list of OCT scans.
        :type theScan: list of :class:`data.OCTscan` or
            single :class:`data.OCTscan`
        :return: True if scans added, False otherwise.
        """
        flagAllOCTscans = True
        if type(theScans) is OCTscan:
            self.__scans.append(theScans)
        elif type(theScans) is list:
            for elem in theScans:
                if type(elem) is not OCTscan:
                    warnMsg = (
                        self.getClassName() + ":addScans: Unexpected scan type. "
                        "Nothing will be added."
                    )
                    flagAllOCTscans = False
                    break
            if flagAllOCTscans:
                self.__scans.extend(theScans)
        return flagAllOCTscans

    def getScans(self, t):
        """
        Retrieves all scans in the volume of type t.
        :param t: Scan type 'A', 'B' or 'C'
        :type t: char
        :return: The set of A scans in the volume
        :rtype: list
        .. seealso:: :func:`data.OCTvolume.getVolume`
        """
        t = t.upper()  # Uppercase
        theScans = list()
        for x in self.__scans:
            if x.scantype == t:
                theScans.append(x)
        return theScans

    def getVolume(self, t):
        """
        Retrieves the (sub-)volume of scans of type t.
        :param t: Scan type 'A', 'B' or 'C'
        :type t: char
        :return: A volume with the set of scans of type t.
        :rtype: :class:`data.OCTvolume`
        .. seealso:: :func:`data.OCTvolume.getScans`
        """
        t = t.upper()  # Uppercase
        theScans = OCTvolume()
        for x in self.__scans:
            if x.scantype == t:
                theScans.addScan(x)
        return theScans

    def getNScans(self, t=None):
        """Get the number of scans of a certain type.
        If type is not given, then the total number of scans is given.
        :param t: Scan type 'A', 'B' or 'C' or None
        :type t: char
        :return: The number of scans of a certain type
        :rtype: int
        """
        if t is None:
            res = len(self.__scans)
        else:
            res = self.getNScans(self.getVolume(t))
        return res


class OCTvolumeSegmentation(object):
    # Sphinx documentation
    """A set of :class:`data.OCTscanSegmentation`
    A set of :class:`data.OCTscanSegmentation`
    .. seealso:: :class:`data.OCTscanSegmentation`
    .. note:: None
    .. todo:: None
    """

    # Class constructor
    def __init__(self):
        """The class constructor.
        The class constructor.
        tmp = OCTvolumeSegmentation() - Creates an empty segmentation volume with no scans.
        """
        self.scanSegmentations = list()
        return

    # Properties getters/setters
    @property
    def scanSegmentations(self):  # scans getter
        """
        The set of OCT scans segmentations.
        :getter: Gets the set of scans segmentations
        :setter: Sets the set of scans segmentations
        :type: list. All scans are of :class:`data.OCTscanSegmentation`
        .. seealso:: :func:`data.OCTvolumeSegmentation.addScan` , :func:`data.OCTvolumeSegmentation.removeScan` , :func:`data.OCTvolumeSegmentation.clear`
        """
        return self.__scanSegmentations

    @scanSegmentations.setter
    def scanSegmentations(self, *args):  # scanSegmentations setter
        tmpScanSet = args[0]
        if type(tmpScanSet) is not list:
            warnMsg = (
                self.getClassName() + ":scanSegmentations: Unexpected type. "
                "Please provide a list of data.OCTscanSegmentation."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        else:
            self.__scanSegmentations = list()
            for x in tmpScanSet:
                if type(x) is not OCTscanSegmentation:
                    warnMsg = (
                        self.getClassName()
                        + ":scanSegmentations: Unexpected scan type "
                        "for object " + x + ". Skipping object."
                    )
                    warnings.warn(warnMsg, SyntaxWarning)
                else:
                    self.__scanSegmentations.append(x)
        return None

    # Public methods
    def clear(self):
        """
        Clears the OCT segmentation volume; Removes all scan segmentations.
        :return: None
        """
        self.__scanSegmentations = list()
        return None

    def addScanSegmentations(self, theScanSegmentations):
        """
        Add one or multiple OCT scan segmentation to the volume.
        :param theScanSegmentations: The list of OCT scan segmentations
        :type theScanSegmentations: list of :class:`data.OCTscanSegmentation` or
            single :class:`data.OCTscanSegmentation`
        :return: None
        """
        flagAllOCTscans = False
        if type(theScanSegmentations) is OCTscanSegmentation:
            flagAllOCTscans = True
            self.__scanSegmentations.append(theScanSegmentations)
        elif type(theScanSegmentations) is list:
            flagAllOCTscans = True
            for elem in theScans:
                if type(elem) is not OCTscanSegmentation:
                    warnMsg = (
                        self.getClassName() + ":addScans: Unexpected scan type. "
                        "Nothing will be added."
                    )
                    flagAllOCTscans = False
                    break
            if flagAllOCTscans:
                self.__scanSegmentations.extend(theScanSegmentations)
        return flagAllOCTscans

    def getScanSegmentations(self, t):
        """
        Retrieves all scans in the volume of type t.
        :param t: Scan type 'A', 'B' or 'C' or scan indexing
        :type t: str, list or int
        :return: The set of scans in the volume of the chosen
        :rtype: list
        .. seealso:: :func:`data.OCTvolumeSegmentation.getVolume`
        """
        theScans = list()
        if type(t) is str:
            t = t.upper()  # Uppercase
            for x in self.__scanSegmentations:
                if x.scantype == t:
                    theScans.append(x)
        elif type(t) is list:
            for elem in t:
                theScans.append(self.__scanSegmentations[elem])
        elif type(t) is int:
            theScans.append(self.__scanSegmentations[elem])
        return theScans

    def getVolume(self, t):
        """
        Retrieves the (sub-)volume of scans of type t.
        :param t: Scan type 'A', 'B' or 'C'
        :type t: char
        :return: A volume with the set of scans of type t.
        :rtype: :class:`data.OCTvolume`
        .. seealso:: :func:`data.OCTvolume.getScans`
        """
        t = t.upper()  # Uppercase
        theScans = OCTvolumeSegmentation()
        for x in self.__scanSegmentations:
            if x.scantype == t:
                theScans.addScanSegmentation(x)
        return theScans

    def getNScans(self):
        """Get the number of scans segmentations.
        :return: The number of scans segmentations
        :rtype: int
        """
        return len(self.__scanSegmentations)


class Document:
    # Sphinx documentation
    """The document class for
    The document class for
    This class represents a document in  A document holds information
    about a study plus some additional metadata.
    Currently, a study is only an OCT image (with several scans) with or without
    segmentation information.
    .. seealso:: None
    .. note:: None
    .. todo:: None
    """

    # Private class attributes shared by all instances
    # Class constructor
    def __init__(self):
        # Call superclass constructor
        # Initialize private attributes unique to this instance
        self.study = OCTvolume()  # The current study.
        # Currently, an OCT volume
        self.segmentation = OCTvolumeSegmentation()  # The current study.
        # Currently, an OCT volumeSegmentation
        self.docsettings = Settings()
        self.docsettings.selectedScan = None
        # shared between study and segmentation
        # Document metadata
        self.folderName = "."  # Folder where the document is currently stored
        self.fileName = "OCTantDocument0001"  # The filename of the document
        self.name = "OCTantDocument0001"  # A logical name for the study
        return

    # Properties getters/setters
    #
    # Remember: Sphinx ignores docstrings on property setters so all
    # documentation for a property must be on the @property method
    @property
    def docsettings(self):  # docsettings getter
        """
        The application settings.
        :getter: Gets the document settings
        :setter: Sets the document settings
        :type: class:`data.Settings`
        """
        return self.__appsettings

    @docsettings.setter
    def docsettings(self, newSettings):  # document setter
        if newSettings is None:
            newSettings = Settings()  # Initialize settings
        if type(newSettings) is Settings:
            self.__appsettings = newSettings
        else:
            warnMsg = self.getClassName() + ":docsettings: Unexpected settings type."
            warnings.warn(warnMsg, SyntaxWarning)
        return None

    @property
    def study(self):  # study getter
        """
        The OCT volume being processed and analysed.
        ..todo: Upgrade to volume. Watch out! This will affect many
            other classess using this method.
        :getter: Gets the OCT volume.
        :setter: Sets the OCT volume.
        :type: :class:`data.OCTvolume`
        """
        return self.__study

    @study.setter
    def study(self, vol):  # study setter
        if vol is None or type(vol) is OCTvolume:
            self.__study = vol
            # ...and reset scan
            self.segmentedScan = None
        elif type(vol) is OCTscan:
            warnMsg = (
                self.getClassName()
                + ":study: OCTvolume expected but OCTscan received. Embedding scan."
            )
            warnings.warn(warnMsg, SyntaxWarning)
            self.__study = OCTvolume()
            self.__study.addScan(vol)
        else:
            warnMsg = self.getClassName() + ":study: Unexpected study type."
            warnings.warn(warnMsg, SyntaxWarning)
        return None

    @property
    def segmentation(self):  # segmentation getter
        """
        The segmentation over the OCT study being processed and analysed.
        :getter: Gets the OCT volume segmentation.
        :setter: Sets the OCT volume segmentation.
        :type: :class:`data.OCTvolumeSegmentation`
        """
        return self.__segmentation

    @segmentation.setter
    def segmentation(self, newSegmentation):  # segmentation setter
        if (newSegmentation is None) or (
            type(newSegmentation) is OCTvolumeSegmentation
        ):
            self.__segmentation = newSegmentation
            if newSegmentation is not None:
                if self.study is None:
                    warnMsg = self.getClassName() + ":segmentation: No reference image."
                    warnings.warn(warnMsg, SyntaxWarning)
                if not (
                    len(newSegmentation.scanSegmentations) == len(self.study.scans)
                ):
                    warnMsg = self.getClassName() + ":segmentation: Unexpected size."
                    warnings.warn(warnMsg, SyntaxWarning)
        elif type(newSegmentation) is OCTscanSegmentation:
            warnMsg = (
                self.getClassName()
                + ":study: OCTvolumeSegmentation expected but OCTscanSegmentation received. Embedding scan."
            )
            warnings.warn(warnMsg, SyntaxWarning)
            self.__segmentation = OCTvolumeSegmentation()
            self.__segmentation.addScanSegmentations(newSegmentation)
        else:
            warnMsg = (
                self.getClassName() + ":segmentation: Unexpected segmented scan type."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        return None

    @property
    def folderName(self):  # folderName getter
        """
        Folder where the document is currently stored.
        ..note: Also retrieve the py:attr:`fileName` to build the full path.
        :getter: Gets the study folder name.
        :setter: Sets the study folder name. If new folder is None,
            the current directory '.' is chosen.
        :type: str
        """
        return self.__folderName

    @folderName.setter
    def folderName(self, d):  # name setter
        if d is None:
            d = "."  # Set to current folder
        if type(d) is str:
            self.__folderName = d
        else:
            warnMsg = self.getClassName() + ":name: Unexpected folderName type."
            warnings.warn(warnMsg, SyntaxWarning)
        return None

    @property
    def fileName(self):  # fileName getter
        """
        The filename of the document.
        ..note: Also retrieve the py:attr:`folderName` to build the full path.
        :getter: Gets the the filename of the document.
        :setter: Sets the The filename of the document. If new name is None,
            a default name is given.
        :type: str
        """
        return self.__folderName

    @fileName.setter
    def fileName(self, newFilename):  # fileName setter
        if newFilename is None:
            newFilename = self.defaultName  # Set to default name
        if type(newFilename) is str:
            self.__fileName = newFilename
        else:
            warnMsg = self.getClassName() + ":name: Unexpected fileName type."
            warnings.warn(warnMsg, SyntaxWarning)
        return None

    @property
    def name(self):  # name getter
        """
        A logical name for the study.
        :getter: Gets the OCT study name.
        :setter: Sets the OCT study name.
        :type: str
        """
        return self.__name

    @name.setter
    def name(self, newName):  # name setter
        if newName is None or type(newName) is str:
            self.__name = newName
        else:
            warnMsg = self.getClassName() + ":name: Unexpected name type."
            warnings.warn(warnMsg, SyntaxWarning)
        return None

    # Private methods
    # Public methods
    def getClassName(self):
        return type(self).__name__

    def getCurrentScan(self):
        """Get the current working OCT scan
        Change the current selection using :func:`pickScan`
        :returns: The current working OCT scan.
        :rtype: :class:`data.OCTscan` or None if the study contains no scans
        """
        if self.docsettings.selectedScan is None:
            self.docsettings.selectedScan = 0
        res = None
        if len(self.__study.scans) > 0:
            res = self.__study.scans[self.docsettings.selectedScan]
        return res

    def setCurrentScan(self, newScan):
        """Sets the current working OCT scan
        Change the current selection using :func:`pickScan`
        :param newScan: An OCT scan to be assigned to the current working OCT scan.
        :type newScan: :class:`data.OCTscan`
        """
        if self.docsettings.selectedScan is None:
            self.docsettings.selectedScan = 0
        if newScan is None:
            if self.__study.getNScans() == 0:
                # do nothing.
                pass
            else:
                warnMsg = (
                    self.getClassName()
                    + ":setcurrentscan: Unexpected scan type NoneType."
                )
                warnings.warn(warnMsg, SyntaxWarning)
        if type(newScan) is OCTscan:
            self.__study.scans[self.docsettings.selectedScan] = newScan
        else:
            warnMsg = self.getClassName() + ":setcurrentscan: Unexpected scan type."
            warnings.warn(warnMsg, SyntaxWarning)
        return None

    def pickScan(self, i):
        """Pick the i-th OCT scan (and its segmentation) for working.
        Sets the docsetting.selectedScan to i checking that it does exist.
        :param i: The selected scan index
        :type i: int
        :return: None
        """
        if type(i) is int and i < self.study.getNScans() and i >= 0:
            self.docsettings.selectedScan = i
        else:
            warnMsg = self.getClassName() + ":pickScan: Selected scan does not exist."
            warnings.warn(warnMsg, SyntaxWarning)
        return None

    def getCurrentScanSegmentation(self):
        """Get the current working OCT scanSegmentation
        Change the current selection using :func:`pickScan`
        :returns: The current working OCT scanSegmentation.
        :rtype: :class:`data.OCTscanSegmentation`
        """
        if self.docsettings.selectedScan is None:
            self.docsettings.selectedScan = 0
        res = None
        if len(self.__segmentation.scanSegmentations) > 0:
            res = self.__segmentation.scanSegmentations[self.docsettings.selectedScan]
        return res

    def setCurrentScanSegmentation(self, newScan):
        """Sets the current working OCT scanSegmentation
        Change the current selection using :func:`pickScan`
        :param newScan: An OCT scan to be assigned to the current working OCT scan.
        :type newScan: :class:`data.OCTscanSegmentation`
        """
        if self.docsettings.selectedScan is None:
            self.docsettings.selectedScan = 0
        if newScan is None:
            if self.__segmentation.getNScans() == 0:
                # do nothing.
                pass
            else:
                warnMsg = (
                    self.getClassName()
                    + ":setcurrentscansegmentation: Unexpected scan type NoneType."
                )
                warnings.warn(warnMsg, SyntaxWarning)
        if type(newScan) is OCTscanSegmentation:
            self.__segmentation.scanSegmentations[
                self.docsettings.selectedScan
            ] = newScan
        else:
            warnMsg = (
                self.getClassName()
                + ":setcurrentscansegmentation: Unexpected scan type."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        return None

    def readFile(self, filename):
        """Reads an OCTant document file.
        This method is currently a sham, and it will be updated
        when serialization is incorporated to  Currently,
        it returns an empty document. Nevertheless, it already
        updates the document, clearing all fields to default values,
        and updates the filename and folder
        The file must exist or an error is generated.
        The file must be in OCTant file format.
        :param fileName: The file name
        :type fileName: str
        :return: This document
        :rtype: :class:`data.Document`
        """
        self = Document()
        self.folderName, self.fileName = os.path.split(filename)
        return self


class Settings(object):
    # Sphinx documentation
    """A class to hold a list of settings.
    A class to hold a list of settings. A list of settings is the classical
    "dictionary" (see note below on python's dict) of pairs key:value
    but with some additional capabilities. In particular, the class provides
    additional file reading and writing capabilities so that settings can
    be read to and from plain text files, as well as some value setting
    checking capabilities.
    The class is intended to behave like a dynamic struct
    where properties of the class, i.e. new settings, can be declared
    "on-the-fly" instead of being predefined.
    Although, creating a dynamic struct class in python itself is trivial
    (see https://stackoverflow.com/questions/1878710/struct-objects-in-python ),
    but because of the additional capabilities, hence the convenience of
    the class.
    .. Background:
    MATLAB's struct allows this "on-the-fly" field declaration on the fly.
    Python's built-in dictionary is not exactly a match because of the
    required syntax i.e. ``mySettingsObj['fieldname']`` instead of
    ``mySettingsObj.fieldname`` and the inability to control value settings.
    .. seealso:: None
    .. note:: None
    .. todo:: None
    """

    # Private class attributes shared by all instances
    # Class constructor
    def __init__(self, **kwargs):
        """The class constructor."""
        # Call superclass constructor
        # Initialize private attributes unique to this instance
        self.__dict__.update(kwargs)  # Permits direct declaration of
        # key:value pairs from declaration, e.g.
        # x = Settings(foo=1, bar=2)
        # See: https://stackoverflow.com/questions/1878710/struct-objects-in-python
        return

    # Properties getters/setters
    #
    # Remember: Sphinx ignores docstrings on property setters so all
    # documentation for a property must be on the @property method
    # Private methods
    def __str__(self):
        """Provides the string representation of the object.
        :returns: The object representation as a string.
        :rtype: str
        """
        s = "<" + self.getClassName() + ": "
        s = s + str(self.__dict__)
        s = s + ">"
        return s

    # Public methods
    def getClassName(self):
        """Retrieves the class name.
        :returns: The class name.
        :rtype: str
        """
        return type(self).__name__

    def read(self, filename):
        """Read settings from JSON file.
        :param filename: The name of the file to be read (including path)
        :type filename: str
        :returns: True if file was sucessfully read. False otherwise.
        :rtype: bool
        """
        with open(filename, "r") as file:
            contentStr = file.read()
        # c=json.loads(jsonminify(content))
        contentDict = json.loads(fastjsonminify(contentStr))
        # c contains a dictionary that has to be
        # traspassed to self
        # Loop over the dictionary
        for fieldName, fieldValue in contentDict.items():
            setattr(self, fieldName, fieldValue)
        return True

    def write(self, filename):
        """Write settings to a JSON file.
        :returns: True if file was sucessfully read. False otherwise.
        :rtype: bool
        """
        contentStr = json.dumps(self.__dict__)
        with open(filename, "w") as file:
            file.write("# \n")
            file.write("# File: " + filename + "\n")
            file.write("# \n")
            file.write("# This is an OCTant settings file.\n")
            file.write("# You can add, edit or remove settings manually here.\n")
            file.write(
                "# File format is in JSON. Although comments are permitted, but they will be lost after resaving because of minification.\n"
            )
            file.write(
                '# If you want your comments to be persistent declared them as "__comment" fields.\n'
            )
            file.write("# \n")
            file.write(
                "# File last saved: "
                + datetime.utcnow().strftime("%d-%b-%Y %H:%M:%S UTC+0")
                + "\n"
            )
            file.write("# \n")
            file.write("# (c) 2019.  Felipe Orihuela-Espina.\n")
            file.write("# \n\n")
            file.write(contentStr)
        return True


class OpScanFlatten(Operation):
    """A flattening operation for :class:`data.OCTscan`.
    A flattening operation for :class:`data.OCTscan`.
    The operation represented by this class rectifies an OCT scan.
    .. seealso:: None
    .. note:: None
    .. todo:: None
    """

    # Private class attributes shared by all instances
    # Class constructor
    def __init__(self):
        # Call superclass constructor
        super().__init__()
        # Set the operation name
        self.name = "Flattening"
        self.__deformationMap = None
        return

    @property
    def deformationMap(self):  # name getter
        """
        A logical name for the study.
        This is a read only property.
        :getter: Gets the deformationMap associated to the last flattening.
        :type: str
        """
        return self.__deformationMap

    # Private methods
    # def __str__(self):
    #     #This not working yet; it gets into an infiite recursion as
    #     #super().__str__() calls self.getClassName() in THIS class.
    #     # s = '<' + self.getClassName() + '([' \
    #     #     + super().__str__() + '])>'
    #     s = '<' + self.getClassName() + '([' \
    #             + str(super()) + '])>'
    #     print(super())
    #     return s
    @staticmethod
    def fittingQuadraticModel(x, a, b, c):
        # quadratic model for curve optimization
        return a * x * x + b * x + c

    # Public methods
    def execute(self, *args, **kwargs):
        """Executes the operation on the :py:attr:`operands`.
        Executes the operation on the :py:attr:`operands` and stores the outcome
        in :py:attr:`result`. Preload operands using
        :func:`Operation.addOperand()`.
        :returns: Result of executing the operation.
        :rtype: :class:`data.OCTscan`
        """
        # print(self._getClasName(),": flattening: Starting flattening")
        # Ensure the operand has been set.
        if len(self.operands) < 1:
            warnMsg = self.getClassName() + ":execute: Operand not set."
            warnings.warn(warnMsg, SyntaxWarning)
            return None
        imgin = self.operands[0]
        if type(imgin) is OCTscan:
            imgin = imgin.data
        # Check whether the image is in RGB (ndim=3) or in grayscale (ndim=2)
        # and convert to grayscale if necessary
        if imgin.ndim == 2:
            # Dimensions are only width and height. The image is already in grayscale.
            I2 = imgin
        elif imgin.ndim == 3:
            # Image is in RGB. Convert.
            I2 = color.rgb2gray(imgin)
        else:  # Unexpected case. Return warning
            print(self._getClasName(), ": Unexpected image shape.")
            self.result = imgin
            return self.result
        aux = np.argmax(I2, axis=0)
        mg = np.mean(aux)
        sdg = np.std(aux)
        markers = []
        remover = []
        x0 = np.arange(len(aux))
        for i in range(0, len(aux)):
            if mg - 3 * sdg <= aux[i] <= mg + 3 * sdg:
                markers += [aux[i]]
            else:
                remover += [i]
        x = np.delete(x0, remover)
        modelCoeffs, pcov = curve_fit(
            self.fittingQuadraticModel, x, markers, method="dogbox", loss="soft_l1"
        )
        a = self.fittingQuadraticModel(x0, *modelCoeffs)
        shift = np.max(a)
        flat = shift - a
        flat = np.round(flat)
        flat = np.ravel(flat).astype(int)
        self.__deformationMap = flat
        newgray = I2
        for i in range(0, len(a)):
            newgray[:, i] = np.roll(I2[:, i], flat[i], axis=0)
        self.result = OCTscan(newgray)
        return self.result

    #    #@deprecated(version='0.2', reason="Deprecated. Use method execute() instead.")
    #    @deprecation.deprecated(deprecated_in="0.2", removed_in="0.3",
    #                        current_version=__version__,
    #                        details="Use method execute() instead.")
    #    def flattening(self,image):
    #        #Encapsulate the image as an OCTscan
    #        tmp=OCTscan(image)
    #        self.clear()
    #        self.addOperand(tmp)
    #        #Execute
    #        self.execute()
    def applyOperation(self, scanA):
        """Apply the current flattening to the given scan.
        Instead of calculating the fitting again needed for the
        flattening, this method applies a known fitted quadratic model to
        the given parameters.
        The result is NOT stored in :py:attr:`result`.
        :param scanA: Image to flatten.
        :type scanA: :class:`data.OCTscan`
        :returns: Result of repeating the last flattening operation onto
             parameter scanA.
        :rtype: :class:`data.OCTscan`
        """
        if type(scanA) is OCTscan:
            scanA = scanA.data
        newgray = scanA
        for i in range(0, len(self.deformationMap)):
            newgray[:, i] = np.roll(scanA[:, i], self.deformationMap[i], axis=0)
        return OCTscan(newgray)


import warnings
from abc import ABC, abstractmethod


class Operation(ABC):
    # Sphinx documentation
    """An abstract base class for operations on :class:`octant.data.OCTvolume` and :class:`octant.data_OCTscan` .
    An abstract base class for operations on :class:`octant.data.OCTvolume` and
    :class:`octant.data.OCTscan`.
    :Example:
        tmp = octant.data.OCTscan(img)
        o = octant.op.OpScanFlattening()
        o.addOperand(tmp)
        o.arity() #return 1
        o.execute() #Flattens the image
    :Known subclasses:
    * :class:`OpScanFlattening`
    * :class:`OpScanMeasureLayerThickness`
    * :class:`OpScanPerfilometer`
    * :class:`OpScanSegment`
    * :class:`OpScanStitch`
    * :class:`OpSegmentationBrush`
    * :class:`OpSegmentationEdit`
    .. seealso:: None
    .. note:: None
    .. todo::
        * Support to hold operand names.
    """

    # Private class attributes shared by all instances
    # Class constructor
    #
    def __init__(self, **kwargs):
        """The class constructor.
        The class constructor. Creates an empty operation
        """
        super().__init__()
        # Initialize attributes (without decorator @property)
        # Initialize properties (with decorator @property)
        self.name = "Operation"  # The operation name
        self.operands = list()  # Operands
        self.parameters = list()  # Parameters
        self.result = None  # Operation outputs (a list in case it is multivalued).
        # None until executed or cleared.
        if kwargs is not None:
            for key, value in kwargs.items():
                if key == "name":
                    self.name = value
        return

    # Properties getters/setters
    #
    # Remember: Sphinx ignores docstrings on property setters so all
    # documentation for a property must be on the @property method
    @property
    def operands(self):  # operands getter
        """
        The list of operands.
        :getter: Gets the list of operands
        :setter: Sets the list of operands.
        :type: list
        """
        return self.__operands

    @operands.setter
    def operands(self, opList):  # operands setter
        # if (not isinstance(opList,(list,))):
        if type(opList) is not list:
            warnMsg = (
                self.getClassName() + ":operands: Unexpected type. "
                "Please provide operands as a list."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        else:
            self.__operands = opList
        return None

    @property
    def name(self):  # name getter
        """
        The operation name
        :getter: Gets the operation name
        :setter: Sets the operation name.
        :type: string
        """
        return self.__name

    @name.setter
    def name(self, opName):  # name setter
        # if (not isinstance(opName,(str,))):
        if type(opName) is not str:
            warnMsg = (
                self.getClassName() + ":name: Unexpected type. "
                "Operations name must be a string."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        else:
            self.__name = opName
        return None

    @property
    def parameters(self):  # operands getter
        """
        The list of parameters.
        :getter: Gets the list of parameters
        :setter: Sets the list of parameters.
        :type: list
        """
        return self.__parameters

    @parameters.setter
    def parameters(self, opList):  # operands setter
        # if (not isinstance(opList,(list,))):
        if type(opList) is not list:
            warnMsg = (
                self.getClassName() + ":parameters: Unexpected type. "
                "Please provide operands as a list."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        else:
            self.__parameters = opList
        return None

    @property
    def result(self):  # result getter
        """
        The list of results.
        This is a read only property. There is no setter method.
        :getter: Gets the list of results
        :setter: Sets the list of results
        :type: list
        """
        return self.__result

    @result.setter
    def result(self, rList):  # result setter
        self.__result = rList
        return None

    # Private methods
    def __str__(self):
        tmp = "["
        for x in self.operands:
            tmp += format(x) + ","
        tmp += "]"
        s = (
            "<"
            + self.getClassName()
            + "(["
            + "name: "
            + self.name
            + ";"
            + " operands: "
            + tmp
            + "])>"
        )
        return s

    # Public methods
    def getClassName(self):
        """Get the class name as a string.
        Get the class name as a string.
        :returns: The class name.
        :rtype: string
        """
        return type(self).__name__

    def addOperand(self, op, i=None):
        """
        Add a new operand.
        :param op: The operand.
        :type op: object
        :param i: (optional) The operand order. If given it may shift the
            order of other operands already set. If not given, the operand
            is appended at the end of the list of operands.
        :type op: int
        :return: None
        """
        if i is None:
            self.__operands.append(op)
        else:
            self.__operands.insert(i, op)
        return None

    def setOperand(self, op, i):
        """
        Set an operand; substitutes an existing operand with a new one.
        Calling setOperand when the :py:attr:`i`-th operand has not been
        previously set will result in an out-of-range error.
        :param op: The new operand.
        :type op: object
        :param i: The operand order. Operand index is zero-base i.e. the
            first operand occupies i=0
        :type op: int
        :return: None
        """
        self.__operands[i] = op
        return None

    def addParameter(self, param, i=None):
        """
        Add a new parameter.
        :param op: The parameter.
        :type op: object
        :param i: (optional) The paremeter order. If given it may shift the
            order of other parameters already set. If not given, the parameter
            is appended at the end of the list of parameters.
        :type op: int
        :return: None
        """
        if i is None:
            self.__parameters.append(op)
        else:
            self.__parameters.insert(i, op)
        return None

    def setParameter(self, op, i):
        """
        Set a parameter; substitutes an existing parameter with a new one.
        Calling setParameter when the :py:attr:`i`-th parameter has not been
        previously set will result in an out-of-range error.
        :param op: The new operand.
        :type op: object
        :param i: The operand order. Operand index is zero-base i.e. the
            first operand occupies i=0
        :type op: int
        :return: None
        """
        self.__operands[i] = op
        return None

    def arity(self):
        """Gets the operation arity (number of operands).
        :return: The operation arity
        :rtype: int
        """
        return len(self.__operands)

    def clear(self):
        """
        Clears the operands; Removes all operands.
        :return: None
        """
        self.__operands = list()
        return None

    # @abstractmethod
    def execute(self, *args, **kwargs):
        """Executes the operation on the operands.
        This is an abstract method. Executes the operation on the .operands
        and stores the outcome in .result
        Operation meta-parameters may be also passed.
        :returns: Result of executing the operation.
        :rtype: Type of result -depends on subclass implementation-.
        """
        pass
<|endoftext|>
[{"Id": 87643680, "ScriptId": 24617507, "ParentScriptVersionId": NaN, "ScriptLanguageId": 9, "AuthorUserId": 359577, "CreationDate": "02/12/2022 13:18:33", "VersionNumber": 1.0, "Title": "OCT segmentation", "EvaluationDate": "02/12/2022", "IsChange": true, "TotalLines": 2855.0, "LinesInsertedFromPrevious": 2855.0, "LinesChangedFromPrevious": 0.0, "LinesUnchangedFromPrevious": 0.0, "LinesInsertedFromFork": NaN, "LinesDeletedFromFork": NaN, "LinesChangedFromFork": NaN, "LinesUnchangedFromFork": NaN, "TotalVotes": 0}]
[{"Id": 117988798, "KernelVersionId": 87643680, "SourceDatasetVersionId": 3175055}]
[{"Id": 3175055, "DatasetId": 1929287, "DatasourceVersionId": 3224540, "CreatorUserId": 359577, "LicenseName": "Unknown", "CreationDate": "02/12/2022 13:17:15", "VersionNumber": 2.0, "Title": "OCTant project", "Slug": "octant-project", "Subtitle": "Retinal imaging data", "Description": NaN, "VersionNotes": "Data Update 2022/02/12", "TotalCompressedBytes": 0.0, "TotalUncompressedBytes": 0.0}]
[{"Id": 1929287, "CreatorUserId": 359577, "OwnerUserId": 359577.0, "OwnerOrganizationId": NaN, "CurrentDatasetVersionId": 3199639.0, "CurrentDatasourceVersionId": 3249374.0, "ForumId": 1952953, "Type": 2, "CreationDate": "02/12/2022 11:00:31", "LastActivityDate": "02/12/2022", "TotalViews": 998, "TotalDownloads": 4, "TotalVotes": 2, "TotalKernels": 1}]
# # Enhancing an open-source OCT segmentation tool with manual segmentation capabilities
# ## Introduction
# Optical coherence tomography (OCT) is a way of capturing images for the diagnosis and monitoring of the retina. The retina is complex tissue composed of ten major layers of distinct cells. An important part of the routine use of OCT retinal images is the segmentation of the tissue layers. The segmentation consists of labelling the image pixels with their corresponding layer and remains a complex problem despite the existence of some approximate solutions.
# Around the central task of segmentation, additional functions need to be performed on the image, for instance to quantify the thickness of the tissue layers, or to stitch several images together to see a wider region of the eye. For ophthalmologists to be able to use the algorithms needed to process the OCT image, these functions are often encapsulated in user friendly applications. But alas, software to segment retinal OCT images is expensive and not always available in developing countries and/or economically disfavored areas.
# OCTant (https://github.com/forihuelaespina/OCTant/) is a python based free open-source software for the analysis of OCT images originally proposed by the National Institute of Astrophysics, Optics and Electronics (INAOE in Spanish) in Mexico but which is still under development. In its current version, OCTant can flatten or mosaic the images, and measure the thickness of each layer.
# Because the performance of even state-of-the-art segmentation algorithms is always suboptimal, ophthalmologists often benefit from being able to make manual adjustments. This project will improve the OCTant tool with the capacity to manually adjust a given segmentation
# The project will teach you:
# * Basic medical image processing.
# * To develop code valid within a collaborative project (i.e. several programming rules ought to be respected so that other collaborators can use and understand you code off-the-shelf. For instance, you will familiarize yourself with GitHub versioning platform.
# * To interact with software in development, where bugs are likely to exist, and errors are likely to appear. For instance, you will learn to report bugs.
# * To produce technical documentation in Sphynx.
# To start off with, let's take a look at an example image:
# import os
# import glob
# import numpy as np
from PIL import Image

# import matplotlib.pyplot as plt
img = Image.open("../input/octant-project/_1323904_.jpeg")
img
# what can we do to the image?
dir(img)
print(img.width, img.height)
spinPicture = img.rotate(30)
spinPicture
from PIL import Image
import matplotlib.pyplot as plt


def getRed(redVal):
    return "#%02x%02x%02x" % (redVal, 0, 0)


def getGreen(greenVal):
    return "#%02x%02x%02x" % (0, greenVal, 0)


def getBlue(blueVal):
    return "#%02x%02x%02x" % (0, 0, blueVal)


# Create an Image with specific RGB value
image = img
# Display the image
image.show()
# Get the color histogram of the image
histogram = image.histogram()
# Take only the Red counts
l1 = histogram[0:256]
# Take only the Blue counts
l2 = histogram[256:512]
# Take only the Green counts
l3 = histogram[512:768]
plt.figure(0)
# R histogram
for i in range(0, 256):
    plt.bar(i, l1[i], color=getRed(i), edgecolor=getRed(i), alpha=0.3)
# G histogram
plt.figure(1)
for i in range(0, 256):
    plt.bar(i, l2[i], color=getGreen(i), edgecolor=getGreen(i), alpha=0.3)
# B histogram
plt.figure(2)
for i in range(0, 256):
    plt.bar(i, l3[i], color=getBlue(i), edgecolor=getBlue(i), alpha=0.3)
plt.show()
# %matplotlib inline
# from matplotlib import pyplot as plt
# from gluoncv import model_zoo, data, utils
# net = model_zoo.get_model('mask_rcnn_fpn_resnet101_v1d_coco', pretrained=True)
# from PIL import Image
# import pandas as pd
# import os
# import numpy as np
# pth = '../input/octant-project/_1323904_.jpeg'
# x, orig_img = data.transforms.presets.rcnn.load_test(pth)
# ids, scores, bboxes, masks = [xx[0].asnumpy() for xx in net(x)]
# width, height = orig_img.shape[1], orig_img.shape[0]
# masks = utils.viz.expand_mask(masks, bboxes, (width, height), scores)
# orig_img = utils.viz.plot_mask(orig_img, masks)
# # identical to Faster RCNN object detection
# fig = plt.figure(figsize=(20, 20))
# ax = fig.add_subplot(1, 1, 1)
# ax = utils.viz.plot_bbox(orig_img, bboxes, scores, ids,
#                          class_names=net.classes, ax=ax)
# plt.show()
from skimage import io
import copy
import octant as oc

fileName = "../input/octant-project/_1323904_.jpeg"
img = io.imread(fileName)
scan = OCTscan(img)
scanSegmentation = OCTscanSegmentation(scan)
print("Creating document.")
study = OCTvolume()  # Initialize a document
study.addScans(scan)
segmentationVol = OCTvolumeSegmentation()  # Initialize a document
segmentationVol.addScanSegmentations(scanSegmentation)
doc = Document()  # Initialize a document
doc.name = fileName
tmp, _ = os.path.split(fileName)
doc.folderName = tmp
doc.fileName = fileName
doc.study = study
doc.segmentation = segmentationVol
# Keep reference image.
print("Replicating image.")
doc2 = copy.deepcopy(doc)
# Flattening
print("-- Flattening.")
flt = OpScanFlatten()
flt.addOperand(doc2.getCurrentScan())
imFlattened = flt.execute()
doc2.setCurrentScan(imFlattened)
# Segmentation
print("-- Segmenting.")
doc3 = copy.deepcopy(doc2)
seg = OpScanSegment()
seg.addOperand(doc3.getCurrentScan())
imSegmented = seg.execute()
doc3.setCurrentScanSegmentation(imSegmented)
# Load colormap
print("-- Plotting.")
appsettings = Settings()
appsettingsfile = "..\\resources\\OCTantApp.config"
appsettings.read(appsettingsfile)
cmap = appsettings.retinallayerscolormap
# hFig = myInitFigure()
# myPaint(hFig,doc) #Plot raw with dummy segmentation
# hFig = myPaint(myInitFigure(),doc) #Plot raw
hFig = myPaint(myInitFigure(), doc2)  # Plot flattened
hFig = myPaint(myInitFigure(), doc3)  # Plot segmented


# imgin = OCTscan(img)
class OCTscan(object):
    # Sphinx documentation
    """A single OCT scan.
    A single OCT scan. A scan is a grayscale image.
    Check .scantype property for the scan type (A, B or C)
    .. seealso:: :class:`data.OCTvolume`
    .. note:: None
    .. todo:: Upgrade to color scans.
    """

    # Class constructor
    def __init__(self, *args):
        """The class constructor.
        The class constructor.
        tmp = OCTscan() - Creates an black scan sized 480x640.
        tmp = OCTScan(img) - Creates scan from the (grayscale) image. Assumed to be an A scan.
        tmp = OCTScan(img,type) - Creates scan from the (grayscale) image.
        :param img: The scan image
        :type img: numpy.ndarray
        :param type: The scan type ('A' -default-, 'B' or 'C')
        :type type: char
        """
        if len(args) > 1:
            warnMsg = (
                self.getClassName() + ":__init__: Unexpected number of input arguments."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        # Initialize attributes (without decorator @property)
        # Initialize properties (with decorator @property)
        self.data = np.zeros(shape=(480, 640), dtype=np.uint8)  # The default scan image
        self.scantype = "A"
        if len(args) > 0:
            self.data = args[0]
        if len(args) > 1:
            self.scantype = args[1]
        return

    # Properties getters/setters
    #
    # Remember: Sphinx ignores docstrings on property setters so all
    # documentation for a property must be on the @property method
    @property
    def data(self):  # data getter
        """
        The (grayscale) scan image. The image is expected to be
        a grayscale image. Colour images will be converted to grayscale.
        :getter: Gets the OCT scan image
        :setter: Sets the OCT scan image.
        :type: numpy.ndarray shaped [width,height]
        """
        return self.__data

    @data.setter
    def data(self, img):  # data setter
        if img is not None:
            # Check whether the image is in RGB (ndim=3) or in grayscale (ndim=2)
            # and convert to grayscale if necessary
            if img.ndim == 2:
                # Dimensions are only width and height. The image is already in grayscale.
                pass
            elif img.ndim == 3:
                # Image is in RGB. Convert.
                img = color.rgb2gray(img)
            else:  # Unexpected case. Return warning
                warnMsg = self.getClassName() + ":data: Unexpected image shape."
                warnings.warn(warnMsg, SyntaxWarning)
            self.__data = img
        return None

    @property
    def shape(self):  # shape getter
        """
        The scan shape [width,height].
        :getter: Gets the scan shape
        :setter: None. This is a read-only property.
        :type: Tuple [width,height]
        """
        return self.__data.shape

    @shape.setter
    def shape(self, *args):  # shape setter
        # Catching attempts to set the shape of the scan
        warnMsg = self.getClassName() + ":shape: shape is a read-only property."
        warnings.warn(warnMsg, UserWarning)
        return

    @property
    def scantype(self):  # scantype getter
        """
        The scan type; 'A', 'B' or 'C'.
        :getter: Gets the scan type.
        :setter: Sets the scan type.
        :type: char 'A', 'B' or 'C'.
        """
        return self.__scantype

    @scantype.setter
    def scantype(self, *args):  # shape setter
        stype = args[0].upper()  # Uppercase
        if not stype in ("A", "B", "C"):
            warnMsg = (
                self.getClassName() + ":scantype: Scan type can only be "
                "A"
                ", "
                "B"
                " or "
                "C"
                "."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        self.__scantype = stype
        return None

    # Private methods
    def __str__(self):
        s = (
            "<"
            + self.getClassName()
            + "(["
            + "  data: "
            + format(self.data)
            + ","
            + "  shape: "
            + format(self.shape)
            + ","
            + "  scantype: "
            + self.scantype
            + "])>"
        )
        return s

    # Public methods
    def getClassName(self):
        """Get the class name as a string.
        Get the class name as a string.
        :returns: The class name.
        :rtype: string
        """
        return type(self).__name__


class OCTscanSegmentation(object):
    # Sphinx documentation
    """A retinal layer segmentation over a :class:`data.OCTscan`
    A retinal layer segmentation over a :class:`data.OCTscan`. A segmentation
    assigns every pixel of the scan a class label.
    Please note that this is a data model class; it keeps the segmentation
    but it is NOT capable of "computing" such segmentation. To compute a
    segmentation please refer to :class:`op.OpScanSegment`.
    The segmentation is sized and shaped equal to its base
    :class:`data.OCTscan`.
    A default segmentation sets the whole segmentation to BACKGROUND.
    .. seealso:: :class:`data.OCTscan`, :class:`op.OpScanSegment`
    .. note:: None
    .. todo:: None
    """
    _BACKGROUND = 0  # The background label identifier

    # Class constructor
    def __init__(self, *args):
        """The class constructor.
        The class constructor.
        tmp = OCTscanSegmentation(theOCTScan) - Creates a default
            segmentation for the given :class:`data.OCTscan`
        :param theOCTScan: The OCT scan to be segmented
        :type img: :class:`data.OCTscan`
        """
        refImage = OCTscan()
        # Dummy reference
        if len(args) == 0:
            warnMsg = (
                self.getClassName()
                + ":__init__: Unexpected number of input arguments. Generating a dummy reference scan."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        else:
            refImage = args[0]
            # if type(refImage) is not OCTscan:
            #     raise ErrorValue #Throw error
        # Initialize attributes (without decorator @property)
        # Initialize properties (with decorator @property)
        self.scan = refImage  # The OCT scan over which the segmentation is made
        self.data = self._BACKGROUND * np.ones(
            refImage.shape
        )  # The segmentation itself
        self.classMap = RetinalLayers().layers  # The map of class labels
        return

    # Properties getters/setters
    #
    # Remember: Sphinx ignores docstrings on property setters so all
    # documentation for a property must be on the @property method
    @property
    def data(self):  # data getter
        """
        The segmentation labels map. Please refer to :py:attr:`classMap` for
        classes.
        ..note: WARNING! This method is not currently checking whether the
            data is sized equal to the scan. This may become a problem later.
            The problem is that trying to check scan.shape will raise an
            error during object creation, when attemting to set the data
            but because the object has not been created yet, it still lacks
            the scan property even if declared in advance.
        :getter: Gets the segmentation map
        :setter: Sets the segmentation map
        :type: numpy.ndarray shaped [width,height]
        """
        return self.__data

    @data.setter
    def data(self, segmentedImg):  # data setter
        self.__data = segmentedImg
        # if segmentedImg is not None:
        #     #Check whether the image is in RGB (ndim=3) or in grayscale (ndim=2)
        #     #and convert to grayscale if necessary
        #     if ((segmentedImg.ndim == 2) & (segmentedImg.shape == self.scan.shape)):
        #         #Dimensions are only width and height, and matches that of
        #         #the scan.
        #         self.__data = segmentedImg;
        #     else: #Unexpected case. Return warning
        #         warnMsg = self.getClassName() + ':data: Unexpected segmentation shape.'
        #         warnings.warn(warnMsg,SyntaxWarning)
        return None

    @property
    def scan(self):  # scan getter
        """
        The base OCT scan. Please refer to :py:attr:`data` for
        the segmentation map.
        :getter: Gets the base OCT scan
        :setter: Sets the base OCT scan
        :type: :class:`data.OCTscan`
        """
        return self.__scan

    @scan.setter
    def scan(self, octScan):  # scan setter
        if octScan is not None:
            # Check whether the image is in RGB (ndim=3) or in grayscale (ndim=2)
            # and convert to grayscale if necessary
            if type(octScan) is OCTscan:
                # Dimensions are only width and height, and matches that of
                # the scan.
                self.__scan = octScan
                self.clear()
            else:  # Unexpected case. Return warning
                warnMsg = self.getClassName() + ":data: Unexpected type for OCT scan."
                warnings.warn(warnMsg, SyntaxWarning)
        return None

    @property
    def shape(self):  # shape getter
        """
        The scan segmentation shape [width,height].
        :getter: Gets the scan segmentation shape
        :setter: None. This is a read-only property.
        :type: Tuple [width,height]
        """
        return self.__data.shape

    @shape.setter
    def shape(self, *args):  # shape setter
        # Catching attempts to set the shape of the scan
        warnMsg = self.getClassName() + ":shape: shape is a read-only property."
        warnings.warn(warnMsg, UserWarning)
        return

    @property
    def classMap(self):  # classMap getter
        """
        The map of classes.
        The map of classes; the list of class names associated to each
        value in the segmentation map.
        ..note: This list does NOT include the BACKGROUND class.
        :getter: Gets the base OCT scan
        :setter: Sets the base OCT scan
        :type: :class:`data.OCTscan`
        """
        return self.__classMap

    @classMap.setter
    def classMap(self, cm):  # classMap setter
        if cm is not None:
            # Check that we are receiving the correct type
            if type(cm) is dict:
                self.__classMap = cm
            else:  # Unexpected case. Return warning
                warnMsg = (
                    self.getClassName() + ":classMap: Unexpected type for classMap."
                )
                warnings.warn(warnMsg, SyntaxWarning)
        return None

    # Private methods
    def __str__(self):
        s = (
            "<"
            + self.getClassName()
            + "(["
            + "  scan: "
            + format(self.scan)
            + ","
            + "  data: "
            + format(self.data)
            + ","
            + "  classMap: "
            + format(self.classMap)
            + "])>"
        )
        return s

    # Public methods
    def getClassName(self):
        """Get the class name as a string.
        Get the class name as a string.
        :returns: The class name.
        :rtype: string
        """
        return type(self).__name__

    def clear(self):
        """Clears/Resets the segmentation map to _BACKGROUND.
        Clears/Resets the segmentation map to _BACKGROUND. All pixels are
        assigned the background label.
        """
        self.data = self._BACKGROUND * np.ones(self.scan.shape)
        return None


class OCTscan(object):
    # Sphinx documentation
    """A single OCT scan.
    A single OCT scan. A scan is a grayscale image.
    Check .scantype property for the scan type (A, B or C)
    .. seealso:: :class:`data.OCTvolume`
    .. note:: None
    .. todo:: Upgrade to color scans.
    """

    # Class constructor
    def __init__(self, *args):
        """The class constructor.
        The class constructor.
        tmp = OCTscan() - Creates an black scan sized 480x640.
        tmp = OCTScan(img) - Creates scan from the (grayscale) image. Assumed to be an A scan.
        tmp = OCTScan(img,type) - Creates scan from the (grayscale) image.
        :param img: The scan image
        :type img: numpy.ndarray
        :param type: The scan type ('A' -default-, 'B' or 'C')
        :type type: char
        """
        if len(args) > 1:
            warnMsg = (
                self.getClassName() + ":__init__: Unexpected number of input arguments."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        # Initialize attributes (without decorator @property)
        # Initialize properties (with decorator @property)
        self.data = np.zeros(shape=(480, 640), dtype=np.uint8)  # The default scan image
        self.scantype = "A"
        if len(args) > 0:
            self.data = args[0]
        if len(args) > 1:
            self.scantype = args[1]
        return

    # Properties getters/setters
    #
    # Remember: Sphinx ignores docstrings on property setters so all
    # documentation for a property must be on the @property method
    @property
    def data(self):  # data getter
        """
        The (grayscale) scan image. The image is expected to be
        a grayscale image. Colour images will be converted to grayscale.
        :getter: Gets the OCT scan image
        :setter: Sets the OCT scan image.
        :type: numpy.ndarray shaped [width,height]
        """
        return self.__data

    @data.setter
    def data(self, img):  # data setter
        if img is not None:
            # Check whether the image is in RGB (ndim=3) or in grayscale (ndim=2)
            # and convert to grayscale if necessary
            if img.ndim == 2:
                # Dimensions are only width and height. The image is already in grayscale.
                pass
            elif img.ndim == 3:
                # Image is in RGB. Convert.
                img = color.rgb2gray(img)
            else:  # Unexpected case. Return warning
                warnMsg = self.getClassName() + ":data: Unexpected image shape."
                warnings.warn(warnMsg, SyntaxWarning)
            self.__data = img
        return None

    @property
    def shape(self):  # shape getter
        """
        The scan shape [width,height].
        :getter: Gets the scan shape
        :setter: None. This is a read-only property.
        :type: Tuple [width,height]
        """
        return self.__data.shape

    @shape.setter
    def shape(self, *args):  # shape setter
        # Catching attempts to set the shape of the scan
        warnMsg = self.getClassName() + ":shape: shape is a read-only property."
        warnings.warn(warnMsg, UserWarning)
        return

    @property
    def scantype(self):  # scantype getter
        """
        The scan type; 'A', 'B' or 'C'.
        :getter: Gets the scan type.
        :setter: Sets the scan type.
        :type: char 'A', 'B' or 'C'.
        """
        return self.__scantype

    @scantype.setter
    def scantype(self, *args):  # shape setter
        stype = args[0].upper()  # Uppercase
        if not stype in ("A", "B", "C"):
            warnMsg = (
                self.getClassName() + ":scantype: Scan type can only be "
                "A"
                ", "
                "B"
                " or "
                "C"
                "."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        self.__scantype = stype
        return None

    # Private methods
    def __str__(self):
        s = (
            "<"
            + self.getClassName()
            + "(["
            + "  data: "
            + format(self.data)
            + ","
            + "  shape: "
            + format(self.shape)
            + ","
            + "  scantype: "
            + self.scantype
            + "])>"
        )
        return s

    # Public methods
    def getClassName(self):
        """Get the class name as a string.
        Get the class name as a string.
        :returns: The class name.
        :rtype: string
        """
        return type(self).__name__


class RetinalLayers(object):
    # Sphinx documentation
    """A collection of constants for identifying retinal layers.
    A collection of constants for identifying retinal layers.
    It is basically a dictionary of pairs key:value
    The following retinal layers are considered:
    * Inner Limiting Membrane (ILM)
    * Nerve Fiber Layer (NFL)
    * Ganglion Cell Layer (GCL)
    * Inner Plexiform Layer (IPL)
    * Inner Nuclear Layer (INL)
    * Outer Plexiform Layer (OPL)
    * Outner Nuclear Layer (ONL)
    * External Limiting Membrane (ELM)
    * Rods and Cones layers (RCL)
    * RetinalPigmentEpithelium (RPE)
    Also the non-retinal layer is indexed:
    * Choroid (CHR)
    For pathological conditions, the following is also encoded:
    * Fluid (FLD)
    .. seealso:: None
    .. note:: None
    .. todo:: None
    """

    # Private class attributes shared by all instances
    # Class constructor
    def __init__(self):
        """The class constructor.
        The class constructor.
        tmp = RetinalLayers() - Creates an RetinalLayers object.
        """
        # Initialize attributes
        # Layer constants
        self.layers = {
            "Inner Limiting Membrane": 1,
            "Nerve Fiber Layer": 2,
            "Ganglion Cell Layer": 3,
            "Inner Plexiform Layer": 4,
            "Inner Nuclear Layer": 5,
            "Outer Plexiform Layer": 6,
            "Outer Nuclear Layer": 7,
            "External Limiting Membrane": 8,
            "Rods and Cones Layers": 9,
            "Retinal Pigment Epithelium": 10,
            "Choroid": 20,
            "Fluid": 1001,
        }
        # Currently:
        # RPE includes Bruch's membrane even though most books will
        # 	consider BM as part of the choroid already.
        # ONL includes Henle’s layer
        # ELM is also known as Inner Segment Layer
        # The RCL includes connecting cilia (CL), outer segment layer (OSL) and Verhoeff membrane (VM)
        # ...and so that one can have different names for each layer
        self.layerNames = {
            "ilm": self.layers["Inner Limiting Membrane"],
            "innerlimitingmembrane": self.layers["Inner Limiting Membrane"],
            "nfl": self.layers["Nerve Fiber Layer"],
            "nervefiberlayer": self.layers["Nerve Fiber Layer"],
            "gcl": self.layers["Ganglion Cell Layer"],
            "ganglioncelllayer": self.layers["Ganglion Cell Layer"],
            "ipl": self.layers["Inner Plexiform Layer"],
            "innerplexiformlayer": self.layers["Inner Plexiform Layer"],
            "inl": self.layers["Inner Nuclear Layer"],
            "innernuclearlayer": self.layers["Inner Nuclear Layer"],
            "opl": self.layers["Outer Plexiform Layer"],
            "outerplexiformlayer": self.layers["Outer Plexiform Layer"],
            "onl": self.layers["Outer Nuclear Layer"],
            "outernuclearlayer": self.layers["Outer Nuclear Layer"],
            "elm": self.layers["External Limiting Membrane"],
            "externallimitingmembrane": self.layers["External Limiting Membrane"],
            "rcl": self.layers["Rods and Cones Layers"],
            "rodsandconeslayers": self.layers["Rods and Cones Layers"],
            "rpe": self.layers["Retinal Pigment Epithelium"],
            "retinalpigmentepithelium": self.layers["Retinal Pigment Epithelium"],
            "chr": self.layers["Choroid"],
            "fld": self.layers["Fluid"],
        }
        # Layer acronyms
        self.layerAcronyms = {
            self.layers["Inner Limiting Membrane"]: "ILM",
            self.layers["Nerve Fiber Layer"]: "NFL",
            self.layers["Ganglion Cell Layer"]: "GCL",
            self.layers["Inner Plexiform Layer"]: "IPL",
            self.layers["Inner Nuclear Layer"]: "INL",
            self.layers["Outer Plexiform Layer"]: "OPL",
            self.layers["Outer Nuclear Layer"]: "ONL",
            self.layers["External Limiting Membrane"]: "ELM",
            self.layers["Rods and Cones Layers"]: "RCL",
            self.layers["Retinal Pigment Epithelium"]: "RPE",
            self.layers["Choroid"]: "Choroid",
            self.layers["Fluid"]: "Fluid",
        }
        return

    # Private methods
    # Public methods
    def getClassName(self):
        """Gets the class name
        return: The class name
        rtype: string
        """
        return type(self).__name__

    def getAllLayersIndexes(self):
        """Gets the list of layer values
        return: The list of layers values
        rtype: list
        """
        return list(self.layers.values())

    def getAllLayersNames(self):
        """Gets the list of layer keys
        return: The list of layers keys
        rtype: list
        """
        # Retrieves a list of layer keys
        return list(self.layers.keys())

    def getLayerAcronym(self, idx):
        """Gets the acronym of the i-th layer
        return: The layer acronym e.g. NFL
        rtype: string
        """
        lacronym = "NaN"
        try:
            lacronym = self.layerAcronyms[idx]
        except:
            lacronym = "Unknown"
            print(
                self.getClassName(),
                ':getLayerAcronym: Unexpected layer index. Returning name "',
                lacronym,
                '"',
            )
        return lacronym

    def getLayerIndex(self, layerName):
        """Retrieve the index of a given layer
        return: The index of the layer
        rtype: int
        """
        r = -1
        try:
            layerName = layerName.replace(" ", "")  # Remove whitespaces
            r = self.layerNames[layerName.lower()]  # Ignore case
        except:
            print(
                self.getClassName(),
                ":getLayerIndex: Unknown layer name. Returning index ",
                r,
            )
        return r

    def getLayerName(self, idx):
        """Retrieve the i-th layer name
        return: The name of the i-th layer
        rtype: string
        """
        lname = "Default"
        try:
            # There is no 'direct' method to access the keys given the value.
            lname = list(self.layers.keys())[list(self.layers.values()).index(idx)]
        except:
            lname = "Unknown"
            print(
                self.getClassName(),
                ':getLayerName: Unexpected layer index. Returning name "',
                lname,
                '"',
            )
        return lname

    def getNumLayers(self):
        """Return the number of known layers.
        Return the number of known layers. Please note that this also
        include known non-retinal layers like the choroid.
        return: Length of property map :func:`data.OCTvolume.layers`
        rtype: int
        """
        return len(self.layers)


## Import
import os
import warnings

# from deprecated import deprecated
import deprecation
import math
import numpy as np
from scipy import signal, ndimage  # Used for 2D convolution
from skimage import feature, color
from functools import reduce
import cv2  # That's OpenCV

# import matlab.engine
import matplotlib.pyplot as plt
import _tkinter
from PyQt5.QtCore import Qt  # Imports constants
from PyQt5.QtWidgets import QProgressDialog

# from matplotlib.backends.qt_compat import QtCore, QtWidgets, is_pyqt5
# from version import __version__
import octant

# from data import OCTscan, OCTscanSegmentation, RetinalLayers
# from util import segmentationUtils
# from .Operation import Operation
from scipy.optimize import curve_fit


class OpScanSegment(Operation):
    # Private class attributes shared by all instances
    # Class constructor
    def __init__(self):
        # Call superclass constructor
        super().__init__()
        # Set the operation name
        self.name = "Segmentation"
        # Initialize private attributes unique to this instance
        # self._imgin = np.zeros(shape = (0,0,0), dtype = np.uint8 ); #Input image
        # self._imgout = np.zeros(shape = (0,0,0), dtype = np.uint8 ); #The segmented image

    # Private methods
    @staticmethod
    def anisotropicDiffusionFilter(img):
        """
        Non-linear anisotropic diffusion filtering
        #Original algorithm in [Perona and Malik, (1990) TPAMI 12(7):629-639]
        #Parameters according to [WangRK2005, Proc. of SPIE 5690:380-385]
        :param img: ndarray. The image to be filtered
        :return: the image filtered
        :rtype: ndarray
        """
        sigma = 1
        GaussMask = generateGaussianMask2D(shape=(3, 3), sigma=sigma)
        # print(GaussMask)
        m = 8
        l = 10  # contrast parameter; structures with s > λ are regarded as edges,
        # while with s < λ are assumed to belong to the interior of a region
        # [Wang2005] uses l=5 for the porcine trachea, but for us that value
        # does not produce as good results as l=10.
        timestep = 0.24
        niterations = 10  # diffusion time (iteration)
        Cm = 3.31488
        img2 = img  # Initialize image
        # Get rid of the "large" noise with morphological closing after opening
        morphkernel = np.ones((5, 5), np.uint8)
        # img2 = cv2.morphologyEx(img2, cv2.MORPH_OPEN, morphkernel)
        # img2 = cv2.morphologyEx(img2, cv2.MORPH_CLOSE, morphkernel)
        if np.max(img2) <= 1:
            # Scale parameters
            # m=m/255
            l = l / 255
            # Cm = Cm/255
        for tau in range(1, niterations):
            # 			#Progress bar
            # 			if progress.wasCanceled():
            # 				break
            # 			progress.setValue(round(100*tau/niterations))
            # Morphological removal of noise
            # Note that this is NOT part of the original diffusion filter
            # but the results are clearly enhanced!.
            img2 = cv2.morphologyEx(img2, cv2.MORPH_OPEN, morphkernel)
            img2 = cv2.morphologyEx(img2, cv2.MORPH_CLOSE, morphkernel)
            # Estimate gradient vector at sigma scale
            gradI = np.gradient(img2)
            tmpGradI = np.add.reduce(gradI)  # For illustrations purposes only
            # Regularize gradient
            # gradientScaleSigma=math.sqrt(2*(tau*timestep))
            gradientScaleSigma = math.sqrt(2 * (tau))
            GaussMask = generateGaussianMask2D(shape=(7, 7), sigma=gradientScaleSigma)
            # Individually convolve the Gaussian filter with each gradient component
            # s = signal.convolve2d(gradI, GaussMask, boundary='symm', mode='same')
            s = [None] * len(gradI)  # Preallocate list
            for dim in range(0, len(gradI)):
                s[dim] = signal.convolve2d(
                    gradI[dim], GaussMask, boundary="symm", mode="same"
                )
            s = np.linalg.norm(s, ord=2, axis=0)  # Norm of the gradient.
            # Calculate diffusivity
            tmp = (s / l) ** m
            tmp = np.divide(
                -Cm, tmp, out=np.zeros(tmp.shape), where=tmp != 0
            )  # Avoid division by 0 when estimating diffusivity
            # D = 1-np.exp(-Cm/((s/l)**m)) #diffusivity or conduction coefficient
            D = 1 - np.exp(tmp)  # diffusivity or conduction coefficient
            # Update image
            img2 = img2 + divergence(np.multiply(D, gradI))  # Update the image
            # Reminder> The divergence of gradient is the Laplacian operator
            # See: https://math.stackexchange.com/questions/690493/what-is-divergence-in-image-processing
        return img2

    @staticmethod
    def divergence(f):
        # """Compute the divergence of n-D SCALAR field `f`."""
        # See: https://stackoverflow.com/questions/11435809/compute-divergence-of-vector-field-using-python
        # return reduce(np.add,np.gradient(f))
        """
        Computes the divergence of the VECTOR field f, corresponding to dFx/dx + dFy/dy + ...
        :param f: List of ndarrays, where every item of the list is one dimension of the vector field
        :return: Single ndarray of the same shape as each of the items in f, which corresponds to a scalar field
        """
        num_dims = len(f)
        return np.ufunc.reduce(
            np.add, [np.gradient(f[i], axis=i) for i in range(num_dims)]
        )

    @staticmethod
    def generateGaussianMask2D(shape=(3, 3), sigma=0.5):
        """
        Generates a 2D gaussian mask
        It should give the same result as MATLAB's
        fspecial('gaussian',[shape],[sigma])
        See: https://stackoverflow.com/questions/17190649/how-to-obtain-a-gaussian-filter-in-python
        """
        m, n = [(ss - 1.0) / 2.0 for ss in shape]
        y, x = np.ogrid[-m : m + 1, -n : n + 1]
        h = np.exp(-(x * x + y * y) / (2.0 * sigma * sigma))
        h[h < np.finfo(h.dtype).eps * h.max()] = 0
        sumh = h.sum()
        if sumh != 0:
            h /= sumh
        return h

    # Public methods
    def execute(self, *args, **kwargs):
        """Executes the operation on the :py:attr:`operands`.
        Executes the operation on the :py:attr:`operands` and stores the outcome
        in :py:attr:`result`. Preload operands using method
        :func:`addOperand()`.
        :returns: Result of executing the operation.
        :rtype: :class:`data.OCTscanSegmentation`
        """
        # Ensure the operand has been set.
        if len(self.operands) < 1:
            warnMsg = self.getClassName() + ":execute: Operand not set."
            warnings.warn(warnMsg, SyntaxWarning)
            return None
        # Establish mode of operation
        # NOT WORKING YET... :(
        # Now, it always detect:
        # 	Terminal for the stdin
        # 	GUI for the stdout
        mode = "terminal"
        if os.isatty(0):
            mode = "gui"
        print("Executing segmentation in mode " + mode + ".")
        if mode == "gui":
            MainWindow = QtWidgets.QWidget()
            progress = QProgressDialog("Segmentation...", "Cancel", 0, 100, MainWindow)
            progress.setWindowModality(Qt.WindowModal)
            progress.setAutoReset(True)
            progress.setAutoClose(True)
            progress.setMinimum(0)
            progress.setMaximum(100)
            progress.setWindowTitle("Automatic segmentation")
            progress.setLabelText("Progress:")
            progress.setMinimumDuration(0)
            progress.resize(500, 100)
            progress.forceShow()
            progress.setValue(0)
        imgin = self.operands[0]
        if type(imgin) is OCTscan:
            imgin = imgin.data
        # Define a default output
        segmentedImage = np.zeros(shape=(0, 0, 0), dtype=np.uint8)
        # Check whether the image is in RGB (ndim=3) or in grayscale (ndim=2)
        # and convert to grayscale if necessary
        # img = cv2.cvtColor(self._imgin, cv2.COLOR_BGR2GRAY)
        if imgin.ndim == 2:
            # Dimensions are only width and height. The image is already in grayscale.
            img = imgin
        elif imgin.ndim == 3:
            # Image is in RGB. Convert.
            img = color.rgb2gray(imgin)
        else:  # Unexpected case. Return warning
            print(self.getClassName(), ": Unexpected image shape.")
            return None
        ## Algoritmo 1: Basado en el código de Matlab CASeReL
        # https://pangyuteng.github.io/caserel/
        #
        # Llama externamente a retSegment.exe
        #
        # Requiere de Matlab Runtime environment
        #  #Remove external noise from image
        #  img = segmentationUtils.ejecuta_close(img,4,4) #Clausura
        #  #segmentationUtils.mostrar_imagen(img)
        #
        #  #Generate temporal intermediate image file to be processed externally
        #  tmpFilename = "paso2.tiff"
        #  cv2.imwrite(tmpFilename, img)
        #
        #  #Rely externally on the matlab algorithm for segmentation
        # # eng = matlab.engine.start_matlab()
        # # eng.retSegment(img)
        #  #segmentationUtils.mostrar_imagen(img)
        #
        #  #Delete the intermediate image file
        #  segmentationUtils.elimina_imagen(tmpFilename)
        ## Algoritmo 3 Felipe
        # Step 1) Non-linear anisotropic diffusion filtering
        # Original algorithm in [Perona and Malik, (1990) TPAMI 12(7):629-639]
        # Parameters according to [WangRK2005, Proc. of SPIE 5690:380-385]
        if mode == "gui":
            progress.setLabelText("Applying diffusion filter:")
        img2 = anisotropicDiffusionFilter(img)
        # The output of the anisotropic filter is a float64 2D array
        # with values between [0 1] (although note that IT MAY not
        # include the segment boundaries 0 and/or 1, so for instance
        # while testing I was getting [0,0.85]). Anyway, in this
        # range, anything below 1% (~0.099) is background, and the
        # rest is tissue.
        if mode == "gui":
            progress.setValue(100)
            if progress.wasCanceled():
                return img  # Return the original image
        # print('Diffusion Filter output range: ['+str(np.amin(img2))+','+str(np.amax(img2))+']')
        # Detect background mask
        BGmask = np.where(img2 <= 0.099, True, False)
        # See above for threshold 0.099 on anisotropic diffusion filter
        # output.
        # Note that there might be "background" pixels
        # found somewhere within the tissue. This is VERY LIKELY fluid!
        # Note as well that the lower part also include a bit of the
        # choroid with this threshold.
        # Finally, the top part should be the vitreous.
        # print(BGmask)
        # Just to see the mask
        # img2 = np.where(BGmask,10,OCTscanSegmentation._BACKGROUND)
        # 		#Plot histogram of image
        # 		binSize = 0.01
        # 		bins=np.arange(0,np.amax(img2)+binSize,binSize)
        # 		hist = np.histogram(img2, bins=bins)
        # 		hfig, ax = plt.subplots(1, 1)
        # 		ax.hist(img2.flatten(),bins)
        # 		ax.set_yscale('log')
        # 		hfig.show()
        # PENDING
        # Assign intratissue BG pixels to fluid.
        # Normalize
        img2 = np.floor(255 * (img2 / np.amax(img2))).astype(int)
        img2 = np.where(BGmask == True, OCTscanSegmentation._BACKGROUND, img2)
        # Assign integers to nRetinalLayers
        r = RetinalLayers()
        nRetinalLayers = r.getLayerIndex("RPE")
        # Note that I'm not counting for choroid or fluid.
        # I can't use pixel intensity alone here for segmenting as several
        # (separated) layers exhibit similar reflectivity, so segmentation
        # criteria has to be a mix of location and intensity.
        # Estimate upper and lower boundary by checking first and last
        # non-BG pixels per column
        # UNFINISHED.
        # upperBoundary =
        # [uVals,uIdx,uInv,uCounts] = np.unique(img2, \
        # 						return_index=True, \
        # 						return_inverse=True, \
        # 						return_counts=True)
        # print(len(uVals))
        # 		## Algoritmo 2: Arlem
        # 		img2 = img
        # 		print(np.amax(img2))
        # 		tmp=np.reshape(img2,(np.size(img2),1))
        # 		print(type(tmp))
        # 		print(np.size(tmp))
        # 		print(np.amax(tmp))
        #
        #
        # 		hfig, ax = plt.subplots(1, 1)
        # 		ax.hist(tmp,255)
        # 		hfig.show()
        #
        # 		#Elimina ruido
        # 		img = segmentationUtils.ejecuta_close(img,4,4)
        #
        # 		#Amplifica capas
        # 		img = segmentationUtils.ejecuta_dilate(img,5,20,1)
        #
        # 		#Tensor
        # 		Axx, Axy, Ayy = feature.structure_tensor(img)
        #
        # 		#Elimina mas ruido
        # 		Ayy = segmentationUtils.ejecuta_close(Ayy,6,1)
        #
        # 		#Resalta las capas que sean mayores a la media
        # 		Ayy = segmentationUtils.resalta_bordes(Ayy,True,0)
        #
        # 		#Elimina aun mas ruido
        # 		Ayy = segmentationUtils.ejecuta_open(Ayy,1,1)
        #
        # 		#Binarizacion
        # 		binary = segmentationUtils.ejecuta_OTSU(Ayy)
        #
        # 		#elimina ruido del posible borde superior
        # 		binary = segmentationUtils.ejecuta_elimina_ruido_extremos(True,0,0,binary)
        #
        # 		#elimina ruido del posible borde inferior
        # 		binary = segmentationUtils.ejecuta_elimina_ruido_extremos(False,0,0,binary)
        #
        # 		#obtiene bordes exteriores
        # 		arraySuperior, arrayInferior = segmentationUtils.obten_bordes_externos(binary)
        #
        # 		#elimina ruido a la imagen original
        # 		img2 = segmentationUtils.elimina_desde_arreglos(img2, arraySuperior, arrayInferior)
        # 		img2 = segmentationUtils.ejecuta_close(img2,2,1)
        # 		img2 = feature.canny(img2,sigma = 2.5)
        # 		img2 = segmentationUtils.elimina_ruido_canny(img2,1)
        # Hide and close progress bar.
        if mode == "gui":
            progress.hide()
        # if isinstance(imgin,(OCTscan,)):
        if type(imgin) is OCTscan:
            self.result = OCTscanSegmentation(imgin)
        else:
            self.result = OCTscanSegmentation(OCTscan(imgin))
        self.result.data = img2
        return self.result


# 	#@deprecated(version='0.2', reason="Deprecated. Use method execute() instead.")
# 	@deprecation.deprecated(deprecated_in="0.2", removed_in="0.3",
# 						current_version=__version__,
# 						details="Use method execute() instead.")
# 	def segmentar(self,image):
# 		#Encapsulate the image as an OCTscan
# 		tmp=OCTscan(image)
# 		self.clear()
# 		self.addOperand(tmp)
# 		#Execute
# 		self.execute()
# 		return None
def anisotropicDiffusionFilter(img):
    """
    Non-linear anisotropic diffusion filtering
    #Original algorithm in [Perona and Malik, (1990) TPAMI 12(7):629-639]
    #Parameters according to [WangRK2005, Proc. of SPIE 5690:380-385]
    :param img: ndarray. The image to be filtered
    :return: the image filtered
    :rtype: ndarray
    """
    sigma = 1
    GaussMask = generateGaussianMask2D(shape=(3, 3), sigma=sigma)
    # print(GaussMask)
    m = 8
    l = 10  # contrast parameter; structures with s > λ are regarded as edges,
    # while with s < λ are assumed to belong to the interior of a region
    # [Wang2005] uses l=5 for the porcine trachea, but for us that value
    # does not produce as good results as l=10.
    timestep = 0.24
    niterations = 10  # diffusion time (iteration)
    Cm = 3.31488
    img2 = img  # Initialize image
    # Get rid of the "large" noise with morphological closing after opening
    morphkernel = np.ones((5, 5), np.uint8)
    # img2 = cv2.morphologyEx(img2, cv2.MORPH_OPEN, morphkernel)
    # img2 = cv2.morphologyEx(img2, cv2.MORPH_CLOSE, morphkernel)
    if np.max(img2) <= 1:
        # Scale parameters
        # m=m/255
        l = l / 255
        # Cm = Cm/255
    for tau in range(1, niterations):
        # 			#Progress bar
        # 			if progress.wasCanceled():
        # 				break
        # 			progress.setValue(round(100*tau/niterations))
        # Morphological removal of noise
        # Note that this is NOT part of the original diffusion filter
        # but the results are clearly enhanced!.
        img2 = cv2.morphologyEx(img2, cv2.MORPH_OPEN, morphkernel)
        img2 = cv2.morphologyEx(img2, cv2.MORPH_CLOSE, morphkernel)
        # Estimate gradient vector at sigma scale
        gradI = np.gradient(img2)
        tmpGradI = np.add.reduce(gradI)  # For illustrations purposes only
        # Regularize gradient
        # gradientScaleSigma=math.sqrt(2*(tau*timestep))
        gradientScaleSigma = math.sqrt(2 * (tau))
        GaussMask = generateGaussianMask2D(shape=(7, 7), sigma=gradientScaleSigma)
        # Individually convolve the Gaussian filter with each gradient component
        # s = signal.convolve2d(gradI, GaussMask, boundary='symm', mode='same')
        s = [None] * len(gradI)  # Preallocate list
        for dim in range(0, len(gradI)):
            s[dim] = signal.convolve2d(
                gradI[dim], GaussMask, boundary="symm", mode="same"
            )
        s = np.linalg.norm(s, ord=2, axis=0)  # Norm of the gradient.
        # Calculate diffusivity
        tmp = (s / l) ** m
        tmp = np.divide(
            -Cm, tmp, out=np.zeros(tmp.shape), where=tmp != 0
        )  # Avoid division by 0 when estimating diffusivity
        # D = 1-np.exp(-Cm/((s/l)**m)) #diffusivity or conduction coefficient
        D = 1 - np.exp(tmp)  # diffusivity or conduction coefficient
        # Update image
        img2 = img2 + divergence(np.multiply(D, gradI))  # Update the image
        # Reminder> The divergence of gradient is the Laplacian operator
        # See: https://math.stackexchange.com/questions/690493/what-is-divergence-in-image-processing
    return img2


segmentedImage = np.zeros(shape=(0, 0, 0), dtype=np.uint8)
img2 = anisotropicDiffusionFilter(img)
BGmask = np.where(img2 <= 0.099, True, False)
img2 = np.floor(255 * (img2 / np.amax(img2))).astype(int)
img2 = np.where(BGmask == True, OCTscanSegmentation._BACKGROUND, img2)
# Assign integers to nRetinalLayers
r = RetinalLayers()
nRetinalLayers = r.getLayerIndex("RPE")
OCTscanSegmentation(imgin)


def divergence(f):
    # """Compute the divergence of n-D SCALAR field `f`."""
    # See: https://stackoverflow.com/questions/11435809/compute-divergence-of-vector-field-using-python
    # return reduce(np.add,np.gradient(f))
    """
    Computes the divergence of the VECTOR field f, corresponding to dFx/dx + dFy/dy + ...
    :param f: List of ndarrays, where every item of the list is one dimension of the vector field
    :return: Single ndarray of the same shape as each of the items in f, which corresponds to a scalar field
    """
    num_dims = len(f)
    return np.ufunc.reduce(np.add, [np.gradient(f[i], axis=i) for i in range(num_dims)])


def generateGaussianMask2D(shape=(3, 3), sigma=0.5):
    """
    Generates a 2D gaussian mask
    It should give the same result as MATLAB's
    fspecial('gaussian',[shape],[sigma])
    See: https://stackoverflow.com/questions/17190649/how-to-obtain-a-gaussian-filter-in-python
    """
    m, n = [(ss - 1.0) / 2.0 for ss in shape]
    y, x = np.ogrid[-m : m + 1, -n : n + 1]
    h = np.exp(-(x * x + y * y) / (2.0 * sigma * sigma))
    h[h < np.finfo(h.dtype).eps * h.max()] = 0
    sumh = h.sum()
    if sumh != 0:
        h /= sumh
    return h


generateGaussianMask2D()


def execute(self, *args, **kwargs):
    """Executes the operation on the :py:attr:`operands`.
    Executes the operation on the :py:attr:`operands` and stores the outcome
    in :py:attr:`result`. Preload operands using method
    :func:`addOperand()`.
    :returns: Result of executing the operation.
    :rtype: :class:`data.OCTscanSegmentation`
    """
    # Ensure the operand has been set.
    if len(self.operands) < 1:
        warnMsg = self.getClassName() + ":execute: Operand not set."
        warnings.warn(warnMsg, SyntaxWarning)
        return None
    # Establish mode of operation
    # NOT WORKING YET... :(
    # Now, it always detect:
    # 	Terminal for the stdin
    # 	GUI for the stdout
    mode = "terminal"
    if os.isatty(0):
        mode = "gui"
    print("Executing segmentation in mode " + mode + ".")
    if mode == "gui":
        MainWindow = QtWidgets.QWidget()
        progress = QProgressDialog("Segmentation...", "Cancel", 0, 100, MainWindow)
        progress.setWindowModality(Qt.WindowModal)
        progress.setAutoReset(True)
        progress.setAutoClose(True)
        progress.setMinimum(0)
        progress.setMaximum(100)
        progress.setWindowTitle("Automatic segmentation")
        progress.setLabelText("Progress:")
        progress.setMinimumDuration(0)
        progress.resize(500, 100)
        progress.forceShow()
        progress.setValue(0)
    imgin = self.operands[0]
    if type(imgin) is OCTscan:
        imgin = imgin.data
    # Define a default output
    segmentedImage = np.zeros(shape=(0, 0, 0), dtype=np.uint8)
    # Check whether the image is in RGB (ndim=3) or in grayscale (ndim=2)
    # and convert to grayscale if necessary
    # img = cv2.cvtColor(self._imgin, cv2.COLOR_BGR2GRAY)
    if imgin.ndim == 2:
        # Dimensions are only width and height. The image is already in grayscale.
        img = imgin
    elif imgin.ndim == 3:
        # Image is in RGB. Convert.
        img = color.rgb2gray(imgin)
    else:  # Unexpected case. Return warning
        print(self.getClassName(), ": Unexpected image shape.")
        return None
    ## Algoritmo 1: Basado en el código de Matlab CASeReL
    # https://pangyuteng.github.io/caserel/
    #
    # Llama externamente a retSegment.exe
    #
    # Requiere de Matlab Runtime environment
    #  #Remove external noise from image
    #  img = segmentationUtils.ejecuta_close(img,4,4) #Clausura
    #  #segmentationUtils.mostrar_imagen(img)
    #
    #  #Generate temporal intermediate image file to be processed externally
    #  tmpFilename = "paso2.tiff"
    #  cv2.imwrite(tmpFilename, img)
    #
    #  #Rely externally on the matlab algorithm for segmentation
    # # eng = matlab.engine.start_matlab()
    # # eng.retSegment(img)
    #  #segmentationUtils.mostrar_imagen(img)
    #
    #  #Delete the intermediate image file
    #  segmentationUtils.elimina_imagen(tmpFilename)
    ## Algoritmo 3 Felipe
    # Step 1) Non-linear anisotropic diffusion filtering
    # Original algorithm in [Perona and Malik, (1990) TPAMI 12(7):629-639]
    # Parameters according to [WangRK2005, Proc. of SPIE 5690:380-385]
    if mode == "gui":
        progress.setLabelText("Applying diffusion filter:")
    img2 = anisotropicDiffusionFilter(img)
    # The output of the anisotropic filter is a float64 2D array
    # with values between [0 1] (although note that IT MAY not
    # include the segment boundaries 0 and/or 1, so for instance
    # while testing I was getting [0,0.85]). Anyway, in this
    # range, anything below 1% (~0.099) is background, and the
    # rest is tissue.
    if mode == "gui":
        progress.setValue(100)
        if progress.wasCanceled():
            return img  # Return the original image
    # print('Diffusion Filter output range: ['+str(np.amin(img2))+','+str(np.amax(img2))+']')
    # Detect background mask
    BGmask = np.where(img2 <= 0.099, True, False)
    # See above for threshold 0.099 on anisotropic diffusion filter
    # output.
    # Note that there might be "background" pixels
    # found somewhere within the tissue. This is VERY LIKELY fluid!
    # Note as well that the lower part also include a bit of the
    # choroid with this threshold.
    # Finally, the top part should be the vitreous.
    # print(BGmask)
    # Just to see the mask
    # img2 = np.where(BGmask,10,OCTscanSegmentation._BACKGROUND)
    # 		#Plot histogram of image
    # 		binSize = 0.01
    # 		bins=np.arange(0,np.amax(img2)+binSize,binSize)
    # 		hist = np.histogram(img2, bins=bins)
    # 		hfig, ax = plt.subplots(1, 1)
    # 		ax.hist(img2.flatten(),bins)
    # 		ax.set_yscale('log')
    # 		hfig.show()
    # PENDING
    # Assign intratissue BG pixels to fluid.
    # Normalize
    img2 = np.floor(255 * (img2 / np.amax(img2))).astype(int)
    img2 = np.where(BGmask == True, OCTscanSegmentation._BACKGROUND, img2)
    # Assign integers to nRetinalLayers
    r = RetinalLayers()
    nRetinalLayers = r.getLayerIndex("RPE")
    # Note that I'm not counting for choroid or fluid.
    # I can't use pixel intensity alone here for segmenting as several
    # (separated) layers exhibit similar reflectivity, so segmentation
    # criteria has to be a mix of location and intensity.
    # Estimate upper and lower boundary by checking first and last
    # non-BG pixels per column
    # UNFINISHED.
    # upperBoundary =
    # [uVals,uIdx,uInv,uCounts] = np.unique(img2, \
    # 						return_index=True, \
    # 						return_inverse=True, \
    # 						return_counts=True)
    # print(len(uVals))
    # 		## Algoritmo 2: Arlem
    # 		img2 = img
    # 		print(np.amax(img2))
    # 		tmp=np.reshape(img2,(np.size(img2),1))
    # 		print(type(tmp))
    # 		print(np.size(tmp))
    # 		print(np.amax(tmp))
    #
    #
    # 		hfig, ax = plt.subplots(1, 1)
    # 		ax.hist(tmp,255)
    # 		hfig.show()
    #
    # 		#Elimina ruido
    # 		img = segmentationUtils.ejecuta_close(img,4,4)
    #
    # 		#Amplifica capas
    # 		img = segmentationUtils.ejecuta_dilate(img,5,20,1)
    #
    # 		#Tensor
    # 		Axx, Axy, Ayy = feature.structure_tensor(img)
    #
    # 		#Elimina mas ruido
    # 		Ayy = segmentationUtils.ejecuta_close(Ayy,6,1)
    #
    # 		#Resalta las capas que sean mayores a la media
    # 		Ayy = segmentationUtils.resalta_bordes(Ayy,True,0)
    #
    # 		#Elimina aun mas ruido
    # 		Ayy = segmentationUtils.ejecuta_open(Ayy,1,1)
    #
    # 		#Binarizacion
    # 		binary = segmentationUtils.ejecuta_OTSU(Ayy)
    #
    # 		#elimina ruido del posible borde superior
    # 		binary = segmentationUtils.ejecuta_elimina_ruido_extremos(True,0,0,binary)
    #
    # 		#elimina ruido del posible borde inferior
    # 		binary = segmentationUtils.ejecuta_elimina_ruido_extremos(False,0,0,binary)
    #
    # 		#obtiene bordes exteriores
    # 		arraySuperior, arrayInferior = segmentationUtils.obten_bordes_externos(binary)
    #
    # 		#elimina ruido a la imagen original
    # 		img2 = segmentationUtils.elimina_desde_arreglos(img2, arraySuperior, arrayInferior)
    # 		img2 = segmentationUtils.ejecuta_close(img2,2,1)
    # 		img2 = feature.canny(img2,sigma = 2.5)
    # 		img2 = segmentationUtils.elimina_ruido_canny(img2,1)
    # Hide and close progress bar.
    if mode == "gui":
        progress.hide()
    # if isinstance(imgin,(OCTscan,)):
    if type(imgin) is OCTscan:
        self.result = OCTscanSegmentation(imgin)
    else:
        self.result = OCTscanSegmentation(OCTscan(imgin))
    self.result.data = img2
    return self.result


class OCTvolume(object):
    # Sphinx documentation
    """A set of :class:`data.OCTscan`
    A set of :class:`data.OCTscan`
    .. seealso:: :class:`data.OCTscan`
    .. note:: None
    .. todo:: None
    """

    # Class constructor
    def __init__(self):
        """The class constructor.
        The class constructor.
        tmp = OCTvolume() - Creates an empty volume with no scans.
        """
        self.scans = list()
        return

    # Properties getters/setters
    @property
    def scans(self):  # scans getter
        """
        The set of OCT scans.
        :getter: Gets the set of scans
        :setter: Sets the set of scans
        :type: list. All scans are of :class:`data.OCTscan`
        .. seealso:: :func:`data.OCTvolume.addScan` , :func:`data.OCTvolume.removeScan` , :func:`data.OCTvolume.clear`
        """
        return self.__scans

    @scans.setter
    def scans(self, *args):  # scans setter
        tmpScanSet = args[0]
        if type(tmpScanSet) is not list:
            warnMsg = (
                self.getClassName() + ":scans: Unexpected type. "
                "Please provide a list of data.OCTscan."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        else:
            self.__scans = list()
            for x in tmpScanSet:
                if type(x) is not OCTscan:
                    warnMsg = (
                        self.getClassName() + ":scans: Unexpected scan type "
                        "for object " + x + ". Skipping object."
                    )
                    warnings.warn(warnMsg, SyntaxWarning)
                else:
                    self.__scans.append(x)
        return None

    # Public methods
    def getClassName(self):
        return type(self).__name__

    def clear(self):
        """
        Clears the OCT volume; Removes all scans.
        :return: None
        """
        self.__scans = list()
        return None

    @deprecation.deprecated(
        deprecated_in="0.3",
        removed_in="1.0",
        #                         current_version=__version__,
        details="Use method addScans() instead.",
    )
    def addScan(self, theScan):
        """
        Add an OCT scan to the volume.
        :param theScan: The OCT scan.
        :type theScan: :class:`data.OCTscan`
        :return: None
        """
        if type(theScan) is not OCTscan:
            warnMsg = (
                self.getClassName() + ":addScan: Unexpected scan type. "
                "Nothing will be added."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        else:
            self.__scans.append(theScan)
        return None

    def addScans(self, theScans):
        """
        Add one or more OCT scans to the volume at once.
        :param theScans: The list of OCT scans.
        :type theScan: list of :class:`data.OCTscan` or
            single :class:`data.OCTscan`
        :return: True if scans added, False otherwise.
        """
        flagAllOCTscans = True
        if type(theScans) is OCTscan:
            self.__scans.append(theScans)
        elif type(theScans) is list:
            for elem in theScans:
                if type(elem) is not OCTscan:
                    warnMsg = (
                        self.getClassName() + ":addScans: Unexpected scan type. "
                        "Nothing will be added."
                    )
                    flagAllOCTscans = False
                    break
            if flagAllOCTscans:
                self.__scans.extend(theScans)
        return flagAllOCTscans

    def getScans(self, t):
        """
        Retrieves all scans in the volume of type t.
        :param t: Scan type 'A', 'B' or 'C'
        :type t: char
        :return: The set of A scans in the volume
        :rtype: list
        .. seealso:: :func:`data.OCTvolume.getVolume`
        """
        t = t.upper()  # Uppercase
        theScans = list()
        for x in self.__scans:
            if x.scantype == t:
                theScans.append(x)
        return theScans

    def getVolume(self, t):
        """
        Retrieves the (sub-)volume of scans of type t.
        :param t: Scan type 'A', 'B' or 'C'
        :type t: char
        :return: A volume with the set of scans of type t.
        :rtype: :class:`data.OCTvolume`
        .. seealso:: :func:`data.OCTvolume.getScans`
        """
        t = t.upper()  # Uppercase
        theScans = OCTvolume()
        for x in self.__scans:
            if x.scantype == t:
                theScans.addScan(x)
        return theScans

    def getNScans(self, t=None):
        """Get the number of scans of a certain type.
        If type is not given, then the total number of scans is given.
        :param t: Scan type 'A', 'B' or 'C' or None
        :type t: char
        :return: The number of scans of a certain type
        :rtype: int
        """
        if t is None:
            res = len(self.__scans)
        else:
            res = self.getNScans(self.getVolume(t))
        return res


class OCTvolumeSegmentation(object):
    # Sphinx documentation
    """A set of :class:`data.OCTscanSegmentation`
    A set of :class:`data.OCTscanSegmentation`
    .. seealso:: :class:`data.OCTscanSegmentation`
    .. note:: None
    .. todo:: None
    """

    # Class constructor
    def __init__(self):
        """The class constructor.
        The class constructor.
        tmp = OCTvolumeSegmentation() - Creates an empty segmentation volume with no scans.
        """
        self.scanSegmentations = list()
        return

    # Properties getters/setters
    @property
    def scanSegmentations(self):  # scans getter
        """
        The set of OCT scans segmentations.
        :getter: Gets the set of scans segmentations
        :setter: Sets the set of scans segmentations
        :type: list. All scans are of :class:`data.OCTscanSegmentation`
        .. seealso:: :func:`data.OCTvolumeSegmentation.addScan` , :func:`data.OCTvolumeSegmentation.removeScan` , :func:`data.OCTvolumeSegmentation.clear`
        """
        return self.__scanSegmentations

    @scanSegmentations.setter
    def scanSegmentations(self, *args):  # scanSegmentations setter
        tmpScanSet = args[0]
        if type(tmpScanSet) is not list:
            warnMsg = (
                self.getClassName() + ":scanSegmentations: Unexpected type. "
                "Please provide a list of data.OCTscanSegmentation."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        else:
            self.__scanSegmentations = list()
            for x in tmpScanSet:
                if type(x) is not OCTscanSegmentation:
                    warnMsg = (
                        self.getClassName()
                        + ":scanSegmentations: Unexpected scan type "
                        "for object " + x + ". Skipping object."
                    )
                    warnings.warn(warnMsg, SyntaxWarning)
                else:
                    self.__scanSegmentations.append(x)
        return None

    # Public methods
    def clear(self):
        """
        Clears the OCT segmentation volume; Removes all scan segmentations.
        :return: None
        """
        self.__scanSegmentations = list()
        return None

    def addScanSegmentations(self, theScanSegmentations):
        """
        Add one or multiple OCT scan segmentation to the volume.
        :param theScanSegmentations: The list of OCT scan segmentations
        :type theScanSegmentations: list of :class:`data.OCTscanSegmentation` or
            single :class:`data.OCTscanSegmentation`
        :return: None
        """
        flagAllOCTscans = False
        if type(theScanSegmentations) is OCTscanSegmentation:
            flagAllOCTscans = True
            self.__scanSegmentations.append(theScanSegmentations)
        elif type(theScanSegmentations) is list:
            flagAllOCTscans = True
            for elem in theScans:
                if type(elem) is not OCTscanSegmentation:
                    warnMsg = (
                        self.getClassName() + ":addScans: Unexpected scan type. "
                        "Nothing will be added."
                    )
                    flagAllOCTscans = False
                    break
            if flagAllOCTscans:
                self.__scanSegmentations.extend(theScanSegmentations)
        return flagAllOCTscans

    def getScanSegmentations(self, t):
        """
        Retrieves all scans in the volume of type t.
        :param t: Scan type 'A', 'B' or 'C' or scan indexing
        :type t: str, list or int
        :return: The set of scans in the volume of the chosen
        :rtype: list
        .. seealso:: :func:`data.OCTvolumeSegmentation.getVolume`
        """
        theScans = list()
        if type(t) is str:
            t = t.upper()  # Uppercase
            for x in self.__scanSegmentations:
                if x.scantype == t:
                    theScans.append(x)
        elif type(t) is list:
            for elem in t:
                theScans.append(self.__scanSegmentations[elem])
        elif type(t) is int:
            theScans.append(self.__scanSegmentations[elem])
        return theScans

    def getVolume(self, t):
        """
        Retrieves the (sub-)volume of scans of type t.
        :param t: Scan type 'A', 'B' or 'C'
        :type t: char
        :return: A volume with the set of scans of type t.
        :rtype: :class:`data.OCTvolume`
        .. seealso:: :func:`data.OCTvolume.getScans`
        """
        t = t.upper()  # Uppercase
        theScans = OCTvolumeSegmentation()
        for x in self.__scanSegmentations:
            if x.scantype == t:
                theScans.addScanSegmentation(x)
        return theScans

    def getNScans(self):
        """Get the number of scans segmentations.
        :return: The number of scans segmentations
        :rtype: int
        """
        return len(self.__scanSegmentations)


class Document:
    # Sphinx documentation
    """The document class for
    The document class for
    This class represents a document in  A document holds information
    about a study plus some additional metadata.
    Currently, a study is only an OCT image (with several scans) with or without
    segmentation information.
    .. seealso:: None
    .. note:: None
    .. todo:: None
    """

    # Private class attributes shared by all instances
    # Class constructor
    def __init__(self):
        # Call superclass constructor
        # Initialize private attributes unique to this instance
        self.study = OCTvolume()  # The current study.
        # Currently, an OCT volume
        self.segmentation = OCTvolumeSegmentation()  # The current study.
        # Currently, an OCT volumeSegmentation
        self.docsettings = Settings()
        self.docsettings.selectedScan = None
        # shared between study and segmentation
        # Document metadata
        self.folderName = "."  # Folder where the document is currently stored
        self.fileName = "OCTantDocument0001"  # The filename of the document
        self.name = "OCTantDocument0001"  # A logical name for the study
        return

    # Properties getters/setters
    #
    # Remember: Sphinx ignores docstrings on property setters so all
    # documentation for a property must be on the @property method
    @property
    def docsettings(self):  # docsettings getter
        """
        The application settings.
        :getter: Gets the document settings
        :setter: Sets the document settings
        :type: class:`data.Settings`
        """
        return self.__appsettings

    @docsettings.setter
    def docsettings(self, newSettings):  # document setter
        if newSettings is None:
            newSettings = Settings()  # Initialize settings
        if type(newSettings) is Settings:
            self.__appsettings = newSettings
        else:
            warnMsg = self.getClassName() + ":docsettings: Unexpected settings type."
            warnings.warn(warnMsg, SyntaxWarning)
        return None

    @property
    def study(self):  # study getter
        """
        The OCT volume being processed and analysed.
        ..todo: Upgrade to volume. Watch out! This will affect many
            other classess using this method.
        :getter: Gets the OCT volume.
        :setter: Sets the OCT volume.
        :type: :class:`data.OCTvolume`
        """
        return self.__study

    @study.setter
    def study(self, vol):  # study setter
        if vol is None or type(vol) is OCTvolume:
            self.__study = vol
            # ...and reset scan
            self.segmentedScan = None
        elif type(vol) is OCTscan:
            warnMsg = (
                self.getClassName()
                + ":study: OCTvolume expected but OCTscan received. Embedding scan."
            )
            warnings.warn(warnMsg, SyntaxWarning)
            self.__study = OCTvolume()
            self.__study.addScan(vol)
        else:
            warnMsg = self.getClassName() + ":study: Unexpected study type."
            warnings.warn(warnMsg, SyntaxWarning)
        return None

    @property
    def segmentation(self):  # segmentation getter
        """
        The segmentation over the OCT study being processed and analysed.
        :getter: Gets the OCT volume segmentation.
        :setter: Sets the OCT volume segmentation.
        :type: :class:`data.OCTvolumeSegmentation`
        """
        return self.__segmentation

    @segmentation.setter
    def segmentation(self, newSegmentation):  # segmentation setter
        if (newSegmentation is None) or (
            type(newSegmentation) is OCTvolumeSegmentation
        ):
            self.__segmentation = newSegmentation
            if newSegmentation is not None:
                if self.study is None:
                    warnMsg = self.getClassName() + ":segmentation: No reference image."
                    warnings.warn(warnMsg, SyntaxWarning)
                if not (
                    len(newSegmentation.scanSegmentations) == len(self.study.scans)
                ):
                    warnMsg = self.getClassName() + ":segmentation: Unexpected size."
                    warnings.warn(warnMsg, SyntaxWarning)
        elif type(newSegmentation) is OCTscanSegmentation:
            warnMsg = (
                self.getClassName()
                + ":study: OCTvolumeSegmentation expected but OCTscanSegmentation received. Embedding scan."
            )
            warnings.warn(warnMsg, SyntaxWarning)
            self.__segmentation = OCTvolumeSegmentation()
            self.__segmentation.addScanSegmentations(newSegmentation)
        else:
            warnMsg = (
                self.getClassName() + ":segmentation: Unexpected segmented scan type."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        return None

    @property
    def folderName(self):  # folderName getter
        """
        Folder where the document is currently stored.
        ..note: Also retrieve the py:attr:`fileName` to build the full path.
        :getter: Gets the study folder name.
        :setter: Sets the study folder name. If new folder is None,
            the current directory '.' is chosen.
        :type: str
        """
        return self.__folderName

    @folderName.setter
    def folderName(self, d):  # name setter
        if d is None:
            d = "."  # Set to current folder
        if type(d) is str:
            self.__folderName = d
        else:
            warnMsg = self.getClassName() + ":name: Unexpected folderName type."
            warnings.warn(warnMsg, SyntaxWarning)
        return None

    @property
    def fileName(self):  # fileName getter
        """
        The filename of the document.
        ..note: Also retrieve the py:attr:`folderName` to build the full path.
        :getter: Gets the the filename of the document.
        :setter: Sets the The filename of the document. If new name is None,
            a default name is given.
        :type: str
        """
        return self.__folderName

    @fileName.setter
    def fileName(self, newFilename):  # fileName setter
        if newFilename is None:
            newFilename = self.defaultName  # Set to default name
        if type(newFilename) is str:
            self.__fileName = newFilename
        else:
            warnMsg = self.getClassName() + ":name: Unexpected fileName type."
            warnings.warn(warnMsg, SyntaxWarning)
        return None

    @property
    def name(self):  # name getter
        """
        A logical name for the study.
        :getter: Gets the OCT study name.
        :setter: Sets the OCT study name.
        :type: str
        """
        return self.__name

    @name.setter
    def name(self, newName):  # name setter
        if newName is None or type(newName) is str:
            self.__name = newName
        else:
            warnMsg = self.getClassName() + ":name: Unexpected name type."
            warnings.warn(warnMsg, SyntaxWarning)
        return None

    # Private methods
    # Public methods
    def getClassName(self):
        return type(self).__name__

    def getCurrentScan(self):
        """Get the current working OCT scan
        Change the current selection using :func:`pickScan`
        :returns: The current working OCT scan.
        :rtype: :class:`data.OCTscan` or None if the study contains no scans
        """
        if self.docsettings.selectedScan is None:
            self.docsettings.selectedScan = 0
        res = None
        if len(self.__study.scans) > 0:
            res = self.__study.scans[self.docsettings.selectedScan]
        return res

    def setCurrentScan(self, newScan):
        """Sets the current working OCT scan
        Change the current selection using :func:`pickScan`
        :param newScan: An OCT scan to be assigned to the current working OCT scan.
        :type newScan: :class:`data.OCTscan`
        """
        if self.docsettings.selectedScan is None:
            self.docsettings.selectedScan = 0
        if newScan is None:
            if self.__study.getNScans() == 0:
                # do nothing.
                pass
            else:
                warnMsg = (
                    self.getClassName()
                    + ":setcurrentscan: Unexpected scan type NoneType."
                )
                warnings.warn(warnMsg, SyntaxWarning)
        if type(newScan) is OCTscan:
            self.__study.scans[self.docsettings.selectedScan] = newScan
        else:
            warnMsg = self.getClassName() + ":setcurrentscan: Unexpected scan type."
            warnings.warn(warnMsg, SyntaxWarning)
        return None

    def pickScan(self, i):
        """Pick the i-th OCT scan (and its segmentation) for working.
        Sets the docsetting.selectedScan to i checking that it does exist.
        :param i: The selected scan index
        :type i: int
        :return: None
        """
        if type(i) is int and i < self.study.getNScans() and i >= 0:
            self.docsettings.selectedScan = i
        else:
            warnMsg = self.getClassName() + ":pickScan: Selected scan does not exist."
            warnings.warn(warnMsg, SyntaxWarning)
        return None

    def getCurrentScanSegmentation(self):
        """Get the current working OCT scanSegmentation
        Change the current selection using :func:`pickScan`
        :returns: The current working OCT scanSegmentation.
        :rtype: :class:`data.OCTscanSegmentation`
        """
        if self.docsettings.selectedScan is None:
            self.docsettings.selectedScan = 0
        res = None
        if len(self.__segmentation.scanSegmentations) > 0:
            res = self.__segmentation.scanSegmentations[self.docsettings.selectedScan]
        return res

    def setCurrentScanSegmentation(self, newScan):
        """Sets the current working OCT scanSegmentation
        Change the current selection using :func:`pickScan`
        :param newScan: An OCT scan to be assigned to the current working OCT scan.
        :type newScan: :class:`data.OCTscanSegmentation`
        """
        if self.docsettings.selectedScan is None:
            self.docsettings.selectedScan = 0
        if newScan is None:
            if self.__segmentation.getNScans() == 0:
                # do nothing.
                pass
            else:
                warnMsg = (
                    self.getClassName()
                    + ":setcurrentscansegmentation: Unexpected scan type NoneType."
                )
                warnings.warn(warnMsg, SyntaxWarning)
        if type(newScan) is OCTscanSegmentation:
            self.__segmentation.scanSegmentations[
                self.docsettings.selectedScan
            ] = newScan
        else:
            warnMsg = (
                self.getClassName()
                + ":setcurrentscansegmentation: Unexpected scan type."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        return None

    def readFile(self, filename):
        """Reads an OCTant document file.
        This method is currently a sham, and it will be updated
        when serialization is incorporated to  Currently,
        it returns an empty document. Nevertheless, it already
        updates the document, clearing all fields to default values,
        and updates the filename and folder
        The file must exist or an error is generated.
        The file must be in OCTant file format.
        :param fileName: The file name
        :type fileName: str
        :return: This document
        :rtype: :class:`data.Document`
        """
        self = Document()
        self.folderName, self.fileName = os.path.split(filename)
        return self


class Settings(object):
    # Sphinx documentation
    """A class to hold a list of settings.
    A class to hold a list of settings. A list of settings is the classical
    "dictionary" (see note below on python's dict) of pairs key:value
    but with some additional capabilities. In particular, the class provides
    additional file reading and writing capabilities so that settings can
    be read to and from plain text files, as well as some value setting
    checking capabilities.
    The class is intended to behave like a dynamic struct
    where properties of the class, i.e. new settings, can be declared
    "on-the-fly" instead of being predefined.
    Although, creating a dynamic struct class in python itself is trivial
    (see https://stackoverflow.com/questions/1878710/struct-objects-in-python ),
    but because of the additional capabilities, hence the convenience of
    the class.
    .. Background:
    MATLAB's struct allows this "on-the-fly" field declaration on the fly.
    Python's built-in dictionary is not exactly a match because of the
    required syntax i.e. ``mySettingsObj['fieldname']`` instead of
    ``mySettingsObj.fieldname`` and the inability to control value settings.
    .. seealso:: None
    .. note:: None
    .. todo:: None
    """

    # Private class attributes shared by all instances
    # Class constructor
    def __init__(self, **kwargs):
        """The class constructor."""
        # Call superclass constructor
        # Initialize private attributes unique to this instance
        self.__dict__.update(kwargs)  # Permits direct declaration of
        # key:value pairs from declaration, e.g.
        # x = Settings(foo=1, bar=2)
        # See: https://stackoverflow.com/questions/1878710/struct-objects-in-python
        return

    # Properties getters/setters
    #
    # Remember: Sphinx ignores docstrings on property setters so all
    # documentation for a property must be on the @property method
    # Private methods
    def __str__(self):
        """Provides the string representation of the object.
        :returns: The object representation as a string.
        :rtype: str
        """
        s = "<" + self.getClassName() + ": "
        s = s + str(self.__dict__)
        s = s + ">"
        return s

    # Public methods
    def getClassName(self):
        """Retrieves the class name.
        :returns: The class name.
        :rtype: str
        """
        return type(self).__name__

    def read(self, filename):
        """Read settings from JSON file.
        :param filename: The name of the file to be read (including path)
        :type filename: str
        :returns: True if file was sucessfully read. False otherwise.
        :rtype: bool
        """
        with open(filename, "r") as file:
            contentStr = file.read()
        # c=json.loads(jsonminify(content))
        contentDict = json.loads(fastjsonminify(contentStr))
        # c contains a dictionary that has to be
        # traspassed to self
        # Loop over the dictionary
        for fieldName, fieldValue in contentDict.items():
            setattr(self, fieldName, fieldValue)
        return True

    def write(self, filename):
        """Write settings to a JSON file.
        :returns: True if file was sucessfully read. False otherwise.
        :rtype: bool
        """
        contentStr = json.dumps(self.__dict__)
        with open(filename, "w") as file:
            file.write("# \n")
            file.write("# File: " + filename + "\n")
            file.write("# \n")
            file.write("# This is an OCTant settings file.\n")
            file.write("# You can add, edit or remove settings manually here.\n")
            file.write(
                "# File format is in JSON. Although comments are permitted, but they will be lost after resaving because of minification.\n"
            )
            file.write(
                '# If you want your comments to be persistent declared them as "__comment" fields.\n'
            )
            file.write("# \n")
            file.write(
                "# File last saved: "
                + datetime.utcnow().strftime("%d-%b-%Y %H:%M:%S UTC+0")
                + "\n"
            )
            file.write("# \n")
            file.write("# (c) 2019.  Felipe Orihuela-Espina.\n")
            file.write("# \n\n")
            file.write(contentStr)
        return True


class OpScanFlatten(Operation):
    """A flattening operation for :class:`data.OCTscan`.
    A flattening operation for :class:`data.OCTscan`.
    The operation represented by this class rectifies an OCT scan.
    .. seealso:: None
    .. note:: None
    .. todo:: None
    """

    # Private class attributes shared by all instances
    # Class constructor
    def __init__(self):
        # Call superclass constructor
        super().__init__()
        # Set the operation name
        self.name = "Flattening"
        self.__deformationMap = None
        return

    @property
    def deformationMap(self):  # name getter
        """
        A logical name for the study.
        This is a read only property.
        :getter: Gets the deformationMap associated to the last flattening.
        :type: str
        """
        return self.__deformationMap

    # Private methods
    # def __str__(self):
    #     #This not working yet; it gets into an infiite recursion as
    #     #super().__str__() calls self.getClassName() in THIS class.
    #     # s = '<' + self.getClassName() + '([' \
    #     #     + super().__str__() + '])>'
    #     s = '<' + self.getClassName() + '([' \
    #             + str(super()) + '])>'
    #     print(super())
    #     return s
    @staticmethod
    def fittingQuadraticModel(x, a, b, c):
        # quadratic model for curve optimization
        return a * x * x + b * x + c

    # Public methods
    def execute(self, *args, **kwargs):
        """Executes the operation on the :py:attr:`operands`.
        Executes the operation on the :py:attr:`operands` and stores the outcome
        in :py:attr:`result`. Preload operands using
        :func:`Operation.addOperand()`.
        :returns: Result of executing the operation.
        :rtype: :class:`data.OCTscan`
        """
        # print(self._getClasName(),": flattening: Starting flattening")
        # Ensure the operand has been set.
        if len(self.operands) < 1:
            warnMsg = self.getClassName() + ":execute: Operand not set."
            warnings.warn(warnMsg, SyntaxWarning)
            return None
        imgin = self.operands[0]
        if type(imgin) is OCTscan:
            imgin = imgin.data
        # Check whether the image is in RGB (ndim=3) or in grayscale (ndim=2)
        # and convert to grayscale if necessary
        if imgin.ndim == 2:
            # Dimensions are only width and height. The image is already in grayscale.
            I2 = imgin
        elif imgin.ndim == 3:
            # Image is in RGB. Convert.
            I2 = color.rgb2gray(imgin)
        else:  # Unexpected case. Return warning
            print(self._getClasName(), ": Unexpected image shape.")
            self.result = imgin
            return self.result
        aux = np.argmax(I2, axis=0)
        mg = np.mean(aux)
        sdg = np.std(aux)
        markers = []
        remover = []
        x0 = np.arange(len(aux))
        for i in range(0, len(aux)):
            if mg - 3 * sdg <= aux[i] <= mg + 3 * sdg:
                markers += [aux[i]]
            else:
                remover += [i]
        x = np.delete(x0, remover)
        modelCoeffs, pcov = curve_fit(
            self.fittingQuadraticModel, x, markers, method="dogbox", loss="soft_l1"
        )
        a = self.fittingQuadraticModel(x0, *modelCoeffs)
        shift = np.max(a)
        flat = shift - a
        flat = np.round(flat)
        flat = np.ravel(flat).astype(int)
        self.__deformationMap = flat
        newgray = I2
        for i in range(0, len(a)):
            newgray[:, i] = np.roll(I2[:, i], flat[i], axis=0)
        self.result = OCTscan(newgray)
        return self.result

    #    #@deprecated(version='0.2', reason="Deprecated. Use method execute() instead.")
    #    @deprecation.deprecated(deprecated_in="0.2", removed_in="0.3",
    #                        current_version=__version__,
    #                        details="Use method execute() instead.")
    #    def flattening(self,image):
    #        #Encapsulate the image as an OCTscan
    #        tmp=OCTscan(image)
    #        self.clear()
    #        self.addOperand(tmp)
    #        #Execute
    #        self.execute()
    def applyOperation(self, scanA):
        """Apply the current flattening to the given scan.
        Instead of calculating the fitting again needed for the
        flattening, this method applies a known fitted quadratic model to
        the given parameters.
        The result is NOT stored in :py:attr:`result`.
        :param scanA: Image to flatten.
        :type scanA: :class:`data.OCTscan`
        :returns: Result of repeating the last flattening operation onto
             parameter scanA.
        :rtype: :class:`data.OCTscan`
        """
        if type(scanA) is OCTscan:
            scanA = scanA.data
        newgray = scanA
        for i in range(0, len(self.deformationMap)):
            newgray[:, i] = np.roll(scanA[:, i], self.deformationMap[i], axis=0)
        return OCTscan(newgray)


import warnings
from abc import ABC, abstractmethod


class Operation(ABC):
    # Sphinx documentation
    """An abstract base class for operations on :class:`octant.data.OCTvolume` and :class:`octant.data_OCTscan` .
    An abstract base class for operations on :class:`octant.data.OCTvolume` and
    :class:`octant.data.OCTscan`.
    :Example:
        tmp = octant.data.OCTscan(img)
        o = octant.op.OpScanFlattening()
        o.addOperand(tmp)
        o.arity() #return 1
        o.execute() #Flattens the image
    :Known subclasses:
    * :class:`OpScanFlattening`
    * :class:`OpScanMeasureLayerThickness`
    * :class:`OpScanPerfilometer`
    * :class:`OpScanSegment`
    * :class:`OpScanStitch`
    * :class:`OpSegmentationBrush`
    * :class:`OpSegmentationEdit`
    .. seealso:: None
    .. note:: None
    .. todo::
        * Support to hold operand names.
    """

    # Private class attributes shared by all instances
    # Class constructor
    #
    def __init__(self, **kwargs):
        """The class constructor.
        The class constructor. Creates an empty operation
        """
        super().__init__()
        # Initialize attributes (without decorator @property)
        # Initialize properties (with decorator @property)
        self.name = "Operation"  # The operation name
        self.operands = list()  # Operands
        self.parameters = list()  # Parameters
        self.result = None  # Operation outputs (a list in case it is multivalued).
        # None until executed or cleared.
        if kwargs is not None:
            for key, value in kwargs.items():
                if key == "name":
                    self.name = value
        return

    # Properties getters/setters
    #
    # Remember: Sphinx ignores docstrings on property setters so all
    # documentation for a property must be on the @property method
    @property
    def operands(self):  # operands getter
        """
        The list of operands.
        :getter: Gets the list of operands
        :setter: Sets the list of operands.
        :type: list
        """
        return self.__operands

    @operands.setter
    def operands(self, opList):  # operands setter
        # if (not isinstance(opList,(list,))):
        if type(opList) is not list:
            warnMsg = (
                self.getClassName() + ":operands: Unexpected type. "
                "Please provide operands as a list."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        else:
            self.__operands = opList
        return None

    @property
    def name(self):  # name getter
        """
        The operation name
        :getter: Gets the operation name
        :setter: Sets the operation name.
        :type: string
        """
        return self.__name

    @name.setter
    def name(self, opName):  # name setter
        # if (not isinstance(opName,(str,))):
        if type(opName) is not str:
            warnMsg = (
                self.getClassName() + ":name: Unexpected type. "
                "Operations name must be a string."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        else:
            self.__name = opName
        return None

    @property
    def parameters(self):  # operands getter
        """
        The list of parameters.
        :getter: Gets the list of parameters
        :setter: Sets the list of parameters.
        :type: list
        """
        return self.__parameters

    @parameters.setter
    def parameters(self, opList):  # operands setter
        # if (not isinstance(opList,(list,))):
        if type(opList) is not list:
            warnMsg = (
                self.getClassName() + ":parameters: Unexpected type. "
                "Please provide operands as a list."
            )
            warnings.warn(warnMsg, SyntaxWarning)
        else:
            self.__parameters = opList
        return None

    @property
    def result(self):  # result getter
        """
        The list of results.
        This is a read only property. There is no setter method.
        :getter: Gets the list of results
        :setter: Sets the list of results
        :type: list
        """
        return self.__result

    @result.setter
    def result(self, rList):  # result setter
        self.__result = rList
        return None

    # Private methods
    def __str__(self):
        tmp = "["
        for x in self.operands:
            tmp += format(x) + ","
        tmp += "]"
        s = (
            "<"
            + self.getClassName()
            + "(["
            + "name: "
            + self.name
            + ";"
            + " operands: "
            + tmp
            + "])>"
        )
        return s

    # Public methods
    def getClassName(self):
        """Get the class name as a string.
        Get the class name as a string.
        :returns: The class name.
        :rtype: string
        """
        return type(self).__name__

    def addOperand(self, op, i=None):
        """
        Add a new operand.
        :param op: The operand.
        :type op: object
        :param i: (optional) The operand order. If given it may shift the
            order of other operands already set. If not given, the operand
            is appended at the end of the list of operands.
        :type op: int
        :return: None
        """
        if i is None:
            self.__operands.append(op)
        else:
            self.__operands.insert(i, op)
        return None

    def setOperand(self, op, i):
        """
        Set an operand; substitutes an existing operand with a new one.
        Calling setOperand when the :py:attr:`i`-th operand has not been
        previously set will result in an out-of-range error.
        :param op: The new operand.
        :type op: object
        :param i: The operand order. Operand index is zero-base i.e. the
            first operand occupies i=0
        :type op: int
        :return: None
        """
        self.__operands[i] = op
        return None

    def addParameter(self, param, i=None):
        """
        Add a new parameter.
        :param op: The parameter.
        :type op: object
        :param i: (optional) The paremeter order. If given it may shift the
            order of other parameters already set. If not given, the parameter
            is appended at the end of the list of parameters.
        :type op: int
        :return: None
        """
        if i is None:
            self.__parameters.append(op)
        else:
            self.__parameters.insert(i, op)
        return None

    def setParameter(self, op, i):
        """
        Set a parameter; substitutes an existing parameter with a new one.
        Calling setParameter when the :py:attr:`i`-th parameter has not been
        previously set will result in an out-of-range error.
        :param op: The new operand.
        :type op: object
        :param i: The operand order. Operand index is zero-base i.e. the
            first operand occupies i=0
        :type op: int
        :return: None
        """
        self.__operands[i] = op
        return None

    def arity(self):
        """Gets the operation arity (number of operands).
        :return: The operation arity
        :rtype: int
        """
        return len(self.__operands)

    def clear(self):
        """
        Clears the operands; Removes all operands.
        :return: None
        """
        self.__operands = list()
        return None

    # @abstractmethod
    def execute(self, *args, **kwargs):
        """Executes the operation on the operands.
        This is an abstract method. Executes the operation on the .operands
        and stores the outcome in .result
        Operation meta-parameters may be also passed.
        :returns: Result of executing the operation.
        :rtype: Type of result -depends on subclass implementation-.
        """
        pass
