genepy3d.obj package

Module contents

Submodules

genepy3d.obj.curves module

Methods for working with Curve objects.

Classes:

Curve(coors)	Curve in 3D.
SimuIntrinsic([seg_lbl, npoints, sigma, ...])	Generate curve in 3D with random intrinsic dimension segments.
RadiusScaleTable(crv, r[, mo, kerlen])	Compute scales table of a curve from given radius of curvature.

Functions:

decompose_intrinsicdim_sequential(p, ...)	Decompose the curve into intrisic segments of 1D line, 2D plane and 3D.
align(target, ref)	Align a target curve to match a reference curve.

class genepy3d.obj.curves.Curve(coors)

Bases: object

Curve in 3D. Also work with 2D curve.

Please check the documentation of each function to see if it supports 2D curve. We assumme that curve is a parametric function g(t) = (x(t),y(t),z(t)). If 2D points are passed then we consider only X and Y coordinates and set Z = 0.

The Curve is assumed open and with no duplicated positions.

coors

list of 3d points.

• if array is given, then each column is the x, y, z coordinates of a specific point.

• if tuple is given, then the items are the three arrays of x, y and z coordinates.

Type:

array | tuple

Examples

import numpy as np
from genepy3d.obj.curves import Curve

# Curve from array input
coors = np.array([[1,2,3],[1,2,3],[1,2,3]])
crv = Curve(coors)

# Curve from tuple input
coors = ([1,1,1],[2,2,2],[3,3,3])
crv = Curve(coors)

Attributes:

size	Number of points on the curve.

Methods:

from_csv(filepath[, column_names, scales, args])	Constructor from csv file.
from_text(filepath)	Read from text-liked file (.txt, .xyz).
compute_derivative(deg[, dt])	Derivative using np.gradient.
compute_length()	Return the length of curve.
compute_direction_change()	Compute the angle (directional change) at a give point and its orientation (CW or CCW).
compute_angles()	Angle between two vectors from two nearby points of a given point.
compute_angles_vector([u])	Angles between vector at a node on the curve and a vector u.
compute_curvature()	Curvatures at every point on the curve.
compute_torsion()	Torsions at every point on the curve.
compute_wiggliness()	Wiggliness at every point on the curve.
compute_tortuosity()	Tortuosity of every point on the curve.
transform([phi, theta, psi, dx, dy, dz, zx, ...])	Apply the same rigid transformation to every point on the curve.
convolve_gaussian(sigma[, mo, kerlen])	Get curve at a specific scale sigma by Gaussian convolution.
resample([unit_length, npoints, ...])	Resample the curve using spline interpolation.
estimate_noise_variance()	Estimate noise variance across each axis x, y and z of the curve.
denoise([sigma_step, max_iter, return_sigma])	Denoise a curve assuming a Gaussian noise.
scale_space(scale_range[, features, ...])	Compute features of a curve across scales.
decompose_intrinsicdim(sig_c[, delta_sig, ...])	Decompose curve into multiscale hierarchichal intrinsic segments at a given sigma sig_c.
plot_decomposed_table(ax[, show_selection, ...])	Display decomposed intrinsic table.
plot_intrinsicdim(ax[, projection, scales, ...])	Display curve superimposed by estimated intrinsic dimensions.
compute_local_3d_scale_sigma(sig_lst[, ...])	Computing local 3d scale of curve from the list of sigma i.e stdev of Gaussian function.
compute_local_3d_scale_radius(r_lst[, ...])	Computing local 3d scale of curve from the list of radius of curvatures.
compute_local_3d_scale(r_lst[, eps_seg_len, ...])	Computing local 3d scale of curve.
main_turns(sig_lst[, search_step, ...])	Compute the main turns of curve.
main_turns_old([nbscales, search_step, ...])	Compute the main turns of curve.
plot_ridge_map(ax[, show_scales])	Display ridge finding from main turns computation.
coors_to_pixel(dx, dy, dz)	Convert coordinates to pixel coordinates.
to_points()	Convert Curve object to Points object.
plot(ax[, projection, scales, show_root, ...])	Plot curve using matplotlib.

property size

Number of points on the curve.

classmethod from_csv(filepath, column_names=['x', 'y', 'z'], scales=(1.0, 1.0, 1.0), args={})

Constructor from csv file.

Parameters:

• filepath (str) – csv file.

• column_names (list of str) – coordinates columns in csv file. Default column names to look for is “x”, “y” and “z”.

• scales (tuple of float) – scales in x, y and z.

• args (dict) – overried parameters of pandas.read_csv().

Returns:

a Curve object.

Examples

from genepy3d.obj.curves import Curve
filepath = 'path/to/csv/file'
crv = Curve.from_csv(filepath)

classmethod from_text(filepath)

Read from text-liked file (.txt, .xyz).

Parameters:

filepath (str) – text file.

Returns:

a Curve object.

Notes

Assuming the first three columns correspond to x, y, z coordinates.

Examples

from genepy3d.obj.curves import Curve
filepath = 'path/to/text/file'
crv = Curve.from_text(filepath)

compute_derivative(deg, dt=1)

Derivative using np.gradient.

Also support 2D. If in 2D, then return 0 for the derivative of Z.

Parameters:

• deg (int) – derivative degree.

• dt (float) – delta t.

Returns:

array of float where each row is the derivative in x, y and z at a specific point.

compute_length()

Return the length of curve.

It is the sum of L2 distances of curve segments.

Also support 2D.

Returns:

a float.

compute_direction_change()

Compute the angle (directional change) at a give point and its orientation (CW or CCW).

Assume three points A => B => C. The direction change at B is the angle between two vectors (A->B) and (B->C). The sign of angle indicates CW (negative) or CCW (positive) rotation (only work in 2D).

Support 2D: YES.

Returns:

Array of angles.

Notes

• I compute the angle using this suggestion: https://stackoverflow.com/questions/14066933/direct-way-of-computing-the-clockwise-angle-between-two-vectors

• The CW/CCW is only worked under 2D case.

Examples

import numpy as np
from genepy3d.obj import curves

# 3D curve from a helix
t = np.arange(50)
a = 1.
b = 1.
x = a * np.cos(t/5)
y = a * np.sin(t/5)
z = b * t
crv = curves.Curve((x,y,z))

# Angles w.r.t. x-vector
crv.compute_direction_change()

compute_angles()

Angle between two vectors from two nearby points of a given point.

The first and last points are not counted.

Assumming three consecutive nodes A, B, C, the angle at node B is the angle between two vector (B=>A) and (B=>C).

Support 2D: YES.

Returns:

Array of angles.

Examples

import numpy as np
from genepy3d.obj import curves

# 3D curve from a helix
t = np.arange(50)
a = 1.
b = 1.
x = a * np.cos(t/5)
y = a * np.sin(t/5)
z = b * t
crv = curves.Curve((x,y,z))

# Angles w.r.t. x-vector
crv.compute_angles()

compute_angles_vector(u=None)

Angles between vector at a node on the curve and a vector u.

Assuming two consecutive nodes A, B, the angle at the node A is the angle between vector (A=>B) and vector u.

Parameters:

u (array (float)) – vector used to compute angles. If None, then return angles with three x, y and z axes.

Returns:

Array of angles.

Examples

import numpy as np
from genepy3d.obj import curves

# 3D curve from a helix
t = np.arange(50)
a = 1.
b = 1.
x = a * np.cos(t/5)
y = a * np.sin(t/5)
z = b * t
crv = curves.Curve((x,y,z))

# Angles w.r.t. x-vector
crv.compute_angles_vector(np.array([1.,0.,0.]))

# Angles w.r.t. three axes x, y and z
# Return array with 3 columns for angles in x, y and z respectively
crv.compute_angles_vector()

compute_curvature()

Curvatures at every point on the curve.

Support 2D: YES.

Returns:

array of curvatures.

Notes

Due to the side effect when computing the derivative, the curvature can be not very accurate at the extremal points of the curve. The number of points affected at each of the two extremes is 2.

compute_torsion()

Torsions at every point on the curve.

Returns:

array of torsions.

Notes

Due to the side effect when computing the derivative, the torsion can be not very accurate at the extremal points of the curve. The number of points affected at each of the two extremes is 3.

compute_wiggliness()

Wiggliness at every point on the curve.

Also support 2D.

Wiggliness = Length / Distance

Returns:

array of wiggliness.

compute_tortuosity()

Tortuosity of every point on the curve.

This function is another name of compute_wiggliness().

Returns:

array of float.

transform(phi=0, theta=0, psi=0, dx=0, dy=0, dz=0, zx=1, zy=1, zz=1)

Apply the same rigid transformation to every point on the curve.

Parameters:

• phi (float) – rotation in x in radian.

• theta (float) – rotation in y in radian.

• psi (float) – rotation in z in radian.

• dx (float) – translation in x.

• dy (float) – translation in y.

• dz (float) – translation in z.

• zx (float) – zoom in x.

• zy (float) – zoom in y.

• zz (float) – zoom in z.

Returns:

a Curve object.

Examples

import numpy as np
from genepy3d.obj import curves
import matplotlib.pyplot as plt

# 3D curve from a helix
t = np.arange(50)
a = 1.
b = 1.
x = a * np.cos(t/5)
y = a * np.sin(t/5)
z = b * t
crv = curves.Curve((x,y,z))

# Plot the curve
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
crv.plot(ax);

# Rotate by 90 degree around x
# Then, shift by 10 in y
crv_transformed = crv.transform(phi=np.pi/2.,dy=10)

# Plot the transformed curve
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
crv_transformed.plot(ax);

convolve_gaussian(sigma, mo='nearest', kerlen=4.0)

Get curve at a specific scale sigma by Gaussian convolution.

Also support 2D.

Parameters:

• sigma (float or list of floats) – scale.

• mo (str) – mode used to treat the border effect (see gaussian_filter1d() in scipy for detail).

• kerlen (float) – specify the length of Gaussian kernel (see gaussian_filter1d() in scipy for detail).

Returns:

Curve object.

Examples

import numpy as np
from genepy3d.obj import curves
import matplotlib.pyplot as plt

# 3D curve from a helix
t = np.arange(50)
a = 1.
b = 1.
x = a * np.cos(t/5)
y = a * np.sin(t/5)
z = b * t
crv = curves.Curve((x,y,z))

# Plot the curve
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
crv.plot(ax);

sigma = 5
crv_gauss = crv.convolve_gaussian(sigma)

# Plot the smoothed curve
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
crv_gauss.plot(ax);

resample(unit_length=None, npoints=None, spline_order=1, return_interp_param=False)

Resample the curve using spline interpolation.

Support 2D: YES.

Parameters:

• unit_length (float) – sampling length. If None, then the curve is resampled from the input npoints.

• npoints (int) – number of resampled points. If None, then the number of resampled points is equal to the current number of points on the curve (i.e. before resampling).

