# Source code for fibtortuosity.fibtortuosity

    # Copyright 2016 Joshua Taillon
    #
    # Licensed under the Apache License, Version 2.0 (the "License");
    # you may not use this file except in compliance with the License.
    # You may obtain a copy of the License at
    #
    #     http://www.apache.org/licenses/LICENSE-2.0
    #
    # Unless required by applicable law or agreed to in writing, software
    # distributed under the License is distributed on an "AS IS" BASIS,
    # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    # See the License for the specific language governing permissions and
    # limitations under the License.

    import sys
    import time
    from datetime import datetime

    from libtiff import TIFFfile, TIFFimage
    import matplotlib
    import numpy as np
    import skfmm
    from matplotlib import pyplot as plt

    # from memory_profiler import profile
    import gc

    __version__ = '0.2.0'

    def run_full_analysis_lsm_ysz(electrolyte_file,
                                  electrolyte_and_pore_file,
                                  electrolyte_and_lsm_file,
                                  electrolyte_and_ysz_file,
                                  date,
                                  phase,
                                  direction,
                                  npzfile=None,
                                  units='nm',
                                  delay=0,
                                  calculate_all=False,
                                  load_from_prev_run=True,
                                  create_hspy_sigs=True,
                                  save_avizo_tiff=True,
                                  tort_profile=True,
                                  save_tort_prof=True,
                                  in_ipython=True,
                                  send_text=False,
                                  text_email=None,
                                  text_number=None,
                                  text_carrier=None):
        """
        Run full tortuosity analysis for LSM/YSZ, calling other functions as
        necessary.

        Parameters
        ----------
        electrolyte_file: str
            Name of bulk electrolyte LabelField
        electrolyte_and_pore_file: str
            Name of LabelField containing bulk electrolyte and pore
        electrolyte_and_lsm_file: str
            Name of LabelField containing bulk electrolyte and LSM
        electrolyte_and_ysz_file: str
            Name of LabelField containing bulk electrolyte and YSZ
        date: str
            date string to be used in output filenames
        phase: str
            phase label to be used in output filenames
        direction: str
            direction of analysis ('x', 'y', or 'z')
        npzfile: str or None
            filename of previously saved data to use (if loading from disk)
        units: str
            units for labeling purpose
        delay: int
            delay to wait before running analysis (useful if running more than
            one at a time)
        calculate_all: bool
            flag to indicate that geodesic distance etc. should be calculated
            from scratch (takes a while)
        load_from_prev_run: bool
            flag to indicate that data should be loaded from disk, using the
            file specified in ``npzfile``
        create_hspy_sigs: bool
            flag to indicate if hyperspy signals should be created (to help
            visualize the results)
        save_avizo_tiff: bool

        tort_profile: bool

        save_tort_prof: bool

        in_ipython: bool
            flag to control if notifications from ipython will be shown.
            Requires the "Notifyme!" `extension <https://gist.github.com/jat255/bde
            d66839b424e20eb8bd1dc2336f841>`_
        send_text: bool
            flag to control if text should be sent (requires |texting|_ and
            |keyring|_ packages)
        text_email: string
            email address to send text from (see |textingTextServer|_)
        text_number: string
            number to send text to (see |textingTextServer|_)
        text_carrier: string
            carrier of number to send to (see |textingTextServer|_)

        Returns
        -------
        results: dict
            results is a dictionary containing all the various parameters that
            were calculated during the analysis

        Notes
        -----
        To authenticate your email server and store the password in your local
        keyring, run the following before trying to send texts:

        >>> import keyring
        ... keyring.set_password('python_TextServer', emailaddress, password)

        Examples
        --------

        >>> res = run_full_analysis_lsm_ysz(electrolyte_file="bulkYSZ.tif",
        ...                                 electrolyte_and_pore_file="bulkYSZandPRE.tif",
        ...                                 electrolyte_and_lsm_file="bulkYSZandLSM.tif",
        ...                                 electrolyte_and_ysz_file="bulkYSZandYSZ.tif",
        ...                                 date='2015-05-06',
        ...                                 phase='Pore',
        ...                                 direction='x',
        ...                                 npzfile=None,
        ...                                 units='nm',
        ...                                 delay=0,
        ...                                 calculate_all=True,
        ...                                 load_from_prev_run=False,
        ...                                 create_hspy_sigs=True,
        ...                                 save_avizo_tiff=True,
        ...                                 tort_profile=True,
        ...                                 save_tort_prof=True,
        ...                                 in_ipython=False,
        ...                                 send_text=True,
        ...                                 text_email="email@server.com",
        ...                                 text_number="8008675309",
        ...                                 text_carrier="verizon")

        .. |texting| replace:: ``texting``
        .. _texting: https://bitbucket.org/jat255/jat255-python-modules/src/09e8acf65230dd3f84df8b25f84000d2eec02755/texting/?at=master
        .. |keyring| replace:: ``keyring``
        .. _keyring: https://pypi.python.org/pypi/keyring
        .. |textingTextServer| replace:: ``texting.TextServer``
        .. _textingTextServer: https://bitbucket.org/jat255/jat255-python-modules/src/09e8acf65230dd3f84df8b25f84000d2eec02755/texting/texting/texting.py?at=master&fileviewer=file-view-default#texting.py-7

        """
        # Dictionary for phase -> filename mapping:
        filedict = {'Pore':electrolyte_and_pore_file,
                    'LSM':electrolyte_and_lsm_file,
                    'YSZ':electrolyte_and_ysz_file}

        # Check to make sure that both calculate and load options are specified
        if calculate_all:
            if load_from_prev_run:
                print("Both loading and calculating were specified. Calculating "
                      "is given preference, so data will not be loaded from disk.")
                sys.stdout.flush()
            load_from_prev_run = False

        # Dictionary to hold all the results
        results = {}

        # Imports
        if in_ipython:
            from IPython import get_ipython
            ipython = get_ipython()

            ipython.magic('load_ext notifyme')
            ipython.magic('notifysound default')

        if send_text:
            try:
                # noinspection PyUnresolvedReferences
                from texting.texting import TextServer
                txt_srv = TextServer(text_email, text_number, text_carrier)
            except ImportError as e:
                print('Could not import texting library. Is it installed?')
                send_text = False

        def get_time():
            return time.time()

        def done_text_send(text_start_time, text_end_time, process):
            """
            Send a text message with descriptive information about calculation
            times using the ``TextServer`` class

            Parameters
            ----------
            text_start_time: datetime.time object
                Starting time of interval
            text_end_time: datetime.time object
                Ending time of interval
            process: str
                Descriptive string to include in the message (usually what the
                name of the process was)
            """
            time_used = text_end_time - text_start_time
            hours_hold = int(time_used / 3600)
            minutes_hold = int(time_used / 60 - int(hours_hold * 60))
            seconds_hold = int(time_used - int(minutes_hold * 60))
            formatted_time = str(hours_hold).zfill(2) + ':' + str(
                minutes_hold).zfill(2) + ':' + str(seconds_hold).zfill(2)
            #
            msg = ('\nDONE\nProcess: ' + str(process) + '\nTime required: ' +
                   str(formatted_time))

            # The actual mail send
            txt_srv.send_message(msg)

        # Delay (if requested):
        import time
        if delay > 0:
            print("Waiting %i seconds before running..." % delay)
            time.sleep(delay)

        # Calculating geodesic, euclidean, and tortuosity
        if calculate_all:
            fn_string = 'tortuosity_from_labels_' + direction
            fn = globals()[fn_string]
            print("Cathode data file is " + filedict[phase])
            sys.stdout.flush()

            start_time = get_time()
            g, e, t, d = fn(electrolyte_file,
                            filedict[phase],
                            phase,
                            units=units,
                            save_output=False,
                            print_output=True)
            results['g'] = g
            results['e'] = e
            results['t'] = t
            results['d'] = d
            print("Data stored in keys \'g\', \'e\', \'t\', and \'d\'")
            print("")
            if in_ipython:
                ipython.magic('notifyme $phase - $direction completed')
            if send_text:
                done_text_send(start_time,
                               get_time(),
                               phase + '-' + direction + ' done')

        # Saving results:
        if calculate_all:
            print("Saving results to disk...")
            sys.stdout.flush()
            start_time = get_time()
            save_results(fname=None,
                         geo_d=g,
                         euc_d=e,
                         tort=t,
                         desc=d,
                         phase=phase,
                         direction=direction)
            if in_ipython:
                ipython.magic('notifyme $phase - $direction saved')
            if send_text:
                done_text_send(start_time,
                               get_time(),
                               phase + '-' + direction + ' saved')

        # Loading data (after, instead of calculating from scratch):
        if load_from_prev_run:
            if npzfile is None:
                print("No .npz file was supplied for loading...")
                raise ValueError('npzfile not defined')
            print("Loading data from previous run...")
            sys.stdout.flush()
            g, e, t, d = load_results(npzfile)
            print("Old data loaded from " + npzfile)
            results['g'] = g
            results['e'] = e
            results['t'] = t
            results['d'] = d
            print("Stored in keys \'g\', \'e\', \'t\', and \'d\'")
            print("")
            sys.stdout.flush()

        # Visualizing results (requires hyperspy):
        if create_hspy_sigs:
            import hyperspy.signals as signals

            g_s = signals.Image(g)
            g_s.metadata.General.title = phase + "_" + direction + " Geodesic"

            e_s = signals.Image(e)
            e_s.metadata.General.title = phase + "_" + direction + " Euclidean"

            t_s = signals.Image(t)
            t_s.metadata.General.title = phase + "_" + direction + " Tortuosity"

            # for i in [g_s, e_s, t_s]:
            #   i.plot()

            results['g_s'] = g_s
            results['e_s'] = e_s
            results['t_s'] = t_s
            print("Hyperspy signals are in values \'g_s\', \'e_s\', and \'t_s\'")
            print("")
            sys.stdout.flush()

        # Saving as tiff for Avizo:
        if save_avizo_tiff:
            save_as_tiff(date + "_" + phase + "_" + direction + "-tort.tif",
                         t,
                         'float32',
                         d)
            print("")
            sys.stdout.flush()

        # Average tortuosity profile:
        if tort_profile:
            print("Calculating average tortuosity profile")
            sys.stdout.flush()

            t_avg, e_avg = tortuosity_profile(t,
                                              e,
                                              axis=direction)
            f = plot_tort_prof(t_avg,
                               e_avg,
                               direction,
                               sns_fontscale=2,
                               sns_color_num=1)
            plt.tight_layout()
            matplotlib.rc('image', cmap='gray')

            results['f'] = f
            results['t_avg'] = t_avg
            results['e_avg'] = e_avg
            print("Average tortuosity plot is in variable \'f\', profiles are "
                  "in \'t_avg\' and \'e_avg\'")
            print("")
            sys.stdout.flush()

        # Save profile to disk:
        if save_tort_prof:
            save_profile_to_csv(date + '_' +
                                phase + '_' +
                                direction + '-tort-profile.csv',
                                e_avg / 1000,  # divide by 1000 to get microns
                                t_avg,
                                direction,
                                phase)

        return results

    def calculate_geodesic_distance(electrolyte_fname,
                                    phase_fname,
                                    units='nm',
                                    print_output=True,
                                    ):
        """
        Calculate the geodesic distance of an indexed tif file.

        .. note:: Deprecated. This function will work, but it is better to use
                  :func:`tortuosity_from_labels_x`,
                  :func:`tortuosity_from_labels_y`,
                  and :func:`tortuosity_from_labels_z`.

        Parameters
        ----------
        electrolyte_fname: str
            Filename of a 3D tiff file that contains only the bulk electrolyte
            labelfield (as 1 values)
        phase_fname: str
            Filename of a 3D tiff file that contains the bulk electrolyte and the
            phase for which the geodesic distance is to be calculated
        units: str
            units of the given dimensions
        print_output: bool
            switch to control whether output is printed

        Returns
        -------
        d: np.array
            Numpy array with geodesic distance from supplied surface
        desc: str
            description string to be used when saving the file (for proper
            reading into Avizo)
        """

        # Load electrolyte data
        electrolyte_tif = TIFFfile(electrolyte_fname)
        electrolyte_array = electrolyte_tif.get_tiff_array()
        electrolyte_data = electrolyte_array[::].astype('int8')
        if print_output:
            print('Loaded electrolyte data.')

        # Load phase to be calculated data
        cathode_tif = TIFFfile(phase_fname)
        cathode_array = cathode_tif.get_tiff_array()
        cathode_data = cathode_array[::].astype('int8')
        if print_output:
            print('Loaded cathode data.')

        # Get information from Avizo-saved tiff file
        ifd = electrolyte_tif.get_first_ifd()
        size_x = ifd.get_value('ImageWidth')
        size_y = ifd.get_value('ImageLength')
        size_z = electrolyte_tif.get_depth()
        image_description = ifd.get_value('ImageDescription')
        if print_output:
            print('ImageDescription is: {}'.format(image_description))

        # Calculate some parameters from the image input
        bb_xmin, bb_xmax, bb_ymin, bb_ymax, bb_zmin, bb_zmax = \
            [eval(i) for i in image_description.split()[1:]]
        bound_x = bb_xmax - bb_xmin
        bound_y = bb_ymax - bb_ymin
        bound_z = bb_zmax - bb_zmin
        if print_output:
            print("Bounding box dimensions are: ({0:.2f}, {1:.2f}, "
                  "{2:.2f})".format(bound_x, bound_y, bound_z) + " " + units)

        vox_x = (bb_xmax - bb_xmin) / float(size_x - 1)
        vox_y = (bb_ymax - bb_ymin) / float(size_y - 1)
        vox_z = (bb_zmax - bb_zmin) / float(size_z - 1)
        if print_output:
            print("Voxel dimensions are: ({0:.2f}, {1:.2f}, "
                  "{2:.2f})".format(vox_x, vox_y, vox_z) + " " + units)

        # Switch electrolyte tag (1) to a negative value to create zero-contour
        cathode_data[electrolyte_data == 1] = -1

        # Mask the voxels where phase is not present:
        mask = cathode_data == 0
        cathode_data = np.ma.MaskedArray(cathode_data, mask)

        # Do geodesic distance calculation
        start_time = datetime.now()
        if print_output:
            print('Starting calculation at: {}'.format(start_time))
        # Calculate geodesic distance:
        d = skfmm.distance(cathode_data, dx=[vox_z, vox_y, vox_x])
        end_time = datetime.now()
        if print_output:
            print('Calculation took: {}'.format(end_time - start_time))

        # Set distances within the bulk electrolyte to zero:
        # d.data[electrolyte_data == 1] = 0

        return d.data, image_description

    def tortuosity_from_labels_x(electrolyte_fname,
                                 phase_fname,
                                 phase,
                                 units='nm',
                                 print_output=True,
                                 save_output=False):
        """
        Calculate the tortuosity (normal to electrolyte interface, x-direction)
        from supplied label fields

        Parameters
        ----------
        electrolyte_fname: str
            Filename of a 3D tiff file that contains only the bulk electrolyte
            labelfield (as 1 values)
        phase_fname: str
            Filename of a 3D tiff file that contains the bulk electrolyte and the
            phase for which the geodesic distance is to be calculated
        phase: str
            Phase that is being calculated. i.e. 'pore', 'LSM', 'YSZ', etc.
        units: str
            units of the given dimensions
        print_output: bool
            switch to control whether output is printed
        save_output: bool
            switch to control whether data is saved to .npy files for later recall

        Returns
        --------
        geo_d: np.array
            Numpy array with geodesic distance from supplied surface
        euc_d: np.array
            Numpy array with euclidean distance from electrolyte
        tort: np.array
            Numpy array with tortuosity at each point
        desc: str
            description string to be used when saving the file (for proper
            reading into Avizo)
        """

        if print_output:
            print('Starting calculation on ' + phase + ' in x direction.')
            sys.stdout.flush()

        global_start_time = datetime.now()

        # Load electrolyte data
        start_time = datetime.now()
        electrolyte_tif = TIFFfile(electrolyte_fname)
        electrolyte_array = electrolyte_tif.get_tiff_array()
        electrolyte_data = electrolyte_array[::].astype('int8')
        end_time = datetime.now()
        if print_output:
            print('Loaded electrolyte data in {} seconds.'.format(end_time -
                                                                  start_time))
            sys.stdout.flush()

        # Load phase to be calculated data
        start_time = datetime.now()
        cathode_tif = TIFFfile(phase_fname)
        cathode_array = cathode_tif.get_tiff_array()
        cathode_data = cathode_array[::].astype('int8')
        end_time = datetime.now()
        if print_output:
            print('Loaded cathode data in {} seconds.'.format(end_time -
                                                              start_time))
            sys.stdout.flush()

        # Get information from Avizo-saved tiff file
        ifd = electrolyte_tif.get_first_ifd()
        size_x = ifd.get_value('ImageWidth')
        size_y = ifd.get_value('ImageLength')
        size_z = electrolyte_tif.get_depth()
        image_description = ifd.get_value('ImageDescription')
        if print_output:
            print('ImageDescription is: {}'.format(image_description))
            sys.stdout.flush()

        # Calculate some parameters from the image input
        bb_xmin, bb_xmax, bb_ymin, bb_ymax, bb_zmin, bb_zmax = \
            [eval(i) for i in image_description.split()[1:]]
        bound_x = bb_xmax - bb_xmin
        bound_y = bb_ymax - bb_ymin
        bound_z = bb_zmax - bb_zmin
        if print_output:
            print("Bounding box dimensions are: ({0:.2f}, {1:.2f}, "
                  "{2:.2f})".format(bound_x, bound_y, bound_z) + " " + units)
            sys.stdout.flush()

        vox_x = (bb_xmax - bb_xmin) / float(size_x - 1)
        vox_y = (bb_ymax - bb_ymin) / float(size_y - 1)
        vox_z = (bb_zmax - bb_zmin) / float(size_z - 1)
        if print_output:
            print("Voxel dimensions are: ({0:.2f}, {1:.2f}, "
                  "{2:.2f})".format(vox_x, vox_y, vox_z) + " " + units)
            sys.stdout.flush()

        # Switch electrolyte tag (1) to a negative value to create zero-contour
        cathode_data[electrolyte_data == 1] = -1

        # Mask the voxels where phase is not present:
        cathode_mask = cathode_data == 0

        cathode_ma = np.ma.MaskedArray(cathode_data, cathode_mask)
        del cathode_mask
        gc.collect()

        # Do geodesic distance calculation
        start_time = datetime.now()
        if print_output:
            print('Starting geodesic calculation at: {}'.format(start_time))
            sys.stdout.flush()
        geo_d = _geo_dist(cathode_ma, vox_x, vox_y, vox_z)
        end_time = datetime.now()
        if print_output:
            print('Geodesic calculation took: {}'.format(end_time -
                                                         start_time))
            sys.stdout.flush()

        # Calculating zero reference for each [x,y] point for euclidean distance
        start_time = datetime.now()
        x_interface = _calc_interface(electrolyte_data)
        end_time = datetime.now()
        if print_output:
            print("Calculating zero distance interface took: "
                  "{}".format(end_time - start_time))
            sys.stdout.flush()

        # Calculate euclidean distance
        start_time = datetime.now()
        # noinspection PyTypeChecker
        euc_d = _calc_euc_x(x_interface, electrolyte_data, vox_x)
        end_time = datetime.now()
        if print_output:
            print("Calculating euclidean distance took: {}".format(end_time -
                                                                   start_time))
            sys.stdout.flush()

        # Calculate tortuosity
        start_time = datetime.now()
        tort = _calc_tort(geo_d, euc_d, electrolyte_data)
        end_time = datetime.now()
        if print_output:
            print("Calculating tortuosity took: {}".format(end_time - start_time))
            sys.stdout.flush()

        start_time = datetime.now()
        if save_output:
            if print_output:
                print("Saving data...")
                sys.stdout.flush()
            save_results(fname=None,
                         direction='x',
                         phase=phase,
                         geo_d=geo_d.data,
                         euc_d=euc_d,
                         tort=tort.data,
                         desc=image_description)
            end_time = datetime.now()
            if print_output:
                print("Saving data took: {}".format(end_time - start_time))
                sys.stdout.flush()

        end_time = datetime.now()
        if print_output:
            print("Total execution time was: {}".format(end_time -
                                                        global_start_time))
            sys.stdout.flush()

        return geo_d.data, euc_d, tort.data, image_description

    # @profile
    def tortuosity_from_labels_y(electrolyte_fname,
                                 phase_fname,
                                 phase,
                                 units='nm',
                                 print_output=True,
                                 save_output=False):
        """
        Calculate the tortuosity (parallel to electrolyte interface, y-direction)
        from supplied label fields

        Parameters
        ----------
        electrolyte_fname: str
            Filename of a 3D tiff file that contains only the bulk electrolyte
            labelfield (as 1 values)
        phase_fname: str
            Filename of a 3D tiff file that contains the bulk electrolyte and the
            phase for which the geodesic distance is to be calculated
        phase: str
            Phase that is being calculated. i.e. 'pore', 'LSM', 'YSZ', etc.
        units: str
            units of the given dimensions
        print_output: bool
            switch to control whether output is printed
        save_output: bool
            switch to control whether data is saved to .npy files for later recall

        Returns
        --------
        geo_d: np.array
            Numpy array with geodesic distance from supplied surface
        euc_d: np.array
            Numpy array with euclidean distance from electrolyte
        tort: np.array
            Numpy array with tortuosity at each point
        desc: str
            description string to be used when saving the file (for proper
            reading into Avizo)
        """

        global_start_time = datetime.now()

        if print_output:
            print('Starting calculation on ' + phase + ' in y direction.')
            sys.stdout.flush()

        # Load electrolyte data
        start_time = datetime.now()
        electrolyte_tif = TIFFfile(electrolyte_fname)
        electrolyte_array = electrolyte_tif.get_tiff_array()
        electrolyte_data = electrolyte_array[::].astype('int8')
        end_time = datetime.now()
        if print_output:
            print('Loaded electrolyte data in {} seconds.'.format(end_time -
                                                                  start_time))
            sys.stdout.flush()

        # Load phase to be calculated data
        start_time = datetime.now()
        cathode_tif = TIFFfile(phase_fname)
        cathode_array = cathode_tif.get_tiff_array()
        cathode_data = cathode_array[::].astype('int8')
        end_time = datetime.now()
        if print_output:
            print('Loaded cathode data in {} seconds.'.format(end_time -
                                                              start_time))
            sys.stdout.flush()

        # Get information from Avizo-saved tiff file
        ifd = electrolyte_tif.get_first_ifd()
        size_x = ifd.get_value('ImageWidth')
        size_y = ifd.get_value('ImageLength')
        size_z = electrolyte_tif.get_depth()
        image_description = ifd.get_value('ImageDescription')
        if print_output:
            print('ImageDescription is: ' + image_description)
            sys.stdout.flush()

        # Calculate some parameters from the image input
        bb_xmin, bb_xmax, bb_ymin, bb_ymax, bb_zmin, bb_zmax = \
            [eval(i) for i in image_description.split()[1:]]
        bound_x = bb_xmax - bb_xmin
        bound_y = bb_ymax - bb_ymin
        bound_z = bb_zmax - bb_zmin
        if print_output:
            print("Bounding box dimensions are: ({0:.2f}, {1:.2f}, "
                  "{2:.2f})".format(bound_x, bound_y, bound_z) + " " + units)
            sys.stdout.flush()

        vox_x = (bb_xmax - bb_xmin) / float(size_x - 1)
        vox_y = (bb_ymax - bb_ymin) / float(size_y - 1)
        vox_z = (bb_zmax - bb_zmin) / float(size_z - 1)
        if print_output:
            print("Voxel dimensions are: ({0:.2f}, {1:.2f}, "
                  "{2:.2f})".format(vox_x, vox_y, vox_z) + " " + units)
            sys.stdout.flush()

        # Set electrolyte area to zero in cathode data
        cathode_data[cathode_data == 1] = 0

        # Mask the voxels where phase is not present:
        cathode_mask = cathode_data == 0

        # Switch XZ zero-plane to a negative value to create zero-contour
        cathode_data[:,0,:] = -1

        cathode_ma = np.ma.MaskedArray(cathode_data, cathode_mask)
        del cathode_mask
        gc.collect()

        # Do geodesic distance calculation
        start_time = datetime.now()
        if print_output:
            print('Starting geodesic calculation at: {}'.format(start_time))
            sys.stdout.flush()
        geo_d = _geo_dist(cathode_ma, vox_x, vox_y, vox_z)
        end_time = datetime.now()
        if print_output:
            print('Geodesic calculation took: {}'.format(end_time -
                                                         start_time))
            sys.stdout.flush()

        # Do not need to calculating zero reference for the y-direction (use
        # just plane)
        # start_time = datetime.now()
        # x_interface = _calc_interface(electrolyte_data)
        # end_time = datetime.now()
        # if print_output:
        #     print("Calculating zero distance interface took: "
        #           "{}".format(end_time - start_time))

        # Calculate euclidean distance
        start_time = datetime.now()

        y_interface = np.zeros_like(electrolyte_data[:,0,:], dtype=np.float32)
        euc_d = np.zeros_like(electrolyte_data, dtype=np.float32)

        for (z,x), value in np.ndenumerate(y_interface):
            # set the value at each [x,z] point to a numpy array that ranges the
            # z dimension times the vox_y size
            # so this value can be compared with that calculated geodesically
            euc_d[z,:,x] = np.arange(len(euc_d[z,:,x])) * vox_y

        euc_d[euc_d < 0] = 0    # set all negative values to zero,
        #                         since they are on the other side of the
        #                         electrolyte (probably don't need this for y-dir)
        end_time = datetime.now()
        if print_output:
            print("Calculating euclidean distance took: {}".format(end_time -
                                                                   start_time))
            sys.stdout.flush()

        # Calculate tortuosity
        start_time = datetime.now()
        tort = _calc_tort(geo_d, euc_d, electrolyte_data)
        end_time = datetime.now()
        if print_output:
            print("Calculating tortuosity took: {}".format(end_time - start_time))
            sys.stdout.flush()

        start_time = datetime.now()
        if save_output:
            if print_output:
                print("Saving data...")
                sys.stdout.flush()
            save_results(fname=None,
                         direction='y',
                         phase=phase,
                         geo_d=geo_d.data,
                         euc_d=euc_d,
                         tort=tort.data,
                         desc=image_description)
            end_time = datetime.now()
            if print_output:
                print("Saving data took: {}".format(end_time - start_time))
                sys.stdout.flush()

        end_time = datetime.now()
        if print_output:
            print("Total execution time was: {}".format(end_time -
                                                        global_start_time))
            sys.stdout.flush()

        return geo_d.data, euc_d, tort.data, image_description

    def tortuosity_from_labels_z(electrolyte_fname,
                                 phase_fname,
                                 phase,
                                 units='nm',
                                 print_output=True,
                                 save_output=False):
        """
        Calculate the tortuosity (parallel to electrolyte interface, z-direction)
        from supplied label fields

        Parameters
        ----------
        electrolyte_fname: str
            Filename of a 3D tiff file that contains only the bulk electrolyte
            labelfield (as 1 values)
        phase_fname: str
            Filename of a 3D tiff file that contains the bulk electrolyte and the
            phase for which the geodesic distance is to be calculated
        phase: str
            Phase that is being calculated. i.e. 'pore', 'LSM', 'YSZ', etc.
        units: str
            units of the given dimensions
        print_output: bool
            switch to control whether output is printed
        save_output: bool
            switch to control whether data is saved to .npy files for later recall

        Returns
        --------
        geo_d: np.array
            Numpy array with geodesic distance from supplied surface
        euc_d: np.array
            Numpy array with euclidean distance from electrolyte
        tort: np.array
            Numpy array with tortuosity at each point
        desc: str
            description string to be used when saving the file (for proper
            reading into Avizo)
        """

        global_start_time = datetime.now()

        if print_output:
            print('Starting calculation on ' + phase + ' in z direction.')
            sys.stdout.flush()

        # Load electrolyte data
        start_time = datetime.now()
        electrolyte_tif = TIFFfile(electrolyte_fname)
        electrolyte_array = electrolyte_tif.get_tiff_array()
        electrolyte_data = electrolyte_array[::].astype('int8')
        end_time = datetime.now()
        if print_output:
            print('Loaded electrolyte data in {} seconds.'.format(end_time -
                                                                  start_time))
            sys.stdout.flush()

        # Load phase to be calculated data
        start_time = datetime.now()
        cathode_tif = TIFFfile(phase_fname)
        cathode_array = cathode_tif.get_tiff_array()
        cathode_data = cathode_array[::].astype('int8')
        end_time = datetime.now()
        if print_output:
            print('Loaded cathode data in {} seconds.'.format(end_time -
                                                              start_time))
            sys.stdout.flush()

        # Get information from Avizo-saved tiff file
        ifd = electrolyte_tif.get_first_ifd()
        size_x = ifd.get_value('ImageWidth')
        size_y = ifd.get_value('ImageLength')
        size_z = electrolyte_tif.get_depth()
        image_description = ifd.get_value('ImageDescription')
        if print_output:
            print('ImageDescription is: ' + image_description)
            sys.stdout.flush()

        # Calculate some parameters from the image input
        bb_xmin, bb_xmax, bb_ymin, bb_ymax, bb_zmin, bb_zmax = \
            [eval(i) for i in image_description.split()[1:]]
        bound_x = bb_xmax - bb_xmin
        bound_y = bb_ymax - bb_ymin
        bound_z = bb_zmax - bb_zmin
        if print_output:
            print("Bounding box dimensions are: ({0:.2f}, {1:.2f}, "
                  "{2:.2f})".format(bound_x, bound_y, bound_z) + " " + units)
            sys.stdout.flush()

        vox_x = (bb_xmax - bb_xmin) / float(size_x - 1)
        vox_y = (bb_ymax - bb_ymin) / float(size_y - 1)
        vox_z = (bb_zmax - bb_zmin) / float(size_z - 1)
        if print_output:
            print("Voxel dimensions are: ({0:.2f}, {1:.2f}, "
                  "{2:.2f})".format(vox_x, vox_y, vox_z) + " " + units)
            sys.stdout.flush()

        # Set electrolyte area to zero in cathode data
        cathode_data[cathode_data == 1] = 0

        # Mask the voxels where phase is not present:
        cathode_mask = cathode_data == 0

        # Switch XY zero-plane to a negative value to create zero-contour
        cathode_data[0,:,:] = -1

        cathode_ma = np.ma.MaskedArray(cathode_data, cathode_mask)
        del cathode_mask
        gc.collect()

        # Do geodesic distance calculation
        start_time = datetime.now()
        if print_output:
            print('Starting geodesic calculation at: {}'.format(start_time))
            sys.stdout.flush()
        geo_d = _geo_dist(cathode_ma, vox_x, vox_y, vox_z)
        end_time = datetime.now()
        if print_output:
            print('Geodesic calculation took: {}'.format(end_time -
                                                         start_time))
            sys.stdout.flush()

        # Do not need to calculating zero reference for the y-direction (use
        # just plane)
        # start_time = datetime.now()
        # x_interface = _calc_interface(electrolyte_data)
        # end_time = datetime.now()
        # if print_output:
        #     print("Calculating zero distance interface took: "
        #           "{}".format(end_time - start_time))

        # Calculate euclidean distance
        start_time = datetime.now()

        z_interface = np.zeros_like(electrolyte_data[0,:,:], dtype=np.float32)
        euc_d = np.zeros_like(electrolyte_data, dtype=np.float32)

        for (y,x), value in np.ndenumerate(z_interface):
            # set the value at each [x,z] point to a numpy array that ranges the
            # z dimension times the vox_z size
            # so this value can be compared with that calculated geodesically
            euc_d[:,y,x] = np.arange(len(euc_d[:,y,x])) * vox_z

        euc_d[euc_d < 0] = 0    # set all negative values to zero,
        #                         since they are on the other side of the
        #                         electrolyte (probably don't need this for y-dir)
        end_time = datetime.now()
        if print_output:
            print("Calculating euclidean distance took: {}".format(end_time -
                                                                   start_time))
            sys.stdout.flush()

        # Calculate tortuosity
        start_time = datetime.now()
        tort = _calc_tort(geo_d, euc_d, electrolyte_data)
        end_time = datetime.now()
        if print_output:
            print("Calculating tortuosity took: {}".format(end_time - start_time))
            sys.stdout.flush()

        start_time = datetime.now()
        if save_output:
            if print_output:
                print("Saving data...")
                sys.stdout.flush()
            save_results(fname=None,
                         direction='z',
                         phase=phase,
                         geo_d=geo_d.data,
                         euc_d=euc_d,
                         tort=tort.data,
                         desc=image_description)
            end_time = datetime.now()
            if print_output:
                print("Saving data took: {}".format(end_time - start_time))
                sys.stdout.flush()

        end_time = datetime.now()
        if print_output:
            print("Total execution time was: {}".format(end_time -
                                                        global_start_time))
            sys.stdout.flush()

        return geo_d.data, euc_d, tort.data, image_description

    def save_results(fname=None,
                     phase="",
                     direction="",
                     geo_d=None,
                     euc_d=None,
                     tort=None,
                     desc=None):
        """
        Saves all the results into a compressed numpy archive

        Parameters
        -----------
        geo_d: numpy array
            Geodesic distance array
        euc_d: numpy array
            Euclidean distance array
        tort: numpy array
            Tortuosity array
        desc: str
            Image description string
        """
        if fname is None:
            fname = time.strftime("%Y-%m-%d_%H.%M.%S_" + phase + '_' + direction +
                                  '_all-results.npz')

        np.savez_compressed(fname,
                            geo_d=geo_d,
                            euc_d=euc_d,
                            tort=tort,
                            image_description=desc
                            )

        print("Saved results to " + fname)

    def load_results(fname):
        """
        Loads the results that have been previously saved with :func:`save_results`

        Parameters
        -----------
        fname: str
            filename to load

        Returns
        --------
        geo_d: numpy array
            Geodesic distance array
        euc_d: numpy array
            Euclidean distance array
        tort: numpy array
            Tortuosity array
        desc: str
            Image description string
        """
        npzfile = np.load(fname)
        geo_d = npzfile['geo_d']
        euc_d = npzfile['euc_d']
        tort = npzfile['tort']
        desc = npzfile['image_description']

        # if description isn't None, cast it as string
        if desc.dtype != 'O':
            desc = str(desc)

        return geo_d, euc_d, tort, desc

    def tortuosity_profile(tort, euc, axis='x'):
        """
        Calculate the average tortuosity along an axis

        Parameters
        ----------
        tort: Numpy array (output of :func:`tortuosity_from_labels_x`...)
            Contains the tortuosity data, with 0 values in the masked areas
        euc: Numpy array (output of :func:`tortuosity_from_labels_x`...)
            Contains the euclidean distance at each point
        axis: str
            Axis along which to calculate profile. Should be 'x', 'y', or 'z'

        Returns
        -------
        tort_prof: 1D Numpy array
            1D profile of average tortuosity at each point in the array
            (ignoring the masked areas)
        euc_prof: 1D Numpy array
            1D profile of the average euclidean distance at each point
        """

        # Set all 0 tortuosity values to nans so they are ignored when calculating
        # the average:
        tort[tort == 0] = np.nan

        axis_dict = {'x':2, 'y':1, 'z':0}

        # Figure out the correct axis:
        axes = [0, 1, 2]
        axes.remove(axis_dict[axis])
        axes = tuple(axes)

        # Calculate averages along relevant axes (ignoring nans):
        euc_prof = np.nanmean(euc, axis=axes)
        tort_prof = np.nanmean(tort, axis=axes)

        return tort_prof, euc_prof

    def plot_tort_prof(tort_prof,
                       euc_prof,
                       direction,
                       units='nm',
                       figsize=None,
                       sns_style='white',
                       sns_context='paper',
                       sns_fontscale=1.5,
                       sns_cmap=None,
                       sns_color_num=0,
                       sns_kw=None,):
        """
        Plot a profile of the average tortuosity.

        Parameters
        ----------
        tort_prof: np.array
            tortuosity profile (output of :func:`tortuosity_profile`)
        euc_prof: np.array
            euclidean profile (output of :func:`tortuosity_profile`)
        direction: str
            direction for labeling
        units: str
            units of distance for labeling
        figsize : tuple of integers, optional, default: None
            width, height in inches. If not provided, defaults to rc
            figure.figsize.
        sns_style: str
            default style for seaborn ('white', 'dark', 'darkgrid', etc.)
        sns_context: str
            context for plotting with seaborn
        sns_fontscale: float
            font scale to pass to seaborn
        sns_cmap: list, str, tuple
            colormap for use with seaborn
            default is [light blue, dark teal, blue, red, green, orange]
            If str or tuple, then these parameters are passed to the
            seaborn.color_palette() function
        sns_color_num: int
            index of color (in ``sns_cmap``) to use
        sns_kw: dict
            additional arguments to be passed to seaborn.set_style
            (see http://goo.gl/WvLdc6 for details)

        Returns
        -------
        matplotlib figure, <matplotlib axes>
            Handle to figure that is plotted (and axes array if subplots)
        """
        # Import and setup seaborn styles
        if sns_kw is None:
            sns_kw = {}
        if sns_cmap is None:
            sns_cmap = ["#4C72B1", "#004040", "#023FA5", "#8E063B",
                        "#098C09", "#EF9708"]
        import seaborn as sns
        matplotlib.rcParams['mathtext.fontset'] = 'stixsans'
        sns.set(rc={'image.cmap': 'Greys_r'})
        sns.set_style(sns_style, sns_kw)
        sns.set_context(sns_context, font_scale=sns_fontscale)
        if type(sns_cmap) is tuple:
            palette = sns.color_palette(*sns_cmap)
        else:
            palette = sns.color_palette(sns_cmap)
        sns.color_palette(palette)

        if units is 'um':
            units = "${\mu}m$"

        fig = plt.figure(figsize=figsize)
        plt.plot(euc_prof, tort_prof, color=palette[sns_color_num])
        plt.xlabel('Euclidean distance (' + direction + " direction, " +
                   units + ")")
        plt.ylabel("Average tortuosity")
        plt.tight_layout()

        return fig

    def save_profile_to_csv(fname,
                            x,
                            y,
                            direction,
                            phase,
                            x_name='Euc_d',
                            y_name='tort',
                            fmt='%10.5f',
                            delimiter=',',
                            ):
        """
        Saves two profiles (x and y) into a text file.

        Parameters
        ----------
        fname: str
            filename to save to (no checks are done before overwriting)
        x: 1D np.array
            X-profile that will be used as first column
            should be a one dimensional numpy row vector
        y: np.array or list of np.arrays
            can be 1D np.array, or a list of them, with data in rows.
            Each row will be transformed to a column in the text files
        direction: str
            direction for labeling
        phase: str
            phase for labeling
        x_name: str
            Name of header in x-column
        y_name: str
            Name of header in y-column
        fmt: str
            Format for decimals to use
        delimiter: str
            Delimiter to use in the csv file
        """
        np.savetxt(fname, np.vstack((x,y)).transpose(),
                   header=phase + " " + direction + "-direction \n" + x_name +
                          delimiter + '\t' + y_name,
                   fmt=fmt,
                   delimiter=delimiter)
        print('Saved', fname)

    def load_profile_from_csv(fname,
                              delimiter=',',
                              skiprows=2):
        """
        Loads columns of a text file into 1D numpy arrays.

        Parameters
        ----------
        fname: str
            filename to load from
        delimiter: str
            character to use as delimiter
        skiprows: int
            number of header rows to skip at beginning of file

        Returns
        --------
        list of 1D np.array objects, 1 for each column in the data
        """
        data = np.loadtxt(fname,
                          delimiter=delimiter,
                          skiprows=skiprows).transpose()
        return [i[0] for i in np.vsplit(data, len(data))]

    def save_as_tiff(fname, data, dtype, desc):
        """
        Save a numpy array as tiff (useful for transferring data back to Avizo)

        Parameters
        ----------
        fname: str
            Filename to save to
        data: np.array
            Data array to be written to image file
        dtype: str
            data type to save as. For some reason, euclidean distance works
            best as 'uint32', while the others are fine as 'float32'.
        desc: str
            Image description that includes bounding box info
        """
        tiff = TIFFimage(data.astype(dtype),
                         description=desc)
        tiff.write_file(fname, compression='lzw')

    def _geo_dist(data, vox_x, vox_y, vox_z):
        """
        Internal helper method to call scikit-fmm distance method

        Parameters
        ----------
        data: masked data array to calculate on
        vox_x: x voxel size
        vox_y: y voxel size
        vox_z: z voxel size

        Returns
        -------
        d: distance transform as calculated by scikit-fmm
        """
        d = skfmm.distance(data, dx=[vox_z, vox_y, vox_x])

        # Convert dtype from float64 to float32 to save memory (almost half):
        d = np.ma.asarray(d, dtype='float32')

        # Set distances within the bulk electrolyte to zero:
        d.data[d.data < 0] = 0
        return d

    def _calc_interface(electrolyte_data):
        """
        Calculates interface between bulk electrolyte and the cathode

        Parameters
        ----------
        electrolyte_data: label data for bulk electrolyte

        Returns
        -------
        x_interface: 2D numpy array
            x values of interface at each y and z locations
        """
        # create a single plane array with same x and y dimensions
        x_interface = np.zeros_like(electrolyte_data[:, :, 0], dtype=int)
        # visit each point on the array
        for (x, y), value in np.ndenumerate(x_interface):
            # set the value at each point to the first place where the electrolyte
            # value is 0 (which is the boundary between electrolyte and cathode)
            x_interface[x, y] = np.where(electrolyte_data[x, y, :] == 0)[0][0]

        return x_interface

    def _calc_tort(geo_d, euc_d, electrolyte_data):
        """
        Calculates tortuosity from geodesic and euclidean distances

        Parameters
        ----------
        geo_d: Geodesic distance (Masked Array)
        euc_d: Euclidean distance (Numpy array)
        electrolyte_data: label data for bulk electrolyte

        Returns
        -------
        tort: tortuosity numpy array
        """

        old = np.seterr(divide='ignore')

        # Calculate tortuosity by dividing geo_d by euc_d
        tort = np.divide(geo_d, euc_d)
        # Clean up the data, removing values < 1 and any infinities
        tort[electrolyte_data == 1] = 0     # 0 so it will be ignored in bulk YSZ
        tort[np.logical_and(tort < 1, tort > 0)] = 1
        tort[np.isinf(tort)] = 1
        tort[np.isnan(tort)] = 1

        # tort should never be very large (greater than 2), so we can
        # safely represent it with a half precision float
        # 53.7MiB before, 44.0MiB after, so saving 18%
        # tort = np.ma.asarray(tort, dtype='float16')

        np.seterr(**old)

        return tort

    def _calc_euc_x(x_interface, electrolyte_data, vox_x):
        """
        Calculates euclidean distance from electrolyte interface

        Parameters
        ----------
        x_interface: Numpy array containing x value of boundary for electrolyte
            for all y and z coordinates
        electrolyte_data: label data for bulk electrolyte
        vox_x: float
            voxel x-size

        Returns
        -------
        euc_d: euclidean distance numpy array
        """
        # create euclidean distance array of the same shape as the others
        euc_d = np.zeros_like(electrolyte_data, dtype=np.float32)
        # visit each [x,y] location in this array:
        for (x,y), value in np.ndenumerate(x_interface):
            # set the value at each [x,y] point to a numpy array that ranges the
            # x dimension minus the zero reference from before times the vox_x size
            # so this value can be compared with that calculated geodesically
            euc_d[x,y] = (np.arange(len(euc_d[x,y])) - x_interface[x,y]) * vox_x

        euc_d[euc_d < 0] = 0    # set all negative values to zero,
        #                         since they are on the other side of the
        #                         electrolyte

        return euc_d

        # How to create multiple processes for jobs:
        ############################################
        #
        # # Setup parallel stuff:
        # manager = multiprocessing.Manager()
        # return_dict = manager.dict()
        # jobs = []
        #
        # Calculate geodesic distance:
        # p = multiprocessing.Process(target=_geo_dist, args=('geo',
        #                                                     return_dict,
        #                                                     cathode_data,
        #                                                     vox_x,
        #                                                     vox_y,
        #                                                     vox_z))
        # jobs.append(p)
        # p.start()
        # if print_output:
        #     print('Started geo calc')
        #
        # Calculate euclidean distance:
        # p = multiprocessing.Process(target=_geo_dist, args=('euc',
        #                                                     return_dict,
        #                                                     all_data,
        #                                                     vox_x,
        #                                                     vox_y,
        #                                                     vox_z))
        # jobs.append(p)
        # p.start()
        # if print_output:
        #     print('Started euc calc')
        #
        # if print_output:
        #     print('Jobs list:' + repr(jobs))
        #
        # for proc in jobs:
        #     proc.join()
        #
        # geo_d = return_dict['geo']
        # euc_d = return_dict['euc']