% For in vivo spine data analyzed using MapManager V_20171028 
% Author: Richard Roth (rhroth@stanford.edu)
% Date: Sept 2023
% Based on MapManager-Matlab-1.0.1-alpha by Robert Cudmore, 20171016

% This script requires a spine pool generated from in vivo imaging data
% that was first analyzed using MapManager V_20171028 and imported to
% Matlab using MapManager-Matlab-1.0.1-alpha. This script adds custom
% analyses to this pool which are required for further analysis, such as of
% spine fluorescence intensity.

%% Add new analysis to all maps in a pool normalize to baseline per spine

%Define baseline session (Up to five, delete/comment out excessive lines)
clear base
base{1,1} = 'day-3'; % Enter name of baseline session 1 (Here: base1)
base{2,1} = 'day-2'; % Enter name of baseline session 1
base{3,1} = 'day-1'; % Enter name of baseline session 1
base{4,1} = 'day0'; % Enter name of baseline session 1
% base{5,1} = 'base5'; % Enter name of baseline session 1
numbase = size(base,1);

% These are for 1 channel images

newStatPair=[{'utsmdtdm'},{1},{'ntbutsmdtdmch1'};...
             {'utssdtds'},{1},{'ntbutssdtdsch1'};...
             {'sLen2d'},{1},{'ntbsLen2d'};...
             {'sLen3d'},{1},{'ntbsLen3d'};...
             {'sbMean'},{1},{'ntbsbMeanch1'};...
             {'sMean'},{1},{'ntbsMeanch1'};...
             ];         
         
                  
for s = 1:size(newStatPair,1)
         
numMaps=length(myMapList);

for i=1:numMaps
    
ps = mmMap.defaultPlotStruct();

ps.stat = newStatPair{s,1}; 
ps.channel  = newStatPair{s,2};
ps = myMapList(i).GetMapValues(ps);
newStatName = newStatPair{s,3};
[m,n] = size(ps.val);
newStatValues = NaN(m,n);

if numbase==1
baseline = myMapList(i).GetMapValuesCond(ps, base{1});
end
if numbase==2
baseline = nanmean([myMapList(i).GetMapValuesCond(ps, base{1}) myMapList(i).GetMapValuesCond(ps, base{2})],2);
end
if numbase==3
baseline = nanmean([myMapList(i).GetMapValuesCond(ps, base{1}) myMapList(i).GetMapValuesCond(ps, base{2}) myMapList(i).GetMapValuesCond(ps, base{3})],2);
end
if numbase==4
baseline = nanmean([myMapList(i).GetMapValuesCond(ps, base{1}) myMapList(i).GetMapValuesCond(ps, base{2}) myMapList(i).GetMapValuesCond(ps, base{3}) myMapList(i).GetMapValuesCond(ps, base{4})],2);
end
if numbase==5
baseline = nanmean([myMapList(i).GetMapValuesCond(ps, base{1}) myMapList(i).GetMapValuesCond(ps, base{2}) myMapList(i).GetMapValuesCond(ps, base{3}) myMapList(i).GetMapValuesCond(ps, base{4}) myMapList(i).GetMapValuesCond(ps, base{5})],2);
end

newStatValues = bsxfun(@rdivide, ps.val, baseline);

myMapList(i).addUserStat(newStatName, newStatValues);

clearvars ps baseline n m newStatValues 
end
end
clearvars s newStatPair 


%% Add new analysis to all maps in a pool background subtracted spine ch1 ./ background subtracted dendrite ch1 #Added Sept 2021

numMaps=length(myMapList);
for i=1:numMaps
    
ps = mmMap.defaultPlotStruct();
ps.stat = 'ubssMean'; 
% ps.stat = 'sMean';
ps.channel  = 1;
ps = myMapList(i).GetMapValues(ps);

ps2 = mmMap.defaultPlotStruct();
ps2.stat = 'ubsdMean'; 
% ps2.stat = 'dMean'; 
ps2.channel  = 1;
ps2 = myMapList(i).GetMapValues(ps2);

newStatName = 'nutsmdtdmch1';
[m,n] = size(ps.val);
newStatValues = NaN(m,n);

newStatValues = ps.val ./ ps2.val;
myMapList(i).addUserStat(newStatName, newStatValues);

clearvars ps ps2
ps = mmMap.defaultPlotStruct();
ps.stat = 'nutsmdtdmch1';
ps = myMapList(i).GetMapValues(ps);

if numbase==1
baseline = myMapList(i).GetMapValuesCond(ps, base{1});
end
if numbase==2
baseline = nanmean([myMapList(i).GetMapValuesCond(ps, base{1}) myMapList(i).GetMapValuesCond(ps, base{2})],2);
end
if numbase==3
baseline = nanmean([myMapList(i).GetMapValuesCond(ps, base{1}) myMapList(i).GetMapValuesCond(ps, base{2}) myMapList(i).GetMapValuesCond(ps, base{3})],2);
end
if numbase==4
baseline = nanmean([myMapList(i).GetMapValuesCond(ps, base{1}) myMapList(i).GetMapValuesCond(ps, base{2}) myMapList(i).GetMapValuesCond(ps, base{3}) myMapList(i).GetMapValuesCond(ps, base{4})],2);
end
if numbase==5
baseline = nanmean([myMapList(i).GetMapValuesCond(ps, base{1}) myMapList(i).GetMapValuesCond(ps, base{2}) myMapList(i).GetMapValuesCond(ps, base{3}) myMapList(i).GetMapValuesCond(ps, base{4}) myMapList(i).GetMapValuesCond(ps, base{5})],2);
end

newStatValues = bsxfun(@rdivide, ps.val, baseline);

newStatName = 'nntbutsmdtdmch1';
myMapList(i).addUserStat(newStatName, newStatValues);

% value at day minus baseline diveded by value at day plus baseline
clearvars ps
ps = mmMap.defaultPlotStruct();
ps.stat = 'nutsmdtdmch1';
ps = myMapList(i).GetMapValues(ps);

if numbase==1
baseline = myMapList(i).GetMapValuesCond(ps, base{1});
end
if numbase==2
baseline = nanmean([myMapList(i).GetMapValuesCond(ps, base{1}) myMapList(i).GetMapValuesCond(ps, base{2})],2);
end
if numbase==3
baseline = nanmean([myMapList(i).GetMapValuesCond(ps, base{1}) myMapList(i).GetMapValuesCond(ps, base{2}) myMapList(i).GetMapValuesCond(ps, base{3})],2);
end
if numbase==4
baseline = nanmean([myMapList(i).GetMapValuesCond(ps, base{1}) myMapList(i).GetMapValuesCond(ps, base{2}) myMapList(i).GetMapValuesCond(ps, base{3}) myMapList(i).GetMapValuesCond(ps, base{4})],2);
end
if numbase==5
baseline = nanmean([myMapList(i).GetMapValuesCond(ps, base{1}) myMapList(i).GetMapValuesCond(ps, base{2}) myMapList(i).GetMapValuesCond(ps, base{3}) myMapList(i).GetMapValuesCond(ps, base{4}) myMapList(i).GetMapValuesCond(ps, base{5})],2);
end

newStatValues = bsxfun(@rdivide, ps.val-baseline, ps.val+baseline);

newStatName = 'nntb2utsmdtdmch1';
myMapList(i).addUserStat(newStatName, newStatValues);

clearvars ps baseline n m newStatValues
end