• spline_order (uint) – if it equal 1, then the linear interpolation is used. Otherwise, spline interpolation is used. Please see interpolate.splprep() in scipy for more detail.

• return_interp_param (bool) – if True then return interpolation parameters.

Returns:

Curve object.

Notes

• Only choose to resample the curve by either the unit_length or npoints.

• The resampling can be failed in some cases due to the spline interpolation constraints. Please see interpolate.splprep() in scipy for more detail.

Examples

import numpy as np
from genepy3d.obj import curves
import matplotlib.pyplot as plt

# 3D curve from a helix
t = np.arange(20)
a = 1.
b = 1.
x = a * np.cos(t/5)
y = a * np.sin(t/5)
z = b * t
crv = curves.Curve((x,y,z))

# Plot the curve
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
crv.plot(ax,point_args={"c":"r"}); # set point_args to display the points on the curve

# Resample by 50 points
crv_resampled = crv.resample(npoints=50)

# Plot the resampled curve
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
crv_resampled.plot(ax,point_args={"c":"r"});

estimate_noise_variance()

Estimate noise variance across each axis x, y and z of the curve. Assuming noise is normal distribution.

Returns:

xsig, ysig and zsig the three estimated sigma for x, y and z.

References

[1]

A posteriori noise estimation in variable data sets - With applications to spectra and light curves. S. Czesla, T. Molle and J. H. M. M. Schmitt. A&A, 609 (2018) A39. DOI: 10.1051/0004-6361/201730618.

denoise(sigma_step=0.25, max_iter=None, return_sigma=False)

Denoise a curve assuming a Gaussian noise.

Work for both 2D and 3D.

Noise variances across x, y and z axes of the curve are estimated using the beta-sigma algorithm implemented in PyAstronomy [1]. The curve is then denoised by Gaussian convolution of an initial sigma = sigma_step. We repeat the denoising process by increasing sigma_step until the stdev between the denoised curve and the noisy curve approaches the estimated noise variances. The process will be also stopped if the number of iteration >= max_iter.

Parameters:

• sigma_step (float) – initial sigma of Gaussian used to smooth the curve.

• max_iter (int) – maximal number of iteration. If None, then max_iter = number of points on the curve.

• return_sigma (bool) – if True, then return the estimated sigmas in x, y and z.

Returns:

a denoised Curve object and list of estimated sigmas if return_sigma is True. If failed, then return None.

References

[1]

A posteriori noise estimation in variable data sets - With applications to spectra and light curves. S. Czesla, T. Molle and J. H. M. M. Schmitt. A&A, 609 (2018) A39. DOI: 10.1051/0004-6361/201730618.

Notes

• We assume noise follows the normal distribution. The denoising can be ineffecient if noise follows another distribution.

• The denoise() cannot guarantee an optimal noise removal, especially when the noise level is high.

Examples

import numpy as np
from genepy3d.obj import curves
import matplotlib.pyplot as plt

# 3D curve from a helix
n = 200 # number of points
t = np.arange(n)
a = 1.
b = 1.
sigma = 0.2 # sigma of noise to add into the curve
x = a * np.cos(t/5) + np.random.normal(0,sigma,n)
y = a * np.sin(t/5) + np.random.normal(0,sigma,n)
z = b * t + np.random.normal(0,sigma,n)
crv = curves.Curve((x,y,z))

# Plot the noisy curve
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
crv.plot(ax,point_args={"c":"r"}); # set point_args to display the points on the curve

# Denoised the curve
crv_denoised, sigma_est = crv.denoise(return_sigma=True)
print(sigma_est)

# Plot the denoised curve
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
crv_denoised.plot(ax,point_args={"c":"r"});

scale_space(scale_range, features={'curvature', 'line', 'planeline', 'ridge', 'torsion', 'valley'}, eps_kappa=0.01, eps_tau=0.01, eps_seg=10, mo='nearest', kerlen=4.0)

Compute features of a curve across scales.

We compute e.g. the curvatures, torsions, etc of a curve convolved by Gaussian of a given sigma. The process repeats over increasing sigma and produce finally a feature-scale space of the curve.

Parameters:

• scale_range (array of float) – range of scales (sigma of Gaussian).

• features (dic) – list of features (detail see below).

• eps_kappa (float) – curvature threshold. (for planeline feature)

• eps_tau (float) – torsion threshold. (for line and planeline features)

• eps_seg (float) – segment length threshold (in number of points, for line and planeline features).

Returns:

dictionary whose each item is a scale space matrix of a specific feature, the rows of matrix are the scale range and columns are the curve indices.

Note

• We support following features

  • curvature: curvature scale space.

  • torsion: torsion scale space.

  • ridge: local maxima of curvature

  • valley: local minima of curvature

  • planeline: plane+line scale space.

  • line: line scale space.

• The line and planeline scale space were used to compute the local 3d scale of the curve. See compute_local_3d_scale_sigma() for detail.

• The curvature scale space was used to compute the main turns of the curve. See main_turns() for detail.

Examples

import numpy as np
from genepy3d.obj import curves
import matplotlib.pyplot as plt

# 3D curve from a helix
t = np.arange(100)
a = 1.
b = 1.
x = a * np.cos(t/5)
y = a * np.sin(t/5)
z = b * t
crv = curves.Curve((x,y,z))

# Plot the curve
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
crv.plot(ax);

scale_range = np.arange(0,50.,1) # sigma ranges
features = crv.scale_space(scale_range,features={'curvature'}) # curvature scale space
kappa_space = features['curvature']

# Plot the curvature scale space
fig = plt.figure()
ax = fig.add_subplot(111)
pl = ax.imshow(kappa_space)
ax.invert_yaxis();
fig.colorbar(pl)

decompose_intrinsicdim(sig_c, delta_sig=None, eps_seg_len=None, eps_crv_len=None, sig_step=1, eps_kappa=0.01, eps_tau=0.01)

Decompose curve into multiscale hierarchichal intrinsic segments at a given sigma sig_c.

The function decomposes the curve into portions of 1D lines, 2D planes or 3D at a given scale sig_c. The decomposition is hierarchichal, i.e. 1D lines can be detected from a 2D plane portion. The detail of algorithm can be found [1].

Parameters:

• sig_c (float) – central scale of the searching region.

• delta_sig (float) – width of the searching interval.

• sig_step (int) – scale step within the searching region.

• eps_kappa (float) – curvature threshold.

• eps_tau (float) – torsion threshold.

• eps_seg_len (float) – segment length threshold.

• eps_crv_len (float) – curve length threshold.

Returns:

list of segment indices specifying line, plane+line.

References

[1]

Phan MS, Matho K, Beaurepaire E, Livet J, Chessel A. nAdder: A scale-space approach for the 3D analysis of neuronal traces. 2022. PLOS Computational Biology 18(7): e1010211. DOI: 10.1371/journal.pcbi.1010211

Notes

The decomposition can be visualized using plot_decomposed_table() and plot_intrinsicdim().

plot_decomposed_table(ax, show_selection=True, show_scales=True, aspect=3, xdiv=50, ydiv=10)

Display decomposed intrinsic table.

It is only used after running decompose_intrinsicdim().

Parameters:

• ax – plot axis.

• show_selection (bool) – if True, then display selection interval.

• show_scales (bool) – if True, display scale range in y axis.

• aspect (int) – aspect ratio between x and y axes.

• xdiv (int) – to calculate xticks.

• ydiv (int) – to calculate yticks.

plot_intrinsicdim(ax, projection='3d', scales=(1.0, 1.0, 1.0), root_args={'c': 'blue', 's': 100}, dim1d_args={'alpha': 1.0, 'c': 'magenta', 'lw': 1}, dim2d_args={'alpha': 0.7, 'c': 'yellow', 'lw': 3}, dim3d_args={'alpha': 0.7, 'c': 'green', 'lw': 3}, plane_color='yellow', overrided_curve=None)

Display curve superimposed by estimated intrinsic dimensions.

It is only used after running decompose_intrinsicdim().

Parameters:

• ax – plot axis.

• projection (str) – support 3d, xy, xz, yz modes.

• scales (tuple of float) – specify x, y, z scales.

• overrided_curve (Curve) – superimpose estimated intrinsic dimensions on the overrided_curve. if None, then a scaled curve computed from decompose_intrinsicdim() is used as default.

compute_local_3d_scale_sigma(sig_lst, delta_sig=0, sig_step=1, eps_seg_len=None, eps_crv_len=None, eps_kappa=0.01, eps_tau=0.01, return_dim_results=False)

Computing local 3d scale of curve from the list of sigma i.e stdev of Gaussian function.

Local 3D scale: scale at which point transforms from 3D to 2D/1D. Detail of the algorithm can be found in [1].

Parameters:

• sig_lst (list (float)) – list of sigma values.

• delta_sig (float) – width of the searching interval.

• sig_step (int|float) – scale step within the searching region.

• eps_seg_len (float) – segment length threshold.

• eps_crv_len (float) – curve length threshold.

• eps_kappa (float) – curvature threshold.

• eps_tau (float) – torsion threshold.

• return_dim_results (bool) – if True, then return intrinsic dim estimation of r_lst

Returns:

list of local 3d scales for every points on curve.

References

[1]

Phan MS, Matho K, Beaurepaire E, Livet J, Chessel A. nAdder: A scale-space approach for the 3D analysis of neuronal traces. 2022. PLOS Computational Biology 18(7): e1010211. DOI: 10.1371/journal.pcbi.1010211

compute_local_3d_scale_radius(r_lst, eps_seg_len=None, eps_crv_len=None, sig_step=1, eps_kappa=0.01, eps_tau=0.01, return_dim_results=False)

Computing local 3d scale of curve from the list of radius of curvatures.

Local 3D scale: scale at which point transforms from 3D to 2D/1D. Detail of the algorithm can be found in [1].

Parameters:

• r_lst (list (float)) – list of radii of curvatures. The radius must be > 0.

• eps_seg_len (float) – segment length threshold.

• eps_crv_len (float) – curve length threshold.

• sig_step (int) – scale step within the searching region.

• eps_kappa (float) – curvature threshold.

• eps_tau (float) – torsion threshold.

• return_dim_results (bool) – if True, then return intrinsic dim estimation of r_lst

Returns:

list of local 3d scales for every points on curve.

References

[1]

Phan MS, Matho K, Beaurepaire E, Livet J, Chessel A. nAdder: A scale-space approach for the 3D analysis of neuronal traces. 2022. PLOS Computational Biology 18(7): e1010211. DOI: 10.1371/journal.pcbi.1010211

compute_local_3d_scale(r_lst, eps_seg_len=None, eps_crv_len=None, sig_step=1, eps_kappa=0.01, eps_tau=0.01, return_dim_results=False)

Computing local 3d scale of curve.

Local 3D scale: scale at which point transforms from 3D to 2D/1D. This function is DEPRECATED. Use compute_local_3d_scale_radius() or compute_local_3d_scale_sigma().

Parameters:

• r_lst (list (float)) – list of radii of curvatures.

