diff --git a/+Stats/collate_tracks.m b/+Stats/collate_tracks.m
index eba30b8e69b9d478a0d55e19290cc46c1b7f63cf..3407f52d5f38980f295114753dc40ef9b62487d2 100644
--- a/+Stats/collate_tracks.m
+++ b/+Stats/collate_tracks.m
@@ -25,6 +25,12 @@ function track_data = collate_tracks(ROOT, filename)
     
     %% Estimate parallelogram mask for first-frame image
     imagePath = CONSTANTS.imageData.imageDir;
+    chkimd = MicroscopeData.ReadMetadataFile(fullfile(imagePath,filesep));
+    if ( isempty(chkimd) )
+        warning(['Unable to load image data for: ' datasetName]);
+        return;
+    end
+    
     im = MicroscopeData.Reader(imagePath, 'timeRange',[1,1]);
     
     parallelMask = Segment.llsm_mask(im);
@@ -52,9 +58,10 @@ function track_data = collate_tracks(ROOT, filename)
         return;
     end
     
-    %% Skip if movie is >= 7fps
+    %% Skip if movie is >= 7fps or < 5
     medTsmp = median(diff(tsmp));
-    if ( medTsmp >= 7.0 )
+    if ( medTsmp >= 7.0 || medTsmp < 5 )
+%     if ( medTsmp >= 7.0 )
         return;
     end
     
@@ -86,15 +93,17 @@ function track_data = collate_tracks(ROOT, filename)
     
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     %% Read per-frame vesicle info
-    bEcto = false(0,1);
-    bSkel = false(0,1);
-    bEndo = false(0,1);
+    cid_max = max([trackInfo{:}]);
+    
+    bEcto = false(cid_max,1);
+    bSkel = false(cid_max,1);
+    bEndo = false(cid_max,1);
     
-    tids = zeros(0,1);
-    coms = zeros(0,3);
-    ctimes = zeros(0,1);
-    csizes = zeros(0,1);
-    cframes = zeros(0,1);
+    tids = zeros(cid_max,1);
+    coms = zeros(cid_max,3);
+    ctimes = zeros(cid_max,1);
+    csizes = zeros(cid_max,1);
+    cframes = zeros(cid_max,1);
     
     imSize = [size(im,1), size(im,2), size(im,3)];
     for t=1:keeplen
@@ -131,11 +140,10 @@ function track_data = collate_tracks(ROOT, filename)
     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);
+        bValidHulls = ((cframes(trackInfo{i}) > 0 ) & cframes(trackInfo{i}) <= keeplen);
+        trackInfo{i} = trackInfo{i}(bValidHulls);
     end
     
     bValidInfo = cellfun(@(x)(length(x) >= 5), trackInfo);
diff --git a/+Stats/export_intensities.m b/+Stats/export_intensities.m
new file mode 100644
index 0000000000000000000000000000000000000000..590a8e1e3fa86b3d19eb6b21289c323651ff8765
--- /dev/null
+++ b/+Stats/export_intensities.m
@@ -0,0 +1,19 @@
+function export_intensities()
+    s = load('ds_intensities.mat');
+    
+    fid = fopen('ds_intensities.csv','wt');
+    fprintf(fid, 'DatasetName,cal_ff_skel,cal_ff_ecto,cal_skel,cal_ecto,fm_ff_skel,fm_ff_ecto,fm_skel,fm_ecto\n');
+    for i=1:length(s.slist)
+        fprintf(fid, '%s,%f,%f,%f,%f,%f,%f,%f,%f\n',...
+                    s.slist(i).datasetName,...
+                    s.slist(i).means(2).ff_skel_mean,...
+                    s.slist(i).means(2).ff_ecto_mean,...
+                    s.slist(i).means(2).skel_mean,...
+                    s.slist(i).means(2).ecto_mean,...
+                    s.slist(i).means(1).ff_skel_mean,...
+                    s.slist(i).means(1).ff_ecto_mean,...
+                    s.slist(i).means(1).skel_mean,...
+                    s.slist(i).means(1).ecto_mean);
+    end
+    fclose(fid);
+end
diff --git a/+Stats/load_ectostats.m b/+Stats/load_ectostats.m
index 9b22f6d78181dd8e3a7fcc6864f8b2c6fc52a067..09fd7f2b026331232372b10e85519b5679e31564 100644
--- a/+Stats/load_ectostats.m
+++ b/+Stats/load_ectostats.m
@@ -123,14 +123,17 @@ function [size_stats, track_stats] = load_ectostats(ROOT, filename)
 
 
     %% Read per-frame vesicle info
