import os
import numpy as np
import pandas as pd
from scipy.spatial.distance import pdist, squareform
import statsmodels.api as sm
from tqdm import tqdm
from .utils import calculate_errors
__all__ = ['bootstrap_tpb_stats',
'read_mv3d',
'write_mv3d',
'crop_tpb_data',
'get_bb_from_data',
'bb_volume',
'get_random_subvolume',
'path_length',
'network_length',
'split_paths_in_network',
'animate_cropped_data',
'get_bb_lines',
'get_box_and_corners',
]
# noinspection PyUnusedLocal
[docs]def bootstrap_tpb_stats(inputs_dict=None,
box_size=4000,
n_volumes=500,
n_bootstrap=100000,
save_output=False,
data_output_fname=None,
err_output_fname=None,
output_avg=False):
"""
Calculate statistics various TPB properties (and their errors) using a
random subvolume sampling method. Total TPB length (for the subvolume),
TPB density, and average length of a TPB path will be calculated.
Parameters
----------
inputs_dict: dict
dictionary of values describing the input data. Keys should be
labels for a particular set of TPB paths, while the values should be
filenames for mv3d skeleton files of the TPB paths.
box_size: float
length of edge of cube used to define the subsampled volumes. The
total volume sampled in each trial will be box_size**3, and the
boxes will be selected randomly throughout the volume of data
n_volumes: int
number of subvolumes to sample from the volume (usually ~500 or so)
n_bootstrap: int
number of bootstrap samples to use when calculating confidence
intervals
save_output: bool
switch to control whether or not the bootstrap data and error output is
written directly to a CSV file in the current directory
data_output_fname: str
filename to use when saving the data output
err_output_fname: str
filename to use when saving the error output
output_avg: bool
switch to control whether "Avg TPB path length" will be calculated
Returns
-------
data_df: :py:class:`~pandas.DataFrame`
Dataframe with data from subvolume statistic calculations
error_df: :py:class:`~pandas.DataFrame`
Dataframe with low and high errors (and std. dev.) calculated using
n_bootstrap samples
Example
-------
>>> from FIBbootstrap.tpb import bootstrap_tpb_stats
>>> inputs_dict = {
... 'active':'smoothActive.savg.mv3d',
... 'inactive':'smoothInactive.savg.mv3d',
... 'unknown':'smoothUnknown.savg.mv3d'}
>>> data_out, \\
... error_out = bootstrap_tpb_stats(inputs_dict,
... n_volumes=500,
... box_size=4000,
... save_output=False,
... data_output_fname='data_N500_s4000.csv',
... err_output_fname='errors_N500_s4000.csv',
... output_avg=False)
"""
# Create descriptive loading bar
loading_bar = tqdm(inputs_dict.items(), desc='Loading skeleton file')
# Empty dataframe for holding data
overall_results = pd.DataFrame()
# Loop through the input dictionary, accessing each file
for k, v in loading_bar:
loading_bar.set_description('Loading skeleton file {}'.format(v))
loading_bar.update()
# Get data and some stats from the mv3d TPB network file
data, num_lines, num_points = read_mv3d(v)
# Infer appropriate bounding box from the data
bb = get_bb_from_data(data)
if output_avg:
data_cols = ["_totL", "_avgL", "_TPBdens"]
else:
data_cols = ["_totL", "_TPBdens"]
# Create temporary data frame to hold data for this key
res_df1 = pd.DataFrame(columns=[k + i for i in data_cols])
# Descriptive progress bar for subvolume sampling
# Loop through n_volumes different subvolumes
subvol_bar = tqdm(range(n_volumes), desc='Calculating on {}'.format(v))
for n in subvol_bar:
# Get a random position within the data, crop the TPB network,
# and try to split outlying large jumps (where the TPB path
# intersects a face of the subvolume)
box_start, box_end = get_random_subvolume(bb, box_size)
cropped_data = crop_tpb_data(data, box_start, box_end)
try:
split_cropped_data = split_paths_in_network(cropped_data)
except ValueError as e:
# If the outlier detection didn't work, just default to the
# regular cropped data
# print("Outlier detection failed for {0} "
# "in iteration {1}".format(v, n))
split_cropped_data = cropped_data
# Gather the data about this cropped network:
lengths = network_length(split_cropped_data)
total_length_nm = lengths[0]
total_length_um = total_length_nm / 1000
avg_length_nm = lengths[1].mean()
avg_length_um = avg_length_nm / 1000
bb_vol_nm = box_size**3
bb_vol_um = bb_vol_nm / (1000**3)
density_nm = total_length_nm / bb_vol_nm
density_um = total_length_um / bb_vol_um
if output_avg:
data_results = [total_length_um, avg_length_um, density_um]
else:
data_results = [total_length_um, density_um]
# Add each of the above values to a temporary dataframe (microns)
df1 = pd.DataFrame(dict(zip([k + i for i in data_cols],
data_results)),
index=[0])
# Append this set of values to the "inner" results list as
# another row. After this loop is finished, there should be
# n_volumes rows in res_df1
res_df1 = pd.concat([res_df1, df1], ignore_index=True)
# Build an "outer" results dataframe by adding the three columns
# from the inner results to the overall results dataframe
overall_results = pd.concat([overall_results, res_df1],
axis=1)
# Calculate actual errors and save data
data_df = overall_results
error_df = calculate_errors(overall_results, n_bootstrap)
if save_output:
data_df.to_csv(path_or_buf=data_output_fname)
print("Data output saved to {}".format(data_output_fname))
error_df.to_csv(path_or_buf=err_output_fname)
print("Error output saved to {}".format(err_output_fname))
return data_df, error_df
# noinspection PyUnusedLocal
[docs]def read_mv3d(filename):
"""
Get number of lines and points, as well as the 3d data contained withing
an MV3D network file
Parameters
----------
filename: str
Name of ``.mv3d`` file to open
Returns
-------
data: :py:class:`~numpy.ndarray`
(N x 4) numpy array containing the index, x, y, and z coordinates of
each point within the network (the last value, ``d`` is discarded)
num_lines: :class:`int`
number of lines contained in the network (read from line 2 of the file)
num_points: :class:`int`
number of points contained in the network (read from line 3 of the
file)
"""
lines = []
# Open file, discard first line, and return the remaining ones into a list
with open(filename, 'r') as f:
_ = f.readline()
lines = f.readlines()
# Remove newlines and tab characters from the text within `lines`
for i, l in enumerate(lines):
lines[i] = l.strip('\n').replace('\t', ' ')
# Get number of lines and points from the file
num_lines = int(lines[0].split(' ')[-1])
num_points = int(lines[1].split(' ')[-1])
# Remove blank lines and lines starting with '#' (header lines)
lines1 = [x for x in lines if x != '' and not x.startswith('#')]
# Split each string into many strings so it can be read into a numpy array
lines2 = [x.split()[:-1] for x in lines1]
# Create numpy data array
data = np.array(lines2,
dtype='float32')
return data, num_lines, num_points
[docs]def get_random_subvolume(bb, size):
"""
Given a particular bounding box and subvolume size, return
two lists with x, y, z coordinates of lower and upper corners of a random
subvolume withing the bounding box. Returned subvolume will be completely
enclosed within bb (i.e. the maximum position for the lower bound of the
cube x value will be bb[0] - size).
Parameters
----------
bb: tuple
tuple of length two, each term should be an :term:`iterable` of length
three with the minimum (position 0) and maximum (position 1) bounding
box coordinate in each dimension from which to take a random volume;
returned volume will be in the range
[min_x, max_x], [min_y, max_y], and [min_z, max_z]
size: :ref:`number <numbers>`
size of cube to return (single edge length, so total volume enclosed
will be size**3)
Returns
-------
box_start: :class:`list`
lower bound corner of the subvolume box (x, y, and z)
box_end: :class:`list`
upper bound corner of the subvolume box (x, y, and z
"""
import random
bb_min, bb_max = bb
box_start = [random.uniform(bb_min[i], x - size) for i, x in enumerate(
bb_max)]
box_end = [x + size for x in box_start]
return box_start, box_end
[docs]def split_paths_in_network(data,
threshold=0.2):
"""
Given an array of data, try to detect the paths that are on the edge of
the volume and split the single path (containing a sudden long jump)
into multiple paths.
Parameters
----------
data: :py:class:`~numpy.ndarray`
array of data, like that returned by
:py:func:`~FIBbootstrap.tpb.read_mv3d` or
:py:func:`~FIBbootstrap.tpb.crop_tpb_data`
threshold: float
weighting threshold in the range [0, 1] to help determine what is an
outlier. All the data in each path is fit by a robust linear model,
so any length values that are significantly different should have
a weight << 1.0. Set this value higher to find more outliers, and thus
split more paths. Set it lower to be more conservative.
Returns
-------
split_data: :py:class:`~numpy.ndarray`
copy of the original data, but with paths split at points that caused
particularly large jumps from point to point
"""
df = pd.DataFrame(data, columns=['id', 'x', 'y', 'z'])
next_id = df['id'].max() + 1
for idx, g_df in df.groupby('id'):
# find pairwise distance and calculate an array of distances between
# subsequent nodes in the path (`step_lengths`)
pdist_vec = pdist(g_df[['x', 'y', 'z']], metric='euclidean')
dist_mat = squareform(pdist_vec)
step_lengths = np.diagonal(dist_mat, offset=1)
# create robust linear model to estimate function describing data
rlm_model = sm.RLM(step_lengths,
range(len(step_lengths)),
M=sm.robust.norms.HuberT())
rlm_results = rlm_model.fit()
# outliers are where the model does not describe the data well
# (weight is below some threshold)
jumps = np.where(np.array(rlm_results.weights) < threshold)[0]
# if there were large jumps in the data:
if len(jumps) > 0:
# split the current grouped dataframe at the jump location
split_df_list = np.split(g_df, [x + 1 for x in jumps])
# change the id of split dataframes (except for the first one) and
# increment the next_id counter
for d in split_df_list[1:]:
d['id'] = int(next_id)
next_id += 1
# remove the current (unsplit) segment from the total dataframe
df = df.query('id != {}'.format(idx))
# append the split segments to the end of the overall dataframe
for d in split_df_list:
df = df.append(d)
split_data = df.as_matrix()
return split_data
[docs]def write_mv3d(fname, data, d=0, overwrite=True):
"""
Output path data to an mv3d file, which can be read by Avizo (and other
software)
Parameters
----------
fname: str
Filename to which to write; will be overwritten if it exists (by
default)
data: :py:class:`~numpy.ndarray`
array containing network (spatial graph) data in the same format as
output by :py:func:`~FIBbootstrap.tpb.read_mv3d` or
:py:func:`~FIBbootstrap.tpb.crop_tpb_data`
d: :py:class:`~numpy.ndarray`
the value to be written in the 'thickness' column of the mv3d file. Can
be used to tag files with a scalar value. If a single number is
given instead of an array, that number will be used for every point.
If a numpy array (the same length as the data), the values will be
specific to each point. Standard pandas/numpy broadcasting rules apply
overwrite: bool
switch to control whether an existing file will be clobbered if it
already exists
"""
if os.path.exists(fname):
if overwrite:
print('Output file exists; overwriting...')
pass
else:
raise IOError('Output file exists;'
'set overwrite=True to clobber the file')
# Put data into a pandas dataframe for processing
df = pd.DataFrame(data, columns=['id', 'x', 'y', 'z'])
# assign 'd' column for dataframe
df['d'] = d
# calculate the number of lines and points (to be written in the mv3d
# header)
num_lines = len(df.groupby('id'))
num_points = len(df)
# Disable some pandas warnings that are not relevant
original_warn_setting = pd.options.mode.chained_assignment
pd.options.mode.chained_assignment = None
# Write mv3d header into fname
with open(fname, 'w') as f:
f.write('# MicroVisu3D file\n')
f.write('# Number of lines {}\n'.format(num_lines))
f.write('# Number of points {}\n'.format(num_points))
f.write('# Number of inter. 0\n')
f.write('#\n')
f.write('# No x y z d\n')
f.write('#\n')
# Loop through the paths defined by 'id', reset their id number to an
# incrementing counter, and then write the data to the file, placing a
# space between each individual path within the network
for i, (idx, group_df) in enumerate(df.groupby('id')):
group_df['id'] = i
with open(fname, 'a') as f:
# write just the data to the file (no header or index info)
group_df.to_csv(f,
float_format='%.6f',
sep='\t',
header=False,
index=False)
f.write('\n')
pd.options.mode.chained_assignment = original_warn_setting
[docs]def crop_tpb_data(data, box_start, box_end):
"""
Crop TPB data to only contain points within the box defined by the corners
box_start and box_end
Parameters
----------
data: :py:class:`~numpy.ndarray`
(N x 4) numpy array with TPB data (as loaded by
:py:func:`~FIBbootstrap.tpb.read_mv3d`),
including the first column containing the index
box_start: list or :py:class:`~numpy.ndarray`
x, y, z coordinates of lower bound corner to crop inside of
box_end: list or :py:class:`~numpy.ndarray`
x, y, z coordinates of upper bound corner to crop inside of
Returns
-------
cropped_data: :py:class:`~numpy.ndarray`
copy of the original data, containing only the points inside of the
crop box
"""
# want to find the rows of the data array to keep:
# find where data is larger than the start, and less than the end
comp = np.logical_and(data[:, 1:] > box_start, data[:, 1:] < box_end)
# Using np.all with the axis=1 option tests just along rows
all_true = np.all(comp, axis=1)
# Find the index of the rows that
true_row_idx = np.where(all_true)[0]
# get only the data inside the box
cropped_data = data[true_row_idx, :]
return cropped_data
[docs]def path_length(path_df):
"""
Calculate the length along a path.
Parameters
----------
path_df: ~pandas.DataFrame
path_df should have 4 columns, 'id', 'x', 'y', 'z'. Using
:py:func:`scipy.spatial.distance.pdist`,
this method will calculate the sum of the distances between successive
rows of (x, y, z) coordinates in the dataframe
Returns
-------
length: float
total sum of path defined by successive points in ``path_df``
"""
# using just the x, y, z coordinates, calculate the condensed distance
# matrix for the series of points (pdist will calculate the distance
# between every point in the list, not just successive points)
pdist_vec = pdist(path_df[['x', 'y', 'z']], metric='euclidean')
# convert condensed distance matrix into a square distance matrix, which
# holds the distance between points i and j at dist_mat[i, j]
dist_mat = squareform(pdist_vec)
# since we are interested in the distance between successive points,
# the diagonal of dist_mat (offset by 1) holds all the individual lengths
# we are interested in. np.trace will sum these for us automatically
length = np.trace(dist_mat, offset=1)
return length
[docs]def network_length(data):
"""
Calculate the total length of a network that has been imported from
the mv3d format
Parameters
----------
data: ~numpy.ndarray
network data, in the format returned by
:py:func:`~FIBbootstrap.tpb.read_mv3d` or
:py:func:`~FIBbootstrap.tpb.crop_tpb_data`
Returns
-------
lengths: :class:`~numpy.ndarray`
lengths of each individual segment within the network
tot_length: :class:`float`
total length of network
"""
# Convert numpy array network data into a pandas Dataframe
df = pd.DataFrame(data, columns=['id', 'x', 'y', 'z'])
# Group the dataframe by segment id, so we can calculate individual lengths
groupby_id = df.groupby(df['id'])
# Apply the path_length function to each group, which will calculate the
# summed Euclidean distance between each point in the segment
lengths_series = groupby_id.apply(path_length)
# Convert the resulting pandas Series into a numpy array
lengths = lengths_series.as_matrix()
# Calculate the total length of all the segments
tot_length = lengths.sum()
return tot_length, lengths
[docs]def get_box_and_corners(box_start, box_end):
"""
Given two opposite corners of a 3D rectangle, find the necessary vectors
that will allow for plotting of the edges of the rectangle and points at
each of the corners using MayaVi
Parameters
----------
box_start: :term:`iterable`
lower bounding corner of the rectangle (i.e. [x, y, z])
box_end: :term:`iterable`
upper bounding corner of the rectangle (i.e. [x, y, z])
Returns
-------
x: :py:class:`~numpy.ndarray`
array of x positions for :py:func:`mayavi.mlab.plot3d` for plotting
edges of 3D rectangle
y: :py:class:`~numpy.ndarray`
array of y positions for :py:func:`~mayavi.mlab.plot3d` for plotting
edges of 3D rectangle
z: :py:class:`~numpy.ndarray`
array of z positions for :py:func:`~mayavi.mlab.plot3d` for plotting
edges of 3D rectangle
x_p: :py:class:`~numpy.ndarray`
array of x positions for :py:func:`~mayavi.mlab.plot3d` for plotting
corners of 3D rectangle
y_p: :py:class:`~numpy.ndarray`
array of y positions for :py:func:`~mayavi.mlab.plot3d` for plotting
corners of 3D rectangle
z_p: :py:class:`~numpy.ndarray`
array of z positions for :py:func:`~mayavi.mlab.plot3d` for plotting
corners of 3D rectangle
"""
c = np.array([box_start, box_end])
x = np.hstack((c[0, 0], c[1, 0], c[1, 0], c[0, 0], c[0, 0],
c[0, 0], c[1, 0], c[1, 0], c[0, 0], c[0, 0],
c[1, 0], c[1, 0], c[1, 0], c[1, 0], c[0, 0],
c[0, 0]))
y = np.hstack((c[0, 1], c[0, 1], c[1, 1], c[1, 1], c[0, 1],
c[0, 1], c[0, 1], c[1, 1], c[1, 1], c[0, 1],
c[0, 1], c[0, 1], c[1, 1], c[1, 1], c[1, 1],
c[1, 1]))
z = np.hstack((c[0, 2].repeat(5), c[1, 2].repeat(5),
c[1, 2], c[0, 2], c[0, 2], c[1, 2], c[1, 2],
c[0, 2]))
x_p = np.hstack((c[0, 0].repeat(4),
c[1, 0].repeat(4)))
y_p = np.hstack((c[0, 1], c[0, 1], c[1, 1], c[1, 1],
c[1, 1], c[0, 1], c[1, 1], c[0, 1]))
z_p = np.hstack((c[0, 2], c[1, 2], c[1, 2], c[0, 2],
c[1, 2], c[0, 2], c[0, 2], c[1, 2]))
return x, y, z, x_p, y_p, z_p
[docs]def get_bb_lines(bb):
"""
Given a bounding box, find the necessary vectors that will allow for
plotting of the edges of the bounding box using MayaVi
Parameters
----------
bb: :class:`tuple` of length 2
tuple of length two, each term should be an :term:`iterable` of length
three with the minimum (position 0) and maximum (position 1) bounding
box coordinate in each dimension from which to take a random volume;
returned volume will be in the range
Returns
-------
x: :py:class:`~numpy.ndarray`
array of x positions for :py:func:`mayavi.mlab.plot3d` for plotting
edges of 3D rectangle
y: :py:class:`~numpy.ndarray`
array of y positions for :py:func:`~mayavi.mlab.plot3d` for plotting
edges of 3D rectangle
z: :py:class:`~numpy.ndarray`
array of z positions for :py:func:`~mayavi.mlab.plot3d` for plotting
edges of 3D rectangle
"""
box_start = list(x for x in bb[0])
box_end = list(x for x in bb[1])
c = np.array([box_start, box_end])
x = np.hstack((c[0, 0], c[1, 0], c[1, 0], c[0, 0], c[0, 0],
c[0, 0], c[1, 0], c[1, 0], c[0, 0], c[0, 0],
c[1, 0], c[1, 0], c[1, 0], c[1, 0], c[0, 0],
c[0, 0]))
y = np.hstack((c[0, 1], c[0, 1], c[1, 1], c[1, 1], c[0, 1],
c[0, 1], c[0, 1], c[1, 1], c[1, 1], c[0, 1],
c[0, 1], c[0, 1], c[1, 1], c[1, 1], c[1, 1],
c[1, 1]))
z = np.hstack((c[0, 2].repeat(5), c[1, 2].repeat(5),
c[1, 2], c[0, 2], c[0, 2], c[1, 2], c[1, 2],
c[0, 2]))
return x, y, z
[docs]def get_bb_from_data(data_list):
"""
Infer from the the existing data the extents of the bounding box
Parameters
----------
data_list: list
list or array of numpy arrays with data (in the format returned by
:py:func:`~FIBbootstrap.tpb.read_mv3d` or
:py:func:`~FIBbootstrap.tpb.crop_tpb_data`.
Expects 4 columns representing ['id', 'x', 'y', 'z'].
Returns
-------
min_bb: :class:`tuple`
list of x, y, z values containing the smallest coordinates in each
dimension found in the data_list
max_bb: :class:`tuple`
list of x, y, z values containing the largest coordinates in each
dimension found in the data_list
"""
max_bb = tuple(np.vstack(data_list).max(axis=0)[1:])
min_bb = tuple(np.vstack(data_list).min(axis=0)[1:])
return min_bb, max_bb
[docs]def bb_volume(bb):
"""
Return the volume enclosed by a bounding box.
Parameters
----------
bb: :class:`tuple` of length 2
tuple of length two, each term should be an :term:`iterable` of length
three with the minimum (position 0) and maximum (position 1) bounding
box coordinate in each dimension
Returns
-------
volume: float
volume enclosed by the bounding box
"""
min_bb, max_bb = bb
volume = 1
for i in range(3):
volume *= (max_bb[i] - min_bb[i])
return volume
# noinspection PyUnusedLocal
[docs]def animate_cropped_data(data_a,
data_i,
data_u,
size=4000,
subsample=1,
bb=None):
"""
Plot a simple animation of the total TPB path network, as well as the
network contained within a randomly cropped volume
Parameters
----------
data_a: :py:class:`~numpy.ndarray`
network data, in the format returned by
:py:func:`~FIBbootstrap.tpb.read_mv3d` or
:py:func:`~FIBbootstrap.tpb.crop_tpb_data`;
will be plotted in green (active)
data_i: :py:class:`~numpy.ndarray`
network data, in the format returned by
:py:func:`~FIBbootstrap.tpb.read_mv3d` or
:py:func:`~FIBbootstrap.tpb.crop_tpb_data`;
will be plotted in red (inactive)
data_u: :py:class:`~numpy.ndarray`
network data, in the format returned by
:py:func:`~FIBbootstrap.tpb.read_mv3d` or
:py:func:`~FIBbootstrap.tpb.crop_tpb_data`;
will be plotted in yellow (unknown)
size: :ref:`number <numbers>`
size of cube to return (single edge length, so total volume enclosed
will be size**3)
subsample: int
factor by which to subsample the total tpb network. If > 1, it will
speed up the plotting (may be helpful on lower-powered CPUs)
bb: :data:`None` or :class:`tuple` of length 2
tuple of length two, each term should be an :term:`iterable` of length
three with the minimum (position 0) and maximum (position 1) bounding
box coordinate in each dimension from which to take a random volume;
returned volume will be in the range
[min_x, max_x], [min_y, max_y], and [min_z, max_z]
If None, the bounding box will be inferred from the supplied data
"""
from mayavi import mlab
green = (0.3, 0.7, 0.4)
yellow = (0.9, 0.9, 0.5)
red = (0.9, 0.2, 0.2)
orange = (255.0 / 255, 165. / 255, 0)
white = (1, 1, 1)
if bb is None:
bb = get_bb_from_data([data_a, data_i, data_u])
min_bb, max_bb = bb
else:
min_bb, max_bb = bb
x_a, y_a, z_a = data_a[::subsample, 1:].T
x_i, y_i, z_i = data_i[::subsample, 1:].T
x_u, y_u, z_u = data_u[::subsample, 1:].T
f1 = mlab.figure(size=(500, 500))
f2 = mlab.figure(size=(500, 500))
# Plotting the subvolume box on both figures
options = {'tube_radius': 50, 'tube_sides': 100}
point_options = {'scale_factor': 200, 'color': (0, 0, 0)}
box_start, box_end = get_random_subvolume(bb, size)
x, y, z, x_p, y_p, z_p = get_box_and_corners(box_start, box_end)
m_lines1 = mlab.plot3d(x, y, z, figure=f1, **options)
m_points1 = mlab.points3d(x_p, y_p, z_p, figure=f1, **point_options)
m_lines2 = mlab.plot3d(x, y, z, figure=f2, **options)
m_points2 = mlab.points3d(x_p, y_p, z_p, figure=f2, **point_options)
# Plotting the tpb points
tpb_options = {'mode': 'point',
'scale_factor': 50,
'opacity': 0.5}
x_tot_a, y_tot_a, z_tot_a = data_a[::subsample, 1:].T
x_tot_i, y_tot_i, z_tot_i = data_i[::subsample, 1:].T
x_tot_u, y_tot_u, z_tot_u = data_u[::subsample, 1:].T
# plot total TPB network in figure 1
m_tot_a = mlab.points3d(x_tot_a, y_tot_a, z_tot_a,
color=green, figure=f1, **tpb_options)
m_tot_i = mlab.points3d(x_tot_i, y_tot_i, z_tot_i,
color=red, figure=f1, **tpb_options)
m_tot_u = mlab.points3d(x_tot_u, y_tot_u, z_tot_u,
color=yellow, figure=f1, **tpb_options)
# plot just the subsampled volume in figure 2
m_a = mlab.points3d(x_a, y_a, z_a,
color=green, figure=f2, **tpb_options)
m_i = mlab.points3d(x_i, y_i, z_i,
color=red, figure=f2, **tpb_options)
m_u = mlab.points3d(x_u, y_u, z_u,
color=yellow, figure=f2, **tpb_options)
# plot overall bounding box
x_bb, y_bb, z_bb = get_bb_lines(bb)
bb_lines1 = mlab.plot3d(x_bb, y_bb, z_bb,
color=orange,
tube_radius=20,
tube_sides=100,
figure=f1)
bb_lines2 = mlab.plot3d(x_bb, y_bb, z_bb,
color=orange,
tube_radius=20,
tube_sides=100,
figure=f2)
# noinspection PyUnusedLocal
@mlab.animate(delay=2000)
def anim():
f = mlab.gcf()
while True:
box_inner_start, box_inner_end = get_random_subvolume(bb, size)
# box_start, box_end = random_box(bb, size)
cropped_a = crop_tpb_data(data_a, box_inner_start, box_inner_end)
cropped_i = crop_tpb_data(data_i, box_inner_start, box_inner_end)
cropped_u = crop_tpb_data(data_u, box_inner_start, box_inner_end)
x_inner_a, y_inner_a, z_inner_a = cropped_a[:, 1:].T
x_inner_i, y_inner_i, z_inner_i = cropped_i[:, 1:].T
x_inner_u, y_inner_u, z_inner_u = cropped_u[:, 1:].T
x_inner, y_inner, z_inner, x_inner_p, y_inner_p, z_inner_p = \
get_box_and_corners(box_inner_start, box_inner_end)
m_lines1.mlab_source.reset(x=x_inner, y=y_inner, z=z_inner)
m_points1.mlab_source.reset(x=x_inner_p, y=y_inner_p, z=z_inner_p)
m_lines2.mlab_source.reset(x=x_inner, y=y_inner, z=z_inner)
m_points2.mlab_source.reset(x=x_inner_p, y=y_inner_p, z=z_inner_p)
m_a.mlab_source.reset(x=x_inner_a, y=y_inner_a, z=z_inner_a)
m_i.mlab_source.reset(x=x_inner_i, y=y_inner_i, z=z_inner_i)
m_u.mlab_source.reset(x=x_inner_u, y=y_inner_u, z=z_inner_u)
mlab.sync_camera(f1, f2)
yield
anim()
mlab.show()