%% MRI 2-exp fit for everything and then apply criteria to delete certain components 

result_2exp_replaced = allResults_1exp.porkbelly_legSensorROI.shorterTE; 

% Create cell array with 999's in case it's a bad fitting and all values need to be replaced with this
replacementArray = 999; % if fit did not work for whatever reason, replace all fit values with 999 (except exitflag)
relaxStdLimit = 1000; 

fn_slices = fieldnames(result_2exp_replaced); 

for i=1:length(fn_slices)
    tempData = result_2exp_replaced.(fn_slices{i}); 
    
    for jj=1:height(tempData)
         
            if tempData.relax1std(jj) > relaxStdLimit
                tempData.relax1(jj) = replacementArray; 
                tempData.relax1std(jj) = replacementArray; 
                tempData.amp1(jj) = replacementArray;
                tempData.amp1std(jj) = replacementArray;
            end
            
            if tempData.relax2std(jj) > relaxStdLimit
                tempData.relax2(jj) = replacementArray; 
                tempData.relax2std(jj) = replacementArray;
                tempData.amp2(jj) = replacementArray;
                tempData.amp2std(jj) = replacementArray;
            end
            
            tempData.ampRatio(jj) = tempData.amp1(jj) ./ tempData.amp2(jj); 
%             
%             if (tempData.ampRatio(jj) > 100) | (tempData.ampRatio(jj) < 1/100)
%                 isAmp1Larger = (tempData.amp1(jj) >= tempData.amp2(jj));
%                 if isAmp1Larger == 1
%                     tempData.amp2(jj) = replacementArray;
%                     tempData.amp2std(jj) = replacementArray;
%                     tempData.relax2(jj) = replacementArray; 
%                     tempData.relax2std(jj) = replacementArray; 
%                     tempData.amp2deleted(jj) = 1; 
%                 end
%                 if isAmp1Larger == 0 
%                     tempData.amp1(jj) = replacementArray; 
%                     tempData.amp1std(jj) = replacementArray; 
%                     tempData.relax1(jj) = replacementArray; 
%                     tempData.relax1std(jj) = replacementArray; 
%                     tempData.amp1deleted(jj) = 1; 
%                 end
%                 clear isAmp1Larger
%             end
    end
    
    result_2exp_replaced.(fn_slices{i}) = tempData; 
end

%% Relax Times: aggregate all relax1 and relax2 from all slices
data = vertcat(result_2exp_replaced.slice1.relax1(:), result_2exp_replaced.slice1.relax2(:), result_2exp_replaced.slice2.relax1(:), result_2exp_replaced.slice2.relax2(:), result_2exp_replaced.slice3.relax1(:), result_2exp_replaced.slice3.relax2(:), result_2exp_replaced.slice4.relax1(:), result_2exp_replaced.slice4.relax2(:)); 
data(data==999) = []; %delete any 999 values
figure; histogram(data,100); xlabel('T2 Relaxation Time (ms)'); ylabel('Pixel Counts')
title({'Relaxation Times from 2-exp Fitting','bin edges not fixed'});   set(gcf, 'color', 'w'); set(gca, 'FontSize', 16);

% Slices
binEdges = [0:5:200];
data_slice1 = vertcat(result_2exp_replaced.slice1.relax1(:), result_2exp_replaced.slice1.relax2(:));
data_slice2 = vertcat(result_2exp_replaced.slice2.relax1(:), result_2exp_replaced.slice2.relax2(:));
data_slice3 = vertcat(result_2exp_replaced.slice3.relax1(:), result_2exp_replaced.slice3.relax2(:));
data_slice4 = vertcat(result_2exp_replaced.slice4.relax1(:), result_2exp_replaced.slice4.relax2(:));
data_slice1(data_slice1==999)=[]; data_slice2(data_slice2==999)=[]; data_slice3(data_slice3==999)=[]; data_slice4(data_slice4==999)=[]; 
figure; histogram(data_slice1,binEdges); hold on; histogram(data_slice2,binEdges); histogram(data_slice3,binEdges); histogram(data_slice4,binEdges)
legend('Slice1', 'Slice 2', 'Slice 3','Slice 4'); legend boxoff
xlabel('T2 Relaxation Time (ms)'); ylabel('Pixel Counts')
title('Relaxation Times from 2-exp Fitting');   set(gcf, 'color', 'w'); set(gca, 'FontSize', 16);

