TubULAR#

Tube-like sUrface Lagrangian Analysis Resource (TubULAR) class

Coordinate Systems#

uv(conformal map onto unit square): Conformally mapping the cylinderCutMesh onto the unit square in the plane results in the instantaneous uv coordinate system. Corners of the unit square are taken directly from the cutMesh, so are liable to include some overall twist.
sphi(proper length x rectified azimuthal coordinate): Cylinder-like coordinate system in which first coordinate is the proper length along the surface and second is a rectified azimuthal coordinate. Rectification means that the surface is rotated at each longitudinal coordinate u so that the surface changes as little as possible over time. Mathematically, v is taken to phi(u), where each u coordinate is offset by a “rotation” about the surface along the v direction (the circumferential direction). This rectification may be based on surface positions in R^3 (geometric) or based on intensity motion in pullback space (material/Lagrangian) inferred through phase-correlation of the tissue at each u value (averaged across v).
sphi_sm(proper length x rectified azimuthal coordinate, smoothed in time): Same as sphi, but coordiantes smoothed over time with a tripulse filter.
APDV(rotated, translated, and scaled coordinate system, like lab frame): Physical coordinates used to visualize dynamics. Aligned meshes are in this 3D coordinate system.
xyz(data coordinate system): Pixel coordinates of the data itself. Raw meshes, cylinder meshes, cutMeshes etc are in this 3D coordinate system.

PIV measurement can be automatically computed in pullback space to determine 3D motion. These measurements can be done in a coordinate system of your choice (uv, sphi, or sphi_sm).

APDV coordinate system and Centerline specification#

QuapSlap allows the user to designate the AP axis, a DV axis, and a centerline for the segmented object. The default APDV coordinate system is automatically determined from the elongation axis of the surface and the data axes.

To customize, you can designate an APDV coordinate system that may be biologically relevant by using iLastik training on the timepoint t0 (which is an attribute of tubi, found as tubi.t0). Train in iLastik on anterior (A), posterior (P), background (B), and dorsal (D) location in different iLastik channels, making small blobs of high probability near the anterior end (A), somewhere along the posterior end for (P) so that AP forms the AP axis, and anywhere along dorsal for D. Then the centers of mass of thresholded contiguous high probability are computed for each to define the coordinate system. By default for this customized option, the ilastik results are read as ch1=A, ch2=P, ch3=B, ch4=D, but you may specify anteriorChannel, posteriorChannel, and dorsalChannel to specify the iLastik training channel that is used for each. Name the h5 file output from iLastik as …_Probabilities_APDVcoords.h5 For example, dorsal for the gut was chosen at the fused site where additional 48YGAL4-expressing muscle-like cells form a seam. Posterior is at the rear of the yolk, where the endoderm closes, for apical surface training. Anterior is at the junction of the midgut with the foregut. Separately, define the AP points for centerline extraction. For most gut data, the posterior point is different in centerline than it is for AP centerline specification, so we use a different ILP to train for A and P for all timepoints of dynamic data.

Properties#

xpImSAnE experiment class instance or struct with fields

expMeta : struct with fields fileMeta : struct with fields

dynamicbool

true if multiple timepoints, false if fixed

timeIntervalnumeric, default=1

increment in time between timepoints with indices differing by 1. For example, if timePoints are [0,1,2,4] and these are [1,1,2] minutes apart, then timeInterval is 1.

timeUnitsstr (default=’min’)

units of the timeInterval (ex ‘min’)

spaceUnitsstr (default = ‘$mu$m’)

units of the embedding space (ex ‘$mu$m’)

imSize2x1 int

size of pullback images to create (default is [a_ratio * 1000, 1000])

dirstruct with fields for different saved data names

directories where QuapSlap data lives

fileNamestruct with fields for different saved data names

Full filenames (with the absolute path) for various data and results

fileBasestruct with fields for different saved data names

Relative fileNames to be populated by timestamp (’…%06d…mat’)

fullFileBasestruct with fields for different saved data names

full path of filenames (like fullfile(tubi.dir.X, tubi.fileBase.X))

ssfactorint

subsampling factor for probabilities

APDVstruct with fields

‘resolution’, float: resolution of data in spaceUnits / pixel
‘rot’3x3 rotation matrix: rotation matrix to transform 3d pts from data frame into APDV frame accordint to v’ = (rot*v+trans)*resolution
‘trans’1x3 float array: translation vector to transform 3d pts from data frame into APDV frame according to v’ = (rot*v+trans)*resolution

flipybool

whether data is mirror image of lab frame coordinates

nVint (default=100)

sampling number along circumferential axis

nUint (default=100)

sampling number along longitudinal axis

uvextenstring, build during instantiation

naming extension with nU and nV like ‘_nU0100_nV0100’

t0numeric

reference time in the experiment, for building rectified sphi coordinates, also default for building pathlines

normalShiftnumeric (default=0)

shift to apply to meshes in pixel space along normal direction

a_fixednumeric (default=1)

aspect ratio for fixed-width pullbacks in uv and sphi coordinates

phiMethodstr specifier (‘3dcurves’ or ‘texture’)

method for stabilizing the v axis: either use the geometry of the surface to minimize surface deformation from timepoint to timepoint, so that v(t) for a given value of u is very near v(t-1) at the same value of u. This is a method for determining Phi map in pullback mesh creation, with the full map from embedding to pullback being [M’=(Phi)o(Z)o(M)], where M is a nearly-conformal mapping to the plane, Z changes the spacing of the longitudinal coordinates to reflect the average proper length traversed from the ‘ring’ of u=const to the ‘ring’ of u=const+du. This string specifier must be ‘3dcurves’ (geometric phi stabilization) or ‘texture’ (optical flow phi stabilization)

endcapOptionsstruct with fields

adistnumeric: distance around anterior point A which is removed from tubular mesh (sliced off)
pdistnumeric: distance around anterior point P which is removed from tubular mesh (sliced off)
trefnumeric: timestamp for reference time used to define th point on the endcap

at which we cut the cylinder mesh into a cylinderCutMesh (a topological disk/square). This “dorsal” point for other timepoints are identified by pointmatching.

Additional fields allowed

plottingstruct with fields

‘preview’bool (default=false): display intermediate results
‘save_ims’bool (default=true): save images along the way
‘xyzlim_um_buff’3x2 float: xyzlimits in um in RS coord sys with buffer
‘xyzlim_raw’3x2 float: xyzlimits in pixels
‘xyzlim_pix’3x2 float: xyzlimits in pixels, in rotated and translated coordinates
‘xyzlim_um’3x2 float: xyzlimits in um in rotated, translated, and scaled coordinate system
‘colors’Nx3 RGB values: color cycle for tubi

apdvPts = struct(‘anteriorPts’, [], …

‘posteriorPts’, [], … ‘antPts_sm’, [], … ‘postPts_sm’, [], … ‘dorsalPts’, [], … ‘antPts_rs’, [], … ‘postPts_rs’, [], … ‘dorsPts_rs’, [])

apdvOptionsoptional struct with fields

options used for finding APDV frame

currentTimenumeric

timestamp, assigned whenever time is set

currentMesh = struct with fields

‘rawMesh’struct with fields: original mesh found by surface detection
‘alignedMesh’: struct with fields: APDV rotated and scaled mesh (raw mesh in APDV coordinates)
‘cylinderMesh’, [],: original mesh with endcaps cut off
‘cylinderMeshClean’, [],: cylinder mesh with “ears” removed (ears give difficulty in mapping to the plane)
‘cutMesh’, [],: cylinder mesh with a seam given by cutPath
‘cutPath’, [],: vertex indices of the cutMesh along which the periodic seam is cut
‘uvcutMesh’, [],: rectilinear cutMesh in (u,v) from Dirichlet map result to rectangle
‘spcutMesh’, [],: rectilinear cutMesh in (s,phi) ‘surface Lagrangian’ parameterization
‘spcutMeshSm’, [],: rectilinear cutMesh in (s,phi) smoothed in time
‘spcutMeshSmRS’, [],: rectilinear cutMesh in (s,phi) smoothed in time with rotated scaled embedding
‘spcutMeshSmRSC’, [],: rectilinear cutMesh as closed cylinder (topological annulus), in (s,phi) smoothed, with rotated scaled embedding
‘ricciMesh’struct with fields: ricci flow result pullback mesh, topological annulus

currentCline = struct(‘mss’, [], …

‘mcline’, [], … ‘avgpts’, []) ;

data = struct(‘adjustlow’, 0, …

‘adjusthigh’, 0, … ‘axisOrder’, [1 2 3], … ‘ilastikOutputAxisOrder’, ‘cxyz’) options for scaling and transposing image intensity data

currentData = struct(‘IV’, [], …

‘adjustlow’, 0, … ‘adjusthigh’, 0 ) image intensity data in 3d and scaling

currentVelocitystruct with field

‘piv3d’ : struct with current timepoint’s PIV information in both 2d and 3d

pivstruct with fields for all timepoints’ PIV information

‘imCoords’str specifier (default=’sp_sme’): image coord system for measuring PIV / optical flow) ;
‘Lx’int: width of image, in pixels (x coordinate)
‘Ly’int: height of image, in pixels (y coordinate)
‘raw’struct: raw PIV results from disk/PIVLab
‘smoothed’ :: smoothed PIV results after gaussian blur
‘smoothing_sigma’ numeric (default=1 ): sigma of gaussian smoothing on PIV, in units of PIV sampling grid pixels

velocityAveragestruct with fields

