Example for using TubULAR just for surface visualization#
Here we visualize surfaces from active droplets. This is an example for extracting DNA droplet from confocal data, courtesy of Alexandra Tayar & the lab of Zvonimir Dogic (UCSB).
%% EXAMPLE MASTER PIPELINE FOR TIME SERIES DATA (3D + time)
% by NPMitchell & Dillon Cislo
% This is a pipeline to analyze a dynamic tube-like interface in 3D data.
% A tube-like surface is one that is either cylindrical in topology or
% elongated and spherical in topology. If the initial mesh is a spherical
% topology, then two 'endcaps' will be truncated in order to transform it
% to a cylindrical topology.
%% Clear workspace ========================================================
% We start by clearing the memory and closing all figures
clear; close all; clc;
cd path/to/data/
dataDir = cd ;
%% ADD PATHS TO THIS ENVIRONMENT ==========================================
origpath = matlab.desktop.editor.getActiveFilename;
cd(fileparts(origpath))
addpath(fileparts(origpath))
addpath(genpath('utility'))
addpath('DEC')
addpath(fullfile('utility','plotting'))
% go back to the data
cd(dataDir)
%% DEFINE NEW MASTER SETTINGS FOR THIS DATASET ============================
if ~exist('./masterSettings.mat', 'file')
% Metadata about the experiment
stackResolution = [1.1 1.1 2.2] ; % resolution in spaceUnits per pixel. For fused lightsheet data, these numbers should all be equal
nChannels = 2 ; % how many channels is the data (ex 2 for GFP + RFP)
channelsUsed = [1,2] ; % which channels are used for analysis
timePoints = 0:54; % timepoints to include in the analysis
ssfactor = 4 ; % subsampling factor
flipy = false ; % whether the data is stored inverted relative to real position in lab frame
timeInterval = 30 ; % physical interval between timepoints
timeUnits = 's' ; % physical unit of time between timepoints
spaceUnits = '$\mu$m' ; % physical unit of time between timepoints
fn = '202203221330_reg2_T%02d'; % filename string pattern, with timestamp and/or channel to be filled in with sprintf()
set_preilastikaxisorder = 'xyzc' ; % data axis order for subsampled h5 data (ilastik input)
swapZT = 0 ; % whether to swap the z and t dimensions
masterSettings = struct('stackResolution', stackResolution, ...
'nChannels', nChannels, ...
'channelsUsed', channelsUsed, ...
'timePoints', timePoints, ...
'ssfactor', ssfactor, ...
'flipy', flipy, ...
'timeInterval', timeInterval, ...
'timeUnits', timeUnits, ...
'spaceUnits', spaceUnits, ...
'fn', fn,...
'swapZT', swapZT, ...
'set_preilastikaxisorder', set_preilastikaxisorder, ...
'nU', 100, ...
'nV', 100);
disp('Saving masterSettings to ./masterSettings.mat')
if exist('./masterSettings.mat', 'file')
ui = input('This will overwrite the masterSettings. Proceed (Y/n)?', 's') ;
if ~isempty(ui) && (strcmp(ui(1), 'Y') || strcmp(ui(1), 'y'))
save('./masterSettings.mat', 'masterSettings')
loadMaster = false ;
else
disp('Loading masterSettings from disk instead of overwriting')
loadMaster = true ;
end
else
save('./masterSettings.mat', 'masterSettings')
loadMaster = false ;
end
else
loadMaster = true ;
end
if loadMaster
% LOAD EXISTING MASTER SETTINGS
disp('Loading masterSettings from ./masterSettings.mat')
load('./masterSettings.mat', 'masterSettings')
% Unpack existing master settings
stackResolution = masterSettings.stackResolution ;
nChannels = masterSettings.nChannels ;
channelsUsed = masterSettings.channelsUsed ;
timePoints = masterSettings.timePoints ;
ssfactor = masterSettings.ssfactor ;
% whether the data is stored inverted relative to real position
flipy = masterSettings.flipy ;
timeInterval = masterSettings.timeInterval ; % physical interval between timepoints
timeUnits = masterSettings.timeUnits ; % physical unit of time between timepoints
spaceUnits = masterSettings.spaceUnits ; % unit of distance of full resolution data pixels ('$\mu$m')
fn = masterSettings.fn ;
set_preilastikaxisorder = masterSettings.set_preilastikaxisorder ;
swapZT = masterSettings.swapZT ;
nU = masterSettings.nU ;
nV = masterSettings.nV ;
end
dir16bit = fullfile(dataDir) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PART 1: Surface detection using TubULAR's getSurface method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Setup a working directory for the project, where extracted surfaces,
% metadata and debugging output will be stored. Also specifiy the
% directory containing the data.
cd(dir16bit)
dataDir = dir16bit ;
projectDir = dataDir ;
%
% The following file metadata information is required:
% * 'directory' , the project directory (full path)
% * 'dataDir' , the data directory (full path)
% * 'filenameFormat' , fprintf type format spec of file name
% * 'timePoints' , list of itmes available stored as a vector
% * 'stackResolution' , stack resolution in microns, e.g. [0.25 0.25 1]
%
% The following file metadata information is optional:
% * 'imageSpace' , bit depth of image, such as uint16 etc., defined
% in Stack class
% * 'stackSize' , size of stack in pixels per dimension
% [xSize ySize zSize]
% * 'swapZT' , set=1 if time is 3rd dimension and z is 4th
% A filename base template - to be used throughout this script
fileMeta = struct();
fileMeta.dataDir = dataDir;
fileMeta.filenameFormat = [fn, '.tif'];
fileMeta.nChannels = nChannels;
fileMeta.timePoints = timePoints ;
fileMeta.stackResolution = stackResolution;
fileMeta.swapZT = masterSettings.swapZT;
% Set required additional information on the experiment.
%
% The following project metadata information is required:
% * 'channelsUsed' , the channels used, e.g. [1 3] for RGB
% * 'channelColor' , mapping from element in channels used to RGB = 123
%
% The following project meta data information is optional:
% * 'description' , string describing the data set set experiments metadata,
% such as a description, and if the surface is dynamic,
% or requires drift correction of the sample.
% first_tp is also required, which sets the tp to do individually.
first_tp = 1 ;
expMeta = struct();
expMeta.channelsUsed = channelsUsed ;
expMeta.channelColor = 1;
expMeta.description = 'microtubule-kinesin gel';
expMeta.fitTime = fileMeta.timePoints(first_tp);
%% SET DETECTION OPTIONS ==================================================
% Load/define the surface detection parameters
msls_detOpts_fn = fullfile(projectDir, 'msls_detectOpts.mat') ;
if exist(msls_detOpts_fn, 'file') && false
load(msls_detOpts_fn, 'detectOptions')
else
outputfilename_ply='mesh_ls_' ;
outputfilename_ls='ls_' ;
outputfilename_smoothply = 'mesh_' ;
ms_scriptDir='/mnt/data/code/morphsnakes_wrapper/morphsnakes_wrapper/' ;
init_ls_fn = 'msls_initguess' ;
meshlabCodeDir = '/mnt/data/code/meshlab_codes/';
mlxprogram = fullfile(meshlabCodeDir, ...
'laplace_surface_rm_resample30k_reconstruct_LS3_1p2pc_ssfactor4.mlx') ;
prob_searchstr = '_stab_Probabilities.h5' ;
preilastikaxisorder = set_preilastikaxisorder; ... % axis order in input to ilastik as h5s. To keep as saved coords use xyzc
ilastikaxisorder= 'cxyz'; ... % axis order as output by ilastik probabilities h5
imsaneaxisorder = 'xyzc'; ... % axis order relative to mesh axis order by which to process the point cloud prediction. To keep as mesh coords, use xyzc
% Name the output mesh directory --------------------------------------
mslsDir = [fullfile(projectDir, 'msls_output') filesep];
% Surface detection parameters ----------------------------------------
detectOptions = struct( 'channel', 1, ...
'ssfactor', 4, ...
'niter', 35,...
'niter0', 160, ...
'pre_pressure', -5, ...
'pre_tension', 0, ...
'pressure', 0, ...
'tension', 0.5, ...
'post_pressure', 2, ...
'post_tension', 3, ...
'exit_thres', 1e-7, ...
'foreGroundChannel', 1, ...
'fileName', fn, ...
'mslsDir', mslsDir, ...
'ofn_ls', outputfilename_ls, ...
'ofn_ply', outputfilename_ply,...
'ms_scriptDir', ms_scriptDir, ...
'timepoint', timePoints(1), ...
'zdim', 2, ...
'ofn_smoothply', outputfilename_smoothply, ...
'mlxprogram', mlxprogram, ...
'init_ls_fn', init_ls_fn, ... % set to none to load prev tp
'run_full_dataset', projectDir,... % projectDir, ... % set to 'none' for single tp
'radius_guess', 40, ...
'dset_name', 'exported_data',...
'save', true, ... % whether to save images of debugging output
'center_guess', '200,75,75',... % xyz of the initial guess sphere ;
'plot_mesh3d', false, ...
'dtype', 'h5',...
'mask', 'none',...
'mesh_from_pointcloud', false, ...
'prob_searchstr', prob_searchstr, ...
'physicalaxisorder', imsaneaxisorder, ...
'preilastikaxisorder', preilastikaxisorder, ...
'ilastikaxisorder', ilastikaxisorder, ...
'include_boundary_faces', true, ...
'smooth_with_matlab', true); % set this to >0 to use matlab laplacian filter instead of meshlab
% save options
if exist(msls_detOpts_fn, 'file')
disp('Overwriting detectOptions --> renaming existing as backup')
backupfn1 = [msls_detOpts_fn '_backup1'] ;
if exist(backupfn1, 'file')
backupfn2 = [msls_detOpts_fn '_backup2'] ;
system(['mv ' backupfn1 ' ' backupfn2])
end
system(['mv ' msls_detOpts_fn ' ' backupfn1])
end
disp('Saving detect Options to disk')
save(msls_detOpts_fn, 'detectOptions') ;
end
% Overwrite certain parameters for script structure
mslsDir = detectOptions.mslsDir ;
%% Collect options into xp experiment struct
xp = struct('expMeta', expMeta, 'fileMeta', fileMeta, ...
'detectOptions', detectOptions) ;
%% Now we have 3d data volumes and surfaces. Define a TubULAR object.
% To visualize data on these surfaces and compute how these surfaces deform
% we now define TubULAR object.
nU = masterSettings.nU ;
nV = masterSettings.nV ;
% Options about the data
opts.meshDir = mslsDir ; % Directory where meshes reside
opts.flipy = flipy ; % Set to true if data volume axes are inverted in chirality wrt physical lab coordinates
opts.timeInterval = timeInterval ; % Spacing between adjacent timepoints in units of timeUnits
opts.timeUnits = timeUnits ; % units of time, so that adjacent timepoints are timeUnits * timeInterval apart
opts.spaceUnits = spaceUnits ; % Units of space in LaTeX, for ex '$mu$m' for micron
% Options for surface parameterization
opts.nU = nU ; % How many points along the longitudinal axis to sample surface
opts.nV = nV ; % How many points along the circumferential axis to sample surface
opts.t0 = xp.fileMeta.timePoints(1) ; % reference timepoint used to define surface-Lagrangian and Lagrangian measurements
opts.normalShift = 10 ; % Additional dilation acting on surface for texture mapping
opts.a_fixed = 2.0 ; % Fixed aspect ratio of pullback images. Setting to 1.0 is most conformal mapping option.
opts.phiMethod = 'curves3d' ; % Method for following surface in surface-Lagrangian mapping [(s,phi) coordinates]
% Options for surface visualization
opts.adjustlow = 1.00 ; % floor for intensity adjustment
opts.adjusthigh = 99.9 ; % ceil for intensity adjustment (clip)
% Options for extraction of tissue kinematics
opts.lambda_mesh = 0.00 ; % Smoothing applied to the mesh before DEC measurements
opts.lambda = 0.0 ; % Smoothing applied to computed values on the surface
opts.lambda_err = 0.0 ; % Additional smoothing parameter, optional
% Now we're ready!
disp('defining TubULAR class instance (tubi= tubular instance)')
tubi = TubULAR(xp, opts) ;
disp('done defining TubULAR instance')
%% CREATE THE SUBSAMPLED H5 FILE FOR INPUT TO ILASTIK =====================
tubi.prepareIlastik()
%% TRAIN NON-STABILIZED DATA IN ILASTIK TO IDENTIFY SURFACE ==============
% Open ilastik, train on h5s until probabilities and uncertainty are
% satisfactory for extracting a mesh. For example, here we train on the
% membrane (channel 1) and the yolk (channel 2), so that a level set will
% enclose the yolk but not escape through the membrane.
%% Create MorphSnakes LevelSets from the Probabilities output of ilastik ==
% Now detect all surfaces
disp('Surfaces are extracted serially using active contours')
overwrite = false ;
tubi.getMeshes(overwrite);
%% Inspect all meshes in 3D
for tp = xp.fileMeta.timePoints
% Load the mesh
meshfn = sprintf( tubi.fullFileBase.mesh, tp ) ;
mesh = read_ply_mod(meshfn) ;
% Plot the mesh in 3d. Color here by Y coordinate
trisurf(mesh.f, mesh.v(:, 1), mesh.v(:, 2), mesh.v(:, 3), ...
mesh.v(:, 3), 'edgecolor', 'none', 'Facealpha', 0.5)
title(['t=' num2str(tp)])
xlabel('x')
ylabel('y')
zlabel('z')
axis equal
view(2)
pause(0.5)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PART 2: TubULAR -- surface parameterization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Obtain APDV coordinates of the surface.
% There are two options for obtaining these coordinates.
% 1. Automatically determine A and P by the extremal points of the
% surface mesh along the elongated axis of the mesh, and define DV as
% pointing perpendicular to this.
% 2. Train in iLastik for an anterior spot in 3d (A), a posterior spot in
% 3d (P), and a spot which is dorsal to the line connecting A and P. Any
% dorsal point is fine, as long as it points dorsal to the AP axis
% defined by A and P. See picture below.
%
% example:
% D
% |
% |
% A -----------------P
%
% Here we use option 2. We must prepare APDV ilastik training first outside
% MATLAB.
% Train on anterior (A), posterior (P), background (B), and
% dorsal anterior (D) location in different iLastik channels by having a
% blob centered on a point that you wish to identify as each label (in 3D).
% anteriorChannel, posteriorChannel, and dorsalChannel specify the iLastik
% training channel that is used for each specification.
% Name the h5 file output from iLastik as ..._Probabilities_apcenterline.h5
% Training for dorsal (D) is only needed at the reference time point, t0,
% because that's the only one that's used.
%
% A dorsal blob for the gut is marked at the site where the gut closes,
% with 48YGAL4-expressing 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.
%% Define global orientation frame (for viewing in canonical frame)
% Compute APDV coordinate system
alignAPDVOpts = struct() ;
alignAPDVOpts.overwrite = false ;
tubi.computeAPDVCoords(alignAPDVOpts) ;
%% PLOT ALL TEXTURED MESHES IN 3D (OPTIONAL: this is SLOW) ================
% Establish texture patch options
metadat = struct() ;
metadat.reorient_faces = false ; % set to true if some mesh normals may be inverted (requires gptoolbox if true)
metadat.normal_shift = tubi.normalShift ; % normal push, in pixels, along normals defined in data XYZ space
metadat.texture_axis_order = [1 2 3] ; % texture space sampling. If the surface and dataspace have axis permutation, enter that here
metadat.blackFigure = true ;
metadat.channel = 2;
metadat.perspective_angle = [70,-50] ;
Options.PSize = 5 ; % Psize is the linear dimension of the grid drawn on each triangular face. Set PSize > 1 for refinement of texture on each triangle of the surface triangulation. Higher numbers are slower but give more detailed images.
Options.numLayers = [10, 2]; % how many layers to MIP over/bundle into stack, as [outward, inward]
Options.layerSpacing = 1 ; % Distance between layers over which we take MIP, in pixels,
metadat.plot_dorsal = false ;
metadat.plot_left = false ;
metadat.plot_right = false ;
metadat.plot_ventral = false ;
metadat.plot_perspective = true ;
metadat.overwrite = true ;
% Plot on surface for all timepoints
tubi.plotSeriesOnSurfaceTexturePatch(metadat, Options)
%%
for tp = tubi.xp.fileMeta.timePoints
tubi.setTime(tp)
tubi.maskCurrentDataWithMesh()
end