% Relax1 vs Relax2
binEdges = [0:10:500];
data_relax1 = vertcat(result_2exp_replaced.slice1.relax1(:), result_2exp_replaced.slice2.relax1(:), result_2exp_replaced.slice3.relax1(:), result_2exp_replaced.slice4.relax1(:));
data_relax2 = vertcat(result_2exp_replaced.slice1.relax2(:), result_2exp_replaced.slice2.relax2(:), result_2exp_replaced.slice3.relax2(:), result_2exp_replaced.slice4.relax2(:));
data_relax1(data_relax1==999)=[];data_relax2(data_relax2==999)=[];
figure; histogram(data_relax1,binEdges); hold on; histogram(data_relax2,binEdges)
legend('Relax 1', 'Relax 2'); legend boxoff
xlabel('T2 Relaxation Time (ms)'); ylabel('Pixel Counts')
title('Relaxation Times from 2-exp Fitting');   set(gcf, 'color', 'w'); set(gca, 'FontSize', 16);


% 1exp vs 2expA vs 2expB
binEdges = [0:10:500];
data_relax1 = vertcat(result_2exp_replaced.slice1.relax1(:), result_2exp_replaced.slice2.relax1(:), result_2exp_replaced.slice3.relax1(:), result_2exp_replaced.slice4.relax1(:));
data_relax2 = vertcat(result_2exp_replaced.slice1.relax2(:), result_2exp_replaced.slice2.relax2(:), result_2exp_replaced.slice3.relax2(:), result_2exp_replaced.slice4.relax2(:));
data_1exp = data_relax1(data_relax2==999); 
data_2expA = data_relax1(data_relax2~=999); 
data_2expB = data_relax2; 
data_1exp(data_1exp==999)=[]; data_2expA(data_2expA==999)=[];data_2expB(data_2expB==999)=[];

figure;  histogram(data_2expA,binEdges); hold on; histogram(data_2expB,binEdges); histogram(data_1exp,binEdges); 
legend(sprintf('1Exp - n=%d', length(data_1exp)), sprintf('2Exp - Relax A, n=%d', length(data_2expA)), sprintf('2Exp - Relax B, n=%d', length(data_2expB))); legend boxoff
xlabel('T2 Relaxation Time (ms)'); ylabel('Pixel Counts')
title('Relaxation Times from 2-exp Fitting');   set(gcf, 'color', 'w'); set(gca, 'FontSize', 16);


%% Relax Stdevs
data_std = vertcat(result_2exp_replaced.slice1.relax1std(:), result_2exp_replaced.slice1.relax2std(:), result_2exp_replaced.slice2.relax1std(:), result_2exp_replaced.slice2.relax2std(:), result_2exp_replaced.slice3.relax1std(:), result_2exp_replaced.slice3.relax2std(:), result_2exp_replaced.slice4.relax1std(:), result_2exp_replaced.slice4.relax2std(:)); 
data_std(data_std==999)=[]; 
figure; histogram(data_std,100); xlabel('95% CI of Relax Times (ms)'); ylabel('Pixel Counts'); title('95% CI of Relaxation Times from 2-exp Fitting');  set(gcf, 'color', 'w'); set(gca, 'FontSize', 16);

data_std_subset = data_std; data_std_subset(data_std_subset > 1000) = []; %remove stdevs > 1e6
figure; histogram(data_std_subset,1000); 
xlabel('95% CI of Relax Times (ms)'); ylabel('Pixel Counts')
title({' 95% CI of Relaxation Times from 2-exp Fitting','Subset w/ bin counts = 1000'});  set(gcf, 'color', 'w'); set(gca, 'FontSize', 16);


% Slices
binEdgesStds = [0:5:200];
data_std_slice1 = vertcat(result_2exp_replaced.slice1.relax1std(:), result_2exp_replaced.slice1.relax2std(:))
data_std_slice2 = vertcat(result_2exp_replaced.slice2.relax1std(:), result_2exp_replaced.slice2.relax2std(:))
data_std_slice3 = vertcat(result_2exp_replaced.slice3.relax1std(:), result_2exp_replaced.slice3.relax2std(:))
data_std_slice4 = vertcat(result_2exp_replaced.slice4.relax1std(:), result_2exp_replaced.slice4.relax2std(:))
data_std_slice1(data_std_slice1==999)=[]; data_std_slice2(data_std_slice2==999)=[]; data_std_slice3(data_std_slice3==999)=[]; data_std_slice4(data_std_slice4==999)=[]; 
figure; histogram(data_std_slice1,binEdgesStds); hold on; histogram(data_std_slice2,binEdgesStds); histogram(data_std_slice3,binEdgesStds); histogram(data_std_slice4,binEdgesStds)
legend('Slice 1', 'Slice 2', 'Slice 3','Slice 4'); legend boxoff
xlabel('95% CI of Relax Times (ms)'); ylabel('Pixel Counts')
title({'95% CI of Relaxation Times from 2-exp Fitting','Subset w/ fixed bin edges'});  set(gcf, 'color', 'w'); set(gca, 'FontSize', 16);


