https://github.com/npmitchell/VisceralOrganMorphogenesisViaCalciumPatternedMuscleConstrictions
Raw File
Tip revision: 3835e0a1faa5853939c2967c867dde5cf49d9887 authored by Noah Mitchell on 25 April 2022, 18:57:35 UTC
added code for figures and dependency code
Tip revision: 3835e0a
figure2_cellSegmentation_and_Tracking.m
%% GUT MASTER PIPELINE FOR TIME SERIES DATA
% by NPMitchell & Dillon Cislo
%
% This is a pipeline to take the surface of the growing Drosophila gut and
% conformally map patches of it to the unit disk

% ============= %
%  CONVENTIONS
% ============= % 
% >> Data and chirality <<
% ------------------------
% xp.stack.image.apply() gives imageLength, imageWidth, imageDepth as XYZ.
% This forms our definition for XYZ. However, this is mirrored wrt the
% physical coordinates in the lab frame. 
% Note that this is different from bfGetReader(fn), which gives yxz -- ie
%   imageWidth, imageLength, imageDepth. In Experiment.m, y and x are 
% swapped to account for this discrepancy after using bfGetReader to build
% the empty array.
%
% Python will read h5s with dimensions backward wrt MATLAB, so if h5s are 
% saved as cxyz in MATLAB, morphsnakes/python will read as zyxc. 
% preilastikaxisorder converts from xyzc to what is given for h5. 
% Typically, ilasik puts the channel at the start in MATLAB (end in python)
% Therefore, set axisorder = 'cxyz' since MATLAB probability h5s are saved 
% as cxyz.
% Note: on 2021-01-30, I set physical = preilastik = ilastik = xyzc, and 
% I get flipy=0 and texture_axis_order = [3 2 1] ; Is this true in general?
%
%
% THIS WORKS BETTER:
% ------------------
% set_preilastikaxisorder = 'xyzc' ;
% 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
%
%    'preilastikaxisorder', preilastikaxisorder, ... 
%    'ilastikaxisorder', ilastikaxisorder, ... 
%    'physicalaxisorder', imsaneaxisorder, ... 
%
% Meshes from integralDetector are stored with normals 'outward'/basal.
% It is important to preserve the triangle orientation for texturepatch,
% but the mesh XY coordinates are flipped with respect to the physical lab
% coordinates. Therefore, we set flipy = true, and the aligned meshes have 
% apical normals. The APDV meshes are black in meshlab when viewed from the
% outside.  flipy indicates that the APDV coordinate system is mirrored 
% across XZ wrt data coordinate system XYZ. 

% >> APDV ilastik training <<
% ---------------------------
% Train on anterior (A), posterior (P), background (B), and 
% dorsal anterior (D) location in different iLastik channels.
% Default order is 1-A, 2-P, 3-bg, 4-D. 
% Train on pre-stabilized H5s, NOT on stabilized H5s. 
% anteriorChannel, posteriorChannel, and dorsalChannel specify the iLastik
% training channel that is used for each specification.
% By default, name the h5 file output from iLastik as 
%    QS.fullFileBase.apdProb --> [QS.fileBase.name '_Probabilities_APDVcoords.h5'] 
% and/or 
%    QS.fullFileBase.apCenterlineProb --> [QS.fileBase.name '_Probabilities_apcenterline.h5'] 
% 
% The APDV coordinate system -- this is where global rot is defined -- is
% computed via 
%    apdvOpts.aProbFileName = QS.fullFileBase.apdProb ;
%    apdvOpts.pProbFileName = QS.fullFileBase.apCenterlineProb ;
%    QS.computeAPDVCoords(alignAPDVOpts) ;
% However, the APD COMs for centerline extraction can be independent. For
% example, if the original segmentation used a third channel which happens
% to be the posterior end, we can specify the posterior COM (P) via:
%    apdvOpts.aProbFileName = QS.fullFileBase.apdProb ; % filename pattern for the apdv training probabilities
%    apdvOpts.pProbFileName = QS.fullFileBase.prob ;
%    apdvOpts.posteriorChannel = 3 ;
%    [acom_sm, pcom_sm] = QS.computeAPDCOMs(apdvOpts) ;
%
% >> Centerline extraction <<
% ---------------------------
% STEP 1: INITIAL CENTERLINE
% Centerlines are first computed using a fast marching procedure on a 
% distance transform (DT) of the mesh data in dataspace. The starting and
% endingpoints of this initial centerline can be specified using APDVcoords
% training (and the resulting COM extraction) from the previous section or
% else can be unique to the centerline extraction. This independence is
% useful, for example, in convoluted organs like the gut, where the 
% centerline of the muscle data curves away from the AP axis in the
% posterior where the midgut turns back on itself to connect to the
% hindgut.
%
% For redundant use of APDV coords for centerline, use the option:
%
% For independent training use the following convention:
% Train on anterior (A), posterior (P), background (B), and 
% dorsal anterior (D) location in different iLastik channels.
% Default order is 1-A, 2-P, 3-bg, 4-D. 
% 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
%
% STEP 2: PRECISE CENTERLINE
% The initial centerline is used to constrain the winding number of the
% orbifold mapping. The second centerline is a series of average 3d points,
% one for each DV ring of the quasi-axisymmetric pullback. 
%

%% Clear workspace ========================================================
% We start by clearing the memory and closing all figures
clear; close all; clc;
% change this path, for convenience
cd /mnt/data/48Ygal4-UAShistRFP/201904031830_great/Time4views_60sec_1p4um_25x_1p0mW_exp0p35_2/data/
% cd /mnt/crunch/48YGal4UasLifeActRuby/201904021800_great/Time6views_60sec_1p4um_25x_1p0mW_exp0p150_3/data/
% cd /mnt/data/48YGal4UasLifeActRuby/201902201200_unusualfolds/Time6views_60sec_1p4um_25x_obis1_exp0p35_3/data/
% cd /mnt/data/48Ygal4UASCAAXmCherry/201902072000_excellent/Time6views_60sec_1p4um_25x_obis1p5_2/data 
% cd /mnt/data/48Ygal4UASCAAXmCherry/201903211930_great/Time6views_60sec_1p4um_25x_1p0mW_exp0p150/
% cd /mnt/crunch/48Ygal4UASsqhGFP/201902271940_excellent_notunpacked/Time6views_60sec_1p4um_25x_1p5mW_exp1p0_3/data/
% cd /mnt/data/mef2GAL4klarUASCAAXmChHiFP/202003151700_1p4um_0p5ms3msexp/data/
% cd /mnt/data/mef2GAL4klarUASCAAXmChHiFP/202003151700_1p4um_0p5ms3msexp/data/
% cd /mnt/data/antpGAL4UASCAAXmChHGFP/202103281352_1p4um_0p15ms0p25ms_1mW1mW_GFPRFP/Time3views_180s/data/
% cd /mnt/data/UbxGAL4UASCAAXmChHGFP/202105132247_UbxG4kCAAXHGFP_1p2um_0p1ms0p2ms_1mW1mW_3v300s/data


% .=========.
% |  VIP10  |
% .=========.
% cd /mnt/crunch/gut/48YGal4UasLifeActRuby/201907311600_48YGal4UasLifeActRuby_60s_exp0p150_1p0mW_25x_1p4um
% cd /mnt/crunch/gut/48YGal4klarUASCAAXmChHiFP/202001221000_60sec_1p4um_25x_1mW_2mW_exp0p25_exp0p7/Time3views_1017/data/
% cd /mnt/crunch/gut/Mef2Gal4klarUASCAAXmChHiFP/202003151700_1p4um_0p5ms3msexp/Time3views_1/data/
% cd /mnt/crunch/gut/Mef2Gal4klarUASCAAXmChHiFP/202007151930_1p4um_0p5msexp/Time3views_25x_60s/data/
% cd /mnt/data/gut/AntpGal4OCRL/202104142135_antpG4kOCRL_1p4um_0p3ms1ms_1mW1mW_GFPRFP/Time2views_60s_RFPRFP_2137_tp0at2132/data/

dataDir = cd ;

%% PATHS ==================================================================
origpath = matlab.desktop.editor.getActiveFilename;
cd(fileparts(origpath))
aux_paths_and_colors
cd(dataDir)

%% Global options
% Decide whether to change previously stored detection Opts, if they exist
overwrite_masterSettings = false ;
overwrite_mips = false ;
overwrite_detOpts = false ;
run_full_dataset_ms = false ;
overwrite_alignAPDVOpts = false ;
overwrite_APDVCOMs = false ;
overwrite_APDVMeshAlignment = false ;
overwrite_alignedMeshIms = false ;
overwrite_centerlines = false ;
overwrite_centerlineIms = false ;
overwrite_TextureMeshOpts = false ;
overwrite_endcapOpts = false ;
overwrite_idAnomClines = false ;
overwrite_cleanCylMesh = false ;

%% DEFINE NEW MASTER SETTINGS
if overwrite_masterSettings || ~exist('./masterSettings.mat', 'file')
    % Metadata about the experiment
    stackResolution = [.2619 .2619 .2619] ;
    nChannels = 1 ;
    channelsUsed = 1 ;
    timePoints = 1:234; %86:211 ;
    ssfactor = 4 ;
    % whether the data is stored inverted relative to real position
    flipy = true ; 
    timeInterval = 1 ;  % physical interval between timepoints
    timeUnits = 'min' ; % physical unit of time between timepoints
    spaceUnits = '$\mu$m' ; % physical unit of time between timepoints
    scale = 0.04 ;      % scale for conversion to 16 bit
    % file32Base = 'TP%d_Ch0_Ill0_Ang0,45,90,135,180,225,270,315.tif'; 
    file32Base = 'TP%d_Ch0_Ill0_Ang0,60,120,180,240,300.tif'; 
    % file32Base = 'TP%d_Ch0_Ill0_Ang0,60,120,180,240,300.tif'; 
    fn = 'Time_%06d_c1_stab';
    fn_prestab = 'Time_%06d_c1.tif';
    set_preilastikaxisorder = 'xyzc' ;
    swapZT = 1 ;
    masterSettings = struct('stackResolution', stackResolution, ...
        'nChannels', nChannels, ...
        'channelsUsed', channelsUsed, ...
        'timePoints', timePoints, ...
        'ssfactor', ssfactor, ...
        'flipy', flipy, ...
        'timeInterval', timeInterval, ...
        'timeUnits', timeUnits, ...
        'spaceUnits', timeUnits, ...
        'scale', scale, ...
        'file32Base', file32Base, ...
        'fn', fn,...
        'fn_prestab', fn_prestab, ...
        'swapZT', swapZT, ...
        'set_preilastikaxisorder', set_preilastikaxisorder, ...
        't0_for_phi0', 110, ... % 40 for mef2 single channel, 110 for CAAX excellent
        'nU', 100, ...  % 150 for mef2 data with posterior midgut loop
        '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')
    scale = masterSettings.scale ;      % scale for conversion to 16 bit
    file32Base = masterSettings.file32Base ; 
    fn = masterSettings.fn ;
    fn_prestab = masterSettings.fn_prestab ;
    set_preilastikaxisorder = masterSettings.set_preilastikaxisorder ;
    swapZT = masterSettings.swapZT ;
    t0_for_phi0 = masterSettings.t0_for_phi0 ;
    nU = masterSettings.nU ;
    nV = masterSettings.nV ;
end
dir32bit = fullfile(dataDir, 'deconvolved_32bit') ;
dir16bit = fullfile(dataDir, 'deconvolved_16bit') ;
dir16bit_prestab = fullfile(dir16bit, 'data_pre_stabilization') ;

%% END OF EXPERIMENT METADATA =============================================
% =========================================================================
% =========================================================================
%% -IIV. make MIPs for 32bit images
% Skip if already done
mipDir = fullfile(dir32bit, 'mips32bit') ;
Options.overwrite_mips = overwrite_mips ;
Options.scale = scale ;
makeMips(timePoints, dir32bit, file32Base, mipDir, Options)

%%  -IV. convert 32 to 16bit images
% Skip if already done
convert32to16bit(timePoints, scale, dir32bit, dir16bit_prestab,...
    file32Base, fn_prestab)