vsmM(#timePoints-1) x (nX*nY) x 3 float array

3d velocities at PIV evaluation coordinates in um/dt rs

vfsmM(#timePoints-1) x (2*nU*(nV-1)) x 3 float array

3d velocities at face barycenters in um/dt rs

vnsmM(#timePoints-1) x (nX*nY) float array

normal velocity at PIV evaluation coordinates in um/dt rs

vvsmM(#timePoints-1) x (nU*nV) x 3 float array

3d velocities at (1x resolution) mesh vertices in um/min rs

v2dsmM(#timePoints-1) x (nX*nY) x 2 float array

2d velocities at PIV evaluation coordinates in pixels/ min

v2dsmMum(#timePoints-1) x (nX*nY) x 2 float array

2d velocities at PIV evaluation coordinates in scaled pix/min, but proportional to um/min (scaled by dilation of map)

‘v3d’: 3D velocities in embedding space [pix/dt]
‘v2d’: 2D tangential velocities in pullback
‘v2dum’: 2D tangential velocity scaled by speed in true embedding space
‘vn’: normal velocity in spaceUnits per timeInterval timeUnits
‘vf’: velocity vielf on face barycenters after Lagrangian avg
‘vv’, []) ;: velocity field on vertices after Lagrangian avg

cleanCntrlines :

centerlines in embedding space after temporal averaging

smoothingstruct with fields

‘lambda’float (default=0.00): diffusion const for field smoothing on mesh, if nonzero
‘lambda_mesh’float (default=0.00): diffusion const for vertex smoothing of mesh itself, if nonzero
‘nmodes’int (default=7): number of low freq modes to keep per DV hoop
‘zwidth’int (default=1): half-width of tripulse filter applied along zeta/z/s/u direction in pullback space, in units of du/dz/ds/dzeta

pathlinesstruct with fields

‘t0’numeric: time at which pathlines are rectilinear (from which time the points are advected). This is a timestamp (not an index) at which pathlines form regular grid in space.
‘refMesh’: reference mesh for pathline advection
‘piv’, []: Lagrangian pathlines from piv coords
‘vertices’, [],: Lagrangian pathlines from mesh vertices
‘vertices3d’,: Lagrangian pathlines from mesh vertices
‘faces’, [],: Lagrangian pathlines from mesh face barycenters
‘beltrami’, struct with fields for beltrami coefficient evaluated along pathlines: ‘mu_material’, [], … ‘mu_material_filtered’, [], … ‘mu_material_vertices’, [], … ‘fitlerOptions’, struct with fields describing filtering

currentStrainstruct with fields

‘pathline’ :: strain from pathlines
struct(‘t0Pathlines’, [], …: t=0 timepoint for pathlines in question
‘strain’ :: strain evaluated along pathlines
‘beltrami’ :: beltrami coefficient for pathlines

Full Contents:#

@TubULAR.measureXYZLims(tubi, options)#

[xyzlim_raw, xyzlim, xyzlim_um, xyzlim_um_buff] = measureXYZLims(tubi, options) Redo calculation of XYZ limits

Parameters:

tubi –
options (struct with fields) –

overwritebool
overwrite previous XYZlimits on disk

NPMitchell 2020

@TubULAR.measureSurfaceAreaVolume(QS, options)#

MEASURESURFACEAREAVOLUME(QS, OPTIONS): Compute the mesh surface area and volume over time

Prerequisites#

QS.alignMeshesAPDV()

param QS:

Current QuapSlap object

type QS:

QuapSlap class instance

param options:

overwritebool, default=false: overwrite results on disk
preivewbool, default=QS.plotting.preview: display intermediate results
thres_fracafloat, default=1.0: threshold fractional change in area above which we say the area and volume are NaN

type options:

struct with fields

returns:

<none>
NPMitchell 2020

@TubULAR.visualizeTracking3D(QS, options)#

Draw cell contours (if method == segmentation) or nuclear track points (if method == nuclei) and possibly connections between pairs of cells (if drawGeodesics == true) in 3D on top of texturePatch image.

Parameters:

options (struct with fields) –

timePoints : timepoints to visualize method : ‘segmentation’ or ‘nuclei’ subdir : ‘labeled_groundTruth’ selectPairs : int (default=0 for no pair selection)

number of cell pairs to select from segmentation for tracking. Selection occurs at t=t0.

drawGeodesics : bool geodesicStyle : ‘pullback’ or ‘geodesic’

’pullback’ draws a curve in 3d that connects the 2d segmentation pts as straight line in pullback space projected into 3d ‘geodesic’ computes geodesic pairs on the surface using CGAL

pairColorsselectPairs x 1 cell array: selectPairs x 1 cell array of colors (either string specifiers or 3x1 float array colors with values between 0-1.0) pairColors{ii} is the color for the pair path connecting tracked cells in pairIDs(ii, :)
viewAngles[-20,20] is default: must match texturePatch angles for overlays to be accurate

t0forPairs : timepoint at which we select which pairs to visualize overlayStyle : ‘mask’, ‘add’

whether to add or mask the segmentation/tracking/geodesics over the texturePatch image

lw_geoPfloat or int: linewidth for the geodesic paths
roiColorstring specifier for color or 3x1 float between 0-1.0: color for the rectangular ROI boundary
overwriteROIbool: if the ROI polygon (projected rectangle at t0) exists on disk, overwrite it with a newly calculated ROI

@TubULAR.measurePIV3d(tubi, options)#

measurePIV3D(tubi, options)

Measure 3d flows from PIV results on smoothed (s,phi) MIP pullbacks. Note that this code accounts for the time interval tubi.timeInterval in the dt between frames.

default options for created MIPS (done before this method) are:: pbOptions = struct() ; pbOptions.overwrite = false ; pbOptions.numLayers = [5, 5] ; pbOptions.layerSpacing = 0.75 ; pbOptions.generate_rsm = true ; pbOptions.generate_spsm = true ; pbOptions.generate_sphi = false ; tubi.data.adjustlow = 1.00 ; tubi.data.adjusthigh = 99.999 ; tubi.generateCurrentPullbacks([], [], [], pbOptions) ;

Parameters:

tubi (TubULAR class instance) –
options (struct with fields) –

overwritebool
overwrite previous results

previewbool
view intermediate results

includeMeshesInPullbackInspectionbool (default=false)
draw an advected mesh and reference mesh on top of pullback images to inspect

timePointsnumeric 1D array
the timepoints to consider for the measurement. For ex, could choose subset of the tubi experiment timePoints

Returns:

m0XY (Px2 float: 2d mesh vertices in pullback image pixel space [pullback pix]” ;)
m0f (#facesx3 float: mesh connectivity list” ;)
m0v3d (Px3 float: 3d mesh coordinates in embedding pixel space [mesh pix / dt]” ;)
x0 (QxR float: x value of 2d velocity evaluation coordinates in pullback image pixel space [pullback pix]” ;)
y0 (QxR float: y value of 2d velocity evaluation coordinates in pullback image pixel space [pullback pix]” ;)
pt0 (Nx3 float: 3d location of evaluation points [mesh pix]” ;)
pt1 (Nx3 float: 3d location of advected point in next mesh [mesh pix]” ;)
v0 (Nx3 float: 3d velocity [mesh pix / dt]” ;)
v0n (Nx1 float: normal velocity [mesh pix / dt]” ;)
v0t (Nx3 float: tangential velocity in 3d [mesh pix / dt]” ;)
facenormals (#faces x 3 float: face normals for all faces [unitless, mesh pix / mesh pix]” ;)
v3dfaces (#faces x 3 float: face-centered velocities in embedding space [mesh pix / dt]” ;)
v3dvertices (#mesh vertices x 3 float: vertex-centered velocities in embedding space [mesh pix / dt]” ;)
v0_rs (Nx3 float: rotated/scaled 3d velocities [um/min]” ;)
v0n_rs (Nx3 float: scaled normal velocities [um/min]” ;)
v0t_rs (Nx3 float: rotated/scaled tangential 3d velocities [um/min]” ;)
normals_rs (Nx3 float: normal vector of field faces in 3d [unitless, um/um]” ;)
v0t2d (Nx2 float: in-plane velocity [pullback pix / dt]” ;)
g_ab (Nx2x2 float: pullback metric” ;)
dilation (Nx1 float: dilation of face from 3d to 2d (A_2d [pullback_pix^2] / A_3d [mesh_pix^2])” ;)
jac (#faces x 1 cell: jacobian of 3d->2d transformation” ;)
fieldfaces (Nx1 int: face indices into spcutMesh where v defined” ;)
NPMitchell 2020

@TubULAR.sliceMeshEndcaps(tubi, opts, methodOpts)#

SLICEMESHENDCAPS() Create cut meshes with endcaps removed Note that a later step involves cutMesh, which slices along AP. Noah Mitchell 2019 This version relies on Gabriel Peyre’s toolbox called toolbox_fast_marching/

Parameters:

tubi (TubULAR class instance) – class instance for which we cut off endcaps of mesh
opts (struct with optional fields) – adist_thres pdist_thres adist_thres2 pdist_thres2 aCapMethod pCapMethod aOffset pOffset aOffset2 pOffset2 aOffsetRate pOffsetRate aOffset2Rate pOffset2Rate aDistRate pDistRate tref custom_aidx custom_pidx
methodOpts (struct with optional fields) –

trefint (timestamp, not index)
reference time stamp for dorsal definition, after and before which we point match to find the “dorsal” boundary vertex on each endcap

quickScanbool
first iterate through all points with large dt between frames to check that settings are good, then go back and compute each frame in order

overwrite save_figs preview timePoints
Prerequisites –
------------- –
pints (alignMeshesAPDV() after identifying APD) –
extractCenterlineSeries() –

First run extract_centerline.m or extractCenterlineSeries() before running this code. Run this code only after identifying anterior (A), posterior (P), and dorsal anterior (D) points via computeAPDpoints().

@TubULAR.alignMeshesAPDV(tubi, opts)#

ALIGNMESHESAPDV(opts) Uses anterior, posterior, and dorsal training in ilastik h5 output to align meshes along APDV coordinate system with global rotation matrix and translation vector. Extracted pts from the segmented training is loaded/saved in h5 file opts.rawapdvname (usually “apdv_pts_from_training.h5”). Smoothed pts from the segmented training –> opts.rawapdvname Smoothed rotated scaled pts –> opts.outapdvname Note that the global rotation matrix for the QuapSlap instance is defined previously in ptputeAPdpts()

This is a function similar to the align_meshes_APDV.m script

Parameters:

tubi (TubULAR class instance) –
opts (struct with fields) –
overwrite : bool overwrite_ims : bool smwindow : float or int

number of timepoints over which we smooth

normal_stepfloat
how far inside to push start/endpoints for centerline extraction (should be about 1 pixel to ensure centerlines can be found with some downsampling)

forceEndpointsInsidebool
push the start/ednpoints for centerline extraction further inside the mesh even if they are already a little bit inside the mesh

Returns:

xyzlim_raw – xyzlimits of raw meshes in units of full resolution pixels (ie not downsampled)
xyzlim – xyzlimits of rotated and translated meshes in units of full resolution pixels (ie not downsampled)
xyzlim_um – xyz limits of rotated and translated meshes in microns
xyzlim_um_buff – xyz limits of rotated and translated meshes in microns, with padding of QS.normalShift * resolution in every dimension
OUTPUTS
——-
xyzlim.txt – xyzlimits of raw meshes in units of full resolution pixels (ie not downsampled)
xyzlim_APDV.txt – xyzlimits of rotated and translated meshes in units of full resolution pixels (ie not downsampled)
xyzlim_APDV_um.txt – xyz limits of rotated and translated meshes in microns
rotation_APDV.txt – rotation matrix to align mesh to APDV frame, saved to fullfile(meshDir, ‘rotation_APDV.txt’) ;
translation_APDV.txt – translation vector to align mesh to APDV frame, saved to fullfile(meshDir, ‘translation_APDV.txt’)
xyzlim.txt – raw bounding box in original frame (not rotated), in full res pixels. Saved to fullfile(meshDir, ‘xyzlim.txt’)
xyzlim_APDV.txt – bounding box in rotated frame, in full resolution pixels. Saved to fullfile(meshDir, ‘xyzlim_APDV.txt’)
xyzlim_APDV_um.txt – bounding box in rotated frame, in microns. Saved to fullfile(meshDir, ‘xyzlim_APDV_um.txt’)
apdv_pts_rs.h5 (outapdvname) – Centers of mass for A, P, and D in microns in rotated, scaled APDV coord system. Note that this coord system is mirrored if flipy==true. Also contains raw apt,ppt,dpt in subsampled pixels. Saved to fullfile(meshDir, ‘centerline/apdv_pts_rs.h5’)
startendpt.h5 – Starting and ending points Saved to fullfile(meshDir, ‘centerline/startendpt_rs.h5’) ;
forceEndpointInside (bool) – push the endpoints for crude (fast marching) centerline extraction further inside the mesh by pushing along vertex normals by normal_step
NPMitchell 2020

@TubULAR.fitPhiOffsetsViaTexture(tubi, uspace_ds_umax, vspace, phi0_init, phi0TextureOpts)#

fit offsets in the DV direction of pullback mesh by measuring phase correlation between current and previous frame

Parameters:

phi0TextureOpts (struct with fields) –

lowerboundy : float lowerboundy : float step_phi0tile : float step_phi0tile : int (pixels)

How far apart each column to correlate, distance in pixels in x dir

width_phi0tilefloat: How wide in pixels to consider each circumferential strip’s correlation
potential_sigmayfloat: Impose quadratic potential cost for motion in phi away from input guess with this prefactor. This penalizes errant large deviations.

NPMitchell 2020

@TubULAR.computePCAoverTime(tubi, options)#

results = computePCAoverTime(tubi, options) Compute Principal Component Analysis of the sequence of deformations over time.

Parameters:

tubi (TubULAR class instance) –
options (struct with fields) –

overwritebool (default = false)
overwrite previous results on disk

overwriteImagesbool (default = false)
overwrite figures of previous results on disk

convert_to_periodbool (default=false)
convert all timestamps to period

Tnumeric (default=NaN)
number of timepoints per period if data is periodic in time

t0int (default=tubi.t0set())
timestamp to mark as t=0

NmodesToViewint (default= 5)
how many modes to plot individually on the reference surface configuration

nTimePoints2RmEndsint (default= 0)
how many timepoints to remove from the beginning and end of the timeseries. This is useful if the first and last few datapoints are noisy

drawArrowsInPCAPlanebool (default= false)
a plotting option, to draw arrows showing trajectory from timepoint to timepoint in 2D PCA projective planes

drawArrowsInPCA3Dbool (default=false)
a plotting option, to draw arrows showing trajectory from timepoint to timepoint in PCA space

meshChoicecell of strings or single string
(default={‘Lagrangian’, ‘sphi’}) the mesh style to use for computing PIV.

axPairsoptional int array (default=[1,2;1,3;2,3])
pairs of mode numbers to plot in 2d projective planes (modes are ordered by rank)

plotArrowsOnModesbool
plot arrows on the surfaces for mode decompositions

nArrowsint
subsampling factor for quiverplot on modes, only used if plotArrowsOnModes == true

displacement_scalefloat
scale factor to multiply pca velocities by before advecting mesh vertices by vn*v_pca –> ie displace vertices for visualization by vn*v_pca*displacement_scale.

displayTypes = {‘normal’, ‘displacement’, ‘rotational’, ‘dilatational’} ;
how to display the mode results (what field to plot on the surfaces)

useCamlightbool
use a camera light for mode images

Returns:

results

Return type:

cell of structs with fields

Saves to disk#

saves fullfile(tubi.dir.PCAoverTime, meshChoice, sprintf(‘pcaResults_%s.mat’, pcaType)) ;: where meshChoice is either ‘Lagrangian’ or ‘sphi’

saves figures to fullfile(tubi.dir.PCAoverTime, meshChoice)

NPMitchell 2022

@TubULAR.plotMetricKinematicsTimePoint(QS, tp, options)#

Plot the metric Kinematics, either for all timepoints or for a single timepoint. This function is called by measureMetricKinematics() TODO: also allow the function to be a standalone method to plot all timepoints if tp is not passed to options.

NPMitchell 2020

@TubULAR.generateCellSegmentation2D(tubi, options)#

segment cell membranes in 2d pullback projections NOTE: Dependency is Nick Noll’s tissueAnalysisSuite: this must be in the current Matlab path.

Parameters:

tubi –
options (struct with fields) – coordSys : coordinate system on which to compute the segmentation

tissueAnalysisSuite fields are different than QuapSlap’s. For tissueAnalysisSuite, we have: vdat:

nverts : neighbor list for vertices ncells : vertxcoord : column of data in which vertex lives vertycoord : row of data in which each vertex lives

cdatcell data: ncells : indices of neighboring cells nverts : vertices that define the cell centroid.coord(1) x position of cell centroid centroid.coord(2) y position of cell centroid
bdat :: nverts ncells pix : linear indices of the pixels associated with that bond

Note on exterior calculus objects d0 is an e x c matrix of exterior derivatives with +1 and -1s at the endpts of each bond d1 is a v x e matrix of exterior derivatives. Upstream is +1, downstream is -1 when moving counterclockwise around a tension plaquette. d0 and d1 are matrices that take derivatives

@TubULAR.plotPathlineStrainRate(tubi, options)#

Plot the strain rate along pathlines as kymographs and correlation plots. These are computed via finite differencing of pathline mesh metrics. That is, vertices are advected along piv flow and projected into 3D from pullback space to embedding, then the metrics g_{ij} of those meshes are measured and differenced.

Parameters:

tubi (QuapSlap class instance) –
options (struct with fields) – plot_kymographs : bool plot_kymographs_cumsum : bool plot_correlations : bool plot_fold_strainRate : bool plot_lobe_strainRate : bool plot_fold_strain : bool plot_lobe_strain : bool

Returns:

<none>
NPMitchell 2020

@TubULAR.measureDxDyStrainFiltered(QS, options)#

Compute strain along bonds that lie along dx and dy in pullback space from integrated pathlines deforming mesh vertices. Filter the results heavily in space. Measurements are taken with respect to fixed Lagrangian frame. Plot results in 2d and/or 3d for each timepoint.

Parameters:

QS (QuapSlap class instance) –
options (struct with fields) – overwrite : bool, overwrite data on disk overwriteImages : bool, overwrite images of results on disk plot_comparison : bool, plot a comparison with DEC traceful dilation median_filter_strainRates : bool
2020 (NPMitchell) –

@TubULAR.plotPathlineVelocities(tubi, options)#

Use pathlines of optical flow in pullback space to query velocities and average along pathlines.

Parameters:

tubi (TubULAR class instance) –
options (struct with fields) –

overwritebool, default=false
overwrite previous results

previewbool, default=false
view intermediate results

t0int, default=tubi.t0set()
timestamp at which the pathlines form a grid onto mesh vertices, mesh face barycenters, or PIV evaluation points

NPMitchell 2020

@TubULAR.plotSeriesOnSurfaceTexturePatch(tubi, options, TexturePatchOptions)#

PLOTDATAONSURFACETEXTUREPATCH(xp, xyzlim, rot, trans): Plot intensity data timeseries on 3d mesh timeseries

Parameters:

tubi (QuapSlap class instance) –
options (struct with fields) –

overwritebool
overwrite options on disk as .mat

texture_shiftfloat (optional, default=0)
shift of texture mesh, but not physical mesh
TexturePatchOptions (struct with fields) –
- Options.PSize: default=5
Special option, defines the image texture size for each individual polygon. A lower number gives a more pixelated texture {64}
- Options.Dilation: default=tubi.APDV.resolution
Special option, defines the units of space (pix2um, for ex)
- Options.Translation: default=tubi.APDV.trans
Special option, defines the translation vector applied after rotation on all surfaces
- Options.Rotation: default=tubi.APDV.rot
Special option, defines the dilation factor applied to all surfaces after rotation and translation
- Options.ApplyAmbientOcclusion: Determines if the texture
colors should be modified by the ambient occlusion of the underlying triangulation {‘false’}
- Options.AmbientOcclusion: #Vx1 list of ambient occlusion
values
- Options.AmbientOcculsionFactor: A scalar between [0,1]
determining how strongly the ambient occlusion modifies the texture mapping {1}
- Options.AmbientOcculsionSamples: The number of samples to use
when computing the ambient occlusion {1000}
- Options.Unoriented: Treat the surface as unoriented when
applying ambient occlusion {‘false’}
- Options.AddLights: Will add lighting to the plot {‘false’}
- Options.ScalarField: A vertex-based or face based scalar field
used to augment texture mapping with external color
- Options.ScalarCLim: The colormap limits used to display the
scalar field. By default it will use the full range of values in the scalar field
- Options.ScalarColorMap: The colormap used to display the scalar
field. Either a a string corresponding to a built-in MATLAB colormap or a user supplied #Cx3 colormap {‘parula’}
- Options.ScalarAlpha: A scalar between [0,1] determining
how the transparency of the overlay of the scalar field mapping {0.1};
- Options.VertexNormals: #Vx3 list of vertex unit normals
- Options.isRGB : bool, whether input arrays are R,G,B channels
- Options.isFalseColor : bool, whether to colorize with unique
color for each channel
- Options.falseColors : #channels x 3 list of colors for each
channel
- Options.Imax: float, maximum value for the data interpolant
object above which we clip the intensity
- Options.Imin: float, minimum value for the data interpolant
object below which we clip the intensity
- Options.extrapolationMethod: ‘nearest’ | ‘linear’ | ‘nearest’ | ‘next’ | ‘previous’ | ‘pchip’ | ‘cubic’ | ‘spline’ | ‘makima’ | ‘none’,
what extrapolation to use outside of the domain of data interpolation values
- Options.smoothIter : int, how many iterations of Laplcaian
smoothing to apply before normal displacement
- Options.numLayers : default=[1, 1]
length 2 array of ints, number of layers in positive, negative directions for MIP option
- Options.layerSpacing : float, default=2
distance between layers in texture space pixel coordinate units

Returns:

<none>
Saves to disk
————-
metadat.mat (mat file with variables ‘metadat’ and ‘Options’) – options for how to plot and compute texturepatches
NPMitchell 2020

@TubULAR.generateCurrentCutMesh(tubi, cutMeshOptions)#

Parameters:

tubi (TubULAR class instance) – The object for which we generate the currentTime’s cutMesh
cutMeshOptions (struct with fields) –

nsegs4path (int, optional, default=2)
How many segments of piecewise geodesics to draw to cut the axisymmetric surface so that the cut does not change winding number with respect to the previous timepoint’s cut around the centerline

maxJitter (float, optional, default=100)
maximum displacement in mesh units to randomly displace nodes of the path in a brute-force search for topologically connected path

maxTwChange (float, optional, default=0.15)
maximum allowed change in the value of the twist of the cutPath with respect to the centerline, as compared to the previous timepoint’s cutPath twist about the centerline. This is a soft proxy for the topology change, but in practice it is superior to the error-prone methods of measuring the (discrete-valued) winding number explored thus far

definePDviaRicci_t0(bool, default=false)
Compute the pullback coords at t0 via Ricci flow, then use Orbifold for other timepoints

useMaxDeviationPtsForCorrectionbool, default = true)
If the twist of a given cutpath is out of bounds, try different paths based on piecewise geodesics that have endpoints chosen based on their distance from the previous cutpath. This makes adjacent timepoints’ cutpaths similar to each other.
2020 (NPMitchell) –

@TubULAR.generateCellSegmentation3D(QS, options)#

In the code below, “aspect” typically means (a/b-1) = ars - 1

Define Q tensor for each cell by: Q = q (n^T n - II/2), where q = (sqrt(I_1/I_2) - 1) is the magnitude of the anisotropy
The Q tensor is related to the aspect ratio of cells via:: abs(norm(Q)) * 2 = |ar| - 1

The mratio

The factor of two arises on the LHS since the norm(n’*n - Identity*0.5) = 0.5 for any unit vector n, and an isotropic cell will have norm(Q) = 0.

NPMitchell 2021

@TubULAR.plotAlignedMeshesPretty(QS, options)#

Plot aligned APDV meshes with pretty lighting for timeseries

Parameters:

QS (QuapSlap class instance) –
options (struct with fields) –

overwritebool
overwrite previous images on disk

ApplyAmbientOcclusion: Determines if the texture
colors should be modified by the ambient occlusion of the underlying triangulation {‘false’}

AmbientOcclusion: #Vx1 list of ambient occlusion
values

AmbientOcculsionFactor: A scalar between [0,1]
determining how strongly the ambient occlusion modifies the texture mapping {1}

AmbientOcculsionSamples: The number of samples to use
when computing the ambient occlusion {1000}

Unoriented: Treat the surface as unoriented when
applying ambient occlusion {‘false’}

AddLights: Will add lighting to the plot {‘false’}

Returns:

<none>
Saves to disk
————-
metadat.mat (mat file with variables ‘metadat’) – options for how plotting was performed
NPMitchell 2021

@TubULAR.generateCurrentPullbacks(tubi, cutMesh, spcutMesh, spcutMeshSm, pbOptions)#

GENERATECURRENTPULLBACKS(tubi, cutMesh, spcutMesh, spcutMeshSm, pbOptions): Generate 2D images of tissue mapped from 3D. If save_as_stack is true, then we generate a 3D stack of 2D images along the normal direction of the surface.

Parameters:

pbOptions (struct) –
overwritebool (default=false)
overwrite existing images on disk

generate_sphibool (default=true)
create pullbacks in sphi coords

generate_relaxedbool (default=false)
create pullbacks in sphi coords, stretched along x to minimize areal distortion (but not anisotropic distortion). So this is not a conformal map.

generate_uvbool (default=false)
create pullbacks in uphi coords. This is a nearly conformal map.

generate_uphibool (default=false)
create pullbacks in uphi coords

generate_spsmbool (default=false)
create pullbacks in uphi coords

generate_rsmbool (default=false)
create pullbacks in uphi coords

generate_pivPathline(s)bool (default=false)
create pullback images from piv pathlines

PSizeint (default=5)
how many interpolation points along each side of a patch

EdgeColorcolorspec (default=’none’)
color of mesh edges in pullback image

YLim1x2 float (default=[0 1])
Y extent of pullback image in v or phi coords

axisorderlength 3 ints (default = [1 2 3 ])
axis permutation for the texture mesh

preTextureLambdafloat (default = 0)
If nonzero, apply laplacian smooth to mesh before rendering and before moving along normal_shift (which occurs before texture mapping)

Additional options as fields passed to texturePatch are
- pbOptions.imSize: The size of the output image
- pbOptions.baseSize: The side length in pixels of the smallest
  side of the output image when using the tight mesh bounding box
- pbOptions.xLim: The x-bounds of the output image
- pbOptions.yLim: The y-bounds of the output image
- pbOptions.pixelSearch: The method of searching for the faces
  
  containing pixel centers
  
  ’AABB’ (requires GPToolBox)
  
  ’Default’ (MATLAB built-ins, faster than AABB)
- pbOptions.numLayers: The number of onion layers to create
  Format is [ (num +), (num -) ]
- pbOptions.layerSpacing: The spacing between adjacent onion layers
  in units of pixels
- pbOptions.smoothIter: Number of iterations of Laplacian mesh
  smoothing to run on the mesh prior to vertex normal displacement (requires GPToolBox) (Default is 0)
- pbOptions.vertexNormal: User supplied vertex unit normals to the
  texture triangulation
- pbOptions.Interpolant: A pre-made texture image volume interpolant
2020 (NPMitchell) –

@TubULAR.prepareIlastik(tubi)#

Prepaere stack for Ilastik segmentation. This outputs an h5 file of subsampled intensity data at isotropic resolution on which to train.

prepareIlastik(stack)

@TubULAR.generateSPCutMeshSmStack(QS, spcutMeshSmStackOptions)#

Create TIFF stacks of normally evolved smoothed meshes in (s, phi) coordinate system.

Parameters:

QS (QuapSlap class instance) –
spcutMeshSmStackOptions (struct with fields) –
- n_outwardint
  number of steps in positive normal direction to sample
- n_inwardint
  number of steps in negative normal direction to sample
- overwritebool
  whether to overwrite spcutMeshSmStack on disk

Return type:

NPMitchell 2020

@TubULAR.plotStrainRate(tubi, options)#

plotStrainRateTimePoint(tubi, tp, options): Plot the traceful and traceless components of the strain rate tensor defined on each face, both for individual timepoints, and also over time as kymographs

Parameters:

tubi (TubULAR class instance) –
tp (int) – timepoint in units of (1/tubi.timeInterval) * tubi.timeUnits
options (struct with fields) –

Return type:

NPMitchell 2020

@TubULAR.timeAverageVelocities(tubi, options)#

measurePullbackStreamlines(tubi, options): Use pathlines of optical flow in pullback space to query velocities and average along Lagrangian pathlines. Default is to weight velocities contributions with a time width 5 weighted by tripulse filter. To avoid any averaging in time, set twidth = 0 ;

Parameters:

tubi (TubULAR class instance) –
options (struct with fields) –
twidth : average over (t-twidth, t+twidth) timepoints XYkernel : float plotOptions : optional struct with fields

vtscale : float vnscale : float

overwritebool
overwrite previous results

previewbool
view intermediate results
disk (Saves to) –
------------- –
vsmM ((#timePoints-1) x (nX*nY) x 3 float array) – 3d velocities at PIV evaluation coordinates in spaceUnits/dt rs
vfsmM ((#timePoints-1) x (2*nU*(nV-1)) x 3 float array) – 3d velocities at face barycenters in um/dt rs
vnsmM ((#timePoints-1) x (nX*nY) float array) – normal velocity at PIV evaluation coordinates in spaceUnits/dt rs
vvsmM ((#timePoints-1) x (nU*nV) x 3 float array) – 3d velocities at (1x resolution) mesh vertices in spaceUnits/dt rs
v2dsmM ((#timePoints-1) x (nX*nY) x 2 float array) – 2d velocities at PIV evaluation coordinates in pixels/dt
v2dsmMum ((#timePoints-1) x (nX*nY) x 2 float array) – 2d velocities at PIV evaluation coordinates in scaled pix/dt, but proportional to spaceUnits/dt (scaled by dilation of map)
2020 (NPMitchell) –

class @TubULAR.TubULAR(xp, opts)#

Bases: handle

Tube-like sUrface Lagrangian Analysis Resource (TubULAR) class. Note that TubULAR is fully documented, with docs accessible via:

doc TubULAR

Coordinate Systems#

uv(conformal map onto unit square, with one fixed vertex at each

endcap) Conformally mapping the cylinderCutMesh onto the unit square in the plane results in the instantaneous uv coordinate system. Corners of the unit square are taken directly from the cutMesh, so are liable to include some overall twist.

sphi(proper length x rectified azimuthal coordinate)

quasi-axisymmetric system in which first coordinate is proper length along surface and second is a rectified azimuthal coordinate. Rectification means that the surface is rotated as v->phi(s), where each s coordinate is offset by a “rotation” about the surface along a direction that is perpendicular to s in pullback space. This rectification may be based on surface positions in R^3 (geometric) or based on intensity motion in pullback space (material/Lagrangian) inferred through phasecorrelation of tissue strips around discretized s values.

r/spr/sphir(proper length x rectified azimuthal coordinate)

same as sphi but with aspect ratio relaxed to minimize isoareal energy cost –> makes triangles more similar in area by scaling the s axis (longitudinal axis) by a scalar factor.

iLastik’s coordinate systemsI recommend using cxyz as the axis

when outputting from iLastik. For our setup, this will mean that the meshes are mirrored with respect to the lab frame, so tubi.flipy = true.

PIV measurements:

‘piv’: principal surface-Lagrangian-frame PIV (sp_sme or up_sme)
–> note that the designation of coordinate system is not explicitly specified in the filenames: tubi.dir.mesh/gridCoords_nU0100_nV0100/piv/piv3d, etc

Properties#

imSize : 2x1 int, size of pullback images to create xyzlim : 3x2 float, mesh limits in full resolution pixels, in data space xyzlim_um : 3x2 float, mesh limits in lab APDV frame in microns resolution : float, resolution of pixels in um rot : 3x3 float, APDV rotation matrix trans : 3x1 float, APDV translation a_fixed : float, aspect ratio for fixed-width pullbacks phiMethod : str, ‘3dcurves’ or ‘texture’ flipy : bool, APDV coord system is mirrored XZ wrt raw data nV : int, sampling number along circumferential axis nU : int, sampling number along longitudinal axis uvexten : str, naming extension with nU and nV like ‘_nU0100_nV0100’ t0 : int, reference timePoint in the experiment velocityAverage : struct with fields

vsmM(#timePoints-1) x (nX*nY) x 3 float array
3d velocities at PIV evaluation coordinates in um/dt rs

vfsmM(#timePoints-1) x (2*nU*(nV-1)) x 3 float array
3d velocities at face barycenters in um/dt rs

vnsmM(#timePoints-1) x (nX*nY) float array
normal velocity at PIV evaluation coordinates in um/dt rs

vvsmM(#timePoints-1) x (nU*nV) x 3 float array
3d velocities at (1x resolution) mesh vertices in um/min rs

v2dsmM(#timePoints-1) x (nX*nY) x 2 float array
2d velocities at PIV evaluation coordinates in pixels/ min

v2dsmMum(#timePoints-1) x (nX*nY) x 2 float array
2d velocities at PIV evaluation coordinates in scaled pix/min, but proportional to um/min (scaled by dilation of map)

readme : struct with this information

APDV coordinate system and Centerline specification#

QuapSlap allows the user to designate the AP axis, a DV axis, and a centerline for the segmented object. To designate APDV coordinate system, use iLastik training on the timepoint t0 (which is an attribute of tubi). In particular, train on anterior (A), posterior (P), background (B), and dorsal (D) location in different iLastik channels. Small blobs of high probability near the anterior end for A, somewhere along the posterior end for P, and anywhere along dorsal for D are best. Then the centers of mass of thresholded contiguous high probability are computed for each to define the coordinate system. By default, the ilastik results are read as ch1=A, ch2=P, ch3=B, ch4=D, but anteriorChannel, posteriorChannel, and dorsalChannel specify the iLastik training channel that is used for each specification. Name the h5 file output from iLastik as …_Probabilities_APDVcoords.h5 For example, dorsal for the gut was chosen at the fused site where additional 48YGAL4-expressing muscle-like cells form a seam. Posterior is at the rear of the yolk, where the endoderm closes, for apical surface training. Anterior is at the junction of the midgut with the foregut. Separately, define the AP points for centerline extraction. For most gut data, the posterior point is different in centerline than it is for AP centerline specification, so we use a different ILP to train for A and P for all timepoints of dynamic data.

Property Summary

APDV#: translation vector to transform data into APDV frame (rot*v+trans)*resolution

a_fixed#: aspect ratio for fixed geometry pullback meshes

apdvOptions#: Options for computing the APDV frame

cleanFMCenterlines#: centerlines from fast marching method (crude centerlines) in embedding space, no temporal filtering

currentData#: image intensity data in 3d and scaling

currentMesh#: ricci flow result pullback mesh, topological annulus

currentSegmentation#: polygonal segmentation of cells in 3D space with manual corrections

currentStrain#: beltrami coefficient for pathlines

currentTime#: int, the timepoint of the current frame under examination (this can be set with setTime and changes during many method calls)

data#: options for scaling and transposing image intensity data

dir#: str, directory where QuapSlap data lives

dynamic#: expMeta : struct with fields fileMeta : struct with fields

endcapOptions#: the full map from embedding to pullback being [M’=(Phi)o()o()]. This string specifier must be ‘3dcurves’ (geometric phi stabilization) or ‘texture’ (optical flow phi stabilization)

features#: #timepoints x #folds float, positional proper length along surface of folds

fileBase#: fileNames to be populated by timestamp (’…%06d…mat’)

fileName#: fileName struct with fields such as rot, trans, etc where measurements/parameters are stored on disk

flipy#: whether data is mirror image of lab frame coordinates

fullFileBase#: full path of filenames (like fullfile(tubi.dir.X, tubi.fileBase.X))

imSize#: size of pullback images to create (default is [a_ratio * 1000, 1000])

nU#: sampling number along longitudinal axis

nV#: sampling number along circumferential axis

normalShift#: shift to apply to meshes in pixel space along normal direction

phiMethod#: method for determining Phi map in pullback mesh creation, with

piv#: sigma of gaussian smoothing on PIV, in units of PIV sampling grid pixels

plotting#

adist : distance around anterior point A which is removed from tubular mesh (sliced off) pdist : distance around anterior point P which is removed from tubular mesh (sliced off) tref : reference time used to define th point on the endcap

at which we cut the cylinder mesh into a cylinderCutMesh (a topological disk/square). This “dorsal” point for other timepoints are identified by pointmatching.

Additional fields allowed :

smoothing#: half-width of tripulse filter applied along zeta/z/s/u direction in pullback space, in units of du/dz/ds/dzeta

spaceUnits#: units of the embedding space (ex ‘$mu$m’)

ssfactor#: subsampling factor for probabilities

t0#: reference time in the experiment

timeInterval#: increment in time between timepoints with

timeStampStringSpec#: String format specifier for the timestamp. For ex, Time_001.tif would be %03d since filename is returned by sprintf(‘Time_%03d.tif’, 3)

timeUnits#: indices differing by 1. For example, if timePoints are [0,1,2,4] and these are [1,1,2] minutes apart, then timeInterval is 1.

uvexten#: naming extension with nU and nV like ‘_nU0100_nV0100’

velocityAverage#: velocity field on vertices after Lagrangian avg

velocityRaw#: velocity field on vertices, no temporal filtering

xp#: struct with fields or ImSAnE experiment class instance

Method Summary

APDV2dxyz(tubi, a)#: ars = xyz2APDV(tubi, a) Transform 3d vectors from APDV coord sys to XYZ data space

APDV2xyz(tubi, a)#: ars = xyz2APDV(tubi, a) Transform 3d coords from APDV coord sys to XYZ data space

static XY2uv(im, XY, doubleCovered, umax, vmax)#

Map pixel positions (1, sizeImX) and (1, sizeImY) to (0, umax) and (0, vmax) of pullback space if singleCover, or y coords are mapped to (-0.5, 1.5)*vmax if doubleCover

NOTE THAT MAP IS [xesz, yesz] = [size(im, 1), size(im, 2)] uv(:, 1) = umax * (XY(:, 1) - 1) / xesz ; uv(:, 2) = vmax * 2.0 * (XY(:, 2) - 1) / yesz - 0.5 ;

NOTE THAT INVERSE MAP IS x–> (xy(:, 1) * (Xsz-1)) / (1*umax) + 1 , … y–> (xy(:, 2) * (Ysz-1)) / (2*vmax) + 1 + (Ysz-1)*0.25 ;

Parameters:

im (NxM numeric array or 2 ints for size) – image in which pixel coordinates are defined or dimensions of the image (pullback image in pixels)
XY (Qx2 numeric array) – pixel coordinates to convert to pullback space
doubleCovered (bool) – the image is a double cover of the pullback (extended/tiled so that the “top” half repeats below the bottom and the “bottom” half repeats above the top. That is, consider im to be a double cover in Y (periodic in Y and covers pullback space twice (-0.5 * Ly, 1.5 * Ly)
umax (float) – extent of pullback mesh coordinates in u direction (X)
vmax (float) – extent of pullback mesh coordinates in v direction (Y) before double covering/tiling

Returns:

uv – pullback coordinates of input pixel coordinates

Return type:

Qx2 numeric array

adjustIV(tubi, IV, adjustlow, adjusthigh, forceValues)#

Parameters:

tubi (QuapSlap class instance) –
IV (3D numeric (optional, loads currentData if empty)) – data to adjust
adjustlow (numeric or #channels x 1 numeric) – minimum intensity or percent intensity to adjust data
adjusthigh (numeric or #channels x 1 numeric) – maximum intensity or percent intensity to adjust data
forceValues (enforce the adjustlow/high to be intensity) –
values –
percentile (not) –
100. (even if <) –

clearTime(tubi)#: clear current timepoint’s data for tubi instance

static clipXY(xx, yy, Lx, Ly)#

Clip x at (1, Lx) and clip Y as periodic (1=Ly, Ly=1), for image that is periodic in Y. Consider Y in [1, Ly]. If the pullback is a doubleCover, Ly = 2*mesh width in pixels

Parameters:

xx –
yy –
Lx (int) – number of pixels along x dimension
Ly (int) – number of pixels along y dimension

Returns:

[xx, yy]

Return type:

x and y coordinates clipped to [1, Lx] and [1, Ly]

static dXY2duv(im, dXY, doubleCovered, umax, vmax)#

XY2uv(im, XY, doubleCovered, umax, vmax): Map difference in pixel positions defined on (1, sizeImX) and (1, sizeImY) to (0, umax) and (0, vmax) of pullback space if singleCover, or y coords are mapped to (-0.5, 1.5)*vmax if doubleCover

NOTE THAT FULL COORDINATE MAP IS [xesz, yesz] = [size(im, 1), size(im, 2)] uv(:, 1) = umax * (XY(:, 1) - 1) / xesz ; uv(:, 2) = vmax * 2.0 * (XY(:, 2) - 1) / yesz - 0.5 ;

NOTE THAT INVERSE MAP IS x–> (xy(:, 1) * (Xsz-1)) / (1*umax) + 1 , … y–> (xy(:, 2) * (Ysz-1)) / (2*vmax) + 1 + (Ysz-1)*0.25 ;

Parameters:

im (NxM numeric array or 2 ints for size) – image in which pixel coordinates are defined or dimensions of the image (pullback image in pixels)
dXY (Qx2 numeric array) – difference in pixel coordinates as vector, to convert to pullback space
doubleCovered (bool) – the image is a double cover of the pullback (extended/tiled so that the “top” half repeats below the bottom and the “bottom” half repeats above the top. That is, consider im to be a double cover in Y (periodic in Y and covers pullback space twice (-0.5 * Ly, 1.5 * Ly)
umax (float) – extent of pullback mesh coordinates in u direction (X)
vmax (float) – extent of pullback mesh coordinates in v direction (Y) before double covering/tiling

Returns:

uv – pullback coordinates of input pixel coordinates

Return type:

Qx2 numeric array

doubleCoverPullbackImages(tubi, options)#

options : struct with fields coordsys : (‘sp’, ‘uv’, ‘up’)

coordinate system to make double cover

overwritebool, default=false: whether to overwrite current images on disk
histeqbool, default=true: perform histogram equilization during pullback extension
ntilesint, default=50: The number of bins in each dimension for histogram equilization for a square original image. That is, the extended image will have (a_fixed * ntiles, 2 * ntiles) bins in (x,y).
a_fixedfloat, default=tubi.a_fixed: The aspect ratio of the pullback image: Lx / Ly

static doubleToSingleCover(XY, Ly)#: detect if XY is passed as a pair of grids

static dvAverageNematic(magnitude, theta)#

[mag_ap, theta_ap] = DVAVERAGENEMATIC(magnitude, theta) Given a nematic field defined on a rectilinear grid, with Q(i,j) being in the ith ap position and jth dv position, average the nematic field over the dv positions (dimension 2)

Parameters:

magnitude (nU x nV numeric array) – magnitude of nematic field in 2D rectilinear grid
theta (nU x nV float array, with values as angles in radians) – angle (mod pi) of nematic field in 2D rectilinear grid

Returns:

mag_ap (nU x 1 float array) – average magnitude of nematic field along the ap dimension
theta_ap (nU x 1 float array) – average angle (mod pi) of nematic field along the ap dimension

dx2APDV(tubi, da)#: dars = dx2APDV(tubi, da) Transform 3d difference vector from XYZ data space to APDV coord sys

getAPpointsSm(tubi)#: Load the anterior and posterior ‘centers of mass’ ie the endpoints of the object’s centerline

getCurrentClineDVhoop(tubi)#: [mss, mcline, avgpts] = getCurrentClineDVhoop(tubi) Compute or load the centerline based on centroids of each circumferential hoop of constant s (mean pathlength along tube)

getCurrentCutMesh(tubi)#: Load/recall the current cutMesh with pullback uv coordinates. This mesh is not a rectilinear mesh topology with vertices in an nU x nV grid, but instead has the topology of the mesh found from getMeshes after cutting.

getCurrentPathlineStrain(tubi, t0Pathlines, varargin)#

Parameters:

tubi (current class instance) –
t0Pathlines (int or empty) – reference timepoint. If empty, set to tubi.t0set()
varargin (strings ('strain', 'beltrami')) – which strain measures to load from disk and attribute to self

Returns:

strain

Return type:

tubi.currentStrain.pathline

getCurrentSegmentation2D(tubi, options)#

Obtain the cell segmentation in 3D pullback space

Parameters:

tubi (current TubULAR class instance) –
options (struct with fields passed to) –
or (generateCellSegmentation() if segmentation is not on disk) –
properties (in current) –

getCurrentSegmentation2DCorrected(tubi, options)#: Decide on coordinate system for corrected pullback binary

getCurrentSegmentation3D(tubi, options)#: Obtain the cell segmentation in 3D pushforward space

getCurrentSegmentation3DCorrected(tubi, options)#: Obtain the cell segmentation in 3D pullback space Note: this is assumed to be coordinate-system independent (ie we do not store what 2D coordsys was used to generate this 3D segmentation.

getCurrentVelocity(tubi, varargin)#

[piv3d, vsm] = getCurrentVelocity(tubi, varargin)

Parameters:

tubi (tubular class instance) –
varargin (cell of strings 'piv3d' and/or 'average') – which velocities to load

Returns:

piv3d (struct of raw, piv-based velocities)
vsm (Lagrangian-frame time-averaged velocities for all timepoints, with)
the averaging done over a chosen/default time window in short
material pathlines.

getFeatures(tubi, varargin)#: GETFEATURES(QS, varargin) Load features of the tubular object (those specied, or all of them). Currently, features include {‘folds’, ‘fold_onset’, ‘ssmax’, ‘ssfold’, ‘rssmax’, ‘rssfold’}.

getPIV(tubi, options)#: Load PIV results and store in tubi.piv if not already loaded

getPathlineStrain(tubi, options)#

[strain, tre, dev, theta] = getPathlineStrain(tubi, options) Compute or load the integrated tissue strain measured in the Lagrangian frame via the deformation of a mesh whose vertices are advected along pathlines.

Parameters:

tubi (TubULAR instance) –
options (optional struct with fields) – t0Pathlines : int (default = t0)

Returns:

strain

Return type:

#tps x nFaces x 2 x 2 float array

getPullbackPathlines(tubi, t0, varargin)#: Discern if we must load pathlines or if already loaded

getRadii(tubi, options)#

Load radii from disk in specified coordinate system

Parameters:

options (struct with fields) –

coordSys: str specifier for how to return radii: ’spcutMeshSmRSC’ –> return nU x (nV-1) array for radii of vertices in spcutMeshSmRSC embedding <add other options here>

getRotTrans(tubi)#: Load the translation to put anterior to origin and AP axis along x axis

getVelocityAverage(tubi, varargin)#: todo: check if all varargin are already loaded

getVelocityRaw(tubi, varargin)#: todo: check if all varargin are already loaded

getVelocitySimpleAverage(tubi, varargin)#: todo: check if all varargin are already loaded

getXYZLims(tubi)#: [raw, pix, um, um_buff] = GETXYZLIMS(tubi) Grab each xyzlim from self, otherwise load from disk full resolution pix

static invertRotation(rot)#: Obtain rotation matrix that undoes the APDV rotation

loadCurrentCutMesh(tubi)#: Load the current cutMesh with pullback uv coordinates. This mesh is not a rectilinear mesh topology with vertices in an nU x nV grid, but instead has the topology of the mesh found from getMeshes after cutting.

loadFeatures(tubi, varargin)#

Load all features stored in QS.features

Parameters:: varargin (optional string list/cell) – which specific features to load

loadPIV(tubi, options)#: Load PIV results from disk and store in tubi.piv

loadPullbackPathlines(tubi, t0, varargin)#

loadVelocityAverage(tubi, varargin)#: Retreive the Velocity on each vertex/face barycenter of the (s,phi) meshes. NOTE: that these are nearly Lagrangian for most applications, but not perfectly Lagrangian. See also: tubi.getPullbackPathlines(t0), which corrects for any residual motion in the pullback plane.

loadVelocityRaw(tubi, varargin)#: Load and pack into struct

loadVelocitySimpleAverage(tubi, varargin)#: Load and pack into struct

makeMIPs(tubi, dim, pages, timePoints, adjustIV)#

Generate maximum intensity projections along the desired axis and for all pages specified.

Parameters:

tubi (TubULAR class instance) –
dim (int) – dimension along which to project
pages (int array) – pages of the tiff stack to project
timePoints (int array) – timePoints for which to take MIPs
adjustIV (bool (default=true)) – apply the intensity adjustment stored in tubi.data.adjustlow/high

manualTrackingAdd(QS, options)#: Add to current tracks on disk

measureLength(tubi, options)#: Measure the length of the centerline curve over time

measureRMSvelocityOverTime(tubi, options)#

Parameters:

options (optional struct with fields) – weights : weights to attribute to each face in mesh coordSys : str specifier (default==’sp’)

Returns:

vrms (#tps x 1 float array) – rms velocity over time
timestamps (#tps x 1 float array) – timestamps associated with each entry

mip(tubi, tp, dim, pages, adjustIV)#: Take maximum intensity projection

plotPathlineBeltramiKymograph(tubi, t0Pathlines, options)#: Example usage for 2021 gut paper: options = struct(‘ylim’, [0, 2] ) tubi.plotPathlineBeltramiKymograph([], options)

static quarterIndicesDV(nV)#: [dorsal, ventral, left, right] = quarterIndicesDV(nV) indices for each quarter of a DV section in grid coordinates

samplePullbackPathlines(tubi, XY0, options)#: [p2d, p3d] = samplePullbackPathlines(tubi, XY0, options) start at XY0, folow flow using barycentric coordinates of PIV pullback pathlines. This is the same as advecting along a fixed Langrangian coordinate.

setDataLimits(tubi, tp, adjustlow_pctile, adjusthigh_pctile)#: Use timepoint (tp) to obtain hard values for intensity limits so that data is rescaled to fixed limits instead of percentile. This is useful to avoid flickering of overall intensity in data in which a few voxels vary a lot in intensity.

setTime(tubi, tt)#

Set the current time of the dataset and clear current data which was associated with the previously considered time

Parameters:: tt (int or float) – timePoint to set to be current, from available times in tubi.xp.fileMeta.timePoints

t0set(tubi, t0)#: t0set(tubi, t0) Set time offset from file or manually. If t0 is saved to disk, it should have a header like: # reference time for dataset 76 <eof> and should be saved at tubi.fileName.t0

uv2APDV(tubi, uv, coordSys, umax, vmax)#: Convert (u,v) pullback coordinates to xyz coordinates in the APDV frame (rotated and scaled frame aligned with AP axis along x and DV axis along z)

static uv2XY(im, uv, doubleCovered, umax, vmax)#

XY = uv2XY(im, uv, doubleCovered, umax, vmax): Map from pullback uv u=(0,1), v=(0,1) to pixel XY

x–> (xy(:, 1) * (size(im, 2)-1)) / (1*umax) + 1 , … y–> (xy(:, 2) * (size(im, 1)-1)) / (2*vmax) + 0.75 + (size(im,1)-1)*0.25

Parameters:

im (NxM numeric array or length(2) int array) – 2D image into whose pixel space to map or size(im)
uv (Q*2 numeric array) – mesh coordinates to convert to pullback pixel space (XY)
doubleCovered (bool) – the image is a double cover of the pullback (extended/tiled so that the “top” half repeats below the bottom and the “bottom” half repeats above the top. That is, consider im to be a double cover in Y (periodic in Y and covers pullback space twice (-0.5 * Ly, 1.5 * Ly)
umax (float) – extent of pullback mesh coordinates in u direction (X)
vmax (float) – extent of pullback mesh coordinates in v direction (Y) before double covering/tiling

Returns:

XY – positions of uv coordinates in pullback pixel space

Return type:

N x 2 float array

xyz2APDV(tubi, a)#: ars = xyz2APDV(tubi, a) Transform 3d coords from XYZ data space to APDV coord sys

@TubULAR.plotAverageVelocitiesTimePoint(QS, tp, options)#

plotTimeAvgVelocitiesTimePoint(QS, tp, options)

Parameters:

QS (QuapSlap class instance) –
options (optional struct with fields) –
vnsm : v2dsmum : overwrite : bool

overwrite previous results

previewbool
view intermediate results

timePointsnumeric 1D array
the timepoints to consider for the measurement. For ex, could choose subset of the QS experiment timePoints

alphaValfloat
the opacity of the heatmap to overlay

invertImagebool
invert the data pullback under the velocity map

NPMitchell 2020

@TubULAR.measurePIV2d(tubi, options)#

Measure 2d PIV on pullbacks in some coordinate system (coordSys), either using PIVLab (default is use_PIVLab=true) or using a simple phasecorrelation with a single pass at a single fixed window size.

Note that this code does NOT use the time interval tubi.timeInterval between frames, since the 2D PIV connects the optical flow of adjacent frames and does not measure velocity in any physical sense.

Parameters:

tubi (TubULAR class instance) –
options (struct with optional fields) –
use_PIVLab : bool (default=True) edgeLength : size of interrogation window in pixels histequilize : bool (default=True)

perform histogram equilization before piv calculation

timePointsNx1 numeric (default is tubi.xp.fileMeta.timePoints)
the timepoints for which to compute PIV

coordSysstring specifier (default is tubi.piv.imCoords)
the image coordinates in which PIV is to be performed

Returns:

none
Saves to disk
————-
tubi.fileName.pivRaw.raw (‘x’, ‘y’, ‘u’, ‘v’, ‘u_filtered’, ‘v_filtered’)
NPMitchell 2022

@TubULAR.measureTwist(tubi, options)#

Measure twist of flow field from pullback: dv_phi / dzeta, where v_phi: is the angular velocity (in radians / timeUnits, NOT in spaceUnits / timeUnits). This is obtained by simply taking the derivative of the 2d velocity field in the (s,phi) coordinate system wrt the longitudinal coordinate.
We might consider projecting onto the instantaneous direction: perpendicular to the longitudinal axis to disambiguate orthogonal contributions, but we have not done this as of 2023.

Parameters:

tubi (TubULAR class instance) –
options (struct with fields) –

overwritebool
overwrite previous results

previewbool
view intermediate results

writheStylestr
style of writhe computation to compare to twist
2020-2023 (NPMitchell) –

@TubULAR.plotClineXSections(QS, options)#

Plot fancy “cross-section” view of centerlines, saved to DVhoop centerline directory

NPMitchell 2021

@TubULAR.measureWrithe(tubi, options)#

MEASUREWRITHE(TUBI, OPTIONS): Compute the centerline length and writhe and save to disk

Parameters:

tubi (TubULAR class instance) – Current TubULAR object
options (struct with fields) –

overwritebool, default=false
overwrite results on disk

overwriteFigsbool, default=false
overwrite figure results on disk

preivewbool, default=QS.plotting.preview
display intermediate results

Wr_stylestr specifier, default=’Levitt’ ;
which style of writhe calculation to plot

black_figsbool
plot figures with black background rather than white

flipy_centerlinebool
flip the centerline data y –> -y (inverts sign of writhe)

Returns:

<none>
NPMitchell 2020

@TubULAR.computeRicciMeshes(QS, options)#

computeRicciMeshes(QS, options): Compute all Ricci meshes (if flow converges) and plot aspect ratio for isothermal PB over time

Parameters:: options (optional struct with fields) – resample : bool (default=true) ;

NPMitchell 2021

@TubULAR.plotMetricKinematics(tubi, options)#

Plot the metric Kinematics as kymographs and correlation plots Out-of-plane motion is v_n * 2H, where v_n is normal velocity and H is mean curvature. In-plane motion considered here is div(v_t) where v_t is tangential velocity on the curved surface. The difference div(v_t) - vn*2H = 1/2 Tr[g^{-1} dot{g}], which is a measure of isotropic metric change over time (dot is the partial derivative wrt time).

Parameters:

tubi (TubULAR class instance) –
options (struct with fields) – plot_kymographs : bool plot_kymographs_cumsum : bool plot_correlations : bool plot_gdot_correlations : bool plot_gdot_decomp : bool plot_Hgdot : bool (default = false)
2020 (NPMitchell) –

@TubULAR.generateCellSegmentationPathlines3D(tubi, options)#

generateCellSegmentation3D(tubi, options)

Load segmentation results from a timepoint t0Pathlines, advect the cell polygons along pathlines from PIV, measure the segmentation properties of the advected “tissue” pattern frozen in the Lagrangian frame as the material deforms.

NPMitchell 2021

@TubULAR.plotTimeAvgVelocities(QS, options)#

plotTimeAvgVelSimple(QS, options): Save images of the velocity field over time that has been “simply” averaged over time in-place in the surface-Lagrangian pullback. That is, the velocity field at location (u,v) has been averaged in time with previous and later timepoints at the same pullback location (u,v). For samplingResolution, all fields plotted in this method depend only on the PIV sampling settings

Parameters:

QS (QuapSlap class instance) –
options (struct with fields) –

overwritebool
overwrite previous results

previewbool
view intermediate results

timePointsnumeric 1D array
the timepoints to consider for the measurement. For ex, could choose subset of the QS experiment timePoints

alphaValfloat
the opacity of the heatmap to overlay

invertImagebool
invert the data pullback under the velocity map

vtscale = 5 ; % um / min vnscale = 2 ; % um / min vscale = 2 ; % um / min alphaVal = 0.7 ; % alpha for normal velocity heatmap washout2d = 0.5 ; % lightening factor for data qsubsample = 5 ; % quiver subsampling in pullback space

NPMitchell 2020

@TubULAR.measureStrainRate(tubi, options)#

Compute epsilon = 1/2 (nabla_i v_j + nabla_j v_i) - vN b_{ij} The result on faces is smoothed insofar as the velocities are smoothed with options.lambda and the mesh (governing b_{ij}) is smoothed with options.labmda_mesh. The result on vertices is additionally smoothed via a quasi-1D spectral filter governed by options.nmodes and options.zwidth.

Parameters:

tubi (TubULAR class object instance) –
options (struct with fields) – lambda : float lambda_mesh : float overwrite : bool overwriteImages : bool preview : bool averagingStyle : str specifier (‘Lagrangian’ or ‘Simple’) samplingResolution : ‘1x’ or ‘2x’ debug : bool

Returns:

(none)
Saves to disk
————-
strainRate results in QS.dir.strainRate.measurements/strainRate_%06d.mat – -> one result per timepoint, with spectral smoothing applied to vertex-based results
NPMitchell 2020

@TubULAR.generateRicciMeshTimePoint(QS, tp, options)#

Compute solution to Ricci flow for each timepoint. If there is only one timepoint, we by default flow the spcutMesh. Otherwise we by default flow the spcutMeshSm.

Parameters:

QS (quapSlap class instance) –
tp (int (default=QS.t0)) – timestamp of timepoint for which to generate Ricci flow mesh
options (optional struct with optional fields) –

resamplebool (default = true)
perform isotropic resampling of the mesh before performing Ricci flow to ensure vertex quality during operation

maxIter : int (default=100) radiusTolerance : float (default=0.01)

maximum allowed fractional deviation of Ricci flow solution inner and outer radius from a true circle with fixed radius (variable inner radius, fixed outer radius of 1)

save_imsbool (default = true)
save images of the ricci flow solution

pathline_computationbool
compute the ricci flow for the pathline mesh of this timepoint

t0Pathlines
only used if pathline_computation, which pathline to look onto for this timpoint (for example load the mesh for t=50 advected along pathlines from t0=1).

Returns:

Saves to disk
————-
Ricci flow solution – ricciFn = sprintf(QS.fullFileBase.ricciSolution, maxIter, tp)
ricciMesh has fields annulus and rectangle – ricciMeshFn = sprintf(QS.fullFileBase.ricciMesh, maxIter, tp) ;
Beltrami coeff mu for this #iterations – mufn = sprintf(QS.fullFileBase.ricciMu, maxIter, tp) ;
Images saved
————
Plot Ricci result in annulus – fullfile(imDir, sprintf(‘%06d_RicciSolution.png’, tp)) ;
plot beltrami for Ricci flow – fullfile(imDir, sprintf(‘%06d_RicciFlowBeltrami.png’, tp)) ;
Plot just the bare triangulation – fullfile(imDir, sprintf(‘%06d_RicciFlowSolution.png’, tp)) ;
Histogram |mu| for each case – fullfile(imDir, sprintf(‘%06d_BeltramiCoefficients.png’, tp)) ;
Image of corrected vertices on inner and outer annulus – fullfile(imDir, sprintf(‘%06d_ricci_InnerCorrection.png’, tp)) ; fullfile(imDir, sprintf(‘%06d_ricci_OuterCorrection.png’, tp)) ;
Orientation and branch cut – fullfile(imDir, sprintf(‘%06d_DrhoDphi_flipped.png’, tp)) ; fullfile(imDir, sprintf(‘%06d_phiOrderInitial.png’, tp)) ; fullfile(imDir, sprintf(‘%06d_phiOrderFinal.png’, tp)) ;

NPMitchell 2020

@TubULAR.visualizeDemoTracks(QS, Options)#

Load demo tracks (segmented or tracked objects) and plot over pullback data.

Parameters:

QS –
Options (optional struct with fields) –

Returns:

<none>
Saves to disk
————-
image sequence of RGB overlays
NPMitchell 2021

@TubULAR.measurePathlineMetricKinematics(tubi, options)#

Query the metric Kinematics along lagrangian pathlines. Plot results as kymographs and correlation plots. Out-of-plane motion is v_n * 2H, where v_n is normal velocity and H is mean curvature. In-plane motion considered here is div(v_t) where v_t is tangential velocity on the curved surface. The difference div(v_t) - vn*2H = Tr[g^{-1} dot{g}], which is a measure of isotropic metric change over time (dot is the partial derivative wrt time).

Parameters:

tubi (TubULAR class instance) –
options (struct with fields) – plot_kymographs : bool plot_kymographs_cumsum : bool plot_correlations : bool plot_gdot_correlations : bool plot_gdot_decomp : bool
2020 (NPMitchell) –

@TubULAR.processCorrectedCellSegmentation2D(QS, options)#

Unfinished code – Lin working on it. ToDo:

tissueAnalysisSuite fields are different than QuapSlap’s. For tissueAnalysisSuite, we have: Vdat:

nverts : neighbor list for vertices ncells : vertxcoord : column of data in which vertex lives vertycoord : row of data in which each vertex lives

Cdatcell data: ncells : indices of neighboring cells nverts : vertices that define the cell centroid.coord(1) x position of cell centroid centroid.coord(2) y position of cell centroid
Bdat :: nverts ncells pix : linear indices of the pixels associated with that bond

@TubULAR.generateRawRicciMeshTimePoint(QS, tp, options)#

Compute solution to Ricci flow for a given timepoint’s cylinderMesh.

Parameters:

QS (quapSlap class instance) –
tp (int (default=QS.t0)) – timestamp of timepoint for which to generate Ricci flow mesh
options (optional struct with optional fields) –
maxIter : int (default=100) radiusTolerance : float (default=0.01)

maximum allowed fractional deviation of Ricci flow solution inner and outer radius from a true circle with fixed radius (variable inner radius, fixed outer radius of 1)

save_imsbool (default = true)
save images of the ricci flow solution

Returns:

Saves to disk
————-
Ricci flow solution – ricciFn = sprintf(QS.fullFileBase.rawRicciSolution, maxIter, tp)
rawRicciMesh has fields annulus and rectangle – rawRicciMeshFn = sprintf(QS.fullFileBase.rawRicciMesh, maxIter, tp) ;
Beltrami coeff mu for this #iterations – mufn = sprintf(QS.fullFileBase.rawRicciMu, maxIter, tp) ;
Images saved
————
Plot Ricci result in annulus – fullfile(imDir, sprintf(‘%06d_RicciSolution.png’, tp)) ;
plot beltrami for Ricci flow – fullfile(imDir, sprintf(‘%06d_RicciFlowBeltrami.png’, tp)) ;
Plot just the bare triangulation – fullfile(imDir, sprintf(‘%06d_RicciFlowSolution.png’, tp)) ;
Histogram |mu| for each case – fullfile(imDir, sprintf(‘%06d_BeltramiCoefficients.png’, tp)) ;
Image of corrected vertices on inner and outer annulus – fullfile(imDir, sprintf(‘%06d_ricci_InnerCorrection.png’, tp)) ; fullfile(imDir, sprintf(‘%06d_ricci_OuterCorrection.png’, tp)) ;
Orientation and branch cut – fullfile(imDir, sprintf(‘%06d_DrhoDphi_flipped.png’, tp)) ; fullfile(imDir, sprintf(‘%06d_phiOrderInitial.png’, tp)) ; fullfile(imDir, sprintf(‘%06d_phiOrderFinal.png’, tp)) ;

NPMitchell 2020

@TubULAR.measureUVPrimePathlineFeatureIDs(QS, options)#

Interactively identify one longitudinal (zeta) position per feature for Lagrangian pathlines using field measurements along Lagrangian pathlines of the coordinate system (u’,v’) [conformally mapped to plane]

One can use, for example, radius of the pathlines on the mesh from the centerline, or the normal velocity, or dz – the ratio of projected (embedding space) proper distance of a unit vector along longitudinal mesh direction to the pullback distance along the longitudinal direction, or dp – similar for circumferential direction, divv – the divergence of the velocity field, etc. Uses two fields to identify the feature positions along the longitudinal dimension.

todo: generalize beyond vP pathline array

See also

measurePathlineFeatureIDs.m

Parameters:

QS (QuapSlap class object instance) –
pathlineType (str specifier ('vertices', 'piv', 'faces')) – whether pathlines in question thread through mesh vertices, piv evaluation coordinates, or face barycenters at t=t0Pathlines
options (optional struct with fields) – overwrite : bool field1 : str specifier (‘radius’, ‘dz’, ‘dp’, ‘divv’, ‘veln’) field2 : str specifier (‘radius’, ‘dz’, ‘dp’, ‘divv’, ‘veln’)

Returns:

featureIDs (#features x 1 int array) – longitudinal pullback coordinate for features, as an index into the pullback pathlines (which are grid-like).
NPMitchell 2020

@TubULAR.measurePathlineMetric(QS, options)#

measurePathlineStrain(QS, options): Integrate the metric strain rate along Lagrangian pathlines. Allow for median filtering along Lagrangian pathlines to avoid spurious spikes in accumulated strain. Plot results in 2d and/or 3d for each timepoint.

Parameters:

QS (QuapSlap class instance) –
options (struct with fields) – overwrite : bool, overwrite data on disk overwriteImages : bool, overwrite images of results on disk plot_comparison : bool, plot a comparison with DEC traceful dilation median_filter_strainRates : bool
2020 (NPMitchell) –

@TubULAR.getMeshes(tubi, overwrite, method)#

Obtain mesh surfaces of volumetric data (like ImSAnE’s surface detection methods), here using level sets activecontours, akin to ImSAnE’s integralDetector method. Note that you can adjust the starting configuration of the level set by using the output of a previous timepoint (which can be modified before evolution using pre_pressure and pre_tension, for example by seeding with a sphere at centerguess.

Parameters:

tubi (current TubULAR instance) –
overwrite (bool) – overwrite previous meshes on disk
method (str) – ‘threshold’ or ‘activecontour’
finding (Note that tubi.xp.detectOptions is unpacked to set parameters for) –
meshes –
preview : bool tension : prepressure : maxIterRelaxMeshSpikes : int >=0 (default=1000)

maximum number of iterations to relax spikes in the mesh

smooth_with_matlabint or float
if positive, uses matlab to smooth mesh with this parameter as the diffusion coefficient. Otherwise, if zero, performs no smoothing. Otherwise, if negative, uses MeshLab to smooth the mesh, using the specified mlxprogram in tubi.xp.detectOptions.mlxprogram
including –
preview : bool tension : prepressure : maxIterRelaxMeshSpikes : int >=0 (default=1000)

maximum number of iterations to relax spikes in the mesh

smooth_with_matlabint or float
if positive, uses matlab to smooth mesh with this parameter as the diffusion coefficient. Otherwise, if zero, performs no smoothing. Otherwise, if negative, uses MeshLab to smooth the mesh, using the specified mlxprogram in tubi.xp.detectOptions.mlxprogram

Returns:

none
Saves to disk
————-
raw meshes as ply files
level sets as .mat files,
smoothed meshes as ply files

NPMitchell 2022

@TubULAR.measureUVPrimePathlines(QS, options)#

measurePullbackPathlines(QS, options): Measure pathlines of optical flow in pullback space. These pathlines can then be used to query velocities or other properties that are parameterized by pullback coordinates. For example, we average velocites along Lagrangian pathlines. For another example, to build a Lagrangian-frame measure of divergence of tangential velocity, div(v_t), we can query div(v_t) at the coords of the PullbackPathlines (via interpolation of div(v_t) defined on pullback mesh vertices). For now assumes all pullbacks are the same size – for ex, 2000 x 2000 pixels.

Parameters:

QS (QuapSlap class instance) –
options (struct with fields) –

overwritebool
overwrite previous results

previewbool
view intermediate results

smoothing_sigmafloat
Gaussian kernel standard deviation if >0 for smoothing in-plane displacements from piv readout
disk (Saves to) –
------------- –
sprintf(QS.fileName.pathlines_uvprime.XY –
; (t0)) –
sprintf(QS.fileName.pathlines_uvprime.vXY –
; –
sprintf(QS.fileName.pathlines_uvprime.fXY –
; –
sprintf(QS.fileName.pathlines_uvprime.XYZ –
; –

See also

timeAverageVelocities, options, pullbackPathlines, options, NPMitchell

@TubULAR.visualizeSegmentationPatch(QS, Options)#

visualizeSegmentationPatch(QS, Options) Load demo tracks (segmented or tracked objects) and plot over pullback data from small patches of the surface mapped to the plane and registered.

Todo: allow for pathline points near periodic boundary.

Parameters:

QS (QuapSlap class instance) –
options (optional struct with fields) –
coordSys : str specifier timePoints : int array

timepoints to visualize, must have segmentation for these timepoints

t0int
reference timepoint to convert to t=0 for timestamp

t0Pathlinesint
timepoint at which pathlines form grid in pullback coords

overwrite : bool buffer : float or int or numeric bufferX : float or int or numeric, trumps buffer for X dim bufferY : float or int or numeric, trumps buffer for Y dim preview : preview (default=false) scaleByMetric : bool (default=false)
rescale parameterization to have det(g) approx 1 on either:
faces queried by pts (if both imW and imH are not supplied)

in window com(pts) +/- [imW*0.5, imH*0.5]
scaleByMetricComponentsbool (default=true)
rescale parameterization to approximate isothermal coordinates

imWint or numeric (default=25)
halfwidth of output image in QS.spaceUnits (approx) (in units of quasi-embedding space –> ie registered/flattened embedding submesh units, should be close to units of QS.spaceUnits if submesh is not too strongly curved in embedding space).

imHint or numeric (default=25)
halfheight of output image in QS.spaceUnits (approx) (in units of quasi-embedding space –> ie registered/flattened embedding submesh units, should be close to units of QS.spaceUnits if submesh is not too strongly curved in embedding space).

Returns:

<none>
Saves to disk
————-
image sequence of RGB overlays

See also

computeLocalSurfacePatch, pts, options, visualizeDemoTracks, options, NPMitchell

@TubULAR.measureCurvatures(QS, options)#

measureCurvatures(QS, options): Measure and plot mean and gaussian curvature for smoothed meshes (QS.currentMesh.spcutMeshSm) for each timepoint

Parameters:

QS (QuapSlap class instance) – Current QuapSlap object
options (struct with fields) –

overwritebool, default=false
overwrite results on disk
2020 (NPMitchell) –

@TubULAR.measureMetricStrainRate(QS, options)#

Advect the current mesh by the current raw velocities and measure differences between metric of advected (future) mesh and current mesh (reference). This gives a measurement of the strain rate tensor with units related to the pullback space.

Parameters:

tubi –
options (struct with fields) –

Returns:

<none>
Saves to disk
————-

NPMitchell 2020

class @TubULAR.Tubular_local_keep_these_changes_incorporate_202207#

Bases: handle

Tube-like sUrface Lagrangian Analysis Resource (TubULAR) class. Note that TubULAR is fully documented, with docs accessible via:

doc TubULAR

Overview Usage#

Instantiation Set experiment metadata and any custom options……………………..xp=struct(…); opts=struct(…); Instantiate TubULAR……………………………………………..tubi = TubULAR(xp, opts)

Surface extraction Create downsampled volumes to use for iLastik & level sets…………..tubi.prepareIlastik() Extract the surfaces using level sets on output of iLastik or raw data..tubi.getMeshes()

Parameterization Set experiment metadata and any custom options……………………..xp=struct(…); opts=struct(…); Instantiate TubULAR……………………………………………..tubi = TubULAR(xp, opts) Define global reference frame (APDV)………………………………tubi.computeAPDVCoords() Define endcap points and one point for a virtual seam……………….tubi.computeAPDpoints() Align meshes into global frame……………………………………tubi.alignMeshesAPDV() Optional: render data on dynamic surfaces………………………….tubi.plotSeriesOnSurfaceTexturePatch() Compute initial centerlines………………………………………tubi.extractCenterlineSeries() Temporally average centerlines and fix any inconsistencies…………..tubi.generateCleanCntrlines() Remove the endcaps………………………………………………tubi.sliceMeshEndcaps() Clean up the mesh endcaps………………………………………..tubi.cleanCylMeshes() Cut the meshes along virtual seams………………………………..tubi.generateCurrentCutMesh() Constrained parameterization of the mesh into (u,v) and (s,phi) frames..tubi.generateCurrentSPCutMesh() Generate pullback images…………………………………………tubi.generateCurrentPullbacks()

Dynamics Initial pass of smoothing embedding over time………………………tubi.smoothDynamicSPhiMeshes() Create pullbacks of smoothed (s,phi) frames………………………..tubi.generateCurrentPullbacks() Tile the images in phi (as ‘double-cover’)…………………………tubi.doubleCoverPullbackImages() Measure any residual material motion in the parameterization…………tubi.measurePIV2d() Translate pullback information into motion in 3D (embedding space)……tubi.measurePIV3d()

Advanced processing of surface dynamics Recommended: smooth the resulting 3D dynamics along material pathlines..tubi.timeAverageVelocities() Compute material pathlines in 3D………………………………….tubi.measurePullbackPathlines() Measure rate of strain over the surface at each timepoint……………tubi.measureStrainRate() Query the rate of strain along pathlines…………………………..tubi.measurePathlineStrainRate() Measure the integrated strain along pathlines………………………tubi.measurePathlineStrain()

Interpretation and mode decomposition Visualize flows…………………………………………………tubi.plotTimeAvgVelocities() Decompose into Divergence and curl………………………………..tubi.helmotzHodge() Compute areal rate of change……………………………………..tubi.measureMetricKinematics() Decompose into Principal Component Analysis………………………..tubi.getPCAoverTime() Decompose into eigenmodes of the Laplace-Beltrami operator…………..tubi.getLBSoverTime()

Coordinate Systems#

uv(conformal map onto unit square)

Conformally mapping the cylinderCutMesh onto the unit square in the plane results in the instantaneous uv coordinate system. Corners of the unit square are taken directly from the cutMesh, so are liable to include some overall twist.

sphi(proper length x rectified azimuthal coordinate)

uvprime(conformal map

[same as uvprime_sm, since uvprime is currently computed via sphi_sm coordinates]

r/spr/sphir(proper length x rectified azimuthal coordinate)

same as sphi but with aspect ratio relaxed to minimize isoareal energy cost –> makes triangles more similar in area by scaling the s axis (longitudinal axis) by a scalar factor.

iLastik’s coordinate systemsI recommend using cxyz as the axis

when outputting from iLastik. For our setup, this will mean that the meshes are mirrored with respect to the lab frame, so tubi.flipy = true.

PIV measurements fall into two classes:

‘piv’: principal surface-Lagrangian-frame PIV (sp_sme or up_sme)
–> note that the designation of coordinate system is not explicitly specified in the filenames: tubi.dir.mesh/gridCoords_nU0100_nV0100/piv/piv3d, etc
‘piv_uvp_sme’: PIV in coordSys for quasiconformal measurements
–> note that these are less Lagrangian than sp_sme or up_sme, so they are treated as independent from the principal pipeline in which we measure velocities in a surface-Lagrangian frame (ie sp_sme or up_sme).

Properties#

vsmM(#timePoints-1) x (nX*nY) x 3 float array
3d velocities at PIV evaluation coordinates in um/dt rs

vfsmM(#timePoints-1) x (2*nU*(nV-1)) x 3 float array
3d velocities at face barycenters in um/dt rs

vnsmM(#timePoints-1) x (nX*nY) float array
normal velocity at PIV evaluation coordinates in um/dt rs

vvsmM(#timePoints-1) x (nU*nV) x 3 float array
3d velocities at (1x resolution) mesh vertices in um/min rs

v2dsmM(#timePoints-1) x (nX*nY) x 2 float array
2d velocities at PIV evaluation coordinates in pixels/ min

v2dsmMum(#timePoints-1) x (nX*nY) x 2 float array
2d velocities at PIV evaluation coordinates in scaled pix/min, but proportional to um/min (scaled by dilation of map)

APDV coordinate system and Centerline specification#

Property Summary

APDV#: translation vector to transform data into APDV frame (rot*v+trans)*resolution

a_fixed#: aspect ratio for fixed geometry pullback meshes

apdvOptions#: Options for computing the APDV frame

cleanCntrlines#: centerlines in embedding space, no temporal filtering

currentData#: image intensity data in 3d and scaling

currentMesh#: ricci flow result pullback mesh, topological annulus

currentStrain#: beltrami coefficient for pathlines

currentTime#: int, the timepoint of the current frame under examination (this can be set with setTime and changes during many method calls)

data#: options for scaling and transposing image intensity data

dir#: str, directory where QuapSlap data lives

dynamic#: expMeta : struct with fields fileMeta : struct with fields

endcapOptions#: the full map from embedding to pullback being [M’=(Phi)o()o()]. This string specifier must be ‘3dcurves’ (geometric phi stabilization) or ‘texture’ (optical flow phi stabilization)

fileBase#: fileNames to be populated by timestamp (’…%06d…mat’)

fileName#: fileName struct with fields such as rot, trans, etc where measurements/parameters are stored on disk

flipy#: whether data is mirror image of lab frame coordinates

fullFileBase#: full path of filenames (like fullfile(tubi.dir.X, tubi.fileBase.X))

imSize#: size of pullback images to create (default is [a_ratio * 1000, 1000])

nU#: sampling number along longitudinal axis

nV#: sampling number along circumferential axis

normalShift#: shift to apply to meshes in pixel space along normal direction

phiMethod#: method for determining Phi map in pullback mesh creation, with

piv#: sigma of gaussian smoothing on PIV, in units of PIV sampling grid pixels

plotting#

at which we cut the cylinder mesh into a cylinderCutMesh (a topological disk/square). This “dorsal” point for other timepoints are identified by pointmatching.

Additional fields allowed :

smoothing#: half-width of tripulse filter applied along zeta/z/s/u direction in pullback space, in units of du/dz/ds/dzeta

spaceUnits#: units of the embedding space (ex ‘$mu$m’)

ssfactor#: subsampling factor for probabilities

t0#: reference time in the experiment

timeInterval#: increment in time between timepoints with

timeUnits#: indices differing by 1. For example, if timePoints are [0,1,2,4] and these are [1,1,2] minutes apart, then timeInterval is 1.

uvexten#: naming extension with nU and nV like ‘_nU0100_nV0100’

velocityAverage#: velocity field on vertices after Lagrangian avg

velocityRaw#: velocity field on vertices, no temporal filtering

xp#: struct with fields or ImSAnE experiment class instance

Method Summary

APDV2dxyz(tubi, a)#: ars = xyz2APDV(tubi, a) Transform 3d vectors from APDV coord sys to XYZ data space

APDV2xyz(tubi, a)#: ars = xyz2APDV(tubi, a) Transform 3d coords from APDV coord sys to XYZ data space

static XY2uv(im, XY, doubleCovered, umax, vmax)#

Map pixel positions (1, sizeImX) and (1, sizeImY) to (0, umax) and (0, vmax) of pullback space if singleCover, or y coords are mapped to (-0.5, 1.5)*vmax if doubleCover

NOTE THAT MAP IS [xesz, yesz] = [size(im, 1), size(im, 2)] uv(:, 1) = umax * (XY(:, 1) - 1) / xesz ; uv(:, 2) = vmax * 2.0 * (XY(:, 2) - 1) / yesz - 0.5 ;

NOTE THAT INVERSE MAP IS x–> (xy(:, 1) * (Xsz-1)) / (1*umax) + 1 , … y–> (xy(:, 2) * (Ysz-1)) / (2*vmax) + 1 + (Ysz-1)*0.25 ;

Parameters:

im (NxM numeric array or 2 ints for size) – image in which pixel coordinates are defined or dimensions of the image (pullback image in pixels)
XY (Qx2 numeric array) – pixel coordinates to convert to pullback space
doubleCovered (bool) – the image is a double cover of the pullback (extended/tiled so that the “top” half repeats below the bottom and the “bottom” half repeats above the top. That is, consider im to be a double cover in Y (periodic in Y and covers pullback space twice (-0.5 * Ly, 1.5 * Ly)
umax (float) – extent of pullback mesh coordinates in u direction (X)
vmax (float) – extent of pullback mesh coordinates in v direction (Y) before double covering/tiling

Returns:

uv – pullback coordinates of input pixel coordinates

Return type:

Qx2 numeric array

adjustIV(tubi, IV, adjustlow, adjusthigh, forceValues)#

Parameters:

tubi (QuapSlap class instance) –
IV (3D numeric (optional, loads currentData if empty)) – data to adjust
adjustlow (numeric or #channels x 1 numeric) – minimum intensity or percent intensity to adjust data
adjusthigh (numeric or #channels x 1 numeric) – maximum intensity or percent intensity to adjust data
forceValues (enforce the adjustlow/high to be intensity) –
values –
percentile (not) –
100. (even if <) –

clearTime(tubi)#: clear current timepoint’s data for tubi instance

static clipXY(xx, yy, Lx, Ly)#

Clip x at (1, Lx) and clip Y as periodic (1=Ly, Ly=1), for image that is periodic in Y. Consider Y in [1, Ly]. If the pullback is a doubleCover, Ly = 2*mesh width in pixels

Parameters:

xx –
yy –
Lx (int) – number of pixels along x dimension
Ly (int) – number of pixels along y dimension

Returns:

[xx, yy]

Return type:

x and y coordinates clipped to [1, Lx] and [1, Ly]

coordSystemDemo(QS, options)#: Image for publication/presentation on method & coordinate system Create coordinate system charts visualization using smoothed meshes options : struct with fields

style : ‘curves’ or ‘surface’

static dXY2duv(im, dXY, doubleCovered, umax, vmax)#

XY2uv(im, XY, doubleCovered, umax, vmax): Map difference in pixel positions defined on (1, sizeImX) and (1, sizeImY) to (0, umax) and (0, vmax) of pullback space if singleCover, or y coords are mapped to (-0.5, 1.5)*vmax if doubleCover

NOTE THAT FULL COORDINATE MAP IS [xesz, yesz] = [size(im, 1), size(im, 2)] uv(:, 1) = umax * (XY(:, 1) - 1) / xesz ; uv(:, 2) = vmax * 2.0 * (XY(:, 2) - 1) / yesz - 0.5 ;

NOTE THAT INVERSE MAP IS x–> (xy(:, 1) * (Xsz-1)) / (1*umax) + 1 , … y–> (xy(:, 2) * (Ysz-1)) / (2*vmax) + 1 + (Ysz-1)*0.25 ;

Parameters:

im (NxM numeric array or 2 ints for size) – image in which pixel coordinates are defined or dimensions of the image (pullback image in pixels)
dXY (Qx2 numeric array) – difference in pixel coordinates as vector, to convert to pullback space
doubleCovered (bool) – the image is a double cover of the pullback (extended/tiled so that the “top” half repeats below the bottom and the “bottom” half repeats above the top. That is, consider im to be a double cover in Y (periodic in Y and covers pullback space twice (-0.5 * Ly, 1.5 * Ly)
umax (float) – extent of pullback mesh coordinates in u direction (X)
vmax (float) – extent of pullback mesh coordinates in v direction (Y) before double covering/tiling

Returns:

uv – pullback coordinates of input pixel coordinates

Return type:

Qx2 numeric array

doubleCoverPullbackImages(tubi, options)#

options : struct with fields coordsys : (‘sp’, ‘uv’, ‘up’)

coordinate system to make double cover

overwritebool, default=false: whether to overwrite current images on disk
histeqbool, default=true: perform histogram equilization during pullback extension
ntilesint, default=50: The number of bins in each dimension for histogram equilization for a square original image. That is, the extended image will have (a_fixed * ntiles, 2 * ntiles) bins in (x,y).
a_fixedfloat, default=tubi.a_fixed: The aspect ratio of the pullback image: Lx / Ly

static doubleToSingleCover(XY, Ly)#: detect if XY is passed as a pair of grids

static dvAverageNematic(magnitude, theta)#

Parameters:

magnitude (nU x nV numeric array) – magnitude of nematic field in 2D rectilinear grid
theta (nU x nV float array, with values as angles in radians) – angle (mod pi) of nematic field in 2D rectilinear grid

Returns:

mag_ap (nU x 1 float array) – average magnitude of nematic field along the ap dimension
theta_ap (nU x 1 float array) – average angle (mod pi) of nematic field along the ap dimension

dx2APDV(tubi, da)#: dars = dx2APDV(tubi, da) Transform 3d difference vector from XYZ data space to APDV coord sys

getAPpointsSm(tubi)#: Load the anterior and posterior ‘centers of mass’ ie the endpoints of the object’s centerline

getCurrentPathlineStrain(tubi, t0Pathlines, varargin)#

Parameters:

tubi (current class instance) –
t0Pathlines (int or empty) – reference timepoint. If empty, set to tubi.t0set()
varargin (strings ('strain', 'beltrami')) – which strain measures to load from disk and attribute to self

Returns:

strain

Return type:

tubi.currentStrain.pathline

getCurrentSegmentation2D(tubi, options)#: Obtain the cell segmentation in 3D pullback space

getCurrentSegmentation2DCorrected(tubi, options)#: Decide on coordinate system for corrected pullback binary

getCurrentSegmentation3D(tubi, options)#: Obtain the cell segmentation in 3D pushforward space

getCurrentSegmentation3DCorrected(tubi, options)#: Obtain the cell segmentation in 3D pullback space

getPIV(tubi, options)#: Load PIV results and store in tubi.piv if not already loaded

getPullbackPathlines(tubi, t0, varargin)#: Discern if we must load pathlines or if already loaded

getRadii(tubi, options)#

Load radii from disk in specified coordinate system

Parameters:

options (struct with fields) –

coordSys: str specifier for how to return radii: ’spcutMeshSmRSC’ –> return nU x (nV-1) array for radii of vertices in spcutMeshSmRSC embedding <add other options here>

getRotTrans(tubi)#: Load the translation to put anterior to origin and AP axis along x axis

getVelocityAverage(tubi, varargin)#: todo: check if all varargin are already loaded

getVelocityRaw(tubi, varargin)#: todo: check if all varargin are already loaded

getVelocitySimpleAverage(tubi, varargin)#: todo: check if all varargin are already loaded

getXYZLims(tubi)#: [raw, pix, um, um_buff] = GETXYZLIMS(tubi) Grab each xyzlim from self, otherwise load from disk full resolution pix

static invertRotation(rot)#: Obtain rotation matrix that undoes the APDV rotation

loadPIV(tubi, options)#: Load PIV results from disk and store in tubi.piv

loadVelocityAverage(tubi, varargin)#: Retreive the Velocity on each vertex/face barycenter of the (s,phi) meshes. NOTE: that these are nearly Lagrangian for most applications, but not perfectly Lagrangian. See also: tubi.getPullbackPathlines(t0), which corrects for any residual motion in the pullback plane.

loadVelocityRaw(tubi, varargin)#: Load and pack into struct

loadVelocitySimpleAverage(tubi, varargin)#: Load and pack into struct

measureRMSvelocityOverTime(tubi, options)#

Parameters:

options (optional struct with fields) – weights : weights to attribute to each face in mesh coordSys : str specifier (default==’sp’)

Returns:

vrms (#tps x 1 float array) – rms velocity over time
timestamps (#tps x 1 float array) – timestamps associated with each entry

plotPathlineBeltramiKymograph(tubi, t0Pathlines, options)#: Example usage for 2021 gut paper: options = struct(‘ylim’, [0, 2] ) tubi.plotPathlineBeltramiKymograph([], options)

static quarterIndicesDV(nV)#: [dorsal, ventral, left, right] = quarterIndicesDV(nV) indices for each quarter of a DV section in grid coordinates

samplePullbackPathlines(tubi, XY0, options)#: [p2d, p3d] = samplePullbackPathlines(tubi, XY0, options) start at XY0, folow flow using barycentric coordinates of PIV pullback pathlines. This is the same as advecting along a fixed Langrangian coordinate.

setDataLimits(tubi, tp, adjustlow_pctile, adjusthigh_pctile)#: Use timepoint (tp) to obtain hard values for intensity limits so that data is rescaled to fixed limits instead of percentile. This is useful to avoid flickering of overall intensity in data in which a few voxels vary a lot in intensity.

setTime(tubi, tt)#

Set the current time of the dataset and clear current data which was associated with the previously considered time

Parameters:: tt (int or float) – timePoint to set to be current, from available times in tubi.xp.fileMeta.timePoints

t0set(tubi, t0)#: t0set(tubi, t0) Set time offset from file or manually

trueTime(tubi, tp, div60)#

Convert timepoint into timestamp in physical units (like minutes, hours, etc) relative to t0

Parameters:

tubi (TubULAR class instance) –
tp (int) – timestamp to convert to “true” time relative to t0 in real units
div60 (bool) – Convert to a larger time unit than tubi.timeUnits, so if tubi.timeUnits is sec(onds), convert to minutes. Or, if tubi.timeUnits is min(utes), convert to hours.

Returns:

tpTrue (float) – the timestamp relative to t0 in true units (tubi.timeUnits)
timestr (char) – string with timestamp relative to t0 and its units

uv2APDV(tubi, uv, coordSys, umax, vmax)#: Convert (u,v) pullback coordinates to xyz coordinates in the APDV frame (rotated and scaled frame aligned with AP axis along x and DV axis along z)

static uv2XY(im, uv, doubleCovered, umax, vmax)#

XY = uv2XY(im, uv, doubleCovered, umax, vmax): Map from pullback uv u=(0,1), v=(0,1) to pixel XY

x–> (xy(:, 1) * (size(im, 2)-1)) / (1*umax) + 1 , … y–> (xy(:, 2) * (size(im, 1)-1)) / (2*vmax) + 0.75 + (size(im,1)-1)*0.25

Parameters:

im (NxM numeric array or length(2) int array) – 2D image into whose pixel space to map or size(im)
uv (Q*2 numeric array) – mesh coordinates to convert to pullback pixel space (XY)
doubleCovered (bool) – the image is a double cover of the pullback (extended/tiled so that the “top” half repeats below the bottom and the “bottom” half repeats above the top. That is, consider im to be a double cover in Y (periodic in Y and covers pullback space twice (-0.5 * Ly, 1.5 * Ly)
umax (float) – extent of pullback mesh coordinates in u direction (X)
vmax (float) – extent of pullback mesh coordinates in v direction (Y) before double covering/tiling

Returns:

XY – positions of uv coordinates in pullback pixel space

Return type:

N x 2 float array

xyz2APDV(tubi, a)#: ars = xyz2APDV(tubi, a) Transform 3d coords from XYZ data space to APDV coord sys

@TubULAR.measurePIV3dMultiChannel(QS, options)#

Measure 3d flows from PIV results on smoothed (s,phi) MIP pullbacks for each channel

default options for created MIPS (done before this method) are:
pbOptions = struct() ; pbOptions.overwrite = false ; pbOptions.numLayers = [5, 5] ; pbOptions.layerSpacing = 0.75 ; pbOptions.generate_rsm = true ; pbOptions.generate_spsm = true ; pbOptions.generate_sphi = false ; QS.data.adjustlow = 1.00 ; QS.data.adjusthigh = 99.999 ; QS.generateCurrentPullbacks([], [], [], pbOptions) ;

Parameters:

QS (QuapSlap class instance) –
options (struct with fields) –

overwritebool
overwrite previous results

previewbool
view intermediate results

timePointsnumeric 1D array
the timepoints to consider for the measurement. For ex, could choose subset of the QS experiment timePoints

Saves to disk#

m0XY : Px2 float: 2d mesh vertices in pullback image pixel space [pullback pix[” ; m0f : #facesx3 float: mesh connectivity list” ; m0v3d : Px3 float: 3d mesh coordinates in embedding pixel space [mesh pix / dt]” ; x0 : QxR float: x value of 2d velocity evaluation coordinates in pullback image pixel space [pullback pix]” ; y0 : QxR float: y value of 2d velocity evaluation coordinates in pullback image pixel space [pullback pix]” ; pt0 : Nx3 float: 3d location of evaluation points [mesh pix]” ; pt1 : Nx3 float: 3d location of advected point in next mesh [mesh pix]” ; v0 : Nx3 float: 3d velocity [mesh pix / dt]” ; v0n : Nx1 float: normal velocity [mesh pix / dt]” ; v0t : Nx3 float: tangential velocity in 3d [mesh pix / dt]” ; facenormals : #faces x 3 float: face normals for all faces [unitless, mesh pix / mesh pix]” ; v3dfaces : #faces x 3 float: face-centered velocities in embedding space [mesh pix / dt]” ; v3dvertices : #mesh vertices x 3 float: vertex-centered velocities in embedding space [mesh pix / dt]” ; v0_rs : Nx3 float: rotated/scaled 3d velocities [um/min]” ; v0n_rs : Nx3 float: scaled normal velocities [um/min]” ; v0t_rs : Nx3 float: rotated/scaled tangential 3d velocities [um/min]” ; normals_rs : Nx3 float: normal vector of field faces in 3d [unitless, um/um]” ; v0t2d : Nx2 float: in-plane velocity [pullback pix / dt]” ; g_ab : Nx2x2 float: pullback metric” ; dilation : Nx1 float: dilation of face from 3d to 2d (A_2d [pullback_pix^2] / A_3d [mesh_pix^2])” ; jac : #faces x 1 cell: jacobian of 3d->2d transformation” ; fieldfaces : Nx1 int: face indices into spcutMesh where v defined” ;

NPMitchell 2020-2021

@TubULAR.computeAPDVCoords(tubi, opts)#

[acom,pcom,dcom, rot,trans] = computeAPDVCoords(tubi, opts) Compute the APDV coordinate system, defined WRT the data coordSys by a rotation, translation, and resolution. This is done by identifying a dorsal point, so that z dim in APDV points from the AP axis to dorsal along the shortest linesegment emanating from the AP axis to the dorsal point. By default, use_iLastik==false, so we just look at the reference time’s mesh in 3D, pick out the data volume axis closest to the long axis of the mesh, and point the dorsal vector perpendicular to that. For ex, if we have an embryo surface elongated along the first dimension of our data volume, then the first dimension becomes our AP axis, and the third dimension is the DV axis.

If acom,pcom,dcom are all inferred from iLastik training, the default channels used to extract the blobs of high probability marking A,P, and a point D that is dorsal to the line connecting AP are:

A : 1 P : 2 D : 4

This would leave channel 3 (and 5+) to be a background channel in the iLastik training.

Chirality: note that the probabilties field is assumed to be meshlike so the first two axis are swapped within com_region().

Todo: back-save acom_for_rot and pcom_for_rot for datasets already processed, inferred from rot and trans

Note that axisOrder is applying upon invoking getCurrentData()

@TubULAR.generatePathlineRicciMeshTimePoint(QS, tp, options)#

generatePathlineRicciRicciMeshTimePoint(QS, tp, options)

Compute Ricci flow pullback for pathline-advected vertices in QS.dir.piv/pathlines/t0_####/quasiconformal/pb2pb/ which is the same as QS.dir.pathlines.quasiconformal/pb2pb

Parameters:

QS (quapSlap class instance) –
tp (int (default=QS.t0)) – timestamp of timepoint for which to generate Ricci flow mesh
options (optional struct with optional fields) –
maxIter : int (default=100) radiusTolerance : float (default=0.01)

maximum allowed fractional deviation of Ricci flow solution inner and outer radius from a true circle with fixed radius (variable inner radius, fixed outer radius of 1)

save_imsbool (default = true)
save images of the ricci flow solution

Returns:

Saves to disk
————-
Ricci flow solution – ricciFn = sprintf(QS.fullFileBase.ricciSolution, maxIter, tp)
ricciMesh has fields annulus and rectangle – ricciMeshFn = sprintf(QS.fullFileBase.ricciMesh, maxIter, tp) ;
Beltrami coeff mu for this #iterations – mufn = sprintf(QS.fullFileBase.ricciMu, maxIter, tp) ;
Images saved
————
Plot Ricci result in annulus – fullfile(imDir, sprintf(‘%06d_RicciSolution.png’, tp)) ;
plot beltrami for Ricci flow – fullfile(imDir, sprintf(‘%06d_RicciFlowBeltrami.png’, tp)) ;
Plot just the bare triangulation – fullfile(imDir, sprintf(‘%06d_RicciFlowSolution.png’, tp)) ;
Histogram |mu| for each case – fullfile(imDir, sprintf(‘%06d_BeltramiCoefficients.png’, tp)) ;
Image of corrected vertices on inner and outer annulus – fullfile(imDir, sprintf(‘%06d_ricci_InnerCorrection.png’, tp)) ; fullfile(imDir, sprintf(‘%06d_ricci_OuterCorrection.png’, tp)) ;
Orientation and branch cut – fullfile(imDir, sprintf(‘%06d_DrhoDphi_flipped.png’, tp)) ; fullfile(imDir, sprintf(‘%06d_phiOrderInitial.png’, tp)) ; fullfile(imDir, sprintf(‘%06d_phiOrderFinal.png’, tp)) ;

NPMitchell 2020

@TubULAR.cleanFastMarchingCenterlines(tubi, idOptions)#

cntrlines = cleanFastMarchingCenterlines(QS, idOptions) Identify anomalous centerlines in a time series, fix them and save.

Parameters:

tubi (TubULAR class instance) –
idOptions (struct with fields) –
- ssr_thresfloat
  threshold sum of squared residuals between adjacent timepoint centerlines to consider the later timepoint centerline to be anomalous & require cleaning
- overwritebool
  overwrite the centerlines on disk

Returns:

QS.cleanCntrlines (length(timePoints) x 1 cell) – populated with cntrlines in XYZ coords (data space, pixels)
QS.fileName.cleanCntrlines
NPMitchell 2020

@TubULAR.plotPathlineStrainTimePoint(tubi, tp, options)#

plotPathlineStrainRateTimePoint(tubi, tp, options): Plot the traceful and traceless components of the strain rate tensor defined at pathlines moving in the domain of parameterization over time. Here, plot the trajectories as immobile in a grid. Uses vertex pathlines.

Parameters:

tubi (TubULAR class instance) –
tp (int) – timepoint in units of (1/tubi.timeInterval) * tubi.timeUnits
options (struct with fields) –

t0Pathlineint
Lagrangian frame is taken to be cutMesh of this timepoint

LagrangianFrameStylestr specifier (‘sphi’ ‘uvprime’)
what kind of cutmesh pullback is strain measured relative to. In other words, if we computed strain with respect to a fixed material frame in (zeta, phi) coordinates defined at timepoint t0Pathline, then this would be ‘zphi’.
2020 (NPMitchell) –

@TubULAR.plotPathlineMetricKinematics(QS, options)#

Plot the metric Kinematics along pathlines as kymographs and correlation plots. Out-of-plane motion is v_n * 2H, where v_n is normal velocity and H is mean curvature. In-plane motion considered here is div(v_t) where v_t is tangential velocity on the curved surface. The difference div(v_t) - vn*2H = 1/2 Tr[g^{-1} dot{g}], which is a measure of isotropic metric change over time (dot is the partial derivative wrt time).

Parameters:

QS (QuapSlap class instance) –
options (struct with fields) – plot_kymographs : bool plot_kymographs_cumsum : bool plot_correlations : bool plot_fold_kinematics : bool plot_lobe_kinematics : bool

Returns:

<none>
NPMitchell 2020

@TubULAR.computeLBSoverTime(tubi, options)#

results = computeLBSoverTime(tubi, options) Compute Laplace-Beltrami spectral decomposition of the sequence of deformations over time.

Parameters:

tubi (TubULAR class instance) –
options (struct with fields) –

overwritebool (default = false)
overwrite previous results on disk

overwriteImagesbool (default = false)
overwrite figures of previous results on disk

convert_to_periodbool (default=false)
convert all timestamps to period

Tnumeric (default=NaN)
number of timepoints per period if data is periodic in time

t0int (default=tubi.t0set())
timestamp to mark as t=0

decompTint (defualt = t0)
The time for which the mesh Laplacian eigenbasis will be constructed

NLBSModesint (default = 200)
The number of Laplacien eigenvectors to calculate

NmodesToViewint (default= 5)
how many modes to plot individually on the reference surface configuration

nTimePoints2RmEndsint (default= 0)
how many timepoints to remove from the beginning and end of the timeseries. This is useful if the first and last few datapoints are noisy

drawArrowsInLBSPlanebool (default= false)
a plotting option, to draw arrows showing trajectory from timepoint to timepoint in 2D LBS projective planes

drawArrowsInLBS3Dbool (default=false)
a plotting option, to draw arrows showing trajectory from timepoint to timepoint in LBS space

meshChoicecell of strings or single string
(default={‘Lagrangian’, ‘sphi’}) the mesh style to use for computing PIV.

axPairsoptional int array (default=[1,2;1,3;2,3])
pairs of mode numbers to plot in 2d projective planes (modes are ordered by rank)

plotArrowsOnModesbool
plot arrows on the surfaces for mode decompositions

nArrowsint
subsampling factor for quiverplot on modes, only used if plotArrowsOnModes == true

Returns:

results

Return type:

cell of structs with fields

Saves to disk#

saves fullfile(tubi.dir.LBSoverTime, meshChoice, sprintf(‘lbsResults_%s.mat’, lbsType)) ;: where meshChoice is either ‘Lagrangian’ or ‘sphi’

saves figures to fullfile(tubi.dir.LBSoverTime, meshChoice)

Dillon Cislo + NPMitchell 2022

@TubULAR.computeAPDpoints(tubi, opts)#

[apts_sm, ppts_sm, dpt] = COMPUTEAPDPOINTS(opts) Compute the anterior, posterior, and dorsal points either from:

Clicking on the points on the mesh from t=t0

moments of inertia of the mesh surface, identifying the endpoints near the long axis of the object at t=t0. (autoAP = true)

iLastik training for CENTERLINE computation. (opts.use_iLastik=true) Note that these are allowed to be different than the APD points for ALIGNMENT computation. Here, we do not require any dorsal point to be trained, as that was defined in the previous step. However, it is useful to have A and P be separable from the previous step. For example, the posterior point might be a point which does NOT form an AP axis with the anteriormost point, as in the illustration of the midgut below. To use this option, set opts.use_iLastik=true and by default the script looks for iLastik files called <<fn>>_Probabilities_apcenterline.h5.

Directly supplying a custom set of anterior and posterior points for each time point

Posterior (distal) pt for centerline

_x_ Dorsal pt

/ / ___x_

/ / /

/ /____/ Anterior (proximal) pt for both centerline and for defining APDV axes

x P for APDV | x ________________/

(ventral here, unlabeled)

The default behavior is to use iLastik training if found, unless opts.use_iLastik is set to false. If use_iLastik is false or no iLastik output is found, the default is to have the user click on the points. If opts.autoAP == true, then we automatically find A and P positions by simply point matching from line intersections onto the mesh.

Parameters:

opts (struct with fields) –

use_iLastik : default=true if training h5s are present on disk
timePoints : timepoints for which to extract APD points
dorsal_thres : float between 0 and 1, threshold for COM extraction
anteriorChannelint, which channel of training to use if training
is used for extracting COM from probability cloud
posteriorChannelint, which channel of training to use if training
is used for extracting COM from probability cloud
dorsalChannelint, which channel of training to use if training
is used for extracting COM from probability cloud
overwrite : bool, overwrite previous results on disk
preview_com : bool, inspect the centers of mass extraction
axorder : length 3 int array, permutation of axes if needed
smwindowfloat or int (optional, default=30)
number of timepoints over which we smooth
preview : bool (optional, default=false)
autoAPbool whether or not to extract A-P points using an automatic
method using the moments of inertia of the mesh (optional, default=false)
custom_aptsa custom set of 3D anterior points
(optional, default=[])
custom_pptsa custom set of 3D posterior points
(optional, default=[])

OUTPUTS#

apdv_pts_for_centerline.h5 (rawapdvname, tubi.fileName.apdv): Raw points (centers of mass if training-based) for A, P, and D in subsampled pixels, in probability data space coordinate system. (note this is downsampled by ssfactor) Saved to fullfile(meshDir, ‘centerline/apdv_pts_for_centerline.h5’)
tubi.fileName.dpt: txt file with dorsal COM for APDV definition

rawapdvmatname=fullfile(tubi.dir.cntrline, ‘apdv_pts_for_centerline.mat’)

NPMitchell 2020

@TubULAR.generateCurrentUVCutMesh(tubi, cutMesh, spcutMeshOptions)#

Create a rectilinear (but periodic in 2nd dim) parameterization of the cutMesh with a nearly-conformal pullback projection uv.

uvcutMesh.f uvcutMesh.uv uvcutMesh.v uvcutMesh.vn

Parameters:

tubi (TubULAR class instance) –

Note that the following properties are used:
tubi.phiMethod = (‘3dcurves’, ‘texture’, ‘combined’) tubi.a_fixed = 2.0
cutMesh (cutMesh struct, optional) – cutMesh with fields
options (struct with fields, optional) –

overwritebool
overwrite previous results

Returns:

uvcutMesh (struct with fields) – uv:
NPMitchell 2020

@TubULAR.plotPathlineVelocitiesTimePoint(QS, tp, options)#

Plot the velocities defined on PIV evaluation points in Lagrangian coordinates.

Parameters:

QS (QuapSlap class instance) –
options (optional struct with fields) –

gridTopologystr
’rectilinear’ or ‘triangulated’ triangulate anew each timepoint or keep grid structure

vnsm : v2dsmum : overwrite : bool

overwrite previous results

previewbool
view intermediate results

timePointsnumeric 1D array
the timepoints to consider for the measurement. For ex, could choose subset of the QS experiment timePoints

alphaValfloat
the opacity of the heatmap to overlay

invertImagebool
invert the data pullback under the velocity map

qscalefloat
scale for quiver arrows

sampleIDxint array (default=[])
vertex indices at which to draw vectors. If not empty, uses ‘custom’ for subsamplingMethod. If empty, computes new farthestPointSearch sampling
2020 (NPMitchell) –

@TubULAR.measureBeltramiCoefficient(tubi, options)#

Previous version measured instantaneous mu as well, which was the Beltrami coefficient to the current (nearly) conformal pullback.

Parameters:

tubi (TubULAR class instance) –
options (struct with fields) –

t0Pathlinesnumeric (default = tubi.t0)
timepoint used to define Lagrangian/material frame for mu measurements to be made with respect to

coordSysstr specifier (‘ricci’, ‘uvprime’ or ‘sp’)
coordinate system in which pathlines are computed

overwrite : bool, overwrite previous results on disk climit : limit for caxis of scalar fields Re(mu) and Im(mu)

Return type:

NPMitchell 2020

@TubULAR.plotPathlineStrain(tubi, options)#

Plot the strain (from piv pathlines) along pathlines as kymographs and correlation plots. These are computed via finite differencing of pathline mesh metrics. That is, vertices are advected along piv flow and projected into 3D from pullback space to embedding, then the metrics g_{ij} of those meshes are measured and differenced.

Parameters:

tubi (QuapSlap class instance) –
options (struct with fields) – plot_kymographs : bool plot_kymographs_cumsum : bool plot_correlations : bool plot_fold_strainRate : bool plot_lobe_strainRate : bool plot_fold_strain : bool plot_lobe_strain : bool

Returns:

<none>
NPMitchell 2020

@TubULAR.pullbackPathlines(QS, x0, y0, t0, options)#

Create paths in pullback space (in pixels, XY) by following optical flow of PIV measured on QS.pivimCoords pullbacks. For non-standard PIV pathline propagation (such as uvprime coords), supply piv, Lx, and Ly in options. Note: I’ve chosen to use spatial smoothing via Gaussian blur rather than temporal smoothing of velocities (as is done in metric kinematics) for this code. This can be changed by setting QS.piv.smoothing_sigma == 0. Ie, if QS.piv.smoothing_sigma > 0, then we advect in the smoothed piv results.

Parameters:

QS (QuapSlap class instance) –
x0 (n*m float array) – x coordinates in pullback pixels to start pathlines at t0
y0 (n*m float array) – y coordinates in pullback pixels to start pathlines at t0
t0 (float) – time at which to begin the pathlines, must be member of QS.xp.fileMeta.timePoints. Note this is not the index of timePoints, but the element of timePoints (ie a timeStamp)
options (struct with fields) –

previewbool
view intermediate results

timePoints1d int array
the timpepoints along which to map pathlines that intersect with (x0, y0) at t0

Returns:

XX (#timePoints x size(x0, 1) x size(x0, 2) float array) – the pullback pixel X coordinates of the pathline spanning timePoints
YY (#timePoints x size(x0, 1) x size(x0, 2) float array) – the pullback pixel Y coordinates of the pathline spanning timePoints
NPMitchell 2020

@TubULAR.plotStrainRateTimePoint(tubi, tp, options)#

Plot the traceful and traceless components of the strain rate tensor defined on each face, using vertex-based results loaded from disk. Note that the vertex-based results are more heavily smoothed (via spectral filtering) than the face-based results on disk.

Parameters:

tubi (Additional optional fields of options that are already specified by) –
tp (int) – timepoint in units of (1/tubi.timeInterval) * tubi.timeUnits
options (struct with fields) –

overwritebool
Overwrite existing results on disk

meshstruct with fields f,v
current Mesh, allowed to be supplied to allow speedup (prevents delay from loading cutMesh from disk)

cutMeshstruct with fields f,v, pathPairs, nU, nV
current cutMesh, allowed to be supplied to allow speedup (prevents delay from loading cutMesh from disk)

clim_trace2x1 numeric
color limits for the colormap/bar of the isotropic strain rate

clim_deviatoric2x1 numeric
color limits for the colormap/bar of the deviatoric strain rate

samplingResolution’1x’ or ‘resampled’
currently only 1x resolution is allowed

averagingStyle’Lagrangian’ or ‘none’
currently we have restricted the velocity averaging style in tubular to either be Lagrangian averaging or none

debugbool
if true, show intermediate results

plot_comparisonbool
plot a comparison between trace of measured strain rate and measured DEC divergence minus H*normal velocity, where H is the mean curvature.
tubi –
are (metadata but can be overwritten by including them in options) – lambda = tubi.smoothing.lambda ; lambda_mesh = tubi.smoothing.lambda_mesh ; nmodes = tubi.smoothing.nmodes ; zwidth = tubi.smoothing.zwidth ;

NPMitchell 2020

@TubULAR.measurePathlineStrainRate(tubi, options)#

Query the metric strain rate along lagrangian pathlines. Plot results as kymographs.

Parameters:

tubi (TubULAR class instance) –
options (struct with fields) – plot_kymographs : bool plot_kymographs_cumsum : bool plot_gdot_correlations : bool plot_gdot_decomp : bool
2020 (NPMitchell) –

@TubULAR.plotStrainRate3DFiltered(tubi, options)#

load spatially-smoothed epsilon = 1/2 (nabla_i v_j + nabla_j v_i) - vN b_{ij} from disk, smooth with time filter and plot in 3d

Parameters:

tubi (QuapSlap class object instance) –
options (struct with fields) – lambda : float lambda_mesh : float overwrite : bool preview : bool averagingStyle : str specifier (‘Lagrangian’ or ‘Simple’) samplingResolution : ‘1x’ or ‘2x’ debug : bool nTimePoints

Returns:

Saves to disk
————-

NPMitchell 2020

@TubULAR.computeStartEndCOMs(QS, opts)#

[acom_sm, pcom_sm, dcom] = COMPUTEAPDCOMS(opts) Compute the anterior, posterior, and dorsal centers of mass from training

Parameters:

opts (struct with fields) –

timePoints
dorsal_thres : float between 0 and 1
anteriorChannel : int
posteriorChannel
dorsalChannel
overwrite : bool
apdvoutdir : string
meshDir : string
preview_com : bool
check_slices : bool
axorder : length 3 int array
smwindowfloat or int (optional, default=30)
number of timepoints over which we smooth
preview : bool (optional, default=false)

OUTPUTS#

apdv_coms_from_training.h5 (rawapdvname, QS.fileName.apdv): Raw centers of mass for A, P, and D in subsampled pixels, in probability data space coordinate system Saved to fullfile(meshDir, ‘centerline/apdv_coms_from_training.h5’)
QS.fileName.dcom: txt file with dorsal COM for APDV definition

rawapdvmatname=fullfile(QS.dir.cntrline, ‘apdv_coms_from_training.mat’)

NPMitchell 2020

@TubULAR.measurePathlineVelocities(tubi, options)#

measurePullbackStreamlines(tubi, options): Use pathlines of optical flow in pullback space to query velocities and average along pathlines.

Parameters:

tubi (TubULAR class instance) –
options (struct with fields) –

overwritebool, default=false
overwrite previous results

previewbool, default=false
view intermediate results

t0int, default=tubi.t0set()
timestamp at which the pathlines form a grid onto mesh vertices, mesh face barycenters, or PIV evaluation points

Returns:

outStruct = struct(‘v2dum’, v2dMum, … – ‘v2d’, v2dM, … ‘vn’, vnM, … ‘v3d’, vM, … ‘vf’, vfM, … ‘vv’, vvM, … ‘v2dsmum’, v2dsmMum, … ‘v2dsm’, v2dsmM, … ‘vnsm’, vnsmM, … ‘v3dsm’, vsmM, … ‘vfsm’, vfsmM, … ‘vvsm’, vvsmM ) ;
Saves to disk
————-
fvsmM = sprintf(tubi.fileName.pathlines.velocities.v3dsm, t0) ;
fvfsmM = sprintf(tubi.fileName.pathlines.velocities.vfsm, t0) ;
fvnsmM = sprintf(tubi.fileName.pathlines.velocities.vnsm, t0) ;
fvvsmM = sprintf(tubi.fileName.pathlines.velocities.vvsm, t0) ;
fv2dsmM = sprintf(tubi.fileName.pathlines.velocities.v2dsm, t0) ;
fv2dsmMum = sprintf(tubi.fileName.pathlines.velocities.v2dsmum, t0) ;
NPMitchell 2020

@TubULAR.helmholtzHodge(QS, options)#

helmholtzHodge(QS, options)

Compute DEC divergence and DEC “curl” on 2d evolving surface in 3d.
Compute laplacian in covariant frame
Compute scalar fields representing rotational and harmonic components
Store all on disk.

Parameters:

QS (QuapSlap class instance) –
options (struct with fields) –

overwritebool
overwrite previous results

previewbool
view intermediate results

averagingStylestr (‘Lagrangian’, ‘simple’)
style in which velocities were averaged over time

alphaValfloat
the opacity of the heatmap to overlay

invertImagebool
invert the data pullback under the velocity map

clipDivfloat (default=5.0)
max allowed divergence measurement

clipRotfloat (default=0.5)
max allowed vorticity measurement

computeLaplacianbool (default = false)
compute the laplacian of the velocity field
2020/2021 (NPMitchell) –

@TubULAR.measurePathlineStrain(tubi, options)#

Compute strain from integrated pathlines deforming mesh vertices. Measurements are taken with respect to fixed Lagrangian frame. Plot results in 2d and/or 3d for each timepoint.

Parameters:

tubi (QuapSlap class instance) –
options (struct with fields) – overwrite : bool, overwrite data on disk overwriteImages : bool, overwrite images of results on disk plot_comparison : bool, plot a comparison with DEC traceful dilation median_filter_strainRates : bool
2020 (NPMitchell) –

@TubULAR.smoothDynamicSPhiMeshes(QS, options)#

SMOOTHDYNAMICSPHIMESHES(QS, options): Using (s,phi) pullback cutmeshes, smooth coord system in time via simple triangular pulse averaging of positions in embedding space while preserving the pullback mesh. Note that the computed geodesic distance along the mesh is preserved in the pullback: mesh extent spcutMeshSm.u(:, 1) varies from (0, GeoLength)

Parameters:

QS (QuapSlap class instance) – Current QuapSlap object
options (struct with fields) –

widthint, default=4
half-width of tripulse smoothing filter on output meshes

overwritebool, default=false
overwrite results on disk

preivewbool, default=QS.plotting.preview
display intermediate results
2020 (NPMitchell) –

@TubULAR.interpolateOntoPullbackXY(tubi, XY, scalar_field, options)#

interpolateOntoPullbackXY(tubi, XY, scalar_field, options): Interpolate scalar field defined on current mesh onto supplied pullback coordinates. Input scalar_field may be defined on vertices or faces. Automatically detects whether field is at 1x mesh resolution or higher. The pixel locations of the pullback are computed from scratch via options.Lx and options.Ly, not passed directly.

Parameters:

tubi (TubULAR class instance) – class for which to interpolate onto the pullback XY
XY (Nx2 numeric array) – positions at which to evaluate scalar field
scalar_field (Nx1 float array) – the field to interpolate, defined on either vertices or faces of the current mesh
options (struct with fields) –

imCoordsstr pullback specifier (default=QS.piv.imCoords)
coordinate system of the pullback space

sfLocation’vertices’ or ‘faces’ (default=’vertices’)
evaluation location of the scalar field on currentMesh

iMethod’linear’, ‘cubic’, or ‘nearest’ (default=’linear’)
interpolation method for scalar field

Lxint or float (optional, saves time if supplied)
extent of pullback image space in horizontal dimension

Lyint or float (optional, saves time if supplied)
extent of pullback image space in vertical dimension

Returns:

sf – scalar field evaluated onto pullback coordinates

Return type:

Nx1 float array

@TubULAR.initializeTubULAR(tubi, xp, opts)#

Hidden method for instantiating TubULAR class

Parameters:

tubi (TubULAR object whose properties to fill in with this function) –
xp (Imsane Experiment class instance belonging to tubi or struct) –
if imsane class instance, tubular uses imsane fields to populate metadata. Otherwise this needs to be a struct with fields

fileMeta expMeta

See docs for entries in these fields and their descriptions.
opts (struct with fields) –

flipybool
Set to true if lab coordinates are mirrored along some axis wrt to data coords.

meshDirstr
path to where meshes are stored and output will be placed

timeUnitsstr
units of time, for ex ‘min’

spaceUnitsstr
units of space, for ex ‘$mu$m’ for microns

nUint
resolution in the longitudinal direction, in number of sampling points per proper length of the mesh

nVint
resolution in the circumferential direction, in number of sampling points per circumference

lambdaoptional float
smoothing applied to fields on a mesh before computation of DEC derivatives

lambda_meshoptional float
smoothing (diffusion constant) applied to mesh vertices before computation of DEC fields

lambda_erroptional float
additional smoothing for fields derived from DEC fields

NPMitchell 2020-2022

@TubULAR.alignMaskedDataAPDV#

function alignMaskedDataAPDV(QS) Given training on midgut cells and spurious (amnioserosa) cells as h5 files in dir16bit/stabilized_h5s/, mask the original data and save 3d volumes. Note: if QS.plotting.preview == True, displays intermediate results

THIS IS AT PRESENT UNFINISHED 2020

NPMitchell 2020

@TubULAR.plotMetric(QS, options)#

Return type:: NPMitchell 2020

@TubULAR.generateSPCutMeshStack(QS, spcutMeshStackOptions)#

generateSPCutMeshSmStack(QS, spcutMeshStackOptions): Create TIFF stacks of normally evolved smoothed meshes in (s, phi) coordinate system. Same as generateSPCutMeshSmStack() but for non-smoothed (in time) pullback images.

Parameters:

QS (QuapSlap class instance) –
spcutMeshStackOptions (struct with fields) –
- n_outwardint
  number of steps in positive normal direction to sample
- n_inwardint
  number of steps in negative normal direction to sample
- overwritebool
  whether to overwrite spcutMeshSmStack on disk

Return type:

NPMitchell 2020

@TubULAR.measurePullbackPathlines(tubi, options)#

Measure pathlines of optical flow in pullback space. These pathlines can then be used to query velocities or other properties that are parameterized by pullback coordinates. For example, we average velocites along Lagrangian pathlines. For another example, to build a Lagrangian-frame measure of divergence of tangential velocity, div(v_t), we can query div(v_t) at the coords of the PullbackPathlines (via interpolation of div(v_t) defined on pullback mesh vertices).

NOTE: Updated 2022-07 so that if the pullback is doubleCover of the surface, we ensure that velocities and pathlines of periodic positions respect the periodicity of the surface. We do this by explicitly AVERAGING the displacments of each cover, then clipping to a single cover cylinderCutMesh.

Parameters:

tubi (TubULAR class instance) –
options (struct with fields) –

overwritebool
overwrite previous results

previewbool
view intermediate results

maxIterint (default=100)
#Ricci steps for refMesh.u_ricci and beltrami calculation
disk (Saves to) –
------------- –
refMesh (sprintf(QS.fileName.pathlines.refMesh, t0)) – reference mesh in 3d and its pullback that is conformal mapping to the plane, created using Ricci flow at time t0
piv_pathlines_v3d (v3d: sprintf(QS.fileName.pathlines.v3d, t0)) – Advected PIV evaluation coordinates in XY pixel space make pathlines in XY.
piv_pathlines_v3d – Advect refMesh vertex coordinates in XY pixel space in optical flow to make pathlines in XY (pixel space).
piv_pathlines_v3d – Advect refMesh face barycenters in XY pixel space to make pathlines
piv_pathlines_v3d – Embedding coordinates of advected PIV evaluation coordinates, pushed into APDV coordinate system
piv_pathlines_v3d – Embedding coordinates of advected refMesh vertices, pushed into APDV coordinate system

See also

timeAverageVelocities, options, pullbackPathlines, options, NPMitchell

@TubULAR.plotSPCutMeshSmSeriesUtility(tubi, coordsys, options)#

Using (s,phi) pullback cutmeshes, smooth coord system in time via simple triangular pulse averaging of positions in embedding space

Parameters:

tubi (TubULAR class instance) – Current TubULAR object
coordsys (str specifier ('spcutMeshSm', 'spcutMeshSmRS', 'spcutMeshSmRSC')) – What kind of mesh to plot. Any string without ‘RS’ or ‘RSC’ directs to plotting ‘spcutMeshSmRSC’ Any string with RS directs to plotting ‘spcutMeshSm’ Any string with RSC directs to plotting ‘spcutMeshSmRS’
options (struct with fields) –

overwritebool, default=false
overwrite results on disk

preivewbool, default=QS.plotting.preview
display intermediate results
2020 (NPMitchell) –

@TubULAR.computeLocalSurfacePatch(QS, pts, options)#

Make local surface parameterization enclosing pts in pullback space (could be used for computing angle between two cells during T1s)

Note: additional points can be found in mapped space via newpts = barycentricMap2d(subm.f, subm.u, subm.V2D_scaled_detg, pts)

Parameters:

QS (QuapSlap class instance) –
options (optional struct with fields) –
coordSys : str specifier overwrite : bool buffer : float or int or numeric bufferX : float or int or numeric, trumps buffer for X dim bufferY : float or int or numeric, trumps buffer for Y dim preview : preview scaleByMetric : bool
rescale parameterization to have det(g) approx 1 on either:
faces queried by pts (if both imW and imH are not supplied)

in window com(pts) +/- [imW*0.5, imH*0.5]
scaleByMetricComponentsbool
rescale parameterization to approximate isothermal coordinates

imWnumeric
if scaleByMetric or scaleByMetricComponents, supplying imW then determines the width in X dimension of V2D parameterization in which to query metric for fixing. If none supplied, uses pointLocation of pts in the submesh to probe which faces are queried. If supplied, probeMetricViaFaceLookup = false ; (in units of quasi-embedding space –> ie registered/flattened embedding submesh units, should be close to units of QS.spaceUnits if submesh is not too strongly curved in embedding space).

imHnumeric
if scaleByMetric or scaleByMetricComponents, supplying imH then determines the width in Y dimension of V2D parameterization in which to query metric for fixing. If none supplied, uses pointLocation of pts in the submesh to probe which faces are queried. If supplied, probeMetricViaFaceLookup = false ; (in units of quasi-embedding space –> ie registered/flattened embedding submesh units, should be close to units of QS.spaceUnits if submesh is not too strongly curved in embedding space).

Returns:

subm – f : (#faces x 3 int) face Connectivity List v : (#vertices x 3 float) vertices in embedding space (must be 3D) u : (#vertices x 2 float) original pullback coordinates of submesh Vcut : (#vertices x 3 float) vertices in embedding space (must be 3D)

This is typically the same as v, in which case it is not returned

Fcut(#faces x 3 int) face Connectivity List of cutMesh submesh: This is typically the same as f, in which case it is not returned
Ucut(#vertices x 2 float) original pullback coordinates of submesh: This is typically the same as u, in which case it is not returned
V2D#parameterized: as-rigid-as-possible re-parameterization of the embedding in a local patch of pullback space

V2D_scaled_detg : (returned if scaleByMetric==true) V2D_scaled_g11g22 : (returned if scaleByMetricComponents==true) connectivityCase : int indicator

submesh is within single Cover (most common, disk-like)
–> no problems here, all triangles are oriented correctly

submesh has one edge on another cover, with some long faces,
but is still a topological disk, so we push those vertices reaching across the branch cut to the other side so that triangles are oriented correctly

submesh has faces on multiple covers (topological disk)

submesh spans the whole single Cover (topological annulus)

Return type:

struct with fields

See also

visualizeSegmentationPatch, options, visualizeDemoTracks, options, NPMitchell

@TubULAR.plotPathlineStrainRateTimePoint(QS, tp, options)#

Plot the traceful and traceless components of the strain rate tensor defined at pathlines moving in the domain of parameterization over time. Here, plot the trajectories as immobile in a grid. Uses vertex pathlines.

Parameters:

QS (QuapSlap class instance) –
tp (int) – timepoint in units of (1/QS.timeInterval) * QS.timeUnits
options (struct with fields) –

NPMitchell 2020

@TubULAR.generateUVPrimeCutMeshes(QS, options)#

Generate pullbacks in the (u’,v’) coordinate system

Parameters:

QS –
options (struct with fields) –

Returns:

none
NPMitchell 2020

@TubULAR.plotCutPath(tubi, cutMesh, cutP)#

Plot the cut path of a cylinderCutMesh on the mesh in 3D and save image

Parameters:

tubi (TubULAR class instance object) –
cutMesh (cutMesh object with fields) – v : vertices f : face connectivity list
cutP (optional 3D path of the cut on the cylinderCutMesh) –
2019 (NPMitchell) –

@TubULAR.generateCurrentSPCutMesh(QS, cutMesh, spcutMeshOptions)#

Compute the parameterization (s(u’),phi) where s is the circumferentially-averaged pathlenth along the surface from u=0 to u’ and phi = v + phi0(u) is the cicumferential coordinate twisted by an amount phi0 relative to the uv conformal parameterization. phi0 is chosen for each discrete value of u based on either:

[spcutMeshOptions.phiMethod==’3dcurves’] the geometric position of the circumferential curve in 3d space relative to the previous timepoint (or next timepoint if tubi.t0 > t, where t is the current timepoint in question). In other words, we rotate each slice of the tube to more closely match the geometric position of the analogous slice in a timepoint closer to t0.

[spcutMeshOptions.phiMethod==’texture’] the optical correspondence of each slice in pullback space with a ‘nearby’ slice in pullback space of a previous timepoint (later timepoint if t<tubi.t0). In other words, we match the material position in the circumferential coordinate of a timepoint closer to t0

[spcutMeshOptions.phiMethod==’combined’] first based on the geometric position of the circumferential curve in 3d space relative to the previous timepoint (or next timepoint if t0 >t, where t is the current timepoint in question), THEN additionally apply an optical matching. This is useful if there is a lot of jittery motion of the tissue. Note that you can pass In other words, we rotate each slice of the tube to more closely match the geometric position of the analogous slice in a timepoint closer to t0.

Note that the only output in APDV (spaceUnits) coordinates are: mss, mcline, avgpts, avgpts_ss

spcutMesh.sphi spcutMesh.v spcutMesh.vn spcutMesh.ringpath_ss spcutMesh.radii_from_mean_uniform_rs % from uniform DV sampling spcutMesh.radii_from_avgpts spcutMesh.mss % from uniform DV sampling, also stored in centerline spcutMesh.mcline % from uniform DV sampling, also stored in centerline spcutMesh.avgpts % from uniform DV sampling, also stored in centerline spcutMesh.avgpts_ss % pathlength from uniform sampling, also

stored in centerline, in QS.spaceUnits

spcutMesh.ar % affine scaling in X that minimizes isoareal energy spcutMesh.phi0s % spcutMesh.phi0_fit %

Note that depending on the relative axisorder between data frame (TIFFs) and the APDV frame, you might need to set phi0_sign to true which would flip the sign of phi0 (added rather than subtracted from v in (u,v)). This can be true even if QS.flipy is false for some nonstandard circumstances.

Parameters:

tubi (TubULAR class instance) –

Note that the following properties are used:
tubi.phiMethod = (‘3dcurves’, ‘texture’, ‘combined’) tubi.a_fixed = 2.0
cutMesh (cutMesh struct, optional) – cutMesh with fields
spcutMesh (spcutMesh struct, optional) –
spcutMeshOptions (struct with fields, optional) –

overwritebool
overwrite previous results

save_phi0patchbool
show the relaxation steps of phi0 determination

iterative_phi0bool
iteratively determine phi0 until convergence

smoothingMethodstr specifier (default=’none’)
method for smoothing phi0 wrt AP axis coordinate (ss)

textureAxisOrder : ‘xyz’, ‘yxz’, etc smoothingWidth : int

width of kernel for smoothing of phi0 that takes v->phi=v-phi0. Must be odd if smoothingMethod==’savgol’, does not matter if smoothingMethod==’none’

phi0TextureOptsstruct with fields
lowerboundy : float lowerboundy : float step_phi0tile : float width_phi0tile : float potential_sigmay :

Returns:

spcutMesh (struct with fields) – sphi: equal dv, resampled 3d points based on spcutMesh.sphi0
NPMitchell 2020

@TubULAR.getCurrentData(tubi, adjustIV, varargin)#

IV = getCurrentData(tubi, adjustIV) Load/return volumetric intensity data for current timepoint Note: axis order permutation is applied here upon loading and assignment to self (current tubi instance). Note that the steps executed by this function depends on whether xp is an ImSAnE experiment class or a simple struct belonging to the TubULAR class instance.

Parameters:

tubi (tubi class instance (self)) –
adjustIV (bool (default=true)) – apply the intensity adjustment stored in tubi.data.adjustlow/high

Returns:

IV – volumetric intensity data of current timepoint Note this is also assigned to self as tubi.currentData.IV, so avoiding a duplicate copy in RAM is useful for large data by not returning an argument and instead calling tubi.currentData.IV when needed.

Return type:

#channels x 1 cell array of X*Y*Z intensities arrays

@TubULAR.measureMetricKinematics(QS, options)#

[gdot_apM, HH_apM, divv_apM, veln_apM] = measureMetricKinematics(QS, options): Measure degree of incompressibility of the flow on the evolving surface Out-of-plane motion is v_n * 2H, where v_n is normal velocity and H is mean curvature. In-plane motion considered here is div(v_t) where v_t is tangential velocity on the curved surface. The difference div(v_t) - vn*2H = Tr[g^{-1} dot{g}], which is a measure of isotropic metric change over time (dot is the partial derivative wrt time).

Parameters:

QS (QuapSlap class instance) –
options (struct with fields) –

overwritebool
overwrite previous results

previewbool
view intermediate results

timePointsnumeric 1D array
the timepoints to consider for the measurement. For ex, could choose subset of the QS experiment timePoints

alphaValfloat
the opacity of the heatmap to overlay

invertImagebool
invert the data pullback under the velocity map

Returns:

Saves to disk
————-
NPMitchell 2020

@TubULAR.measureBeltramiCoefficientPullbackToPullback(QS, options)#

Measure the Beltrami using the Ricci flow pullback of a (advected) mesh: at some timepoint t1 relative to the pullback of the same (non-advected) mesh at t0. Since Ricci flow is a conformal map, these Beltramis should equal those of the 3D embedding to 2D pullback comparison performed in QS.measureBeltramiCoefficient().

Parameters:

QS (QuapSlap class instance) –
options (struct with fields) –

t0Pathlinesnumeric (default = QS.t0)
timepoint used to define Lagrangian/material frame for mu measurements to be made with respect to

coordSysstr specifier (‘ricci’, ‘uvprime’ or ‘sp’)
coordinate system in which pathlines are computed

maxIterint
only used if coordSys == ‘ricci’, #ricci iterations

overwrite : bool, overwrite previous results on disk climit : limit for caxis of scalar fields Re(mu) and Im(mu)

See also

QS.measureBeltramiCoefficient, NPMitchell

@TubULAR.generateFastMarchingCenterlines(tubi, cntrlineOptions)#

Extract centerlines for all timepoints in the QuapSlap instance

Parameters:

tubi –

uses in particular
timePoints : meshDir flipy : bool (optional, default is true)

APDV coordinate system is mirrored across XZ wrt data coordinate system XYZ.
Options (struct with fields) –
- overwritebool (optional, default is false)
  whether to overwrite previous results on disk
- exponentfloat (optional, default is 1.0)
  exponent of the distance transform to use as speed
- resfloat (optional, default is 1.0)
  resolution with which to sample data volume to build shortest path
- xx, yy, zzNx1 float or int arrays (optional)
  coordinates of the data volume in which to find the shortest path
- xyzlim3x1 or 3x2 or 1x3 float (required if xx,yy,zz not given)
  extrema (maxima or maxima & minima) of each dimension for mesh extent
- meshAPDVFileNamestr (optional)
  If supplied, plot rotated and scaled centerline with the previously saved APDV aligned meshes matching this filename pattern
- reorient_facesbool (optional, default is true)
  ensure proper face orientation on each mesh (slower but rigorous)
Outputs –
------- –
curves (centerline) – [fullfile(outdir, name) ‘_centerline’ extenstr] ;
timepoint (one for each) – [fullfile(outdir, name) ‘_centerline’ extenstr] ;
in (stored) – [fullfile(outdir, name) ‘_centerline’ extenstr] ;
APDV (centerline curves in RS frame (rotated and scaled to align with) – coordsys), stored in: skel_rs_outfn = [fullfile(outdir, name) ‘_centerline_scaled’ extenstr ]
outDirs (2 x 1 cell array) – The 2 output directories: centerlineOutDir, outdir for images
2020 (NPMitchell) –

@TubULAR.cleanCylMeshes(tubi, cleanCylOptions)#

Parameters:

tubi (TubULAR class instance) –
cleanCylOptions (optional struct with optional fields) – overwrite : bool save_ims : bool

Returns:

Saves to disk
————-
tubi.fileName.aBoundaryDorsalPtsClean ;
tubi.fileName.pBoundaryDorsalPtsClean
[tubi.fileName.aBoundaryDorsalPtsClean(1 (end-3) ‘.txt’] ;)
[tubi.fileName.pBoundaryDorsalPtsClean(1 (end-3) ‘.txt’] ;)
NPMitchell 2020

@TubULAR.coordinateSystemDemo(tubi, options)#

COORDINATESYSTEMDEMO(tubi): Draw coordinate system for currentTime stamp, presentation/publication style. Image for publication/presentation on method & coordinate system Create coordinate system charts visualization using smoothed meshes

Parameters:

tubi (TubULAR class instance) –
options (struct with fields) –
coordSys : char (default=’sphism’, ‘sphi’, ‘uv’) style : ‘curves’ or ‘surface’ exten : ‘.png’ or ‘.pdf’ or ‘.jpg’ interpreter : ‘latex’ or ‘default’

whether to use the latex interpreter for axes labels

normLongitudinal : bool includeCenterline : bool fillHoops : bool

Returns:

<none>
Saves to disk
————-
figure saved to fullfile(tubi.dir.spcutMesh, … – ‘coordinate_system_demo’, ‘coordinate_system_demo_zeta’) ;
figure saved to fullfile(tubi.dir.spcutMesh, … – ‘coordinate_system_demo’, ‘coordinate_system_demo_phi’) ;
NPMitchell 2020-2022

TubULAR#

Coordinate Systems#

APDV coordinate system and Centerline specification#

Properties#

Full Contents:#

Prerequisites#

Saves to disk#

Coordinate Systems#

Properties#

APDV coordinate system and Centerline specification#

Overview Usage#

Coordinate Systems#

Properties#

APDV coordinate system and Centerline specification#

Saves to disk#

Saves to disk#

OUTPUTS#

OUTPUTS#

Indices and tables#