%% Amplitude: aggregate all amp1 and amp2 from all slices
data_amp = vertcat(result_2exp_replaced.slice1.amp1(:), result_2exp_replaced.slice1.amp2(:), result_2exp_replaced.slice2.amp1(:), result_2exp_replaced.slice2.amp2(:), result_2exp_replaced.slice3.amp1(:), result_2exp_replaced.slice3.amp2(:),result_2exp_replaced.slice4.amp1(:), result_2exp_replaced.slice4.amp2(:)); 
data_amp(data_amp==999) = []; %delete any 999 values
figure; histogram(data_amp); xlabel('Amplitude (a.u.)'); ylabel('Pixel Counts')
title({'Amplitudes from 2-exp Fitting','bin edges not fixed'});   set(gcf, 'color', 'w'); set(gca, 'FontSize', 16);

% Slices
binEdgesAmp = [0:100:5000];
data_slice1 = vertcat(result_2exp_replaced.slice1.amp1(:), result_2exp_replaced.slice1.amp2(:))
data_slice2 = vertcat(result_2exp_replaced.slice2.amp1(:), result_2exp_replaced.slice2.amp2(:))
data_slice3 = vertcat(result_2exp_replaced.slice3.amp1(:), result_2exp_replaced.slice3.amp2(:))
% data_slice3a = vertcat(result_2exp_replaced.slice3.amp1(:))
% data_slice3b = vertcat(result_2exp_replaced.slice3.amp2(:))
data_slice4 = vertcat(result_2exp_replaced.slice4.amp1(:), result_2exp_replaced.slice4.amp2(:))
% data_slice3a(data_slice3a==999)=[]; data_slice3b(data_slice3b==999)=[];
data_slice1(data_slice1==999)=[]; data_slice2(data_slice2==999)=[]; data_slice3(data_slice3==999)=[]; data_slice4(data_slice4==999)=[]; 
figure; histogram(data_slice1,binEdgesAmp); hold on; 
histogram(data_slice2,binEdgesAmp);
% figure; histogram(data_slice3a,binEdgesAmp); hold on; histogram(data_slice3b,binEdgesAmp); legend('amp1', 'amp2')
histogram(data_slice3,binEdgesAmp); 
histogram(data_slice4,binEdgesAmp)
legend('Slice1', 'Slice 2', 'Slice 3','Slice 4'); legend boxoff
xlabel('Amplitude (a.u.)'); ylabel('Pixel Counts')
title('Amplitudes from 2-exp Fitting');   set(gcf, 'color', 'w'); set(gca, 'FontSize', 16);



%% Amplitude Std
data_amp_std = vertcat(result_2exp_replaced.slice1.amp1std(:), result_2exp_replaced.slice1.amp2std(:), result_2exp_replaced.slice2.amp1std(:), result_2exp_replaced.slice2.amp2std(:), result_2exp_replaced.slice3.amp1std(:), result_2exp_replaced.slice3.amp2std(:),result_2exp_replaced.slice4.amp1std(:), result_2exp_replaced.slice4.amp2std(:)); 
data_amp_std(data_amp_std==999) = []; %delete any 999 values
figure; histogram(data_amp_std,500); xlabel('Amplitude Std (a.u.)'); ylabel('Pixel Counts')
title({'Amplitude Std from 2-exp Fitting','default binning'});   set(gcf, 'color', 'w'); set(gca, 'FontSize', 16);

