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'


[docs]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
[docs]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
[docs]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
[docs]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
[docs]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
[docs]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)
[docs]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
[docs]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
[docs]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
[docs]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)
[docs]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))]
[docs]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']