-    bEcto = false(0,1);
-    bSkel = false(0,1);
-    bEndo = false(0,1);
+    cid_max = max([trackInfo{:}]);
     
-    coms = zeros(0,3);
-    ctimes = zeros(0,1);
-    csizes = zeros(0,1);
-    cframes = zeros(0,1);
+    bEcto = false(cid_max,1);
+    bSkel = false(cid_max,1);
+    bEndo = false(cid_max,1);
+    
+    tids = zeros(cid_max,1);
+    coms = zeros(cid_max,3);
+    ctimes = zeros(cid_max,1);
+    csizes = zeros(cid_max,1);
+    cframes = zeros(cid_max,1);
     
     imSize = [size(im,1), size(im,2), size(im,3)];
     for t=1:keeplen
@@ -175,9 +178,10 @@ function [size_stats, track_stats] = load_ectostats(ROOT, filename)
     size_stats.endo.com = coms(bEndo,:);
     
     %% Tracking stats
-    %% Skip if movie is >= 7fps
+    %% Skip if movie is >= 7fps or < 5
     medTsmp = median(diff(tsmp));
-    if ( medTsmp >= 7.0 )
+    if ( medTsmp >= 7.0 || medTsmp < 5 )
+%     if ( medTsmp >= 7.0 )
         return;
     end
     
@@ -187,11 +191,10 @@ function [size_stats, track_stats] = load_ectostats(ROOT, filename)
     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);
+        bValidHulls = ((cframes(trackInfo{i}) > 0 ) & cframes(trackInfo{i}) <= keeplen);
+        trackInfo{i} = trackInfo{i}(bValidHulls);
     end
     
     bValidInfo = cellfun(@(x)(length(x) >= 5), trackInfo);
@@ -245,12 +248,109 @@ function [size_stats, track_stats] = load_ectostats(ROOT, filename)
         %% Compute per-track info
         cids = trackInfo{i};
         
+        tframe = cframes(cids);
         ttime = ctimes(cids);
         tsize = csizes(cids);
         tpos = coms(cids,:);
         
         tdt = diff(ttime);
+        tdfr = diff(tframe);
+        
+        tdx = diff(tpos,1,1);
+        tsd = sum(tdx.^2, 2);
+        tssd = cumsum(tsd);
+        
+        n_msd = floor(length(tdfr)/2);
+        mx_msd = sum(tdfr);
         
+%         ndf = (1:mx_msd).';
+%         cdf = zeros(mx_msd,1);
+%         csd = zeros(mx_msd,1);
+%         for k=1:n_msd
+%             idxk = (1:ceil(length(tdt)/k)).';
+%             idxsum = reshape(repmat(idxk,1,k).',[],1);
+%             idxsum = idxsum(1:length(tdt));
+%             
+%             chk_df = accumarray(idxsum,tdfr);
+%             chk_sd = accumarray(idxsum,tsd);
+%             
+%             cdf = cdf + accumarray(chk_df,ones(length(chk_sd),1), [mx_msd,1]);
+%             csd = csd + accumarray(chk_df,chk_sd, [mx_msd,1]);
+%         end
+%         
+%         bvdf = (cdf > 0);
+%         ndf = ndf(bvdf);
+%         cdf = cdf(bvdf);
+%         csd = csd(bvdf);
+%         
+%         if ( length(csd) >= 6 )
+%             ndf = ndf(1:ceil(end/2));
+%             cdf = cdf(1:ceil(end/2));
+%             csd = csd(1:ceil(end/2));
+%         end
+%         
+%         % Compute Diffusion coefficient (with bias)
+%         tmsd = (csd ./ cdf);
+%         tmdt = median(tdt ./ tdfr);
+%         
+%         x1_rd = ndf * tmdt;
+%         x2_rd = ones(length(ndf),1);
+%         
+%         [b_rd,~,~,~,rstat_rd] = regress(tmsd, [x1_rd x2_rd]);
+%         D_rd = b_rd(1)/(2*msd_gamma);
+        
+        % Use single delta-t diffusoin computation if (D_rd < 0) this
+        % generally indicates a tracking error on a short track
+%         if ( D_rd < 0 )
+            x1 = cumsum(tdt);
+            x2 = ones(size(tdt,1),1);
+            [b,~,~,~,rstat] = regress(tssd, [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));
+            
+%             msde = mean(tssd ./ x1);
+%             rsqe = 1 - mean((tssd - msde*x1).^2) / var(tssd);
+%             
+%             track_msd_D(i) = msde/(2*msd_gamma);
+%             track_msd_R(i) = sqrt(rsqe);
+%         else
+%             track_msd_D(i) = D_rd;
+%             track_msd_sig(i) = b_rd(2)/(2*msd_gamma);
+%             track_msd_R(i) = sqrt(rstat_rd(1));
+%         end
+
+%         if ( length(tssd) > 30 && rstat(1) > 0.9 )
+%             xl = [0,max(x1)+1];
+%             yl = [0,max(tssd)+(max(tssd)/length(tssd))];
+%             xx = [xl(1) xl(2)+1];
+%             yy = b(1)*xx + b(2);
+%             
+%             m2 = mean(tssd ./ x1);
+%             r2 = 1 - mean((tssd - m2*x1).^2) / var(tssd);
+%             
+%             yy = m2*xx;
+%             
+%             cmap = lines(7);
+%             figure;hold on
+%             plot(x1,tssd,'.', 'Color',cmap(7,:), 'MarkerSize',12);
+%             plot(xx,yy, '-', 'Color',cmap(1,:), 'LineWidth',2);
+%             xlim(xl);
+%             ylim(yl);
+%         end
+        
+        
+%         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);
+
+        % Handle spicule stats if has defined spicule
         if ( ~isempty(sp_struct) )
             T = sp_struct.root .* pxSize;
             tds = spicule_distances(cs,Rs,rs,als, T,tpos);