• eps_seg_len (float) – segment length threshold.

• eps_crv_len (float) – curve length threshold.

• sig_step (int) – scale step within the searching region.

• eps_kappa (float) – curvature threshold.

• eps_tau (float) – torsion threshold.

• return_dim_results (bool) – if True, then return intrinsic dim estimation of r_lst

Returns:

list of local 3d scales for every points on curve.

main_turns(sig_lst, search_step=5, eps_kappa=None, ridgelength_thr=0.1, angle_thr=20, min_dist=None)

Compute the main turns of curve.

Main turns are the positions on the curve showing important changes of orientation. We calculate the curvature scale space of the curve across a range of sigmas (Gaussian func). Then, track the change of curvatures of points on the curve within the curvature scale space. Finally, select the points as the main turns based on some criteria.

Parameters:

• sig_lst (list (float)) – list of sigma values used to compute curvature scale space.

• search_step (int) – number of neighboring points that are grouped to find the main turns.

• eps_kappa (float) – curvature threshold. The points whose curvatures are larger than the eps_kappa are taken as the candidates to find the main turns. If None, then it is set relatively to the maximal curvature of the curve.

• ridgelength_thr (float) – ridge length threshold (in % compared to the longest ridge). ridge is a track of a point in the curvature scale space. The point whose ridge length is larger than the ridgelength_thr is considered as main turn.

• angle_thr (float) – angle threshold in degree. Remove position whose angle is smaller than the angle_thr.

• min_dist (float) – minimal distance between two main turns. Make sure that the distance between two main turns is not smaller than min_dist.

Returns:

indices of main turns.

main_turns_old(nbscales=None, search_step=5, eps_kappa=None, ridgelength_thr=0.1, angle_thr=20, min_dist=None)

Compute the main turns of curve.

Main turns are the positions on the curve showing important changes of orientation. We calculate the curvature scale space of the curve across a range of sigmas (Gaussian func). Then, track the change of curvatures of points on the curve within the curvature scale space. Finally, select the points as the main turns based on some criteria.