% Slices
binEdgesAmpStd = [0:50:2000];
data_slice1_std = vertcat(result_2exp_replaced.slice1.amp1std(:), result_2exp_replaced.slice1.amp2std(:))
data_slice2_std = vertcat(result_2exp_replaced.slice2.amp1std(:), result_2exp_replaced.slice2.amp2std(:))
data_slice3_std = vertcat(result_2exp_replaced.slice3.amp1std(:), result_2exp_replaced.slice3.amp2std(:))
% data_slice3a = vertcat(result_2exp_replaced.slice3.amp1(:))
% data_slice3b = vertcat(result_2exp_replaced.slice3.amp2(:))
data_slice4_std = vertcat(result_2exp_replaced.slice4.amp1std(:), result_2exp_replaced.slice4.amp2std(:))
% data_slice3a(data_slice3a==999)=[]; data_slice3b(data_slice3b==999)=[];
data_slice1_std(data_slice1_std==999)=[]; data_slice2_std(data_slice2_std==999)=[]; data_slice3_std(data_slice3_std==999)=[]; data_slice4_std(data_slice4_std==999)=[]; 
figure; histogram(data_slice1_std,binEdgesAmpStd); hold on; 
histogram(data_slice2_std,binEdgesAmpStd);
% figure; histogram(data_slice3a,binEdgesAmp); hold on; histogram(data_slice3b,binEdgesAmp); legend('amp1', 'amp2')
histogram(data_slice3_std,binEdgesAmpStd); 
histogram(data_slice4_std,binEdgesAmpStd)
legend('Slice1', 'Slice 2', 'Slice 3','Slice 4'); legend boxoff
xlabel('Amplitude Std (a.u.)'); ylabel('Pixel Counts')
title('Amplitude Std from 2-exp Fitting');   set(gcf, 'color', 'w'); set(gca, 'FontSize', 16);


%% Generate Decay from Pixel-by-Pixel Fitting Results
%  INPUT: result_2exp_replaced with the bad components already deleted out

fn_slices = fieldnames(result_2exp_replaced);  
time = 8:8:8*32; 
for ii=1:length(fn_slices)
        
    tempTable = result_2exp_replaced.(fn_slices{ii}); 
    tempTable = standardizeMissing(tempTable, 999); % Replace 999's with NaN
    
    for jj=1:height(tempTable)
        if isnan(tempTable.amp2(jj)) | isnan(tempTable.relax2(jj)) %if 2nd component is NaN, then do a 1-exp fit 
%             tempDecay(jj,:) = tempTable.amp1(jj).*exp(-time./tempTable.relax1(jj));
            tempDecay(jj,:) = exp(-time./tempTable.relax1(jj));

        else
%             tempDecay(jj,:) = tempTable.amp1(jj).*exp(-time./tempTable.relax1(jj)) + tempTable.amp2(jj).*exp(-time./tempTable.relax2(jj));
            tempDecay(jj,:) = 0.5.*exp(-time./tempTable.relax1(jj)) + 0.5.*exp(-time./tempTable.relax2(jj));

        end
    end
    
    decay(ii,:) = nansum(tempDecay,1); 
end 
avgDecay = nansum(decay,1); 


figure; plot(mri_ROI_means.legSensorROI.shorterTE.slice1(2:end)./max(mri_ROI_means.legSensorROI.shorterTE.slice1(2:end)), 'k-o'); 
hold on; plot(mri_ROI_means.legSensorROI.shorterTE.slice2(2:end)./max(mri_ROI_means.legSensorROI.shorterTE.slice2(2:end)), 'k-s'); 
plot(mri_ROI_means.legSensorROI.shorterTE.slice3(2:end)./max(mri_ROI_means.legSensorROI.shorterTE.slice3(2:end)), 'k-+');
plot(mri_ROI_means.legSensorROI.shorterTE.slice4(2:end)./max(mri_ROI_means.legSensorROI.shorterTE.slice4(2:end)), 'k-^');

hold on; plot(decay(1,:)./max(decay(1,:)), 'r-o'); plot(decay(2,:)./max(decay(2,:)), 'r-s'); plot(decay(3,:)./max(decay(3,:)), 'r-+'); plot(decay(4,:)./max(decay(4,:)), 'r-^')
hold on; plot(decay(1,2:end)./max(decay(1,2:end)), 'r-o'); plot(decay(2,2:end)./max(decay(2,2:end)), 'r-s'); plot(decay(3,2:end)./max(decay(3,2:end)), 'r-+'); plot(decay(4,2:end)./max(decay(4,2:end)), 'r-^')


% Add lines to histogram
hold on; 



line([19.44 19.44], [0 25], 'LineWidth',2, 'Color', 'b') %mri whole ROI avg 2exp relax 1
line([109 109], [0 25], 'LineWidth',2, 'Color', 'b') %mri whole ROI avg 2exp relax2

line([25.78 25.78], [0 25], 'LineWidth',2, 'Color', 'k') %leg sensor AM wholePorkbelly 0 dummy 2exp fit
line([98 98], [0 25], 'LineWidth',2, 'Color', 'k') %leg sensor AM wholePorkbelly 0 dummy 2exp fit

