function ExtractMoreData4(optmov,pathexcel,pathfigs);
%optmov = option to make movie (set optmov=1)
%v2: correcting nodulation zone cutoff for lateral roots: not keeping 1st
%root below if nodule is on primary root, instead, record length of primary root up to lowest nodule - or actually up to bottom of nodulation zone.
if nargin<1
    optmov=0;
end;
%Load data files: nodules, roots, whole plants data

close all
% Locate data etc
% if ~exist('mpd')
%     mpd=uigetdir(pwd,['Select the excel files folder']);
% end
% cd(mpd);


% Import data on nodules

N='NoduleCoordinates';
%import_as_column_vectors([mpd '\' N '.xls'])
Mnod=xlsread([pathexcel '\' N '.xls']);
delete([pathexcel '\' N '.xls']);
PlantID=Mnod(:,1);
Age=Mnod(:,2);
PlaneID=Mnod(:,3);
RootID=Mnod(:,4);
DistAlongPri=Mnod(:,5);
DistAlongLateral=Mnod(:,6);
PlaneRank=Mnod(:,7);

%resave without plane rank
%N='NoduleCoordinates';
H={'PlantID','Age','PlaneID','RootID','DistAlongPri','DistAlongLateral'};
Mnew=Mnod(:,1:6);
xlswrite(N,H,'Sheet1','A1');
xlswrite(N,Mnew,'Sheet1','A2')
xlswrite(N,{'DistAlongPri: distance along primary root (from top) to attachment point of the root that the nodule is on; DistAlongLateral: distance along lateral root (from attachment point) where nodule formed'},'Notes','A1');


%M=[PlantID Age PlaneID RootID DistAlongPri DistAlongLateral PlaneRank]; % nodules data

N2='RootData'; % data for lateral roots recorded in planes 1-3, primary root recorded as plane 0.
M=xlsread([pathexcel '\' N2 '.xls']);
%M=[PlantID Age PlaneID RootID RootLength RootdistAlongPri NodulePerRoot Nodule lateral range Planerank Innermost Outermost]
M2=M;%to make new root file based only on nodulation zone (does not extend below).
M2=cat(2,M2,ones(size(M2,1),2));%adding two columns to tag for nodulation zone for the plant, and nodulation zone for the dataset.

%added for V2
MnodzonedepthP=zeros(size(M2,1),1);%for plant
MnodzonedepthD=zeros(size(M2,1),1);%for dataset
MnodzoneTopP=zeros(size(M2,1),1);%for plant
MnodzonetopD=zeros(size(M2,1),1);%for dataset


[vald,indxlo]=max(DistAlongPri);%val = distance of lowest nodule for the whole dataset
[valdhi,indxhi]=min(DistAlongPri);%val = distance of highest nodule for the whole dataset
%tag lateral roots lower than lower nodule for the dataset
indx2=find(M2(:,6)>vald & M2(:,3)>0);
M2(indx2,end)=0;%roots below nodulation zone for dataset are tagged with 0

%THIS BELOW REMOVED FOR V2
% if RootID(indxlo)==0
%     Mtmp=cat(2,M2,(1:size(M2,1))');
%     tmp=Mtmp(indx2,:);
%     tmp=sortrows(tmp,6);%sorting root by distance from top of primary - the highest (1st) one should be kept
%     indxkeep=tmp(1,end);
%     M2(indxkeep,end)=1;
%     clear Mtmp tmp
% end 


Plts=unique(PlantID);
%for each plant, find lowest nodule
for i=1:length(Plts)
    indx=find(PlantID==Plts(i));%get current plant
    %find lowest nodule for the plant
    [valphi,indxphi]=min(DistAlongPri(indx));%indxphi = index of highest nodule, valphi = distance of highest nodule for the plant
    [val,indxlo]=max(DistAlongPri(indx));%indxlo = index of lowest nodule, val = distance of lowest nodule for the plant
    %tag lateral roots lower than lower nodule for the plant
    indx2=find(M2(:,1)==Plts(i) & M2(:,6)>val & M2(:,3)>0);
    M2(indx2,end-1)=0;%roots below nodulation zone for current plant are tagged with 0

%THIS BELOW REMOVED FOR V2
%     %find if lowest nodule is primary root nodule, if so keep root just
%     %below as tag 1
%     if RootID(indx(indxlo))==0 % in nodule file, primary root labeled as 0
%         Mtmp=cat(2,M2,(1:size(M2,1))');
%         tmp=Mtmp(indx2,:);
%         tmp=sortrows(tmp,6);%sorting root by distance from top of primary - the highest (1st) one should be kept
%         indxkeep=tmp(1,end);
%         M2(indxkeep,end-1)=1;
%         clear Mtmp tmp
%     end

%changed to:
indx=find(M2(:,1)==Plts(i) & M2(:,3)==0);%get current plant, primary root
MnodzonedepthP(indx)=val;%depth of nod zone for the plant
MnodzonedepthD(indx)=vald;%depth of nod zone for dataset
MnodzonehiP(indx)=valphi;%depth of nod zone for the plant
MnodzonehiD(indx)=valdhi;%depth of nod zone for dataset
    
    
end


%RESAVE ROOT DATA with added columns - remove plane rank as not used in
%software
%rootfile
Ns='RootData';
M2(:,9)=[];
Hs={'PlantID','Age','PlaneID','RootID' 'RootLength', 'RootDistAlongPri', 'NodulePerRoot', 'NoduleLateralRange', 'Innermost','Outermost','PltNZ','DatasetNZ'};
xlswrite(Ns,Hs,'Sheet1','A1');
xlswrite(Ns,M2,'Sheet1','A2')


%Prepare variables to save to whole plant file
SumL1=zeros(length(Plts),1);
NbL1=zeros(length(Plts),1);
SumL2=zeros(length(Plts),1);
NbL2=zeros(length(Plts),1);
NbLN1=zeros(length(Plts),1);
NbLN2=zeros(length(Plts),1);
NodYlo=zeros(length(Plts),1);
NodYhi=zeros(length(Plts),1);
NodXmin=zeros(length(Plts),1);
NodXmax=zeros(length(Plts),1);
NodvertRange=zeros(length(Plts),1);
NodhorzRange=zeros(length(Plts),1);
NodHull2=zeros(length(Plts),1);
RootHull21=zeros(length(Plts),1);
RootHull22=zeros(length(Plts),1);
RootHull23=zeros(length(Plts),1);
NodHull3=zeros(length(Plts),1);
RootHull31=zeros(length(Plts),1);
RootHull32=zeros(length(Plts),1);
RootHull33=zeros(length(Plts),1);
for i=1:length(Plts)
    
    % TO ADD:
    %- Total length (SumL1) and number (NbL1) of lateral roots within nodulation zone defined
    % per plant
    %- Total length (SumL2) and number (NbL2) of lateral roots within nodulation zone
    % defined per dataset,
    %-Number of lateral roots bearing nodules (NbLN1) 
    indx=find(M2(:,1)==Plts(i) & M2(:,end-1)==1 & M2(:,3)>0);%get lateral roots for current plant that are in nodulation zone defined per plant
    SumL1(i)=sum(M2(indx,5));
    NbL1(i)=length(indx);
    indx=find(M2(:,1)==Plts(i) & M2(:,end-1)==1 & M2(:,3)>0 & M2(:,7)>0);
    NbLN1(i)=length(indx);
    %length of lateral roots bearing nodules:
    SumLLnod(i)=sum(M2(indx,5));%sum of length of nodulated lateral roots
    %Tag if primary root is bearing noduls: 0 no, 1 yes
    indx=find(M2(:,1)==Plts(i) & M2(:,3)==0 & M2(:,7)>0);
    if isempty (indx)
        Prim_nod(i)=0;
    else
        Prim_nod(i)=1;
    end
    
    indx=find(M2(:,1)==Plts(i) & M2(:,end)==1 & M2(:,3)>0);%get lateral roots for current plant that are in nodulation zone defined per dataset
    SumL2(i)=sum(M2(indx,5));
    NbL2(i)=length(indx);
    indx=find(M2(:,1)==Plts(i) & M2(:,end)==1 & M2(:,3)>0 & M2(:,7)>0);
    
    
    %[PlantID Age PlaneID RootID DistAlongPri DistAlongLateral PlaneRank]; % nodules data
    %RSE 2D Data
    %-Y distance of Lowest nodule on plant
    %-Y distance of highest
    %-x distance min lat
    %-x distance max lat
    %-vertrange
    %-horzrange
    indx=find(PlantID==Plts(i));%get current plant
    NodYlo(i)=max(DistAlongPri(indx));
    NodYhi(i)=min(DistAlongPri(indx));
    NodXmin(i)=min(DistAlongLateral(indx));
    NodXmax(i)=max(DistAlongLateral(indx));
    NodvertRange(i)=abs(NodYlo(i)-NodYhi(i));
    NodhorzRange(i)=abs(NodXmax(i)-NodXmin(i));
    
    
    %-Hull of 2D RSE for nodules
    %NOTE: for now for hulls we are only using lateral roots, but could use
    %primary too (no sense for now it seems as primary very long and not all roots recorded)
    %-Hull of 2D RSE for lateral roots (all kept) RootHullArea1
    %-Hull of 2D RSE for lateral roots in dataset nodulation zone RootHullArea2
    %-Hull of 2D RSE for lateral roots in plant nodulation zone RootHullArea3
    
    indx=find(PlantID==Plts(i));
    Xnod=DistAlongLateral(indx);
    Ynod=DistAlongPri(indx);
    K1 = convhull(Xnod,Ynod);
    NodHull2(i)=polyarea(Xnod(K1),Ynod(K1));
    %M2=[PlantID Age PlaneID RootID RootLength RootdistAlongPri NodulePerRoot Nodule lateral range Planerank Innermost Outermost tagNZplt tagNZdataset]
    indx=find(M2(:,1)==Plts(i) & M2(:,3)>0);
    Xroot=M2(indx,5);
    Yroot=M2(indx,6);
    
    XX=Xroot;
    YY=Yroot;
    %for plot
    figure;hold on
    
    for ii=1:length(Xroot)
        h=line([0;Xroot(ii)],[Yroot(ii);Yroot(ii)]);
        set(h,'Color',[0.5 0.5 0.5]);
    end
    
    Xroot=cat(1,Xroot,zeros(size(Xroot,1),1));
    Yroot=cat(1,Yroot,Yroot);
    K2 = convhull(Xroot,Yroot);
    RootHull21(i)=polyarea(Xroot(K2),Yroot(K2));
    %h=line(Xroot(K2),Yroot(K2));set(h,'Color',[0 0 1])
    %patch(Xroot(K2),Yroot(K2),'b','FaceAlpha',.1,'EdgeColor',[0 0 1])
    
    indx=find(M2(:,1)==Plts(i) & M2(:,3)>0 & M2(:,end)==1);%find lateral roots in dataset nodulation zone
        
    Xroot=M2(indx,5);
    Yroot=M2(indx,6);
    Xroot=cat(1,Xroot,zeros(size(Xroot,1),1));
    Yroot=cat(1,Yroot,Yroot);
    
    %adding primary root up to bottonm of dataset nodulation zone
    indp1=find(M2(:,1)==Plts(i) & M2(:,3)==0);
    Xroot=cat(1,Xroot,[0;0]);
    Yroot=cat(1,Yroot,[MnodzonehiD(indp1); MnodzonedepthD(indp1)]);
    
    K4 = convhull(Xroot,Yroot);
    RootHull22(i)=polyarea(Xroot(K4),Yroot(K4));
    hh2=line(Xroot(K4),Yroot(K4));set(hh2,'Color',[0 1 0])
    %patch(Xroot(K4),Yroot(K4),'g','FaceAlpha',.1,'EdgeColor',[0 1 0])
    
    indx=find(M2(:,1)==Plts(i) & M2(:,3)>0 & M2(:,end-1)==1);%find lateral roots in plant nodulation zone
    Xroot=M2(indx,5);
    Yroot=M2(indx,6);
    Xroot=cat(1,Xroot,zeros(size(Xroot,1),1));
    Yroot=cat(1,Yroot,Yroot);
    
    %adding primary root up to bottonm of plant nodulation zone
    indp1=find(M2(:,1)==Plts(i) & M2(:,3)==0);
    Xroot=cat(1,Xroot,[0;0]);
    Yroot=cat(1,Yroot,[MnodzonehiP(indp1); MnodzonedepthP(indp1)]);
    
    K3 = convhull(Xroot,Yroot);
    RootHull23(i)=polyarea(Xroot(K3),Yroot(K3));
    hh3=line(Xroot(K3),Yroot(K3));set(hh3,'Color',[1 0 0])
    %patch(Xroot(K3),Yroot(K3),'r','FaceAlpha',.1,'EdgeColor',[1 0 0])
    
    
    h=plot(Xnod,Ynod,'.');set(h,'Color',[0 0 0]);
    h=line(Xnod(K1),Ynod(K1));set(h,'Color',[0 0 0]);
    axis ij
    axis equal
    
    %save figure
    figsave=fullfile(pathfigs,['Hull2D_RSEfig' num2str(Plts(i)) '.fig']);
    saveas(gcf,figsave);
    close gcf
    %3D Hulls: Nodules: for each plane, assign X coord, from x, Z coord from y info,
    %then assign Y using polar coordinates in XY plane
    %[PlantID Age PlaneID RootID DistAlongPri DistAlongLateral PlaneRank]; % nodules data
    
%     %-Hull of 3D RSE for nodules
%     theta=[0 120 240]*pi/180;
%     Xnodtot=[];
%     Ynodtot=[];
%     Znodtot=[];
%     for j=1:3
%         
%         indx=find(PlantID==Plts(i) & PlaneID==j);
%         Xnod=DistAlongLateral(indx);
%         Znod=DistAlongPri(indx);
%         [Xnod,Ynod]=pol2cart(ones(size(Xnod,1),1)*theta(j),Xnod);
%         Xnodtot=cat(1,Xnodtot,Xnod);
%         Ynodtot=cat(1,Ynodtot,Ynod);
%         Znodtot=cat(1,Znodtot,Znod);
%         
%     end
%     figure;
%     plot3(Xnodtot,Ynodtot,-Znodtot,'.');
%     hold on
%     
%     [K31,NodHull3(i)]=convhull(Xnodtot,Ynodtot,-Znodtot,'simplify',true);
%     h=trisurf(K31,Xnodtot,Ynodtot,-Znodtot);
%     set(h,'EdgeColor',[0.7 0 0],'EdgeAlpha',0.1,'FaceColor',[0.5 0.5 0],'FaceAlpha',0.1)
%     
%     
%     %-Hull of 3D RSE for roots (all kept)
%     [Xroottot,Yroottot,Zroottot,KR1,RootHull31(i)]=getroots3D(M2,Plts(i),MnodzonedepthD,MnodzonedepthP,MnodzonehiD,MnodzonehiP,1);
%     h=trisurf(KR1,Xroottot,Yroottot,-Zroottot);
%     set(h,'EdgeColor',[0 0 0.7],'EdgeAlpha',0.1,'FaceColor',[0 0 1],'FaceAlpha',0.05)
%     
%     %-Hull of 3D RSE for roots in dataset nodulation zone
%     [Xroottot,Yroottot,Zroottot,KR2,RootHull32(i)]=getroots3D(M2,Plts(i),MnodzonedepthD,MnodzonedepthP,MnodzonehiD,MnodzonehiP,2);
%     h=trisurf(KR2,Xroottot,Yroottot,-Zroottot);
%     set(h,'EdgeColor',[0 0.7 0],'EdgeAlpha',0.1,'FaceColor',[0 1 0],'FaceAlpha',0.05)
%     
%     %-Hull of 3D RSE for roots in plant nodulation zone
%     [Xroottot,Yroottot,Zroottot,KR3,RootHull33(i)]=getroots3D(M2,Plts(i),MnodzonedepthD,MnodzonedepthP,MnodzonehiD,MnodzonehiP,3);
%     h=trisurf(KR3,Xroottot,Yroottot,-Zroottot);
%     set(h,'EdgeColor',[0.7 0 0],'EdgeAlpha',0.1,'FaceColor',[1 0 0],'FaceAlpha',0.05)
%     
%     %save figure
%     figsave=fullfile(pathfigs,['Hull3fig' num2str(Plts(i)) '.fig']);
%     saveas(gcf,figsave);
%     close gcf
%     %Make movie
%     disp('')
%     if optmov==1,
%         axis vis3d
%         fna=['moviehull3PLTnew2' num2str(Plts(i))];
%         nbf=20;
%         aviobj=avifile([fna '.avi'],'FPS',nbf,'Compression','Cinepak');
%         
%        axis([-25 25 -20 20 -25 1]) 
%        for cc=1:360
%            view(cc,0)
%             axis([-25 25 -20 20 -25 1]) 
%             %camorbit(5,0,'camera')
%             drawnow
%             %set(gcf,'Position',[360 500 420 315])
%             set(gcf,'Position',[1281 1 1280 950])
% 
%             %axis off;
%            % set(gcf, 'Renderer', 'zbuffer');
%             frame=gcf;
%             aviobj=addframe(aviobj,frame);
%         end
%         aviobj=close(aviobj);
%         close all;
%     end;
    
    %Save new variables in whole plant data:
    
end

%ADDING VARIABLES TO WHOLE PLANT FILE
N3='WholePlantData'; % data for whole plants
P=xlsread([pathexcel '\' N3 '.xls']);
%Increase measures extracted at the RSE level if we are to promote RSE as a
%way to standardize analyses (derive a wide range of measures for RSE).
%P=[PlantID Age LengthIryRoot TotalRootLength NbNod NbLatNod NbIryNod];



Ns='WholePlantData';
Hs={'PlantID' 'Age' 'LengthIryRoot' 'TotalRootLength' 'NbNod' 'NbLatNod' 'NbIryNod' 'NodYlo' 'NodYhi' 'NodXmin' 'NodXmax' 'NodYrange' 'NodXrange' 'RootL_NZP' 'RootNb_NZP' 'RootL_NZD' 'RootNb_NZD' 'NbNodulated_Roots' 'NodHull2A' ...
    'RootHull21A' 'RootHull22A' 'RootHull23A' 'LengthIrydat','totLength_nodulatedlatroots','Iryhasnodules'};
P=cat(2,P,NodYlo,NodYhi,NodXmin,NodXmax,NodvertRange,NodhorzRange,SumL1,NbL1,SumL2,NbL2,NbLN1,NodHull2,RootHull21,RootHull22,RootHull23,zeros(length(Plts),1)+vald,SumLLnod',Prim_nod');
xlswrite(Ns,Hs,'Sheet1','A1');
xlswrite(Ns,P,'Sheet1','A2')

function [Xroottot,Yroottot,Zroottot,K,V]=getroots3D(M2,pltnum,MnodzonedepthD,MnodzonedepthP,MnodzonehiD,MnodzonehiP,opt)

hold on
Xroottot=[];
Yroottot=[];
Zroottot=[];

%plot primary root
indp=find(M2(:,1)==pltnum & M2(:,3)==0);
h=line([0;0],[0;0],[0 -M2(indp,5)]);
set(h,'Color',[0 0 0.8]);

% root data: M2=[PlantID Age PlaneID RootID RootLength RootdistAlongPri NodulePerRoot Nodule lateral range Planerank Innermost Outermost tagNZplt tagNZdataset]
theta=[0 120 240]*pi/180;
for i=1:3 %going through planes
    if opt==1
        indx=find(M2(:,3)==i & M2(:,1)==pltnum & M2(:,3)>0);%keeping all roots
    elseif opt==2
        indx=find(M2(:,3)==i & M2(:,1)==pltnum & M2(:,3)>0 & M2(:,end)==1);%keep roots in nodulation zone for the plant
    elseif opt==3
        indx=find(M2(:,3)==i & M2(:,1)==pltnum & M2(:,3)>0 & M2(:,end-1)==1);%keep roots in nodulation zone for dataset
        
    end
    Xroot=M2(indx,5);
    Zroot=M2(indx,6);
    [Xroot,Yroot]=pol2cart(ones(size(Xroot,1),1)*theta(i),Xroot);
    
    for ii=1:length(Xroot)
        h=line([0;Xroot(ii)],[0;Yroot(ii)],-[Zroot(ii);Zroot(ii)]);
        set(h,'Color',[0 0 0.8]);
    end
    
    Xroot=cat(1,Xroot,zeros(size(Xroot,1),1));
    Yroot=cat(1,Yroot,zeros(size(Yroot,1),1));
    Zroot=cat(1,Zroot,Zroot);
    
    Xroottot=cat(1,Xroottot,Xroot);
    Yroottot=cat(1,Yroottot,Yroot);
    Zroottot=cat(1,Zroottot,Zroot);
    
end

%adding primary root up to bottonm of dataset nodulation zone
 indp1=find(M2(:,1)==pltnum & M2(:,3)==0);
 
if opt==1
    vall=[];
elseif opt==2
   vall=MnodzonedepthD(indp1);
   vallh=MnodzonehiD(indp1);
elseif opt==3
    vall=MnodzonedepthP(indp1);
    vallh=MnodzonehiP(indp1);
end
    
if ~isempty(vall)
    Xroottot=cat(1,Xroottot,[0;0]);
    Yroottot=cat(1,Yroottot,[0;0]);
    %Zroottot=cat(1,Zroottot,[vallh;vall]);
    Zroottot=cat(1,Zroottot,[0;vall]);
else
    Xroottot=cat(1,Xroottot,[0]);
    Yroottot=cat(1,Yroottot,[0]);
    Zroottot=cat(1,Zroottot,[0]);
end

[K,V] = convhull(Xroottot,Yroottot,Zroottot);