The sigmas are defined from ``nbscales` input. sigmas = [0, 1, …, nbscales]

This function is DEPRECATED. Please refer to main_turns().

Parameters:

• nbscales (int) – number of scales used to compute the curvature scale space. If None, then it is set as halft of number of points of the curve.

• search_step (int) – number of neighboring points that are grouped to find the main turns.

• eps_kappa (float) – curvature threshold. The points whose curvatures are larger than the eps_kappa are taken as the candidates to find the main turns. If None, then it is set relatively to the maximal curvature of the curve.

• ridgelength_thr (float) – ridge length threshold (in % compared to the longest ridge). ridge is a track of a point in the curvature scale space. The point whose ridge length is larger than the ridgelength_thr is considered as main turn.

• angle_thr (float) – angle threshold in degree. Remove position whose angle is smaller than the angle_thr.

• min_dist (float) – minimal distance between two main turns. Make sure that the distance between two main turns is not smaller than min_dist.

Returns:

indices of main turns.

plot_ridge_map(ax, show_scales=True)

Display ridge finding from main turns computation.

This is only used after running main_turns().

Parameters:

• ax – plot axis.

• show_scales (bool) – if True, display scale range in y axis.

coors_to_pixel(dx, dy, dz)

Convert coordinates to pixel coordinates.

Parameters:

• dx (float) – pixel size in x.

• dy (float) – pixel size in y.

• dz (float) – pixel size in z.

Returns:

array of pixel coordinates.

to_points()

Convert Curve object to Points object.

Returns:

Points object.

plot(ax, projection='3d', scales=(1.0, 1.0, 1.0), show_root=True, root_args={'color': 'red'}, line_args={'c': 'blue', 'lw': 1}, point_args=None)

Plot curve using matplotlib.

Parameters:

• ax – matplotlib axis.

• projection (str) – support ‘3d’, ‘xy’, ‘xz’ and ‘yz’.

• scales (tuple) – scale in x, y and z.

• show_root (bool) – If True, then show the root position (the first point).

• root_args (dic) – Pass this to the scatter() in matplotlib.

• line_args (dic) – Pass this to the plot() in matplotlib.

• point_args (dic) – Pass this to the scatter() in matplotlib.

Examples

import numpy as np
from genepy3d.obj import curves
import matplotlib.pyplot as plt

# 3D curve from a helix
n = 200 # number of points
t = np.arange(n)
a = 1.
b = 1.
x = a * np.cos(t/5)
y = a * np.sin(t/5)
z = b * t
crv = curves.Curve((x,y,z))

# Plot the noisy curve
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
crv.plot(ax,point_args={"c":"r"}); # set point_args to display the points on the curve

class genepy3d.obj.curves.SimuIntrinsic(seg_lbl=[0, 1, 2], npoints=1000, sigma=1.0, random_state=None, params={}, max_seg=None, nb_seg=None, remove_zero_replica=True, spline_order=1)

Bases: object

Generate curve in 3D with random intrinsic dimension segments.

The process is simulated based on 2D and 3D Brownian motion. See active_brownian_2d() and active_brownian_3d() in genepy3d.util.geo module for the generation of Brownian motion.

seg_lbl

list of generated curve segments, 0 for 3d, 1 for 2d, 2 for 1d.

Type:

array of int

npoints

number of points on simulated curve.

Type:

int

sigma

for adding Gaussian noise.

Type:

float

params

simulation parameters.

Type:

array of dic | dic

random_state

seed number used to memorize the simulated curve.

Type:

int

max_seg

maximal number of segments. if this is set, then number of simulated intrinsic segments will be <= max_seg.

Type:

int

nb_seg

number of segments. if this is set, then number of simulated intrinsic segments will be fixed at nb_seg.

Type:

int

remove_zero_replica

if True, then replace the ‘0,0,..’ sequence into only one ‘0’.

Type:

bool

Notes

specify either seg_lbl, max_seg or nb_seg whose priority order is seg_lbl first, then max_seg and nb_seg.

Examples

from genepy3d.obj import curves
import matplotlib.pyplot as plt

# Create a simulator to generate curve consists of 3D=>1D=>3D=>2D=>1D=>3D segments.
simulator = curves.SimuIntrinsic(seg_lbl=[0,2,0,1,2,0]).gen()

print("Number of points:",simulator.npoints)
print("Line indices:",simulator.line_ids)
print("Plane indices:",simulator.plane_ids)

# Plot the simulated curve
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
simulator.plot(ax)

Methods:

gen()	Generate a 3D curve with specified intrinsic lines, planes and 3D.
add_noise(sigma)	Add Gaussian noise to curve.
get_curve([noisy])	Return curve with/without noise.
evaluate(pred_line_ids, pred_plane_ids[, metric])	Evaluate accuracies of intrinsic dimension estimation.
plot(ax[, projection, noisy, scales, ...])	Plot simulated curve with its intrinsic segments.

gen()

Generate a 3D curve with specified intrinsic lines, planes and 3D.

Returns:

itself.

Notes

generated curve for noiseless case can contain duplicata points.

add_noise(sigma)

Add Gaussian noise to curve.

Parameters:

sigma (float) – std of Gaussian kernel.

Returns:

itself.

get_curve(noisy=True)

Return curve with/without noise.

Parameters:

noisy (bool) – if True, then return noisy curve, and False otherwise.

Returns:

Curve object.

evaluate(pred_line_ids, pred_plane_ids, metric='f1')

Evaluate accuracies of intrinsic dimension estimation.

Parameters:

• pred_line_ids (list) – list of predicted lines indices.

• pred_plane_ids (list) – list of predicted planes indices.

• metric (str) – accuracy metric, we support: f1, jaccard, balanced as in scikit learn.

Returns:

A tuple including

• [line_acc, line_pre, line_rec]: line accuracy, precision and recall.

• [plane_acc, plane_pre, plane_rec]: plane accuracy, precision and recall.

• [nonplane_acc, nonplane_pre, nonplane_rec]: 3D accuracy, precision and recall.

plot(ax, projection='3d', noisy=True, scales=(1.0, 1.0, 1.0), show_interpoints=False)

Plot simulated curve with its intrinsic segments.

Parameters:

• ax – plot axis.

• projection (str) – support 3d, xy, xz, yz modes.

• noisy (bool) – plot noisy curve or noiseless curve.

• scales (tuple of float) – specify x, y and z scales.

• show_interpoints (bool) – if True, points in between intrinsic segments are displayed.

class genepy3d.obj.curves.RadiusScaleTable(crv, r, mo=None, kerlen=None)

Bases: object

Compute scales table of a curve from given radius of curvature.

This is the support class for the function compute_local_3d_scale_radius().

crv

curve to compute radius scale table.

Type:

Curve

r

reference radius of curvature.

Type:

float

mo

method used in gaussian convolution to treat the extreme values.

Type:

str

kerlen

gaussian kernel length.

Type:

float

Methods:

build_scales_from_radius([max_iter])	Compute scales range from given reference radius.
compute_scales_from_radius([r])	Compute scales range from given reference radius.

build_scales_from_radius(max_iter=None)

Compute scales range from given reference radius.

Support func of intrinsic dimensions decomposition.

Parameters:

max_iter (int) – number of convolution iteration. If None, it will be set = haft of number of points.

Retuns:

kappa_tbl (array (float)): curvature table.

compute_scales_from_radius(r=None)

Compute scales range from given reference radius.

Support func of intrinsic dimensions decomposition.

Parameters:

r (float) – reference radius.

Retuns:

sig_c (float): central value of sigma delta_sig (float): std value of sigma

genepy3d.obj.curves.decompose_intrinsicdim_sequential(p, eps_line, eps_plane, eps_seg)

Decompose the curve into intrisic segments of 1D line, 2D plane and 3D.

This decomposition method is implemented from the paper of J. Yang et al. [1]. This func is only used to compare with our algorithm [2].

Parameters:

• p (array of float) – list of 3D points on curve.

• eps_line (float) – line threshold.

• eps_plane (float) – plane_threshold.

• eps_seg (float) – segment length threshold.

Returns:

list of line and plane indices.

References

[1]

Yang J, Yuan J, Li Y. Parsing 3D motion trajectory for gesture recognition. Journal of Visual Communication and Image Representation. 2016. Pages 627-640. DOI: 10.1016/j.jvcir.2016.04.010.

[2]

Phan MS, Matho K, Beaurepaire E, Livet J, Chessel A. nAdder: A scale-space approach for the 3D analysis of neuronal traces. 2022. PLOS Computational Biology 18(7): e1010211. DOI: 10.1371/journal.pcbi.1010211

genepy3d.obj.curves.align(target, ref)

Align a target curve to match a reference curve.

This func is used to align two curves that minimize their spatial differences.

The func is useful as a preparation before computing the distance between the two curves. E.g. we provided emd() in the Points class to compute the Earth Mover’s Distance between two point clouds.

Parameters:

• target (Curve) – target curve.

• ref (Curve) – reference curve.

Returns:

aligned target curve.

genepy3d.obj.points module

Methods for working with Points objects.

Classes:

Points(_coors)

Points cloud in 3D.

Functions:

emd(ps1, ps2[, return_flows])	Compute the Earth Mover's Distance (EMD) [1] between two points clouds ps1, ps2.
gen_ellipsoid([axes_length, n])	Generate random uniform points cloud on the surface of an ellipsoid.
gen_sphere([r, n])	Generate random uniform points cloud on a sphere.

class genepy3d.obj.points.Points(_coors)

Bases: object

Points cloud in 3D.

coors

list of 3d points.

• if array is given, then each column is the x, y, z coordinates.

• if tuple is given, then the items are the three arrays of x, y and z coordinates.

Type:

array | tuple

Examples

from genepy3d.obj import points
from sklearn.datasets import make_blobs
import matplotlib.pyplot as plt

# Generate two clusters of 3D points
coors, _ = make_blobs(n_features=3, centers = [(0, 0, 0), (5, 5, 5)])

# Create Points object
pnts = points.Points(coors)

# Plot
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
pnts.plot(ax)

Attributes:

size	Number of points.

Methods:

from_csv(filepath[, column_names, scales, args])	Create Points object from x, y and z coordinates in the csv file.
from_text(filepath)	Create Points object from x, y, z coordinates in the text file (e.g. txt, xyz).
center()	Center the cloud to the mean of each coordinates.
append(new_points)	Append a new points cloud to the current points cloud.
transform([phi, theta, psi, dx, dy, dz, zx, ...])	Transform the points cloud.
bbox()	Return bounding box covered the point cloud
fit_plane()	Fit a hyperplane to the points cloud.
pca()	Compute principal components analysis from the points cloud.
test_isotropy([n_iter])	Test if the points distribution is isotropic using bootstraping.
export_to_VTK(filepath)	Export to VTK file format.
to_curve()	Convert to Curve under the assumption that the points are listed in order
to_surface_qhull()	Convert to Surface by computing the convex hull of the points cloud.
to_surface_alphashape([alpha])	Convert to Surface by computing the alpha shape, i.e. the 'concave up to concavities of size alpha' hull.
plot(ax[, projection, point_args, equal_aspect])	Plot points cloud using matplotlib.

property size

Number of points.

classmethod from_csv(filepath, column_names=['x', 'y', 'z'], scales=(1.0, 1.0, 1.0), args={})

Create Points object from x, y and z coordinates in the csv file.

Parameters:

• filepath (str) – path of the csv file.

• column_names (list of str) – coordinates columns in the csv file.

• scales (tuple of float) – scales in x, y and z.

• args (dict) – overried parameters of pandas.read_csv().

Returns:

A Points object.

Examples

from genepy3d.obj import points

filepath = "path/to/your/csv/file.csv"

# Create Points object from column "x", "y", "z" in the given csv file
pnts = points.Points.from_csv(filepath)

classmethod from_text(filepath)

Create Points object from x, y, z coordinates in the text file (e.g. txt, xyz).

Parameters:

filepath (str) – path to the text file.

Returns:

A Points object.

Notes

We assume that the first three columns correspond to the x, y, z coordinates.

Examples

from genepy3d.obj import points

filepath = "path/to/your/text/file.txt"

# Create Points object from  the given text file
pnts = points.Points.from_text(filepath)

center()

Center the cloud to the mean of each coordinates.

The coordinates of points will be updated. Nothing is returned.

append(new_points)

Append a new points cloud to the current points cloud.

Parameters:

new_points (Points) – a Points object to be appended.

Notes

The coordinates of the current Points object will be updated.

Examples

from genepy3d.obj import points
from sklearn.datasets import make_blobs
import matplotlib.pyplot as plt

# Create Points object from blob cluster centered at (0,0,0)
coors, _ = make_blobs(n_features=3, centers = [(0, 0, 0)])
pnts_1 = points.Points(coors)

# Plot the pnts_1
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
pnts_1.plot(ax)

# Create Points object from blob cluster centered at (5,5,5)
coors, _ = make_blobs(n_features=3, centers = [(5, 5, 5)])
pnts_2 = points.Points(coors)

# Append pnts_2 into pnts_1
pnts_1.append(pnts_2)

# Plot the pnts_1 after appending the pnts_2
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
pnts_1.plot(ax)

transform(phi=0, theta=0, psi=0, dx=0, dy=0, dz=0, zx=1, zy=1, zz=1)

Transform the points cloud.

Parameters:

• phi (float) – rotation in x (in rad).

• theta (float) – rotation in y (in rad).

• psi (float) – rotation in z (in rad).

• dx (float) – translation in x.

• dy (float) – translation in y.

• dz (float) – translation in z.

• zx (float) – zoom in x.

• zy (float) – zoom in y.

• zz (float) – zoom in z.

Returns:

A Points object.

Examples

import numpy as np
from genepy3d.obj import points
from sklearn.datasets import make_blobs
import matplotlib.pyplot as plt

# Generate 3D points cloud contains two clusters
coors, _ = make_blobs(n_features=3, centers = [(0, 0, 0), (5, 5, 5)])
pnts = points.Points(coors)

# Plot the points cloud
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
pnts.plot(ax)

# Rotate points cloud around z by 90 degrees.
pnts_transformed = pnts.transform(psi=np.pi/2.)

# Plot the rotated points
pnts_transformed.plot(ax,point_args={'color':'r'})

bbox()

Return bounding box covered the point cloud

Returns:

top-left x y0 (float): top-left y z0 (float): top-left z xlen (float): length in x ylen (float): length in y zlen (float): length in z

Return type:

x0 (float)

fit_plane()

Fit a hyperplane to the points cloud.

Returns:

Tuple containing

• c (float): intercept scalar.

• normal (1d array): normal vector.

Examples

from genepy3d.obj import points
from sklearn.datasets import make_blobs
import matplotlib.pyplot as plt

# Generate 3D points cloud contains two clusters
coors, _ = make_blobs(n_features=3, centers = [(0, 0, 0), (5, 5, 5)])
pnts = points.Points(coors)

# Plot
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
pnts.plot(ax)

# Compute the z-intercept and the normal vector of the hyperplane that fit the point cloud
z0, normal = pnts.fit_plane()

pca()

Compute principal components analysis from the points cloud.

Returns:

Tuple containing

• (1d array): empirical mean of points cloud.

• (2d array): component vectors sorted by explained variances.

• (2d array): explained variances.

Examples

from genepy3d.obj import points
import matplotlib.pyplot as plt

# Generate random points on the surface of a 3D ellipsoid
pnts = points.gen_ellipsoid(axes_length=(1.5,1.,0.5),n=500)

# PCA
means, components, variances = pnts.pca()

# Plot
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
pnts.plot(ax,point_args={"alpha":0.3})
ax.quiver(means[0], means[1], means[2], components[0,0], components[0,1], components[0,2], length=2*variances[0], color='red')
ax.quiver(means[0], means[1], means[2], components[1,0], components[1,1], components[1,2], length=2*variances[1], color='red')
ax.quiver(means[0], means[1], means[2], components[2,0], components[2,1], components[2,2], length=2*variances[2], color='red')

test_isotropy(n_iter=500)

Test if the points distribution is isotropic using bootstraping.

Parameters:

n_iter (int) – number of bootstraping iterations.

Returns:

p_value.

export_to_VTK(filepath)

Export to VTK file format.

This is the wrapper of the pointsToVTK() in pyevtk package.

Parameters:

filepath (str) – full path of file to save, without extention.

Returns:

A VTK object.

to_curve()

Convert to Curve under the assumption that the points are listed in order

Returns:

A Curve object.

to_surface_qhull()

Convert to Surface by computing the convex hull of the points cloud.

Returns:

A Surface object.

to_surface_alphashape(alpha=None)

Convert to Surface by computing the alpha shape, i.e. the ‘concave up to concavities of size alpha’ hull. Assume a relatively homogeneous repartition of the points.

Please see surfaces.Surface.from_points_alpha_shape().

Returns:

a Surface object.

plot(ax, projection='3d', point_args={}, equal_aspect=True)

Plot points cloud using matplotlib.

Parameters:

• ax – axis to be plotted.

• projection (str) – support ‘3d’|’xy’|’xz’|’yz’ plot.

• point_args (dic) – matplotlib arguments of scatter().

• equal_aspect (bool) – make equal aspect for both axes.

genepy3d.obj.points.emd(ps1, ps2, return_flows=False)

Compute the Earth Mover’s Distance (EMD) [1] between two points clouds ps1, ps2.

Parameters:

• ps1 (Points) – point cloud.

• ps2 (Points) – point cloud.

• return_flows (bool) – if True then return maching flows between ps1 and ps2.

Returns:

EMD distance between ps1 and ps2 and if return_flows is True, then return array of matching flows.

References

[1]

https://en.wikipedia.org/wiki/Earth_mover%27s_distance

Examples

import numpy as np
from genepy3d.obj import points
import matplotlib.pyplot as plt

# Generate two points clouds
x = np.arange(10)
y = x
z = np.zeros(10)
pnt1 = points.Points((x,y,z))
pnt2 = pnt1.transform(dx=20) # translate in x

# EMD
dist, flows = points.emd(pnt1,pnt2,return_flows=True)
print(dist)
print(flows)

# Plot
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
pnt1.plot(ax)
pnt2.plot(ax)

genepy3d.obj.points.gen_ellipsoid(axes_length=[1.0, 1.0, 1.0], n=100)

Generate random uniform points cloud on the surface of an ellipsoid.

Parameters:

• axes_length (list of float) – half lengths of the major, median and minor axis.

• n (int) – number of points to be generated.

Returns:

A Points object.

Examples

from genepy3d.obj import points
import matplotlib.pyplot as plt

# Generate random points on the surface of 3D ellipsoid
pnts = points.gen_ellipsoid(axes_length=(1.5,1.,0.5),n=500)

# Plot
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
pnts.plot(ax)

genepy3d.obj.points.gen_sphere(r=1, n=1000)

Generate random uniform points cloud on a sphere.

Parameters:

• r (float) – sphere’s radius.

• n (int) – number of points.

Returns:

A Points object.

Examples

from genepy3d.obj import points
import matplotlib.pyplot as plt

# Generate random points on the surface of 3D sphere
pnts = points.gen_sphere(r=1,n=500)

# Plot
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
pnts.plot(ax)

genepy3d.obj.surfaces module

Methods for working with Surface objects.

Classes:

Surface([vertices, faces, face_normals, ...])

Triangulated surface in 3D; inherits from Trimesh and thus is triangulated upon creation if the mesh is not a triangulation.

Functions:

available_formats()

Return file formats that support by Trimesh for loading the surface.

class genepy3d.obj.surfaces.Surface(vertices: _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None = None, faces: _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None = None, face_normals: _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None = None, vertex_normals: _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None = None, face_colors: _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None = None, vertex_colors: _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None = None, face_attributes: dict[str, _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]] | None = None, vertex_attributes: dict[str, _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]] | None = None, metadata: dict[str, Any] | None = None, process: bool = True, validate: bool = False, merge_tex: bool | None = None, merge_norm: bool | None = None, use_embree: bool = True, initial_cache: dict[str, ndarray] | None = None, visual: ColorVisuals | TextureVisuals | None = None, **kwargs)

Bases: Trimesh

Triangulated surface in 3D; inherits from Trimesh and thus is triangulated upon creation if the mesh is not a triangulation.

The Trimesh library is here: https://trimsh.org/trimesh.html.

vertices

vertices coordinates.

Type:

ndarray (float)

faces

index of the vertices of each of the simplices.

Type:

ndarray (float)

Methods:

from_points_qhull(points)	Create surface object from point cloud via the convex hull algorithms.
from_points_alpha_shape(points[, alpha])	Create surface object from point cloud via the alpha-shape algorithms.
from_stl(afile[, xy2zratio, center])	Create surface object from STL file.
from_ply(afile[, xy2zratio, center])	Create surface object from PLY file.
from_file(afile[, xy2zratio, center])	Wrapper of trimesh.load() to load a surface from file.
from_volume(vol, lbl[, spacing, step_size, ...])	Create isosurface of a 3D voxels volume from image stack using marching cubes.
get_qhull()	Return the Convex Hull of the surface.
get_bbox()	Return the bounding box of the surface.
get_obbox()	Return the oriented bounding box of the surface.
get_obbox_axes()	Return the three basic axis-vectors covered the oriented bounding box.
compute_sphericity()	Compute the sphericity index of the surface.
get_paraboloid_centroid()	Find the centroid lying on the paraboloid.
split()	Override split() by casting Trimesh to Surface.
to_points()	Convert to Points object.
export_to_VTK(filepath[, kind])	Export surface in VTK fileformat, either as PolyData or unstructuredGrid.
plot(ax, **kwds)	Plot outline in 3D or 2D-projection (xy, xz and yz) using matplotlib.

classmethod from_points_qhull(points)

Create surface object from point cloud via the convex hull algorithms.

Parameters:

points (nx3 array of float) – coordinates of point cloud.

Returns:

a Surface object.

Examples

import numpy as np
from genepy3d.obj import surfaces
import matplotlib.pyplot as plt

# coors = np.array(...): array whose columns are x, y and z coordinates

surf = surfaces.Surface.from_points_qhull(coors) # surface from convex hull

fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
surf.plot(ax);

classmethod from_points_alpha_shape(points, alpha=None)

Create surface object from point cloud via the alpha-shape algorithms.

The result is the set triangles part of tetrahedrals with a circumsphere of radius less that alpha that face outward.

Parameters:

• points (nx3 array of float) – coordinates of point cloud.

• alpha (float) – value of alpha.

Returns:

a Surface instace.

Notes

There is no guarantee on the topology of the resulting mesh.

Examples

import numpy as np
from genepy3d.obj import surfaces
from genepy3d.util import geo
import matplotlib.pyplot as plt

# coors = np.array(...): array whose columns are x, y and z coordinates

# Diagonal distance of the bounding box of the point cloud
# Alpha value can be set relatively to this distance
# Larger value of alpha produces surface closer to the convex hull
min_coors = np.min(coors,axis=0)
max_coors = np.max(coors,axis=0)
diadst = geo.l2(min_coors,max_coors) # diagonal distance
alpha = diadst/3. # alpha value

# Surface from the alpha shape
surf = surfaces.Surface.from_points_alpha_shape(coors,alpha)

fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
surf.plot(ax)

classmethod from_stl(afile, xy2zratio=None, center=False)

Create surface object from STL file.

This function is DEPRECATED. use Surface.from_file() instead.

Parameters:

• afile (string) – filepath to load the surface.

• xy2zratio (float) – in case of anisotromic measurement, ratio to apply to z coordinates w.r.t. x and y.

• center (bool) – center the mesh.

Returns:

a Surface object.

classmethod from_ply(afile, xy2zratio=None, center=False)

Create surface object from PLY file.

This function is DEPRECATED. use Surface.from_file() instead.

Parameters:

• afile (string) – filepath to load the surface.

• xy2zratio (float) – in case of anisotromic measurement, ratio to apply to z coordinates w.r.t. x and y.

• center (bool) – center the mesh.

Returns:

a Surface object.

classmethod from_file(afile, xy2zratio=None, center=False)

Wrapper of trimesh.load() to load a surface from file.

To see which file formats are supported. Use surfaces.available_formats().

Parameters:

• afile (string) – filepath to load the surface.

• xy2zratio (float) – in case of anisotromic measurement, ratio to apply to z coordinates w.r.t. x and y.

• center (bool) – center the mesh.

Returns:

a Surface object.

Examples

from genepy3d.obj import surfaces
import matplotlib.pyplot as plt

filepath = "path/to/file"
surf = surfaces.Surface.from_file(filepath)

fig = plt.figure()
ax = fig.add_subplot(111,projection="3d")
surf.plot(ax)

classmethod from_volume(vol, lbl, spacing=(1.0, 1.0, 1.0), step_size=1, level=0.5, opening=False, fillholes=False)

Create isosurface of a 3D voxels volume from image stack using marching cubes.

We assume that the image stack contains multiple volumic objects specified by integer labels. This function extracts the volume of a given label, then estimates the outlined surface using marching cubes from scikit-image.

Parameters:

• vol (3D ndarray) – 3D volume.

• lbl (int) – voxel label of the volume to be extracted from the image stack.

• spacing (tuple (float)) – voxel spacing using in marching cubes.

• step_size (int) – parameter defining the mesh resolution (higher step_size yields coarser mesh).

• level (float) – between 0 and 1, determine proportion between background (0.) and foreground (1.) in volume.

• opening (bool) – preprocess volumne by morphology opening.

• fillholes (bool) – preprocess volumne by morphology fillholes.

Returns:

a Surface object.

Examples

import numpy as np
from genepy3d.obj import surfaces
import matplotlib.pyplot as plt

stack = np.zeros((100,100,100)) # an image stack of size 100x100x100
stack[20:40,20:40,20:40] = 1 # a volume of label = 1

surf = surfaces.Surface.from_volume(stack,1) # outline surface of volume

fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
surf.plot(ax);

get_qhull()

Return the Convex Hull of the surface.

Returns:

a Surface object.

get_bbox()

Return the bounding box of the surface.

Returns:

a Surface object.

get_obbox()

Return the oriented bounding box of the surface.

Returns:

a Surface object.

get_obbox_axes()

Return the three basic axis-vectors covered the oriented bounding box.

Returns:

• coordinates of the three axes.

• lengths of the three axes.

Return type:

Tuple contains

compute_sphericity()

Compute the sphericity index of the surface.

The sphericity formula is from this paper: https://www.sciencedirect.com/science/article/pii/S1877750318304757#bib0205

Sphericity = ((minor**2) / (median*major))**(1/3),

where minor, median and major are lengths of the elliposoid fitting the surface. These lengths can be found from the orient bounding box of the surface.

Returns:

a float.

get_paraboloid_centroid()

Find the centroid lying on the paraboloid.

Assuming the surface has a paraboloid form. This function estimate the centroid lying on the surface. It’s not similar to the centroid computed by averaging faces positions.

split()

Override split() by casting Trimesh to Surface.

If the surface consists of distinct subsurfaces, then the function split() returns these subsurfaces.

Returns:

a list of distinct subsurfaces.

to_points()

Convert to Points object.

Returns:

a Points object, consisting of the vertices of the mesh.

export_to_VTK(filepath, kind='PolyData')

Export surface in VTK fileformat, either as PolyData or unstructuredGrid.

Parameters:

• filepath (string) – path to export the file to, with no extention

• kind (string) – one of ‘PolyData’ or ‘UnstructuredGrid’ (default:’PolyData’).

Returns:

a VTK fileformat.

plot(ax, **kwds)

Plot outline in 3D or 2D-projection (xy, xz and yz) using matplotlib.

Parameters:

• ax – axis to be plotted.

• projection (str) – support ‘3d’|’xy’|’xz’|’yz’.

• scales – (tuples [float]): scales in x y and z.

• args_3d (dic) – matplotlib arguments to plot 3d boundary (set of 3d polygons).

• args_2d (dic) – matplotlib arguments to plot 2d boundary.

• equal_aspect (bool) – make equal aspect for both axes.

genepy3d.obj.surfaces.available_formats()

Return file formats that support by Trimesh for loading the surface.

Returns:

list of support file formats.

genepy3d.obj.trees module

Methods for working with Tree objects.

Classes:

Tree(_nodes[, _id, _name])

Tree in 3D.

class genepy3d.obj.trees.Tree(_nodes, _id=0, _name='GeNePy3D')

Bases: object

Tree in 3D. The tree structure is constructed based on anytree package.

We assume there are no duplicated nodes on the Tree.

nodes

structure of tree nodes.

Type:

anytree

id

tree ID (optional).

Type:

int

name

tree name (optional).

Type:

str

Methods:

from_table(nodes_tbl[, connectors_tbl, ...])	Build tree from tables.
from_swc(filepath)	Build tree from a SWC file.
from_eswc(filepath)	Build tree from ESWC file.
from_ims(filepath)	Build tree from Imaris file.
from_catmaid_server(catmaid_host, token, ...)	Build tree from CATMAID [1].
get_root()	Return root node ID.
get_parent(nodeid)	Return parent id of node id.
get_children(nodeid)	Return children ids of nodeid.
get_preorder_nodes([rootid])	Return list of nodes in tree.
get_number_nodes([rootid])	Return the total number of nodes in tree.
get_leaves([rootid])	Return list of leaf node IDs.
get_branchingnodes([rootid])	Return list of branching node IDs.
get_connectors([rootid])	Return list of connector node IDs.
get_coordinates([nodeid, rootid])	Return coordinates of given nodes.
get_features([features, nodeid, rootid])	Return list of features from given nodes.
get_bbox([nodeid, rootid])	Return bounding box coordinates covering a list of nodes from the tree.
compute_length([nodeid, rootid])	Return the total length of the tree or subpart of the tree.
compute_spine([rootid])	Return spine node IDs (longest branch from the root to a leaf).
compute_strahler_order([nodeid, rootid])	Return Strahler order for given node IDs [1].
path(target[, source, rootid])	Find list of nodes going from source node to target node.
extract_subtrees(nodeid[, to_children, ...])	Extract a sub tree from the given node.
copy([rootid])	Make a copy of the tree.
prune_leaves([nodeid, length, rootid])	Prune leaf branches based on its length.
decompose_segments([rootid])	Decompose tree in segments separating by the branching nodes.
decompose_leaves([rootid])	Decompose tree into segments starting from root to every leaf.
decompose_spines([rootid])	Decompose tree into list of spines (longest branches) starting from the root node.
resample([sampling_length, rootid, ...])	Resample the tree by a given sampling length.
compute_angles([decomposed_method, sigma, ...])	Return angle of each tree node.
compute_angles_vector([decomposed_method, ...])	Return angles of each tree node with x, y and z axes.
compute_local_3d_scale_sigma(sig_lst[, ...])	Compute local 3d scale of neuron from a list of sigma of Gaussian.
compute_local_3d_scale_radius(r_lst[, ...])	Compute local 3d scale of neuron from a list of radius of curvatures.
summary()	Return a brief summary of the tree.
to_curve([nodeid])	Convert a segment of tree given by the list of node ID to a Curve object.
to_points([nodeid])	Convert a segment of tree given by the list of node ID to a Points object.
plot(ax[, projection, rootid, spine_only, ...])	Plot tree using matplotlib.

classmethod from_table(nodes_tbl, connectors_tbl=None, tree_id=0, tree_name='genepy3d')

Build tree from tables.

Parameters:

• nodes_tbl (pandas dataframe) – tree nodes table.

• connectors_tbl (pandas dataframe) – connector table (optional).

• tree_id (int) – tree ID (optional).

• tree_name (str) – tree name (optional).

Returns:

A Tree object.

Notes

• It needs at least nodes_tbl to create the tree object. The nodes_tbl stores information of nodes and their linked parent-nodes. It is the same as in the SWC file specification. Please see this [1] for more detail. The column names must include:

  • treenode_id: ID of node

  • structure_id: ID of node structure (see the column Structure Identifier in [1] for the meaning).

  • x, y, z: coordinates of node

  • r: radius of node

  • parent_treenode_id: ID of parent node

• The table connectors_tbl is optional. It is used to describe the synaptic connection between neuronal trees. It is the same as connector table in CATMAID [2]. It must consist of columns:

  • connector_id: ID of a connector

  • treenode_id: ID of a treenode from a neuronal tree object

  • relation_id: connectomic relation, must be ‘presynaptic_to’ or ‘postsynaptic_to’ (see [3] to understand about connectomic relation in CATMAID)

References

[1]

http://www.neuronland.org/NLMorphologyConverter/MorphologyFormats/SWC/Spec.html

[2]

https://catmaid.readthedocs.io/en/stable/index.html

[3]

https://catmaid.readthedocs.io/en/stable/tracing_neurons.html

Examples

import pandas as pd
from genepy3d.obj import trees
import matplotlib.pyplot as plt

# Generate a dummy nodes data
data = [[0,0,0.,0.,0.,1.,-1],
        [1,0,1.,2.,3.,1.,0],
        [2,0,4.,5.,6.,1.,1]]

# Dataframe from nodes data
node_tbl = pd.DataFrame(data,columns=['treenode_id','structure_id','x','y','z','r','parent_treenode_id'])

# Create a dummy neuron from node_tbl
neuron = trees.Tree.from_table(node_tbl)

# Plot it
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
neuron.plot(ax)

classmethod from_swc(filepath)

Build tree from a SWC file.

See this [1] for the description of SWC format.

Parameters:

filepath (str) – path to swc file.

Returns:

A Tree object.

References

[1]

http://www.neuronland.org/NLMorphologyConverter/MorphologyFormats/SWC/Spec.html

Examples

from genepy3d.obj import trees
import matplotlib.pyplot as plt

# An example of swc file can be downloaded from:
# https://neuromorpho.org/neuron_info.jsp?neuron_name=AD1202-PMC-L-4201-fiber129
swc_path = "file/to/swc/file"
neu = trees.Tree.from_swc(swc_path)

fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
neu.plot(ax)

classmethod from_eswc(filepath)

Build tree from ESWC file. The ESWC is used e.g. in Vaa3D-Neuron software [1].

ESWC is an extension of SWC with additional columns ‘segment_id’, ‘level’, ‘mode’, ‘timestamp’, ‘featureval’.

Parameters:

filepath (str) – path to eswc file.

Returns:

A Tree object.

Notes

ESWC has extended columns compared to default SWC file, but we don’t pass these extended columns into the Tree object. We only pass the similar columns to the Tree object as in the SWC file.

References

[1]

Vaa3D/Vaa3D_Wiki

Examples

Reading ESWC should be the same as reading SWC. Please see example from from_swc().

classmethod from_ims(filepath)

Build tree from Imaris file.

Parameters:

filepath (str) – path to ims file.

Returns:

A Tree object.

classmethod from_catmaid_server(catmaid_host, token, project_id, neuron_id)

Build tree from CATMAID [1].

This function queries from the CATMAID server a neuronal trace given by a project_id and neuron_id. The neuronal trace data is fetched using CATMAID API [2].

Parameters:

• catmaid_host (str) – address of CATMAID server

• token (str) – authenticated string

• project_id (int) – project ID

• neuron_id (int) – neuron ID

Returns:

A Tree object.

References

[1]

https://catmaid.readthedocs.io/en/stable/index.html

[2]

https://catmaid.readthedocs.io/en/stable/api.html

Examples

from genepy3d.obj import trees

# Replace these fake values with your own
catmaid_host = 'https://your-catmaid-host.com'
token = "9944b09199c62bcf9418ad846dd0e4bbdfc6ee4b"
project_id = 1
neuron_id = 1

# Retreive neuronal tree from from CATMAID host
neuron = trees.Tree.from_catmaid_server(catmaid_host,token,project_id,neuron_id)

get_root()

Return root node ID.

Returns:

List of root ID (list[int]).

get_parent(nodeid)

Return parent id of node id.

Parameters:

nodeid (int) – a node ID.

Returns:

Node ID (int).

get_children(nodeid)

Return children ids of nodeid.

Parameters:

nodeid (int) – a node ID.

Returns:

List of node ID (list[int]).

get_preorder_nodes(rootid=None)

Return list of nodes in tree. The order of nodes in the list follows the preorder traversal [1].

Parameters:

rootid (int) – root node ID. If None, then the first rootid will be taken.

Returns:

List of node ID (list[int]).

References

[1]

https://en.wikipedia.org/wiki/Tree_traversal

get_number_nodes(rootid=None)

Return the total number of nodes in tree.

Parameters:

rootid (int) – root node ID. If None, then the first rootid will be taken.

Returns:

Total number of nodes (int).

get_leaves(rootid=None)

Return list of leaf node IDs.

Parameters:

rootid (int) – root node ID. If None, then the first rootid will be taken.

Returns:

list of int.

get_branchingnodes(rootid=None)

Return list of branching node IDs.

Parameters:

rootid (int) – root node ID. If None, then the first rootid will be taken.

Returns:

list of int.

get_connectors(rootid=None)

Return list of connector node IDs.

Parameters:

rootid (int) – root node ID. If None, then the first rootid will be taken.

Returns:

Pandas dataframe whose indices are node IDs and values are corresponding connector types, connector ids.

get_coordinates(nodeid=None, rootid=None)

Return coordinates of given nodes.

Parameters:

• nodes (int | list of int) – list of node IDs. If None, then return the coordinates of all nodes from the given root ID.

• rootid (int) – root node ID. If None, then the first rootid will be taken.

Returns:

Pandas dataframe whose indices are node IDs and values are corresponding coordinates.

get_features(features=['x', 'y', 'z'], nodeid=None, rootid=None)

Return list of features from given nodes.

Parameters:

• features (list of str) – list of node-wise features, e.g. “x”, “y”, “z”, “r”, “structure_id”, “connector_relation”. If None, then return all features.

• nodes (int | list of int) – list of node IDs. If None, then return the coordinates of all nodes from the given root ID.

• rootid (int) – root node ID. If None, then the first rootid will be taken.

Returns:

Pandas dataframe whose indices are node IDs and values are corresponding features.

get_bbox(nodeid=None, rootid=None)

Return bounding box coordinates covering a list of nodes from the tree.

Parameters:

• nodeid (list[int]) – list of node IDs. If None, then return bounding box of the entire tree.

• rootid (int) – root node ID. If None, then the first rootid will be taken.

Returns:

List of [(Xmin,Xmax),(Ymin,Ymax),(Zmin,Zmax)]

compute_length(nodeid=None, rootid=None)

Return the total length of the tree or subpart of the tree.

Parameters:

• nodeid (list[int]) – list of node IDs. If None, then return the total length of the entire tree.

• rootid (int) – root node ID. If None, then the first rootid will be taken.

Returns:

int.

Notes

If you provide a list of nodes. We consider it as a curve, then we return the curve length.

Examples

import pandas as pd
from genepy3d.obj import trees
import matplotlib.pyplot as plt

# Generate a dummy nodes data
data = [[0,0,0.,0.,0.,1.,-1],
        [1,0,1.,2.,3.,1.,0],
        [2,0,4.,5.,6.,1.,1],
        [3,0,5.,5.,5.,1.,1]]

# Dataframe from nodes data
node_tbl = pd.DataFrame(data,columns=['treenode_id','structure_id','x','y','z','r','parent_treenode_id'])

# Create a dummy neuron from node_tbl
neuron = trees.Tree.from_table(node_tbl)

# Plot it
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
neuron.plot(ax)
plt.tight_layout();

# Total length of the entire neuron
total_length = neuron.compute_length()
print(total_length)

# Get a segment from root to the first branching point
segment_lst = neuron.decompose_segments()
print(segment_lst)
segment = segment_lst['0_1']

# Length of that segment
segment_length = neuron.compute_length(segment)
print(segment_length)

compute_spine(rootid=None)

Return spine node IDs (longest branch from the root to a leaf).

Parameters:

rootid (int) – root node ID. If None, then the first rootid will be taken.

Returns:

list of int.

compute_strahler_order(nodeid=None, rootid=None)

Return Strahler order for given node IDs [1].

Parameters:

• nodeid (int|list[int]) – list of node IDs. If None, return all the Strahler order for all nodes from the given root ID.

• rootid (int) – root node ID. If None, then the first rootid will be taken.

Returns:

Pandas serie whose indices are node IDs and values are corresponding strahler orders.

References

[1]

https://en.wikipedia.org/wiki/Strahler_number

path(target, source=None, rootid=None)

Find list of nodes going from source node to target node.

Parameters:

• target (int) – target node ID.

• source (int) – source node ID. If None, then we take root as the source.

• rootid (int) – root node ID. If None, then the first rootid will be taken.

Returns:

list of traveled node IDs included source and target node IDs.

Examples

import pandas as pd
from genepy3d.obj import trees
import matplotlib.pyplot as plt

# Generate a dummy nodes data
data = [[0,0,0.,0.,0.,1.,-1],
        [1,0,1.,2.,3.,1.,0],
        [2,0,4.,5.,6.,1.,1],
        [3,0,5.,5.,5.,1.,1]]

# Dataframe from nodes data
node_tbl = pd.DataFrame(data,columns=['treenode_id','structure_id','x','y','z','r','parent_treenode_id'])

# Create a dummy neuron from node_tbl
neuron = trees.Tree.from_table(node_tbl)

# Get ID of a leaf
leaf_node = neuron.get_leaves()[0]
print(leaf_node)

# Path from root to the leaf
print(neuron.path(leaf_node))

extract_subtrees(nodeid, to_children=True, separate_children=False, rootid=None)

Extract a sub tree from the given node.

Parameters:

• nodeid (int) – node ID from that we extract the subtree. It becomes new root node in the subtree.

• to_children (bool) – if False, then extract upper subtree (toward parent). Otherwise, lower subtrees are extracted (toward children).

• separate_children (bool) – if True, then each subtree is extracted for each child. This parameter only works when setting to_children = True.

• rootid (int) – root node ID. If None, then the first rootid will be taken.

Returns:

list of Tree objects.

Examples

import pandas as pd
from genepy3d.obj import trees
import matplotlib.pyplot as plt

# Generate a dummy nodes data
data = [[0,0,0.,0.,0.,1.,-1],
        [1,0,1.,2.,3.,1.,0],
        [2,0,4.,5.,6.,1.,1],
        [3,0,7.,2.,2.,1.,1],
        [4,0,6.,6.,6.,1.,2]]

# Dataframe from nodes data
node_tbl = pd.DataFrame(data,columns=['treenode_id','structure_id','x','y','z','r','parent_treenode_id'])

# Create a dummy neuron from node_tbl
neuron = trees.Tree.from_table(node_tbl)

# Plot it
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
neuron.plot(ax)
plt.tight_layout();

# Get a branching node
brnode = neuron.get_branchingnodes()[0]
print(brnode)

# Extract a subtree from the branching node toward its children
subtree = neuron.extract_subtrees(brnode,to_children=True,separate_children=False)

# Plot it
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
subtree.plot(ax)
plt.tight_layout();

# Extract list of subtress from the branching node toward its children
# Each subtree is extracted from each child
subtree_lst = neuron.extract_subtrees(brnode,to_children=True,separate_children=True)

# Plot first subtree
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
subtree_lst[0].plot(ax)
plt.tight_layout();

# Plot second subtree
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
subtree_lst[1].plot(ax)
plt.tight_layout();

copy(rootid=None)

Make a copy of the tree.

Parameters:

rootid (int) – root node ID. If None, then the first rootid will be taken.

Returns:

A Tree object.

prune_leaves(nodeid=None, length=None, rootid=None)

Prune leaf branches based on its length.

Parameters:

• nodeid (list [int]) – list of leaves IDs used to examine the pruning. If None, then take all leaves into account.

• length (float) – length threshold for pruning the leaf branches from the leaves given by nodeid. If None, then we eliminate all the leaf branches.

• rootid (int) – root node ID. If None, then the first rootid will be taken.

Returns:

A new Tree object after pruning.

Examples

import pandas as pd
from genepy3d.obj import trees
import matplotlib.pyplot as plt

# Generate a dummy nodes data
data = [[0,0,0.,0.,0.,1.,-1],
        [1,0,1.,2.,3.,1.,0],
        [2,0,4.,5.,6.,1.,1],
        [3,0,7.,2.,2.,1.,1],
        [4,0,6.,6.,6.,1.,2],
        [5,0,3.,7.,6.,1.,2]]

# Dataframe from nodes data
node_tbl = pd.DataFrame(data,columns=['treenode_id','structure_id','x','y','z','r','parent_treenode_id'])

# Create a dummy neuron from node_tbl
neuron = trees.Tree.from_table(node_tbl)

# Plot it
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
neuron.plot(ax)
plt.tight_layout();

# Prune all the leaf branches
neuron_pruned = neuron.prune_leaves()

# Plot pruned tree
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
neuron_pruned.plot(ax)
plt.tight_layout();

# Prune only leaf branches whose length < 5 microns
neuron_pruned = neuron.prune_leaves(length=5)

# Plot pruned tree
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
neuron_pruned.plot(ax)
plt.tight_layout();

decompose_segments(rootid=None)

Decompose tree in segments separating by the branching nodes.

Parameters:

rootid (int) – root node ID. If None, then the first rootid will be taken.

Returns:

A dictionary whose each item is a decomposed segment, where

• key is a string represents the first node and the last node of the decomposed segment.

• value is a list of nodes IDs of the decomposed segment.

Examples

import pandas as pd
from genepy3d.obj import trees
import matplotlib.pyplot as plt

# Generate a dummy nodes data
data = [[0,0,0.,0.,0.,1.,-1],
        [1,0,1.,2.,3.,1.,0],
        [2,0,4.,5.,6.,1.,1],
        [3,0,7.,2.,2.,1.,1],
        [4,0,6.,6.,6.,1.,2],
        [5,0,3.,7.,6.,1.,2],
        [6,0,5.,5.,9.,1.,5]]

# Dataframe from nodes data
node_tbl = pd.DataFrame(data,columns=['treenode_id','structure_id','x','y','z','r','parent_treenode_id'])

# Create a dummy neuron from node_tbl
neuron = trees.Tree.from_table(node_tbl)

# Plot neuron
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
neuron.plot(ax)

# Plot node ID
coors = neuron.get_coordinates()
for nodeid in coors.index:
    tmp = coors.loc[nodeid]
    ax.text(tmp.x-0.5,tmp.y+0.1,tmp.z+0.1,nodeid,fontsize=12)

ax.axis('off');
plt.tight_layout();

# Decompose into segments separating by branching nodes
print(neuron.decompose_segments())

# The output should be as follows:
# {'2_4': [2, 4], '2_6': [2, 5, 6], '1_3': [1, 3], '0_1': [0, 1], '1_2': [1, 2]}

decompose_leaves(rootid=None)

Decompose tree into segments starting from root to every leaf.

Parameters:

rootid (int) – root node ID. If None, then the first rootid will be taken.

Returns:

A dictionary whose each item is a decomposed segment, where

• key is a string represents the leaf node of the decomposed segment.

• value is a list of nodes IDs of the decomposed segment.

Examples

import pandas as pd
from genepy3d.obj import trees
import matplotlib.pyplot as plt

# Generate a dummy nodes data
data = [[0,0,0.,0.,0.,1.,-1],
        [1,0,1.,2.,3.,1.,0],
        [2,0,4.,5.,6.,1.,1],
        [3,0,7.,2.,2.,1.,1],
        [4,0,6.,6.,6.,1.,2],
        [5,0,3.,7.,6.,1.,2],
        [6,0,5.,5.,9.,1.,5]]

# Dataframe from nodes data
node_tbl = pd.DataFrame(data,columns=['treenode_id','structure_id','x','y','z','r','parent_treenode_id'])

# Create a dummy neuron from node_tbl
neuron = trees.Tree.from_table(node_tbl)

# Plot neuron
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
neuron.plot(ax)

# Plot node ID
coors = neuron.get_coordinates()
for nodeid in coors.index:
    tmp = coors.loc[nodeid]
    ax.text(tmp.x-0.5,tmp.y+0.1,tmp.z+0.1,nodeid,fontsize=12)

ax.axis('off');
plt.tight_layout();

# Decompose into segments from root to every leaves
print(neuron.decompose_leaves())

# The output should be as follows:
# {4: [0, 1, 2, 4], 6: [0, 1, 2, 5, 6], 3: [0, 1, 3]}

decompose_spines(rootid=None)

Decompose tree into list of spines (longest branches) starting from the root node.

We first extract the longest branch from the root to a leaf. Then for every branching node on that longest branch, we extract the subtree from that, and repeat the process by extracting the longest branch of that subtree. The process is repeated until no branch can be extracted.

Parameters:

rootid (int) – root node ID. If None, then the first rootid will be taken.

Returns:

A dictionary whose each item is a decomposed segment, where

• key is a string represents the first node and the last node of the decomposed segment.

• value is a list of nodes IDs of the decomposed segment.

Examples

import pandas as pd
from genepy3d.obj import trees
import matplotlib.pyplot as plt

# Generate a dummy nodes data
data = [[0,0,0.,0.,0.,1.,-1],
        [1,0,1.,2.,3.,1.,0],
        [2,0,4.,5.,6.,1.,1],
        [3,0,7.,2.,2.,1.,1],
        [4,0,6.,6.,6.,1.,2],
        [5,0,3.,7.,6.,1.,2],
        [6,0,5.,5.,9.,1.,5]]

# Dataframe from nodes data
node_tbl = pd.DataFrame(data,columns=['treenode_id','structure_id','x','y','z','r','parent_treenode_id'])

# Create a dummy neuron from node_tbl
neuron = trees.Tree.from_table(node_tbl)

# Plot neuron
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
neuron.plot(ax)

# Plot node ID
coors = neuron.get_coordinates()
for nodeid in coors.index:
    tmp = coors.loc[nodeid]
    ax.text(tmp.x-0.5,tmp.y+0.1,tmp.z+0.1,nodeid,fontsize=12)

ax.axis('off');
plt.tight_layout();

# Decompose into segments of longest branches
print(neuron.decompose_spines())

# The output should be as follows:
# {'0_6': [0, 1, 2, 5, 6], '1_3': [1, 3], '2_4': [2, 4]}

resample(sampling_length=None, rootid=None, spline_order=1, spline_smoothness=0, decompose_method='branching')

Resample the tree by a given sampling length. Spline interpolation is used for resampling.

The tree is first decomposed into segments either by “branching” (decompose_segments()) or by “spine” (decompose_spines()). Then, each segment is resampled by the given sampling length.

Parameters:

• sampling_length (float) – sampling length. If None, then we resample the tree by the same number of nodes.

• rootid (int) – root node ID. If None, then the first rootid will be taken.

• spline_order (uint) – Spline order. value of 1 means linear interpolation.

• spline_smoothness (float or None) – Parameter to control the smoothness of spline interpolation. The parameter is similar to the s parameter from this [1].

• decompose_method (str) – choose the mode of decomposition: “branching” or “spine”.

Returns:

A new Tree object with resampled nodes.

References

[1]

https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.splprep.html

Examples

import pandas as pd
from genepy3d.obj import trees
import matplotlib.pyplot as plt

# Generate a dummy nodes data
data = [[0,0,0.,0.,0.,1.,-1],
        [1,0,1.,2.,3.,1.,0],
        [2,0,4.,5.,6.,1.,1],
        [3,0,7.,2.,2.,1.,1],
        [4,0,6.,6.,6.,1.,2],
        [5,0,3.,7.,6.,1.,2],
        [6,0,5.,5.,9.,1.,5]]

# Dataframe from nodes data
node_tbl = pd.DataFrame(data,columns=['treenode_id','structure_id','x','y','z','r','parent_treenode_id'])

# Create a dummy neuron from node_tbl
neuron = trees.Tree.from_table(node_tbl)

# Plot neuron, only showing the node points
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
neuron.plot(ax,show_nodes=True,show_leaves=False,show_root=False,show_branchingnodes=False)

print("Nb. of nodes before resampling:",neuron.get_number_nodes())

# Resampling by sampling length = 1 micron
neuron_resampled = neuron.resample(sampling_length=1)

# Plot resampled neuron, only showing the node points
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
neuron_resampled.plot(ax,show_nodes=True,show_leaves=False,show_root=False,show_branchingnodes=False)

print("Nb. of nodes after resampling:",neuron_resampled.get_number_nodes())

compute_angles(decomposed_method='branching', sigma=0, rootid=None)

Return angle of each tree node.

The tree is decomposed into segments. Then, the angle of node N is the angle between two vectors (N => N-1) and (N => N+1). One node can have multiple angle values depending on the decomposition mode.

Parameters:

• decomposed_method (str) – “leaf”, “spine” or “branching”, method for decomposing the neuronal trace.

• sigma (float) – sigma used in Gaussian function for smoothing the tree branches.

• rootid (int) – root node ID. If None, then the first rootid will be taken.

Returns:

Pandas dataframe containing angles at every nodes.

Examples

import pandas as pd
from genepy3d.obj import trees
import matplotlib.pyplot as plt

# Generate a dummy nodes data
data = [[0,0,0.,0.,0.,1.,-1],
        [1,0,1.,2.,3.,1.,0],
        [2,0,4.,5.,6.,1.,1],
        [3,0,7.,2.,2.,1.,1],
        [4,0,6.,6.,6.,1.,2],
        [5,0,3.,7.,6.,1.,2],
        [6,0,5.,5.,9.,1.,5]]

# Dataframe from nodes data
node_tbl = pd.DataFrame(data,columns=['treenode_id','structure_id','x','y','z','r','parent_treenode_id'])

# Create a dummy neuron from node_tbl
neuron = trees.Tree.from_table(node_tbl)

# Plot neuron, only showing the node points
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
neuron.plot(ax,show_nodes=True,show_leaves=False,show_root=False,show_branchingnodes=False)

# Plot node ID
coors = neuron.get_coordinates()
for nodeid in coors.index:
    tmp = coors.loc[nodeid]
    ax.text(tmp.x-0.5,tmp.y+0.1,tmp.z+0.1,nodeid,fontsize=12)

plt.tight_layout();

# Compute the node angles with x, y and z axes
# The default decomposition is "branching"
neuron.compute_angle()

compute_angles_vector(decomposed_method='branching', sigma=0, rootid=None)

Return angles of each tree node with x, y and z axes.

Angle at at a node N is the angle between vector (N => N + 1) and three unit vectors of x, y and z axes. The tree is decomposed into segments. Then, the angles of nodes on each segment are computed. One node can have multiple angle values depending on the decomposition mode.

Parameters:

• decomposed_method (str) – “leaf”, “spine” or “branching”, method for decomposing the neuronal trace.

• sigma (float) – sigma used in Gaussian function for smoothing the tree branches.

• rootid (int) – root node ID. If None, then the first rootid will be taken.

Returns:

Pandas dataframe containing angles at every nodes.

Examples

import pandas as pd
from genepy3d.obj import trees
import matplotlib.pyplot as plt

# Generate a dummy nodes data
data = [[0,0,0.,0.,0.,1.,-1],
        [1,0,1.,2.,3.,1.,0],
        [2,0,4.,5.,6.,1.,1],
        [3,0,7.,2.,2.,1.,1],
        [4,0,6.,6.,6.,1.,2],
        [5,0,3.,7.,6.,1.,2],
        [6,0,5.,5.,9.,1.,5]]

# Dataframe from nodes data
node_tbl = pd.DataFrame(data,columns=['treenode_id','structure_id','x','y','z','r','parent_treenode_id'])

# Create a dummy neuron from node_tbl
neuron = trees.Tree.from_table(node_tbl)

# Plot neuron, only showing the node points
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
neuron.plot(ax,show_nodes=True,show_leaves=False,show_root=False,show_branchingnodes=False)

# Plot node ID
coors = neuron.get_coordinates()
for nodeid in coors.index:
    tmp = coors.loc[nodeid]
    ax.text(tmp.x-0.5,tmp.y+0.1,tmp.z+0.1,nodeid,fontsize=12)

plt.tight_layout();

# Compute the node angles with x, y and z axes
# The default decomposition is "branching"
neuron.compute_angle_axes()

compute_local_3d_scale_sigma(sig_lst, dim_param=None, decomposed_method='leaf', rootid=None, postprocess='mean')

Compute local 3d scale of neuron from a list of sigma of Gaussian.

The tree is decomposed into segments. Then, local 3d scale is computed for each segment. Each segment is considered as a Curve object. Please see compute_local_3d_scale_sigma() from Curve class for more detail. One node can have multiple local 3d scales due to the decomposition. You can specify the postprocess to return e.g. the mean local 3d scale of the node.

Parameters:

• sig_lst (list of float) – list of sigma values of Gaussian.

• dim_param (dic) – parameters for dimension decomposition. If None, default parameters are used. See parameters of compute_local_3d_scale_sigma() from Curve class for the meaning.

• decomposed_method (str) – “leaf”, “spine” or “branching”, method for decomposing the neuronal tree.

• rootid (int) – root node ID. If None, then the first rootid will be taken.

• postprocess (str) – support “mean”, “std”, “max”, “min” or None.

Returns:

1D, 2D, 3D flags, curvature and torsion.

Return type:

Pandas dataframe containing local 3D scale and additional information

Notes

If postprocess = None, then a Pandas dataframe containing local 3d scale and additional features (curvature, torsion, etc) of each decomposed segment is returned. Otherwise, a Pandas serie containing only the local 3d scale of each node is returned.

Examples

import pandas as pd
from genepy3d.obj import trees
import matplotlib.pyplot as plt

# Generate a dummy nodes data
data = [[0,0,0.,0.,0.,1.,-1],
        [1,0,1.,2.,3.,1.,0],
        [2,0,4.,5.,6.,1.,1],
        [3,0,7.,2.,2.,1.,1],
        [4,0,6.,6.,6.,1.,2],
        [5,0,3.,7.,6.,1.,2],
        [6,0,5.,5.,9.,1.,5]]

# Dataframe from nodes data
node_tbl = pd.DataFrame(data,columns=['treenode_id','structure_id','x','y','z','r','parent_treenode_id'])

# Create a dummy neuron from node_tbl
neuron = trees.Tree.from_table(node_tbl)

# Resample the neuron with sampling length = 1 simply to get more points
neuron = neuron.resample(sampling_length=1)

# Plot neuron, only showing the node points
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
neuron.plot(ax,show_nodes=True,show_leaves=False,show_root=False,show_branchingnodes=False)
plt.tight_layout();

# Averaged local 3D scale of each node for sigma from 0 => 100
l3ds_tbl = neuron.compute_local_3d_scale_sigma(sig_lst=np.arange(100),postprocess='mean');

# Plot neuron whose nodes are colored by the local 3d scale
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
neuron.plot(ax,weights=l3ds_tbl,
            show_cbar=True,point_args={"cmap":"rainbow"},
            show_nodes=False,show_leaves=False,show_root=False,show_branchingnodes=False)
plt.tight_layout();

compute_local_3d_scale_radius(r_lst, dim_param=None, decomposed_method='leaf', rootid=None, postprocess='mean')

Compute local 3d scale of neuron from a list of radius of curvatures.

This is the advanced approach to compute the local 3d scale. Please see this [1] for more detail.

The tree is decomposed into segments. Then, local 3d scale is computed for each segment. Each segment is considered as a Curve object. Please see compute_local_3d_scale_radius() from Curve class for more detail. One node can have multiple local 3d scales due to the decomposition. You can specify the postprocess to return e.g. the mean local 3d scale of the node.

Parameters:

• r_lst (list of float) – list of radius of curvatures. The radius must be > 0.

• dim_param (dic) – parameters for dimension decomposition. If None, default parameters are used. See parameters of compute_local_3d_scale_radius() from Curve class for the meaning.

• decomposed_method (str) – “leaf”, “spine” or “branching”, method for decomposing the neuronal tree.

• rootid (int) – root node ID. If None, then the first rootid will be taken.

• postprocess (str) – support “mean”, “std”, “max”, “min” or None.

Returns:

1D, 2D, 3D flags, curvature and torsion.

Return type:

Pandas dataframe containing local 3D scale and additional information

Notes

• If postprocess = None, then a Pandas dataframe containing local 3d scale and additional features (curvature, torsion, etc) of each decomposed segment is returned. Otherwise, a Pandas serie containing only the local 3d scale of each node is returned.

• Running this method is much longer than compute_local_3d_scale_sigma().

References

[1]

Phan MS, Matho K, Beaurepaire E, Livet J, Chessel A. nAdder: A scale-space approach for the 3D analysis of neuronal traces. 2022. PLOS Computational Biology 18(7): e1010211. DOI: 10.1371/journal.pcbi.1010211

summary()

Return a brief summary of the tree.

Returns:

A Pandas Serie contains

• number of nodes

• number of leaves

• number of branching nodes

• number of connectors

to_curve(nodeid=None)

Convert a segment of tree given by the list of node ID to a Curve object.

Parameters:

nodeid (list of int) – List of node ID. If None, then consider all nodes of the tree.

Returns:

A Curve object.

to_points(nodeid=None)

Convert a segment of tree given by the list of node ID to a Points object.

Parameters:

nodeid (list of int) – List of node ID. If None, then consider all nodes of the tree.

Returns:

A Points object.

plot(ax, projection='3d', rootid=None, spine_only=False, show_root=True, show_leaves=True, show_branchingnodes=True, show_connectors=True, show_nodes=False, root_args={}, leaves_args={}, branchingnodes_args={}, connectors_args={}, weights=None, weights_display_type='c', point_args={}, show_cbar=False, cbar_args={}, line_args={}, scales=(1.0, 1.0, 1.0), equal_axis=True)

Plot tree using matplotlib.

Parameters:

• ax – plot axis.

• projection (str) – we support 3d, xy, xz, yz modes.

• spine_only (bool) – if True, then only plot tree spine.

• show_root (bool) – if True, then display root node.

• show_leaves (bool) – if True, then display leaves nodes.

• show_branchingnodes (bool) – if True, then display branching nodes.

• show_connectors (bool) – if True, then display connector nodes.

• show_nodes (bool) – if True, then display all nodes.

• root_args (dic) – plot params for root node.

• leaves_args (dic) – plot params for leaves nodes.

• branchingnodes_args (dic) – plot params for branching nodes.

• connectors_args (dic) – plot params for connectors nodes.

• weights (pandas serie) – serie indexed by node ID and contains corresponding weights

• weights_display_type (str) – “s” for point size or “c” for point color.

• point_args (dic) – point plot params.

• line_args (dic) – line plot params.

• scales (tuple of float) – use to set x, y and z scales.

• equal_axis (bool) – fix equal axes.