line([26 26], [0 25], 'LineWidth',2, 'Color', 'r') %leg sensor AM wholePorkbelly 3 dummy 2exp fit
line([114.1 114.1], [0 25], 'LineWidth',2, 'Color', 'r') %leg sensor AM wholePorkbelly 3 dummy 2exp fit

line([16.44 16.44], [0 25], 'LineWidth',2, 'Color', 'c') %MRI weighted with Bo_p0_avg
line([89.8 89.3], [0 25], 'LineWidth',2, 'Color', 'c') %

line([17.72 17.72], [0 25], 'LineWidth',2, 'Color', 'g') %MRI weighted with Bo_p0_avg
line([88.88 88.88], [0 25], 'LineWidth',2, 'Color', 'g') %

line([19.42 19.41], [0 25], 'LineWidth',2, 'Color', 'm') %MRI weighted with custom (green diamond)
line([112.5 112.5], [0 25], 'LineWidth',2, 'Color', 'm') %

%% Add Weightings to MRI Data

%
legWeightsUse = table(legSensorWeights.dist_mm, legSensorWeights.Bo_p2_avg);%USER: pick which weighting (2nd input)
legWeightsUse(isnan(legSensorWeights.Bo_p0_avg),:) = [];
legWeightsUse.Properties.VariableNames = {'dist','weight'};

lookupTable = table( [157, 156, 155, 154, 153, 152, 151]', [1,2,3,4,5,6,7]'); %USER INPUTS
lookupTable.Properties.VariableNames = {'col', 'dist'};

% ======= CUSTOM ========
lookupTable = table( [157, 156, 155, 154, 153, 152, 151]', [1,2,3,4,5,6,7]',[.4, 1,0.4,.8,.8,.8,.8]'); %USER INPUTS
lookupTable.Properties.VariableNames = {'col', 'dist','weight'};
% ========================

master = join(pixelData.porkbelly_legSensorROI.shorterTE.slice3, lookupTable);%USER: which pixelData
master = join(master, legWeightsUse);

for ii=1:height(master)
    tempDecay(ii,:) = master.weight(ii).*master.decay{ii}; 
end
decay = sum(tempDecay,1); 

figure; hold on; plot(decay, 'd-')

tempTable = allResults_2exp.porkbelly_legSensorROI.shorterTE.slice3; 
for ii=1:height(tempTable)
    imageTest_relax1(tempTable.rowID(ii), tempTable.colID(ii)) = tempTable.relax1(ii); 
    imageTest_relax2(tempTable.rowID(ii), tempTable.colID(ii)) = tempTable.relax2(ii);
end

%select just column 
master = pixelData.porkbelly_legSensorROI.shorterTE.slice3; 
rows = (master.col ==157)
vars = {'row', 'col', 'decay'}
master = master(rows, vars);


for ii=1 %:height(master)
    decay(ii,:) = master.decay{ii}; 
end
decay = mean(decay,1); 

%% ILT the weighted decay
weightedDecay1 = decay; 
time = [8:8:8*32]'; 
fieldnames_timePts = {sprintf('col%d_row%d',master.col(ii),master.row(ii))}; 

% NNLS parameters 
num_guesses = 500; %500; %precision of T2 guesses
t2Min = 0.01; %.01; %minimum possible value for T2
t2Max = 1500; %1500;
ignore = 1; 
% Smoothing Factor
alpha_values = logspace(-5,0.4771,40); %logspace(-5,0.4771,40); %, 0.1, 0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9, 1, 2]; %[linspace(0,0.0098,50) logspace(-2,1,25)];%[0,0.0001,0.0003,0.0005,0.001,0.003,0.005, 0.007,0.01];
name = 'LAC';

% Create fieldnames_alpha (a string array) 
fieldnames_alpha = strrep(strrep(strcat('alpha_', strrep(cellstr(num2str(alpha_values(:))),'.','_')), ' ',''),'-','_'); 

