Figure_YMMT_Substracted.m

function Figure_YMMT_Substracted
% Create figures illustrating the OR, optic tract, FA on the these tracts,
% and possibly the LGN position
%
% Repository dependencies
%    VISTASOFT
%    AFQ
%    LHON2
%
% SO Vista lab, 2014

%% Identify the directories and subject types in the study
% The full call can be
[homeDir,subDir] = Tama_subj;
[homeDir3,subDir3] = Tama_subj3;

%% load fiber groups (fg) and ROI files% select the subject{i}

% This selects a specific
whichSubject = 16; % MT =16, YM = 17;

% These directories are where we keep the data at Stanford.
% The pointers must be directed to any site.
SubDir=fullfile(homeDir,subDir{whichSubject});

% dt6 files
Pre_dt6 =fullfile(SubDir,'/dwi_2nd/dt6.mat');
switch whichSubject
    case 16
        Post_dt6 =fullfile(homeDir3,subDir3{11},'/dwi_2nd/dt6.mat');
    case 17
        Post_dt6 =fullfile(homeDir3,subDir3{13},'/dwi_2nd/dt6.mat');
end


% Load fiber groups

% define fg names
Fg = {'LOTD3L2_1206.pdb','ROTD3L2_1206.pdb','LOR1206_D4L4.pdb','ROR1206_D4L4.pdb'};
% Optic radiation
fgrOR = fgRead(fullfile(SubDir,'dwi_2nd/fibers',Fg{4}));
fglOR = fgRead(fullfile(SubDir,'dwi_2nd/fibers',Fg{3}));

% Optic tract
fglOT = fgRead(fullfile(SubDir,'dwi_2nd/fibers',Fg{1}));
fgrOT = fgRead(fullfile(SubDir,'dwi_2nd/fibers',Fg{2}));

FG = {fgrOR , fglOR, fglOT, fgrOT};

% Load dt6
PreDt = dtiLoadDt6(Pre_dt6);
PostDt = dtiLoadDt6(Post_dt6);
t1 = niftiRead(PostDt.files.t1);

% load ROIs for Figure 3, but not Figure 4
% dirROI = fullfile(SubDir,'dwi_2nd','ROIs');
% roiList = {'Optic-Chiasm.mat','Rt-LGN4.mat','Lt-LGN4.mat','lh_V1_smooth3mm_NOT.mat','rh_V1_smooth3mm_NOT.mat'};

%% draw visual pathway figure using AFQ_render
% Fig4 A

mrvNewGraphWin;hold on;
% OR
AFQ_RenderFibers(fgrOR, 'newfig', [0],'numfibers', 1000 ,'color', [0.854,0.65,0.125],'radius',[0.5,2]); %fg() = fg
AFQ_RenderFibers(fglOR, 'newfig', [0],'numfibers', 1000 ,'color', [0.854,0.65,0.125],'radius',[0.5,2]); %fg() = fg

% Occipital callosal fibers
% AFQ_RenderFibers(fg3, 'newfig', [0],'numfibers', 100 ,'color', [0.4, 0.7, 0.7],'radius',[0.5,2]); %fg() = fg

%Optic Tract
AFQ_RenderFibers(fglOT, 'newfig', [0],'numfibers', 50 ,'color', [0.67,0.27,0.51],'radius',[0.5,2]); %fg() = fg
AFQ_RenderFibers(fgrOT, 'newfig', [0],'numfibers', 50 ,'color', [0.67,0.27,0.51],'radius',[0.5,2]); %fg() = fg

% T1w
AFQ_AddImageTo3dPlot(t1, [0, 0, -30]);

% adjust figure
view(0 ,89);
set(gcf,'Color',[1 1 1])
set(gca,'Color',[1 1 1])
axis image
axis off
camlight('headlight');
title('Optic tract and optic radiation')
hold off;

%% Fig4 B
% Overlay the FA for this subject on the OR and OT

% Set the value for overlay to 'fa'
if(~exist('valName','var') || isempty(valName))
    valName = 'fa'; % 'fa','md','ad', and 'rd' are available. Change value range
end

mrvNewGraphWin;hold on;
for kk = 1:length(FG)
    % Get the FA values from the fiber group
    PostVals = dtiGetValFromFibers(PostDt.dt6,FG{kk},inv(PostDt.xformToAcpc),valName);
    PreVals = dtiGetValFromFibers(PreDt.dt6,FG{kk},inv(PreDt.xformToAcpc),valName);
    
    % substract Post from Pre values
    vals = PreVals;
    for ii =1:length(PreVals)
        k = PostVals{1,ii}(:)-PreVals{1,ii}(:);
        vals{1,ii}(:)  = k;
    end
    
    rgb = vals2colormap(vals);
    switch kk
        case {1,2}
            AFQ_RenderFibers(FG{kk},'color',rgb,'crange',[-0.2 0.2],'newfig',0,'numfibers',1000);
        case {3,4}
            AFQ_RenderFibers(FG{kk},'color',rgb,'crange',[-0.2 0.2],'newfig',0);
    end
end

% Put T1w
t1 = niftiRead(PostDt.files.t1);
AFQ_AddImageTo3dPlot(t1, [0, 0, -30]);

axis image, axis off, view(0,89)
cb = colorbar('location','eastoutside');
caxis([-0.2 0.2])
T = get(cb,'Title'); set(T,'string','FA','FontSize',14);
title('Substracted Post - Pre')

%% Fig4 C
% Fa values averaged along tracts

mrvNewGraphWin; hold on;

for kk = 1:length(FG)
    % TractProfile along current tract
    Pre_TP = SO_FiberValsInTractProfiles(FG{kk},PreDt,'AP',100,1);
    Post_TP = SO_FiberValsInTractProfiles(FG{kk},PostDt,'AP',100,1);
    
    % params
    radius = 3;
    subdivs = 20;
    crange = [-0.2 0.2];
    cmap    = 'jet';
    newfig = 0;
    
    % render core fiber with averaged FA values 
    AFQ_RenderTractProfile(Pre_TP.coords.acpc, radius, Post_TP.vals.fa - Pre_TP.vals.fa,...
        subdivs, cmap, crange, newfig)
end
% Put T1w
AFQ_AddImageTo3dPlot(t1, [0, 0, -30]);

axis image, axis off, view(0,89)
camlight('headlight');
title('Changed FA value')
% colorbar('off')
colorbar('location','eastoutside');
caxis(crange)

%% End