@@ -264,29 +364,10 @@ function [size_stats, track_stats] = load_ectostats(ROOT, filename)
             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]);
+        [b,~,~,~,rstat] = regress(log(tssd), [x1 x2]);
         
         track_nimsd_alpha(i) = b(1);
         track_nimsd_D(i) = exp(b(2))/(2*msd_gamma);
diff --git a/+Stats/plot_ecto_stats.m b/+Stats/plot_ecto_stats.m
index de36c1418532a650bf06940ffbd772ffb8c7d04b..f67325b0d3ef559ff4521bd12ccc88035824bfe8 100644
--- a/+Stats/plot_ecto_stats.m
+++ b/+Stats/plot_ecto_stats.m
@@ -242,37 +242,37 @@ 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');
+% figure;
+% boxplot(skel_wdiridx, exp_skel_grp, 'notch','on');
+% xtickangle(30);
+% zoom('reset');
 % ylim([0.3,0.6]);
-title('Directionality Index - Endoderm');
+% 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_msdR = real(vertcat(ecto_track_info.msd_R));
 ecto_msdsig = vertcat(ecto_track_info.msd_sig);
 
 bBB94_ecto = cellfun(@(x)(strcmpi(x,'BB94')), exp_ecto_grp);
@@ -394,6 +394,7 @@ 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_diff_R = [vertcat(ecto_track_info.msd_R); vertcat(skel_track_info.msd_R)];
 
 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);
@@ -410,7 +411,7 @@ all_direct_len = [vertcat(ecto_track_info.direct_len); vertcat(skel_track_info.d
 all_path_ratio = all_direct_len ./ all_path_len;
 
 %%
-[~,ctr] = hist(all_mean_ivel(bSkelDMSO), ctr);
+[~,ctr] = hist(all_mean_ivel(bSkelDMSO), 100);
 figure;
 hist(all_mean_ivel(bSkelDMSO), ctr);
 xlim([0,0.3]);
@@ -466,11 +467,17 @@ 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)]);
+[b,~,r,~,rstDMSO] = regress(all_mean_ivel(bSkelDMSO), [ones(size(all_mean_ivel(bSkelDMSO))), all_mean_sizes(bSkelDMSO)]);
+xx = 0:0.1:1.3;
+yy = b(1) + xx*b(2);
+plot(xx,yy, '-', 'Color',cmap(2,:), 'LineWidth',2);
+
+[b,~,r,~,rstAxt] = regress(all_mean_ivel(bSkelAxt), [ones(size(all_mean_ivel(bSkelAxt))), all_mean_sizes(bSkelAxt)]);
 xx = 0:0.1:1.3;
 yy = b(1) + xx*b(2);
-plot(xx,yy, '-k', 'LineWidth',2);
-legend('Mesoderm: Control','Mesoderm: VEGFR Inhib.');
+plot(xx,yy, '-', 'Color',cmap(1,:),'LineWidth',2);
+
+legend('Skeletogenic Cells: Control','Skeletogenic Cells: VEGFR Inhib.');
 saveas(gcf,'C:\Users\mwinter\Desktop\figs\sup2_meso_speed_v_size.png');
 
 % figure;plot(all_mean_sizes(bEctoDMSO|bSkelDMSO), r, '.r');