tic
for i= 1 %2:(nTimePts-1); %[1, nTimePts];  
%     disp(sprintf('Analyzing time point %d of %d',i, nTimePts)) %which time point are we on? 
    amp = weightedDecay1(:);     
        for j = 1:length(alpha_values)
            alpha = alpha_values(j) %which alpha value are we on? 

            % Add ILT data for this time point and this alpha value
            %ilt(j).data(i) = struct(sensorTimeString{i},'hi')
            ilt.(fieldnames_timePts{i}).(fieldnames_alpha{j}) = struct('alpha',alpha);  % OR sensorTimeString{i} IF I DON'T WANT THE 'HD' IN THERE
                       
            % Ignore initial points
            tempAmp = amp(1+ignore:end,i); 
            tempTime = time(1+ignore:end,i); 

            % Construct array of values to guess for T2
            xSpec = logspace(log10(t2Min),log10(t2Max),num_guesses);
            k = exp(-tempTime*(1./xSpec));
            k_2 = vertcat(k, eye(num_guesses) * alpha);
            amp_2 = vertcat(tempAmp,zeros(num_guesses,1));

            tic            
            %%%%%%%%% PERFORM THE INVERSION %%%%%%%%%%%%%

            % no smoothing factor
            %[ySpec,resnorm,residual,exitflag,output,lambda] = lsqnonneg(k,amp');

            % with smoothing factor
            [ySpec, resnorm, residual, exitflag, output, lambda] = lsqnonneg(k_2,amp_2);

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%            
            toc
            
            % Calculate normalized ySpec
            ySpec_norm = ySpec ./ trapz(xSpec, ySpec);
            % Calculate cumtrapz 
            integral_ySpec_norm = cumtrapz(xSpec, ySpec);
            % Calculate SNR of this time point
%             SNR = max(tempAmp)/std(tempAmp(end-100:end)); % SNR = (Max of T2 decay curve) / (Standard deviation of last 100 points, hopefully noise floor, in decay)
            % Optimal alpha calculations
            Lx = norm(k_2 * ySpec - amp_2); %save
            Ly = norm(ySpec); %save
            for kk =1:length(ySpec)
                fitArray(:,kk) = (ySpec(kk) .* exp(-tempTime./xSpec(kk)))';
            end
            fit = sum(fitArray,2); 
%             noise = std(tempAmp(end-100:end));
%             chiSquared = sum(((tempAmp-fit)/noise).^2); %save
            
            % Add to ILT structure
            ilt.(fieldnames_timePts{i}).(fieldnames_alpha{j}).xSpec = xSpec'; % = struct('ySpec',ySpec);             
            ilt.(fieldnames_timePts{i}).(fieldnames_alpha{j}).ySpec = ySpec; % = struct('ySpec',ySpec); 
            ilt.(fieldnames_timePts{i}).(fieldnames_alpha{j}).ySpec_norm = ySpec_norm;
            ilt.(fieldnames_timePts{i}).(fieldnames_alpha{j}).interal_ySpec_norm = integral_ySpec_norm;
            ilt.(fieldnames_timePts{i}).(fieldnames_alpha{j}).resnorm = resnorm;
%             ilt.(fieldnames_timePts{i}).(fieldnames_alpha{j}).snr = SNR;
            ilt.(fieldnames_timePts{i}).(fieldnames_alpha{j}).residual = residual; 
            ilt.(fieldnames_timePts{i}).(fieldnames_alpha{j}).exitflag = exitflag; 
            ilt.(fieldnames_timePts{i}).(fieldnames_alpha{j}).output = output; 
            ilt.(fieldnames_timePts{i}).(fieldnames_alpha{j}).lambda = lambda;
            ilt.(fieldnames_timePts{i}).(fieldnames_alpha{j}).Lx = Lx; 
            ilt.(fieldnames_timePts{i}).(fieldnames_alpha{j}).Ly = Ly; 
%             ilt.(fieldnames_timePts{i}).(fieldnames_alpha{j}).chiSquared = chiSquared; 
%             ilt.(fieldnames_timePts{i}).(fieldnames_alpha{j}).nTrials = sum((trialTE(:,1)~=0)); %trials(i);
%             ilt.(fieldnames_timePts{i}).(fieldnames_alpha{j}).te = trialTE(1,i); 
            ilt.(fieldnames_timePts{i}).(fieldnames_alpha{j}).ignore = ignore; % # points to ignore at start of T2 decay of data
%             ilt.(fieldnames_timePts{i}).(fieldnames_alpha{j}).nDummy = nDummies; 
%             ilt.(fieldnames_timePts{i}).(fieldnames_alpha{j}).stdev = median(stdev(:,i)); % median of stdev between all trials in this time point            
%             ilt.(fieldnames_timePts{i}).(fieldnames_alpha{j}).sensor = whichSensor; 
            ilt.(fieldnames_timePts{i}).(fieldnames_alpha{j}).user = name;
            ilt.(fieldnames_timePts{i}).(fieldnames_alpha{j}).date = datestr(now, 'yyyy-mm-dd HH:MM:SS'); 

        end 
end
toc 
 

%% Visualize ILT of the weighted decay
whichScan1 = 'porkbelly_legSensorROI_slice3'; % MRI AM
whichScan2 = 'customWeight1'; % MRI PM
whichScanLeg1 = 'porkbellyfat_whole_skinDown_cpmg_3dummy'; % Leg 1
% whichScanLeg2 = 'porkbellyfat_whole_skinDown_cpmg_3dummy'; % Leg 2

ilt1 = mri_ROI_ilt; % Original MRI 
ilt_leg1 = leg_ilt; %Leg Sensor
ilt2 = ilt; % Weighted MRI 

fieldnames_trial1 = fieldnames(ilt1.(whichScan1));
fieldnames_trial_leg1 = fieldnames(ilt_leg1.(whichScanLeg1));
fieldnames_trial_leg1(end) = []; 
figure
for j= 1:length(fieldnames_trial1)
    plot(ilt1.(whichScan1).(fieldnames_trial1{j}).xSpec, ilt1.(whichScan1).(fieldnames_trial1{j}).ySpec_norm, 'c'); hold on; 
    plot(ilt2.(whichScan2).(fieldnames_trial1{j}).xSpec, ilt2.(whichScan2).(fieldnames_trial1{j}).ySpec_norm, 'k'); 
    try 
        plot(ilt_leg1.(whichScanLeg1).(fieldnames_trial_leg1{j}).xSpec, ilt_leg1.(whichScanLeg1).(fieldnames_trial_leg1{j}).ySpec_norm,'r');
% %         plot(ilt_leg1.(whichScanLeg2).(fieldnames_trial_leg1{j}).xSpec, ilt_leg1.(whichScanLeg2).(fieldnames_trial_leg1{j}).ySpec_norm,'r-.');
%         alpha_values_leg1(j) = ilt_leg1.(whichScanLeg1).(fieldnames_trial_leg1{j}).alpha; 
%         alpha_values_leg2(j) = ilt_leg1.(whichScanLeg2).(fieldnames_trial_leg1{j}).alpha; 
    catch
    end
end
legend('MRI Original', 'MRI Weighted w/ CustomWeight3', 'Leg Sensor'); legend boxoff
%legend('cyan MRI Original', 'blue MRI col 151','red Leg Sensor','magenta MRI col 152','black MRI col 153','yellow MRI col 154','green MRI col 155','blue .. MRI col 156','magenta .. MRI col 157'); legend boxoff
set(gcf, 'color', 'w'); set(gca, 'fontsize',16) %, 'xscale', 'log') 
xlabel('T2 Relaxation Time (ms)'); 
title({sprintf('MRI ILT %s',whichScan1),'vs',sprintf('Leg Sensor %s',whichScanLeg1)},'Interpreter', 'none')