% % Rename stab to prestab -- hack
% % fns = fullfile('./deconvolved_16bit/Time*stab')
% % for qq=1:length(fns)
% %     command = ['mv ' fullfile(fns.folder, fns.name) fullfile(dir16bit, fns.name
% % end

% -III. make MIPs for 16bit images
% Skip if already done
mipDir = fullfile(dir16bit_prestab, 'mips') ;
Options.overwrite_mips = false ;
Options.scale = -1 ; % do NOT rescale intensities during intensity projection
makeMips(timePoints, dir16bit_prestab, fn_prestab, mipDir, Options)

%  -II. stabilize images, based on script stabilizeImagesCorrect.m
% Skip if already done
% name of directory to check the stabilization of mips
mips_stab_check = fullfile(mipDir, 'stab_check') ;
mipoutdir = fullfile(mipDir, 'mips_stab') ;
% Choose bit depth as typename
typename = 'uint16' ;
% Give file names for I/O
fileNameIn = fullfile(dir16bit_prestab, fn_prestab) ;
fileNameOut = fullfile(dir16bit, [fn '.tif']) ;
rgbName = [fn '.png'] ;
typename = 'uint16' ; 
Options = struct(); 
Options.overwrite_mips = false ;
Options.overwrite_tiffs = false ;
Options.im_intensity = 1 ; % 0.01 ;
Options.imref_intensity = 1 ; % 0.005 ;
stabilizeImages(fileNameIn, fileNameOut, rgbName, typename, ...
    timePoints, timePoints, timePoints(50), ...
    mipDir, mipoutdir, mips_stab_check, Options)


%%   -I. master_gut_timeseries_prestab_for_training.m
% Skip if already done
cd(dir16bit)
dataDir = cd ;
masterSettings.dir16bit_prestab = dir16bit_prestab ;
makeH5SeriesPrestabForTraining(masterSettings)

%% 
cd(dir16bit)

%% I. INITIALIZE ImSAnE PROJECT ===========================================
% 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 = cd ;
projectDir = dataDir ;
% [ projectDir, ~, ~ ] = fileparts(matlab.desktop.editor.getActiveFilename); 
cd(projectDir);
if projectDir(end) ~= '/'
    projectDir = [projectDir '/'];
end

% Start by creating an experiment object, optionally pass on the project
% directory (otherwise it will ask), and change into the directory of the
% data.  This serves as a front-end for data loading, detection, fitting
% etc.
xp = project.Experiment(projectDir, dataDir);

% Set file and experiment meta data
%
% We assume on individual image stack for each time point, labeled by time.
%  To be able to load the stack, we need to tell the project wehre the data
%  is, what convention is assumed for the file names, available time
%  points, and the stack resolution.  Options for modules in ImSAnE are
%  organized in MATLAB structures, i.e a pair of field names and values are
%  provided for each option.
%
% 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. A verbal data set
% description, Jitter correct by translating  the sample, which time point
% to use for fitting, etc.
%
% 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
% * 'dynamicSurface'    , Not implemented yet, future plan: boolean, false: static surface
% * 'detectorType'      , name of detector class, e.g. radielEdgeDetector
%                         ,(user threshholded), fastCylinderDetector
% * 'fitterType'        , name of fitter class
%
% 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.
% * 'jitterCorrection', Boolean, false: No fft based jitter correction 

% 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         = 'Drosophila gut';
expMeta.dynamicSurface      = 1;
expMeta.jitterCorrection    = 0;  % 1: Correct for sample translation
expMeta.fitTime             = fileMeta.timePoints(first_tp);
expMeta.detectorType        = 'surfaceDetection.morphsnakesDetector';
expMeta.fitterType          = 'surfaceFitting.meshWrapper';

% Now set the meta data in the experiment.
xp.setFileMeta(fileMeta);
xp.setExpMeta(expMeta);
xp.initNew();

clear fileMeta expMeta

%% LOAD THE FIRST TIME POINT ==============================================
xp.setTime(xp.fileMeta.timePoints(1)) ;
% xp.loadTime(xp.fileMeta.timePoints(first_tp));
% xp.rescaleStackToUnitAspect();

%% SET DETECT OPTIONS =====================================================
% Must run this section for later functionality.
% Mesh extraction options
if run_full_dataset_ms
    run_full_dataset = projectDir ; 
else
    run_full_dataset = 'none' ;
end

% Load/define the surface detection parameters
msls_detOpts_fn = fullfile(projectDir, 'msls_detectOpts.mat') ;
if exist(msls_detOpts_fn, 'file') && ~overwrite_detOpts
    load(msls_detOpts_fn, 'detectOptions')
else
    channel = 1;
    foreGroundChannel = 1;
    % ssfactor = 4;
    % niter = 15 ;
    % niter0 = 15 ;
    % ofn_ply = 'mesh_' ; % mesh_apical_ms_stab_' ; 
    % ofn_ls = 'msls_' ; % mesh_apical_stab_' ;
    % ofn_smoothply = 'mesh_apical_stab_' ; % mesh_apical_stab_' ;
    % lambda1 = 1 ;
    % lambda2 = 1 ;
    % exit_thres = 0.0001 ;
    % if contains(projectDir, 'LifeActRuby') || contains(projectDir, 'CAAXmCherry') 
    %     nu = 4 ;  % For volumetric
    %     smoothing = 0.20 ;
    % else
    %     nu = 0.00 ;
    %     smoothing = 0.10 ;
    % end
    % pre_nu = -4 ;
    % pre_smoothing = 0 ;
    % post_nu = 3 ;
    % post_smoothing = 1 ;
    zdim = 2 ;
    
    % for caax_great
    ssfactor=4
    niter=35
    niter0=35;
    ofn_ply='mesh_ms_'
    ofn_ls='msls_'
    ofn_smoothply = 'mesh_stab_' ; % mesh_apical_stab_' ;
    ms_scriptDir='/mnt/data/code/morphsnakes_wrapper/morphsnakes_wrapper/'
    pre_nu=-5
    pre_smoothing=0
    lambda1=1
    lambda2=1
    exit_thres=0.000005
    smoothing=0.20
    nu=0
    post_nu=2
    post_smoothing=5
    
    init_ls_fn = 'msls_initguess.h5';
    % mlxprogram = fullfile(meshlabCodeDir, ...
    %     'laplace_refine_HCLaplace_LaplaceSPreserve_QuadEdgeCollapse60kfaces.mlx') ;
    mlxprogram = fullfile(meshlabCodeDir, ...
        'laplace_surface_rm_resample30k_reconstruct_LS3_1p2pc_ssfactor4.mlx') ;
    radius_guess = 40 ;
    center_guess = '200,75,75' ;
    dtype = 'h5' ;
    mask = 'none' ;
    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
    include_boundary_faces = true ;
    smooth_with_matlab = -1;
    
    % Name the output mesh directory ------------------------------------------
    % msls_exten = ['_prnu' strrep(strrep(num2str(pre_nu, '%d'), '.', 'p'), '-', 'n')];
    % msls_exten = [msls_exten '_prs' strrep(num2str(pre_smoothing, '%d'), '.', 'p') ];
    % msls_exten = [msls_exten '_nu' strrep(num2str(nu, '%0.2f'), '.', 'p') ];
    % msls_exten = [msls_exten '_s' strrep(num2str(smoothing, '%0.2f'), '.', 'p') ];
    % msls_exten = [msls_exten '_pn' num2str(post_nu, '%d') '_ps',...
    %     num2str(post_smoothing)];
    % msls_exten = [msls_exten '_l' num2str(lambda1) '_l' num2str(lambda2) ];
    meshDir = [projectDir 'msls_output' filesep];
    % meshDir = [meshDir msls_exten '/'] ;

    % Surface detection parameters --------------------------------------------
    detectOptions = struct( 'channel', channel, ...
        'ssfactor', ssfactor, ...
        'niter', niter,...
        'niter0', niter0, ...
        'lambda1', lambda1, ...
        'lambda2', lambda2, ...
        'nu', nu, ...
        'smoothing', smoothing, ...
        'post_nu', post_nu, ...
        'post_smoothing', post_smoothing, ...
        'exit_thres', exit_thres, ...
        'foreGroundChannel', foreGroundChannel, ...
        'fileName', sprintf( fn, xp.currentTime ), ...
        'mslsDir', meshDir, ...
        'ofn_ls', ofn_ls, ...
        'ofn_ply', ofn_ply,...
        'ms_scriptDir', ms_scriptDir, ...
        'timepoint', xp.currentTime, ...
        'zdim', zdim, ...
        'ofn_smoothply', ofn_smoothply, ...
        'pre_nu', pre_nu, ...
        'pre_smoothing', pre_smoothing, ...
        'mlxprogram', mlxprogram, ...
        'init_ls_fn', init_ls_fn, ... % set to none to load prev tp
        'run_full_dataset', run_full_dataset,... % projectDir, ... % set to 'none' for single tp
        'radius_guess', radius_guess, ...
        'dset_name', 'exported_data',...
        'center_guess', center_guess,... % xyz of the initial guess sphere ;
        'save', true, ... % whether to save images of debugging output
        'plot_mesh3d', false, ...
        'dtype', dtype,...
        'mask', mask,...
        'mesh_from_pointcloud', false, ...
        'prob_searchstr', prob_searchstr, ...
        'preilastikaxisorder', preilastikaxisorder, ... 
        'ilastikaxisorder', ilastikaxisorder, ... 
        'physicalaxisorder', imsaneaxisorder, ... 
        'include_boundary_faces', include_boundary_faces, ...
        'smooth_with_matlab', smooth_with_matlab, ...
        'pythonVersion', 2) ;

    % 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
detectOptions.fileName = sprintf( fn, xp.currentTime ) ;
detectOptions.run_full_dataset = run_full_dataset ;
detectOptions.ms_scriptDir = ms_scriptDir ;
meshDir = detectOptions.mslsDir ;
% These are now in QS.
% meshFileBase = [ofn_smoothply '%06d'] ;
% alignedMeshBase = [ofn_smoothply '%06d_APDV_um'] ;

% Set detect options ------------------------------------------------------
xp.setDetectOptions( detectOptions );
disp('done')

%% Dictionary for updating detectOptions
% detectOptions.pressure = detectOptions.nu ;
% detectOptions.tension = detectOptions.smoothing ;
% detectOptions.pre_pressure = detectOptions.pre_nu ;
% detectOptions.post_pressure = detectOptions.post_nu ;
% detectOptions.pre_tension = detectOptions.pre_smoothing ;
% detectOptions.post_tension = detectOptions.post_smoothing ;
% 
% detectOptions = rmfield(detectOptions, 'nu') ;
% detectOptions = rmfield(detectOptions, 'pre_nu') ; 
% detectOptions = rmfield(detectOptions, 'post_nu') ;
% detectOptions = rmfield(detectOptions, 'smoothing') ;
% detectOptions = rmfield(detectOptions, 'pre_smoothing') ;
% detectOptions = rmfield(detectOptions, 'post_smoothing') ;

%% CREATE THE SUBSAMPLED H5 FILE FOR INPUT TO ILASTIK =====================
% skip if already done

for tt = xp.fileMeta.timePoints
    if ~exist(fullfile(projectDir, 'stabilized_h5s', [sprintf(fn, tt) '.h5']), 'file')
        disp(['Did not find file: ', fullfile(projectDir, 'stabilized_h5s', [sprintf(fn, tt) '.h5'])])
        xp.loadTime(tt);
        xp.rescaleStackToUnitAspect();
        % make a copy of the detectOptions and change the fileName
        detectOpts2 = detectOptions ;
        detectOpts2.fileName = sprintf( fn, xp.currentTime ) ;
        xp.setDetectOptions( detectOpts2 );
        xp.detector.prepareIlastik(xp.stack);
        disp(['done outputting downsampled data h5: tp=' num2str(tt) ' for surface detection'])
    else
        disp(['h5 ' num2str(tt) ' was already output, skipping...'])
    end
end    
disp('Open with ilastik if not already done')


%if ~exist([projectDir sprintf(fn, xp.currentTime) '.h5'], 'file')
 %   xp.detector.prepareIlastik(xp.stack);
  %  disp('done outputting downsampled data h5 for surface detection')
%else
 %   disp('h5 was already output, skipping...')
%end
%disp('Open with ilastik if not already done')

%% TRAIN NON-STABILIZED DATA IN ILASTIK TO IDENTIFY APICAL/YOLK ===========
% Skip if already done.
% Open ilastik, train pre-stab h5s until probabilities and uncertainty are 
% satisfactory, then run on stab images.

%% Note that old h5's used to be different order. To convert, do
% Skip 
% if false
%     tmp = h5read(fullfile(meshDir, init_ls_fn), '/implicit_levelset');
%     % OLD ORDER: yxz for implicit levelset
%     tmp2 = permute(tmp, [3,2,1]);
%     h5create(fullfile(meshDir, init_ls_fn), '/implicit_levelset', size(tmp2), 'Datatype', 'int8')
%     h5write(fullfile(meshDir, init_ls_fn), '/implicit_levelset', tmp2)
% end

%% Create MorphSnakesLevelSet from the Probabilities from ilastik ========
% Skip if already done
% Now detect all surfaces
detectOptions.run_full_dataset = 'none' ;  % projectDir ; % 'none' ;  % override here
if strcmp(detectOptions.run_full_dataset, projectDir)
    % assert(run_full_dataset_ms)
    disp('Running dataset mode')
    xp.setTime(xp.fileMeta.timePoints(1));
    detectOpts2 = detectOptions ;
    detectOpts2.fileName = sprintf( fn, xp.currentTime ) ;
    detectOpts2.nu = 4 ;
    detectOpts2.niter0 = 5 ;
    xp.setDetectOptions( detectOpts2 );
    xp.detectSurface();
else
    assert(~run_full_dataset_ms)
    assert(strcmp(detectOptions.run_full_dataset, 'none'))
    % Morphosnakes for all remaining timepoints INDIVIDUALLY ==============
    for tp = xp.fileMeta.timePoints(100:end)
        try
            xp.setTime(tp);
            % xp.loadTime(tp) ;
            % xp.rescaleStackToUnitAspect();

            % make a copy of the detectOptions and change the fileName
            detectOpts2 = detectOptions ;
            % detectOpts2.post_smoothing = 1 ;
            detectOpts2.timepoint = xp.currentTime ;
            detectOpts2.fileName = sprintf( fn, xp.currentTime );
            % detectOpts2.mlxprogram = fullfile(meshlabCodeDir, ...
            %      'surface_rm_resample30k_reconstruct_LS3_1p2pc_ssfactor4') ;
            %  _octree12.mlx') ;
            % detectOpts2.mlxprogram = fullfile(meshlabCodeDir, ...
            %     'laplace_surface_rm_resample25k_reconstruct_LS3_1p2pc_ssfactor4_vip8test.mlx') ;
            detectOpts2.mlxprogram = fullfile(meshlabCodeDir, ...
                 'laplace_surface_rm_resample25k_reconstruct_LS3_wu13_ssfactor4.mlx') ;
            xp.setDetectOptions( detectOpts2 );
            xp.detectSurface();
            % For next time, use the output mesh as an initial mesh
            detectOpts2.init_ls_fn = 'none' ;
        catch
        %    error('here')
            disp('Could not create mesh -- skipping for now')
        %     % On next timepoint, use the tp previous to current time
        %     detectOptions.init_ls_fn = [detectOptions.ofn_ls, ...
        %             num2str(tp - 1, '%06d' ) '.' detectOptions.dtype] ;
        end
    end
end

%% Define QuapSlap object
nU = masterSettings.nU ;
nV = masterSettings.nV ;
opts.meshDir = meshDir ;
opts.flipy = flipy ;
opts.timeInterval = timeInterval ;
opts.timeUnits = timeUnits ;
opts.spaceUnits = spaceUnits ;
opts.nU = nU ;
opts.nV = nV ;
opts.normalShift = 10 ;
opts.a_fixed = 2.0 ;
opts.adjustlow = 1.00 ;         % floor for intensity adjustment
opts.adjusthigh = 99.9 ;        % ceil for intensity adjustment (clip)
opts.phiMethod = 'curves3d' ;
% for caax
opts.lambda_mesh = 0.00 ;
opts.lambda = 0.0 ;
opts.lambda_err = 0.0 ;

% for histRFP:
opts.nmodes = 7 ;
opts.zwidth = 1 ;
opts.lambda_mesh = 0.00 ;
opts.lambda = 0.01; % should use 0.01 for fig2/3 eLife, I think
opts.lambda_err = 0.01; %shoudl use 0.01 for fig2/3 eLife ;

disp('defining QS')
QS = QuapSlap(xp, opts) ;
disp('done')

%% Make some mips of shallow stacks
adjustIV = false ; 
% QS.makeMIPs(1, {400:490, 510:550, 570:630, 800:880}, [], adjustIV)
stacks = {400:490, 510:550, 570:630, 800:880} ;
QS.makeMIPs(1, stacks, [], adjustIV)

%% HACK -- transpose meshes
% % xyzc -> zyxc
% for tp = QS.xp.fileMeta.timePoints
% 
%     meshfn = sprintf(QS.fullFileBase.mesh, tp) ;  
%     msh = read_ply_mod(meshfn) ;
%     msh.v = msh.v(:, [2,3,1]) ;
%     msh.vn = msh.vn(:, [2,3,1]) ;
%     plywrite_with_normals(meshfn, msh.f, msh.v, msh.vn)
% end


%% Inspect a single mesh
% Skip if already done
tp = 100 ;
meshfn = sprintf(QS.fullFileBase.mesh, tp) ;  
QS.setTime(tp)
QS.getCurrentData() ;
IV = QS.currentData.IV ;
mesh = read_ply_mod(meshfn) ;
for leaf=1:40:size(IV{1}, 1)
    tmp = mesh.v ;
    inds = find(abs(tmp(:, 1) - leaf) < 5) ;
    imshow(imadjust(squeeze(IV{1}(leaf, :, :))'))
    if any(inds)
        hold on;
        % plot(mesh.v(inds, 2), mesh.v(inds, 3), 'co')
        plot(tmp(inds, 2), tmp(inds, 3), 'co')
    end
    pause(0.2)
    clf
end
close all

%% 2-color cross-section inspection
% Skip if already done

% Define your two-channel colors
color1 = [0.5 0 0.5] ;
color2 = [0 0.5 0.5] ;
% which timepoint
tp = 50 ;
% grab mesh filename
meshfn = sprintf(QS.fullFileBase.mesh, tp) ;
% Load the raw data via imsane Experiment instance xp
QS.setTime(tp)
QS.getCurrentData()
IV = QS.currentData.IV ;
% which page do you want to look at in cross section?
leaf = 500 ;
% unpack two channel data
if any(size(IV) > 1)
    IV1 = IV{1} ;
    IV2 = IV{2} ;
    % make an rgb image cyan/magenta
    red = IV1(leaf, :, :) * color1(1) + IV2(leaf, :, :) * color2(1) ;
    grn = IV1(leaf, :, :) * color1(2) + IV2(leaf, :, :) * color2(2) ;
    blu = IV1(leaf, :, :) * color1(3) + IV2(leaf, :, :) * color2(3) ;
    im = cat(3, red, grn, blu) ;
else
    IV = IV{1} ;
    im = squeeze(IV(leaf, :, :)) ;
    im = cat(3, im, im, im) ;
end
% Load up the mesh
mesh = read_ply_mod(meshfn) ;
% Make this number larger to sample more of the nearby mesh
width = 5 ;
% Show the cross-section
inds = find(abs(mesh.v(:, 1) - leaf) < width) ;
imshow(permute(im, [2, 1, 3]))
if any(inds)
    hold on;
    plot(mesh.v(inds, 2), mesh.v(inds, 3), 'co')
end
                
%% Make sure all meshes exist
for tp = xp.fileMeta.timePoints
    % Load the mesh
    meshfn = sprintf( QS.fullFileBase.mesh, tp ) ;    
    if ~exist(meshfn, 'file')
        error(['mesh does not exist: ' meshfn ])
    end
end

%% Inspect all meshes in 3D
% Skip if already done

% Make an output directory for the quick-and-dirty inspection
for tp = xp.fileMeta.timePoints
    % Load the mesh
    meshfn = sprintf( QS.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)
    % saveas(gcf, fullfile(outputdir, sprintf('inspect_%04d.png', tp)))
    title(['t=' num2str(tp)])
    axis equal
    view(2)
    pause(0.1)
end
                    

%% adjust all meshes by 0.5 --- now this is done in integralDetector
% Skip
% for tp = fileMeta.timePoints
%     disp(num2str(tp))
%     meshfn = fullfile(meshDir, sprintf([ofn_smoothply '%06d.ply'], tp)) ;
%     outfn = fullfile(projectDir, sprintf([ofn_smoothply '%06d.ply'], tp)) ;
%     if ~exist(outfn, 'file')
%         mesh = read_ply_mod(meshfn) ;
%         mesh.v = mesh.v + 0.5 * [1 1 1] ;
%         plywrite_with_normals(outfn, mesh.f, mesh.v, mesh.vn)
%     else
%         disp('skip')
%     end
%     
%     meshfn = fullfile(meshDir, sprintf([ofn_ply '%06d.ply'], tp)) ;
%     outfn = fullfile(projectDir, sprintf([ofn_ply '%06d.ply'], tp)) ;
%     if ~exist(outfn, 'file')
%         mesh = read_ply_mod(meshfn) ;
%         mesh.v = mesh.v + 0.5 * [1 1 1] ;
%         plywrite(outfn, mesh.f, mesh.v)
%     else
%         disp('skip')
%     end
% end
% disp('done')


%% Check that all have been created
% Skip if already done
for tp = xp.fileMeta.timePoints
    meshfn = fullfile(sprintf(QS.fullFileBase.mesh, tp)) ;
       
    % Check that smoothply exists
    if ~exist(meshfn, 'file') 
        xp.setTime(tp);
        % xp.loadTime(tp) ;
        % xp.rescaleStackToUnitAspect();

        % make a copy of the detectOptions and change the fileName
        detectOpts2 = detectOptions ;
        detectOpts2.timepoint = xp.currentTime ;
        detectOpts2.fileName = sprintf( fn, xp.currentTime );
        if tp > xp.fileMeta.timePoints(1)
            detectOpts2.init_ls_fn = 'none' ;
        end
        xp.setDetectOptions( detectOpts2 );
        xp.detectSurface();
    else
        disp('exists')
    end
end
disp('done')

clearvars lambda1 lambda2 ilastikaxisorder mlxprogram ms_scriptDir
clearvars msls_exten msls_detOpts_fn niter niter0 nu 


%% Fix ssfactor error
% % Skip -- only do if needed
% for tp = xp.fileMeta.timePoints
%     meshfn = fullfile(sprintf(QS.fullFileBase.mesh, tp)) ;
%        
%     % Check that smoothply exists
%     mesh = read_ply_mod(meshfn) ;
%     mesh.v = mesh.v * ssfactor ;
%     plywrite_with_normals(meshfn,mesh.f,mesh.v,mesh.vn)
% end
% %%
% for tp = xp.fileMeta.timePoints
%     mesh2fn = fullfile(QS.dir.mesh, sprintf('mesh_%06d.ply', tp)) ;
%     mesh2 = read_ply_mod(mesh2fn) ;
%     disp(['begin writing - ' num2str(tp)])
%     plywrite(mesh2fn, mesh2.f, mesh2.v * ssfactor) ;
%     disp('done writing')
% end
% disp('done')

%% APDV ilastik training
% Train on anterior (A), posterior (P), background (B), and 
% dorsal anterior (D) location in different iLastik channels.
% Perform on pre-stabilized H5s, NOT on stabilized H5s. 
% 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
% Train for anterior dorsal (D) only at the first time point, because
% that's the only one that's used. 
% Dorsal anterior for the gut is 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.

%% 3. align_meshes_APDV or load transformations if already done
% Skip if already done.

% Try to load the results to test if we must do calculation 
try
    [rot, trans] = QS.getRotTrans() ;
    [xyzlim_raw, xyzlim, xyzlim_um, xyzlim_buff] = getXYZLims(QS) ;
    assert(exist(QS.fileName.startendPt, 'file') ~= 0)
    redo_alignmesh = false ;
    
    % check all timepoints
    tidx = 1; 
    while ~redo_alignmesh && tidx <= length(xp.fileMeta.timePoints)
        tt = xp.fileMeta.timePoints(tidx) ;
        searchfn = sprintf(QS.fullFileBase.alignedMesh, tt);
        if ~strcmp(searchfn(end-3:end), '.ply')
            searchfn = [searchfn '.ply'] ;
        end
        if ~exist(searchfn, 'file')
            redo_alignmesh = true ;
            disp(['did not find aligned mesh for tp = ' num2str(tt)])
        end
        tidx = tidx + 1;
    end
    
catch
    redo_alignmesh = true ;
end

if redo_alignmesh || overwrite_APDVMeshAlignment || overwrite_APDVCOMs
    disp('Calculating/Overwriting APDVMeshAlignment')
    % OUTPUTS OF THIS SECTION
    % -------
    % 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 (APDV coord sys) in microns
    % rotation_APDV.txt
    %   rotation matrix to align mesh to APDV frame
    % translation_APDV.txt
    %   translation vector to align mesh to APDV frame
    % xyzlim.txt 
    %   raw bounding box in original frame (not rotated), in full res pixels
    % xyzlim_APDV.txt
    %   bounding box in rotated frame, in full resolution pixels
    % xyzlim_APDV_um.txt
    %   bounding box in rotated frame, in microns
    % apdv_coms_rs.h5
    %   Centers of mass for A,aux_paths_and_colors P, and D in microns in rotated APDV coord system
    % 
    % Notes
    % -----
    % vertices are in units of pixels (at full resolution)
    % To take mesh to rotated + translated mesh in physical units, apply:
    %         xs = mesh.vertex.z ;
    %         ys = mesh.vertex.y ;
    %         zs = mesh.vertex.x ;
    %         vtx_rs = (rot * vtx' + trans)' * resolution
    %         
    save_figs = true ;  % save images of cntrline, etc, along the way
    preview = false ;  % display intermediate results, for debugging
    dorsal_thres = 0.5 ;  % threshold for extracting Dorsal probability cloud 
    buffer = 5 ;  % extra space in meshgrid of centerline extraction, to ensure mesh contained in volume
    plot_buffer = 40;  % extra space in plots, in um
    weight = 0.1;  % for speedup of centerline extraction. Larger is less precise
    normal_step = 0.5 ;  % how far to move normally from ptmatched vtx if a/pcom is not inside mesh
    eps = 0.01 ;  % value for DT outside of mesh in centerline extraction
    meshorder = 'xyz' ;  % ordering of axes in loaded mesh wrt iLastik output
    anteriorChannel = 1;  % which channel of APD training is anterior
    posteriorChannel = 2;  % which channel of APD training is posterior 
    dorsalChannel = 4 ;  % which channel of APD training is dorsal

    clearvars opts
    optsfn = fullfile(projectDir, 'alignAPDV_Opts.mat') ;
    if exist(optsfn, 'file') && ~overwrite_alignAPDVOpts
        disp('Loading options from disk')
        load(optsfn, 'alignAPDVOpts', 'apdvOpts')
    else
        disp('No alignAPDV_Opts on disk or overwriting, defining')
        apdvOpts.smwindow = 30 ;
        apdvOpts.dorsal_thres = dorsal_thres ;
        apdvOpts.buffer = buffer ;  
        apdvOpts.plot_buffer = plot_buffer ;
        apdvOpts.anteriorChannel = anteriorChannel ;
        apdvOpts.posteriorChannel = posteriorChannel ;
        apdvOpts.dorsalChannel = dorsalChannel ;% filename pattern for the apdv training probabilities
        apdvOpts.tpref = t0_for_phi0 ;
        
        alignAPDVOpts.weight = weight ;
        alignAPDVOpts.normal_step = normal_step ;
        alignAPDVOpts.eps = eps ;
        alignAPDVOpts.meshorder = meshorder ;
        alignAPDVOpts.anteriorChannel = anteriorChannel ;
        alignAPDVOpts.posteriorChannel = posteriorChannel ;
        alignAPDVOpts.dorsalChannel = dorsalChannel ;
        alignAPDVOpts.dorsal_thres = 0.5 ;
        alignAPDVOpts.tref = t0_for_phi0 ;
        
        % alignAPDVOpts.fn = fn ;  % filename base

        % save the optionsc
        save(optsfn, 'alignAPDVOpts', 'apdvOpts')
    end

    % Add script-instance-specific options
    apdvOpts.overwrite = overwrite_APDVCOMs ;  % recompute APDV coms from training
    apdvOpts.check_slices = false ;
    apdvOpts.preview = preview ;
    apdvOpts.preview_com = false ;
    
    %% Compute APDV coordinate system
    alignAPDVOpts.overwrite = false ;
    QS.computeAPDVCoords(alignAPDVOpts) ;
    
    % Compute the APD COMs
    [acom_sm, pcom_sm] = QS.computeAPDCOMs(apdvOpts) ;
    
    % Align the meshes APDV & plot them
    alignAPDVOpts.overwrite_ims = overwrite_alignedMeshIms ;  % overwrite images even if centerlines are not overwritten
    alignAPDVOpts.overwrite = overwrite_APDVCOMs || overwrite_APDVMeshAlignment ; % recompute APDV rotation, translation
    QS.alignMeshesAPDV(alignAPDVOpts) ;
else
    % Display APDV COMS over time
    acom_sm = h5read(QS.fileName.apdv, '/acom_sm') ;
    pcom_sm = h5read(QS.fileName.apdv, '/pcom_sm') ;
    acoms = h5read(QS.fileName.apdv, '/acom') ;
    pcoms = h5read(QS.fileName.apdv, '/pcom') ;
    dcom = dlmread(QS.fileName.dcom) ;
    for tidx = 1:length(timePoints)
        tp = timePoints(tidx) ;
        % Plot the APDV points
        clf
        plot3(acom_sm(tidx, 1), acom_sm(tidx, 2), acom_sm(tidx, 3), 'ro')
        hold on;
        plot3(acoms(tidx, 1), acoms(tidx, 2), acoms(tidx, 3), 'r.')
        plot3(pcom_sm(tidx, 1), pcom_sm(tidx, 2), pcom_sm(tidx, 3), 'b^')
        plot3(pcoms(tidx, 1), pcoms(tidx, 2), pcoms(tidx, 3), 'b.')
        plot3(dcom(1, 1), dcom(1, 2), dcom(1, 3), 'cs')
        axis equal
        title(['t = ', num2str(tp)]) 
        pause(0.01)
    end
    disp('Already done')
end
disp('done')
clearvars normal_step 


%% Fix vertex normals in alignedMeshes (hack)
% for tt = QS.xp.fileMeta.timePoints ;
%     alignedmeshfn = fullfile(QS.dir.alignedMesh, ...
%         sprintf([QS.fileBase.alignedMesh '.ply'], tt)) ;
%     mm = read_ply_mod(alignedmeshfn) ;
%     xyzrs = mm.v ;
%     vn_rs = mm.vn ;
%     vn_rs(:, 1) = -vn_rs(:, 1) ;
%     vn_rs(:, 3) = -vn_rs(:, 3) ;
%     
%     % Save it
%     if overwrite || ~exist(alignedmeshfn, 'file')
%         disp('Saving the aligned mesh...')
%         disp([' --> ' alignedmeshfn])
%         plywrite_with_normals(alignedmeshfn, mm.f, xyzrs, vn_rs)
%     end
% end

%% MAKE MASKED DATA FOR PRETTY VIDEO ======================================
% Skip if already done
% Generate masks for isolating intensity data of the gut alone, for making
% pretty videos without fiducial markers and other spurious fluorescent
% bits. 
% Save output h5s trained on stabilized h5s from iLastik as 
%   -->   <QS.fileBase.name>_Probabilities_mask3d.h5
% and, optionally, a probability cloud to exclude 
% (applied post extraction of previous mask)
%   -->   <QS.fileBase.name>_Probabilities_maskDorsal.h5
options = struct() ;
options.maskMask = false ;
QS.generateMaskedData(options)

%% MAKE ORIENTED MASKED DATA FOR PRETTY VIDEO =============================
% Skip if already done
QS.alignMaskedDataAPDV()

%% PLOT ALL TEXTURED MESHES IN 3D =========================================
% Skip if already done
overwrite = false ;

% Get limits and create output dir
% Establish texture patch options
metafn = fullfile(QS.dir.texturePatchIm, 'metadat.mat') ;
if ~exist(metafn, 'file') || overwrite_TextureMeshOpts
    [~,~,~,xyzbuff] = QS.getXYZLims() ;
    xyzbuff(:, 1) = xyzbuff(:, 1) - 20 ; 
    xyzbuff(:, 2) = xyzbuff(:, 2) + 20 ; 
    % Define & Save metadata
    metadat.xyzlim = xyzbuff ;                  % xyzlimits
    metadat.reorient_faces = false ;            % if some normals are inverted
    metadat.normal_shift = QS.normalShift ;             % normal push, in pixels, along normals defined in data XYZ space
    metadat.texture_axis_order = QS.data.axisOrder ;    % texture space sampling
        
    % Psize is the linear dimension of the grid drawn on each triangular face
    Options.PSize = 5 ;
    Options.EdgeColor = 'none';
    QS.getRotTrans() ;
    Options.Rotation = QS.APDV.rot ;
    Options.Translation = QS.APDV.trans ;
    Options.Dilation = QS.APDV.resolution ;
    % outward, inward --> reasonable values are ~[3,3]
    Options.numLayers = [1, -5];  % at layerSpacing 2, 2 marches ~0.5 um 
    Options.layerSpacing = 2 ;
    
    % Save it
    save(metafn, 'metadat', 'Options')
else
    load(metafn, 'metadat', 'Options')
end

% Use first timepoint's intensity limits throughout
QS.setDataLimits(QS.xp.fileMeta.timePoints(101), 1.0, 99.9)
% QS.setDataLimits(QS.xp.fileMeta.timePoints(101), 1.0, 99.95)

%% Plot on surface for all TP 
options = metadat ;
options.overwrite = false ;
options.plot_dorsal = false ;
options.plot_ventral = false ;
options.plot_right = false ;
options.plot_left = false ;
options.plot_perspective = true ;
options.blackFigure = false ;
options.timePoints = 0:169 ;
QS.plotSeriesOnSurfaceTexturePatch(options, Options)
clearvars Options

%% Morphsnakes demo evolution visualizeMeshEvolution orthoviews data planes
options = struct() ;
options.plot_evolution = true ;
options.plot_growth = false ;
options.plot_texture = false ;
options.plot_meshOnly = false ;
options.brighten = 100; 
options.growth_t0 = 85 ;
options.y0 = 20 ;
options.overwrite = true ;
% options.viewAngles = [-0.75, 1, 0.7] ;
QS.visualizeMeshEvolution(options)

%% Highligh junction for BWF fund application
amesh = QS.getCurrentAlignedMesh() ;
xwidth = 32 ; % cm
ywidth = 20 ; % cm
% scaled mesh vertices
amesh.v = laplacian_smooth(amesh.v, amesh.f, 'cotan', [], 0.001,'implicit',amesh.v,1000) ;
xa = amesh.v(:, 1)  ;
ya = amesh.v(:, 2)  ;
za = amesh.v(:, 3)  ;
sky    = colors(6, :) ;
yellow = colors(3, :) ;
lscolors = sky ;


close all            
fig = figure('units', 'centimeters', ...
    'outerposition', [0 0 xwidth xwidth], 'visible', 'off') ;
h = trimesh(amesh.f, xa, ya, za, 'EdgeColor', 'none', ...
    'facecolor', lscolors) ;
axis equal
hold on;
axis off

lighting gouraud    % preferred method for lighting curved surfaces
material dull    % set material to be dull, no specular highlights

% Lighting and view 
% view(-20, 20)
view(viewAngles) 
lgt = camlight('headlight') ;
set(gcf,'color','w');

%% Symmetry figure
% meshL = read_ply_mod(sprintf(...
%     [QS.fullFileBase.alignedMesh(1:end-4) '_dense_extrasmooth.ply'], 130)) ;
% trisurf(triangulation(meshL.f, meshL.v), meshL.vn(:, 2), 'edgecolor', 'none')
% meshR = meshL ;
% meshR.v(:, 2) = - meshL.v(:, 2) + 200 ;
% clf
% trisurf(triangulation(meshL.f, meshL.v), -meshL.vn(:, 2), 'edgecolor', 'none')
% hold on;
% trisurf(triangulation(meshR.f, meshR.v), meshL.vn(:, 2), 'edgecolor', 'none')
% view(2)
% axis equal
% xlims = get(gca, 'xlim') ;
% xlim([xlims(1)-10, xlims(2) + 10])
% grid off
% set(gca, 'color', 'k', 'xcol', 'w', 'ycol', 'w', 'zcol', 'w')
% set(gcf, 'InvertHardCopy', 'off');
% set(gcf, 'Color', 'k')
% set(gcf, 'color', 'k')
% export_fig(fullfile(meshDir, 'symmetry2.png'), '-nocrop', '-r400')

%% EXTRACT CENTERLINES
% Skip if already done
% Note: these just need to be 'reasonable' centerlines for topological
% checks on the orbifold cuts.
exponent = 1.0 ;
res = 2.0 ; 
cntrlineOpts.overwrite = overwrite_centerlines ;     % overwrite previous results
cntrlineOpts.overwrite_ims = overwrite_centerlineIms ;     % overwrite previous results
cntrlineOpts.weight = 0.1;               % for speedup of centerline extraction. Larger is less precise
cntrlineOpts.exponent = exponent ;       % how heavily to scale distance transform for speed through voxel
cntrlineOpts.res = res ;                 % resolution of distance tranform grid in which to compute centerlines
cntrlineOpts.preview = false ;           % preview intermediate results
cntrlineOpts.reorient_faces = false ;    % not needed for our well-constructed meshes
cntrlineOpts.dilation = 0 ;              % how many voxels to dilate the segmentation inside/outside before path computation

% Note: this takes about 400s per timepoint for res=2.0
%
QS.extractCenterlineSeries(cntrlineOpts)
disp('done with centerlines')


%% Fix flip in Y for centerlines
% aux_fix_flip

%% Identify anomalies in centerline data
% Skip if already done
idOptions.ssr_thres = 15 ;  % distance of sum squared residuals in um as threshold
idOptions.overwrite = overwrite_idAnomClines ;
QS.generateCleanCntrlines(idOptions) ;


%% Cylinder cut mesh
% Skip if already done
overwrite_endcapOpts = true ;
if overwrite_endcapOpts || ~exist(QS.fileName.endcapOptions, 'file')
    endcapOpts = struct( 'adist_thres', 20, ...  % 20, distance threshold for cutting off anterior in pix
                'pdist_thres', 20, ...  % 15-20, distance threshold for cutting off posterior in pix
                'tref', 110) ;  % reference timepoint at which time dorsal-most endcap vertices are defined
    QS.setEndcapOptions(endcapOpts) ;
    % Save the options to disk
    QS.saveEndcapOptions() ;
else
    % load endcapOpts
    QS.loadEndcapOptions() ;
    endcapOpts = QS.endcapOptions ;
end

clearvars methodOpts
methodOpts.overwrite = overwrite_endcapOpts ;  % recompute sliced endcaps
methodOpts.save_figs = true ;   % save images of cutMeshes along the way
methodOpts.preview = false  ;     % display intermediate results
QS.sliceMeshEndcaps(endcapOpts, methodOpts) ;

%% Clean Cylinder Meshes
% May skip if already done
cleanCylOptions.overwrite = overwrite_cleanCylMesh ;
cleanCylOptions.save_ims = true ;
QS.cleanCylMeshes(cleanCylOptions)
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ORBIFOLD -> begin populating Qs.dir.mesh/gridCoords_nUXXXX_nVXXXX/ 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Parameters
% Overwriting options
overwrite_pullbacks = false ;
overwrite_cutMesh = false ;
overwrite_spcutMesh = false ;
% Plotting params
washout2d = 0.5 ;

%% Iterate Through Time Points to Create Pullbacks ========================
% Skip if already done
% outcutfn = fullfile(cutFolder, 'cutPaths_%06d.txt') ;
overwrite_spcutMesh = false ;
overwrite_cutMesh = false ;
tp2do = [t0_for_phi0:max(xp.fileMeta.timePoints), ...
    fliplr(min(xp.fileMeta.timePoints):(t0_for_phi0-1))] ;

for tt = tp2do
    disp(['NOW PROCESSING TIME POINT ', num2str(tt)]);
    tidx = xp.tIdx(tt);
    
    % Load the data for the current time point ------------------------
    QS.setTime(tt) ;
    
    %----------------------------------------------------------------------
    % Create the Cut Mesh
    %----------------------------------------------------------------------
    cutMeshfn = sprintf(QS.fullFileBase.cutMesh, tt) ;
    cutPathfn = sprintf(QS.fullFileBase.cutPath, tt) ;
    if ~exist(cutMeshfn, 'file') || ~exist(cutPathfn, 'file') || overwrite_cutMesh
        if exist(cutMeshfn, 'file')
            disp('Overwriting cutMesh...') ;
        else
            disp('cutMesh not found on disk. Generating cutMesh... ');
        end
        options = struct() ;
        options.preview = false ;
        options.t0 = t0_for_phi0 ;
        QS.generateCurrentCutMesh(options)
        disp('Saving cutP image')
        % Plot the cutPath (cutP) in 3D
        QS.plotCutPath(QS.currentMesh.cutMesh, QS.currentMesh.cutPath)
        compute_pullback = true ;
    else
        fprintf('Loading Cut Mesh from disk... ');
        QS.loadCurrentCutMesh()
        compute_pullback = ~isempty(QS.currentMesh.cutPath) ;
    end
    
    spcutMeshOptions = struct() ;
    spcutMeshOptions.t0_for_phi0 = t0_for_phi0 ;
    spcutMeshOptions.overwrite = overwrite_spcutMesh ;
    spcutMeshOptions.save_phi0patch = false ;
    spcutMeshOptions.iterative_phi0 = true ;
    spcutMeshOptions.smoothingMethod = 'none' ;
    QS.plotting.preview = false ;
    QS.generateCurrentSPCutMesh([], spcutMeshOptions) ;
    
    % Compute the pullback if the cutMesh is ok
    if compute_pullback || ~exist(sprintf(QS.fullFileBase.im_sp, tt), 'file')
        pbOptions.overwrite = false ;
        pbOptions.generate_uv = false ;
        pbOptions.generate_sphi = true ;
        pbOptions.generate_uphi = false ;
        pbOptions.generate_relaxed = false ;
        QS.generateCurrentPullbacks([], [], [], pbOptions) ;
    else
        disp('Skipping computation of pullback')
    end
    clear Options IV
        
    % Save SMArr2D (vertex positions in the 2D pullback) -----------------
    % disp(['Saving meshStack to disk: ' mstckfn])
    % save(mstckfn, 'meshStack') ;
    % 
    % %% Save SMArr2D (vertex positions in the 2D pullback) -----------------
    % disp(['Saving spmeshStack to disk: ' spmstckfn])
    % if generate_sphi_coord
    %     save(spmstckfn, 'spmeshStack') ;
    % end
end
clearvars t cutP spcutMesh spvn3d ss pp uv tileCount slin plin plotfn
clearvars IVloaded IV uphi sphi
disp('Done with generating spcutMeshes and cutMeshes')

%% FLIP Y for spcutMesh mcline
% DO ONLY ONCE IF NEEDED (no longer needed for datasets processed after
% 12-03-2020
% for tidx = 1:length(QS.xp.fileMeta.timePoints)
%     tp = QS.xp.fileMeta.timePoints(tidx) ;
%     load(sprintf(QS.fullFileBase.spcutMesh, tp), 'spcutMesh') ;
%     spcutMesh.mcline(:, 2) = - spcutMesh.mcline(:, 2) ;
%     spcutMesh.avgpts(:, 2) = - spcutMesh.avgpts(:, 2) ;
%     save(sprintf(QS.fullFileBase.spcutMesh, tp), 'spcutMesh') ;
% end


%% Preview results ========================================================
check = false ;
if check
    aux_preview_results
end

%% FIND THE FOLDS SEPARATING COMPARTMENTS =================================
% Skip if already done
options = struct() ;
options.overwrite = true ;
options.preview = true ;
options.first_tp_allowed = [-1, 120, -1] ;  % enforce that no folds before this tp
options.guess123 = [.29 .55 .80] ;
options.max_wander = 5 ;
QS.identifyFolds(options)
disp('done')

%% Circularlity of constrictions -- todo?
% options = struct() ;
% QS.measureFoldRadiiVariance(options)

%% COMPUTE MESH SURFACE AREA AND VOLUME ===================================
% Skip if already done
% Note: doing this after fold identification so that t0 is defined for
% plotting purposes
options = struct() ;
options.overwrite = true ;
QS.measureSurfaceAreaVolume(options)
disp('done')

% RECOMPUTE WRITHE OF MEANCURVE CENTERLINES ==============================
% Skip if already done
options = struct() ;
options.overwrite = true ;
QS.measureWrithe(options)
disp('done')

%% Plot fancy "cross-section" view of centerlines
options = struct() ;
QS.plotClineXSections(options)

%% Compute surface area and volume for each compartment ===================
% Skip if already done
options = struct() ;
options.overwrite = true ;
QS.measureLobeDynamics(options) ;

%% plot length, area, and volume for each lobe ============================
% Skip if already done
options = struct() ;
options.overwrite = true ;
QS.plotLobes(options)

%% Plot motion of avgpts & DVhoops at folds in yz plane over time =========
% Skip if already done
opts = struct() ;
opts.overwrite = true ;
opts.timePoints = [123:30:256] ;
QS.plotConstrictionDynamics(opts) ;
disp('done')

%% Tubular figure dec rotation smoothed

outdir = fullfile([QS.dir.piv.avgDEC.rot3d '_smoothed']) ;
mkdir(outdir)
[~,~,~,xyzlims] = QS.getXYZLims() ;
for tt = QS.t0set():10:max(QS.xp.fileMeta.timePoints)-1
    
    % smooth dec in time
    first = true ;
    twidth = 15 ;
    tIdxs2Add = -twidth:twidth ;
    tpulse = tripulseFunction(twidth) ;
    for qq = 1:length(tIdxs2Add)
        tadd = tIdxs2Add(qq) ;
        QS.setTime(tt+tadd)
        dec = QS.getCurrentDEC() ;
        rotv_qq = dec.rots.rotv(1:end-QS.nU) ;
        mesh = QS.getCurrentSPCutMeshSmRSC() ;

        % Silence endcaps
        rotv_qq(1:QS.nU:end) = 0 ;
        rotv_qq(QS.nU:QS.nU:end) = 0 ;
        
        if first
            rotv = tpulse(qq) * rotv_qq ;
            first = false ;
        else
            rotv = rotv + tpulse(qq) * rotv_qq ;
        end
    end
    
    % get current mesh
    QS.setTime(tt)
    mesh = QS.getCurrentSPCutMeshSmRSC() ;

    % Smooth result in space
    rotv = laplacian_smooth(mesh.v, mesh.f, 'cotan',[],0.05,'implicit',rotv,1000) ;
    opts = struct('nmodesY', 7, 'widthX', 3) ;
    rotv = reshape(rotv, [QS.nU, QS.nV-1]) ;
    rotv = modeFilterQuasi1D(rotv, opts) ;
    
    % Plot it
    subtightplot(2, 1, 1)
    trisurf(triangulation(mesh.f, mesh.v), rotv(:), ...
        'edgecolor', 'none') ;
    % colormap(brewermap(256, '*RdBu'))
    colormap(bwr)
    cmax = 0.1 ;
    caxis(cmax*[-1,1]) ;
    axis equal
    grid off
    set(gcf, 'color', 'w')
    view(0,0)
    xlim(xyzlims(1, :))
    ylim(xyzlims(2, :))
    zlim(xyzlims(3, :))
    axis off
    
    % flat pullback with aspect 1:1
    subtightplot(2, 1, 2)
    rotvM = rotv ;
    rotvM(:,QS.nV) = rotv(:, 1) ;
    imagesc(linspace(0, 1.31953235, 100), linspace(0,1,100), rotvM')
    axis equal
    axis tight
    axis off
    caxis(cmax*[-1,1]) ;
    outfn = fullfile(outdir, sprintf('rot_sm_%06d.png',tt)) ;
    saveas(gcf, outfn) 
end

%% SMOOTH MEAN CENTERLINE RADIUS ==========================================
% Skip if already done
% todo: rework this section so that it takes place after smoothing meshes
% overwrite = false ;
% aux_smooth_avgptcline_radius_before_mesh_smoothing(overwrite, ...
%     QS.xp.fileMeta.timePoints, QS.fullFileBase.spcutMesh, ...
%     QS.fullFileBase.clineDVhoop, radiusImDir, ...
%     QS.APDV.rot, QS.APDV.trans, QS.APDV.resolution, ...
%     QS.plotting.xyzlim_um, QS.nU, QS.nV)

%% Smooth the sphi grid meshes in time ====================================
% Skip if already done
options = struct() ;
options.overwrite = false ;
options.width = 4 ;  % before 2020-12-09, this was 6.
QS.smoothDynamicSPhiMeshes(options) ;

%% Fix one field of spcutMeshSm / spcutMeshSmRS / spcutMeshSmRSC -- 2020-12-21
% for tidx = 1:length(QS.xp.fileMeta.timePoints)
%     tt = QS.xp.fileMeta.timePoints(tidx) ;
%     
%     spcutMeshBase = QS.fullFileBase.spcutMesh ;
%     spcutMeshSmBase = QS.fullFileBase.spcutMeshSm ;
%     spcutMeshSmRSBase = QS.fullFileBase.spcutMeshSmRS ;
%     spcutMeshSmRSCBase = QS.fullFileBase.spcutMeshSmRSC ;
% 
%     % Load the files
%     load(sprintf(spcutMeshBase, tt), 'spcutMesh') ;
%     load(sprintf(spcutMeshSmBase, tt), 'spcutMeshSm') ;
%     load(sprintf(spcutMeshSmRSBase, tt), 'spcutMeshSmRS') ;
%     load(sprintf(spcutMeshSmRSCBase, tt), 'spcutMeshSmRSC') ;
%     
%     spcutMeshSm.u = spcutMesh.sphi ;
%     spcutMeshSmRS.u = spcutMesh.sphi ;
%     spcutMeshSmRSC.u = spcutMesh.sphi(1:nU*(nV-1), :) ;
%     assert(size(spcutMeshSmRSC.u, 1) == size(spcutMeshSmRSC.v, 1))
%     
%     % Resave the files
%     disp(['Saving ' sprintf(spcutMeshSmBase, tt)])
%     save(sprintf(spcutMeshSmBase, tt), 'spcutMeshSm') ;
%     save(sprintf(spcutMeshSmRSBase, tt), 'spcutMeshSmRS') ;
%     save(sprintf(spcutMeshSmRSCBase, tt), 'spcutMeshSmRSC') ;
%             
% end

%% Plot the time-smoothed meshes
% Skip if already done
options.overwrite = false ;
QS.plotSPCutMeshSmRS(options) ;

% Images for publication/presentation on method & coordinate system
% Skip if already done
% Create coordinate system charts visualization using smoothed meshes

QS.coordSystemDemo()
QS.setTime(QS.t0 + 90)
mesh = QS.loadCurrentSPCutMeshSmRS() ;
fig = figure('units', 'centimeters', 'position', [0, 0, 13, 13]) ;
uu = mesh.u(:, 1) ;
vv = mesh.u(:, 2) ;
xx = (mesh.v(:, 1)) ;
yy = (mesh.v(:, 2)) ;
zz = (mesh.v(:, 3)) ;
colorsV = viridis(mesh.nV) ;
colorsU = viridis(mesh.nU) ;

% LONGITUDE 3D
subplot(2, 3, [1,2])
hold off
for qq = 1:mesh.nV  
    plot3(xx(qq:nU:end), yy(qq:nU:end), zz(qq:nU:end), '-', ...
        'color', colorsV(qq, :));
    hold on ;
end
axis equal
xlabel('ap position [$\mu$m]', 'interpreter', 'latex')
ylabel('lateral position [$\mu$m]', 'interpreter', 'latex')
zlabel('dv position [$\mu$m]', 'interpreter', 'latex')

% AZIMUTH 3D
subplot(2, 3, [4,5])
hold off
for qq = 1:2:mesh.nU
    inds = (qq-1)*nU+1:qq*nU ;
    plot3(xx(inds), yy(inds), zz(inds), '-', ...
        'color', colorsV(qq, :));
    hold on ;
end
axis equal
xlabel('ap position [$\mu$m]', 'interpreter', 'latex')
ylabel('lateral position [$\mu$m]', 'interpreter', 'latex')
zlabel('dv position [$\mu$m]', 'interpreter', 'latex')

% LONGITUDE
subplot(2, 3, 3)
hold off
for qq = 1:mesh.nV  
    plot(uu(qq:nU:end), vv(qq:nU:end), '-', ...
        'color', colorsV(qq, :));
    hold on ;
end
axis square
xlabel('$s$ [$\mu$m]', 'interpreter', 'latex')
ylabel('$\phi$ [1/$2\pi$]', 'interpreter', 'latex')

% AZIMUTH
subplot(2, 3, 6)
hold off
for qq = 1:mesh.nV  
    inds = (qq-1)*nU+1:qq*nU ;
    plot(uu(inds), vv(inds), '-', ...
        'color', colorsV(qq, :));
    hold on ;
end
axis square
xlabel('$s$ [$\mu$m]', 'interpreter', 'latex')
ylabel('$\phi$ [1/$2\pi$]', 'interpreter', 'latex')

set(gcf, 'color', 'w')
export_fig(fullfile(QS.dir.uvCoord, 'coordSystemDemo.png'), '-nocrop','-r600')

%% COMPUTE MEAN AND GAUSSIAN CURVATURES OF SMOOTHED MESHES
% Skip if already done
options = struct() ;
options.overwrite = false ;
QS.measureCurvatures(options)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Redo Pullbacks with time-smoothed meshes ===============================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Skip if already done
disp('Create pullback using S,Phi coords with time-averaged Meshes')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for tt = QS.t0set() % , QS.xp.fileMeta.timePoints]
    disp(['NOW PROCESSING TIME POINT ', num2str(tt)]);
    tidx = QS.xp.tIdx(tt);
    
    % Load the data for the current time point ------------------------
    QS.setTime(tt) ;
    
    % OPTION 1: Keep constant luminosity throughout, modify default 
    % intensity limits.
    % if tidx == 1        
    %     % Use first timepoint's intensity limits throughout
    %     QS.setDataLimits(QS.xp.fileMeta.timePoints(1), 1.0, 99.995)
    % end
    % QS.data.adjustlow 
    % QS.data.adjusthigh
    
    % OPTION 2 : Adjust intensity to scale from timepoint to timepoint
    adjustlow = 1.00 ;         % floor for intensity adjustment
    adjusthigh = 99.9 ;        % ceil for intensity adjustment (clip)
    QS.data.adjustlow = adjustlow ;
    QS.data.adjusthigh = adjusthigh ;
    
    % Establish custom Options for MIP
    pbOptions = struct() ;
    pbOptions.overwrite = true ;
    pbOptions.numLayers = [7 15] ; % previously [7, 7] ;  % previously [5,5]
    pbOptions.layerSpacing = 0.75 ;
    pbOptions.generate_rsm = false ;
    pbOptions.generate_spsm = false ;
    pbOptions.generate_sphi = false ;
    pbOptions.generate_uvprime = false ;
    pbOptions.generate_ruvprime = false ;
    pbOptions.generate_ricci = true ;
    pbOptions.imSize = 2000 ;
    pbOptions.normal_shift = -4 ;
    QS.data.adjustlow = adjustlow ;
    QS.data.adjusthigh = adjusthigh ;
    pbOptions.save_as_stack = true ;
    QS.generateCurrentPullbacks([], [], [], pbOptions) ;
end


%% TILE/EXTEND SMOOTHED IMAGES IN Y AND RESAVE =======================================
% Skip if already done
options = struct() ;
options.overwrite = false ;
% options.coordsys = 'spsm' ;
% QS.doubleCoverPullbackImages(options)
options.coordsys = 'rsm' ;
QS.doubleCoverPullbackImages(options)
disp('done')

%% Compute whether pullback is isothermal -> metric images
options = struct() ;
options.overwrite = false ;
options.makeRawMetricComponentFigures = false ;
options.lambda_mesh = 0.002 ;
options.coordSys = 'ricci' ; % 'spsm_rs'
QS.plotMetric(options) ;
% !!! todo: continue here for diagnostic
% Note: demo_FundFormNematic

%% Compare 2nd fund form director to radon of image
options = struct() ;
options.method = 'pullback' ; 
options.coordSys = 'spsmre' ;
QS.measurePolarity(options) ;
% inside measurePolarity: comparePolarityHopfDifferential()

%% Add radii to spcutMeshSm's in post
% 
% for tidx = 1:length(QS.xp.fileMeta.timePoints)
%     tp = QS.xp.fileMeta.timePoints(tidx) ;
%     % Load this timepoint spcutMeshSm
%     QS.setTime(tp)
%     spcutMeshSm = QS.getCurrentSPCutMeshSm() ;
%     spcutMeshSmRS = QS.getCurrentSPCutMeshSmRS() ;
%     spcutMeshSmRSC = QS.getCurrentSPCutMeshSmRSC() ;
%     
%     % Make avgpts in pixel space (not RS)
%     fprintf('Resampling uvgrid3d curves in pix...\n')
%     nU = spcutMeshSm.nU ;
%     nV = spcutMeshSm.nV ;
%     curves3d_pix = reshape(spcutMeshSm.v, [nU, nV, 3]) ;
%     c3d_dsv_pix = zeros(size(curves3d_pix)) ;  % in units of pix
%     avgpts_pix = zeros(nU, 3) ;
%     radius_pix = zeros(nU, nV) ;
%     for i=1:nU
%         % Note: no need to add the first point to the curve
%         % since the endpoints already match exactly in 3d and
%         % curvspace gives a curve with points on either
%         % endpoint (corresponding to the same 3d location).
%         c3d_dsv_pix(i, :, :) = resampleCurvReplaceNaNs(squeeze(curves3d_pix(i, :, :)), nV, true) ;
%         if vecnorm(squeeze(c3d_dsv_pix(i, 1, :)) - squeeze(c3d_dsv_pix(i, end, :))) > 1e-7
%             error('endpoints do not join! Exiting')
%         end
%         % Drop the final endpoint in the mean pt determination
%         avgpts_pix(i, :) = mean(squeeze(c3d_dsv_pix(i, 1:end-1, :)), 1) ; 
%         radius_pix(i, :) = vecnorm(squeeze(curves3d_pix(i, :, :)) - avgpts_pix(i, :), 2, 2) ;
%     end
% 
%     % uvpcutMesh.raw.avgpts_pix = avgpts_pix ;
%     % uvpcutMesh.raw.radius_pix = radius_pix ;
%     spcutMeshSm.avgpts_um = QS.xyz2APDV(avgpts_pix) ;
%     spcutMeshSm.radius_um = radius_pix * QS.APDV.resolution ;
%     
%     % Add to RS and RSC versions
%     spcutMeshSmRS.avgpts_um = spcutMeshSm.avgpts_um ;
%     spcutMeshSmRS.radius_um = spcutMeshSm.radius_um ;
%     % RSC
%     spcutMeshSmRSC.avgpts_um = spcutMeshSm.avgpts_um ;
%     spcutMeshSmRSC.radius_um = spcutMeshSm.radius_um(:, 1:end-1) ;
%     
%     % Save spcutMeshSm / RS / RSC
%     save(sprintf(QS.fullFileBase.spcutMeshSm, QS.currentTime), ...
%         'spcutMeshSm')
%     save(sprintf(QS.fullFileBase.spcutMeshSmRS, QS.currentTime), ...
%         'spcutMeshSmRS')
%     save(sprintf(QS.fullFileBase.spcutMeshSmRSC, QS.currentTime), ...
%         'spcutMeshSmRSC')
%     
%     if mod(tidx, 10) == 0 
%         clf; set(gcf, 'visible', 'on')
%         trisurf(triangulation(spcutMeshSmRSC.f, spcutMeshSmRSC.v), 'edgecolor', 'none')
%         axis equal 
%         title(num2str(tidx))
%         pause(1)
%         cla
%         trisurf(triangulation(spcutMeshSmRS.f, spcutMeshSmRS.v), 'edgecolor', 'k')
%         axis equal
%         title([num2str(tidx) ' open'])
%         pause(0.0001)
%     end
% end

%% Cell Segmentation

% Automatic segmentation gives a rough guess
error('only do this once')
options = struct() ;
options.overwrite = true ;
options.overwriteImages = false;
options.timePoints = [96:10:206] ;
QS.generateCellSegmentation2D(options) 

%% PERFORM Manual corrections (in GIMP, for ex) and re-process results
options = struct() ;
options.overwrite = false ;
options.overwriteImages = false;
options.timePoints = [123,126] ; %[96:10:206] ;
QS.processCorrectedCellSegmentation2D(options) 

%% Project into 3D
options = struct() ;
options.timePoints = [96:10:206] ;
options.overwrite = false ;
options.overwriteImages = false ;
options.correctedSegmentation = true ;
QS.generateCellSegmentation3D(options) 

%% Report how many cells segmented
coordSys = 'spsme' ;
for tii = 1:length(options.timePoints) 
    tp = options.timePoints(tii) ;
    segfn = sprintf(...
        QS.fullFileBase.segmentation2dCorrected, ...
        coordSys, tp) ;
    %% Convert to simpler format
    tmp = load(segfn, 'seg2d') ;
    nseg = size(tmp.seg2d.cdat.centroid, 1) ;
    disp(['n= ' num2str(nseg)])
end

%% Summarize statistics of segmentation results (in 3d)
options = struct() ;
options.correctedSegmentation = true ;
options.timePoints = [96:10:206] ;
options.overwrite = false  ;
QS.plotSegmentationStatisticsLobes(options)

%% Compare segmentation to an advected segmentation
options = struct() ;
% options.timePoints = [93:15:263] ;
options.timePoints = [96:10:206] ;
options.overwrite = false ;
options.overwritePathlines = false ;
options.overwriteImages = false ;

% Optional: load pathlines before call
cellVertexPathlineFn = fullfile(QS.dir.segmentation, 'pathlines', ...
    sprintf('cellVertexPathlines_%06dt0.mat', QS.t0set())) ;
% load(cellVertexPathlineFn, 'segVertexPathlines2D', ...
%         'segVertexPathlines3D', 'cellIDs')
% options.segmentationPathlines = struct() ;
% options.segmentationPathlines.segVertexPathlines2D = segVertexPathlines2D ;
% options.segmentationPathlines.segVertexPathlines3D = segVertexPathlines3D ;
% options.segmentationPathlines.cellIDs = cellIDs ;
QS.generateCellSegmentationPathlines3D(options)
% QS.estimateIntercalationRate(options)


%% Visualize kymograph of trace of strain
% TODO HERE


%% Visualize demoTracks in 2D ARAP submesh parameterization patch
tracks2demo = dir(fullfile(QS.dir.tracking, 'demoTracks', 'demoTracks*0.mat')) ;
options = struct() ;
options.tracks2demo = tracks2demo ;
options.overwrite = true ;
options.scaleByMetric = false ;
options.scaleByMetricComponents = true ;
QS.visualizeDemoTracks(options)

%% Visualize tracking in 3D
options = struct() ;
options.overwrite = false ;
options.method = 'segmentation'; 
options.subdir = 'labeled_groundTruth' ;
options.timePoints = [96:2:206] ;
QS.visualizeTracking3D(options)

%% Dump labeled_groundTruth into manualTracks
% load each
dmyk = 1 ;
first = true ;
tracks = cell(178, 1) ;
for qq = QS.xp.fileMeta.timePoints
    % Load this track timepoint
    fn = fullfile(QS.dir.tracking, 'labeled_groundTruth', ...
        sprintf('seg2d_%06d.mat', qq)) ;
    
    % Check against image
    segIm = fullfile(QS.dir.tracking, 'labeled_groundTruth', ...
        sprintf('tracks_label_%06d.mat', qq)) ;
    
    if exist(fn, 'file')
        disp(['loading ' num2str(qq)])
        seg2d = load(fn) ;
        for cellID = 1:length(seg2d.seg2d.cdat.centroid)
            if first
                tracks{cellID} = ...
                    nan(length(QS.xp.fileMeta.timePoints), 3);
            end
            tracks{cellID}(dmyk, 1:2) = seg2d.seg2d.cdat.centroid(cellID, :) - [0, 500] ;
            tracks{cellID}(dmyk, 3) = 0 ;
        end
        first = false ;
    else
        disp(['skipping ' num2str(qq)])
    end
    dmyk = dmyk + 1 ;
end
tracks{176} = tracksAppend{2} ;
tracks{177} = tracksAppend{3} ;
save([options.trackOutfn '_backup'], 'tracks')

%% Visualize Segmentation in a true-scale patch
options = struct() ;
options.demoPatchName = 'demoPatch001' ;
options.timePoints = [96, 126, 156, 166, 186, 206] ;
options.scaleByMetricComponents = true ;
options.overwrite = false ;
QS.visualizeSegmentationPatch(options)

%% Manual track and compute geodesics between tracked cells
options = struct() ;
options.trackOutfn = fullfile(QS.dir.tracking,...
    'manualTracking', 'manualTracks.mat') ;
options.tracks2Add = [125] ;
QS.manualTrackingAdd(options) 

%%
options = struct() ;
options.overwrite = true ;
options.coordSys = 'spsm' ;
options.method = 'nuclei'; 
options.subdir = 'manualTracking' ;
options.trackOutfn = fullfile(QS.dir.tracking, ...
    options.subdir, ...
    'manualTracks.mat') ;
options.drawPairRectangle = false ;
options.selectPairs = 2 ;
options.t0forPairs = 123 ;
options.overlayStyle = 'mask' ;
options.nPairs = 2 ;
options.overwriteROI = false ;
% options.pairIDs = [121, 181, NaN, NaN] ;
options.pairIDs = [101, 85; ...
    189,125] ;
options.pairColors = {[255,158,2]/255,...
    [0,140,255]/255.}
% options.pairIDs = [] ;
options.timePoints = [123,126,166,176,206, 124:2:216] ;
options.viewAngles = [-20, 20] ;
options.normal_shift = 2 ;
QS.visualizeTracking3D(options)



%% Measure Cell density
% Skip if already done
options = struct() ;
options.overwrite = false ;
options.preview = false ;
QS.measureCellDensity('nuclei', options)
QS.plotCellDensity(options) ;

%% Cell density kymograph
% Skip if already done
options = struct() ;
options.overwrite = false ;
options.timePoints = 0:85 ;
QS.plotCellDensityKymograph(options)

%% CREATE PULLBACK STACKS =================================================
% Skip if already done
disp('Create pullback stack using S,Phi coords with time-averaged Meshes');
% Load options
overwrite = true ;
optionfn = fullfile(QS.dir.im_r_sme_stack, 'spcutMeshSmStackOptions.mat') ;
if ~exist(optionfn, 'file') || overwrite
    spcutMeshSmStackOptions.layer_spacing = 0.5 / QS.APDV.resolution ; % pixel resolution roughly matches xy
    spcutMeshSmStackOptions.n_outward = 20 ;
    spcutMeshSmStackOptions.n_inward = 40 ;
    spcutMeshSmStackOptions.smoothIter = 0 ;
    spcutMeshSmStackOptions.preSmoothIter = 35 ;
    spcutMeshSmStackOptions.imSize = 500 ;
    % Save options
    save(optionfn, 'spcutMeshSmStackOptions')
else
    load(optionfn, 'smSPCutMeshStackOptions')
end
spcutMeshSmStackOptions.overwrite = overwrite ;
QS.generateSPCutMeshSmStack(spcutMeshSmStackOptions)

%% TRAIN IN ILASTIK ON MIDGUT TISSUE TO GET THICKNESS
% Read thickness training output, train with first channel fg, second bg,
% third optional folded-over tissue. Name output <filename>_thickness.h5

%% MEASURE THICKNESS ======================================================
% Measure thickness from images of imFolder_sp_r_sme_stack (extendedLUT_smoothed)
thicknessOptions = struct() ;
QS.measureThickness(thicknessOptions)

%% DUMP OR LOAD HERE [break]
clearvars fold_ofn fig1exist fig2exist id idx mcline nsmM prevcline tmp
clearvars fig1exist fig1fn fig2exist fp1fn fp2fn fp3fn fp4fn IV 
clearvars n_inward n_outward TV2D TV3D 
clearvars layer_spacing 
dumpfn = fullfile(meshDir, 'orbifold_dump_before_piv.mat') ;
save(dumpfn)
load(dumpfn)
clearvars dumpfn


%% PERFORM PIV ON PULLBACK MIPS ===========================================
% % Compute PIV in PIVLab
% % ---------------------
% % Open PIVLab
% % Select all frames in meshDir/PullbackImages_010step_sphi/smoothed_extended/
% % Select Sequencing style 1-2, 2-3, ... 
% % Image Preprocessing (used to select all, but now:)
% %  --> Enable CLAHE with 20 pix
% %  --> DO NOT Enable highpass with 15 pix
% %  --> DO NOT Enable Intensity capping
% %  --> Wiener2 denoise filter with 3 pix
% %  --> DO NOT Auto constrast stretch
% % PIV settings: 
% %  --> 128 (32 step), 64 (32 step), 32 (16 step), 16 (8 step)
% %  --> Linear window deformation interpolator
% %  --> 5x repeated correlation 
% %  --> Disable auto-correlation
% % Post-processing
% %  --> Standard deviation filter: 7 stdev
% %  --> Local median filter: thres=5, eps=0.1
% %  --> Interpolate missing data
% % Export 
% %  --> File > Save > MAT file
disp('Loading PIV results...')
tmp = load(fullfile(QS.dir.piv, 'piv_results.mat')) ;


%% Measure velocities =============================================
disp('Making map from pixel to xyz to compute velocities in 3d for smoothed meshes...')
options = struct() ;
options.overwrite = false ;
options.preview = false ;
options.show_v3d_on_data = false ;
options.save_ims = true ;
QS.measurePIV3d(options) ;


%% AUTOCORRELATIONS
% overwrite_autocorrelations = false ;
% do_acorr = false ;
% redo_acorr = ~exist(fullfile(pivSimAvgDir, 'autocorr_velocities.png'), 'file') ;
% if (redo_acorr || overwrite_autocorrelations) && do_acorr
%     aux_autocorrelations
% end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% First do very simpleminded averaging of velocities
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options.overwrite = false ;
% options.plot_vxyz = true ;
% QS.timeAverageVelocitiesSimple('1x', options)
% %% VELOCITY PLOTS
% options.overwrite = true ;
% options.plot_vxyz = false ;
% options.invertImage = true ;
% QS.plotTimeAvgVelSimple(options)
% %% Divergence and Curl (Helmholtz-Hodge)
% options = struct() ;
% options.overwrite = false ;
% options.samplingResolution = '1x' ; 
% options.averagingStyle = 'simple' ;
% QS.helmholtzHodge(options) ;
% %% Measure Compressibility (div(v), 2*vn*H, and gdot)
% options = struct() ;
% options.overwrite = false ;
% options.plot_Hgdot = false ;
% options.plot_flows = true ;
% options.plot_factors = true ;
% options.plot_kymographs = true ;
% options.plot_kymographs_cumsum = true ;
% options.plot_correlations = true ;
% options.plot_gdot_correlations = false ;
% options.plot_gdot_decomp = true ;
% options.lambda_mesh = 0.002 ;
% options.lambda = 0.02 ;
% options.lambda_err = 0.03 ;
% options.samplingResolution = '1x'; 
% QS.measureMetricKinematics(options)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% DOUBLE RESOLUTION Simple average
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Measure velocities at 2*nU x 2*nV resolution
% options = struct() ;
% options.overwrite = false ;
% options.preview = false ;
% options.show_v3d_on_data = false ;
% options.save_ims = true ;
% QS.generateSPCutMeshSm2x(options.overwrite) ;
% QS.measurePIV3d2x(options) ;
% %% time average the velocities
% options.overwrite = false ;
% options.plot_vxyz = true ;
% QS.timeAverageVelocitiesSimple('2x', options) ;
% %% Velocity plots for 2x
% options.overwrite = false ;
% options.plot_vxyz = false ;
% options.invertImage = true ;
% options.samplingResolution = '2x'; 
% QS.plotTimeAvgVelSimple(options)
% %% Divergence and Curl (Helmholtz-Hodge) for 2x
% options = struct() ;
% options.overwrite = true ;
% options.samplingResolution = '2x' ;
% QS.helmholtzHodgeSimple(options) ;
% %% doubleResolution compressibility
% options = struct() ;
% options.overwrite = false ;
% options.plot_Hgdot = false ;
% options.plot_flows = true ;
% options.plot_factors = true ;
% options.plot_kymographs = true ;
% options.plot_kymographs_cumsum = true ;
% options.plot_correlations = true ;
% options.plot_gdot_correlations = false ;
% options.plot_gdot_decomp = true ;
% options.samplingResolution = '2x'; 
% QS.measureMetricKinematics(options)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Lagrangian dynamics
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Pullback pathline time averaging of velocities
options = struct() ;
options.overwrite = false ;
options.preview = false ;
QS.timeAverageVelocities(options)
%% Velocity plots for pathline time averaging 
options.overwrite = false ;
options.plot_vxyz = false ;
options.invertImage = true ;
options.averagingStyle = 'Lagrangian'; 
options.samplingResolution = '1x'; 
QS.plotTimeAvgVelocities(options)
%% Divergence and Curl (Helmholtz-Hodge) for Lagrangian
options = struct() ;
options.overwrite = false ;
options.samplingResolution = '1x' ;
options.averagingStyle = 'Lagrangian' ;
options.lambda = 0 ;
options.lambda_mesh = 0 ; 
QS.helmholtzHodge(options) ;

%% Compressibility & kinematics for Lagrangian
options = struct() ;
options.overwrite = false ;
options.samplingResolution = '1x'; 
options.lambda_mesh = 0 ;
options.lambda_err = 0.00 ;
options.lambda = 0.00 ;
options.nmodes = 7 ;  %% bandwidth filtering
options.zwidth = 2 ; 
QS.measureMetricKinematics(options)

%% Contraction rate to compare with GCaMP
options = struct() ;
QS.plotMetricKinematicsEulerianFrame(options)

%% Metric Kinematics Kymographs & Correlations -- Bandwidth Filtered
options = struct() ;
options.overwrite = false ;
options.overwrite_timePoints = false ;
options.plot_Hgdot = false ;
options.plot_flows = true ;
options.plot_factors = false ;
options.plot_kymographs = false ;
options.plot_kymographs_cumsum = false ;
options.plot_kymographs_cumprod = false ;
options.plot_correlations = false ;
options.plot_spaceMaps = true ;
options.plot_raw_correlations = false ;
options.plot_gdot_correlations = false ;
options.plot_gdot_decomp = false ;
options.lambda = 0.00 ;
options.lambda_err = 0.00 ;
options.lambda_mesh = 0 ;
options.nmodes = 7 ;  %% bandwidth filtering
options.zwidth = 1 ; 
options.climit = 0.3 ;
QS.plotMetricKinematics(options)

%% Measure shear stresses against gradients in pressure
% eta \nabla_i d^i_j = \nalba_j p, 
% where d^i_j = 1/2(\nabla_i v_j + \nalba_j v_i) - b_ij vn
options = struct() ;
options.overwrite = true ;
options.zwidth = 0 ;
QS.measureStokesForces(options) ;
options.plot_kymographs = true ;
QS.plotStokesForces(options) ;  % !!! continue here



%% Pullback pathlines connecting Lagrangian grids
options = struct() ;
options.overwrite = false ;
options.preview = false ;
options.debug = false ; 
QS.measurePullbackPathlines(options)

%% Query velocities along pathlines
options = struct() ;
options.overwrite = false ;
options.preview = false ;
QS.measurePathlineVelocities(options)
% plot the pathline velocities 
options = struct() ;
options.overwrite = false ;
options.gridTopology = 'triangulated' ;
QS.plotPathlineVelocities(options)

%% Measure Pathline Kinematics
% todo here
options = struct() ;
options.overwrite = false ;
options.lambda = 0.01 ;
options.lambda_error = 0.01 ;
options.lambda_mesh = 0 ;
options.nmodes = 7 ;
options.zwidth = 1 ; 
QS.measurePathlineMetricKinematics(options)
% todo: resume this BLOCK !!!

%% Plot Pathline Kinematics
options = struct() ;
options.overwrite = true ;
options.plot_kymographs = true ;
options.plot_kymographs_cumsum = false ;
options.plot_kymographs_cumprod = false ;
options.plot_correlations = true ;
options.plot_fold_kinematics = true ;
options.plot_lobe_kinematics = true ;
options.climit = 0.30 ;
QS.plotPathlineMetricKinematics(options)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Create ricci mesh at t0 to measure Beltrami coefficient in pathlines
options = struct() ;
options.overwrite = false ;
options.climit = 1 ;
options.coordSys = 'ricci' ;
QS.measureBeltramiCoefficient(options) ;

%% Generate all Beltramis from all Riccis & plot aspect ratio over time
options = struct() ;
options.overwrite = false ;
options.resample = false ;
QS.computeRicciMeshes(options)
% return to here

%% Generate all Beltramis from all Riccis & Plot aspect ratio for isothermal PB over time
% Note: this really shouldn't be necessary, as we show in detail now
options = struct() ;
QS.measureBeltramiCoefficientPullbackToPullback(options) ;

%% Check mu by computing for all times and get the same thing as 
% embedding back to same Reference tp

for tp = QS.xp.fileMeta.timePoints
    disp(['t = ', num2str(tp)])
    QS.setTime(tp)
    try
        rmesh = QS.loadCurrentRicciMesh() ;
    catch
        disp(['could not load Ricci mesh for t=' num2str(tp)])
    end
    
    % Compute mu from mapping from PB(0) to PB(1)
    
end


%% Compare to linearized description
options = struct() ;
options.overwrite = false ;
options.overwriteImages = false ;
options.coordSys = 'ricci' ;
QS.compareBeltramiToLinearizedConstriction(options) ;

%% Compare to nonlinear isothermal description
options = struct() ;
options.overwrite = true ;
options.overwriteImages = false ;
options.coordSys = 'ricci' ;
QS.compareBeltramiToConstriction(options) ;

%% UVPrime option (ignore)
master_uvprime_module
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Measure twist (d v_phi / d zeta)
options = struct() ;
options.overwrite = false ;
QS.measureTwist(options)

%% Measure required stress pattern
options = struct() ;
options.overwrite = false ;
QS.measureStressPattern(options) ;

%% Measure metric strain rate (epsilon=gdot/2)
% option = struct() ;
% options.overwrite = true ;
% options.preview = false ;
% QS.measureMetricStrainRate(options) 

%% Strain rate (epsilon = 1/2 (djvi+divj) -vn bij)
options = struct() ;
options.overwrite = false ;
options.overwriteImages = false ;
options.preview = false ;
options.lambda = 0.000 ;
options.lambda_error = 0 ;
options.lambda_mesh = 0 ;
options.nmodes = 7 ;
options.zwidth = 1 ; 
QS.measureStrainRate(options) 

%% Plot time-averaged strain rates in 3d on mesh
options = struct() ;
% error('This does not look great -- tweak how the averaging is done to be Lagrangian?')
options.overwrite = false ;
options.preview = false ;
QS.plotStrainRate3DFiltered(options) 

%% Kymograph strain rates
options = struct() ;
options.overwrite = false ;
options.skipTimePoint = true ;
options.clim_trace = 0.05 ;
options.clim_deviatoric = 0.05 ;
QS.plotStrainRate(options)

% Measure strain rate along pathlines
options = struct() ;
options.overwrite = false ;
options.overwriteImages = false ;
options.plot_dzdp = false ;
QS.measurePathlineStrainRate(options)
% Integrate strain rates along pathlines (noisy)
% options = struct() ;
% options.overwrite = true ;
% options.overwriteImages = true ;
% options.plot_dzdp = false ;
% options.median_filter_strainRates = false ;
% options.climitInitial = 0.05 ;
% options.climitRamp = 0.01 ;
% options.climitRatio = 1 ;
% QS.measurePathlineIntegratedStrain(options)

% Pathline strain rate plots
options = struct() ;
options.overwrite = false ;
options.plot_kymographs = false ;
options.plot_fold_strainRate = true ;
options.plot_lobe_strainRate = true ;
options.climit = 0.05 ;
options.climitWide = 1.0 ;
QS.plotPathlineStrainRate(options)

%% (REPEAT FROM ABOVE) Pullback pathlines connecting Lagrangian grids
% options = struct() ;
% options.overwrite = false ;
% options.preview = false ;
% options.debug = false ; 
% QS.measurePullbackPathlines(options)
%% Measure strain along pathlines -- note this is from pathlines, not integrating rates
options = struct() ;
options.overwrite = true ;
options.overwriteImages = false ;
options.plot_dzdp = false ;
options.climitInitial = 0.05 ;
options.climitRamp = 0.01 ;
options.climitRatio = 1 ;
QS.measurePathlineStrain(options)
QS.plotPathlineStrain(options)

%% Measure coarse-grained bond contraction and dilation in zeta=s/L, phi
% todo: consider putting this in Ricci map frame instead of (s,phi) frame
options = struct() ;
options.overwrite = true ;
options.overwriteImages = true ;
QS.measureDxDyStrainFiltered(options) ;

%% Output gg and bb in fixed frame (t0 Lagrangian zeta phi frame) -- this is done in measurePathlineStrain()
% options = struct() ;
% options.overwrite = false ;
% measurePathlineMetric(QS, options)

%% Simulate this experiment
options = struct() ;
QS.simulateNES(options)

%% Measure surface area growh of lobes and folds in Eulerian frame
options = struct() ;
options.overwrite = false ;
options.preview = false ;
QS.measureEulerianMetricDynamics(options)

%% Plot texture map in 3d with velocities
for i = 1:size(vM, 1)
    t = time(i) ;
    vtqfn0 = fullfile(pivSimAvgImQDir, [sprintf('%04d', time(i)) '_perspective.png']) ;
    vtqfn1 = fullfile(pivSimAvgImQDir, [sprintf('%04d', time(i)) '_dorsal.png']) ;
    vtqfn2 = fullfile(pivSimAvgImQDir, [sprintf('%04d', time(i)) '_ventral.png']) ;
    vtqfn3 = fullfile(pivSimAvgImQDir, [sprintf('%04d', time(i)) '_lateralL.png']) ;
    vtqfn4 = fullfile(pivSimAvgImQDir, [sprintf('%04d', time(i)) '_lateralR.png']) ;

    % Plot the tangential velocity as quiver on top of the embedding,
    % colored by normal velocity (all in one in embedding!)
    if ~exist(vtqfn0, 'file') || overwrite_vsm_plots
        disp(['Creating ' vtqfn])
        tic

        v2dsmum_ii = squeeze(v2dsmMum(i, :, :)) ;
        vnsm_ii = squeeze(vnsmM(i, :, :)) ;

        % Load the data for the current time point ------------------------
        xp.setT=ime(t) ;
        % Load 3D data for coloring mesh pullback
        xp.loadTime(t);
        xp.rescaleStackToUnitAspect();

        fig = figure('units', 'normalized', ...
                'outerposition', [0 0 1 1], 'visible', 'off') ;

        % Psize is the linear dimension of the grid drawn on each triangular face
        Options.PSize = 5;
        Options.EdgeColor = 'none';
        Options.Rotation = rot ;
        Options.FaceScalarField = vn_interp ;

        % Raw stack data
        IV = xp.stack.image.apply();
        IV = imadjustn(IV{1});

        % First args are physical vertices, then texture faces (same number as 
        % physical faces, but connectivity can be different), then texture
        % vertices, which can also be different. The vertices are in row, column, 
        % page format, so x and y get flipped. IV is the texture volume.
        % Options.PSize 
        [ patchIm, imref, zeroID, MIP, SIP ] = ...
            texture_patch_3d( mesh.f, mesh.v, ...
            mesh.f, mesh.v(:, [2 1 3]), IV, Options );

        % Add quiver to image
        % Downsample by a factor of qsubsample
        xx = piv.x{i}(1, :)' ;
        yy = piv.y{i}(:, 1) ;
        ww = length(xx) ;
        hh = length(yy) ;
        vx = reshape(v2dsmum_ii(:, 1), [ww, hh]) ;
        vy = reshape(v2dsmum_ii(:, 2), [ww, hh]) ;
        QX = imresize(vx, [ww / qsubsample, hh / qsubsample], 'bicubic') ;
        QY = imresize(vy, [ww / qsubsample, hh / qsubsample], 'bicubic') ;
        xq = 1:qsubsample:ww ;
        yq = 1:qsubsample:hh ;
        [xg, yg] = meshgrid(xq, yq) ;
        h2 = quiver(xg(:), yg(:), QX(:), QY(:), 0) ;

        % Format the figure and axis
        xlim(xyzmin(1, :))
        ylim(xyzmin(2, :))
        zlim(xyzmin(3, :))
        title(['$t = $' num2str(t) ' min'], 'Interpreter', 'Latex', 'Color', 'white') 
        xlabel('AP position [$\mu$m]', 'Interpreter', 'Latex', 'Color', 'white')
        ylabel('DV position [$\mu$m]', 'Interpreter', 'Latex', 'Color', 'white')
        zlabel('lateral position [$\mu$m]', 'Interpreter', 'Latex', 'Color', 'white')
        set(fig, 'PaperUnits', 'centimeters');
        set(fig, 'PaperPosition', [0 0 xwidth ywidth]);

        % Make background black
        set(gca,'Color','k')
        set(gcf, 'InvertHardCopy', 'off');
        set(gcf, 'Color', 'k')

        % Make tick labels white
        labeltype = {'XTickLabel', 'YTickLabel', 'ZTickLabel'} ;
        for q = 1:2
            % get the current tick labeks
            ticklabels = get(gca, labeltype{q});
            % prepend a color for each tick label
            ticklabels_new = cell(size(ticklabels));
            for i = 1:length(ticklabels)
                ticklabels_new{i} = ['\color{white} ' ticklabels{i}];
            end
            % set the tick labels
            set(gca, labeltype{q}, ticklabels_new);
        end

        % Capture all five views (perspective + DVLR
        disp(['saving figure...' num2str(t, '%06d')])
        saveas(fig, vtqfn0)
        view(0, 90)
        % DORSAL
        saveas(fig, vtqfn1)
        view(0, 270)
        % VENTRAL
        saveas(fig, vtqfn2)

        % Lateral views
        view(0, 0)
        % make white z tick labels
        q = 3;
        % get the current tick labeks
        ticklabels = get(gca, labeltype{q});
        % prepend a color for each tick label
        ticklabels_new = cell(size(ticklabels));
        for i = 1:length(ticklabels)
            ticklabels_new{i} = ['\color{white} ' ticklabels{i}];
        end
        % set the tick labels
        set(gca, labeltype{q}, ticklabels_new);

        % LEFT
        saveas(fig, vtqfn3)
        view(0, 180)
        % RIGHT
        saveas(fig, vtqfn4)
        close all
        toc
        clear Options IV
    end

end

% Check the orientation of the phasebar
% imshow(im)
% hold on;
clf
xsize = 2000 ;
ysize = 2000 ;
step = 100 ;
imsz = [ xsize ysize ] ;
[xx, yy] = meshgrid(1:step:xsize, 1:step:ysize) ;
ucheck = xx - xsize * 0.5 ;
vcheck = yy - ysize * 0.5 ;
vangle = mod(atan2(vcheck, -ucheck), 2* pi) ;
imshow(im * washout2d + max(im) * (1-washout2d)) ;
hold on ;
h2 = imagesc(xx(:), yy(:), vangle) ;
hold on ;
quiver(xx, yy, ucheck, vcheck, 5) ;
colormap phasemap
phasebar
caxis([0, 2*pi])
axis on
set(gcf, 'visible', 'on')
waitfor(gcf)


%% Simpleminded streamlines from velocity scaled by dilation
% % Build 3d grid of positions and velocities
% % Assume that x0, y0 are fixed for all time
% ntimes = length(piv.x) ;
% x0 = piv.x{i} ;
% y0 = piv.y{i} ;
% xlen = size(v2dsmM, 1) ;
% ylen = size(v2dsmM, 2) ;
% vdat = v2dsmM(:) ;
% vxdat = reshape(vdat(:, 1), [ntimes, xlen, ylen]) ;
% vydat = reshape(vdat(:, 2), [ntimes, xlen, ylen]) ;
% vzdat = ones(size(vydat)) ; % we march through time at 1 index / timestep
% % Define positions we track through the streamlines
% startx = x0(1:200:end) ;
% starty = y0(1:200:end) ;
% startz = zeros(size(starty)) ;
% streamline(x0, y0, z0, vxdat, vydat, vzdat, startx, starty, startz)
% view(2)
% save(gcf, fullfile(pivDir, 'streamlines.png')) 
% error('break')

back to top