%function NodPlotData


%% Locate data etc
if ~exist('mpd')
    mpd=uigetdir(pwd,['Select the main project folder']);
end
if ~exist([mpd '\Excel files']);
    disp('Excel files folder does not exist')
end; %cd([mpd '\Excel files'])

if ~exist([mpd '\Figures']);
    mkdir([mpd '\Figures']);   
end; %cd([mpd '\Figures'])


%% SET FIGURE PARAMETERS
w=3; h=3;
fmts={'fig','eps','png'};
figpth=[mpd '\Figures'];
xlab='Lateral root axis (cm)';
ylab='Primary root axis (cm)';
fsz=10;

% 
% %% NODULE DISTRIBUTION PLOT
% N='NodulecoordinatesAPR';
% clearvars -except mpd N w h fsz fmts figpth xlab ylab  ax xlab ylab
% import_as_column_vectors([mpd '\Excel files\' N '.xls'])
% ax=[0  floor(max(max(DistAlongLateral))) 0 floor(max(max(DistAlongPri)))];
% figure
% circle6([DistAlongLateral DistAlongPri],.15,16,'k',1/length(DistAlongPri)*200);
% SaveFigs(w,h,fsz,xlab,ylab,ax,fmts,figpth,N)
% disp(['Total number of nodules from all plants all planes: ' num2str(length(DistAlongPri))])
% 
% T='Nodulation zones per plane';
% 
% %Note: adding hull of root system for that plane: for each plane we would
% %have a ratio of nodule occupation based on hulls
% N2='RootDataAPR';
% MatRoot=xlsread([mpd '\Excel files\' N2 '.xls']); %AG added apr 22 2013 % 
% RootL=MatRoot(:,5);
% RootDPri=MatRoot(:,6);
% 
% figure; hold on
% Plts=unique(PlantID);
% ax=[0  floor(max(max(DistAlongLateral))) 0 floor(max(max(DistAlongPri)))];
% A=[];
% for plti=1:length(Plts)
%     for planei=1:3
%         NodPlaX=DistAlongLateral(PlantID==Plts(plti)&PlaneID==planei);
%         NodPlaY=DistAlongPri(PlantID==Plts(plti)&PlaneID==planei);
%         
%         %AG
%         RPlaX=RootL(MatRoot(:,1)==Plts(plti)&MatRoot(:,3)==planei);
%         RPlaX=cat(1,RPlaX,zeros(size(RPlaX,1),1));
%         RPlaY=RootDPri(MatRoot(:,1)==Plts(plti)&MatRoot(:,3)==planei);
%         RPlaY=cat(1,RPlaY,RPlaY);
%         %
%         if length(NodPlaY)>2
%         K = convhull(NodPlaX,NodPlaY);      
%         K2 = convhull(RPlaX,RPlaY);%AG 22 apr 2013
%          patch(NodPlaX(K),NodPlaY(K),'k','FaceAlpha',.1)
%          plot(NodPlaX(K),NodPlaY(K))
%          
%          
%         hold all
%         A=cat(1,A,[polyarea(NodPlaX(K),NodPlaY(K)) polyarea(RPlaX(K2),RPlaY(K2))]);%recording hull of nodules and roots
%         end
%     end
% end
% SaveFigs(w,h,fsz,xlab,ylab,ax,fmts,figpth,T)
% disp(['Mean nodulation zone area per plane: ' num2str(mean(A(:,1)))])
% disp(['SD of nodulation zone area per plane: ' num2str(std(A(:,1)))])
% 
% T='Nodulation zones per plant';
% figure; hold on
% Plts=unique(PlantID);
% ax=[0  floor(max(max(DistAlongLateral))) 0 floor(max(max(DistAlongPri)))];
% A=[];
% col=jet(length(Plts));
% for plti=1:length(Plts)
% %     for planei=1:3
%         NodPlaX=DistAlongLateral(PlantID==Plts(plti));%&PlaneID==planei);
%         NodPlaY=DistAlongPri(PlantID==Plts(plti));%&PlaneID==planei);
%         
%         %AG
%         RPlaX=RootL(MatRoot(:,1)==Plts(plti)&MatRoot(:,3)>0);
%         RPlaX=cat(1,RPlaX,zeros(size(RPlaX,1),1));
%         RPlaY=RootDPri(MatRoot(:,1)==Plts(plti)&MatRoot(:,3)>0);
%         RPlaY=cat(1,RPlaY,RPlaY);
%         %
%         
%         if length(NodPlaY)>2
%         K = convhull(NodPlaX,NodPlaY);        
%         K2 = convhull(RPlaX,RPlaY);%AG 22 apr 2013
%          h=line(NodPlaX(K),NodPlaY(K));set(h,'Color',col(plti,:),'LineWidth',[1.5]);
%          %h=plot(NodPlaX,NodPlaY,'.');set(h,'Color',[0.5 0.5 0.5]);
%          hh=line(RPlaX(K2),RPlaY(K2));set(hh,'Color',col(plti,:),'LineWidth',[2],'LineStyle','--');
%         hold all
%         %A=[A;polyarea(NodPlaX(K),NodPlaY(K))];
%         A=cat(1,A,[polyarea(NodPlaX(K),NodPlaY(K)) polyarea(RPlaX(K2),RPlaY(K2))]);%recording hull of nodules and roots
% %         end
%     end
% end
% %SaveFigs(w,h,fsz,xlab,ylab,ax,fmts,figpth,T)
% disp(['Mean nodulation zone area per plant: ' num2str(mean(A(:,1)))])
% disp(['SD of nodulation zone area per plant: ' num2str(std(A(:,1)))])
% 
% %% MEAN NODULE COUNTS
% 
% % %%Nodules per plane REMOVED AG APR 22
% % N='Nodules per plane';
% % clearvars -except mpd N w h fsz fmts figpth xlab ylab  ax
% % import_as_column_vectors([mpd '\Excel files\' N '.xls'])
% % disp(['Mean lateral nodules per plane: ' num2str(mean(LateralNodulesPerPlane))])
% % disp(['SD of lateral nodules per plane: ' num2str(std(LateralNodulesPerPlane))])
% % disp(['Mean nodule vertical range: ' num2str(mean(VerticalRange))])
% % disp(['SD of nodule vertical range: ' num2str(std(VerticalRange))])
% 
% %%Nodules per plant
% N='Nodules per plant';
% clearvars -except mpd N w h fsz fmts figpth xlab ylab  ax
% import_as_column_vectors([mpd '\Excel files\' N '.xls'])
% disp(['Mean lateral nodules per plant: ' num2str(mean(LateralNodules))])
% disp(['SD of lateral nodules per plant ' num2str(std(LateralNodules))])
% disp(['Mean total nodules per plant: ' num2str(mean(TotalNodules))])
% disp(['SD of total nodules per plant ' num2str(std(TotalNodules))])



%% NODULE AND ROOT DENSITY MAPS
N='Nodulecoordinates';
clearvars -except mpd N w h fsz fmts figpth xlab ylab  ax
import_as_column_vectors([mpd '\Excel files\' N '.xls'])
ax=[0  floor(max(max(DistAlongLateral))) 0 floor(max(max(DistAlongPri)))];
nodRtID=RootID(RootID>0);
xpos=DistAlongLateral(RootID>0);
ypos=DistAlongPri(RootID>0);
nodPltID=PlantID(RootID>0);
w=3; h=3;

N='RootDataext';
import_as_column_vectors([mpd '\Excel files\' N '.xls'])
RootDistAlongPri=RootDistAlongPri+.00001;

inc=1; disp('Densities are computed per cm2 per plant')
xg=0:inc:ax(2)+inc;
yg=0:inc:ax(4)+inc;
NodIn=[]; RtIn=[]; NodRtRatio=[];
Plts=unique(PlantID);
for plti=1:length(Plts)
    for xi=1:length(xg)-1
        for yi=1:length(yg)-1
            NodIn(yi,xi,plti)=sum(inpolygon(xpos(nodPltID==Plts(plti)),ypos(nodPltID==Plts(plti)),[xg(xi) xg(xi+1) xg(xi+1) xg(xi)],[yg(yi) yg(yi) yg(yi+1) yg(yi+1)]));
            RtIn(yi,xi,plti)=sum(RootLength>xg(xi) & RootDistAlongPri>yg(yi) & RootDistAlongPri<=yg(yi+1) & PlantID==Plts(plti));
            %AG adding info on root length in square
            indx=find(RootLength>xg(xi) & RootDistAlongPri>yg(yi) & RootDistAlongPri<=yg(yi+1) & PlantID==Plts(plti));
            RtL(yi,xi,plti)=0;
            for cp=1:length(indx)
                if RootLength(indx(cp))>=xg(xi+1)
                    RtL(yi,xi,plti)=RtL(yi,xi,plti)+1;
                else
                    RtL(yi,xi,plti)=RtL(yi,xi,plti)+RootLength(indx(cp))-xg(xi);
                end
            end
            NodRtnbRatio(yi,xi,plti)=NodIn(yi,xi,plti)./RtIn(yi,xi,plti);
            NodpermmRoot(yi,xi,plti)=NodIn(yi,xi,plti)./RtL(yi,xi,plti);   
        end
    end
end
[xx yy]=meshgrid(0:inc:ax(2),0:inc:ax(4));

uchoosex=input(['Choose length (cm) of lateral root (x) axis for nodulation zone maps (min ' num2str(ax(2)+inc) '): ']);
uchoosey=input(['Choose length (cm) of primary root (y) axis for nodulation zone maps (min ' num2str(ax(4)+inc) '): ']);
ax1=[0 uchoosex 0 uchoosey]; 
T='Map mean Nodule density grid';
M=mean(NodIn,3);
cmax=max(max(M));

draw_gridmap(M,cmax,size(M,1),size(M,2),0);
SaveFigs2(w,h,fsz,xlab,ylab,ax1,fmts,figpth,T)

xlswrite([mpd '\Excel files\' T],xx,'Sheet1','A1');
xlswrite([mpd '\Excel files\' T],yy','Sheet1','A2');
xlswrite([mpd '\Excel files\' T],M,'Sheet1','B2');

T='Map SE Nodule density grid';
M=std(NodIn,1,3)/sqrt(size(NodIn,3));
draw_gridmap(M,cmax,size(M,1),size(M,2),0);
SaveFigs2(w,h,fsz,xlab,ylab,ax1,fmts,figpth,T)

xlswrite([mpd '\Excel files\' T],xg,'Sheet1','A1');
xlswrite([mpd '\Excel files\' T],yg','Sheet1','A2');
xlswrite([mpd '\Excel files\' T],M,'Sheet1','B2');

T='Map mean Root number per cm2 grid';
M=mean(RtIn,3);
cmax=max(max(M));
draw_gridmap(M,cmax,size(M,1),size(M,2),0);
SaveFigs2(w,h,fsz,xlab,ylab,ax1,fmts,figpth,T)

xlswrite([mpd '\Excel files\' T],xg,'Sheet1','A1');
xlswrite([mpd '\Excel files\' T],yg','Sheet1','A2');
xlswrite([mpd '\Excel files\' T],M,'Sheet1','B2');

T='Map SE Root number per cm2 grid';
M=std(RtIn,1,3)/sqrt(size(NodIn,3));
draw_gridmap(M,cmax,size(M,1),size(M,2),0);
SaveFigs2(w,h,fsz,xlab,ylab,ax1,fmts,figpth,T)

xlswrite([mpd '\Excel files\' T],xg,'Sheet1','A1');
xlswrite([mpd '\Excel files\' T],yg','Sheet1','A2');
xlswrite([mpd '\Excel files\' T],M,'Sheet1','B2');

T='Map mean Nodule-to-Root number ratio grid';
NodRtnbRatio(isnan(NodRtnbRatio))=0;
M=mean(NodRtnbRatio,3);
cmax=max(max(M));
draw_gridmap(M,cmax,size(M,1),size(M,2),0);
SaveFigs2(w,h,fsz,xlab,ylab,ax1,fmts,figpth,T)

xlswrite([mpd '\Excel files\' T],xg,'Sheet1','A1');
xlswrite([mpd '\Excel files\' T],yg','Sheet1','A2');
xlswrite([mpd '\Excel files\' T],M,'Sheet1','B2');


T='Map SE Nodule-to-Root number ratio grid';
M=std(NodRtnbRatio,1,3)/sqrt(size(NodIn,3));

draw_gridmap(M,cmax,size(M,1),size(M,2),0);
SaveFigs2(w,h,fsz,xlab,ylab,ax1,fmts,figpth,T)

xlswrite([mpd '\Excel files\' T],xg,'Sheet1','A1');
xlswrite([mpd '\Excel files\' T],yg','Sheet1','A2');
xlswrite([mpd '\Excel files\' T],M,'Sheet1','B2');


T='Map mean Root length per cm2 grid';
M=mean(RtL,3);
cmax=max(max(M));
draw_gridmap(M,cmax,size(M,1),size(M,2),0);
SaveFigs2(w,h,fsz,xlab,ylab,ax1,fmts,figpth,T)

xlswrite([mpd '\Excel files\' T],xg,'Sheet1','A1');
xlswrite([mpd '\Excel files\' T],yg','Sheet1','A2');
xlswrite([mpd '\Excel files\' T],M,'Sheet1','B2');


T='Map SE Root length per cm2 grid';
M=std(RtL,1,3)/sqrt(size(NodIn,3));
draw_gridmap(M,cmax,size(M,1),size(M,2),0);
SaveFigs2(w,h,fsz,xlab,ylab,ax1,fmts,figpth,T)

xlswrite([mpd '\Excel files\' T],xg,'Sheet1','A1');
xlswrite([mpd '\Excel files\' T],yg','Sheet1','A2');
xlswrite([mpd '\Excel files\' T],M,'Sheet1','B2');


T='Map mean nb nodule per mm of Root per cm2 grid';
NodpermmRoot(isnan(NodpermmRoot))=0;
M=mean(NodpermmRoot,3);
cmax=max(max(M));
draw_gridmap(M,cmax,size(M,1),size(M,2),0);
SaveFigs2(w,h,fsz,xlab,ylab,ax1,fmts,figpth,T)

xlswrite([mpd '\Excel files\' T],xg,'Sheet1','A1');
xlswrite([mpd '\Excel files\' T],yg','Sheet1','A2');
xlswrite([mpd '\Excel files\' T],M,'Sheet1','B2');

T='Map SE nb nodule per mm of Root per cm2 grid';
M=std(NodpermmRoot,1,3)/sqrt(size(NodIn,3));
draw_gridmap(M,cmax,size(M,1),size(M,2),0);
SaveFigs2(w,h,fsz,xlab,ylab,ax1,fmts,figpth,T)

xlswrite([mpd '\Excel files\' T],xg,'Sheet1','A1');
xlswrite([mpd '\Excel files\' T],yg','Sheet1','A2');
xlswrite([mpd '\Excel files\' T],M,'Sheet1','B2');



%% LATERAL ROOT PLOT
N='RootDataext';
clearvars -except mpd N w h fsz fmts figpth xlab ylab  ax
import_as_column_vectors([mpd '\Excel files\' N '.xls'])
RootDistAlongPri=RootDistAlongPri+.001;
w=6; h=6; ax=[0  floor(max(max(RootLength))) 0 floor(max(max(RootDistAlongPri)))];
figure; hold on
colmatch=1:max(RootLength)/64:max(RootLength);
matrgb=jet;
for k=1:length(RootLength)
    [v colind]=min(abs(RootLength(k)-colmatch));
    plot([0 RootLength(k)],[RootDistAlongPri(k) RootDistAlongPri(k)],'color', matrgb(colind,:))
end
T='Lateral root plot';
SaveFigs(w,h,fsz,xlab,ylab,ax,fmts,figpth,T)
disp(['Total number of lateral roots from all plants all planes: ' num2str(length(RootDistAlongPri(RootID>0)))])


%% MEAN ROOT MEASUREMENTS AG removed Apr 22 2013
% N='Longest lateral root per plane';
% clearvars -except mpd N w h fsz fmts figpth xlab ylab  ax
% import_as_column_vectors([mpd '\Excel files\' N '.xls'])
% disp(['Mean longest lateral root length: ' num2str(mean(RootLength))])
% disp(['SD of longest lateral root length: ' num2str(std(RootLength))])
% 
% N='Primary root lengths';
% clearvars -except mpd N w h fsz fmts figpth xlab ylab  ax
% import_as_column_vectors([mpd '\Excel files\' N '.xls'])
% disp(['Mean longest primary root length: ' num2str(mean(Length))])
% disp(['SD of longest primary root length: ' num2str(std(Length))])
% 
% N='Total root length per plane';
% clearvars -except mpd N w h fsz fmts figpth xlab ylab  ax
% import_as_column_vectors([mpd '\Excel files\' N '.xls'])
% disp(['Mean total lateral root length per plane: ' num2str(mean(TotalLateralRootLength))])
% disp(['SD of total lateral root length per plane ' num2str(std(TotalLateralRootLength))])
% 
% N='Total root length per plant';
% clearvars -except mpd N w h fsz fmts figpth xlab ylab  ax
% import_as_column_vectors([mpd '\Excel files\' N '.xls'])
% disp(['Mean total lateral root length per plant: ' num2str(mean(TotalRootLength))])
% disp(['SD of total lateral root length per plant ' num2str(std(TotalRootLength))])


%% NODULE SPACING % AG Removed Apr 22 2013
% N='Nodule spacing';
% clearvars -except mpd N w h fsz fmts figpth xlab ylab  ax
% import_as_column_vectors([mpd '\Excel files\' N '.xls'])
% ax=[0  floor(max(max(DistAlongLateral))) 0 floor(max(max(RootDistAlongPri)))];
% 
% %Scatter plot of nodule spacing
% figure
% x=DistAlongLateral;%midpoint distance of nodule pair along the lateral root
% y=RootDistAlongPri;%distance of lateral root along primary
% priInd=isnan(y);
% y(priInd)=x(priInd);
% x(priInd)=0;
% z=NoduleSpacing;%nodule spacing
% scatter(x,y,ones(size(x))*5,z,'filled')
% SaveFigs(w,h,fsz,xlab,ylab,ax,fmts,figpth,N)
% 
% %Nodule spacing map
% figure
% inc=1;
% [xg yg]=meshgrid(0:inc:ax(2),0:inc:ax(4));
% f=TriScatteredInterp(x,y,z);
% zg=f(xg,yg);
% % surf(xg,yg,zg,'EdgeAlpha',0)
% pcolor(xg,yg,zg); shading flat
% T='Nodule spacing (cm)';
% SaveFigs(w,h,fsz,xlab,ylab,ax,fmts,figpth,T)
% 
% %Nodule spacing data
% disp(['Mean nodule spacing: ' num2str(mean(NoduleSpacing))])
% disp(['SD of nodule spacing: ' num2str(std(NoduleSpacing))])