clear peaks1 allPeakLocations1 allPeakNormalizedAUCs1 peaks2 allPeakLocations2 allPeakNormalizedAUCs2 peaks_leg1 allPeakLocations_leg1 allPeakNormalizedAUCs_leg1 peaks_leg2 allPeakLocations_leg2 allPeakNormalizedAUCs_leg2
% Find Peaks & Calculate AUCs
doILTPlot = 'N'; 
for i=1:length(fieldnames_trial1)
        peaks1.(fieldnames_trial1{i}) = peakDetectionIntegration(ilt1.(whichScan1).(fieldnames_trial1{i}).xSpec(:), ilt1.(whichScan1).(fieldnames_trial1{i}).ySpec(:), ilt1.(whichScan1).(fieldnames_trial1{i}).ySpec_norm(:), doILTPlot);
        peaks2.(fieldnames_trial1{i}) = peakDetectionIntegration(ilt2.(whichScan2).(fieldnames_trial1{i}).xSpec(:), ilt2.(whichScan2).(fieldnames_trial1{i}).ySpec(:), ilt2.(whichScan2).(fieldnames_trial1{i}).ySpec_norm(:), doILTPlot);
        try
            peaks_leg1.(fieldnames_trial_leg1{i}) = peakDetectionIntegration(ilt_leg1.(whichScanLeg1).(fieldnames_trial_leg1{i}).xSpec(:), ilt_leg1.(whichScanLeg1).(fieldnames_trial_leg1{i}).ySpec(:), ilt_leg1.(whichScanLeg1).(fieldnames_trial_leg1{i}).ySpec_norm(:), doILTPlot);
%             peaks_leg2.(fieldnames_trial_leg1{i}) = peakDetectionIntegration(ilt_leg1.(whichScanLeg2).(fieldnames_trial_leg1{i}).xSpec(:), ilt_leg1.(whichScanLeg2).(fieldnames_trial_leg1{i}).ySpec(:), ilt_leg1.(whichScanLeg2).(fieldnames_trial_leg1{i}).ySpec_norm(:), doILTPlot);
        catch
        end