@@ -483,10 +490,16 @@ 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)]);
+[b,~,r,~,rstDMSO] = regress(all_mean_ivel(bEctoDMSO), [ones(size(all_mean_ivel(bEctoDMSO))), all_mean_sizes(bEctoDMSO)]);
+xx = 0:0.1:1.3;
+yy = b(1) + xx*b(2);
+plot(xx,yy, '-', 'Color',cmap(2,:), 'LineWidth',2);
+
+[b,~,r,~,rstAxt] = regress(all_mean_ivel(bEctoAxt), [ones(size(all_mean_ivel(bEctoAxt))), all_mean_sizes(bEctoAxt)]);
 xx = 0:0.1:1.3;
 yy = b(1) + xx*b(2);
-plot(xx,yy, '-k', 'LineWidth',2);
+plot(xx,yy, '-', 'Color',cmap(1,:),'LineWidth',2);
+
 legend('Ectoderm: Control','Ectoderm: VEGFR Inhib.');
 saveas(gcf,'C:\Users\mwinter\Desktop\figs\sup2_ecto_speed_v_size.png');
 
@@ -500,11 +513,18 @@ 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)]);
+% [b,~,r,~,rst] = regress(all_diff_coeff(bSkelDMSO|bSkelAxt), [ones(size(all_diff_coeff(bSkelDMSO|bSkelAxt))), all_mean_sizes(bSkelDMSO|bSkelAxt)]);
+[b,~,r,~,rstDMSO] = regress(all_diff_coeff(bSkelDMSO), [ones(size(all_diff_coeff(bSkelDMSO))), all_mean_sizes(bSkelDMSO)]);
+xx = 0:0.1:1.3;
+yy = b(1) + xx*b(2);
+plot(xx,yy, '-', 'Color',cmap(2,:), 'LineWidth',2);
+
+[b,~,r,~,rstAxt] = regress(all_diff_coeff(bSkelAxt), [ones(size(all_diff_coeff(bSkelAxt))), all_mean_sizes(bSkelAxt)]);
 xx = 0:0.1:1.3;
 yy = b(1) + xx*b(2);
-plot(xx,yy, '-k', 'LineWidth',2);
-legend('Mesoderm: Control','Mesoderm: VEGFR Inhib.');
+plot(xx,yy, '-', 'Color',cmap(1,:),'LineWidth',2);
+
+legend('Skeletogenic Cells: Control','Skeletogenic Cells: VEGFR Inhib.');
 saveas(gcf,'C:\Users\mwinter\Desktop\figs\sup3_meso_diff_v_size.png');
 
 % figure;plot(all_mean_sizes(bEctoDMSO|bSkelDMSO), r, '.r');
@@ -517,10 +537,17 @@ 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)]);
+% [b,~,r,~,rst] = regress(all_diff_coeff(bEctoAxt|bEctoDMSO), [ones(size(all_diff_coeff(bEctoAxt|bEctoDMSO))), all_mean_sizes(bEctoAxt|bEctoDMSO)]);
+[b,~,r,~,rstDMSO] = regress(all_diff_coeff(bEctoDMSO), [ones(size(all_diff_coeff(bEctoDMSO))), all_mean_sizes(bEctoDMSO)]);
 xx = 0:0.1:1.3;
 yy = b(1) + xx*b(2);
-plot(xx,yy, '-k', 'LineWidth',2);
+plot(xx,yy, '-', 'Color',cmap(2,:), 'LineWidth',2);
+
+[b,~,r,~,rstAxt] = regress(all_diff_coeff(bEctoAxt), [ones(size(all_diff_coeff(bEctoAxt))), all_mean_sizes(bEctoAxt)]);
+xx = 0:0.1:1.3;
+yy = b(1) + xx*b(2);
+plot(xx,yy, '-', 'Color',cmap(1,:),'LineWidth',2);
+
 legend('Ectoderm: Control','Ectoderm: VEGFR Inhib.');
 saveas(gcf,'C:\Users\mwinter\Desktop\figs\sup3_ecto_diff_v_size.png');
 
diff --git a/+Stats/plot_vel_histograms.m b/+Stats/plot_vel_histograms.m
index 5d0dcc0661e43641befe03854b74fcd5642b79ac..55e12644081c8d525c84f45fdee4164f47f7eac4 100644
--- a/+Stats/plot_vel_histograms.m
+++ b/+Stats/plot_vel_histograms.m
@@ -69,7 +69,6 @@ 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);
@@ -78,12 +77,11 @@ 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');
+saveas(gcf, 'C:\Users\mwinter\Desktop\figs\sup_inst_speed_hist_fix.png');
 
