diff --git a/+Ecto/compare_manual_ecto.m b/+Ecto/compare_manual_ecto.m
new file mode 100644
index 0000000000000000000000000000000000000000..e4fde8d63e52a9fc3b4c72c2c230e87e9aba4983
--- /dev/null
+++ b/+Ecto/compare_manual_ecto.m
@@ -0,0 +1,84 @@
+function compare_manual_ecto(datasetName, imageRoot)
+ if ( ~exist('Data','dir') )
+ return;
+ end
+
+ maskInfo = struct('name',{'ecto','endo'}, 'file',{'',''});
+
+ maskInfo(1).file = fullfile('Data','ectomasks',[datasetName '_ectoregion.mat']);
+ maskInfo(2).file = fullfile('Data','ectomasks',[datasetName '_endoregion.mat']);
+
+ if ( ~exist(maskInfo(1).file,'file') )
+ warning([datasetName ': Ectoderm mask doesn''t exist -- Skipping']);
+ return;
+ end
+
+ ncvDir = fullfile('Data', 'ectodeltas');
+ if ( ~exist(ncvDir, 'dir') )
+ mkdir(ncvDir);
+ end
+
+ jsonFile = fullfile(ncvDir,[datasetName '.json']);
+ if ( exist(jsonFile ,'file') )
+ return;
+ end
+
+ maskDeltas = struct('type',{}, 'deltas',{}, 'ncvs',{});
+ for i=1:length(maskInfo)
+ if ( ~exist(maskInfo(i).file,'file') )
+ continue;
+ end
+
+ imPath = getImPath(imageRoot, datasetName, {'matlab_decon','CPPdecon'});
+ im1 = MicroscopeData.Reader(imPath, 'timeRange',[1,1], 'chanList',[1]);
+
+ s = load(maskInfo(i).file);
+ mask = s.regionInfo.bwMask(:,:,:,1);
+
+ try
+ [deltas, ncvs] = getNcvDeltas(im1,mask, imPath, 50);
+ catch me
+ warning([datasetName ': Skipping -- ' me.message]);
+ return;
+ end
+
+ maskDeltas(i).type = maskInfo(i).name;
+ maskDeltas(i).deltas = deltas;
+ maskDeltas(i).ncvs = ncvs;
+ end
+
+ jsonData = Utils.CreateJSON(maskDeltas);
+ fid = fopen(jsonFile, 'wt');
+ fprintf(fid, '%s\n', jsonData);
+ fclose(fid);
+end
+
+function imPath = getImPath(imageRoot, datasetName, imSubdirs)
+ imPath = fullfile(imageRoot,datasetName);
+ for i=1:length(imSubdirs)
+ if ( exist(fullfile(imPath,imSubdirs{i}), 'dir') )
+ imPath = fullfile(imPath,imSubdirs{i});
+ return;
+ end
+ end
+end
+
+function [deltas, ncvs] = getNcvDeltas(im, mask, imPath, maxFrame)
+ imD = MicroscopeData.ReadMetadataFile(fullfile(imPath,filesep));
+
+ maxFrame = min(maxFrame, imD.NumberOfFrames);
+
+ prg = Utils.CmdlnProgress(maxFrame, true, [imD.DatasetName ' -- Computing Deltas']);
+
+ deltas = NaN(maxFrame,3);
+ ncvs = NaN(maxFrame,1);
+ for t=1:maxFrame
+ imt = MicroscopeData.Reader(imPath, 'timeRange',[t,t], 'chanList',[1]);
+
+ [chkDelta, chkNCV] = Ecto.get_masked_ncv(im,mask, imt);
+ deltas(t,:) = chkDelta;
+ ncvs(t) = chkNCV;
+
+ prg.PrintProgress(t);
+ end
+end
diff --git a/+Ecto/ectoSeg.m b/+Ecto/ectoSeg.m
new file mode 100644
index 0000000000000000000000000000000000000000..2a457a2698019189027b8f725853395cf0df80cc
--- /dev/null
+++ b/+Ecto/ectoSeg.m
@@ -0,0 +1,86 @@
+% [im,imD] = MicroscopeData.Reader('H:/Smadar/Raw/22-06-2016_TimeLapse3_DMSO_20hrs_Calcein_FM464/matlab_decon', 'timeRange',[1,1]);
+[im,imD] = MicroscopeData.Reader('C:/Users/mwinter/Documents/Lab/Smadar/Raw/22-06-2016_TimeLapse3_DMSO_20hrs_Calcein_FM464/matlab_decon', 'timeRange',[1,1]);
+
+% addpath('C:\Users\mwinter\Documents\MATLAB\external\AOSLevelsetSegmentationToolboxM');
+pxScale = imD.PixelPhysicalSize;
+
+% mul = 20;
+bwMask = llsm_mask(im);
+
+se = strel('square', 31);
+bwShrink = false(size(bwMask));
+for z=1:size(bwMask,3)
+ bwShrink(:,:,z) = imerode(bwMask(:,:,z), se);
+end
+
+imFM = mat2gray(im(:,:,:,1));
+
+alpha = 0.95;
+thresh = alpha*graythresh(imFM);
+bwFM = imFM >= thresh;
+bwFM(~bwShrink) = false;
+
+% Approximate ectoderm thickness
+ectoDist = 15.0;
+idxPts = find(bwFM);
+coords_xy = Utils.SwapXY_RC(Utils.IndToCoord(size(bwFM), idxPts));
+coords_um = coords_xy.*pxScale;
+
+shp = alphaShape(coords_um(:,1),coords_um(:,2),coords_um(:,3), ectoDist);
+[f,pts_xyz] = boundaryFacets(shp);
+
+unscalePts = pts_xyz ./ pxScale;
+
+poly = D3d.Polygon.MakeEmptyStruct();
+nrms = D3d.Polygon.CalcNorms(unscalePts, f);
+poly(1).index = 1;
+poly(1).frame = 1;
+poly(1).label = '1';
+poly(1).color = [1,0,0];
+poly(1).faces = f;
+poly(1).verts = unscalePts - 0.5;
+poly(1).norms = nrms;
+poly(1).CenterOfMass = mean(unscalePts);
+
+D3d.Viewer.AddPolygons(poly);
+
+% coords = Utils.IndToCoord(size(bwMask), find(bwMask));
+%
+% init_phi = -ones(size(imFM));
+% ctr = round(mean(coords));
+% padding = 20;
+% padctr = {ctr(1)-padding:ctr(1)+padding, ctr(2)-padding:ctr(2)+padding, ctr(3)-padding:ctr(3)+padding};
+% init_phi(padctr{:}) = 1;
+%
+% init_phi = ac_reinit(init_phi);
+%
+% smooth_weight = 3;
+% image_weight = 1e-3;
+% delta_t = 2;
+% n_iters = 100;
+% show_result = 1;
+%
+% phi = ac_ChanVese_model(imFM, init_phi, smooth_weight, image_weight, delta_t, n_iters, show_result);
+
+
+% bwFM = mask;
+% updatePolys(bwFM);
+%
+% for i=1:500
+% bwFM = activecontour(imBMask, bwFM, 2, 'edge', 'ContractionBias', -20.0*mul, 'SmoothFactor',mul);
+%
+% updatePolys(bwFM);
+% end
+
+function updatePolys(bwFM)
+ rp = regionprops(bwFM, 'PixelList');
+ polys = [];
+ for i=1:length(rp)
+ polys = [polys; D3d.Polygon.Make(rp(i).PixelList, i, ['Endo_' num2str(i)], 1, [1 0 0], 1)];
+ end
+
+ D3d.Viewer.DeleteAllPolygons();
+ D3d.Viewer.AddPolygons(polys);
+
+ D3d.Update();
+end
diff --git a/+Ecto/get_masked_ncv.m b/+Ecto/get_masked_ncv.m
new file mode 100644
index 0000000000000000000000000000000000000000..3097226ba6a59e906d51efdf3a1221bd358e4096
--- /dev/null
+++ b/+Ecto/get_masked_ncv.m
@@ -0,0 +1,8 @@
+function [delta, ncv] = get_masked_ncv(frame1Im, frame1Ecto, frametIm)
+ frametMask = Segment.llsm_mask(frametIm);
+
+ [dX,dY,dZ,ncv,overlapSize,ncvMatrixROI] = Registration.FFT.RegisterTwoImages(frame1Im,[],frametIm,[],...
+ [],[],[],[], false, frame1Ecto, frametMask);
+
+ delta = [dX;dY;dZ];
+end
diff --git a/+Export/get_subdir_list.m b/+Export/get_subdir_list.m
new file mode 100644
index 0000000000000000000000000000000000000000..72b28812d50d47bba3424b56ac6d31ca4ebdf315
--- /dev/null
+++ b/+Export/get_subdir_list.m
@@ -0,0 +1,6 @@
+function subdirs = get_subdir_list(rootDir)
+ flist = dir(rootDir);
+ bValidDir = arrayfun(@(x)((x.isdir>0) && ~strncmpi(x.name,'.',1) && ~strncmpi(x.name,'$',1)), flist);
+
+ subdirs = {flist(bValidDir).name}.';
+end
\ No newline at end of file
diff --git a/+Export/smadar_export_datasets.m b/+Export/smadar_export_datasets.m
new file mode 100644
index 0000000000000000000000000000000000000000..d1fccbd48ff43a5ad82148339d2c7a3b74d8b5a7
--- /dev/null
+++ b/+Export/smadar_export_datasets.m
@@ -0,0 +1,14 @@
+IN_ROOT = 'E:/';
+OUT_ROOT = 'H:/Smadar/Raw';
+
+dateDirs = Export.get_subdir_list(IN_ROOT);
+for i=1:length(dateDirs)
+ dataDirs = Export.get_subdir_list(fullfile(IN_ROOT,dateDirs{i}));
+ for j=1:length(dataDirs)
+ outName = [dateDirs{i} '_' dataDirs{j}];
+
+ dirIn = fullfile(IN_ROOT, dateDirs{i}, dataDirs{j});
+ dirOut = fullfile(OUT_ROOT, outName);
+ LLSM.ConvertToKLB(dirIn, dirOut, {'matlab_decon', 'CPPdecon', 'Deskewed', 'deskewed'}, outName);
+ end
+end
diff --git a/+Fixup/add_default_center_param.m b/+Fixup/add_default_center_param.m
new file mode 100644
index 0000000000000000000000000000000000000000..7392d64a4c128dd45aef8d7b369ea7bd87b6a46c
--- /dev/null
+++ b/+Fixup/add_default_center_param.m
@@ -0,0 +1,39 @@
+function add_default_center_param(ROOT)
+ if ( ~exist('ROOT','var') )
+ ROOT = pwd();
+ end
+
+ flist = dir(fullfile(ROOT,'*params.json'));
+ for i=1:length(flist)
+ filepath = fullfile(ROOT,flist(i).name);
+ jsonparams = fileread(filepath);
+
+ s = jsondecode(jsonparams);
+ fieldList = fieldnames(s);
+ if ( ~all(ismember({'render','view'}, fieldList)) )
+ continue;
+ end
+
+ if ( ~ismember({'center'},fieldnames(s.view)) )
+ s.view = add_center(s.view);
+ end
+
+ jsonData = Utils.CreateJSON(s, true);
+
+% tok = regexp(flist(i).name, '(.*?)_params\.json', 'once','tokens');
+% movieName = tok{1};
+
+% fid = fopen(fullfile(ROOT,[movieName '_params.json']), 'wt');
+ fid = fopen(filepath,'wt');
+ fprintf(fid, '%s\n', jsonData);
+ fclose(fid);
+ end
+end
+
+
+function newview = add_center(view)
+ newview = view;
+ newview.center = [0.0;0.0;0.0];
+
+ newview = orderfields(newview, {'zoom','pos','center','worldRot','bClip','clipPlane'});
+end
diff --git a/+Fixup/fixup_im_paths.m b/+Fixup/fixup_im_paths.m
new file mode 100644
index 0000000000000000000000000000000000000000..051e2e8cf62fd2853986d6e0d55a92317005ff38
--- /dev/null
+++ b/+Fixup/fixup_im_paths.m
@@ -0,0 +1,36 @@
+% Load lever files and retarget imageData.imageDir and ui.rootImageFolder
+% NOTE: The special variable <DATASETNAME> is replaced in patterns
+function fixup_im_paths(leverPath, uiImPatterns, rawImPatterns)
+ flist = dir(fullfile(leverPath, '*.LEVER'));
+ for i=1:length(flist)
+ updateConstants(fullfile(leverPath,flist(i).name), uiImPatterns, rawImPatterns);
+ end
+end
+
+function updateConstants(leverFile, uiImPatterns, rawImPatterns)
+ conn = database(leverFile, '','', 'org.sqlite.JDBC', 'jdbc:sqlite:');
+ CONSTANTS = Read.getConstants(conn);
+
+ datasetName = CONSTANTS.imageData.DatasetName;
+
+ chkUIPaths = strrep(uiImPatterns, '<DATASETNAME>', datasetName);
+ bUIExist = cellfun(@(x)(exist(x,'dir') > 0), chkUIPaths);
+
+ validUIPaths = chkUIPaths(bUIExist);
+ if ( isempty(validUIPaths) )
+ error(['No UI image patterns match real directories for dataset: ' datasetName]);
+ end
+
+ chkRawPaths = strrep(rawImPatterns, '<DATASETNAME>', datasetName);
+ bRawExist = cellfun(@(x)(exist(x,'dir') > 0), chkRawPaths);
+
+ validRawPaths = chkRawPaths(bRawExist);
+ if ( isempty(validRawPaths) )
+ error(['No raw image patterns match real directories for dataset: ' datasetName]);
+ end
+
+ CONSTANTS.imageData.imageDir = validRawPaths{1};
+ CONSTANTS.ui.rootImageFolder = validUIPaths{1};
+
+ Write.updateConstants(conn, CONSTANTS);
+end
diff --git a/+Fixup/fixup_movie_params.m b/+Fixup/fixup_movie_params.m
new file mode 100644
index 0000000000000000000000000000000000000000..59e3bb6bc1d57fb9c109616e9af209ecb9063570
--- /dev/null
+++ b/+Fixup/fixup_movie_params.m
@@ -0,0 +1,38 @@
+function fixup_movie_params(ROOT)
+ if ( ~exist('ROOT','var') )
+ ROOT = pwd();
+ end
+
+ flist = dir(fullfile(ROOT,'*.mat'));
+ for i=1:length(flist)
+ filepath = fullfile(ROOT,flist(i).name);
+ vlist = whos('-file', filepath);
+
+ if ( ~all(ismember({'ui','renderparams','viewparams'}, {vlist.name})) )
+ continue;
+ end
+
+ s = load(filepath);
+ if ( ~ismember({'clipPlane'}, fieldnames(s.viewparams)) )
+ s.viewparams = add_clipinfo(s.viewparams);
+ end
+
+ outparams = struct('render',{s.renderparams}, 'view',{s.viewparams});
+
+ jsonData = Utils.CreateJSON(outparams, true);
+
+ [~,movieName,ext] = fileparts(flist(i).name);
+ fid = fopen(fullfile(ROOT,[movieName '_params.json']), 'wt');
+ fprintf(fid, '%s\n', jsonData);
+ fclose(fid);
+ end
+end
+
+
+function newview = add_clipinfo(view)
+ newview = view;
+ newview.bClip = false;
+ newview.clipPlane = 0.0;
+
+ newview = orderfields(newview, {'zoom','pos','worldRot','bClip','clipPlane'});
+end
diff --git a/+Fixup/fixup_smadar_impaths.m b/+Fixup/fixup_smadar_impaths.m
new file mode 100644
index 0000000000000000000000000000000000000000..749a12bb67d40306de0e6a2f1243366095a9f963
--- /dev/null
+++ b/+Fixup/fixup_smadar_impaths.m
@@ -0,0 +1,6 @@
+rawPatterns = {'E:/Smadar/Raw/<DATASETNAME>/CPPdecon';
+ 'E:/Smadar/Raw/<DATASETNAME>/matlab_decon'};
+
+uiPatterns = {'E:/Smadar/LEVERJS/<DATASETNAME>'};
+
+fixup_im_paths('E:/Smadar/LEVERJS', uiPatterns, rawPatterns);
diff --git a/+Fixup/rename_all_frames.m b/+Fixup/rename_all_frames.m
new file mode 100644
index 0000000000000000000000000000000000000000..b8fa4f9503ee688a4a1af40d3131e986e89c63bf
--- /dev/null
+++ b/+Fixup/rename_all_frames.m
@@ -0,0 +1,11 @@
+ROOT='Data';
+
+flist = dir(fullfile(ROOT,'*_frames'));
+
+bValidDirs = arrayfun(@(x)(x.isdir>0), flist);
+frameDirs = flist(bValidDirs);
+
+for i=1:length(frameDirs)
+ movieTok = regexp(frameDirs(i).name, '([-\w]+?)_frames', 'tokens', 'once');
+ rename_frame_files(movieTok{1});
+end
diff --git a/+Fixup/rename_frame_files.m b/+Fixup/rename_frame_files.m
new file mode 100644
index 0000000000000000000000000000000000000000..d25c199c2ac6f35ade67963f8f85103d364d5ab2
--- /dev/null
+++ b/+Fixup/rename_frame_files.m
@@ -0,0 +1,35 @@
+function rename_frame_files(movieName)
+ dataDir = 'Data';
+
+ %% Setup movie framce-cap dir
+ frameDir = fullfile(dataDir, [movieName '_frames']);
+ if ( ~exist(frameDir, 'dir') )
+ mkdir(frameDir);
+ end
+
+ flist = dir(fullfile(frameDir,'*.png'));
+ fnames = {flist.name};
+
+ tok = regexpi(fnames, '([-\w]+?)_t(\d+)\.png','tokens', 'once');
+
+ bValid = cellfun(@(x)(~isempty(x)), tok);
+
+ validNames = fnames(bValid);
+ validTok = tok(bValid);
+
+ ftimes = cellfun(@(x)(str2double(x{2})), validTok);
+ firstT = min(ftimes);
+
+ if ( firstT == 0 )
+ fprintf('---- %s: Already frame-indexed: skipping\n', frameDir);
+ return;
+ end
+
+ fprintf('---- %s: Renaming frames\n', frameDir);
+ for i=1:length(validNames)
+ newName = [validTok{i}{1} '_t' num2str(ftimes(i)-firstT,'%04d') '.png'];
+ fprintf('Renaming %s - %s\n', validNames{i},newName);
+
+ movefile(fullfile(frameDir,validNames{i}), fullfile(frameDir,newName));
+ end
+end
\ No newline at end of file
diff --git a/+LLSM/ConvertToKLB.m b/+LLSM/ConvertToKLB.m
new file mode 100644
index 0000000000000000000000000000000000000000..38691486551e5d85c45bbdeda26beccbf70c8197
--- /dev/null
+++ b/+LLSM/ConvertToKLB.m
@@ -0,0 +1,171 @@
+function ConvertToKLB(dirIn,dirOut, subfolders, forceDatasetName)
+ settingsList = dir(fullfile(dirIn,'*_Settings.txt'));
+ firstMetaSettings = LLSM.ParseSettingsFile(fullfile(dirIn,settingsList(1).name));
+
+ if ( ~exist('forceDatasetName','var') )
+ forceDatasetName = '';
+ end
+
+ %% Cleanup function to try and avoid partial file-writes
+ fileMap = containers.Map();
+ if ( ~isempty(dirOut) )
+ cleanupObj = onCleanup(@()(cleanupFunc(fileMap)));
+ end
+
+ for s = 1:length(subfolders)
+ fullOutDir = fullfile(dirOut,subfolders{s});
+ if (~exist(fullfile(dirIn, subfolders{s}),'dir'))
+ warning('Cannot find %s.\nSkipping\n',fullfile(dirIn,subfolders{s}));
+ continue
+ end
+
+ fprintf('Converting %s\n', fullfile(dirIn,subfolders{s}));
+ if ( ~exist(fullOutDir, 'dir') )
+ mkdir(fullOutDir);
+ end
+
+ flist = dir(fullfile(dirIn,subfolders{s},'*.tif'));
+ parseStruct = LLSM.ParseFileNames(flist,'*.tif');
+
+ imDT = MicroscopeData.ReadMetadataFile([fullOutDir filesep]);
+ if ( ~isempty(imDT) )
+ fprintf('Found metadata file for %s... Skipping\n', fullOutDir);
+ continue;
+ end
+
+ %% TODO: Make this a little more generalized
+ chanMap = unique(vertcat(parseStruct.chan));
+ camMap = unique({parseStruct.cam}.');
+
+ frames = unique(vertcat(parseStruct.stack));
+ if ( length(frames) == 1 && ~isempty(parseStruct.iter) )
+ frames = unique(parseStruct.iter);
+ end
+
+ tiffName = LLSM.MatchFilename(parseStruct, camMap(1), frames(1), chanMap(1));
+ chkIm = LLSM.loadtiff(fullfile(dirIn,subfolders{s},[tiffName{1} '.tif']));
+
+ imSize = size(chkIm);
+
+ datasetName = parseStruct(1).datasetName;
+ if ( ~isempty(forceDatasetName) )
+ datasetName = forceDatasetName;
+ end
+
+ %%
+ imD = MicroscopeData.GetEmptyMetadata();
+ imD.DatasetName = datasetName;
+
+ imD.Dimensions = imSize([2,1,3]);
+ imD.NumberOfChannels = length(camMap);
+ imD.NumberOfFrames = length(frames);
+ imD.PixelFormat = class(chkIm);
+
+ imD.ChannelColors = [1,1,1; 0,1,0];
+ imD.ChannelNames = {'FM 4-64','Calcein'};
+ imD.PixelPhysicalSize = [0.104, 0.104, firstMetaSettings.zOffset];
+ imD.StartCaptureDate = firstMetaSettings.startCaptureDate;
+
+ %% Some KLB Helpers for getting block sizing
+ bytes = MicroscopeData.Helper.GetBytesPerVoxel(chkIm);
+
+ blockSize_xyzct = [64, 64, imSize(3), 1, 1];
+ blockMem = prod(blockSize_xyzct)*bytes;
+ blockSize_xyzct(5) = max(1,floor((1024^2) / blockMem)); % try to get close to 1MB block size
+
+ myCluster = parcluster('local');
+ threads = myCluster.NumWorkers;
+
+ %%
+% poolobj = gcp('nocreate');
+% if isempty(poolobj)
+% parpool(2);
+% end
+
+ prg = Utils.CmdlnProgress(length(frames), true, sprintf('Reading %s', imD.DatasetName));
+ for t=1:size(frames,1)
+ frameIm = zeros([imSize, length(camMap)]);
+
+ % Try to skip completed output files
+ fileList = getKLBFiles(imD.DatasetName,fullOutDir,1:length(camMap),[t,t],true,false);
+ bFilesExist = cellfun(@(x)(exist(x,'file') > 0), fileList);
+ if ( all(bFilesExist) )
+ continue;
+ end
+
+ bReadErr = false;
+ for c=1:length(camMap)
+ tiffName = LLSM.MatchFilename(parseStruct, camMap(c), frames(t), chanMap(1));
+ try
+ frameIm(:,:,:,c) = LLSM.loadtiff(fullfile(dirIn,subfolders{s},[tiffName{1} '.tif']));
+ catch me
+ bReadErr = true;
+ warning(['ERROR reading ''' tiffName{1} '.tif' ''' - skipping: ' me.message]);
+ break;
+ end
+ end
+ if ( bReadErr )
+ % Skip frame-write if there was an error reading one of the tiffs
+ continue;
+ end
+
+ fileNameT = sprintf('%s_t%04d.klb', fullfile(fullOutDir,imD.DatasetName), t);
+
+ wrtFunc = @()(MicroscopeData.KLB.writeKLBstack(frameIm, fileNameT, threads, [0.104 0.104, firstMetaSettings.zOffset, 1,1], blockSize_xyzct));
+ safeWrite(fileMap, fileList, wrtFunc);
+
+ prg.PrintProgress(t);
+ end
+
+ MicroscopeData.CreateMetadata(fullOutDir, imD);
+ end
+end
+
+function fileList = getKLBFiles(datasetName,outDir,chanList,timeRange,filePerT,filePerC)
+ fileList = {};
+
+ fileName = fullfile(outDir,datasetName);
+ if (filePerT)
+ if (filePerC)
+ for t=1:length(timeRange(1):timeRange(2))
+ for c=1:length(chanList)
+ fileList = [fileList; {sprintf('%s_c%d_t%04d.klb',fileName,chanList(c),(t-1)+timeRange(1))}];
+ end
+ end
+ else
+ for t=1:length(timeRange(1):timeRange(2))
+ fileList = [fileList; {sprintf('%s_t%04d.klb',fileName,(t-1)+timeRange(1))}];
+ end
+ end
+ elseif (filePerC)
+ for c=1:length(chanList)
+ fileList = [fileList; {sprintf('%s_c%d.klb',fileName,chanList(c))}];
+ end
+ else
+ fileList = {[fileName,'.klb']};
+ end
+end
+
+
+function safeWrite(fileMap, fileList, writeFunc)
+ for i=1:length(fileList)
+ fileMap(fileList{i}) = false;
+ end
+
+ writeFunc();
+
+ for i=1:length(fileList)
+ remove(fileMap,fileList{i});
+ end
+end
+
+function cleanupFunc(fileMap)
+ k = keys(fileMap);
+
+ for i=1:length(k)
+ if ( exist(k{i}, 'file') )
+ delete(k{i});
+ end
+ end
+end
+
diff --git a/+LLSM/FixupMetadata.m b/+LLSM/FixupMetadata.m
new file mode 100644
index 0000000000000000000000000000000000000000..817fa093051e0cb68a38aa01674a710cce8d020f
--- /dev/null
+++ b/+LLSM/FixupMetadata.m
@@ -0,0 +1,38 @@
+function FixupMetadata(metadataPath, rawNames)
+ [imD, jsonDir] = MicroscopeData.ReadMetadataFile(metadataPath);
+ if ( isempty(imD) )
+ return
+ end
+
+ parseStruct = LLSM.ParseFileNames(rawNames);
+
+ % TODO: This is very specific to older LLSM data names
+ chanMap = unique(vertcat(parseStruct.chan));
+ camMap = unique(vertcat(parseStruct.cam));
+
+ numFrames = nnz(vertcat(parseStruct.cam)==camMap(1) & vertcat(parseStruct.chan)==chanMap(1));
+
+ msecDeltas = zeros(numFrames, length(camMap));
+ msecAbs = zeros(numFrames, length(camMap));
+ for c=1:length(camMap)
+ bSelect = (vertcat(parseStruct.cam)==camMap(c) & vertcat(parseStruct.chan)==chanMap(1));
+ msecDeltas(:,c) = vertcat(parseStruct(bSelect).msec);
+ msecAbs(:,c) = vertcat(parseStruct(bSelect).msecAbs);
+ end
+
+ % TODO: Won't work for actual iterations datasets
+ if ( nnz(msecDeltas(:,1)) < size(msecDeltas,1)-1 )
+ msecDeltas = msecAbs - msecAbs(1,:);
+ end
+
+ % TODO: Need to handle timestamps better
+ tsDeltas = permute(msecDeltas / 1000, [3,2,1]);
+ imD.StartCaptureDate = datestr(imD.StartCaptureDate, 'yyyy-mm-dd HH:MM:SS');
+ imD.TimeStampDelta = tsDeltas(:,:,1:imD.NumberOfFrames);
+
+ if ( imD.NumberOfFrames ~= size(tsDeltas,3) )
+ warning(sprintf('%s: Mismatched frames/timestamps truncated timestamps (%d ~= %d)', imD.DatasetName, imD.NumberOfFrames, size(tsDeltas,3)));
+ end
+
+ MicroscopeData.CreateMetadata(jsonDir, imD);
+end
diff --git a/+LLSM/FixupSmadarData.m b/+LLSM/FixupSmadarData.m
new file mode 100644
index 0000000000000000000000000000000000000000..2bf358f1ec7f955491d48674e0c2d37c78d7eda3
--- /dev/null
+++ b/+LLSM/FixupSmadarData.m
@@ -0,0 +1,42 @@
+function FixupSmadarData(rootDir, infoDir)
+ fList = dir(infoDir);
+ bValidDir = arrayfun(@(x)(~strncmp(x.name, '.',1)&&x.isdir>0), fList);
+
+ expList = fList(bValidDir);
+ for i=1:length(expList)
+ infoFiles = dir(fullfile(infoDir, expList(i).name, '*.txt'));
+ for j=1:length(infoFiles)
+ tok = regexp(infoFiles(j).name, '(.*)_flist\.txt', 'tokens','once');
+ if ( isempty(tok) )
+ continue;
+ end
+ datasetName = [expList(i).name '_' tok{1}];
+
+ if ( ~exist(fullfile(rootDir, datasetName), 'dir') )
+ continue;
+ end
+
+ updateSubdirs(fullfile(rootDir, datasetName), datasetName, fullfile(infoDir, expList(i).name, infoFiles(j).name));
+ end
+ end
+end
+
+function updateSubdirs(dataDir, datasetName, infoFile)
+ fprintf('Updating %s ...\n', dataDir);
+ orgFilenames = importdata(infoFile);
+
+ fList = dir(dataDir);
+ bValidDir = arrayfun(@(x)(~strncmp(x.name, '.',1)&&x.isdir>0), fList);
+
+ dataSubdirs = fList(bValidDir);
+ for i=1:length(dataSubdirs)
+ imDir = [fullfile(dataDir,dataSubdirs(i).name) filesep];
+ LLSM.FixupMetadata(imDir, orgFilenames);
+% LLSM.RenameDataset(imDir, datasetName);
+%
+% if ( endsWith(dataSubdirs(i).name,'KLB') )
+% dropKLB = dataSubdirs(i).name(1:end-3);
+% movefile(fullfile(dataDir,dataSubdirs(i).name),fullfile(dataDir,dropKLB));
+% end
+ end
+end
diff --git a/+LLSM/MatchFilename.m b/+LLSM/MatchFilename.m
new file mode 100644
index 0000000000000000000000000000000000000000..d7a73835605ac4c4ba22d32050d119b6df3247a4
--- /dev/null
+++ b/+LLSM/MatchFilename.m
@@ -0,0 +1,70 @@
+function outNames = MatchFilename(parseStruct, camera,frame,channel)
+ outNames = {};
+
+ bValidStacks = checkNonemptyField(parseStruct, 'stack');
+ if (bValidStacks)
+ useStacks = true;
+ numStacks = max(vertcat(parseStruct.stack))+1;
+ else
+ useStacks = false;
+ numStacks = [];
+ end
+
+ bValidIter = checkNonemptyField(parseStruct, 'iter');
+ if ( bValidIter )
+ useIter = true;
+ else
+ useIter = false;
+ end
+
+ bValidCam = checkNonemptyField(parseStruct, 'cam');
+ if ( bValidCam )
+ unqCams = unique({parseStruct.cam}.');
+ if (length(unqCams)>1)
+ useCams = true;
+ else
+ if ( ~isempty(camera) && all(unqCams~=camera) )
+ return
+ end
+ useCams = false;
+ end
+ else
+ if ( ~isempty(camera) )
+ return
+ end
+ useCams = false;
+ end
+
+ if ( useIter )
+ if ( useStacks )
+ iterIdx = floor((frame) / numStacks);
+ stackIdx = frame - iterIdx*numStacks;
+ timeMask = (vertcat(parseStruct.stack) == stackIdx) & (vertcat(parseStruct.iter) == iterIdx);
+ else
+ timeMask = (vertcat(parseStruct.iter) == frame);
+ end
+ else
+ timeMask = (vertcat(parseStruct.stack) == frame);
+ end
+
+ name = {};
+ if ( ~isempty(camera) && useCams )
+ camChanMask = ( strcmpi(camera,{parseStruct.cam}.') & (vertcat(parseStruct.chan)==channel) );
+ if (any(camChanMask))
+ name = {parseStruct(timeMask & camChanMask).name}.';
+ end
+ elseif (~isempty(channel))
+ chanMask = (vertcat(parseStruct.chans) == channel);
+ if (any(chanMask))
+ name = {parseStruct(timeMask & chanMask).name}.';
+ end
+ else
+ name ={parseStruct(timeMask).name}.';
+ end
+
+ outNames = name;
+end
+
+function bNonempty = checkNonemptyField(parseStruct, fieldname)
+ bNonempty = all(arrayfun(@(x)(~isempty(x.(fieldname))), parseStruct));
+end
diff --git a/+LLSM/RenameDataset.m b/+LLSM/RenameDataset.m
new file mode 100644
index 0000000000000000000000000000000000000000..f2ee420a541f979430d4d19dd6cfc5658f7b2a0c
--- /dev/null
+++ b/+LLSM/RenameDataset.m
@@ -0,0 +1,28 @@
+function RenameDataset(imPath, datasetName)
+ [imD, jsonDir, jsonFile] = MicroscopeData.ReadMetadataFile(imPath);
+ if ( isempty(imD) )
+ return
+ end
+
+ oldDatasetName = imD.DatasetName;
+
+ oldFile = fullfile(jsonDir,jsonFile);
+ delete(oldFile);
+
+ imD.DatasetName = datasetName;
+ MicroscopeData.CreateMetadata(imPath, imD);
+
+ flist = dir(fullfile(imPath, [oldDatasetName '*.klb']));
+
+ inNames = {flist.name}.';
+ guardName = regexptranslate('escape', oldDatasetName);
+ outNames = regexprep(inNames, [guardName '(.*)'], [datasetName '$1']);
+
+ for i=1:length(outNames)
+ if ( strcmp(inNames{i},outNames{i}) )
+ continue;
+ end
+% fprintf('%s --> %s\n', inNames{i}, outNames{i});
+ movefile(fullfile(imPath,inNames{i}), fullfile(imPath,outNames{i}));
+ end
+end
diff --git a/+Movie/movie_enable_uitoolbars.m b/+Movie/movie_enable_uitoolbars.m
new file mode 100644
index 0000000000000000000000000000000000000000..5064e432dffc9cd65352579d41f7a11641188426
--- /dev/null
+++ b/+Movie/movie_enable_uitoolbars.m
@@ -0,0 +1,20 @@
+function movie_enable_uitoolbars(bEnable, uiserver_info)
+ if ( ~exist('uiserver_info','var') )
+ web_opts = weboptions('Timeout',60);
+ uiserver_info = struct('host','http://localhost:4444', ...
+ 'opts',web_opts);
+ end
+ uiserver_info.opts.MediaType = 'application/json';
+
+ ui = webread([uiserver_info.host '/ui'],uiserver_info.opts);
+
+ toolbarDisplayMode = 'none';
+ if ( bEnable )
+ toolbarDisplayMode = 'block';
+ end
+
+ ui.sidebar = toolbarDisplayMode;
+ ui.webToolbar = toolbarDisplayMode;
+
+ webwrite([uiserver_info.host '/ui'],uiserver_info.opts, ui);
+end
\ No newline at end of file
diff --git a/+Movie/movie_hist_match.m b/+Movie/movie_hist_match.m
new file mode 100644
index 0000000000000000000000000000000000000000..5dbbd87dcea90eb60485383c4552254a9ca87b0b
--- /dev/null
+++ b/+Movie/movie_hist_match.m
@@ -0,0 +1,332 @@
+function movie_hist_match(movieName, matchFrame, startFrame, endFrame, match_type)
+ dataDir = 'Data';
+
+ if ( ~any(strcmpi(match_type, {'ks','mdpa'})) )
+ error('Invalid matching type');
+ end
+
+ %% Setup movie frame-cap dir
+ frameDir = fullfile(dataDir, [movieName '_frames']);
+ if ( ~exist(frameDir, 'dir') )
+ mkdir(frameDir);
+ end
+
+ %% Setup server info and options to pass around
+ web_opts = weboptions('Timeout',60, 'MediaType','application/json');
+ uiserver_info = struct('host','http://localhost:4444', ...
+ 'opts',web_opts);
+
+ %% Load constant info
+ CONSTANTS = webread([uiserver_info.host '/CONSTANTS'], uiserver_info.opts);
+ imD = CONSTANTS.imageData;
+ imD.Dimensions = imD.Dimensions';
+
+ %% Put frame start/end and matching info in the histogram params file (for reference)
+ movieInfo = struct('startFrame',{startFrame}, 'endFrame',{endFrame}, 'matchFrame',{matchFrame});
+
+ %% Assume current transfer function is target (and get params)
+ renderparams = webread([uiserver_info.host '/renderparams'], uiserver_info.opts);
+ viewparams = webread([uiserver_info.host '/view'], uiserver_info.opts);
+
+ mvparams = struct('render',{renderparams}, 'view',{viewparams});
+ mvparams.movie = movieInfo;
+
+ %% Load match-frame and make transformed histogram
+ im = MicroscopeData.Reader('imageData',imD, 'timeRange',[matchFrame,matchFrame]);
+% [rawHist,centers] = computeRawHist(im);
+% matchhsts = computeTfmHistIm(im, mvparams.render);
+
+ [rhst,hstc] = computeRawHist(im);
+ thsts = transformHist(rhst,hstc, mvparams.render);
+% thsts = transformHistGauss(rhst,hstc, mvparams.render);
+
+ %% Load all other frames and try to match histogram
+ mvparams.frametfm = cell(endFrame,1);
+ prg = Utils.CmdlnProgress(endFrame-startFrame+1,[],'Matching Histogram');
+ for t=startFrame:endFrame
+ if ( t == matchFrame )
+ continue;
+ end
+
+ im = MicroscopeData.Reader('imageData',imD, 'timeRange',[t,t]);
+ matchparams = search_params_match(im, thsts, mvparams.render, match_type);
+ mvparams.frametfm{t} = matchparams;
+
+ prg.PrintProgress(t-startFrame+1);
+ end
+
+ %% Save matched params
+ fid = fopen(fullfile(frameDir, [movieName '_histp_' match_type '.json']), 'wt');
+ fprintf(fid, '%s\n', Utils.CreateJSON(mvparams, true));
+ fclose(fid);
+end
+
+function matchparams = search_params_match(im, matchhsts, initparams, match_type)
+ matchparams = [];
+
+ if ( strcmpi(match_type,'ks') )
+ cmp_func = @ksdist;
+ elseif ( strcmpi(match_type,'mdpa') )
+ cmp_func = @mdpa;
+ end
+
+ [rhst,hstc] = computeRawHist(im);
+
+ stepsz = 0.05;
+ xmin = (0:stepsz:(1.0-stepsz)).';
+ xmax = (stepsz:stepsz:1.0).';
+ ymid = (stepsz:stepsz:(1.0-stepsz)).';
+ for c=1:size(matchhsts,1)
+% ymid = clamp(ymid, 0.01,0.99);
+
+ [xxmin,yymid,xxmax] = meshgrid(xmin,ymid,xmax);
+
+ xxmin = xxmin(:);
+ xxmax = xxmax(:);
+
+ bValid = (xxmax-xxmin) >= stepsz;
+ xxmin = xxmin(bValid);
+ xxmax = xxmax(bValid);
+
+ xxmid = (xxmin+xxmax)/2;
+
+ yymid = yymid(bValid);
+
+ d = Inf(size(xxmin,1), 1);
+ paramInfo = [];
+ for i=1:size(xxmin,1)
+ A = [xxmin(i)^2 xxmin(i) 1; xxmid(i)^2 xxmid(i) 1; xxmax(i)^2 xxmax(i) 1];
+ p = A \ [0; yymid(i); 1];
+
+ pa = p(1);
+ pb = p(2);
+ pc = p(3);
+
+ pmin = xxmin(i);
+ pmax = xxmax(i);
+
+ if (pa > 0)
+ pmin = max(-pb/(2*pa), pmin);
+ elseif (pa < 0)
+ pmax = min(-pb/(2*pa), pmax);
+ end
+
+ params = struct('a',{pa}, 'b',{pb}, 'c',{pc}, 'minRange',{pmin}, 'maxRange',{pmax});
+ %% TODO - Fix this!!!!!!
+% thst = resampleTfmHistGauss(rhst(c,:), hstc, pa,pb,pc,pmin,pmax);
+ thst = resampleTfmHist(rhst(c,:), hstc, pa,pb,pc,pmin,pmax);
+% thst = computeTfmHistIm(im, params, c);
+ d(i) = cmp_func(thst, matchhsts(c,:));
+% d(i) = mdpa(thst, matchhsts(c,:));
+% d(i) = ksdist(thst, matchhsts(c,:));
+ paramInfo = [paramInfo; params];
+ end
+
+ [~,bestIdx] = min(d);
+ bestParams = paramInfo(bestIdx);
+ matchparams = [matchparams; bestParams];
+ end
+end
+
+function d = ksdist(h1,h2)
+ F1 = cumsum(h1);
+ F2 = cumsum(h2);
+
+ d = max(abs(F2-F1));
+end
+
+% Mimimum distance of pair assignments (Cha and Srihari, Pattern Recognition 2002)
+function d = mdpa(h1,h2)
+ dH = cumsum(h1-h2);
+ d = sum(abs(dH));
+end
+
+function thst = computeTfmHistIm(data, tfmParams, chans)
+ tfmBins = 256;
+ if ( ~exist('chans','var') )
+ chans = 1:size(data,4);
+ end
+
+ tfmIdx = 1:length(chans);
+ if ( length(tfmParams) == 1 )
+ tfmIdx = ones(1,length(chans));
+ end
+
+ thst = zeros(length(chans), tfmBins);
+ for i = 1:length(chans)
+ ci = chans(i);
+ ti = tfmIdx(i);
+
+ tmin = tfmParams(ti).minRange;
+ tmax = tfmParams(ti).maxRange;
+
+ a = tfmParams(ti).a;
+ b = tfmParams(ti).b;
+ c = tfmParams(ti).c;
+
+ d = mat2gray(data(:,:,:,ci));
+
+ tv = clamp(d, tmin,tmax);
+ tvals = clamp(a*tv.^2 + b*tv + c, 0,1);
+
+ thst(i,:) = histcounts(tvals, tfmBins);
+ end
+end
+
+function [hst,c] = computeRawHist(data)
+ startBins = 8192;
+ numCh = size(data,4);
+ hst = zeros(numCh, startBins);
+ for i = 1:size(data,4)
+ d = mat2gray(data(:));
+ hst(i,:) = histcounts(d, startBins);
+ end
+
+ spc = 1 / startBins;
+ c = (0.5*spc):spc:1.0;
+end
+
+function thst = transformHist(inhst, hstc, tfmParams)
+ thst = [];
+ for i=1:size(inhst,1)
+ xmin = tfmParams(i).minRange;
+ xmax = tfmParams(i).maxRange;
+
+ a = tfmParams(i).a;
+ b = tfmParams(i).b;
+ c = tfmParams(i).c;
+
+ chst = resampleTfmHist(inhst(i,:), hstc, a,b,c,xmin,xmax);
+ thst = [thst; chst];
+ end
+end
+
+function thst = transformHistGauss(inhst, hstc, tfmParams)
+ thst = [];
+ for i=1:size(inhst,1)
+ xmin = tfmParams(i).minRange;
+ xmax = tfmParams(i).maxRange;
+
+ a = tfmParams(i).a;
+ b = tfmParams(i).b;
+ c = tfmParams(i).c;
+
+ chst = resampleTfmHistGauss(inhst(i,:), hstc, a,b,c,xmin,xmax);
+ thst = [thst; chst];
+ end
+end
+
+% function thst = computeTfmHist(rawHist, hctrs, tfmParams)
+% tfmBins = 256;
+%
+% tspc = 1 / tfmBins;
+% ledge = (0:(tfmBins-1)) * tspc;
+% redge = (1:tfmBins) * tspc;
+%
+% % HACK - Makes sure 0 or 1 value maps into output histogram space
+% ledge(1) = -1;
+% redge(end) = 2;
+%
+% thst = zeros(size(rawHist,1), tfmBins);
+% for i=1:size(rawHist,1)
+% tmin = tfmParams(i).minRange;
+% tmax = tfmParams(i).maxRange;
+%
+% a = tfmParams(i).a;
+% b = tfmParams(i).b;
+% c = tfmParams(i).c;
+%
+% tc = clamp(hctrs, tmin,tmax);
+% tctrs = clamp(a*tc.^2 + b*tc + c, 0,1);
+%
+% cmt = repmat(tctrs, tfmBins,1);
+% bbins = (cmt > ledge.') & (cmt <= redge.');
+%
+% cnts = repmat(rawHist(i,:), tfmBins,1);
+% n = zeros(size(cnts));
+% n(bbins) = cnts(bbins);
+%
+% thst(i,:) = sum(n,2).';
+% end
+% end
+
+% function hst = resampleTfmHist(inhst, hstc, a,b,c,xmin,xmax)
+% tfmBins = 256;
+%
+% tspc = 1 / tfmBins;
+% ledge = ((0:(tfmBins-1)) * tspc).';
+% redge = ((1:tfmBins) * tspc).';
+%
+% if ( abs(a) < 1e-6 )
+% iledge = (ledge - c) / b;
+% iredge = (redge - c) / b;
+% else
+% iledge = (-b + sqrt(b^2 - 4*a*(c-ledge))) / (2*a);
+% iredge = (-b + sqrt(b^2 - 4*a*(c-redge))) / (2*a);
+% end
+%
+% %% Hack to make far left/right bins eat everything else
+% iledge(1) = -1;
+% iredge(end) = 2;
+%
+% repc = repmat(hstc, tfmBins,1);
+% bbins = (repc > iledge) & (repc <= iredge);
+%
+% reph = repmat(inhst, tfmBins,1);
+% n = zeros(size(reph));
+% n(bbins) = reph(bbins);
+%
+% hst = (sum(n,2).');
+% end
+
+function hst = resampleTfmHist(inhst, hstc, a,b,c,xmin,xmax)
+ tfmBins = 256;
+
+ tspc = 1 / tfmBins;
+ ledge = ((0:(tfmBins-1)) * tspc).';
+ redge = ((1:tfmBins) * tspc).';
+
+ %% Hack to make far left/right bins eat everything else
+ ledge(1) = -1;
+ redge(end) = 2;
+
+ hcc = clamp(hstc, xmin,xmax);
+ tfmc = clamp(a*hcc.^2 + b*hcc + c, 0,1);
+
+ repc = repmat(tfmc, tfmBins,1);
+ bbins = (repc > ledge) & (repc <= redge);
+
+ reph = repmat(inhst, tfmBins,1);
+ n = zeros(size(reph));
+ n(bbins) = reph(bbins);
+
+ hst = (sum(n,2).');
+end
+
+function hst = resampleTfmHistGauss(inhst, hstc, a,b,c,xmin,xmax)
+ tfmBins = 256;
+
+ tspc = 1 / tfmBins;
+ ledge = ((0:(tfmBins-1)) * tspc).';
+ redge = ((1:tfmBins) * tspc).';
+
+ smean = (ledge + redge) / 2;
+ ssig = tspc;
+
+ hcc = clamp(hstc, xmin,xmax);
+ tfmc = clamp(a*hcc.^2 + b*hcc + c, 0,1);
+
+ repc = repmat(tfmc, tfmBins,1);
+ reph = repmat(inhst, tfmBins,1);
+
+ sval = exp(-0.5*(repc-smean).^2 / ssig^2);
+ sw = sum(sval(:)) / tfmBins;
+
+ sval = ((length(inhst) / tfmBins) / sw) * sval;
+
+ hst = (sum(reph .* sval, 2)).';
+end
+
+function cv = clamp(v, minv,maxv)
+ cv = max(min(v, maxv), minv);
+end
diff --git a/+Movie/render_matched_movie.m b/+Movie/render_matched_movie.m
new file mode 100644
index 0000000000000000000000000000000000000000..ecd812c3c13880d4f17af71aa45be155db1f5415
--- /dev/null
+++ b/+Movie/render_matched_movie.m
@@ -0,0 +1,135 @@
+function render_matched_movie(movieName, match_type)
+ dataDir = 'Data';
+
+ %% Check for match_type info json
+ prvFrameDir = fullfile(dataDir, [movieName '_frames']);
+ jsonFile = fullfile(prvFrameDir, [movieName '_histp_' match_type '.json'] );
+ if ( ~exist(prvFrameDir, 'dir') )
+ warning([jsonFile ' json not found']);
+ return;
+ end
+
+ if ( ~exist(jsonFile, 'file') )
+ warning([jsonFile ' json not found']);
+ return;
+ end
+
+
+ %% Setup output frame-cap dir
+ frameDir = fullfile(dataDir, [movieName '_' match_type '_frames']);
+ if ( ~exist(frameDir, 'dir') )
+ mkdir(frameDir);
+ end
+
+ movieInfo = jsondecode(fileread(jsonFile));
+ movieInfo.movieName = movieName;
+ movieInfo.frameDir = frameDir;
+
+ %% Setup server info and options to pass around
+ web_opts = weboptions('Timeout',60, 'MediaType','application/json');
+ uiserver_info = struct('host','http://localhost:4444', 'opts',web_opts);
+
+ %% Render all frames to folder
+ renderframes(movieInfo, uiserver_info);
+
+ Movie.movie_enable_uitoolbars(true, uiserver_info);
+
+% %% Load constant info
+% CONSTANTS = webread([uiserver_info.host '/CONSTANTS'], uiserver_info.opts);
+%
+% %% TODO: FFMPEG encode
+% samplerate = diff(squeeze(CONSTANTS.imageData.TimeStampDelta(1,1,:)));
+% asr = median(samplerate);
+%
+% fps = 30;
+% minrate = 5;
+%
+% framerate = 12 * (6 / (asr));
+% framerate = max(framerate, minrate);
+%
+% inpattern = fullfile(frameDir,[movieName '_t%04d.png']);
+% movieFile = [movieName '_fps' num2str(round(framerate)) '.mp4'];
+%
+% ffmpegStr = sprintf('ffmpeg -y -framerate %4.2f -i %s -pix_fmt yuv420p -r %d %s', framerate, inpattern, fps, movieFile);
+% status = system(ffmpegStr);
+end
+
+function renderframes(movieInfo, uiserver_info)
+ movieName = movieInfo.movieName;
+ frameDir = movieInfo.frameDir;
+
+ renderInfo = movieInfo.render;
+ frametfm = movieInfo.frametfm;
+
+ %% Disable leverjs toolbars and set movie rendering parameters
+ Movie.movie_enable_uitoolbars(false, uiserver_info);
+ pause(0.5);
+
+ %% Try to wait for last initial-render to complete
+ if ( ~awaitRenderDone(300, uiserver_info) )
+ error('Full rendering timed-out after 5 minutes!');
+ end
+
+ pause(0.5);
+
+ baseRender = renderInfo;
+ viewInfo = movieInfo.view;
+
+ maxFrame = length(frametfm);
+ for i=1:maxFrame
+ renderInfo = baseRender;
+
+ if ( ~isempty(frametfm{i}) )
+ for c=1:length(renderInfo)
+ renderInfo(c).a = frametfm{i}(c).a;
+ renderInfo(c).b = frametfm{i}(c).b;
+ renderInfo(c).c = frametfm{i}(c).c;
+ renderInfo(c).minRange = frametfm{i}(c).minRange;
+ renderInfo(c).maxRange = frametfm{i}(c).maxRange;
+ end
+ end
+
+ set_all_params(i, renderInfo, viewInfo, uiserver_info);
+
+ %% Wait for complete volume render
+ if ( ~awaitRenderDone(300, uiserver_info) )
+ error('Full rendering timed-out after 5 minutes!');
+ end
+
+ %% Write frame to directory
+ imFrame = webread([uiserver_info.host '/screenCap'], uiserver_info.opts);
+
+ frameSize = [size(imFrame,1) size(imFrame,2)];
+
+ framePad = mod(frameSize,2);
+ imPad = zeros([frameSize + framePad, 3], 'like',imFrame);
+
+ imPad(1:frameSize(1),1:frameSize(2),:) = imFrame;
+
+ frameFile = fullfile(frameDir, [movieName '_t' num2str(i-1, '%04d') '.png']);
+ imwrite(imPad, frameFile, 'png');
+
+ end
+end
+
+function set_all_params(frame, renderInfo, viewInfo, uiserver_info)
+ webwrite([uiserver_info.host '/renderparams'],uiserver_info.opts, renderInfo);
+ webwrite([uiserver_info.host '/view'],uiserver_info.opts, viewInfo);
+ webwrite([uiserver_info.host '/time/' num2str(frame)] ,'');
+end
+
+function bCompleted = awaitRenderDone(timeout, uiserver_info)
+ bCompleted=false;
+
+ elapsed = 0;
+ chk_timer = tic();
+ while ( ~bCompleted && (timeout <= 0 || elapsed <= timeout) )
+ bCompleted = webread([uiserver_info.host '/drawComplete'], uiserver_info.opts);
+ if ( bCompleted )
+ return;
+ end
+
+ pause(0.5);
+ elapsed = toc(chk_timer);
+ end
+end
diff --git a/+Movie/render_small_movie.m b/+Movie/render_small_movie.m
new file mode 100644
index 0000000000000000000000000000000000000000..33fa6c1c209c2b8bd1c5f45646e404758b2f9939
--- /dev/null
+++ b/+Movie/render_small_movie.m
@@ -0,0 +1,32 @@
+function render_small_movie(movieName, asr, size_target_mb, resx)
+% samplerate = diff(squeeze(CONSTANTS.imageData.TimeStampDelta(1,1,:)));
+% asr = median(samplerate);
+ frameDir = fullfile('Data', [movieName '_frames']);
+ nframes = length(dir(fullfile(frameDir,[movieName '*.png'])));
+
+ if ( nframes <= 0 )
+ error('No video frames found!');
+ end
+
+ fps = 30;
+ minrate = 5;
+
+ framerate = 12 * (6 / (asr));
+ framerate = max(framerate, minrate);
+
+ duration = fps * nframes;
+ bitrate = round(8192 * size_target_mb / duration);
+
+ inpattern = fullfile(frameDir,[movieName '_t%04d.png']);
+ movieFile = fullfile('final_movies', 'small_movies', [movieName '_small_fps' num2str(round(framerate)) '.mp4']);
+
+% ffmpegP1Str = sprintf('ffmpeg -y -framerate %4.2f -i %s -vf scale=%d:-1 -c:v libx264 -b:v %dk -pass 1 -pix_fmt yuv420p -r %d -an -f mp4 NUL', framerate, inpattern, resx, bitrate, fps);
+% status = system(ffmpegP1Str);
+%
+% ffmpegP2Str = sprintf('ffmpeg -y -framerate %4.2f -i %s -vf scale=%d:-1 -c:v libx264 -b:v %dk -pass 2 -pix_fmt yuv420p -r %d %s', framerate, inpattern, resx, bitrate, fps, movieFile);
+% % ffmpegP2Str = sprintf('ffmpeg -y -i %s -c:v libx264 -b:v %dk -pass 2 %s', inpattern, bitrate, movieFile);
+% status = system(ffmpegP2Str);
+
+ ffmpegStr = sprintf('ffmpeg -y -framerate %4.2f -i %s -vf scale=%d:-1 -preset veryslow -crf 28 -pix_fmt yuv420p -r %d %s', framerate, inpattern, resx, fps, movieFile);
+ status = system(ffmpegStr);
+end
diff --git a/+Movie/rendermovie.m b/+Movie/rendermovie.m
new file mode 100644
index 0000000000000000000000000000000000000000..b6de1a6868beb998832581cfc44513d7468e5cbc
--- /dev/null
+++ b/+Movie/rendermovie.m
@@ -0,0 +1,138 @@
+function rendermovie(movieName, startFrame, endFrame)
+ dataDir = 'Data';
+
+ %% Setup movie framce-cap dir
+ frameDir = fullfile(dataDir, [movieName '_frames']);
+ if ( ~exist(frameDir, 'dir') )
+ mkdir(frameDir);
+ end
+
+ %% Setup server info and options to pass around
+ web_opts = weboptions('Timeout',60, 'MediaType','application/json');
+ uiserver_info = struct('host','http://localhost:4444', ...
+ 'opts',web_opts);
+
+
+
+ %% Disable leverjs toolbars and set movie rendering parameters
+ Movie.movie_enable_uitoolbars(false, uiserver_info);
+ Movie.set_movie_params(movieName, uiserver_info);
+
+ pause(0.5);
+
+ %% Try to wait for last initial-render to complete
+ if ( ~awaitRenderDone(300, uiserver_info) )
+ error('Full rendering timed-out after 5 minutes!');
+ end
+
+ pause(0.5);
+
+ %% Load constant info
+ CONSTANTS = webread([uiserver_info.host '/CONSTANTS'], uiserver_info.opts);
+
+ samplerate = diff(squeeze(CONSTANTS.imageData.TimeStampDelta(1,1,:)));
+ asr = median(samplerate);
+
+ %% Put frame start/end info in the movie-params file (for reference)
+ movieInfo = struct('startFrame',{startFrame}, 'endFrame',{endFrame}, 'avgSampleRate',{asr});
+
+ mvParams = jsondecode(fileread(fullfile(dataDir,[movieName '_params.json'])));
+ mvParams.movie = movieInfo;
+
+ fid = fopen(fullfile(dataDir, [movieName '_params.json']), 'wt');
+ fprintf(fid, '%s\n', Utils.CreateJSON(mvParams, true));
+ fclose(fid);
+
+ %% Compare current params to file-params
+% if ( ~cmp_movie_params(movieName) )
+% % Delete all old frames if movie params have changed
+% delete(fullfile(frameDir,'*.*'));
+% end
+ paramfname = [movieName '_params.json'];
+ if ( exist(fullfile(frameDir,paramfname),'file') )
+ movieJSON = strtrim(fileread(fullfile(frameDir,paramfname)));
+ paramJSON = strtrim(fileread(fullfile(dataDir,paramfname)));
+ if ( ~strcmpi(movieJSON,paramJSON) )
+ button = questdlg('','','Continue','Cancel','Cancel');
+ if ( strcmpi(button,'Cancel') )
+ return;
+ else
+ % Delete all old frames if movie params have changed
+ delete(fullfile(frameDir,'*.*'));
+ copyfile(fullfile(dataDir,paramfname), fullfile(frameDir,paramfname));
+ end
+ end
+ else
+ copyfile(fullfile(dataDir,paramfname), fullfile(frameDir,paramfname));
+ end
+
+ %% Render movie frames
+ startFrame = max(1, startFrame);
+ endFrame = min(CONSTANTS.imageData.NumberOfFrames, endFrame);
+
+ for time=startFrame:endFrame
+ % Use frame-index for image output (0-based)
+ frmIdx = time - startFrame;
+ %% Skip frames that have already been written
+ frameFile = fullfile(frameDir, [movieName '_t' num2str(frmIdx, '%04d') '.png']);
+ if ( exist(frameFile,'file') )
+ continue;
+ end
+
+ %% Set frame time
+ webwrite([uiserver_info.host '/time/' num2str(time)] ,'');
+
+ %% Wait for complete volume render
+ if ( ~awaitRenderDone(300, uiserver_info) )
+ error('Full rendering timed-out after 5 minutes!');
+ end
+
+ %% Write frame to directory
+ imFrame = webread([uiserver_info.host '/screenCap'], web_opts);
+
+ frameSize = [size(imFrame,1) size(imFrame,2)];
+
+ framePad = mod(frameSize,2);
+ imPad = zeros([frameSize + framePad, 3], 'like',imFrame);
+
+ imPad(1:frameSize(1),1:frameSize(2),:) = imFrame;
+
+ imwrite(imPad, frameFile, 'png');
+ end
+
+ Movie.movie_enable_uitoolbars(true, uiserver_info);
+
+ %% TODO: FFMPEG encode
+ fps = 30;
+ minrate = 5;
+
+ framerate = 12 * (6 / (asr));
+ framerate = max(framerate, minrate);
+
+ inpattern = fullfile(frameDir,[movieName '_t%04d.png']);
+ movieFile = fullfile('final_movies', [movieName '_fps' num2str(round(framerate)) '.mp4']);
+
+ ffmpegStr = sprintf('ffmpeg -y -framerate %4.2f -i %s -pix_fmt yuv420p -r %d %s', framerate, inpattern, fps, movieFile);
+ status = system(ffmpegStr);
+
+% %% Remove frame dir (if ffmpeg returns no errors)
+% if ( status == 0 )
+% rmdir(frameDir,'s');
+% end
+end
+
+function bCompleted = awaitRenderDone(timeout, uiserver_info)
+ bCompleted=false;
+
+ elapsed = 0;
+ chk_timer = tic();
+ while ( ~bCompleted && (timeout <= 0 || elapsed <= timeout) )
+ bCompleted = webread([uiserver_info.host '/drawComplete'], uiserver_info.opts);
+ if ( bCompleted )
+ return;
+ end
+
+ pause(0.5);
+ elapsed = toc(chk_timer);
+ end
+end
diff --git a/+Movie/save_movie_params.m b/+Movie/save_movie_params.m
new file mode 100644
index 0000000000000000000000000000000000000000..3b6bd690a3fc7562e6b91b1c5827bf3aa6051e95
--- /dev/null
+++ b/+Movie/save_movie_params.m
@@ -0,0 +1,21 @@
+function save_movie_params(movieName, uiserver_info)
+ if ( ~exist('uiserver_info','var') )
+ web_opts = weboptions('Timeout',60);
+ uiserver_info = struct('host','http://localhost:4444', ...
+ 'opts',web_opts);
+ end
+ uiserver_info.opts.MediaType = 'application/json';
+
+
+ renderparams = webread([uiserver_info.host '/renderparams'], uiserver_info.opts);
+ viewparams = webread([uiserver_info.host '/view'], uiserver_info.opts);
+
+ outparams = struct('render',{renderparams}, 'view',{viewparams});
+
+ jsonData = Utils.CreateJSON(outparams, true);
+
+ dataDir = 'Data';
+ fid = fopen(fullfile(dataDir, [movieName '_params.json']), 'wt');
+ fprintf(fid, '%s\n', jsonData);
+ fclose(fid);
+end
diff --git a/+Movie/set_movie_params.m b/+Movie/set_movie_params.m
new file mode 100644
index 0000000000000000000000000000000000000000..7d98b7d73ab74049b4cd5bb70b925c38c0e3e42b
--- /dev/null
+++ b/+Movie/set_movie_params.m
@@ -0,0 +1,21 @@
+function set_movie_params(movieName, uiserver_info)
+ if ( ~exist('uiserver_info','var') )
+ web_opts = weboptions('Timeout',60);
+ uiserver_info = struct('host','http://localhost:4444', ...
+ 'opts',web_opts);
+ end
+ uiserver_info.opts.MediaType = 'application/json';
+
+ dataDir = 'Data';
+ jsonData = fileread(fullfile(dataDir, [movieName '_params.json']));
+ inParams = jsondecode(jsonData);
+
+ if ( ~isfield(inParams.view, 'volColor') )
+ inParams.view.volColor = [0,0,0];
+ end
+
+ inParams.view.bgColor = [0,0,0,1];
+
+ webwrite([uiserver_info.host '/renderparams'],uiserver_info.opts, inParams.render);
+ webwrite([uiserver_info.host '/view'],uiserver_info.opts, inParams.view);
+end
diff --git a/+Segment/add_fake_sphere_seg.m b/+Segment/add_fake_sphere_seg.m
new file mode 100644
index 0000000000000000000000000000000000000000..1a3d6128736853b5fe249ad4a9015ef9871e5ab8
--- /dev/null
+++ b/+Segment/add_fake_sphere_seg.m
@@ -0,0 +1,22 @@
+function add_fake_sphere_seg(leverFile, radius)
+ conn = database(leverFile, '','', 'org.sqlite.JDBC', 'jdbc:sqlite:');
+ CONSTANTS = Read.getConstants(conn);
+
+ imDims = Utils.SwapXY_RC(CONSTANTS.imageData.Dimensions);
+ pxSize = CONSTANTS.imageData.PixelPhysicalSize;
+
+ [f,e,v,n] = Segment.make_fake_sphere(radius, imDims, pxSize);
+
+ [X,Y,Z] = meshgrid(1:imDims(1),1:imDims(2),1:imDims(3));
+ coords = [X(:),Y(:),Z(:)];
+
+ com = imDims / 2;
+ bInside = (sum((coords-com).^2, 2) <= radius^2);
+ pts = uint16(coords(bInside,:));
+
+ segStruct = struct('time',{[1]}, 'channel',{[2]}, 'maxRadius',{[radius]}, 'area',{size(pts,1)}, ...
+ 'centroid',{imDims/2}, 'pts',{coords(bInside)}, 'segCC',{[1]}, ...
+ 'verts', {v}, 'edges',{e}, 'normals',{n}, 'faces',{f});
+
+ Write.CreateCells_3D(conn, segStruct);
+end
diff --git a/+Segment/hipdog.m b/+Segment/hipdog.m
new file mode 100644
index 0000000000000000000000000000000000000000..8aa10222f2d66ac10bd695ad6dc5b3f7ed734ba8
--- /dev/null
+++ b/+Segment/hipdog.m
@@ -0,0 +1,6 @@
+function imDoG = hipdog(imIn, sigma)
+ s2 = sigma * sqrt(2);
+ s1 = sigma / sqrt(2);
+
+ imDoG = HIP.Gaussian(imIn, s2,[],[]) - HIP.Gaussian(imIn, s1,[],[]);
+end
diff --git a/+Segment/leverjs_seg_dataset.m b/+Segment/leverjs_seg_dataset.m
new file mode 100644
index 0000000000000000000000000000000000000000..69015a08d0291e06afb73b154f7f7224042203de
--- /dev/null
+++ b/+Segment/leverjs_seg_dataset.m
@@ -0,0 +1,38 @@
+function leverjs_seg_dataset(leverFile, thresh)
+ conn = database(leverFile, '','', 'org.sqlite.JDBC', 'jdbc:sqlite:');
+ CONSTANTS = Read.getConstants(conn);
+
+ maxFrame = CONSTANTS.imageData.NumberOfFrames;
+ chkDeltas = squeeze(CONSTANTS.imageData.TimeStampDelta(1,1,:));
+
+ medDelta = median(diff(chkDeltas,1));
+ if ( medDelta < 7.0 )
+ % Cut high-framerates at 20-minutes
+ maxTime = 20*60;
+ dropFrame = find(chkDeltas >= maxTime ,1,'first');
+ if ( ~isempty(dropFrame) )
+ maxFrame = dropFrame;
+ end
+ end
+
+
+ prg = Utils.CmdlnProgress(maxFrame, true, ['Segmenting: ' CONSTANTS.imageData.DatasetName]);
+ for t=1:maxFrame
+ segs = Segment.segFrame(t, CONSTANTS, thresh, 0.17);
+ if ( isempty(segs) )
+ continue;
+ end
+
+ Write.CreateCells_3D(conn, segs);
+
+ prg.PrintProgress(t);
+ end
+
+ prg = Utils.CmdlnProgress(maxFrame, true, ['Computing distances: ' CONSTANTS.imageData.DatasetName]);
+ for t=1:maxFrame
+ Distance.cellDistance(conn,CONSTANTS, t, [(t-2):(t-1)]);
+
+ prg.PrintProgress(t);
+ end
+ fprintf('Finished segmenting: %s\n\n', CONSTANTS.imageData.DatasetName);
+end
diff --git a/+Segment/llsm_mask.m b/+Segment/llsm_mask.m
new file mode 100644
index 0000000000000000000000000000000000000000..286a32a2844a4279a634d365f1cc3547c02cc29b
--- /dev/null
+++ b/+Segment/llsm_mask.m
@@ -0,0 +1,10 @@
+function bwMask = llsm_mask(im)
+ %% Try to find the "parallelogram" mask of LLSM data
+ bwLLSM = any(im > 0, 4);
+
+ % Fill holes in each slice to get the full mask
+ bwMask = false(size(bwLLSM));
+ for z=1:size(bwLLSM,3)
+ bwMask(:,:,z) = imfill(bwLLSM(:,:,z), 'holes');
+ end
+end
diff --git a/+Segment/logSeg.m b/+Segment/logSeg.m
new file mode 100644
index 0000000000000000000000000000000000000000..67f7100124acd8ab300daa95d03888b368830ada
--- /dev/null
+++ b/+Segment/logSeg.m
@@ -0,0 +1,30 @@
+function [bwIm,imNg] = logSeg(im, imD, logThresh, vesicle_rad_um)
+ %%
+% logThresh = 0.003;
+ scale_factor = 1.0;
+ vesicle_rad_px = vesicle_rad_um ./ imD.PixelPhysicalSize;
+
+ imCalcium = single(mat2gray(im));
+% szFilter = round(size(imCalcium)/3);
+% imLP = HIP.Gaussian(imCalcium, szFilter, 1, []);
+% imHP = imCalcium - imLP;
+% % radKern = HIP.MakeEllipsoidMask([2 2 1]);
+% % imx = HIP.MedianFilter(imHP, radKern, 1, []);
+% imx = imHP;
+
+ %%
+% imDog = hipdog(imCalcium, 2*scale_factor*vesicle_rad_px.*[1 1 1]);
+ imDog = HIP.LoG(imCalcium, 2*scale_factor*vesicle_rad_px.*[1 1 1], []);
+
+ bwMask = false(size(imDog));
+ bwMask(imDog<0) = true;
+
+ imNg = zeros(size(imDog));
+ imNg(bwMask) = -imDog(bwMask);
+
+ % TODO: Figure out how to pick a better threshold
+ bw = imNg >= logThresh;
+
+ se = HIP.MakeEllipsoidMask([2,2,1]);
+ bwIm = HIP.Opener(bw, se, [], []);
+end
diff --git a/+Segment/makePoly.m b/+Segment/makePoly.m
new file mode 100644
index 0000000000000000000000000000000000000000..5bc0d70580e7ad92b37922ed269ac2a73da69dc4
--- /dev/null
+++ b/+Segment/makePoly.m
@@ -0,0 +1,64 @@
+function [faces,verts,normals,edges] = makePoly(pixels_xyz, imSize_rc, maxScale)
+ faces = [];
+ verts = [];
+ normals = [];
+ edges = [];
+
+ if ( ~exist('maxScale','var') )
+ maxScale = 8;
+ end
+
+ PADDING = 2 * maxScale;
+
+ scale_rc = max(imSize_rc) ./ imSize_rc / maxScale;
+ scale_rc = max(min(scale_rc, 0.5), 0.01);
+
+ startCoords_rcz = Utils.SwapXY_RC(min(round(pixels_xyz),[],1));
+ endCoords_rcz = Utils.SwapXY_RC(max(round(pixels_xyz),[],1));
+
+ startPadded = startCoords_rcz-PADDING;
+ endPadded = endCoords_rcz+PADDING;
+
+ roiIm = zeros(ceil(endPadded-startPadded+1));
+
+ paddedPixelList_rcz = Utils.SwapXY_RC(pixels_xyz)-startPadded+1;
+ indList = Utils.CoordToInd(size(roiIm),paddedPixelList_rcz);
+ roiIm(indList) = 1;
+
+
+ if ( all(scale_rc == 1) )
+ smImSize = size(roiIm);
+ imR = roiIm;
+ else
+ smImSize = round(scale_rc .* size(roiIm));
+ imR = imresize3(roiIm, smImSize);
+ end
+
+ [f,v] = isosurface(imR, graythresh(imR));
+ if ( size(v,1) < 10 )
+ smImSize = round(2 * scale_rc .* size(roiIm));
+ imR = imresize3(roiIm, smImSize);
+ [f,v] = isosurface(imR, graythresh(imR));
+ end
+
+ if ( isempty(v) )
+ return;
+ end
+
+ % rescale verts back to
+ unscale_rc = (size(roiIm) ./ size(imR));
+ % put unscale on x,y instead of r,c
+ unscale_xy = Utils.SwapXY_RC(unscale_rc);
+
+ verts = (v-0.5) .* unscale_xy + Utils.SwapXY_RC(startPadded);
+ normals = D3d.Polygon.CalcNorms(verts, f);
+
+ % make sure face/edge indices are 0-based
+ faces = f-1;
+ edges = makeEdges(f)-1;
+end
+
+function edges = makeEdges(faces)
+ fullEdges = [faces(:,[1 2]); faces(:,[2 3]); faces(:,[3 1])];
+ edges = unique(sort(fullEdges, 2), 'rows','stable');
+end
diff --git a/+Segment/makeSegStruct.m b/+Segment/makeSegStruct.m
new file mode 100644
index 0000000000000000000000000000000000000000..74f994a2f5c8ea0cec55b812a917f1ea47edbf6a
--- /dev/null
+++ b/+Segment/makeSegStruct.m
@@ -0,0 +1,46 @@
+function segs = makeSegStruct(time,chan, rp, imSize, isoScale, orgL)
+ if ( ~exist('orgL','var') )
+ orgL = [];
+ end
+
+ % Pre-allocate segs structure for frame
+ times = num2cell(repmat(time, length(rp),1));
+ chans = num2cell(repmat(chan, length(rp),1));
+ segs = struct('time',times, 'channel',chans, 'maxRadius',{[]}, 'area',{[]}, ...
+ 'centroid',{[]}, 'pts',{[]}, 'segCC',{[]}, ...
+ 'verts', {[]}, 'edges',{[]}, 'normals',{[]}, 'faces',{[]});
+
+ for i=1:length(rp)
+ ns = segs(i);
+
+ ns.area = size(rp(i).PixelList,1);
+ ns.centroid = rp(i).Centroid;
+ ns.pts = uint16(rp(i).PixelList);
+
+ if ( isempty(orgL) )
+ ns.segCC = i;
+ else
+ ns.segCC = findOrgCC(rp(i).PixelIdxList, orgL);
+ end
+
+ [ns.faces,ns.verts,ns.normals,ns.edges] = Segment.makePoly(rp(i).PixelList, imSize, isoScale);
+ if ( isempty(ns.faces) )
+ continue;
+ end
+
+ ns.maxRadius = sqrt(max(sum((ns.verts - ns.centroid).^2, 2)));
+
+ segs(i) = ns;
+ end
+
+ %% Drop empty segmentations (e.g. failed isosurfacing)
+ bValidIso = arrayfun(@(x)(~isempty(x.faces)), segs);
+ segs = segs(bValidIso);
+end
+
+function orgCC = findOrgCC(idxList, orgL)
+ orgCC = unique(orgL(idxList));
+ if ( length(orgCC) > 1 )
+ orgCC = orgCC(find(orgCC,1,'first'));
+ end
+end
diff --git a/+Segment/make_fake_sphere.m b/+Segment/make_fake_sphere.m
new file mode 100644
index 0000000000000000000000000000000000000000..569641235c0af86957a0fe83571b05df428d3f52
--- /dev/null
+++ b/+Segment/make_fake_sphere.m
@@ -0,0 +1,141 @@
+function [f,e,v,n] = make_fake_sphere(radius, imDims, pxSize)
+ [f,e,v,n] = make_sphere(radius);
+
+ scl = repmat(min(pxSize) ./ pxSize, size(v,1),1);
+ v = v .* scl;
+
+ n = calcNormals(v,f);
+% n = n ./ scl;
+% n = n ./ sqrt(sum(n.^2,2));
+
+% figure;hold on;
+% trisurf(f, v(:,1),v(:,2),v(:,3), 'EdgeColor','none');
+%
+% ppX = [v(e(:,1),1) v(e(:,2),1)].';
+% ppY = [v(e(:,1),2) v(e(:,2),2)].';
+% ppZ = [v(e(:,1),3) v(e(:,2),3)].';
+% plot3(ppX, ppY, ppZ, '-k');
+%
+% npX = [v(:,1) v(:,1)+0.2*n(:,1)].';
+% npY = [v(:,2) v(:,2)+0.2*n(:,2)].';
+% npZ = [v(:,3) v(:,3)+0.2*n(:,3)].';
+% plot3(npX,npY,npZ, '-r');
+%
+% v = v + (imDims ./ 2);
+end
+
+function [f,e,v,n] = make_sphere(r)
+ nsegs = 20;
+ theta = linspace(-180,180, 2*nsegs);
+ phi = linspace(-90,90, nsegs);
+
+ [ttheta,tphi] = meshgrid(theta(1:end-1), phi(2:end-1));
+
+ uy = sind(tphi);
+ ux = cosd(tphi) .* cosd(ttheta);
+ uz = cosd(tphi) .* sind(ttheta);
+
+ v = r * [0,-1,0;
+ ux(:),uy(:),uz(:);
+ 0,1,0];
+
+ n = v ./ sqrt(sum(v.^2, 2));
+
+ f = [];
+ e = [];
+
+ vq_offset = 1;
+ %% Make traingulated quad faces (and quad edges)
+ for i=1:size(uy,1)-1
+ for j=1:size(uy,2)
+ jnxt = mod(j, size(uy,2)) + 1;
+
+ bl = sub2ind(size(uy), i,j) + vq_offset;
+ br = sub2ind(size(uy), i,jnxt) + vq_offset;
+ tl = sub2ind(size(uy), i+1,j) + vq_offset;
+ tr = sub2ind(size(uy), i+1,jnxt) + vq_offset;
+
+ qt1 = [bl, br, tl];
+ qt2 = [br, tr, tl];
+
+ qel = [tl, bl];
+ qeb = [bl, br];
+
+ f = [f; qt1;qt2];
+ e = [e; qel;qeb];
+ end
+ end
+
+ %% Make bottom triangle cap (and edges)
+ bcf = [];
+ bce = [];
+ for j=1:size(uy,2)
+ jnxt = mod(j, size(uy,2)) + 1;
+
+ tl = sub2ind(size(uy), 1,j) + vq_offset;
+ tr = sub2ind(size(uy), 1,jnxt) + vq_offset;
+
+ bcf = [bcf; 1,tr,tl];
+ bce = [bce; tl,1];
+ end
+
+ %% Make top triangle cap (and edges)
+ tcf = [];
+ tce = [];
+ for j=1:size(uy,2)
+ jnxt = mod(j, size(uy,2)) + 1;
+
+ bl = sub2ind(size(uy), size(uy,1),j) + vq_offset;
+ br = sub2ind(size(uy), size(uy,1),jnxt) + vq_offset;
+
+ tel = [size(v,1), bl];
+ teb = [bl, br];
+
+ tcf = [tcf; size(v,1),bl,br];
+ tce = [tce; tel;teb];
+ end
+
+ f = [bcf;f;tcf];
+ e = [bce;e;tce];
+
+% figure;hold on;
+% trisurf(f, v(:,1),v(:,2),v(:,3), 'EdgeColor','none');
+%
+% ppX = [v(e(:,1),1) v(e(:,2),1)].';
+% ppY = [v(e(:,1),2) v(e(:,2),2)].';
+% ppZ = [v(e(:,1),3) v(e(:,2),3)].';
+% plot3(ppX, ppY, ppZ, '-k');
+%
+% npX = [v(:,1) v(:,1)+0.2*n(:,1)].';
+% npY = [v(:,2) v(:,2)+0.2*n(:,2)].';
+% npZ = [v(:,3) v(:,3)+0.2*n(:,3)].';
+% plot3(npX,npY,npZ, '-r');
+end
+
+function vertNormals = calcNormals(vertList, faceList)
+ vertNormals = zeros(size(vertList,1),3);
+ for i=1:size(faceList,1)
+ face = faceList(i,:);
+ edges = [vertList(face(2),:)-vertList(face(1),:);
+ vertList(face(3),:)-vertList(face(1),:);
+ vertList(face(3),:)-vertList(face(2),:)];
+
+ edgeLen = sqrt(sum(edges.^2,2));
+ edges = edges ./ repmat(edgeLen,1,3);
+
+ faceDir = [dot(edges(1,:),edges(2,:));
+ dot(edges(3,:),-edges(1,:));
+ dot(edges(2,:),edges(3,:))];
+
+ faceAngles = acos(faceDir);
+
+ faceVec = cross(edges(1,:), edges(2,:));
+ faceNorm = faceVec / norm(faceVec);
+
+ for k=1:3
+ vertNormals(face(k),:) = vertNormals(face(k),:) + faceAngles(k)*faceNorm;
+ end
+ end
+ vertLengths = sqrt(sum(vertNormals.^2,2));
+ vertNormals = -(vertNormals ./ repmat(vertLengths,1,3));
+end
diff --git a/+Segment/segFrame.m b/+Segment/segFrame.m
new file mode 100644
index 0000000000000000000000000000000000000000..42d98e062eb6c67c715edddcaf5a3a45997ef4e0
--- /dev/null
+++ b/+Segment/segFrame.m
@@ -0,0 +1,86 @@
+function [segs,imNg] = segFrame(t, CONSTANTS, thresh, rad_um)
+ [im,imD] = MicroscopeData.Reader('imageData',CONSTANTS.imageData, 'timeRange',[t,t]);
+
+ validVolume = 20;
+
+ fmSegs = [];
+% %% TODO: Use FM 4-64 to try and identify ectoderm as a whole
+% imFM = mat2gray(im(:,:,:,1));
+% thEct = 0.8*graythresh(imFM);
+%
+% bwFM = (imFM >= thEct);
+%
+% closeKern = HIP.MakeEllipsoidMask([3,3,1]);
+% bwEct = HIP.Closure(bwFM, closeKern, [],[]);
+%
+% rp = regionprops(bwEct, 'Area', 'PixelIdxList', 'PixelList', 'Centroid');
+% [mxArea,maxIdx] = max([rp.Area]);
+%
+% if (mxArea < 9000)
+% warning(['Couldn''t find ectoderm in frame: ' t]);
+% else
+% fmSegs = makeSegStruct(t, 1, rp, size(bwEct), 1);
+% end
+
+
+ %% Guess spicule brightness from graythresh
+ % TODO: Try to autoscale threshold somehow?
+
+ %% TODO: Deal with bright spicule first
+ [bw,imNg] = Segment.logSeg(im(:,:,:,2), imD, thresh, rad_um);
+
+ %% Only keep large CCs
+ rp = regionprops(bw, 'Area', 'PixelIdxList');
+ bValid = arrayfun(@(x)(x.Area > validVolume), rp);
+
+ idxList = vertcat(rp(bValid).PixelIdxList);
+
+ bwValid = false(size(bw));
+ bwValid(idxList) = true;
+
+ %% Get initial labeling
+ orgL = bwlabeln(bwValid);
+
+ %% Apply a small erosion to try and separate nearby components
+ axScale = min(CONSTANTS.imageData.PixelPhysicalSize) ./ CONSTANTS.imageData.PixelPhysicalSize;
+
+ erodeKern = HIP.MakeEllipsoidMask([3,3,1]);
+ bwEr = HIP.MinFilter(bwValid, erodeKern, [],[]);
+
+ %% Add back disappeared labels
+ allL = unique(orgL(orgL>0));
+ keptL = orgL(bwEr);
+
+ missingL = setdiff(allL, keptL);
+ bAddBack = ismember(orgL(:), missingL);
+
+ bwEr(bAddBack) = true;
+
+ %% Use distance tfm to backfill
+ L = bwlabeln(bwEr);
+
+ [~, idxFg] = bwdist(bwEr);
+ L(bwValid) = L(idxFg(bwValid));
+
+
+% %% Use Watershed to backfill
+% D = bwdist(bwEr);
+% D(~bw) = -Inf;
+%
+% L = watershed(D);
+% bwReassign = (bw & (L <= 1));
+%
+% % Assign the border/bg pixels that watershed drops
+% [~,idxFg] = bwdist(L>1);
+% L(bwReassign) = idxFg(bwReassign);
+%
+% % Reorder the labels
+% L(L<=1) = 1;
+% L = L-1;
+
+ %% Build seg structure
+ rp = regionprops(L, 'PixelIdxList', 'PixelList', 'Centroid');
+ caSegs = Segment.makeSegStruct(t, 2, rp, size(orgL), 1, orgL);
+
+ segs = [fmSegs; caSegs];
+end
diff --git a/+Stats/abbreviateName.m b/+Stats/abbreviateName.m
new file mode 100644
index 0000000000000000000000000000000000000000..e9955a4beaaadc66723bd3f3c986cc8524876f76
--- /dev/null
+++ b/+Stats/abbreviateName.m
@@ -0,0 +1,22 @@
+function abrv = abbreviateName(name)
+ n_date = name(1:5);
+ type = regexpi(name, '(?:DMSO)|(?:Axtinib)|(?:Axitinib)|(?:BB94)', 'match','once');
+ tlid = regexpi(name, 'TimeLapse(\d+)', 'tokens','once');
+ posid = regexpi(name, 'Pos(\d+)', 'tokens','once');
+
+ % Handle spelling difference
+ type = strrep(type, 'Axitinib','Axtinib');
+
+ wash = regexpi(name, '_wash_', 'match','once');
+ if ( ~isempty(wash) )
+ type = [type '_w'];
+ end
+
+ abrv = [type '_' n_date];
+ if ( ~isempty(tlid) )
+ abrv = [abrv '_T' tlid{:}];
+ end
+ if ( ~isempty(posid) )
+ abrv = [abrv '_P' posid{:}];
+ end
+end
diff --git a/+Stats/check_bb94_stats.m b/+Stats/check_bb94_stats.m
new file mode 100644
index 0000000000000000000000000000000000000000..ccc6c6adc5528ff6a58c0db771ffda5d5989b320
--- /dev/null
+++ b/+Stats/check_bb94_stats.m
@@ -0,0 +1,12 @@
+%%
+Stats.export_bb94_stats
+
+%% Create tables for easier data access
+size_tbl = readtable('ecto_size_stats_bb94.csv');
+
+%% Create csv with comparisons
+fid = fopen('mmpl_sigs.csv', 'wt');
+
+%% Write size comparison
+sum_tbl = Stats.write_summary(fid, size_tbl, 'size', 'Vesicle Size', {'Control';'Metalloproteinase Inhibition'});
+Stats.write_compare(fid, size_tbl, sum_tbl, 'size', 'Vesicle Size Comparison', {'Control';'Metalloproteinase Inhibition'});
diff --git a/+Stats/check_ecto_stats.m b/+Stats/check_ecto_stats.m
new file mode 100644
index 0000000000000000000000000000000000000000..3124e56da361ffdb8b67ac73ff108a4f57d6c255
--- /dev/null
+++ b/+Stats/check_ecto_stats.m
@@ -0,0 +1,57 @@
+%%
+Stats.export_ecto_stats;
+Stats.export_spicule_stats;
+
+%% Create tables for easier data access
+size_tbl = readtable('ecto_size_stats_axtinib.csv');
+vesicle_tbl = readtable('ecto_track_stats_axtinib.csv');
+spc_tbl = readtable('spicule_stats.csv');
+
+%% Create csv with comparisons
+fid = fopen('ecto_sigs.csv', 'wt');
+
+%% Write size comparison
+sum_tbl = Stats.write_summary(fid, size_tbl, 'size', 'Vesicle Size');
+Stats.write_compare(fid, size_tbl, sum_tbl, 'size', 'Vesicle Size Comparison');
+
+%% Speed comparison
+sum_tbl = Stats.write_summary(fid, vesicle_tbl, 'avg_speed', 'Vesicle Speed');
+Stats.write_compare(fid, vesicle_tbl, sum_tbl, 'avg_speed', 'Vesicle Speed Comparison');
+
+%% Directionality comparison
+sum_tbl = Stats.write_summary(fid, vesicle_tbl, 'dir_index', 'Vesicle Directionality');
+Stats.write_compare(fid, vesicle_tbl, sum_tbl, 'dir_index', 'Vesicle Directionality Comparison');
+
+%% Diffusion comparison
+sum_tbl = Stats.write_summary(fid, vesicle_tbl, 'msd_D', 'Vesicle Diffusion');
+Stats.write_compare(fid, vesicle_tbl, sum_tbl, 'msd_D', 'Vesicle Diffusion Comparison');
+
+%% Spicule motion stats
+sum_tbl = Stats.write_summary_spicule(fid, spc_tbl, 'abs_vel', 'Instantaneous Speed');
+[pKW_s,pM_s] = Stats.write_compare_spicule(fid, spc_tbl, sum_tbl, 'abs_vel', 'Instantaneous Speed Comparison');
+
+%% Spicule vel stats
+Stats.write_summary_spicule(fid, spc_tbl, 'inst_vel', 'Instantaneous Velocity');
+[pKW_v,pM_v] = Stats.write_compare_spicule(fid, spc_tbl, sum_tbl, 'inst_vel', 'Velocity Toward Spicul Comparison');
+
+%% Plot pM as a "matrix"
+cmaplh = hot(10);
+cmap = flip(cmaplh(1:9,:),1);
+clims = log10([min(pM_s(:)) max(pM_s(:))]);
+
+figure;imagesc(log10(pM_s),[-9,0]);colormap(cmap);
+hold on;
+for i=1:10
+ plot([i+0.5,i+0.5],[0.5,10.5], '-k');
+ plot([0.5,10.5],[i+0.5,i+0.5], '-k');
+end
+labels = flip(['1', arrayfun(@(x)(['10^{' num2str(x) '}']), (-1:-1:-8), 'UniformOutput',false)]);
+colorbar('Ticks',flip(0:-1:-8), 'TickLabels',labels)
+set(gca,'XAxisLocation','top')
+axis tight
+xlim([0.5,10.5]);
+ylim([0.5,10.5]);
+
+
+%%
+fclose(fid);
diff --git a/+Stats/collate_all_tracks.m b/+Stats/collate_all_tracks.m
new file mode 100644
index 0000000000000000000000000000000000000000..26f8e74a1b1935e53c8876aa25167d931208eb94
--- /dev/null
+++ b/+Stats/collate_all_tracks.m
@@ -0,0 +1,36 @@
+ROOT = 'E:\Smadar\LeverJS';
+
+partial_prefix = 'dataset_collate_tracks_p';
+
+flist = dir(fullfile(ROOT,'*.LEVER'));
+parfor i=1:length(flist)
+ track_data = Stats.collate_tracks(ROOT, flist(i).name);
+ if ( isempty(track_data) )
+ continue;
+ end
+
+ Stats.parsave(partial_prefix, i, {'track_data'}, {track_data});
+end
+
+% flist = dir(fullfile(ROOT,'*.LEVER'));
+% for i=1:length(flist)
+% track_data = Stats.collate_tracks(ROOT, flist(i).name);
+% if ( isempty(track_data) )
+% continue;
+% end
+%
+% Stats.parsave(partial_prefix, i, {'track_data'}, {track_data});
+% end
+
+date_suffix = ['-' datestr(now(),'yyyy-mm-dd')];
+
+track_data = [];
+matList = dir([partial_prefix '*.mat']);
+for i=1:length(matList)
+ s = load(matList(i).name);
+
+ track_data = [track_data; s.track_data];
+
+ outname = ['dataset_collate_tracks' date_suffix '.mat'];
+ save(outname, 'track_data');
+end
diff --git a/+Stats/collate_tracks.m b/+Stats/collate_tracks.m
new file mode 100644
index 0000000000000000000000000000000000000000..eba30b8e69b9d478a0d55e19290cc46c1b7f63cf
--- /dev/null
+++ b/+Stats/collate_tracks.m
@@ -0,0 +1,183 @@
+function track_data = collate_tracks(ROOT, filename)
+ leverFile = fullfile(ROOT,filename);
+
+ track_data = [];
+
+ %% Load constants
+ conn = database(leverFile, '','', 'org.sqlite.JDBC', 'jdbc:sqlite:');
+ CONSTANTS = Read.getConstants(conn);
+
+ %% Setup time and pixel resolution info
+ cutoff = 60*15;
+ bKeep = CONSTANTS.imageData.TimeStampDelta(1,1,:) <= cutoff;
+
+ pxSize = CONSTANTS.imageData.PixelPhysicalSize;
+ pxVol = prod(pxSize);
+
+ tsmp = squeeze(CONSTANTS.imageData.TimeStampDelta(1,1,:));
+
+ %% Set number of frames to consider
+ keeplen = nnz(bKeep);
+
+ %% Create abbreviation for dataset
+ datasetName = CONSTANTS.imageData.DatasetName;
+ abrv = Stats.abbreviateName(datasetName);
+
+ %% Estimate parallelogram mask for first-frame image
+ imagePath = CONSTANTS.imageData.imageDir;
+ im = MicroscopeData.Reader(imagePath, 'timeRange',[1,1]);
+
+ parallelMask = Segment.llsm_mask(im);
+ % HACK: Drop messy top/bottom slices in mask
+ parallelMask(:,:,1) = false;
+ parallelMask(:,:,end) = false;
+
+ %% Get ecto/endo regions
+ dataPath = 'Data';
+ maskPath = fullfile(dataPath, 'ectoMasks');
+
+ ectoPath = fullfile(maskPath, [datasetName '_ectoregion.mat']);
+ endoPath = fullfile(maskPath, [datasetName '_endoregion.mat']);
+
+ bHasEcto = exist(ectoPath, 'file');
+ bHasEndo = exist(endoPath, 'file');
+
+ endoInfo = [];
+ bwEndo = false(size(parallelMask));
+
+ %%%%%%%%%%%%%%%%%%%%
+ %% Exit if there's not ectoderm region info
+ size_stats = [];
+ if ( ~bHasEcto )
+ return;
+ end
+
+ %% Skip if movie is >= 7fps
+ medTsmp = median(diff(tsmp));
+ if ( medTsmp >= 7.0 )
+ return;
+ end
+
+ %% Set frame cutoff should be no more than 5 minutes
+ cutoff = 60 * 5;
+ bKeep = CONSTANTS.imageData.TimeStampDelta(1,1,:) <= cutoff;
+ keeplen = nnz(bKeep);
+
+ %% Load edge/track info
+ trackInfoPath = fullfile(dataPath, [datasetName '_trackinfo.mat']);
+ if ( ~exist(trackInfoPath, 'file') )
+ return;
+ end
+
+ load(fullfile(dataPath, [datasetName '_edgeinfo.mat']));
+ load(trackInfoPath);
+
+ %% Load ectoderm and endoderm region info
+ ectoInfo = load(ectoPath);
+ bwEcto = (parallelMask & ectoInfo.regionInfo.bwMask(:,:,:,1));
+
+ if ( bHasEndo )
+ endoInfo = load(endoPath);
+ bwEndo = (parallelMask & endoInfo.regionInfo.bwMask(:,:,:,1));
+ end
+
+ bwSkel = (parallelMask & ~bwEcto & ~bwEndo);
+
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %% Read per-frame vesicle info
+ bEcto = false(0,1);
+ bSkel = false(0,1);
+ bEndo = false(0,1);
+
+ tids = zeros(0,1);
+ coms = zeros(0,3);
+ ctimes = zeros(0,1);
+ csizes = zeros(0,1);
+ cframes = zeros(0,1);
+
+ imSize = [size(im,1), size(im,2), size(im,3)];
+ for t=1:keeplen
+ Cells = Read.getCells3DTime(conn,t);
+ if ( isempty(Cells) )
+ return;
+ end
+
+ cellsIdx = cellfun(@(x)(Utils.CoordToInd(imSize, double(Utils.SwapXY_RC(x)))), {Cells.pts}.', 'UniformOutput',false);
+ cids = vertcat(Cells.cellID);
+
+ %% Assign vesicles to regions by max overlap
+ nEcto = cellfun(@(x)(nnz(bwEcto(x))), cellsIdx);
+ nSkel = cellfun(@(x)(nnz(bwSkel(x))), cellsIdx);
+ nEndo = cellfun(@(x)(nnz(bwEndo(x))), cellsIdx);
+
+ bEcto(cids) = (nEcto >= nSkel & nEcto >= nEndo);
+ bSkel(cids) = (nSkel >= nEndo & ~bEcto(cids));
+ bEndo(cids) = (~bEcto(cids) & ~bSkel(cids));
+
+ tids(cids,:) = vertcat(Cells.trackID);
+ coms(cids,:) = horzcat(Cells.centroid).' .* pxSize;
+ csizes(cids) = vertcat(Cells.area) .* pxVol;
+ ctimes(cids) = tsmp(vertcat(Cells.time));
+ cframes(cids) = vertcat(Cells.time);
+ end
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %% Tracking info
+
+ %% Get tracks to keep
+ minTrackLen = 7;
+
+ tlens = cellfun(@(x)(length(x)), trackInfo);
+ bKeep = (tlens >= minTrackLen);
+
+ lastId = length(csizes);
+ trackInfo = trackInfo(bKeep);
+ for i=1:length(trackInfo)
+ bValid = (trackInfo{i} <= lastId);
+ trackInfo{i} = trackInfo{i}(bValid);
+ end
+
+ bValidInfo = cellfun(@(x)(length(x) >= 5), trackInfo);
+ trackInfo = trackInfo(bValidInfo);
+
+
+ %% Assign tracks to region by majority vote
+ tnEcto = cellfun(@(x)(sum(bEcto(x))), trackInfo);
+ tnSkel = cellfun(@(x)(sum(bSkel(x))), trackInfo);
+ tnEndo = cellfun(@(x)(sum(bEndo(x))), trackInfo);
+
+ btEcto = (tnEcto >= tnSkel & tnEcto >= tnEndo);
+ btSkel = (tnSkel >= tnEndo & ~btEcto);
+ btEndo = (~btEcto & ~btSkel);
+
+ %%
+ %% Load per-track data
+ typeCell = regexpi(abrv,'(\w*?)_','tokens','once');
+ type = typeCell{1};
+
+ regions = cell(length(trackInfo),1);
+ regions(btEcto) = {'ecto'};
+ regions(btSkel) = {'skel'};
+ regions(btEndo) = {'endo'};
+
+ track_data = [];
+ for i=1:length(trackInfo)
+ %% Compute per-track info
+ cids = trackInfo{i};
+
+ tid = tids(cids(1));
+ ttime = ctimes(cids);
+ tsize = csizes(cids);
+ tpos = coms(cids,:);
+
+ tdt = diff(ttime);
+ tdx = diff(tpos,1,1);
+ tmsd = sum(tdx.^2, 2);
+
+ nts = struct('name',{datasetName}, 'abrv',{abrv}, 'type',{type}, 'region',{regions{i}}, 'srate',{medTsmp}, 'tid',{tid}, 'times',{ttime}, 'pos',{tpos}, 'sizes',{tsize}, 'dt',{tdt}, 'dx',{tdx}, 'msd',{tmsd});
+ track_data = [track_data; nts];
+ end
+
+% track_data = track_data(btEcto | btSkel);
+end
diff --git a/+Stats/export_bb94_stats.m b/+Stats/export_bb94_stats.m
new file mode 100644
index 0000000000000000000000000000000000000000..8f617893c4cb1093b6556a02c9d6f5f63a21d310
--- /dev/null
+++ b/+Stats/export_bb94_stats.m
@@ -0,0 +1,154 @@
+%%
+dmso_order = {'23-06-2016_TimeLapse1_DMSO_20hrs_Calcein_FM464';
+ '28-06-2016_TimeLapse1_Pos2_DMSO_18hrs_Calcein_FM464';
+ '28-06-2016_TimeLapse1_Pos1_DMSO_18hrs_Calcein_FM464';
+ '23-06-2016_TimeLapse2_DMSO_20hrs_Calcein_AzepRH';
+ '28-06-2016_TimeLapse2_DMSO_22hrs_Calcein_FM464';
+ '22-06-2016_TimeLapse1_DMSO_24hrs_Calcein_FM464';
+ '22-06-2016_TimeLapse2_DMSO_24hrs_Calcein_FM464';
+ '22-06-2016_TimeLapse3_DMSO_20hrs_Calcein_FM464';
+ '28-06-2016_TimeLapse3_Pos1_DMSO_25hrs_Calcein_FM464';
+ '28-06-2016_TimeLapse3_Pos3_DMSO_25hrs_Calcein_FM464';
+ '28-06-2016_TimeLapse3_Pos4_DMSO_25hrs_Calcein_FM464';
+ '24-06-2016_Timelapse4_DMSO_24hrs_Calcein_FM464';
+ '01-07-2016_TimeLapse3_DMSO_21hrs_Calcein_FM464';
+ '01-07-2016_TimeLapse4_DMSO_21hrs_Calcein_FM464';
+ '01-07-2016_TimeLapse5_DMSO_23hrs_Calcein_FM464'};
+
+axt_order = {'01-07-2016_TimeLapse1_Pos0_Axtinib_150_18hrs_Calcein_FM464';
+ '01-07-2016_TimeLapse1_Pos1_Axtinib_150_18hrs_Calcein_FM464';
+ '01-07-2016_TimeLapse2_Pos2_Axtinib_150_19hrs_Calcein_FM464';
+ '01-07-2016_TimeLapse2_Pos3_Axtinib_150_19hrs_Calcein_FM464';
+ '24-06-2016_TimeLapse1_Axtinib_150_19hrs_Calcein_FM464';
+ '24-06-2016_TimeLapse2_Axtinib_150_20hrs_Calcein_FM464';
+ '24-06-2016_TimeLapse3_Axtinib_150_21hrs_Calcein_FM464';
+ '30-06-2016_TimeLapse1_Pos0_Axtinib_150_19hrs_Calcein_FM464';
+ '30-06-2016_TimeLapse1_Pos1_Axtinib_150_19hrs_Calcein_FM464';
+ '30-06-2016_TimeLapse1_Pos2_Axtinib_150_19hrs_Calcein_FM464';
+ '23-06-2016_Timelapse_3_Axitinib_36hr_Calcein_FM464'};
+
+bb94_order = {'29-06-2016_TimeLapse1_Pos0_BB94_17hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse1_Pos1_BB94_17hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse1_Pos2_BB94_17hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse2_Pos0_BB94_21hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse2_Pos1_BB94_21hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse2_Pos2_BB94_21hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse3_Pos0_BB94_23hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse3_Pos2_BB94_23hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse3_Pos3_BB94_23hrs_Calcein_FM464';
+ '30-06-2016_TimeLapse2_Pos0_BB94_21hrs_Calcein_FM464';
+ '30-06-2016_TimeLapse2_Pos1_BB94_21hrs_Calcein_FM464';
+ '30-06-2016_TimeLapse2_Pos2_BB94_21hrs_Calcein_FM464';
+ '30-06-2016_TimeLapse2_Pos3_BB94_21hrs_Calcein_FM464'};
+
+all_order = [dmso_order; axt_order; bb94_order];
+
+flist = dir('dataset_ecto_stats-*.mat');
+fnames = {flist.name}.';
+
+load(fnames{end});
+
+%% Sort and only keep matched tracks for now
+matchCells = arrayfun(@(x)(find(strcmpi(all_order,x.name))), ecto_size_stats, 'UniformOutput',false);
+matchTracks = arrayfun(@(x)(find(strcmpi(all_order,x.name))), ecto_track_stats, 'UniformOutput',false);
+
+bMatchedS = cellfun(@(x)(~isempty(x)), matchCells);
+bMatchedT = cellfun(@(x)(~isempty(x)), matchTracks);
+
+chk_size_stats = ecto_size_stats(bMatchedS);
+chk_track_stats = ecto_track_stats(bMatchedT);
+
+[~,invMatch] = sort(vertcat(matchCells{bMatchedS}));
+chk_size_stats = chk_size_stats(invMatch);
+
+[~,invMatch] = sort(vertcat(matchTracks{bMatchedT}));
+chk_track_stats = chk_track_stats(invMatch);
+
+%% Drop empty data
+bValidEctoS = arrayfun(@(x)(~isempty(x.ecto.vsizes)), chk_size_stats);
+bValidSkelS = arrayfun(@(x)(~isempty(x.skel.vsizes)), chk_size_stats);
+
+bValidEctoT = arrayfun(@(x)(~isempty(x.ecto.mean_sizes)), chk_track_stats);
+bValidSkelT = arrayfun(@(x)(~isempty(x.skel.mean_sizes)), chk_track_stats);
+
+bError = (bValidEctoS & ~bValidSkelS);
+err_names = arrayfun(@(x)(x.abrv), chk_size_stats(bError), 'UniformOutput',false);
+if ( any(bError) )
+ fprintf('WARNING: Invalid data detected: No skeletogenic segmentations: (%s)\n', strjoin(err_names,','));
+end
+
+chk_size_stats = chk_size_stats(bValidEctoS);
+
+bError = (bValidEctoT & ~bValidSkelT);
+err_names = arrayfun(@(x)(x.abrv), chk_size_stats(bError), 'UniformOutput',false);
+if ( any(bError) )
+ fprintf('WARNING: Invalid data detected: No skeletogenic segmentations: (%s)\n', strjoin(err_names,','));
+end
+
+chk_track_stats = chk_track_stats(bValidEctoT);
+
+%%
+exp_types = arrayfun(@(x)(regexpi(x.abrv,'(\w*?)_','tokens','once')), chk_size_stats, 'UniformOutput',false);
+exp_types = vertcat(exp_types{:});
+
+ecto_size_info = vertcat(chk_size_stats.ecto);
+skel_size_info = vertcat(chk_size_stats.skel);
+endo_size_info = vertcat(chk_size_stats.endo);
+
+exp_types_track = arrayfun(@(x)(regexpi(x.abrv,'(\w*?)_','tokens','once')), chk_track_stats, 'UniformOutput',false);
+exp_types_track = vertcat(exp_types_track{:});
+
+ecto_track_info = vertcat(chk_track_stats.ecto);
+skel_track_info = vertcat(chk_track_stats.skel);
+endo_track_info = vertcat(chk_track_stats.endo);
+
+%%
+outFile = 'ecto_size_stats_bb94.csv';
+fid = fopen(outFile, 'wt');
+
+% heading = {'experiment','exp_short','type','vid',};
+heading = {'experiment','exp_short','type','region','vid','size'};
+
+hstr = strjoin(heading,',');
+fprintf(fid, '%s\n', hstr);
+
+for i=1:length(chk_size_stats)
+ exp_name = chk_size_stats(i).name;
+ exp_short = chk_size_stats(i).abrv;
+
+ type_cell = regexpi(chk_size_stats(i).abrv,'(\w*?)_','tokens','once');
+ exp_type = type_cell{1};
+
+ if ( strcmpi(exp_type,'Axtinib') )
+ continue;
+ end
+
+ exp_type = replace(exp_type, 'DMSO','Control');
+ exp_type = replace(exp_type, 'Axtinib','VEGFR Inhibition');
+ exp_type = replace(exp_type, 'BB94','Metalloproteinase Inhibition');
+
+ vid = 1;
+ for j=1:length(chk_size_stats(i).skel.vsizes)
+ fprintf(fid, '%s,%s,%s,%s,', exp_name, exp_short, exp_type, 'Mesoderm');
+ fprintf(fid, '%d,', vid);
+
+ fprintf(fid, '%f', chk_size_stats(i).skel.vsizes(j));
+
+ fprintf(fid, '\n');
+
+ vid = vid + 1;
+ end
+
+ for j=1:length(chk_size_stats(i).ecto.vsizes)
+ fprintf(fid, '%s,%s,%s,%s,', exp_name, exp_short, exp_type, 'Ectoderm');
+ fprintf(fid, '%d,', vid);
+
+ fprintf(fid, '%f', chk_size_stats(i).ecto.vsizes(j));
+
+ fprintf(fid, '\n');
+
+ vid = vid + 1;
+ end
+end
+
+fclose(fid);
diff --git a/+Stats/export_ecto_stats.m b/+Stats/export_ecto_stats.m
new file mode 100644
index 0000000000000000000000000000000000000000..eab5347bca2b4d87ed189c3853ffdf20e041e22d
--- /dev/null
+++ b/+Stats/export_ecto_stats.m
@@ -0,0 +1,235 @@
+%%
+dmso_order = {'23-06-2016_TimeLapse1_DMSO_20hrs_Calcein_FM464';
+ '28-06-2016_TimeLapse1_Pos2_DMSO_18hrs_Calcein_FM464';
+ '28-06-2016_TimeLapse1_Pos1_DMSO_18hrs_Calcein_FM464';
+ '23-06-2016_TimeLapse2_DMSO_20hrs_Calcein_AzepRH';
+ '28-06-2016_TimeLapse2_DMSO_22hrs_Calcein_FM464';
+ '22-06-2016_TimeLapse1_DMSO_24hrs_Calcein_FM464';
+ '22-06-2016_TimeLapse2_DMSO_24hrs_Calcein_FM464';
+ '22-06-2016_TimeLapse3_DMSO_20hrs_Calcein_FM464';
+ '28-06-2016_TimeLapse3_Pos1_DMSO_25hrs_Calcein_FM464';
+ '28-06-2016_TimeLapse3_Pos3_DMSO_25hrs_Calcein_FM464';
+ '28-06-2016_TimeLapse3_Pos4_DMSO_25hrs_Calcein_FM464';
+ '24-06-2016_Timelapse4_DMSO_24hrs_Calcein_FM464';
+ '01-07-2016_TimeLapse3_DMSO_21hrs_Calcein_FM464';
+ '01-07-2016_TimeLapse4_DMSO_21hrs_Calcein_FM464';
+ '01-07-2016_TimeLapse5_DMSO_23hrs_Calcein_FM464'};
+
+axt_order = {'01-07-2016_TimeLapse1_Pos0_Axtinib_150_18hrs_Calcein_FM464';
+ '01-07-2016_TimeLapse1_Pos1_Axtinib_150_18hrs_Calcein_FM464';
+ '01-07-2016_TimeLapse2_Pos2_Axtinib_150_19hrs_Calcein_FM464';
+ '01-07-2016_TimeLapse2_Pos3_Axtinib_150_19hrs_Calcein_FM464';
+ '24-06-2016_TimeLapse1_Axtinib_150_19hrs_Calcein_FM464';
+ '24-06-2016_TimeLapse2_Axtinib_150_20hrs_Calcein_FM464';
+ '24-06-2016_TimeLapse3_Axtinib_150_21hrs_Calcein_FM464';
+ '30-06-2016_TimeLapse1_Pos0_Axtinib_150_19hrs_Calcein_FM464';
+ '30-06-2016_TimeLapse1_Pos1_Axtinib_150_19hrs_Calcein_FM464';
+ '30-06-2016_TimeLapse1_Pos2_Axtinib_150_19hrs_Calcein_FM464';
+ '23-06-2016_Timelapse_3_Axitinib_36hr_Calcein_FM464'};
+
+bb94_order = {'29-06-2016_TimeLapse1_Pos0_BB94_17hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse1_Pos1_BB94_17hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse1_Pos2_BB94_17hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse2_Pos0_BB94_21hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse2_Pos1_BB94_21hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse2_Pos2_BB94_21hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse3_Pos0_BB94_23hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse3_Pos2_BB94_23hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse3_Pos3_BB94_23hrs_Calcein_FM464';
+ '30-06-2016_TimeLapse2_Pos0_BB94_21hrs_Calcein_FM464';
+ '30-06-2016_TimeLapse2_Pos1_BB94_21hrs_Calcein_FM464';
+ '30-06-2016_TimeLapse2_Pos2_BB94_21hrs_Calcein_FM464';
+ '30-06-2016_TimeLapse2_Pos3_BB94_21hrs_Calcein_FM464'};
+
+all_order = [dmso_order; axt_order; bb94_order];
+
+flist = dir('dataset_ecto_stats-*.mat');
+fnames = {flist.name}.';
+
+load(fnames{end});
+
+%% Sort and only keep matched tracks for now
+matchCells = arrayfun(@(x)(find(strcmpi(all_order,x.name))), ecto_size_stats, 'UniformOutput',false);
+matchTracks = arrayfun(@(x)(find(strcmpi(all_order,x.name))), ecto_track_stats, 'UniformOutput',false);
+
+bMatchedS = cellfun(@(x)(~isempty(x)), matchCells);
+bMatchedT = cellfun(@(x)(~isempty(x)), matchTracks);
+
+chk_size_stats = ecto_size_stats(bMatchedS);
+chk_track_stats = ecto_track_stats(bMatchedT);
+
+[~,invMatch] = sort(vertcat(matchCells{bMatchedS}));
+chk_size_stats = chk_size_stats(invMatch);
+
+[~,invMatch] = sort(vertcat(matchTracks{bMatchedT}));
+chk_track_stats = chk_track_stats(invMatch);
+
+%% Drop empty data
+bValidEctoS = arrayfun(@(x)(~isempty(x.ecto.vsizes)), chk_size_stats);
+bValidSkelS = arrayfun(@(x)(~isempty(x.skel.vsizes)), chk_size_stats);
+
+bValidEctoT = arrayfun(@(x)(~isempty(x.ecto.mean_sizes)), chk_track_stats);
+bValidSkelT = arrayfun(@(x)(~isempty(x.skel.mean_sizes)), chk_track_stats);
+
+bError = (bValidEctoS & ~bValidSkelS);
+err_names = arrayfun(@(x)(x.abrv), chk_size_stats(bError), 'UniformOutput',false);
+if ( any(bError) )
+ fprintf('WARNING: Invalid data detected: No skeletogenic segmentations: (%s)\n', strjoin(err_names,','));
+end
+
+chk_size_stats = chk_size_stats(bValidEctoS);
+
+bError = (bValidEctoT & ~bValidSkelT);
+err_names = arrayfun(@(x)(x.abrv), chk_size_stats(bError), 'UniformOutput',false);
+if ( any(bError) )
+ fprintf('WARNING: Invalid data detected: No skeletogenic segmentations: (%s)\n', strjoin(err_names,','));
+end
+
+chk_track_stats = chk_track_stats(bValidEctoT);
+
+%%
+exp_types = arrayfun(@(x)(regexpi(x.abrv,'(\w*?)_','tokens','once')), chk_size_stats, 'UniformOutput',false);
+exp_types = vertcat(exp_types{:});
+
+ecto_size_info = vertcat(chk_size_stats.ecto);
+skel_size_info = vertcat(chk_size_stats.skel);
+endo_size_info = vertcat(chk_size_stats.endo);
+
+exp_types_track = arrayfun(@(x)(regexpi(x.abrv,'(\w*?)_','tokens','once')), chk_track_stats, 'UniformOutput',false);
+exp_types_track = vertcat(exp_types_track{:});
+
+ecto_track_info = vertcat(chk_track_stats.ecto);
+skel_track_info = vertcat(chk_track_stats.skel);
+endo_track_info = vertcat(chk_track_stats.endo);
+
+%%
+outFile = 'ecto_size_stats_axtinib.csv';
+fid = fopen(outFile, 'wt');
+
+% heading = {'experiment','exp_short','type','vid',};
+heading = {'experiment','exp_short','type','region','vid','size'};
+
+hstr = strjoin(heading,',');
+fprintf(fid, '%s\n', hstr);
+
+for i=1:length(chk_size_stats)
+ exp_name = chk_size_stats(i).name;
+ exp_short = chk_size_stats(i).abrv;
+
+ type_cell = regexpi(chk_size_stats(i).abrv,'(\w*?)_','tokens','once');
+ exp_type = type_cell{1};
+
+ if ( strcmpi(exp_type,'BB94') )
+ continue;
+ end
+
+ exp_type = replace(exp_type, 'DMSO','Control');
+ exp_type = replace(exp_type, 'Axtinib','VEGFR Inhibition');
+ exp_type = replace(exp_type, 'BB94','Metalloproteinase Inhibition');
+
+ vid = 1;
+ for j=1:length(chk_size_stats(i).skel.vsizes)
+ fprintf(fid, '%s,%s,%s,%s,', exp_name, exp_short, exp_type, 'Mesoderm');
+ fprintf(fid, '%d,', vid);
+
+ fprintf(fid, '%f', chk_size_stats(i).skel.vsizes(j));
+
+ fprintf(fid, '\n');
+
+ vid = vid + 1;
+ end
+
+ for j=1:length(chk_size_stats(i).ecto.vsizes)
+ fprintf(fid, '%s,%s,%s,%s,', exp_name, exp_short, exp_type, 'Ectoderm');
+ fprintf(fid, '%d,', vid);
+
+ fprintf(fid, '%f', chk_size_stats(i).ecto.vsizes(j));
+
+ fprintf(fid, '\n');
+
+ vid = vid + 1;
+ end
+end
+
+fclose(fid);
+
+%%
+outFile = 'ecto_track_stats_axtinib.csv';
+fid = fopen(outFile, 'wt');
+
+% heading = {'experiment','exp_short','type','vid',};
+heading = {'experiment','exp_short','type','region','tid',...
+ 'avg_size','avg_speed','path_ratio',...
+ 'win_path_ratio','dir_index','dir_avg_index',...
+ 'msd_D','msd_R',...
+ 'nimsd_D','nimsd_alpha','nimsd_R'};
+
+hstr = strjoin(heading,',');
+fprintf(fid, '%s\n', hstr);
+
+didx_idx = 9;
+
+for i=1:length(chk_track_stats)
+ exp_name = chk_track_stats(i).name;
+ exp_short = chk_track_stats(i).abrv;
+
+ type_cell = regexpi(chk_track_stats(i).abrv,'(\w*?)_','tokens','once');
+ exp_type = type_cell{1};
+
+ exp_type = replace(exp_type, 'DMSO','Control');
+ exp_type = replace(exp_type, 'Axtinib','VEGFR Inhibition');
+ exp_type = replace(exp_type, 'BB94','Metalloproteinase Inhibition');
+
+ tid = 1;
+ for j=1:length(chk_track_stats(i).skel.mean_sizes)
+ fprintf(fid, '%s,%s,%s,%s,', exp_name, exp_short, exp_type, 'Mesoderm');
+ fprintf(fid, '%d,', tid);
+
+ path_rat = chk_track_stats(i).skel.direct_len(j) / chk_track_stats(i).skel.path_len(j);
+ fprintf(fid, '%f,', chk_track_stats(i).skel.mean_sizes(j));
+ fprintf(fid, '%f,', chk_track_stats(i).skel.mean_inst_vel(j));
+ fprintf(fid, '%f,', path_rat);
+
+ fprintf(fid, '%f,', chk_track_stats(i).skel.avg_wpath_ratio{j}(didx_idx));
+ fprintf(fid, '%f,', chk_track_stats(i).skel.avg_wdir_index{j}(didx_idx));
+ fprintf(fid, '%f,', chk_track_stats(i).skel.avg_wdir_mindex{j}(didx_idx));
+
+ fprintf(fid, '%f,', chk_track_stats(i).skel.msd_D(j));
+ fprintf(fid, '%f,', chk_track_stats(i).skel.msd_R(j));
+
+ fprintf(fid, '%f,', chk_track_stats(i).skel.nimsd_D(j));
+ fprintf(fid, '%f,', chk_track_stats(i).skel.nimsd_alpha(j));
+ fprintf(fid, '%f', chk_track_stats(i).skel.nimsd_R(j));
+
+ fprintf(fid, '\n');
+
+ tid = tid + 1;
+ end
+
+ for j=1:length(chk_track_stats(i).ecto.mean_sizes)
+ fprintf(fid, '%s,%s,%s,%s,', exp_name, exp_short, exp_type, 'Ectoderm');
+ fprintf(fid, '%d,', tid);
+
+ path_rat = chk_track_stats(i).ecto.direct_len(j) / chk_track_stats(i).ecto.path_len(j);
+ fprintf(fid, '%f,', chk_track_stats(i).ecto.mean_sizes(j));
+ fprintf(fid, '%f,', chk_track_stats(i).ecto.mean_inst_vel(j));
+ fprintf(fid, '%f,', path_rat);
+
+ fprintf(fid, '%f,', chk_track_stats(i).ecto.avg_wpath_ratio{j}(didx_idx));
+ fprintf(fid, '%f,', chk_track_stats(i).ecto.avg_wdir_index{j}(didx_idx));
+ fprintf(fid, '%f,', chk_track_stats(i).ecto.avg_wdir_mindex{j}(didx_idx));
+
+ fprintf(fid, '%f,', chk_track_stats(i).ecto.msd_D(j));
+ fprintf(fid, '%f,', chk_track_stats(i).ecto.msd_R(j));
+
+ fprintf(fid, '%f,', chk_track_stats(i).ecto.nimsd_D(j));
+ fprintf(fid, '%f,', chk_track_stats(i).ecto.nimsd_alpha(j));
+ fprintf(fid, '%f', chk_track_stats(i).ecto.nimsd_R(j));
+
+ fprintf(fid, '\n');
+
+ tid = tid + 1;
+ end
+end
+
+fclose(fid);
diff --git a/+Stats/export_spicule_stats.m b/+Stats/export_spicule_stats.m
new file mode 100644
index 0000000000000000000000000000000000000000..bbfd6322068fe316e845d8a64ab0e94578052937
--- /dev/null
+++ b/+Stats/export_spicule_stats.m
@@ -0,0 +1,51 @@
+load spicule_stats.mat
+
+dists = [];
+abs_vel = [];
+
+abs_dist = [];
+inst_vel = [];
+for i=1:length(spicule_stats)
+ dists = [dists; spicule_stats(i).skel.sp_mean_dist];
+ abs_vel = [abs_vel; spicule_stats(i).skel.mean_inst_vel];
+
+ abs_dist = [abs_dist; spicule_stats(i).skel.sp_absidelta];
+ inst_vel = [inst_vel; -spicule_stats(i).skel.sp_idelta];
+end
+
+ledges = 0:9;
+redges = 1:10;
+
+dist_grps = zeros(length(dists),1);
+for i=1:length(dists)
+ bIdx = dists(i) > ledges & dists(i) <= redges;
+ if ( ~any(bIdx) )
+ continue;
+ end
+
+ dist_grps(i) = redges(bIdx);
+end
+
+bDropIdx = (dist_grps == 0);
+
+dist_grps = dist_grps(~bDropIdx);
+
+abs_vel_grp = abs_vel(~bDropIdx);
+abs_dist_grp = abs_dist(~bDropIdx);
+inst_vel_grp = inst_vel(~bDropIdx);
+
+%%
+outFile = 'spicule_stats.csv';
+fid = fopen(outFile, 'wt');
+
+heading = {'sp_dist','abs_vel','abs_dist','inst_vel'};
+
+hstr = strjoin(heading,',');
+fprintf(fid, '%s\n', hstr);
+
+for i=1:length(dist_grps)
+ fprintf(fid, '%f,%f,%f,%f\n', dist_grps(i), abs_vel_grp(i), abs_dist_grp(i), inst_vel_grp(i));
+end
+
+fclose(fid);
+
diff --git a/+Stats/load_all_ectostats.m b/+Stats/load_all_ectostats.m
new file mode 100644
index 0000000000000000000000000000000000000000..f8604d3c4b8d384f52e34439e387e38a5a36862f
--- /dev/null
+++ b/+Stats/load_all_ectostats.m
@@ -0,0 +1,39 @@
+ROOT = 'E:\Smadar\LeverJS';
+
+ecto_size_stats = [];
+ecto_track_stats = [];
+
+partial_prefix = 'dataset_ecto_p';
+
+flist = dir(fullfile(ROOT,'*.LEVER'));
+parfor i=1:length(flist)
+ [ecto_sizes, ecto_tracks] = Stats.load_ectostats(ROOT, flist(i).name);
+ if ( isempty(ecto_sizes) )
+ continue;
+ end
+
+ Stats.parsave(partial_prefix, i, {'ecto_sizes','ecto_tracks'}, {ecto_sizes, ecto_tracks});
+end
+
+% flist = dir(fullfile(ROOT,'*.LEVER'));
+% for i=1:length(flist)
+% [ecto_sizes, ecto_tracks] = Stats.load_ectostats(ROOT, flist(i).name);
+% if ( isempty(ecto_sizes) )
+% continue;
+% end
+%
+% Stats.parsave('dataset_ecto_p', i, {'ecto_sizes','ecto_tracks'}, {ecto_sizes, ecto_tracks});
+% end
+
+date_suffix = ['-' datestr(now(),'yyyy-mm-dd')];
+
+matList = dir([partial_prefix '*.mat']);
+for i=1:length(matList)
+ s = load(matList(i).name);
+
+ ecto_size_stats = [ecto_size_stats; s.ecto_sizes];
+ ecto_track_stats = [ecto_track_stats; s.ecto_tracks];
+
+ outname = ['dataset_ecto_stats' date_suffix '.mat'];
+ save(outname, 'ecto_size_stats','ecto_track_stats');
+end
diff --git a/+Stats/load_all_stats.m b/+Stats/load_all_stats.m
new file mode 100644
index 0000000000000000000000000000000000000000..40cc6275ade423f2bab9305671adcb02362fb72e
--- /dev/null
+++ b/+Stats/load_all_stats.m
@@ -0,0 +1,22 @@
+ROOT = 'E:\Smadar\LeverJS';
+
+% if ( ~exist('dataset_stats.mat','file') )
+ size_stats = [];
+ track_stats = [];
+% else
+% load dataset_stats.mat
+% end
+
+flist = dir(fullfile(ROOT,'*.LEVER'));
+for i=1:length(flist)
+% if strncmpi(flist(i).name,'27',2)
+% break;
+% end
+
+ [ds_sizes,ds_tracks] = load_stats(ROOT, flist(i).name);
+
+ size_stats = [size_stats; ds_sizes];
+ track_stats = [track_stats; ds_tracks];
+
+ save('dataset_stats.mat', 'size_stats','track_stats');
+end
diff --git a/+Stats/load_ectostats.m b/+Stats/load_ectostats.m
new file mode 100644
index 0000000000000000000000000000000000000000..9b22f6d78181dd8e3a7fcc6864f8b2c6fc52a067
--- /dev/null
+++ b/+Stats/load_ectostats.m
@@ -0,0 +1,420 @@
+function [size_stats, track_stats] = load_ectostats(ROOT, filename)
+ leverFile = fullfile(ROOT,filename);
+
+ size_stats = [];
+ track_stats = [];
+
+ %% Load constants
+ conn = database(leverFile, '','', 'org.sqlite.JDBC', 'jdbc:sqlite:');
+ CONSTANTS = Read.getConstants(conn);
+
+ %% Setup time and pixel resolution info
+ cutoff = 60*15;
+ bKeep = CONSTANTS.imageData.TimeStampDelta(1,1,:) <= cutoff;
+
+ pxSize = CONSTANTS.imageData.PixelPhysicalSize;
+ pxVol = prod(pxSize);
+
+ tsmp = squeeze(CONSTANTS.imageData.TimeStampDelta(1,1,:));
+
+ %% Set number of frames to consider
+ keeplen = nnz(bKeep);
+
+ %% Create abbreviation for dataset
+ datasetName = CONSTANTS.imageData.DatasetName;
+ abrv = Stats.abbreviateName(datasetName);
+
+ %% Estimate parallelogram mask for first-frame image
+ imagePath = CONSTANTS.imageData.imageDir;
+ im = MicroscopeData.Reader(imagePath, 'timeRange',[1,1]);
+
+ parallelMask = Segment.llsm_mask(im);
+ % HACK: Drop messy top/bottom slices in mask
+ parallelMask(:,:,1) = false;
+ parallelMask(:,:,end) = false;
+
+ %% Get ecto/endo regions
+ dataPath = 'Data';
+ maskPath = fullfile(dataPath, 'ectoMasks');
+
+ ectoPath = fullfile(maskPath, [datasetName '_ectoregion.mat']);
+ endoPath = fullfile(maskPath, [datasetName '_endoregion.mat']);
+
+ bHasEcto = exist(ectoPath, 'file');
+ bHasEndo = exist(endoPath, 'file');
+
+ endoInfo = [];
+ bwEndo = false(size(parallelMask));
+
+ size_stats = [];
+ if ( ~bHasEcto )
+ return;
+ end
+
+ %% Set frame cutoff should be no more than 5 minutes
+ cutoff = 60 * 5;
+ bKeep = CONSTANTS.imageData.TimeStampDelta(1,1,:) <= cutoff;
+ keeplen = nnz(bKeep);
+
+ %% Load edge/track info
+ trackInfoPath = fullfile(dataPath, [datasetName '_trackinfo.mat']);
+ if ( ~exist(trackInfoPath, 'file') )
+ return;
+ end
+
+ load(fullfile(dataPath, [datasetName '_edgeinfo.mat']));
+ load(trackInfoPath);
+
+ %% Load spicule if it exists
+ sp_struct = [];
+ if ( exist(fullfile(dataPath, [datasetName '_spicule.json']), 'file') )
+ sp_json = fileread(fullfile(dataPath, [datasetName '_spicule.json']));
+ sp_struct = jsondecode(sp_json);
+ [cs,rs,als,Rs] = calc_spicule_coords(sp_struct, pxSize);
+ end
+
+ %% Load ectoderm and endoderm region info
+ ectoInfo = load(ectoPath);
+ bwEcto = (parallelMask & ectoInfo.regionInfo.bwMask(:,:,:,1));
+
+ if ( bHasEndo )
+ endoInfo = load(endoPath);
+ bwEndo = (parallelMask & endoInfo.regionInfo.bwMask(:,:,:,1));
+ end
+
+ bwSkel = (parallelMask & ~bwEcto & ~bwEndo);
+
+ %% Initialize size stats structures
+ size_stats = struct('name',{datasetName}, 'abrv',{abrv}, 'tsmp',{tsmp}, 'ecto',{[]}, 'skel',{[]}, 'endo',{[]});
+ size_stats.ecto = struct('rsize',{[]}, 'vsizes',{[]}, 'count',{[]}, 'com',{[]});
+ size_stats.skel = struct('rsize',{[]}, 'vsizes',{[]}, 'count',{[]}, 'com',{[]});
+ size_stats.endo = struct('rsize',{[]}, 'vsizes',{[]}, 'count',{[]}, 'com',{[]});
+
+ size_stats.ecto.rsize = nnz(bwEcto) * pxVol;
+ size_stats.skel.rsize = nnz(bwSkel) * pxVol;
+ size_stats.endo.rsize = nnz(bwEndo) * pxVol;
+
+ %% Initialize track stat structures
+ track_stats = struct('name',{datasetName}, 'abrv',{abrv}, 'tsmp',{tsmp}, 'ecto',{[]}, 'skel',{[]}, 'endo',{[]});
+ track_stats.ecto = struct('mean_sizes',{[]}, 'track_lens',{[]}, 'mean_inst_vel',{[]}, 'path_len',{[]}, ...
+ 'brownian_var',{[]}, 'brownian_cov',{[]}, 'direct_len',{[]}, 'start_pos',{[]},...
+ 'ivel_sm',{[]}, 'isize_sm',{[]},...
+ 'msd_D',{[]}, 'msd_sig',{[]}, 'msd_R',{[]},...
+ 'nimsd_D',{[]}, 'nimsd_alpha',{[]}, 'nimsd_R',{[]},...
+ 'wlen',{[]}, 'avg_path_wlen',{[]}, 'avg_max_wdist',{[]}, 'avg_wdir_index',{[]},...
+ 'avg_wdir_mindex',{[]}, 'avg_wpath_ratio',{[]},...
+ 'sp_mean_dist',{[]}, 'sp_idelta',{[]}, 'sp_tdelta',{[]}, 'sp_dstd',{[]}, 'sp_absidelta',{[]});
+ track_stats.skel = struct('mean_sizes',{[]}, 'track_lens',{[]}, 'mean_inst_vel',{[]}, 'path_len',{[]}, ...
+ 'brownian_var',{[]}, 'brownian_cov',{[]}, 'direct_len',{[]}, 'start_pos',{[]},...
+ 'ivel_sm',{[]}, 'isize_sm',{[]},...
+ 'msd_D',{[]}, 'msd_sig',{[]}, 'msd_R',{[]},...
+ 'nimsd_D',{[]}, 'nimsd_alpha',{[]}, 'nimsd_R',{[]},...
+ 'wlen',{[]}, 'avg_path_wlen',{[]}, 'avg_max_wdist',{[]}, 'avg_wdir_index',{[]},...
+ 'avg_wdir_mindex',{[]}, 'avg_wpath_ratio',{[]},...
+ 'sp_mean_dist',{[]}, 'sp_idelta',{[]}, 'sp_tdelta',{[]}, 'sp_dstd',{[]}, 'sp_absidelta',{[]});
+ track_stats.endo = struct('mean_sizes',{[]}, 'track_lens',{[]}, 'mean_inst_vel',{[]}, 'path_len',{[]}, ...
+ 'brownian_var',{[]}, 'brownian_cov',{[]}, 'direct_len',{[]}, 'start_pos',{[]},...
+ 'ivel_sm',{[]}, 'isize_sm',{[]},...
+ 'msd_D',{[]}, 'msd_sig',{[]}, 'msd_R',{[]},...
+ 'nimsd_D',{[]}, 'nimsd_alpha',{[]}, 'nimsd_R',{[]},...
+ 'wlen',{[]}, 'avg_path_wlen',{[]}, 'avg_max_wdist',{[]}, 'avg_wdir_index',{[]},...
+ 'avg_wdir_mindex',{[]}, 'avg_wpath_ratio',{[]},...
+ 'sp_mean_dist',{[]}, 'sp_idelta',{[]}, 'sp_tdelta',{[]}, 'sp_dstd',{[]}, 'sp_absidelta',{[]});
+
+
+ %% Read per-frame vesicle info
+ bEcto = false(0,1);
+ bSkel = false(0,1);
+ bEndo = false(0,1);
+
+ coms = zeros(0,3);
+ ctimes = zeros(0,1);
+ csizes = zeros(0,1);
+ cframes = zeros(0,1);
+
+ imSize = [size(im,1), size(im,2), size(im,3)];
+ for t=1:keeplen
+ Cells = Read.getCells3DTime(conn,t);
+ if ( isempty(Cells) )
+ return;
+ end
+
+ cellsIdx = cellfun(@(x)(Utils.CoordToInd(imSize, double(Utils.SwapXY_RC(x)))), {Cells.pts}.', 'UniformOutput',false);
+ cids = vertcat(Cells.cellID);
+
+ %% Assign vesicles to regions by max overlap
+ nEcto = cellfun(@(x)(nnz(bwEcto(x))), cellsIdx);
+ nSkel = cellfun(@(x)(nnz(bwSkel(x))), cellsIdx);
+ nEndo = cellfun(@(x)(nnz(bwEndo(x))), cellsIdx);
+
+ bEcto(cids) = (nEcto >= nSkel & nEcto >= nEndo);
+ bSkel(cids) = (nSkel >= nEndo & ~bEcto(cids));
+ bEndo(cids) = (~bEcto(cids) & ~bSkel(cids));
+
+ coms(cids,:) = horzcat(Cells.centroid).' .* pxSize;
+ csizes(cids) = vertcat(Cells.area) .* pxVol;
+ ctimes(cids) = tsmp(vertcat(Cells.time));
+ cframes(cids) = vertcat(Cells.time);
+ end
+
+ bFrame1 = (cframes == 1);
+
+ %% Ecto stats
+ size_stats.ecto.vsizes = csizes(bEcto & bFrame1);
+ size_stats.ecto.count = nnz(bEcto & bFrame1);
+ size_stats.ecto.com = coms(bEcto,:);
+
+ %% Skel stats
+ size_stats.skel.vsizes = csizes(bSkel & bFrame1);
+ size_stats.skel.count = nnz(bSkel & bFrame1);
+ size_stats.skel.com = coms(bSkel,:);
+
+ %% Endo stats
+ size_stats.endo.vsizes = csizes(bEndo & bFrame1);
+ size_stats.endo.count = nnz(bEndo & bFrame1);
+ size_stats.endo.com = coms(bEndo,:);
+
+ %% Tracking stats
+ %% Skip if movie is >= 7fps
+ medTsmp = median(diff(tsmp));
+ if ( medTsmp >= 7.0 )
+ return;
+ end
+
+ %% Get tracks to keep
+ minTrackLen = 7;
+
+ tlens = cellfun(@(x)(length(x)), trackInfo);
+ bKeep = (tlens >= minTrackLen);
+
+ lastId = length(csizes);
+ trackInfo = trackInfo(bKeep);
+ for i=1:length(trackInfo)
+ bValid = (trackInfo{i} <= lastId);
+ trackInfo{i} = trackInfo{i}(bValid);
+ end
+
+ bValidInfo = cellfun(@(x)(length(x) >= 5), trackInfo);
+ trackInfo = trackInfo(bValidInfo);
+
+ %% Assign tracks to region by majority vote
+ tnEcto = cellfun(@(x)(sum(bEcto(x))), trackInfo);
+ tnSkel = cellfun(@(x)(sum(bSkel(x))), trackInfo);
+ tnEndo = cellfun(@(x)(sum(bEndo(x))), trackInfo);
+
+ btEcto = (tnEcto >= tnSkel & tnEcto >= tnEndo);
+ btSkel = (tnSkel >= tnEndo & ~btEcto);
+ btEndo = (~btEcto & ~btSkel);
+
+ % Setup track info to put into groups
+ track_lens = zeros(length(trackInfo),1);
+ track_inst_vel = zeros(length(trackInfo),1);
+ track_mean_size = zeros(length(trackInfo),1);
+ track_path_len = zeros(length(trackInfo),1);
+ track_brn_var = zeros(length(trackInfo),1);
+ track_brn_cov = zeros(3,3,length(trackInfo));
+ track_direct_len = zeros(length(trackInfo),1);
+ track_start_pos = zeros(length(trackInfo),3);
+
+ track_sp_dists = zeros(length(trackInfo),1);
+ track_sp_idelta = zeros(length(trackInfo),1);
+ track_sp_tdelta = zeros(length(trackInfo),1);
+ track_sp_dstd = zeros(length(trackInfo),1);
+ track_sp_iabsdelta = zeros(length(trackInfo),1);
+ track_sp_ivel = zeros(length(trackInfo),1);
+
+ track_wsec = (15:9:105);
+ track_wlen = round(track_wsec ./ medTsmp);
+ track_avg_path_wlen = cell(length(trackInfo),1);
+ track_avg_max_wdisp = cell(length(trackInfo),1);
+ track_avg_wdir_index = cell(length(trackInfo),1);
+ track_avg_wdir_mindex = cell(length(trackInfo),1);
+ track_avg_wpath_rat = cell(length(trackInfo),1);
+
+ %% TODO: Setup and compute MSD fit for each track
+ msd_gamma = 3;
+ track_msd_D = zeros(length(trackInfo),1);
+ track_msd_sig = zeros(length(trackInfo),1);
+ track_msd_R = zeros(length(trackInfo),1);
+
+ track_nimsd_D = zeros(length(trackInfo),1);
+ track_nimsd_alpha = zeros(length(trackInfo),1);
+ track_nimsd_R = zeros(length(trackInfo),1);
+
+ for i=1:length(trackInfo)
+ %% Compute per-track info
+ cids = trackInfo{i};
+
+ ttime = ctimes(cids);
+ tsize = csizes(cids);
+ tpos = coms(cids,:);
+
+ tdt = diff(ttime);
+
+ if ( ~isempty(sp_struct) )
+ T = sp_struct.root .* pxSize;
+ tds = spicule_distances(cs,Rs,rs,als, T,tpos);
+ tds_delta = diff(tds,1);
+
+ track_sp_dists(i) = mean(tds);
+ track_sp_idelta(i) = mean(tds_delta);
+ track_sp_tdelta(i) = tds(end)-tds(1);
+ track_sp_dstd(i) = std(tds);
+ track_sp_iabsdelta(i) = mean(abs(tds_delta));
+ track_sp_ivel(i) = mean(tds_delta ./ tdt);
+ end
+
+ tdx = diff(tpos,1,1);
+ tmsd = cumsum(sum(tdx.^2, 2));
+
+ tmot = sqrt(sum(tdx.^2, 2));
+ tvel = tmot ./ tdt;
+
+ nrmtdx = tdx ./ sqrt(tdt);
+ brntvel = sqrt(sum(nrmtdx.^2, 2));
+% brntvel = tmot ./ sqrt(dt);
+
+ % Compute Diffusion coefficient (with bias)
+ x1 = cumsum(tdt);
+ x2 = ones(size(tdt,1),1);
+ [b,~,~,~,rstat] = regress(tmsd, [x1 x2]);
+
+ track_msd_D(i) = b(1)/(2*msd_gamma);
+ track_msd_sig(i) = b(2)/(2*msd_gamma);
+ track_msd_R(i) = sqrt(rstat(1));
+
+ % Compute non-ideal diffusion coeffient
+ x1 = log(cumsum(tdt));
+ x2 = ones(size(tdt,1),1);
+ [b,~,~,~,rstat] = regress(log(tmsd), [x1 x2]);
+
+ track_nimsd_alpha(i) = b(1);
+ track_nimsd_D(i) = exp(b(2))/(2*msd_gamma);
+ track_nimsd_R(i) = sqrt(rstat(1));
+
+ track_start_pos(i,:) = tpos(1,:);
+ track_lens(i) = length(trackInfo);
+
+ track_inst_vel(i) = mean(tvel);
+ track_mean_size(i) = mean(tsize);
+
+ track_path_len(i) = sum(tmot, 1);
+ track_direct_len(i) = sqrt(sum((tpos(end,:)-tpos(1,:)).^2, 2));
+
+ track_brn_var(i) = var(brntvel);
+ track_brn_cov(:,:,i) = cov(nrmtdx);
+
+ for j=1:length(track_wlen)
+ wlen = track_wlen(j);
+% if ( wlen > length(ttime) )
+% break;
+% end
+
+ nwin = max(length(ttime) - wlen + 1, 1);
+ max_wdisp = zeros(nwin,1);
+ mean_wdisp = zeros(nwin,1);
+ path_wlen = zeros(nwin,1);
+ se_wdisp = zeros(nwin,1);
+ for k=1:nwin
+ idxs = k:min(k+wlen-1, length(ttime));
+ dsps = pdist(tpos(idxs,:));
+
+ path_wlen(k) = sum(tmot(idxs(1:end-1)));
+ max_wdisp(k) = max(dsps);
+ se_wdisp(k) = sqrt(sum((tpos(idxs(end))-tpos(idxs(1))).^2, 2));
+ end
+
+ wdir_index = max_wdisp ./ path_wlen;
+ wdir_mindex = mean_wdisp ./ path_wlen;
+ wpath_rat = se_wdisp ./ path_wlen;
+
+ track_avg_max_wdisp{i} = [track_avg_max_wdisp{i} mean(max_wdisp)];
+ track_avg_path_wlen{i} = [track_avg_path_wlen{i} mean(path_wlen)];
+ track_avg_wdir_index{i} = [track_avg_wdir_index{i} mean(wdir_index)];
+ track_avg_wdir_mindex{i} = [track_avg_wdir_mindex{i} mean(wdir_mindex)];
+ track_avg_wpath_rat{i} = [track_avg_wpath_rat{i} mean(wpath_rat)];
+ end
+ end
+
+ %% Ecto stats
+ track_stats.ecto.mean_sizes = track_mean_size(btEcto);
+ track_stats.ecto.track_lens = track_lens(btEcto);
+ track_stats.ecto.mean_inst_vel = track_inst_vel(btEcto);
+ track_stats.ecto.path_len = track_path_len(btEcto);
+ track_stats.ecto.direct_len = track_direct_len(btEcto);
+ track_stats.ecto.brownian_var = track_brn_var(btEcto);
+ track_stats.ecto.brownian_cov = track_brn_cov(:,:,btEcto);
+ track_stats.ecto.start_pos = track_start_pos(btEcto,:);
+ track_stats.ecto.sp_mean_dist = track_sp_dists(btEcto);
+ track_stats.ecto.sp_idelta = track_sp_ivel(btEcto);
+ track_stats.ecto.sp_tdelta = track_sp_tdelta(btEcto);
+ track_stats.ecto.sp_dstd = track_sp_dstd(btEcto);
+ track_stats.ecto.sp_absidelta = track_sp_iabsdelta(btEcto);
+
+ track_stats.ecto.msd_D = track_msd_D(btEcto);
+ track_stats.ecto.msd_sig = track_msd_sig(btEcto);
+ track_stats.ecto.msd_R = track_msd_R(btEcto);
+ track_stats.ecto.nimsd_D = track_nimsd_D(btEcto);
+ track_stats.ecto.nimsd_alpha = track_nimsd_alpha(btEcto);
+ track_stats.ecto.nimsd_R = track_nimsd_R(btEcto);
+
+% track_stats.ecto.wlen = track_wlen;
+ track_stats.ecto.avg_path_wlen = track_avg_path_wlen(btEcto);
+ track_stats.ecto.avg_max_wdist = track_avg_max_wdisp(btEcto);
+ track_stats.ecto.avg_wdir_index = track_avg_wdir_index(btEcto);
+ track_stats.ecto.avg_wdir_mindex = track_avg_wdir_mindex(btEcto);
+ track_stats.ecto.avg_wpath_ratio = track_avg_wpath_rat(btEcto);
+
+
+ %% Skel stats
+ track_stats.skel.mean_sizes = track_mean_size(btSkel);
+ track_stats.skel.track_lens = track_lens(btSkel);
+ track_stats.skel.mean_inst_vel = track_inst_vel(btSkel);
+ track_stats.skel.path_len = track_path_len(btSkel);
+ track_stats.skel.direct_len = track_direct_len(btSkel);
+ track_stats.skel.brownian_var = track_brn_var(btSkel);
+ track_stats.skel.brownian_cov = track_brn_cov(:,:,btSkel);
+ track_stats.skel.start_pos = track_start_pos(btSkel,:);
+ track_stats.skel.sp_mean_dist = track_sp_dists(btSkel);
+ track_stats.skel.sp_idelta = track_sp_ivel(btSkel);
+ track_stats.skel.sp_tdelta = track_sp_tdelta(btSkel);
+ track_stats.skel.sp_dstd = track_sp_dstd(btSkel);
+ track_stats.skel.sp_absidelta = track_sp_iabsdelta(btSkel);
+
+ track_stats.skel.msd_D = track_msd_D(btSkel);
+ track_stats.skel.msd_sig = track_msd_sig(btSkel);
+ track_stats.skel.msd_R = track_msd_R(btSkel);
+ track_stats.skel.nimsd_D = track_nimsd_D(btSkel);
+ track_stats.skel.nimsd_alpha = track_nimsd_alpha(btSkel);
+ track_stats.skel.nimsd_R = track_nimsd_R(btSkel);
+
+ % track_stats.ecto.wlen = track_wlen;
+ track_stats.skel.avg_path_wlen = track_avg_path_wlen(btSkel);
+ track_stats.skel.avg_max_wdist = track_avg_max_wdisp(btSkel);
+ track_stats.skel.avg_wdir_index = track_avg_wdir_index(btSkel);
+ track_stats.skel.avg_wdir_mindex = track_avg_wdir_mindex(btSkel);
+ track_stats.skel.avg_wpath_ratio = track_avg_wpath_rat(btSkel);
+
+ %% Endo stats
+ track_stats.endo.mean_sizes = track_mean_size(btEndo);
+ track_stats.endo.track_lens = track_lens(btEndo);
+ track_stats.endo.mean_inst_vel = track_inst_vel(btEndo);
+ track_stats.endo.path_len = track_path_len(btEndo);
+ track_stats.endo.direct_len = track_direct_len(btEndo);
+ track_stats.endo.brownian_var = track_brn_var(btEndo);
+ track_stats.endo.brownian_cov = track_brn_cov(:,:,btEndo);
+ track_stats.endo.start_pos = track_start_pos(btEndo,:);
+ track_stats.endo.sp_mean_dist = track_sp_dists(btEndo);
+ track_stats.endo.sp_idelta = track_sp_ivel(btEndo);
+ track_stats.endo.sp_tdelta = track_sp_tdelta(btEndo);
+ track_stats.endo.sp_dstd = track_sp_dstd(btEndo);
+ track_stats.endo.sp_absidelta = track_sp_iabsdelta(btEndo);
+
+ track_stats.endo.msd_D = track_msd_D(btEndo);
+ track_stats.endo.msd_sig = track_msd_sig(btEndo);
+ track_stats.endo.msd_R = track_msd_R(btEndo);
+ track_stats.endo.nimsd_D = track_nimsd_D(btEndo);
+ track_stats.endo.nimsd_alpha = track_nimsd_alpha(btEndo);
+ track_stats.endo.nimsd_R = track_nimsd_R(btEndo);
+
+end
diff --git a/+Stats/load_stats.m b/+Stats/load_stats.m
new file mode 100644
index 0000000000000000000000000000000000000000..0ae244124ed9c4686ea06ede814b2f8c1e9dfc1a
--- /dev/null
+++ b/+Stats/load_stats.m
@@ -0,0 +1,135 @@
+function [size_stats,track_stats] = load_stats(ROOT,filename)
+ leverFile = fullfile(ROOT,filename);
+
+ conn = database(leverFile, '','', 'org.sqlite.JDBC', 'jdbc:sqlite:');
+ CONSTANTS = Read.getConstants(conn);
+
+ cutoff = 60*15;
+ bKeep = CONSTANTS.imageData.TimeStampDelta(1,1,:) <= cutoff;
+
+ pixSize = CONSTANTS.imageData.PixelPhysicalSize;
+ tsmp = squeeze(CONSTANTS.imageData.TimeStampDelta(1,1,:));
+
+ keeplen = nnz(bKeep);
+
+ datasetName = CONSTANTS.imageData.DatasetName;
+ abrv = Stats.abbreviateName(datasetName);
+ size_stats = struct('name',{datasetName}, 'abrv',{abrv}, 'sizes',{[]}, 'cpf',{[]}, 'vpf',{[]}, 'mvpf',{[]});
+ track_stats = struct('name',{datasetName}, 'abrv',{abrv}, 'sizes',{[]}, 'tlen',{[]}, 'ivel',{[]}, ...
+ 'tvel',{[]}, 'slvec',{[]}, 'slvel',{[]}, 'slmot',{[]}, 'slpos',{[]},...
+ 'brnVar',{[]}, 'brnCov',{[]});
+
+ %% Load edge/track info
+ if ( ~exist([datasetName '_trackinfo.mat'], 'file') )
+ return;
+ end
+
+ load([datasetName '_edgeinfo.mat']);
+ load([datasetName '_trackinfo.mat']);
+
+ %% TODO: What is an appropriate minimum track length
+ minTrackLen = 7;
+
+ tlens = cellfun(@(x)(length(x)), trackInfo);
+ bKeep = (tlens >= minTrackLen);
+
+ trackInfo = trackInfo(bKeep);
+
+ numCells = max(cellfun(@(x)(max(x)), trackInfo));
+ coms = zeros(numCells,3);
+ vols = zeros(numCells,1);
+ times = zeros(numCells,1);
+
+ %% Get cell info and size stats
+ for t=1:keeplen
+ Cells = Read.getCells3DTime(conn,t);
+
+ size_stats.cpf = [size_stats.cpf; 0];
+ size_stats.vpf = [size_stats.vpf; 0];
+ size_stats.mvpf = [size_stats.mvpf; 0];
+
+ if ( isempty(Cells) )
+ continue;
+ end
+
+ ids = vertcat(Cells.cellID);
+ cellSizes = vertcat(Cells.area) .* prod(pixSize);
+
+ totalSize = sum(cellSizes);
+ meanSize = mean(cellSizes);
+ countVes = length(cellSizes);
+
+ size_stats.cpf(end) = countVes;
+ size_stats.vpf(end) = totalSize;
+ size_stats.mvpf(end) = meanSize;
+ size_stats.sizes = [size_stats.sizes; cellSizes];
+
+ coms(ids,:) = horzcat(Cells.centroid).' .* pixSize;
+ vols(ids,:) = cellSizes;
+ times(ids,:) = t;
+ end
+
+ %% Track length
+ track_stats.tlen = zeros(length(trackInfo),1);
+ for i=1:length(trackInfo)
+ track_stats.tlen(i) = length(trackInfo{i});
+ end
+
+ %% Average vesicle sizes
+ track_stats.sizes = zeros(length(trackInfo),1);
+ for i=1:length(trackInfo)
+ hids = trackInfo{i};
+ track_stats.sizes(i) = mean(vols(hids));
+ end
+
+ %% Vesicle average instantaneous velocity
+ track_stats.ivel = zeros(length(trackInfo),1);
+ track_stats.tvel = zeros(length(trackInfo),1);
+ for i=1:length(trackInfo)
+ hids = trackInfo{i};
+
+ dx = diff(coms(hids,:), 1);
+ dt = diff(tsmp(times(hids)), 1);
+
+ vel = sqrt(sum(dx.^2, 2)) ./ dt;
+
+ track_stats.ivel(i) = mean(vel);
+ track_stats.tvel(i) = sum(vel);
+ end
+
+ %% Wiener proces variance (\sigma^2)
+ track_stats.brnVar = zeros(length(trackInfo),1);
+ track_stats.brnCov = zeros(3,3, length(trackInfo));
+ for i=1:length(trackInfo)
+ hids = trackInfo{i};
+
+ dx = diff(coms(hids,:), 1);
+ dt = diff(tsmp(times(hids)), 1);
+
+ vel = sum(dx.^2, 2) ./ sqrt(dt);
+ nrmdx = dx ./ repmat(sqrt(dt), 1,3);
+
+ track_stats.brnVar(i) = var(vel);
+ track_stats.brnCov(:,:,i) = cov(nrmdx);
+ end
+
+ %% Vesicle total motion
+ track_stats.slvec = zeros(length(trackInfo),3);
+ track_stats.slvel = zeros(length(trackInfo),1);
+ track_stats.slpos = zeros(length(trackInfo),3);
+ track_stats.slmot = zeros(length(trackInfo),1);
+ for i=1:length(trackInfo)
+ hids = trackInfo{i};
+
+ dvec = (coms(hids(end),:) - coms(hids(1),:));
+
+ dx = sqrt(sum(dvec.^2, 2));
+ dt = (tsmp(times(hids(end))) - tsmp(times(hids(1))));
+
+ track_stats.slmot(i) = dx;
+ track_stats.slvel(i) = dx/dt;
+
+ track_stats.slvec(i,:) = dvec;
+ track_stats.slpos(i,:) = coms(hids(1),:);
+ end
+end
diff --git a/+Stats/make_all_trackinfo.m b/+Stats/make_all_trackinfo.m
new file mode 100644
index 0000000000000000000000000000000000000000..1fd73ad792dc68e1b7e57f025a4ef4def99d9ce4
--- /dev/null
+++ b/+Stats/make_all_trackinfo.m
@@ -0,0 +1,6 @@
+ROOT = 'E:\Smadar\LeverJS';
+
+flist = dir(fullfile(ROOT,'*.LEVER'));
+for i=1:length(flist)
+ Stats.make_trackinfo(ROOT,flist(i).name);
+end
diff --git a/+Stats/make_movie_table.m b/+Stats/make_movie_table.m
new file mode 100644
index 0000000000000000000000000000000000000000..2d8d2934956ed56d235e70eeddba06b3c4a7f8de
--- /dev/null
+++ b/+Stats/make_movie_table.m
@@ -0,0 +1,47 @@
+flist = dir('dataset_ecto_stats-*.mat');
+fnames = {flist.name}.';
+
+load(fnames{end});
+
+exp_types = arrayfun(@(x)(regexpi(x.abrv,'(\w*?)_','tokens','once')), ecto_size_stats, 'UniformOutput',false);
+exp_types = vertcat(exp_types{:});
+
+tsmp = arrayfun(@(x)(median(diff(x.tsmp))), ecto_size_stats);
+
+bExcludeMotion = (tsmp >= 7.0);
+exclMotSrt = repmat({'A'}, length(bExcludeMotion),1);
+exclMotSrt(bExcludeMotion) = {'B'};
+
+inclMotLbl = repmat({'Y'}, length(bExcludeMotion),1);
+inclMotLbl(bExcludeMotion) = {'N'};
+
+chkrows = arrayfun(@(x,y)([x,y]), exp_types, exclMotSrt, 'UniformOutput',false);
+chkrows = vertcat(chkrows{:});
+
+[~,srtidx] = sortrows(chkrows);
+
+fid = fopen('sup_table_info.csv','wt');
+fprintf(fid, 'Dataset,Treatment,dt (sec),First Frame Vesicles (Mesoderm),Motion Statistics,First Frame Vesicles (Ectoderm),Total Vesicles (Mesoderm),Total Vesicles (Ectoderm)\n');
+for i=1:length(srtidx)
+ idx = srtidx(i);
+
+ datasetName = ecto_size_stats(idx).name;
+
+ ectoVesicleCount = ecto_size_stats(idx).ecto.count;
+ skelVesicleCount = ecto_size_stats(idx).skel.count;
+
+ if ( ~bExcludeMotion(idx) )
+ ectoTotalCount = length(ecto_track_stats(idx).ecto.mean_sizes);
+ skelTotalCount = length(ecto_track_stats(idx).skel.mean_sizes);
+
+ fprintf(fid, '%s,%s,%.2f,%s,%d,%d,%d,%d\n', datasetName, exp_types{idx}, tsmp(idx), inclMotLbl{idx},...
+ skelVesicleCount, ectoVesicleCount, skelTotalCount, ectoTotalCount);
+ else
+ fprintf(fid, '%s,%s,%.2f,%s,%d,%d\n', datasetName, exp_types{idx}, tsmp(idx), inclMotLbl{idx}, skelVesicleCount, ectoVesicleCount);
+ end
+end
+
+fclose(fid);
+
+%
+%
\ No newline at end of file
diff --git a/+Stats/make_trackinfo.m b/+Stats/make_trackinfo.m
new file mode 100644
index 0000000000000000000000000000000000000000..474dd3c07be2e352f483ecf1fe57d191c83cc53c
--- /dev/null
+++ b/+Stats/make_trackinfo.m
@@ -0,0 +1,92 @@
+function make_trackinfo(ROOT,filename)
+ leverFile = fullfile(ROOT,filename);
+
+ conn = database(leverFile, '','', 'org.sqlite.JDBC', 'jdbc:sqlite:');
+ CONSTANTS = Read.getConstants(conn);
+
+ cutoff = 60*15;
+ bKeep = CONSTANTS.imageData.TimeStampDelta(1,1,:) <= cutoff;
+
+ keeplen = nnz(bKeep);
+ datasetName = CONSTANTS.imageData.DatasetName;
+
+ %% Load track information from cells db table
+ fprintf('Loading tracking info... ');
+ cellTrackinfo = getCellTracks(conn);
+ if ( isempty(cellTrackinfo) )
+ return;
+ end
+
+ szCells = max(cellTrackinfo.cellID);
+ fprintf('done\n');
+
+ %% Load edge information from db
+ fprintf('Loading edge info... ');
+ edgeInfo = getEdgeInfo(conn, szCells);
+ save([datasetName '_edgeinfo.mat'], 'edgeInfo');
+ fprintf('done\n');
+
+ %% Track lengths
+ trackLengths = zeros(1,max([cellTrackinfo.tid]));
+ for i=1:length(cellTrackinfo.tid)
+ trackLengths(cellTrackinfo.tid(i)) = trackLengths(cellTrackinfo.tid(i))+1;
+ end
+
+ validTracks = find(trackLengths > 5);
+ trackMap = containers.Map(validTracks,1:length(validTracks));
+
+ prg = Utils.CmdlnProgress(length(cellTrackinfo.tid));
+ trackInfo = cell(1,length(validTracks));
+ trackHash = cell(1,keeplen);
+ for i=1:length(cellTrackinfo.tid)
+ if ( cellTrackinfo.time(i) > keeplen )
+ prg.PrintProgress(i);
+ continue;
+ end
+
+ if ( ~isKey(trackMap,cellTrackinfo.tid(i)) )
+ prg.PrintProgress(i);
+ continue;
+ end
+
+ idx = trackMap(cellTrackinfo.tid(i));
+
+ trackInfo{idx} = [trackInfo{idx} cellTrackinfo.cellID(i)];
+ trackHash{cellTrackinfo.time(i)} = [trackHash{cellTrackinfo.time(i)}; idx];
+
+ prg.PrintProgress(i);
+ end
+
+ prg.ClearProgress(true);
+ save([datasetName '_trackinfo.mat'], 'trackInfo', 'trackHash', 'validTracks', 'trackMap');
+end
+
+function trackInfo = getCellTracks(conn)
+ trackInfo = [];
+
+ cmd=['SELECT cellID,time,trackID FROM tblCells ORDER BY time,cellID'];
+ Q=ljsFetch(conn,cmd);
+
+ if ( isempty(Q) )
+ return;
+ end
+
+ trackInfo.cellID = vertcat(Q{:,1});
+ trackInfo.time = vertcat(Q{:,2});
+ trackInfo.tid = vertcat(Q{:,3});
+end
+
+function edgeInfo = getEdgeInfo(conn, szCells)
+ edgeInfo = [];
+
+ cmd=['SELECT cellID_src,cellID_dst,cost FROM tblCosts ORDER BY cellID_src'];
+ Q=ljsFetch(conn,cmd);
+
+ edgeInfo.costs = sparse([Q{:,1}],[Q{:,2}],[Q{:,3}], szCells,szCells);
+
+
+ cmd=['SELECT cellID_src,cellID_dst,cost FROM tblDistCC ORDER BY cellID_src'];
+ Q=ljsFetch(conn,cmd);
+ edgeInfo.dists = sparse([Q{:,1}],[Q{:,2}],[Q{:,3}], szCells,szCells);
+end
+
diff --git a/+Stats/parsave.m b/+Stats/parsave.m
new file mode 100644
index 0000000000000000000000000000000000000000..e8e2de7265d6bf9c8b8052e3b1f41d06eb7bddbd
--- /dev/null
+++ b/+Stats/parsave.m
@@ -0,0 +1,9 @@
+function parsave(prefix, didx, names, vals)
+ sv_struct = [];
+
+ for i=1:length(names)
+ sv_struct.(names{i}) = vals{i};
+ end
+
+ save([prefix num2str(didx,'%02d') '.mat'], '-struct','sv_struct');
+end
diff --git a/+Stats/plot_directional_stats.m b/+Stats/plot_directional_stats.m
new file mode 100644
index 0000000000000000000000000000000000000000..47014013c2f7b18be765909fdee68242bd04c7cd
--- /dev/null
+++ b/+Stats/plot_directional_stats.m
@@ -0,0 +1,43 @@
+flist = dir('dataset_collate_tracks-*.mat');
+fnames = {flist.name}.';
+
+s = load(fnames{end});
+
+da = cell(length(s.track_data),1);
+for i=1:length(s.track_data)
+ dx = s.track_data(i).dx;
+ ldx = sqrt(sum(dx.^2, 2));
+
+ bValid = ldx > 1e-8;
+ dx = dx(bValid,:);
+ ldx = ldx(bValid);
+
+ ndx = dx ./ sqrt(ldx);
+
+ dirc = dot(ndx(1:end-1,:),ndx(2:end,:), 2);
+ na = cross(ndx(1:end-1,:),ndx(2:end,:), 2);
+
+ nna = na ./ sqrt(sum(na.^2, 2));
+ dsgn = [1; sign(dot(nna(1:end-1,:),nna(2:end,:), 2))];
+ dsgn(dsgn==0) = 1;
+
+ da{i} = acosd(dirc) .* dsgn;
+end
+
+de = [];
+bEcto = [];
+bDMSO = [];
+
+figure;
+abins = linspace(-180,180,6);
+for i=1:length(da)
+ if ( length(da{i}) < 20 )
+ continue;
+ end
+
+ bEcto = [bEcto; strcmpi(s.track_data(i).region, 'ecto')];
+ bDMSO = [bDMSO; strcmpi(s.track_data(i).type, 'DMSO')];
+
+% n = histcounts(da{i}, abins);
+ hist(da{i});
+end
diff --git a/+Stats/plot_ecto_stats.m b/+Stats/plot_ecto_stats.m
new file mode 100644
index 0000000000000000000000000000000000000000..de36c1418532a650bf06940ffbd772ffb8c7d04b
--- /dev/null
+++ b/+Stats/plot_ecto_stats.m
@@ -0,0 +1,918 @@
+%%
+dmso_order = {'23-06-2016_TimeLapse1_DMSO_20hrs_Calcein_FM464';
+ '28-06-2016_TimeLapse1_Pos2_DMSO_18hrs_Calcein_FM464';
+ '28-06-2016_TimeLapse1_Pos1_DMSO_18hrs_Calcein_FM464';
+ '23-06-2016_TimeLapse2_DMSO_20hrs_Calcein_AzepRH';
+ '28-06-2016_TimeLapse2_DMSO_22hrs_Calcein_FM464';
+ '22-06-2016_TimeLapse1_DMSO_24hrs_Calcein_FM464';
+ '22-06-2016_TimeLapse2_DMSO_24hrs_Calcein_FM464';
+ '22-06-2016_TimeLapse3_DMSO_20hrs_Calcein_FM464';
+ '28-06-2016_TimeLapse3_Pos1_DMSO_25hrs_Calcein_FM464';
+ '28-06-2016_TimeLapse3_Pos3_DMSO_25hrs_Calcein_FM464';
+ '28-06-2016_TimeLapse3_Pos4_DMSO_25hrs_Calcein_FM464';
+ '24-06-2016_Timelapse4_DMSO_24hrs_Calcein_FM464';
+ '01-07-2016_TimeLapse3_DMSO_21hrs_Calcein_FM464';
+ '01-07-2016_TimeLapse4_DMSO_21hrs_Calcein_FM464';
+ '01-07-2016_TimeLapse5_DMSO_23hrs_Calcein_FM464'};
+
+axt_order = {'01-07-2016_TimeLapse1_Pos0_Axtinib_150_18hrs_Calcein_FM464';
+ '01-07-2016_TimeLapse1_Pos1_Axtinib_150_18hrs_Calcein_FM464';
+ '01-07-2016_TimeLapse2_Pos2_Axtinib_150_19hrs_Calcein_FM464';
+ '01-07-2016_TimeLapse2_Pos3_Axtinib_150_19hrs_Calcein_FM464';
+ '24-06-2016_TimeLapse1_Axtinib_150_19hrs_Calcein_FM464';
+ '24-06-2016_TimeLapse2_Axtinib_150_20hrs_Calcein_FM464';
+ '24-06-2016_TimeLapse3_Axtinib_150_21hrs_Calcein_FM464';
+ '30-06-2016_TimeLapse1_Pos0_Axtinib_150_19hrs_Calcein_FM464';
+ '30-06-2016_TimeLapse1_Pos1_Axtinib_150_19hrs_Calcein_FM464';
+ '30-06-2016_TimeLapse1_Pos2_Axtinib_150_19hrs_Calcein_FM464';
+ '23-06-2016_Timelapse_3_Axitinib_36hr_Calcein_FM464'};
+
+bb94_order = {'29-06-2016_TimeLapse1_Pos0_BB94_17hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse1_Pos1_BB94_17hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse1_Pos2_BB94_17hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse2_Pos0_BB94_21hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse2_Pos1_BB94_21hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse2_Pos2_BB94_21hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse3_Pos0_BB94_23hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse3_Pos2_BB94_23hrs_Calcein_FM464';
+ '29-06-2016_TimeLapse3_Pos3_BB94_23hrs_Calcein_FM464';
+ '30-06-2016_TimeLapse2_Pos0_BB94_21hrs_Calcein_FM464';
+ '30-06-2016_TimeLapse2_Pos1_BB94_21hrs_Calcein_FM464';
+ '30-06-2016_TimeLapse2_Pos2_BB94_21hrs_Calcein_FM464';
+ '30-06-2016_TimeLapse2_Pos3_BB94_21hrs_Calcein_FM464'};
+
+all_order = [dmso_order; axt_order; bb94_order];
+
+flist = dir('dataset_ecto_stats-*.mat');
+fnames = {flist.name}.';
+
+load(fnames{end});
+
+%% Sort and only keep matched tracks for now
+matchCells = arrayfun(@(x)(find(strcmpi(all_order,x.name))), ecto_size_stats, 'UniformOutput',false);
+matchTracks = arrayfun(@(x)(find(strcmpi(all_order,x.name))), ecto_track_stats, 'UniformOutput',false);
+
+bMatchedS = cellfun(@(x)(~isempty(x)), matchCells);
+bMatchedT = cellfun(@(x)(~isempty(x)), matchTracks);
+
+chk_size_stats = ecto_size_stats(bMatchedS);
+chk_track_stats = ecto_track_stats(bMatchedT);
+
+[~,invMatch] = sort(vertcat(matchCells{bMatchedS}));
+chk_size_stats = chk_size_stats(invMatch);
+
+[~,invMatch] = sort(vertcat(matchTracks{bMatchedT}));
+chk_track_stats = chk_track_stats(invMatch);
+
+%% Drop empty data
+bValidEctoS = arrayfun(@(x)(~isempty(x.ecto.vsizes)), chk_size_stats);
+bValidSkelS = arrayfun(@(x)(~isempty(x.skel.vsizes)), chk_size_stats);
+
+bValidEctoT = arrayfun(@(x)(~isempty(x.ecto.mean_sizes)), chk_track_stats);
+bValidSkelT = arrayfun(@(x)(~isempty(x.skel.mean_sizes)), chk_track_stats);
+
+bError = (bValidEctoS & ~bValidSkelS);
+err_names = arrayfun(@(x)(x.abrv), chk_size_stats(bError), 'UniformOutput',false);
+if ( any(bError) )
+ fprintf('WARNING: Invalid data detected: No skeletogenic segmentations: (%s)\n', strjoin(err_names,','));
+end
+
+chk_size_stats = chk_size_stats(bValidEctoS);
+
+bError = (bValidEctoT & ~bValidSkelT);
+err_names = arrayfun(@(x)(x.abrv), chk_size_stats(bError), 'UniformOutput',false);
+if ( any(bError) )
+ fprintf('WARNING: Invalid data detected: No skeletogenic segmentations: (%s)\n', strjoin(err_names,','));
+end
+
+chk_track_stats = chk_track_stats(bValidEctoT);
+
+%%
+exp_types = arrayfun(@(x)(regexpi(x.abrv,'(\w*?)_','tokens','once')), chk_size_stats, 'UniformOutput',false);
+exp_types = vertcat(exp_types{:});
+
+ecto_size_info = vertcat(chk_size_stats.ecto);
+skel_size_info = vertcat(chk_size_stats.skel);
+endo_size_info = vertcat(chk_size_stats.endo);
+
+exp_types_track = arrayfun(@(x)(regexpi(x.abrv,'(\w*?)_','tokens','once')), chk_track_stats, 'UniformOutput',false);
+exp_types_track = vertcat(exp_types_track{:});
+
+ecto_track_info = vertcat(chk_track_stats.ecto);
+skel_track_info = vertcat(chk_track_stats.skel);
+endo_track_info = vertcat(chk_track_stats.endo);
+
+%% Interleave ecto/endo and skel
+inter_info = arrayfun(@(x)([x.ecto;x.skel;x.endo]), chk_size_stats, 'UniformOutput',false);
+inter_info = vertcat(inter_info{:});
+
+inter_grp_names = arrayfun(@(x)({[x.abrv '_ecto'];[x.abrv '_skel'];[x.abrv '_ecto']}), chk_size_stats, 'UniformOutput',false);
+inter_grp_names = vertcat(inter_grp_names{:});
+
+inter_sizes = vertcat(inter_info.vsizes);
+inter_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.vsizes,1),1)), inter_grp_names, inter_info, 'UniformOutput',false);
+inter_grp = vertcat(inter_grp_cell{:});
+
+
+%% Make sizes/group names for boxplots
+ecto_sizes = vertcat(ecto_size_info.vsizes);
+skel_sizes = vertcat(skel_size_info.vsizes);
+endo_sizes = vertcat(endo_size_info.vsizes);
+
+ecto_grp_cell = arrayfun(@(x,y)(repmat({x.abrv},size(y.vsizes,1),1)), chk_size_stats, ecto_size_info, 'UniformOutput',false);
+ecto_grp = vertcat(ecto_grp_cell{:});
+
+skel_grp_cell = arrayfun(@(x,y)(repmat({x.abrv},size(y.vsizes,1),1)), chk_size_stats, skel_size_info, 'UniformOutput',false);
+skel_grp = vertcat(skel_grp_cell{:});
+
+endo_grp_cell = arrayfun(@(x,y)(repmat({x.abrv},size(y.vsizes,1),1)), chk_size_stats, endo_size_info, 'UniformOutput',false);
+endo_grp = vertcat(endo_grp_cell{:});
+
+exp_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.vsizes,1),1)), [exp_types;exp_types], [ecto_size_info;skel_size_info], 'UniformOutput',false);
+exp_grp = vertcat(exp_grp_cell{:});
+
+%% Make groups for track stats
+ecto_track_sz = vertcat(ecto_track_info.mean_sizes);
+skel_track_sz = vertcat(skel_track_info.mean_sizes);
+endo_track_sz = vertcat(endo_track_info.mean_sizes);
+
+ecto_grp_cell = arrayfun(@(x,y)(repmat({x.abrv},size(y.mean_sizes,1),1)), chk_track_stats, ecto_track_info, 'UniformOutput',false);
+ecto_track_grp = vertcat(ecto_grp_cell{:});
+
+skel_grp_cell = arrayfun(@(x,y)(repmat({x.abrv},size(y.mean_sizes,1),1)), chk_track_stats, skel_track_info, 'UniformOutput',false);
+skel_track_grp = vertcat(skel_grp_cell{:});
+
+endo_grp_cell = arrayfun(@(x,y)(repmat({x.abrv},size(y.mean_sizes,1),1)), chk_track_stats, endo_track_info, 'UniformOutput',false);
+endo_track_grp = vertcat(endo_grp_cell{:});
+
+exp_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.mean_sizes,1),1)), [exp_types_track;exp_types_track], [ecto_track_info;skel_track_info], 'UniformOutput',false);
+exp_track_grp = vertcat(exp_grp_cell{:});
+
+
+%%
+close all
+
+
+
+% %%
+% figure;
+% boxplot(inter_sizes, inter_grp, 'notch','on');
+% xtickangle(90);
+% zoom('reset');
+% ylim([0.1,1.0]);
+% title('Vesicle size [um^3]');
+%
+% % draw alternating patches
+% pcol = [0,0,0];
+% palpha = 0.1;
+% pwidth = 2.0;
+% pheight = 1.0;
+%
+% xstart = 0.5;
+% for i=1:length(inter_info)
+% if ( mod(i-1,2) == 0 )
+% pcol = [0,1,0];
+% else
+% pcol = [0,0,0];
+% end
+%
+% x0 = xstart;
+% x1 = xstart + pwidth;
+%
+% y0 = 0;
+% y1 = pheight;
+%
+% patch([x0,x1,x1,x0],[y0,y0,y1,y1],pcol, 'FaceColor',pcol, 'FaceAlpha',palpha, 'EdgeColor','none');
+%
+% xstart = xstart + pwidth;
+% end
+
+%%
+wdir_grp_names = arrayfun(@(x)(num2str(x+2)), (1:length(ecto_track_info(1).avg_wdir_index{1})).', 'UniformOutput',false);
+
+ecto_wdiridx = [];
+ecto_pl_diridx = [];
+ecto_wdirgrp = {};
+
+ecto_all_didx = vertcat(ecto_track_info.avg_wdir_index);
+ecto_all_plen = vertcat(ecto_track_info.avg_path_wlen);
+ecto_all_mdisp = vertcat(ecto_track_info.avg_max_wdist);
+
+ecto_all_lens = cellfun(@(x)(length(x)), ecto_all_didx);
+for i=1:length(wdir_grp_names)
+ bValid = (ecto_all_lens >= i);
+
+ didx = cellfun(@(x)(x(i)), ecto_all_didx(bValid));
+
+ ecto_wdiridx = [ecto_wdiridx; didx];
+ ecto_wdirgrp = [ecto_wdirgrp; repmat(wdir_grp_names(i), size(didx,1),1)];
+end
+
+figure;
+boxplot(ecto_wdiridx,ecto_wdirgrp, 'notch','on');
+
+%%
+didx_idx = 8;
+
+exp_ecto_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.mean_sizes,1),1)), exp_types_track, ecto_track_info, 'UniformOutput',false);
+exp_ecto_grp = vertcat(exp_ecto_grp_cell{:});
+
+ecto_wdir_cell = vertcat(ecto_track_info.avg_wdir_index);
+ecto_wdiridx = cellfun(@(x)(x(didx_idx)), ecto_wdir_cell);
+
+bBB94_ecto = cellfun(@(x)(strcmpi(x,'BB94')), exp_ecto_grp);
+exp_ecto_grp = exp_ecto_grp(~bBB94_ecto);
+ecto_wdiridx = ecto_wdiridx(~bBB94_ecto);
+
+figure;
+boxplot(ecto_wdiridx, exp_ecto_grp, 'notch','on');
+xtickangle(30);
+zoom('reset');
+% ylim([0.3,0.6]);
+title('Directionality Index - Ectoderm');
+
+%
+exp_skel_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.mean_sizes,1),1)), exp_types_track, skel_track_info, 'UniformOutput',false);
+exp_skel_grp = vertcat(exp_skel_grp_cell{:});
+
+skel_wdir_cell = vertcat(skel_track_info.avg_wdir_index);
+skel_wdiridx = cellfun(@(x)(x(didx_idx)), skel_wdir_cell);
+
+bBB94_skel = cellfun(@(x)(strcmpi(x,'BB94')), exp_skel_grp);
+exp_skel_grp = exp_skel_grp(~bBB94_skel);
+skel_wdiridx = skel_wdiridx(~bBB94_skel);
+
+figure;
+boxplot(skel_wdiridx, exp_skel_grp, 'notch','on');
+xtickangle(30);
+zoom('reset');
+% ylim([0.3,0.6]);
+title('Directionality Index - Mesoderm');
+
+%
+exp_endo_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.mean_sizes,1),1)), exp_types_track, endo_track_info, 'UniformOutput',false);
+exp_endo_grp = vertcat(exp_endo_grp_cell{:});
+
+endo_wdir_cell = vertcat(endo_track_info.avg_wdir_index);
+endo_wdiridx = cellfun(@(x)(x(didx_idx)), endo_wdir_cell);
+
+bBB94_ecto = cellfun(@(x)(strcmpi(x,'BB94')), exp_endo_grp);
+exp_endo_grp = exp_endo_grp(~bBB94_ecto);
+endo_wdiridx = endo_wdiridx(~bBB94_ecto);
+
+figure;
+boxplot(endo_wdiridx, exp_endo_grp, 'notch','on');
+xtickangle(30);
+zoom('reset');
+% ylim([0.3,0.6]);
+title('Directionality Index - Endoderm');
+
+%%
+exp_ecto_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.mean_sizes,1),1)), exp_types_track, ecto_track_info, 'UniformOutput',false);
+exp_ecto_grp = vertcat(exp_ecto_grp_cell{:});
+
+ecto_msdD = vertcat(ecto_track_info.msd_D);
+ecto_msdR = vertcat(ecto_track_info.msd_R);
+ecto_msdsig = vertcat(ecto_track_info.msd_sig);
+
+bBB94_ecto = cellfun(@(x)(strcmpi(x,'BB94')), exp_ecto_grp);
+exp_ecto_grp = exp_ecto_grp(~bBB94_ecto);
+ecto_msdD = ecto_msdD(~bBB94_ecto);
+ecto_msdR = ecto_msdR(~bBB94_ecto);
+ecto_msdsig = ecto_msdsig(~bBB94_ecto);
+
+figure;
+boxplot(ecto_msdD, exp_ecto_grp, 'notch','on');
+xtickangle(30);
+zoom('reset');
+ylim([0,0.04]);
+title('Ideal Diffusion Coefficient [um^2/s] - Ectoderm');
+
+% figure;
+% boxplot(ecto_msdsig, exp_ecto_grp, 'notch','on');
+% xtickangle(30);
+% zoom('reset');
+% ylim([0,2.5]);
+% title('Ideal Diffusion Variance - Ectoderm');
+
+figure;
+hist(ecto_msdR,100);
+title('Ideal Diffusion Fit (Pearson) - Ectoderm')
+
+%
+exp_skel_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.mean_sizes,1),1)), exp_types_track, skel_track_info, 'UniformOutput',false);
+exp_skel_grp = vertcat(exp_skel_grp_cell{:});
+
+skel_msdD = vertcat(skel_track_info.msd_D);
+skel_msdR = vertcat(skel_track_info.msd_R);
+skel_msdsig = vertcat(skel_track_info.msd_sig);
+
+bBB94_ecto = cellfun(@(x)(strcmpi(x,'BB94')), exp_skel_grp);
+exp_skel_grp = exp_skel_grp(~bBB94_ecto);
+skel_msdD = skel_msdD(~bBB94_ecto);
+skel_msdR = skel_msdR(~bBB94_ecto);
+skel_msdsig = skel_msdsig(~bBB94_ecto);
+
+figure;
+boxplot(skel_msdD, exp_skel_grp, 'notch','on');
+xtickangle(30);
+zoom('reset');
+ylim([0,0.04]);
+title('Ideal Diffusion Coefficient [um^2/s] - Mesoderm');
+
+% figure;
+% boxplot(skel_msdsig, exp_skel_grp, 'notch','on');
+% xtickangle(30);
+% zoom('reset');
+% ylim([0,2.5]);
+% title('Ideal Diffusion Variance - Mesoderm');
+
+figure;
+hist(skel_msdR.^2,100);
+title('Ideal Diffusion Fit (R^2) - Mesoderm')
+
+bDMSO_skel = cellfun(@(x)(strcmpi(x,'DMSO')), exp_skel_grp);
+
+figure;hold on;
+plot(skel_msdD(bDMSO_skel),skel_msdR(bDMSO_skel).^2, '.k');
+plot(skel_msdD(~bDMSO_skel),skel_msdR(~bDMSO_skel).^2, '.r');
+
+test = 1;
+
+% %%
+% didx_idx = 8;
+%
+% exp_ecto_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.mean_sizes,1),1)), exp_types_track, ecto_track_info, 'UniformOutput',false);
+% exp_ecto_grp = vertcat(exp_ecto_grp_cell{:});
+%
+% ecto_wpathrat_cell = vertcat(ecto_track_info.avg_wpath_ratio);
+% ecto_wpathrat = cellfun(@(x)(x(didx_idx)), ecto_wpathrat_cell);
+%
+% bBB94_ecto = cellfun(@(x)(strcmpi(x,'BB94')), exp_ecto_grp);
+% exp_ecto_grp = exp_ecto_grp(~bBB94_ecto);
+% ecto_wpathrat = ecto_wpathrat(~bBB94_ecto);
+%
+% figure;
+% boxplot(ecto_wpathrat, exp_ecto_grp, 'notch','on');
+% xtickangle(30);
+% zoom('reset');
+% % ylim([0.3,0.6]);
+% title('Windowed Path Ratio - Ectoderm');
+%
+% %
+% exp_skel_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.mean_sizes,1),1)), exp_types_track, skel_track_info, 'UniformOutput',false);
+% exp_skel_grp = vertcat(exp_skel_grp_cell{:});
+%
+% skel_wpathrat_cell = vertcat(skel_track_info.avg_wpath_ratio);
+% skel_wpathrat = cellfun(@(x)(x(didx_idx)), skel_wdir_cell);
+%
+% bBB94_skel = cellfun(@(x)(strcmpi(x,'BB94')), exp_skel_grp);
+% exp_skel_grp = exp_skel_grp(~bBB94_skel);
+% skel_wpathrat = skel_wpathrat(~bBB94_skel);
+%
+% figure;
+% boxplot(skel_wpathrat, exp_skel_grp, 'notch','on');
+% xtickangle(30);
+% zoom('reset');
+% % ylim([0.3,0.6]);
+% title('Windowed Path Ratio - Mesoderm');
+
+%%
+exp_ecto_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.mean_sizes,1),1)), exp_types_track, ecto_track_info, 'UniformOutput',false);
+exp_ecto_grp = vertcat(exp_ecto_grp_cell{:});
+exp_skel_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.mean_sizes,1),1)), exp_types_track, skel_track_info, 'UniformOutput',false);
+exp_skel_grp = vertcat(exp_skel_grp_cell{:});
+
+cmap = lines(4);
+bEctoDMSO = [cellfun(@(x)(strcmpi(x,'DMSO')), exp_ecto_grp); false(size(exp_skel_grp))];
+bEctoAxt = [cellfun(@(x)(strcmpi(x,'Axtinib')), exp_ecto_grp); false(size(exp_skel_grp))];
+bSkelDMSO = [false(size(exp_ecto_grp)); cellfun(@(x)(strcmpi(x,'DMSO')), exp_skel_grp)];
+bSkelAxt = [false(size(exp_ecto_grp)); cellfun(@(x)(strcmpi(x,'Axtinib')), exp_skel_grp)];
+
+didx_idx = 8;
+
+all_mean_sizes = [vertcat(ecto_track_info.mean_sizes); vertcat(skel_track_info.mean_sizes)];
+all_mean_ivel = [vertcat(ecto_track_info.mean_inst_vel); vertcat(skel_track_info.mean_inst_vel)];
+all_diff_coeff = [vertcat(ecto_track_info.msd_D); vertcat(skel_track_info.msd_D)];
+
+all_mean_didx = [vertcat(ecto_track_info.avg_wdir_index); vertcat(skel_track_info.avg_wdir_index)];
+all_mean_didx = cellfun(@(x)(x(didx_idx)), all_mean_didx);
+
+all_mean_mdidx = [vertcat(ecto_track_info.avg_wdir_mindex); vertcat(skel_track_info.avg_wdir_mindex)];
+all_mean_mdidx = cellfun(@(x)(x(didx_idx)), all_mean_mdidx);
+
+all_wpath_ratio = [vertcat(ecto_track_info.avg_wpath_ratio); vertcat(skel_track_info.avg_wpath_ratio)];
+all_wpath_ratio = cellfun(@(x)(x(didx_idx)), all_wpath_ratio);
+
+all_path_len = [vertcat(ecto_track_info.path_len); vertcat(skel_track_info.path_len)];
+all_direct_len = [vertcat(ecto_track_info.direct_len); vertcat(skel_track_info.direct_len)];
+
+all_path_ratio = all_direct_len ./ all_path_len;
+
+%%
+[~,ctr] = hist(all_mean_ivel(bSkelDMSO), ctr);
+figure;
+hist(all_mean_ivel(bSkelDMSO), ctr);
+xlim([0,0.3]);
+title('Avg. Speed [um/s]: Control - Mesoderm');
+saveas(gcf, 'C:/Users/mwinter/Desktop/figs/speed_hist_control_mesoderm.png','png');
+%
+figure;
+hist(all_mean_ivel(bSkelAxt), ctr);
+xlim([0,0.3]);
+title('Avg. Speed [um/s]: VEGFR Inhib. - Mesoderm');
+saveas(gcf, 'C:/Users/mwinter/Desktop/figs/speed_hist_axt_mesoderm.png','png');
+%
+figure;
+hist(all_mean_ivel(bEctoDMSO), ctr);
+xlim([0,0.3]);
+title('Avg. Speed [um/s]: Control - Ectoderm');
+saveas(gcf, 'C:/Users/mwinter/Desktop/figs/speed_hist_control_ectoderm.png','png');
+%
+figure;
+hist(all_mean_ivel(bEctoAxt), ctr);
+xlim([0,0.3]);
+title('Avg. Speed [um/s]: VEGFR Inhib. - Ectoderm');
+saveas(gcf, 'C:/Users/mwinter/Desktop/figs/speed_hist_axt_ectoderm.png','png');
+
+% %%
+% figure;
+% hist(log(all_mean_ivel(bSkelDMSO)), 25);
+% title('Log-Avg. Speed: Control - Mesoderm');
+%
+% figure;
+% hist(log(all_mean_ivel(bEctoDMSO)), 25);
+% title('Log-Avg. Speed: Control - Ectoderm');
+%
+% figure;
+% hist(log(all_mean_ivel(bSkelAxt)), 25);
+% title('Log-Avg. Speed: VEGFR Inhib. - Mesoderm');
+%
+% figure;
+% hist(log(all_mean_ivel(bEctoAxt)), 25);
+% title('Log-Avg. Speed: VEGFR Inhib. - Ectoderm');
+
+% %%
+% figure;
+% hist(log(all_mean_ivel(bEctoDMSO)), 25);
+% title('Log-Avg. Speed [um/s]: VEGFR Inhib. - Ectoderm');
+
+%% Plot size vs. average instantaneous speed
+figure;hold on;
+plot(all_mean_sizes(bSkelDMSO),all_mean_ivel(bSkelDMSO), '.', 'Color',cmap(2,:));
+% title('Avg. Speed [um/s] vs. Avg. Size [um]');
+xlim([0.1,1.3]);
+ylim([0,0.2]);
+%
+plot(all_mean_sizes(bSkelAxt),all_mean_ivel(bSkelAxt), '.', 'Color',cmap(1,:));
+
+[b,~,r,~,rst] = regress(all_mean_ivel(bSkelAxt|bSkelDMSO), [ones(size(all_mean_ivel(bSkelAxt|bSkelDMSO))), all_mean_sizes(bSkelAxt|bSkelDMSO)]);
+xx = 0:0.1:1.3;
+yy = b(1) + xx*b(2);
+plot(xx,yy, '-k', 'LineWidth',2);
+legend('Mesoderm: Control','Mesoderm: VEGFR Inhib.');
+saveas(gcf,'C:\Users\mwinter\Desktop\figs\sup2_meso_speed_v_size.png');
+
+% figure;plot(all_mean_sizes(bEctoDMSO|bSkelDMSO), r, '.r');
+%%
+figure;hold on;
+plot(all_mean_sizes(bEctoDMSO),all_mean_ivel(bEctoDMSO), '.', 'Color',cmap(2,:));
+% title('Avg. Speed [um/s] vs. Avg. Size [um]');
+xlim([0.1,1.3]);
+ylim([0,0.2]);
+%
+plot(all_mean_sizes(bEctoAxt),all_mean_ivel(bEctoAxt), '.', 'Color',cmap(1,:));
+
+[b,~,r,~,rst] = regress(all_mean_ivel(bEctoAxt|bEctoDMSO), [ones(size(all_mean_ivel(bEctoAxt|bEctoDMSO))), all_mean_sizes(bEctoAxt|bEctoDMSO)]);
+xx = 0:0.1:1.3;
+yy = b(1) + xx*b(2);
+plot(xx,yy, '-k', 'LineWidth',2);
+legend('Ectoderm: Control','Ectoderm: VEGFR Inhib.');
+saveas(gcf,'C:\Users\mwinter\Desktop\figs\sup2_ecto_speed_v_size.png');
+
+
+%% Mean size vs. diffusion coeff
+figure;hold on;
+plot(all_mean_sizes(bSkelDMSO),all_diff_coeff(bSkelDMSO), '.', 'Color',cmap(2,:));
+% title('Avg. Speed [um/s] vs. Avg. Size [um]');
+xlim([0.1,1.3]);
+ylim([0,0.05]);
+%
+plot(all_mean_sizes(bSkelAxt),all_diff_coeff(bSkelAxt), '.', 'Color',cmap(1,:));
+
+[b,~,r,~,rst] = regress(all_diff_coeff(bSkelAxt|bSkelDMSO), [ones(size(all_diff_coeff(bSkelAxt|bSkelDMSO))), all_mean_sizes(bSkelAxt|bSkelDMSO)]);
+xx = 0:0.1:1.3;
+yy = b(1) + xx*b(2);
+plot(xx,yy, '-k', 'LineWidth',2);
+legend('Mesoderm: Control','Mesoderm: VEGFR Inhib.');
+saveas(gcf,'C:\Users\mwinter\Desktop\figs\sup3_meso_diff_v_size.png');
+
+% figure;plot(all_mean_sizes(bEctoDMSO|bSkelDMSO), r, '.r');
+%%
+figure;hold on;
+plot(all_mean_sizes(bEctoDMSO),all_diff_coeff(bEctoDMSO), '.', 'Color',cmap(2,:));
+% title('Avg. Speed [um/s] vs. Avg. Size [um]');
+xlim([0.1,1.3]);
+ylim([0,0.05]);
+%
+plot(all_mean_sizes(bEctoAxt),all_diff_coeff(bEctoAxt), '.', 'Color',cmap(1,:));
+
+[b,~,r,~,rst] = regress(all_diff_coeff(bEctoAxt|bEctoDMSO), [ones(size(all_diff_coeff(bEctoAxt|bEctoDMSO))), all_mean_sizes(bEctoAxt|bEctoDMSO)]);
+xx = 0:0.1:1.3;
+yy = b(1) + xx*b(2);
+plot(xx,yy, '-k', 'LineWidth',2);
+legend('Ectoderm: Control','Ectoderm: VEGFR Inhib.');
+saveas(gcf,'C:\Users\mwinter\Desktop\figs\sup3_ecto_diff_v_size.png');
+
+
+
+%%
+figure;hold on;
+plot(all_mean_sizes(bEctoDMSO),all_path_ratio(bEctoDMSO), '.', 'Color',cmap(1,:));
+plot(all_mean_sizes(bEctoAxt),all_path_ratio(bEctoAxt), '.', 'Color',cmap(2,:));
+plot(all_mean_sizes(bSkelDMSO),all_path_ratio(bSkelDMSO), '.', 'Color',cmap(3,:));
+plot(all_mean_sizes(bSkelAxt),all_path_ratio(bSkelAxt), '.', 'Color',cmap(4,:));
+legend('DMSO - Ecto','Axtinib - Ecto','DMSO - Meso','Axtinib - Meso');
+title('Path Ratio vs. Avg. Size [um]');
+xlim([0.1,1.3]);
+
+%%
+figure;hold on;
+plot(all_mean_sizes(bEctoDMSO),all_mean_didx(bEctoDMSO), '.', 'Color',cmap(1,:));
+plot(all_mean_sizes(bEctoAxt),all_mean_didx(bEctoAxt), '.', 'Color',cmap(2,:));
+plot(all_mean_sizes(bSkelDMSO),all_mean_didx(bSkelDMSO), '.', 'Color',cmap(3,:));
+plot(all_mean_sizes(bSkelAxt),all_mean_didx(bSkelAxt), '.', 'Color',cmap(4,:));
+legend('DMSO - Ecto','Axtinib - Ecto','DMSO - Meso','Axtinib - Meso');
+title('Directed Index vs. Avg. Size [um]');
+xlim([0.1,1.3]);
+
+% %%
+% figure;hold on;
+% plot(all_mean_sizes(bEctoDMSO),all_mean_mdidx(bEctoDMSO), '.', 'Color',cmap(1,:));
+% plot(all_mean_sizes(bEctoAxt),all_mean_mdidx(bEctoAxt), '.', 'Color',cmap(2,:));
+% plot(all_mean_sizes(bSkelDMSO),all_mean_mdidx(bSkelDMSO), '.', 'Color',cmap(3,:));
+% plot(all_mean_sizes(bSkelAxt),all_mean_mdidx(bSkelAxt), '.', 'Color',cmap(4,:));
+% legend('DMSO - Ecto','Axtinib - Ecto','DMSO - Meso','Axtinib - Meso');
+% title('Directed Index (MeanDisp) vs. Avg. Size [um]');
+% xlim([0.1,1.3]);
+
+%%
+figure;hold on;
+plot(all_mean_sizes(bEctoDMSO),all_wpath_ratio(bEctoDMSO), '.', 'Color',cmap(1,:));
+plot(all_mean_sizes(bEctoAxt),all_wpath_ratio(bEctoAxt), '.', 'Color',cmap(2,:));
+plot(all_mean_sizes(bSkelDMSO),all_wpath_ratio(bSkelDMSO), '.', 'Color',cmap(3,:));
+plot(all_mean_sizes(bSkelAxt),all_wpath_ratio(bSkelAxt), '.', 'Color',cmap(4,:));
+legend('DMSO - Ecto','Axtinib - Ecto','DMSO - Meso','Axtinib - Meso');
+title('Avg. Windowed Path Ratio vs. Avg. Size [um]');
+xlim([0.1,1.3]);
+
+
+test = 1;
+
+%%
+exp_ecto_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.vsizes,1),1)), exp_types, ecto_size_info, 'UniformOutput',false);
+exp_ecto_grp = vertcat(exp_ecto_grp_cell{:});
+
+figure;
+boxplot(ecto_sizes, exp_ecto_grp, 'notch','on');
+xtickangle(30);
+zoom('reset');
+ylim([0.0,1.0]);
+title('Vesicle Size - Ectoderm (N=3057,7663,6050) [um^3]');
+
+%
+exp_skel_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.vsizes,1),1)), exp_types, skel_size_info, 'UniformOutput',false);
+exp_skel_grp = vertcat(exp_skel_grp_cell{:});
+
+figure;
+boxplot(skel_sizes, exp_skel_grp, 'notch','on');
+xtickangle(30);
+zoom('reset');
+ylim([0.0,1.0]);
+title('Vesicle Size - Mesoderm (N=1443,2125,1832) [um^3]');
+
+%
+all_sizes = [ecto_sizes;skel_sizes];
+
+exp_all_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.vsizes,1),1)), [exp_types;exp_types], [ecto_size_info;skel_size_info], 'UniformOutput',false);
+exp_all_grp = vertcat(exp_all_grp_cell{:});
+
+figure;
+boxplot(all_sizes, exp_all_grp, 'notch','on');
+xtickangle(30);
+zoom('reset');
+ylim([0.0,1.0]);
+title('Vesicle Size - All (N=4500,9788,7882) [um^3]');
+
+%%
+figure;
+boxplot(skel_sizes, skel_grp, 'notch','on');
+xtickangle(30);
+zoom('reset');
+ylim([0.1,1.0]);
+title('Vesicle size (mesoderm) [um^3]');
+
+%%
+figure;
+boxplot(ecto_sizes, ecto_grp, 'notch','on');
+xtickangle(30);
+zoom('reset');
+ylim([0.1,1.0]);
+title('Vesicle size (ectoderm) [um^3]');
+
+%%
+figure;
+stem(vertcat(skel_size_info.count),'k');
+set(gca, 'TickLabelInterpreter','none');
+xticks(1:length(chk_size_stats));
+xticklabels('manual');
+xticklabels({chk_size_stats.abrv});
+xtickangle(30)
+title('Vesicle counts (mesoderm) [um^3]');
+
+%%
+figure;
+stem(vertcat(skel_size_info.count)./vertcat(ecto_size_info.rsize),'k');
+set(gca, 'TickLabelInterpreter','none');
+xticks(1:length(chk_size_stats));
+xticklabels('manual');
+xticklabels({chk_size_stats.abrv});
+xtickangle(30)
+title('Volume-normalized vesicle counts (mesoderm)');
+
+%%
+figure;
+stem(vertcat(ecto_size_info.count),'k');
+set(gca, 'TickLabelInterpreter','none');
+xticks(1:length(chk_size_stats));
+xticklabels('manual');
+xticklabels({chk_size_stats.abrv});
+xtickangle(30)
+title('Vesicle counts (ectoderm) [um^3]');
+
+%%
+figure;
+stem(vertcat(ecto_size_info.count)./vertcat(ecto_size_info.rsize),'k');
+set(gca, 'TickLabelInterpreter','none');
+xticks(1:length(chk_size_stats));
+xticklabels('manual');
+xticklabels({chk_size_stats.abrv});
+xtickangle(30)
+title('Volume-normalized vesicle counts (ectoderm)');
+
+% %%
+% figure;
+% stem(vertcat(ecto_info.rsize),'k');
+% set(gca, 'TickLabelInterpreter','none');
+% xticks(1:length(chk_size_stats));
+% xticklabels('manual');
+% xticklabels({chk_size_stats.abrv});
+% xtickangle(30)
+% title('Volume (ectoderm) [um^3]');
+
+% %%
+% figure;
+% stem(vertcat(skel_size_info.count)./vertcat(ecto_size_info.count),'k');
+% set(gca, 'TickLabelInterpreter','none');
+% xticks(1:length(chk_size_stats));
+% xticklabels('manual');
+% xticklabels({chk_size_stats.abrv});
+% xtickangle(30)
+% title('Vesicle count ratio (mesoderm/ectoderm)');
+%
+% %%
+% figure;
+% boxplot(vertcat(skel_size_info.count)./vertcat(ecto_size_info.count),exp_types);
+% xtickangle(30)
+% title('Vesicle count ratio (mesoderm/ectoderm)');
+
+% %%
+% ecto_vol = arrayfun(@(x)(sum(x.vsizes)), ecto_size_info);
+% skel_vol = arrayfun(@(x)(sum(x.vsizes)), skel_size_info);
+%
+% figure;
+% boxplot(skel_vol ./ ecto_vol, exp_types);
+% xtickangle(30)
+% title('Vesicle volume ratio (mesoderm/ectoderm)');
+%
+% %%
+% figure;
+% boxplot(skel_vol, exp_types);
+% xtickangle(30)
+% title('Vesicle volume (mesoderm)');
+%
+% %%
+% figure;
+% boxplot(ecto_vol, exp_types);
+% xtickangle(30)
+% title('Vesicle volume (ectoderm)');
+%
+% %%
+% figure;
+% boxplot(vertcat(ecto_size_info.rsize), exp_types);
+% xtickangle(30)
+% title('Ectoderm volume');
+
+% %%
+% figure;
+% stem(vertcat(skel_info.rsize)./vertcat(ecto_info.rsize),'k');
+% set(gca, 'TickLabelInterpreter','none');
+% xticks(1:length(chk_size_stats));
+% xticklabels('manual');
+% xticklabels({chk_size_stats.abrv});
+% xtickangle(30)
+% title('Volume ratio (mesoderm/ectoderm)');
+%
+% %%
+% figure;
+% stem((vertcat(skel_size_info.count)./vertcat(ecto_size_info.count)) ./ (vertcat(skel_size_info.rsize)./vertcat(ecto_size_info.rsize)),'k');
+% set(gca, 'TickLabelInterpreter','none');
+% xticks(1:length(chk_size_stats));
+% xticklabels('manual');
+% xticklabels({chk_size_stats.abrv});
+% xtickangle(30)
+% title('Count per volume ratio (mesoderm/ectoderm)');
+
+%% Vesicle path ratio
+%%
+skel_path_rat = vertcat(skel_track_info.direct_len) ./ vertcat(skel_track_info.path_len);
+
+figure;
+boxplot(skel_path_rat, skel_track_grp, 'notch','on');
+xtickangle(30);
+zoom('reset');
+ylim([0.1,1.0]);
+title('Vesicle path ratio (mesoderm)');
+
+%%
+ecto_path_rat = vertcat(ecto_track_info.direct_len) ./ vertcat(ecto_track_info.path_len);
+
+figure;
+boxplot(ecto_path_rat, ecto_track_grp, 'notch','on');
+xtickangle(30);
+zoom('reset');
+ylim([0.1,1.0]);
+title('Vesicle path ratio (ectoderm)');
+
+%% Vesicle path ratio -- bulk groups
+%%
+exp_ecto_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.mean_sizes,1),1)), exp_types_track, ecto_track_info, 'UniformOutput',false);
+exp_ecto_grp = vertcat(exp_ecto_grp_cell{:});
+
+bBB94_ecto = cellfun(@(x)(strcmpi(x,'BB94')), exp_ecto_grp);
+exp_ecto_grp = exp_ecto_grp(~bBB94_ecto);
+ecto_path_rat = ecto_path_rat(~bBB94_ecto);
+
+figure;
+boxplot(ecto_path_rat, exp_ecto_grp, 'notch','on');
+xtickangle(30);
+zoom('reset');
+ylim([0.0,1.0]);
+title('Path Length Ratio - Ectoderm (N=5539,7369)');
+
+%
+exp_skel_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.mean_sizes,1),1)), exp_types_track, skel_track_info, 'UniformOutput',false);
+exp_skel_grp = vertcat(exp_skel_grp_cell{:});
+
+bBB94_skel = cellfun(@(x)(strcmpi(x,'BB94')), exp_skel_grp);
+exp_skel_grp = exp_skel_grp(~bBB94_skel);
+skel_path_rat = skel_path_rat(~bBB94_skel);
+
+figure;
+boxplot(skel_path_rat, exp_skel_grp, 'notch','on');
+xtickangle(30);
+zoom('reset');
+ylim([0.0,1.0]);
+title('Path Length Ratio - Mesoderm (N=2606,2096)');
+
+%
+all_path_rat = [vertcat(ecto_track_info.direct_len);vertcat(skel_track_info.direct_len)] ./ [vertcat(ecto_track_info.path_len);vertcat(skel_track_info.path_len)];
+
+exp_all_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.mean_sizes,1),1)), [exp_types_track;exp_types_track], [ecto_track_info;skel_track_info], 'UniformOutput',false);
+exp_all_grp = vertcat(exp_all_grp_cell{:});
+
+bBB94_all = cellfun(@(x)(strcmpi(x,'BB94')), exp_all_grp);
+exp_all_grp = exp_all_grp(~bBB94_all);
+all_path_rat = all_path_rat(~bBB94_all);
+
+figure;
+boxplot(all_path_rat, exp_all_grp, 'notch','on');
+xtickangle(30);
+zoom('reset');
+ylim([0.0,1.0]);
+title('Path Length Ratio - All (N=8145,9465)');
+
+%% Instantaneous Velocity -- bulk groups
+%%
+exp_ecto_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.mean_sizes,1),1)), exp_types_track, ecto_track_info, 'UniformOutput',false);
+exp_ecto_grp = vertcat(exp_ecto_grp_cell{:});
+
+ecto_inst_vel = vertcat(ecto_track_info.mean_inst_vel);
+
+bBB94_ecto = cellfun(@(x)(strcmpi(x,'BB94')), exp_ecto_grp);
+exp_ecto_grp = exp_ecto_grp(~bBB94_ecto);
+ecto_inst_vel = ecto_inst_vel(~bBB94_ecto);
+
+figure;
+boxplot(ecto_inst_vel, exp_ecto_grp, 'notch','on');
+xtickangle(30);
+zoom('reset');
+ylim([0.0,0.2]);
+title('Instantaneous Velocity - Ectoderm (N=5539,7369)[um/s]');
+
+%
+exp_skel_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.mean_sizes,1),1)), exp_types_track, skel_track_info, 'UniformOutput',false);
+exp_skel_grp = vertcat(exp_skel_grp_cell{:});
+
+skel_inst_vel = vertcat(skel_track_info.mean_inst_vel);
+
+bBB94_skel = cellfun(@(x)(strcmpi(x,'BB94')), exp_skel_grp);
+exp_skel_grp = exp_skel_grp(~bBB94_skel);
+skel_inst_vel = skel_inst_vel(~bBB94_skel);
+
+figure;
+boxplot(skel_inst_vel, exp_skel_grp, 'notch','on');
+xtickangle(30);
+zoom('reset');
+ylim([0.0,0.2]);
+title('Instantaneous Velocity - Mesoderm (N=2606,2096)[um/s]');
+
+%
+all_inst_vel = [vertcat(ecto_track_info.mean_inst_vel); vertcat(skel_track_info.mean_inst_vel)];
+
+exp_all_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.mean_sizes,1),1)), [exp_types_track;exp_types_track], [ecto_track_info;skel_track_info], 'UniformOutput',false);
+exp_all_grp = vertcat(exp_all_grp_cell{:});
+
+bBB94_all = cellfun(@(x)(strcmpi(x,'BB94')), exp_all_grp);
+exp_all_grp = exp_all_grp(~bBB94_all);
+all_inst_vel = all_inst_vel(~bBB94_all);
+
+figure;
+boxplot(all_inst_vel, exp_all_grp, 'notch','on');
+xtickangle(30);
+zoom('reset');
+ylim([0.0,0.2]);
+title('Instantaneous Velocity - All (N=8145,9465)[um/s]');
+
+
+%%
+all_inst_vel = [vertcat(ecto_track_info.mean_inst_vel); vertcat(skel_track_info.mean_inst_vel)];
+
+figure;
+boxplot(all_inst_vel, [ecto_track_grp;skel_track_grp], 'notch','on');
+xtickangle(30);
+zoom('reset');
+ylim([0.0,0.15]);
+title('Vesicle average instantaneous velocity (all) [um/sec]');
+
+%% Instantaneous velocity comparison
+%%
+exp_ecto_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.mean_sizes,1),1)), exp_types_track, ecto_track_info, 'UniformOutput',false);
+exp_ecto_grp = vertcat(exp_ecto_grp_cell{:});
+
+figure;
+boxplot(vertcat(ecto_track_info.mean_inst_vel), ecto_track_grp, 'notch','on');
+xtickangle(30);
+zoom('reset');
+ylim([0.0,0.15]);
+title('Vesicle average instantaneous velocity (ectoderm) [um/sec]');
+
+%
+exp_skel_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.mean_sizes,1),1)), exp_types_track, skel_track_info, 'UniformOutput',false);
+exp_skel_grp = vertcat(exp_skel_grp_cell{:});
+
+figure;
+boxplot(vertcat(skel_track_info.mean_inst_vel), skel_track_grp, 'notch','on');
+xtickangle(30);
+zoom('reset');
+ylim([0.0,0.15]);
+title('Vesicle average instantaneous velocity (mesoderm) [um/sec]');
+
+
+%% Average instantaneous velocity
+%%
+exp_grp_cell = arrayfun(@(x,y)(repmat(x,size(y.mean_sizes,1),1)), exp_types_track, skel_track_info, 'UniformOutput',false);
+exp_grp = vertcat(exp_skel_grp_cell{:});
+
+
+%% Vesicle brownian variance
+%%
+skel_brn_var = vertcat(skel_track_info.brownian_var);
+
+figure;
+boxplot(skel_brn_var, skel_track_grp, 'notch','on');
+xtickangle(30);
+zoom('reset');
+ylim([0.1,1.0]);
+title('Vesicle brownian variance (mesoderm)');
+
+%%
+ecto_brn_var = vertcat(ecto_track_info.brownian_var);
+
+figure;
+boxplot(ecto_brn_var, ecto_track_grp, 'notch','on');
+xtickangle(30);
+zoom('reset');
+ylim([0.1,1.0]);
+title('Vesicle brownian variance (ectoderm)');
+
diff --git a/+Stats/plot_motion_size.m b/+Stats/plot_motion_size.m
new file mode 100644
index 0000000000000000000000000000000000000000..84bd69a80e034b6e27291e18538c1055bba0941a
--- /dev/null
+++ b/+Stats/plot_motion_size.m
@@ -0,0 +1,100 @@
+flist = dir('dataset_collate_tracks-*.mat');
+fnames = {flist.name}.';
+
+s = load(fnames{end});
+
+% figure;hold on;
+% for i=1:length(s.track_data)
+% szs = s.track_data(i).sizes;
+% h = ones(5,1);
+% h = h / sum(h);
+%
+% smsz = conv(szs, h, 'same');
+% nrm = conv(ones(size(szs,1),1), h, 'same');
+%
+% smsz = smsz ./ nrm;
+% plot(smsz, '-k');
+% end
+
+window_size = 2*60;
+
+
+avg_length = 15;
+
+bDMSO = strcmpi({s.track_data.type}.', 'DMSO');
+bSkel = strcmpi({s.track_data.region}.', 'skel');
+
+s.track_data = s.track_data(bDMSO);
+
+bSlowing = false(length(s.track_data),1);
+slowSz = zeros(length(s.track_data),1);
+for i=1:length(s.track_data)
+ if ( length(s.track_data(i).times) < avg_length )
+ continue;
+ end
+
+ mlen = min(avg_length, length(s.track_data(i).msd)-2);
+
+ bIdx = 1:mlen;
+ eIdx = (length(s.track_data(i).msd)-mlen+1):length(s.track_data(i).msd);
+
+ avgBeg = mean(sqrt(s.track_data(i).msd(bIdx)) ./ s.track_data(i).dt(bIdx));
+ avgEnd = mean(sqrt(s.track_data(i).msd(eIdx)) ./ s.track_data(i).dt(eIdx));
+
+ bSlowing(i) = (avgEnd <= 0.5*avgBeg);
+ slowSz(i) = avgEnd-avgBeg;
+end
+
+bShrinking = false(length(s.track_data),1);
+shrinkSz = zeros(length(s.track_data),1);
+for i=1:length(s.track_data)
+ mlen = min(avg_length, length(s.track_data(i).sizes)-2);
+
+ bIdx = 1:mlen;
+ eIdx = (length(s.track_data(i).sizes)-mlen+1):length(s.track_data(i).sizes);
+
+ avgBeg = mean(s.track_data(i).sizes(bIdx));
+ avgEnd = mean(s.track_data(i).sizes(eIdx));
+
+ bShrinking(i) = (avgEnd <= 0.5*avgBeg);
+ shrinkSz(i) = avgEnd-avgBeg;
+end
+
+figure;hist(shrinkSz(bShrinking&bSlowing),25);
+figure;hist(slowSz(bShrinking&bSlowing),25);
+
+cmap = lines(10);
+
+figure;hold on;
+shsl = s.track_data(bShrinking&bSlowing);
+for i=1:length(shsl)
+ fprintf('%s:%d\n', shsl(i).name, shsl(i).tid);
+ mlen = min(avg_length, length(shsl(i).msd)-2);
+
+ ivels = sqrt(shsl(i).msd) ./ shsl(i).dt;
+ h = ones(mlen,1);
+ h = h / sum(h);
+
+ smvel = conv(ivels, h, 'same');
+ nrm = conv(ones(size(ivels,1),1), h, 'same');
+
+ smvel = smvel ./ nrm;
+ plot(smvel, '-', 'Color',cmap(mod(i-1,10)+1,:));
+ title('Smoothed instantaneous speed');
+end
+
+figure;hold on;
+for i=1:length(shsl)
+ mlen = min(avg_length, length(shsl(i).sizes)-2);
+
+ isz = shsl(i).sizes;
+ h = ones(mlen,1);
+ h = h / sum(h);
+
+ smsz = conv(isz, h, 'same');
+ nrm = conv(ones(size(isz,1),1), h, 'same');
+
+ smsz = smsz ./ nrm;
+ plot(smsz, '-', 'Color',cmap(mod(i-1,10)+1,:));
+ title('Smoothed size');
+end
diff --git a/+Stats/plot_vel_histograms.m b/+Stats/plot_vel_histograms.m
new file mode 100644
index 0000000000000000000000000000000000000000..5d0dcc0661e43641befe03854b74fcd5642b79ac
--- /dev/null
+++ b/+Stats/plot_vel_histograms.m
@@ -0,0 +1,89 @@
+flist = dir('dataset_collate_tracks-*.mat');
+fnames = {flist.name}.';
+
+s = load(fnames{end});
+
+bSkel = arrayfun(@(x)(strcmpi(x.region,'skel')), s.track_data);
+bEcto = arrayfun(@(x)(strcmpi(x.region,'ecto')), s.track_data);
+
+bDMSO = arrayfun(@(x)(strcmpi(x.type,'DMSO')), s.track_data);
+bAxt = arrayfun(@(x)(strcmpi(x.type,'Axtinib')), s.track_data);
+
+dmso_skel_data = s.track_data(bDMSO&bSkel);
+dmso_ecto_data = s.track_data(bDMSO&bEcto);
+axt_skel_data = s.track_data(bAxt&bSkel);
+axt_ecto_data = s.track_data(bAxt&bEcto);
+
+dmso_skel_ivel_c = arrayfun(@(x)(sum(x.dx.^2,2) ./ x.dt), dmso_skel_data, 'UniformOutput',false);
+dmso_ecto_ivel_c = arrayfun(@(x)(sum(x.dx.^2,2) ./ x.dt), dmso_ecto_data, 'UniformOutput',false);
+axt_skel_ivel_c = arrayfun(@(x)(sum(x.dx.^2,2) ./ x.dt), axt_skel_data, 'UniformOutput',false);
+axt_ecto_ivel_c = arrayfun(@(x)(sum(x.dx.^2,2) ./ x.dt), axt_ecto_data, 'UniformOutput',false);
+
+dmso_skel_ivel = vertcat(dmso_skel_ivel_c{:});
+dmso_ecto_ivel = vertcat(dmso_ecto_ivel_c{:});
+axt_skel_ivel = vertcat(axt_skel_ivel_c{:});
+axt_ecto_ivel = vertcat(axt_ecto_ivel_c{:});
+
+nbins = 100;
+
+figure;
+subplot(2,2,1)
+hist((dmso_skel_ivel), 100*nbins);
+xlim([0,0.1]);
+title('Instantaneous Velocities - Control - Mesoderm')
+
+subplot(2,2,2)
+hist((dmso_ecto_ivel), 100*nbins);
+xlim([0,0.1]);
+title('Instantaneous Velocities - Control - Ectoderm')
+
+subplot(2,2,3)
+hist((axt_skel_ivel), 100*nbins);
+xlim([0,0.1]);
+title('Instantaneous Velocities - VEGFR Inhib. - Mesoderm')
+
+subplot(2,2,4)
+hist((axt_ecto_ivel), 100*nbins);
+xlim([0,0.1]);
+title('Instantaneous Velocities - VEGFR Inhib. - Ectoderm')
+
+
+%%
+figure;
+%%
+subplot(2,2,1)
+hist(log10(dmso_skel_ivel), nbins);
+oldL = get(gca, 'XTickLabel');
+newL = cellfun(@(x)(['10^{' x '}']), oldL,'UniformOutput',false);
+set(gca, 'XTickLabel',newL);
+title('Instantaneous Speed - Control - Mesoderm')
+xlabel('Vesicle Speed [\mum/s]')
+
+subplot(2,2,2)
+hist(log10(dmso_ecto_ivel), nbins);
+oldL = get(gca, 'XTickLabel');
+newL = cellfun(@(x)(['10^{' x '}']), oldL,'UniformOutput',false);
+set(gca, 'XTickLabel',newL);
+title('Instantaneous Speed - Control - Ectoderm')
+xlabel('Vesicle Speed [\mum/s]')
+
+subplot(2,2,3)
+hist(log10(axt_skel_ivel), nbins);
+hist(log10(dmso_ecto_ivel), nbins);
+oldL = get(gca, 'XTickLabel');
+newL = cellfun(@(x)(['10^{' x '}']), oldL,'UniformOutput',false);
+set(gca, 'XTickLabel',newL);
+title('Instantaneous Speed - VEGFR Inhib. - Mesoderm')
+xlabel('Vesicle Speed [\mum/s]')
+
+subplot(2,2,4)
+hist(log10(axt_ecto_ivel), nbins);
+hist(log10(dmso_ecto_ivel), nbins);
+oldL = get(gca, 'XTickLabel');
+newL = cellfun(@(x)(['10^{' x '}']), oldL,'UniformOutput',false);
+set(gca, 'XTickLabel',newL);
+title('Instantaneous Speed - VEGFR Inhib. - Ectoderm')
+xlabel('Vesicle Speed [\mum/s]')
+
+saveas(gcf, 'C:\Users\mwinter\Desktop\figs\sup_figure2.png');
+
diff --git a/+Stats/write_compare.m b/+Stats/write_compare.m
new file mode 100644
index 0000000000000000000000000000000000000000..1d603d1e72f10b65dc706e72a2f68286603d73b0
--- /dev/null
+++ b/+Stats/write_compare.m
@@ -0,0 +1,55 @@
+function tbl = write_compare(fid, stats_tbl, sum_tbl, vname, title, types)
+ sgrp = cellfun(@(x,y)([x ':' y]), stats_tbl.region, stats_tbl.type, 'UniformOutput',false);
+ [p,tbl,stats] = kruskalwallis(stats_tbl.(vname), sgrp, 'off');
+ [c,~,~,gn] = multcompare(stats, 'CType','dunn-sidak', 'Display','off');
+
+ m = size(c,1);
+ pt = zeros(m,1);
+ for i=1:m
+ bga = strcmp(sgrp, gn{c(i,1)});
+ bgb = strcmp(sgrp, gn{c(i,2)});
+
+ pt(i) = 1 - (1-ranksum(stats_tbl.(vname)(bga),stats_tbl.(vname)(bgb)))^m;
+ end
+
+ gt = {};
+ midx = [];
+
+ %% Assumed regions/types (to get proper sorting)
+ regions = {'Mesoderm';'Ectoderm'};
+% types = {'Control';'VEGFR Inhibition'};
+ for i=1:length(regions)
+ for j=1:(length(types)-1)
+ t2 = types((j+1):end);
+ np = cellfun(@(x)({[regions{i} ':' types{j}], [regions{i} ':' x]}), t2, 'UniformOutput',false);
+
+ gt = [gt; vertcat(np{:})];
+ end
+ end
+
+ for i=1:length(types)
+ for j=1:(length(regions)-1)
+ r2 = regions((j+1):end);
+ np = cellfun(@(x)({[regions{j} ':' types{i}], [x ':' types{i}]}), r2, 'UniformOutput',false);
+
+ gt = [gt; vertcat(np{:})];
+ end
+ end
+
+ fprintf(fid, '%s\n', title);
+ fprintf(fid, 'Kruskal-Wallis, %g\n\n', p);
+
+% fprintf(fid, '\n');
+ fprintf(fid, 'Pair-wise (Dunn-Sidak), Difference, P-Value\n');
+ for i=1:size(gt,1)
+ ia = find(strcmpi(gn,gt{i,1}));
+ ib = find(strcmpi(gn,gt{i,2}));
+
+ ma = sum_tbl.median(ia);
+ mb = sum_tbl.median(ib);
+
+ bridx = all((c(:,1:2)==[ia,ib]), 2);
+ fprintf(fid, '%s to %s,%f,%g\n', gt{i,1},gt{i,2}, mb-ma, c(bridx,6));
+ end
+ fprintf(fid, '\n\n');
+end
diff --git a/+Stats/write_compare_spicule.m b/+Stats/write_compare_spicule.m
new file mode 100644
index 0000000000000000000000000000000000000000..41a834a1077bbc2e443ebcaa7ce95962a30f7a77
--- /dev/null
+++ b/+Stats/write_compare_spicule.m
@@ -0,0 +1,32 @@
+function [pKW,pM] = write_compare_spicule(fid, spc_tbl, sum_tbl, fieldname, cmp_title)
+ [pKW,~,stats] = kruskalwallis(spc_tbl.(fieldname), spc_tbl.sp_dist,'off');
+ [c,~,~,gnames] = multcompare(stats, 'CType','dunn-sidak', 'Display','off');
+
+ pM = ones(length(gnames),length(gnames));
+ for i=1:size(c,1)
+ pM(c(i,1),c(i,2)) = c(i,6);
+ pM(c(i,2),c(i,1)) = c(i,6);
+ end
+
+ fprintf(fid, '%s\n', cmp_title);
+ fprintf(fid, 'Kruskal-Wallis, %g\n\n', pKW);
+
+ fprintf(fid, 'Pair-Wise Difference\n');
+ fprintf(fid, 'Median (Col - Row),%s\n', strjoin(gnames,','));
+ for i=1:size(gnames,1)
+ diff_row = (sum_tbl.median - sum_tbl.median(i));
+
+ nstr = arrayfun(@(x)(num2str(x,'%f')), diff_row, 'UniformOutput',false);
+ fprintf(fid, '%s,%s\n', gnames{i}, strjoin(nstr,','));
+ end
+ fprintf(fid, '\n\n');
+
+
+ fprintf(fid, 'Pair-wise (Dunn-Sidak)\n');
+ fprintf(fid, 'P-Value,%s\n', strjoin(gnames,','));
+ for i=1:size(gnames,1)
+ nstr = arrayfun(@(x)(num2str(x,'%f')), pM(i,:), 'UniformOutput',false);
+ fprintf(fid, '%s,%s\n', gnames{i}, strjoin(nstr,','));
+ end
+ fprintf(fid, '\n\n');
+end
diff --git a/+Stats/write_summary.m b/+Stats/write_summary.m
new file mode 100644
index 0000000000000000000000000000000000000000..07aa6aa0c322f852c3461bc0eb7059bceee96575
--- /dev/null
+++ b/+Stats/write_summary.m
@@ -0,0 +1,40 @@
+function stbl = write_summary(fid, stats_tbl, vname, title, types)
+ %% Assumed regions/types (to get proper sorting)
+ regions = {'Mesoderm';'Ectoderm'};
+% types = {'Control';'VEGFR Inhibition'};
+
+ tnames = cell(length(regions),length(types));
+ rnames = cell(length(regions),length(types));
+ stats_means = zeros(length(regions),length(types));
+ stats_meds = zeros(length(regions),length(types));
+ for i=1:length(regions)
+ for j=1:length(types)
+ bReg = strcmpi(stats_tbl.region,regions{i});
+ bType = strcmpi(stats_tbl.type,types{j});
+ chk_stats = stats_tbl.(vname);
+
+ tnames{i,j} = types{j};
+ rnames{i,j} = regions{i};
+ stats_means(i,j) = mean(chk_stats(bReg&bType));
+ stats_meds(i,j) = median(chk_stats(bReg&bType));
+ end
+ end
+
+ stbl = table(tnames(:),rnames(:),stats_means(:),stats_meds(:), 'VariableNames',{'type','region','mean','median'});
+
+ typeStr = strjoin(types, ',');
+ fprintf(fid, 'Average %s,%s\n', title, typeStr);
+ for i=1:length(regions)
+ nstr = arrayfun(@(x)(num2str(x,'%f')), stats_means(i,:), 'UniformOutput',false);
+ fprintf(fid, '%s,%s\n', regions{i}, strjoin(nstr,','));
+ end
+ fprintf(fid, '\n');
+
+ fprintf(fid, 'Median %s,%s\n', title, typeStr);
+ for i=1:length(regions)
+ nstr = arrayfun(@(x)(num2str(x,'%f')), stats_meds(i,:), 'UniformOutput',false);
+ fprintf(fid, '%s,%s\n', regions{i}, strjoin(nstr,','));
+ end
+ fprintf(fid, '\n');
+
+end
diff --git a/+Stats/write_summary_spicule.m b/+Stats/write_summary_spicule.m
new file mode 100644
index 0000000000000000000000000000000000000000..8607da8f9847038605d1f970dee91898dea3a149
--- /dev/null
+++ b/+Stats/write_summary_spicule.m
@@ -0,0 +1,27 @@
+function stbl = write_summary_spicule(fid, spc_tbl, fieldname, cmp_title)
+ distances = unique(spc_tbl.sp_dist);
+
+ stats_means = zeros(length(distances),1);
+ stats_meds = zeros(length(distances),1);
+ for i=1:length(distances)
+ bDist = (spc_tbl.sp_dist == distances(i));
+
+ stats_means(i) = mean(spc_tbl.(fieldname)(bDist));
+ stats_meds(i) = median(spc_tbl.(fieldname)(bDist));
+ end
+
+ stbl = table(distances(:),stats_means(:),stats_meds(:), 'VariableNames',{'sp_distance','mean','median'});
+
+ groupStr = 'Average Spicule Distance';
+ fprintf(fid, '%s,Average %s\n', groupStr, cmp_title);
+ for i=1:length(distances)
+ fprintf(fid, '%d,%f\n', distances(i), stats_means(i));
+ end
+ fprintf(fid, '\n');
+
+ fprintf(fid, '%s,Median %s\n', groupStr, cmp_title);
+ for i=1:length(distances)
+ fprintf(fid, '%d,%f\n', distances(i), stats_meds(i));
+ end
+ fprintf(fid, '\n');
+end
diff --git a/plot_axt.py b/plot_axt.py
new file mode 100644
index 0000000000000000000000000000000000000000..51d63cf450c16418dee5dbcf031ecdd30b569406
--- /dev/null
+++ b/plot_axt.py
@@ -0,0 +1,166 @@
+import pandas as pd
+import matplotlib
+import matplotlib.pyplot as plt
+import seaborn as sns
+# import statannot as sta
+
+data = pd.read_csv('ecto_size_stats_axtinib.csv')
+
+##
+plt.figure(figsize=(10,5))
+# matplotlib.rcParams['figure.dpi'] = 200
+ax = sns.boxplot(x="region", y="size", hue="type", showmeans=True, data=data, palette='Set1',
+ meanprops={"marker":"s","markerfacecolor":"white", "markeredgecolor":"black"})
+ax.set_ylim(0,1.2)
+ax.set_ylabel('')
+ax.set_xlabel('')
+# ax.set_title('Vesicle Average Volume [um^3]', fontsize=20)
+ax.set_title('')
+ax.legend().set_title('')
+ax.legend(loc='upper left')
+ax.tick_params(labelsize=16)
+# ax.set(xticklabels=['Control','VEGFR Inhibition'])
+# ax.yaxis.set_major_locator(matplotlib.ticker.MultipleLocator(5))
+# ax.yaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
+plt.setp(ax.get_legend().get_texts(), fontsize=16)
+plt.tight_layout()
+plt.savefig('C:/Users/mwinter/Desktop/figs/size_axt.svg')
+
+
+#####################
+data = pd.read_csv('ecto_track_stats_axtinib.csv')
+
+##
+plt.figure(figsize=(5,4.5))
+# matplotlib.rcParams['figure.dpi'] = 200
+ax = sns.boxplot(x="region", y="avg_speed", hue="type", showmeans=True, data=data, palette='Set1',
+ meanprops={"marker":"s","markerfacecolor":"white", "markeredgecolor":"black"})
+ax.set_ylim(0.001,0.2)
+ax.set_ylabel('')
+ax.set_xlabel('')
+ax.yaxis.set_major_locator(plt.MaxNLocator(4))
+# ax.set_title('Vesicle Average Speed [um/s]', fontsize=20)
+ax.set_title('')
+ax.legend().set_title('')
+ax.legend(loc='upper left')
+ax.tick_params(labelsize=16)
+# ax.set(xticklabels=['Control','VEGFR Inhibition'])
+# ax.yaxis.set_major_locator(matplotlib.ticker.MultipleLocator(5))
+# ax.yaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
+plt.setp(ax.get_legend().get_texts(), fontsize=16)
+plt.tight_layout()
+plt.savefig('C:/Users/mwinter/Desktop/figs/track_speed.svg')
+
+# ##
+# plt.figure()
+# ax = sns.boxplot(x="type", y="path_ratio", hue="region", data=data, palette='Set3')
+# # ax.set_ylim(0,1.0)
+# ax.set_ylabel('')
+# ax.set_xlabel('')
+# ax.set_title('Vesicle Path Ratio', fontsize=20)
+# ax.legend().set_title('')
+# ax.tick_params(labelsize=16)
+# # ax.yaxis.set_major_locator(matplotlib.ticker.MultipleLocator(5))
+# # ax.yaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
+# plt.setp(ax.get_legend().get_texts(), fontsize=18)
+
+# ##
+# plt.figure()
+# ax = sns.boxplot(x="type", y="win_path_ratio", hue="region", data=data, palette='Set1')
+# # ax.set_ylim(0,1.0)
+# ax.set_ylabel('')
+# ax.set_xlabel('')
+# ax.set_title('Vesicle Windowed Path Ratio', fontsize=20)
+# ax.legend().set_title('')
+# ax.tick_params(labelsize=16)
+# # ax.yaxis.set_major_locator(matplotlib.ticker.MultipleLocator(5))
+# # ax.yaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
+# plt.setp(ax.get_legend().get_texts(), fontsize=18)
+
+##
+plt.figure(figsize=(5,4.5))
+ax = sns.boxplot(x="region", y="dir_index", hue="type", showmeans=True, data=data, palette='Set1',
+ meanprops={"marker":"s","markerfacecolor":"white", "markeredgecolor":"black"})
+# ax.set_ylim(0,0.2)
+ax.set_ylabel('')
+ax.set_xlabel('')
+# ax.set_title('Vesicle Directionality Index', fontsize=20)
+ax.set_title('')
+ax.legend().set_title('')
+ax.legend(loc='upper left')
+ax.tick_params(labelsize=16)
+# ax.set(xticklabels=['Control','VEGFR Inhibition'])
+# ax.yaxis.set_major_locator(matplotlib.ticker.MultipleLocator(5))
+# ax.yaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
+plt.setp(ax.get_legend().get_texts(), fontsize=16)
+plt.tight_layout()
+plt.savefig('C:/Users/mwinter/Desktop/figs/track_di.svg')
+
+##
+plt.figure(figsize=(5,4.5))
+# matplotlib.rcParams['figure.dpi'] = 200
+ax = sns.boxplot(x="region", y="nimsd_D", hue="type", showmeans=True, data=data, palette='Set1',
+ meanprops={"marker":"s","markerfacecolor":"white", "markeredgecolor":"black"})
+ax.set_ylim(0.0,0.05)
+ax.set_ylabel('')
+ax.set_xlabel('')
+ax.yaxis.set_major_locator(plt.MaxNLocator(4))
+# ax.set_title('Vesicle Average Speed [um/s]', fontsize=20)
+# ax.set_title('Non-Ideal Diffusion Coefficient [um^2/s]')
+ax.legend().set_title('')
+ax.legend(loc='upper left')
+ax.tick_params(labelsize=16)
+# ax.set(xticklabels=['Control','VEGFR Inhibition'])
+# ax.yaxis.set_major_locator(matplotlib.ticker.MultipleLocator(5))
+# ax.yaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
+plt.setp(ax.get_legend().get_texts(), fontsize=16)
+plt.tight_layout()
+plt.savefig('C:/Users/mwinter/Desktop/figs/track_nidiff_D.svg')
+
+##
+plt.figure(figsize=(5,4.5))
+# matplotlib.rcParams['figure.dpi'] = 200
+ax = sns.boxplot(x="region", y="nimsd_alpha", hue="type", showmeans=True, data=data, palette='Set1',
+ meanprops={"marker":"s","markerfacecolor":"white", "markeredgecolor":"black"})
+ax.set_ylim(0.5,2.5)
+ax.set_ylabel('')
+ax.set_xlabel('')
+ax.yaxis.set_major_locator(plt.MaxNLocator(4))
+ax.legend().set_title('')
+ax.legend(loc='upper left')
+ax.tick_params(labelsize=16)
+# ax.set(xticklabels=['Control','VEGFR Inhibition'])
+# ax.yaxis.set_major_locator(matplotlib.ticker.MultipleLocator(5))
+# ax.yaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
+plt.setp(ax.get_legend().get_texts(), fontsize=16)
+plt.tight_layout()
+plt.savefig('C:/Users/mwinter/Desktop/figs/track_nidiff_alpha.svg')
+
+##
+plt.figure(figsize=(5,4.5))
+# matplotlib.rcParams['figure.dpi'] = 200
+ax = sns.boxplot(x="region", y="msd_D", hue="type", showmeans=True, data=data, palette='Set1',
+ meanprops={"marker":"s","markerfacecolor":"white", "markeredgecolor":"black"})
+ax.set_ylim(0,0.03)
+ax.set_ylabel('')
+ax.set_xlabel('')
+ax.yaxis.set_major_locator(plt.MaxNLocator(4))
+ax.legend().set_title('')
+ax.legend(loc='upper left')
+ax.tick_params(labelsize=16)
+# ax.set(xticklabels=['Control','VEGFR Inhibition'])
+# ax.yaxis.set_major_locator(matplotlib.ticker.MultipleLocator(5))
+# ax.yaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
+plt.setp(ax.get_legend().get_texts(), fontsize=16)
+plt.tight_layout()
+plt.savefig('C:/Users/mwinter/Desktop/figs/track_diff_D.svg')
+
+##
+plt.figure(figsize=(5,4.5))
+# matplotlib.rcParams['figure.dpi'] = 200
+msdr2 = [i**2 for i in data['msd_R']]
+ax = sns.distplot(msdr2, kde=False, norm_hist=False, color=[0.07,0.30,0.59], hist_kws={'alpha':1})
+ax.set_ylabel('')
+ax.set_xlabel('')
+plt.tight_layout()
+plt.savefig('C:/Users/mwinter/Desktop/figs/track_diff_R_hist.svg')
diff --git a/plot_bb94.py b/plot_bb94.py
new file mode 100644
index 0000000000000000000000000000000000000000..f9c7c466418692f612944e319cc6f2c569334dbb
--- /dev/null
+++ b/plot_bb94.py
@@ -0,0 +1,25 @@
+import pandas as pd
+import matplotlib
+import matplotlib.pyplot as plt
+import seaborn as sns
+
+data = pd.read_csv('ecto_size_stats_bb94.csv')
+
+##
+plt.figure(figsize=(10,5))
+# matplotlib.rcParams['figure.dpi'] = 200
+ax = sns.boxplot(x="region", y="size", hue="type", data=data, palette='Set1',
+ showmeans = True, meanprops={"marker":"s","markerfacecolor":"white", "markeredgecolor":"black"})
+ax.set_ylim(0,1.0)
+ax.set_ylabel('')
+ax.set_xlabel('')
+# ax.set_title('Vesicle Average Volume [um^3]', fontsize=20)
+ax.set_title('')
+ax.legend().set_title('')
+ax.legend(loc='upper left')
+ax.tick_params(labelsize=14)
+# ax.set(xticklabels=['Control','Metalloproteinase Inhibition'])
+# ax.yaxis.set_major_locator(matplotlib.ticker.MultipleLocator(5))
+# ax.yaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
+plt.setp(ax.get_legend().get_texts(), fontsize=14)
+plt.savefig('C:/Users/mwinter/Desktop/figs/mmpl/size_bb94.svg')
diff --git a/plot_spicule_stats.py b/plot_spicule_stats.py
new file mode 100644
index 0000000000000000000000000000000000000000..78ef25726e8bba1bba392ed3590773b88673959b
--- /dev/null
+++ b/plot_spicule_stats.py
@@ -0,0 +1,49 @@
+import pandas as pd
+import matplotlib
+import matplotlib.pyplot as plt
+import seaborn as sns
+
+#####################
+data = pd.read_csv('spicule_stats.csv')
+
+##
+plt.figure(figsize=(4.8,5))
+ax = sns.boxplot(x="sp_dist", y="abs_vel", data=data, color='cornflowerblue', showmeans=True,meanprops={"marker":"s","markerfacecolor":"white", "markeredgecolor":"black"})
+# ax.set_ylim(0,0.2)
+ax.set_ylabel('')
+ax.set_xlabel('')
+# ax.set_title('Vesicle Directionality Index', fontsize=20)
+ax.set_title('')
+# ax.legend().set_title('')
+# ax.legend(loc='upper left')
+ax.tick_params(labelsize=16)
+ax.set_xticklabels(['1','2','3','4','5','6','7','8','9','10'])
+# ax.set(xticklabels=['Control','VEGFR Inhibition'])
+# ax.yaxis.set_major_locator(matplotlib.ticker.MultipleLocator(5))
+# ax.yaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
+# plt.setp(ax.get_legend().get_texts(), fontsize=16)
+plt.tight_layout()
+plt.setp(ax.artists, edgecolor = 'k')
+plt.setp(ax.lines, color='k')
+plt.savefig('C:/Users/mwinter/Desktop/figs/sp_dist_v_vel.svg')
+
+##
+plt.figure(figsize=(4.8,5))
+ax = sns.boxplot(x="sp_dist", y="inst_vel", data=data, color='cornflowerblue', showmeans=True,meanprops={"marker":"s","markerfacecolor":"white", "markeredgecolor":"black"})
+ax.set_ylim(-0.4,0.4)
+ax.set_ylabel('')
+ax.set_xlabel('')
+# ax.set_title('Vesicle Directionality Index', fontsize=20)
+ax.set_title('')
+# ax.legend().set_title('')
+# ax.legend(loc='upper left')
+ax.tick_params(labelsize=16)
+ax.set_xticklabels(['1','2','3','4','5','6','7','8','9','10'])
+# ax.set(xticklabels=['Control','VEGFR Inhibition'])
+# ax.yaxis.set_major_locator(matplotlib.ticker.MultipleLocator(5))
+# ax.yaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
+# plt.setp(ax.get_legend().get_texts(), fontsize=16)
+plt.tight_layout()
+plt.setp(ax.artists, edgecolor = 'k')
+plt.setp(ax.lines, color='k')
+plt.savefig('C:/Users/mwinter/Desktop/figs/sp_dist_v_prjdist.svg')