end


% Extract peaks, areas, etc. 
for kk=1:length(fieldnames_trial1)
    % MRI 1
    tempPeaks1 = peaks1.(fieldnames_trial1{kk});
    allPeakLocations1(kk,1:height(tempPeaks1)) = tempPeaks1.peakLocation; 
    allPeakNormalizedAUCs1(kk,1:height(tempPeaks1)) = tempPeaks1.peakAUC_percentage; 
    allPeakLocations1(allPeakLocations1==0)=nan; %replace any 0's with NaN
    allPeakNormalizedAUCs1(allPeakNormalizedAUCs1==0)=nan; 

    % MRI 2
    tempPeaks2 = peaks2.(fieldnames_trial1{kk});
    allPeakLocations2(kk,1:height(tempPeaks2)) = tempPeaks2.peakLocation; 
    allPeakNormalizedAUCs2(kk,1:height(tempPeaks2)) = tempPeaks2.peakAUC_percentage; 
    allPeakLocations2(allPeakLocations2==0)=nan; %replace any 0's with NaN
    allPeakNormalizedAUCs2(allPeakNormalizedAUCs2==0)=nan; 
    
    try 
    %Leg 1
    tempPeaks_leg1 = peaks_leg1.(fieldnames_trial_leg1{kk});
    allPeakLocations_leg1(kk,1:height(tempPeaks_leg1)) = tempPeaks_leg1.peakLocation; 
    allPeakNormalizedAUCs_leg1(kk,1:height(tempPeaks_leg1)) = tempPeaks_leg1.peakAUC_percentage; 
    allPeakLocations_leg1(allPeakLocations_leg1==0)=nan; %replace any 0's with NaN
    allPeakNormalizedAUCs_leg1(allPeakNormalizedAUCs_leg1==0)=nan; 
    
%     %Leg 2
%     tempPeaks_leg2 = peaks_leg2.(fieldnames_trial_leg1{kk});
%     allPeakLocations_leg2(kk,1:height(tempPeaks_leg2)) = tempPeaks_leg2.peakLocation; 
%     allPeakNormalizedAUCs_leg2(kk,1:height(tempPeaks_leg2)) = tempPeaks_leg2.peakAUC_percentage; 
%     allPeakLocations_leg2(allPeakLocations_leg2==0)=nan; %replace any 0's with NaN
%     allPeakNormalizedAUCs_leg2(allPeakNormalizedAUCs_leg2==0)=nan; 
    catch 
    end

end


% Visualize all Peaks/Parameters vs. Smoothing Factor
alpha_values1 = logspace(-5,0.4771,40);
alpha_values_leg1 = alpha_values_leg1 + .000001; % if smoothing factor is 0, then scatter plot w/ alpha as size of marker doesn't work
figure
for ii=1:max([ size(allPeakLocations1,2), size(allPeakLocations_leg1,2)]) %max([ size(allPeakLocations1,2), size(allPeakLocations2,2), size(allPeakLocations_leg1,2), size(allPeakLocations_leg2,2)])
    try
        scatter(allPeakLocations1(:,ii),allPeakNormalizedAUCs1(:,ii),alpha_values1.*500,'MarkerEdgeColor','c'); hold on %MRI 1
        scatter(allPeakLocations2(:,ii),allPeakNormalizedAUCs2(:,ii),alpha_values1.*500,'MarkerEdgeColor','r'); hold on %MRI 2

    catch 
    end
    try
        scatter(allPeakLocations_leg1(:,ii),allPeakNormalizedAUCs_leg1(:,ii),alpha_values_leg1.*500,'o','MarkerEdgeColor','k'); %Leg Sensor 1
%         scatter(allPeakLocations_leg2(:,ii),allPeakNormalizedAUCs_leg2(:,ii),alpha_values_leg1.*500,'o','MarkerEdgeColor','k'); %Leg Sensor 2

    catch
    end
end
set(gcf, 'color', 'w'); set(gca, 'fontsize',16) %, 'xscale', 'log') 
xlabel('Peak Locations (ms)'); ylabel('Normalized AUC (%)')
title({sprintf('MRI ILT %s',whichScan1),'vs',sprintf('Leg Sensor %s',whichScanLeg1)},'Interpreter', 'none')
legend('MRI Original','MRI Weighted', 'Leg Sensor 3 dummy'); legend boxoff

    
    