diff --git a/+Stats/region_avg_intensity.m b/+Stats/region_avg_intensity.m
new file mode 100644
index 0000000000000000000000000000000000000000..922d488ecfcfecf92aef90d76114ee2b25dc9f93
--- /dev/null
+++ b/+Stats/region_avg_intensity.m
@@ -0,0 +1,128 @@
+function region_avg_intensity()
+    dslist = {'23-06-2016_TimeLapse1_DMSO_20hrs_Calcein_FM464'
+            '22-06-2016_TimeLapse1_DMSO_24hrs_Calcein_FM464'
+            '22-06-2016_TimeLapse3_DMSO_20hrs_Calcein_FM464'
+            '22-06-2016_TimeLapse2_DMSO_24hrs_Calcein_FM464'
+            '01-07-2016_TimeLapse4_DMSO_21hrs_Calcein_FM464'
+            '01-07-2016_TimeLapse3_DMSO_21hrs_Calcein_FM464'
+            '01-07-2016_TimeLapse5_DMSO_23hrs_Calcein_FM464'
+            '28-06-2016_TimeLapse1_Pos1_DMSO_18hrs_Calcein_FM464'
+            '28-06-2016_TimeLapse1_Pos2_DMSO_18hrs_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'
+            '23-06-2016_Timelapse_3_Axitinib_36hr_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'
+            '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'};
+
+	ROOT = 'E:\Smadar\LeverJS';
+    
+    slist = [];
+    for i=1:length(dslist)
+        leverFile = fullfile(ROOT,[dslist{i} '.LEVER']);
+        [datasetName,int_means] = get_region_intensities(leverFile);
+        if ( isempty(int_means) )
+            continue;
+        end
+        
+        s = struct('datasetName',{datasetName}, 'means',{int_means});
+        slist = [slist; s];
+    end
+    
+    save('ds_intensities.mat', 'slist');
+end
+
+function [datasetName,int_means] = get_region_intensities(leverFile)
+    int_means = [];
+
+    %% Load constants
+    conn = database(leverFile, '','', 'org.sqlite.JDBC', 'jdbc:sqlite:');
+    CONSTANTS = Read.getConstants(conn);
+    
+    %% Estimate parallelogram mask for first-frame image
+    datasetName = CONSTANTS.imageData.DatasetName;
+    
+    imagePath = CONSTANTS.imageData.imageDir;
+    chkimd = MicroscopeData.ReadMetadataFile(fullfile(imagePath,filesep));
+    if ( isempty(chkimd) )
+        warning(['Unable to load image data for: ' datasetName]);
+        return;
+    end
+    
+    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
+    if ( ~bHasEcto )
+        warning(['Unable to load ectoderm region info for: ' datasetName]);
+        return;
+    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);
+    
+    %% Get average intensities for first and first 50 frames
+    int_means = struct('ff_ecto_mean',{},'ff_skel_mean',{},'ecto_mean',{},'skel_mean',{});
+    
+    mxFrame = min(50, CONSTANTS.imageData.NumberOfFrames);
+    
+    ecto_m = zeros(mxFrame,size(im,4));
+    skel_m = zeros(mxFrame,size(im,4));
+    parfor t=1:min(mxFrame)
+        im = MicroscopeData.Reader(imagePath, 'timeRange',[t,t]);
+        % For simplicity we use mean of means for intensity
+        
+        ecto_mt = zeros(1,size(im,4));
+        skel_mt = zeros(1,size(im,4));
+        for c = 1:size(im,4)
+            imC = im(:,:,:,c);
+            
+            ecto_mean = mean(imC(bwEcto));
+            skel_mean = mean(imC(bwSkel));
+            
+            ecto_mt(1,c) = ecto_mean;
+            skel_mt(1,c) = skel_mean;
+        end
+        
+        ecto_m(t,:) = ecto_mt;
+        skel_m(t,:) = skel_mt;
+    end
+    
+    for c=1:size(im,4)
+        int_means(c).ff_ecto_mean = ecto_m(1,c);
+        int_means(c).ff_skel_mean = skel_m(1,c);
+        int_means(c).ecto_mean = mean(ecto_m(:,c));
+        int_means(c).skel_mean = mean(skel_m(:,c));
+    end
+end