% function NoduMaps8
%To make nodulation maps

%Load dataset lists
if ~exist('mpd','var') 
    mpd=uigetdir(pwd,['Select the main project folder']);
end
load([mpd '\dataset.mat']);

% %Determine where data for a group ends
% for s=1:size(dataset,1)-1;
%     grps(s,1)=strcmp(dataset(s,1:3),dataset(s+1,1:3));
% end
% grps=[grps;0];
PriL=[]; longestLat=[]; NodXY=[]; NodSpc=[]; NpRt=[]; NpPlane=[]; NpPlt=[]; TotRtLen=[]; TotRtLenpPlt=[];

for s=1:size(imgpth,1); %FOR EACH PLANE
    c_imgpth=imgpth(s,:);
    c_resfol=resfol(s,:);
    c_data=dataset(s,:);
    genotypeID=c_data(1:3);
    if str2num(c_data(end))>0
        disp(['Processing ' c_data '...']);
        load([c_resfol '/picinfo.mat'],'pltID','age','plane');
        load([c_resfol '/nodattchG.mat'],'noddat');
        load([c_resfol '/individrootsG.mat'],'attg','xygg');
        %         load([c_resfol '/LatRootsOrderAlongPrimart.mat'],'ordered','ranked','attdistp');
        load([c_resfol '/imageresolution.mat']);
        load([c_resfol '/primaryS.mat'],'prixyS');
        %lat=figure; imshow(imread(c_imgpth)); hold on
        
        %%Add primary root to all roots data
        prix=flipud(prixyS(:,1));
        priy=flipud(prixyS(:,2));
        xygg=[flipud(prixyS) zeros(length(prixyS),1);xygg];
        
        for rtID=0:max(xygg(:,3)); %% FOR EACH ROOT
            
            %%Coordinates of this root
            rtInd=find(xygg(:,3)==rtID);
            rxy=xygg(rtInd,1:2); %coordinates of this root
            rxy(isnan(rxy(:,2)),:)=[];
            if rtID>0 %add the attachment point
                rxy=[attg(rtID,:);rxy];
            end
            
            
            %%ROOT LENGTH
            dRt=0; %Distance of each point along the root
            for j=2:length(rxy)
                dRt(j,1)=dRt(j-1)+distance(rxy(j,1:2),rxy(j-1,1:2));
            end
            RtLength=dRt(end)*res; %Rt length
            %text(attg(k,1),attg(k,2)-50,num2str(ranked(k)),'color','w')
            %scatter(rxy(1:5:end,1), rxy(1:5:end,2),ones(size(rxy(1:5:end,2)))*3,RtLength+zeros(size(rxy(1:5:end,2))),'filled');
            caxis([0 20])
            
            
            %%DISTANCE TO ROOT ATTACHMENT POINT
            if rtID==0
                dpri=dRt;
                distpri=nan;
                
                %%RECORDING PRIMARY ROOT LENGTH
                if str2num(c_data(end))==3
                    PriLength=RtLength;
                    %T={'PlantID', 'Age',  'Length'};
                    PriL=[PriL; pltID age PriLength];
                end
                
            elseif rtID>0
                %locate which point of the primary corresponds to the attachmment point
                %i=find(round(prix)==attg(rtID,1)&round(priy)==attg(rtID,2));
                i=find(round(prix)==round(attg(rtID,1)) & round(priy)== round(attg(rtID,2)));
                %record the distance from the primary root top
                distpri=dpri(i)*res;
            end
            
            
            %%NUMBER OF NODULES PER ROOT
            %%Find all nodules on this root
            nodInd=find(noddat(:,end)==rtID);
            nxy=noddat(nodInd,3:4);
            NodpRt=length(nodInd);%%Number of nodules on this root
            
            if ~isempty(nodInd) %If there are nodules on the root
                
                %% NODULE DISTANCE ALONG ROOT
                attdistlat_all=[];
                for j=1:size(nxy,1)
                    i=find(rxy(:,1)==nxy(j,1)&rxy(:,2)==nxy(j,2)); %locate which point of the root corresponds to the nodule attachmment point
                    attdistlat=dRt(i)*res;%record distance of each nodule along the lateral root
                    
                    %% RECORD NODULE POSITIONS
                    %HeadingsNodXY={'PlantID','Age','PlaneID','RootID','RootDistAlongPri','NodDistAlongLateral','PlaneRank'};
                    if rtID>0
                        NodXY=[NodXY; pltID age plane rtID distpri attdistlat nan];
                    else
                        NodXY=[NodXY; pltID age plane rtID attdistlat 0 nan];
                    end
                    attdistlat_all=[attdistlat_all;attdistlat];
                end
                
                
                %% Get proper order of nodules and compute SPACING between
                ordNodD=sort(attdistlat_all);
                for j=1:length(ordNodD)-1
                    midpt=mean(ordNodD(j:j+1));
                    spc=ordNodD(j+1)-ordNodD(j);
                    % HeadingsNodSpc={'1-PlantID','2-Age','3-PlaneID','4-RootID','5-RootLength','6-RootDistAlongPri','7-NodDistAlongLateral','8-NoduleSpacing','9-PlaneRank'};
                    NodSpc=[NodSpc; pltID age plane rtID RtLength distpri midpt spc nan];
                end
                
                
                %% LATERAL RANGE of nodules per root
                NodLatRange=max(attdistlat_all)-min(attdistlat_all);
                NodInner=min(attdistlat_all);
                NodOuter=max(attdistlat_all);
            else
                NodLatRange=nan;
                NodInner=nan;
                NodOuter=nan;
            end
            
            %% Record data per root
            % HeadingsNpRt={'1-PlantID','2-Age','3-PlaneID','4-RootID' '5-RootLength', '6-RootDistAlongPri', '7-NodulePerRoot', '8-NoduleLateralRange', '9-PlaneRank', 'Innermost','Outermost'};
            NpRt=[NpRt; pltID age plane rtID RtLength distpri NodpRt NodLatRange nan NodInner NodOuter];
            
        end %END FOR EACH ROOT
        
        %% FOR EACH PLANE
        thisPltPlane=NpRt(:,1)==pltID & NpRt(:,3)==plane;
        
        %% Re-record proper ORDER of roots along primary
        [B,ordered]=sort(NpRt(thisPltPlane,6));
        [B,ranked]=sort(ordered);
        [v indPri]=max(ranked);
        ranked(indPri)=0;
        NpRt(thisPltPlane,4)=ranked;
        %save([c_resfol '/LatRootsOrderAlongPrimart.mat'],'ordered','attdistp','ranked');
        
        %% Record longest lateral root length
        %longestLat={'PlantID','Age','PlaneID','RootLength'};
        longestLat=[longestLat;pltID age plane max(NpRt(NpRt(:,1)==pltID & NpRt(:,3)==plane & NpRt(:,4)>0,5)) nan];
        
        %% Get NODULES PER PLANE
        %HeadingsNpPlane={'PlantID','Age','Plane','NodulesPerPlane','LateralNodulesPerPlane','PlaneRank'};
        thisPlaneTot=sum(NpRt(thisPltPlane,7));
        hi=max(NpRt(thisPltPlane,6));
        lo=min(NpRt(thisPltPlane,6));
        inner=min(NpRt(thisPltPlane,10));
        outer=max(NpRt(thisPltPlane,11));
        thisPlaneLat=sum(NpRt(thisPltPlane & NpRt(:,4)>0,7));
        NpPlane=[NpPlane; pltID age plane thisPlaneTot thisPlaneLat nan hi lo hi-lo inner outer outer-inner (outer-inner)/longestLat(end,4)];
        
        %% TOTAL ROOT LENGTH PER PLANE
        %HeadingsTotRtLen={'PlantID','Age','Plane','TotalLateralRootLength','PlaneRank'};
        TotRtLen=[TotRtLen; pltID age plane sum(NpRt(thisPltPlane&NpRt(:,4)>0,5)) nan];
    end
    
    if str2num(c_data(end))==3 % if on last plane
        
        %% TOTAL NODULES PER PLANT
        thisPltTot=sum(NpRt(NpRt(:,1)==pltID,7));
        thisPltLat=sum(NpRt(NpRt(:,1)==pltID & NpRt(:,4)>0,7));
        thisPltPri=sum(NpRt(NpRt(:,1)==pltID & NpRt(:,4)==0,7));
        HeadingsNpPlt={'PlantID','Age','TotalNodules','LateralNodules','PrimarNodules'};
        NpPlt=[NpPlt; pltID age thisPltTot thisPltLat thisPltPri];
        
        %% TOTAL ROOT LENGTH PER PLANT
        %HeadingsTotRtLenpPlt={'PlantID','Age','TotalRootLength'};
        TotRtLenpPlt=[TotRtLenpPlt; pltID age sum(NpRt(NpRt(:,1)==pltID & NpRt(:,4)>0,5))];
        
        %% PRIMARY ROOT NODULES
        % HeadingsNpRt={'1-PlantID','2-Age','3-PlaneID','4-RootID' '5-RootLength', '6-RootDistAlongPri', '7-NodulePerRoot', '8-NoduleLateralRange', '9-PlaneRank'};
        indPri=find(NpRt(:,1)==pltID & NpRt(:,4)==0);
        nodPri=sum(NpRt(indPri,7));
        if nodPri>0
            NpRt=[NpRt;pltID age 0 0 PriLength NaN nodPri max(NpRt(indPri,8)) nan nan nan];
        else
            NpRt=[NpRt;pltID age 0 0 PriLength NaN nodPri 0 nan nan nan];
        end
        NpRt(indPri,:)=[];
        
        %% Label planes based on nodule count rank
        indplt=find(NpPlane(:,1)==pltID & NpPlane(:,3)>0); %Find indices of plant's planes 1-3 (not 0, the primary)
        for k=1:3
            [v i]=max(NpPlane(indplt,4)); %i=which plane has the most nodules
            OriginalPlaneID=NpPlane(indplt(i),3);
            
            r=find(TotRtLen(:,1)==pltID & TotRtLen(:,3)==OriginalPlaneID);
            TotRtLen(r,end)=k;
            
            r=find(NpPlane(:,1)==pltID & NpPlane(:,3)==OriginalPlaneID);
            NpPlane(r,6)=k;
            
            r=find(NodXY(:,1)==pltID & NodXY(:,3)==OriginalPlaneID);
            NodXY(r,end)=k;
            
            r=find(NodSpc(:,1)==pltID & NodSpc(:,3)==OriginalPlaneID);
            NodSpc(r,end)=k;
            
            r=find(NpRt(:,1)==pltID & NpRt(:,3)==OriginalPlaneID);
            NpRt(r,9)=k;
            
            indplt(i)=[];
        end
    end
end


%% Save data to excel and mat files
if ~exist([mpd '\Excel files']);
    mkdir([mpd '\Excel files']);
    %mkdir([mainfol '\Figures\Nodule spacing maps']);
    %mkdir([mainfol '\Figures\Nodule distribution maps']);
end; 
pathexcel=[mpd '\Excel files'];

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

cd(pathexcel);
%plantfile- OK1 - NOTE: change into a RSE file.
%add:longest lateral root, vert range and lat range for the RSE, ADD Number
%of roots
N='WholePlantData';
H={ 'PlantID', 'Age','LengthIryRoot','TotalRootLength','NbNod','NbLatNod','NbIryNod'};
M=cat(2,PriL,TotRtLenpPlt(:,3),NpPlt(:,3:5));
xlswrite(N,H,'Sheet1','A1');
xlswrite(N,M,'Sheet1','A2')


%planefile - OK1
% N='PlaneData';
% H={'PlantID','Age','PlaneID','PlaneRank','TotalLateralRootLength','LgestLatRootLen','NodulesPerPlane','LateralNodulesPerPlane','LowestNodule','HighestNodule','VerticalRange', 'Innermost','Outermost','HorizontalRange','HNZR'};
% M=cat(2,TotRtLen(:,1:3),TotRtLen(:,5),TotRtLen(:,4),longestLat(:,4),NpPlane(:,[4:5,7:13]));
% xlswrite(N,H,'Sheet1','A1');
% xlswrite(N,M,'Sheet1','A2')

%nodulefile
N='NoduleCoordinates';
H={'PlantID','Age','PlaneID','RootID','DistAlongPri','DistAlongLateral','PlaneRank'};
M=NodXY;
xlswrite(N,H,'Sheet1','A1');
xlswrite(N,M,'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');

%rootfile
N='RootData';
H={'PlantID','Age','PlaneID','RootID' 'RootLength', 'RootDistAlongPri', 'NodulePerRoot', 'NoduleLateralRange', 'PlaneRank', 'Innermost','Outermost'};
M=NpRt;
xlswrite(N,H,'Sheet1','A1');
xlswrite(N,M,'Sheet1','A2')

%root file?? CHECK - NO: nodule as pairs of nodules and pos and dist
%between pairs (like 'segments' for which we would have seg os and length)
% N='NoduleSpacing';
% H={'PlantID','Age','PlaneID','RootID','RootLength','RootDistAlongPri','DistAlongLateral','NoduleSpacing','PlaneRank'};
% M=NodSpc;
% xlswrite(N,H,'Sheet1','A1');
% xlswrite(N,M,'Sheet1','A2')
% xlswrite(N,{'DistAlongLateral: midpoint of the distancs along lateral root for the pair of nodules'},'Notes','A1');

optmov=0;%note: set to 1 to generate movies of 3D hulls
ExtractMoreData4_2D(optmov,pathexcel,pathfigs);
NodPlotResults17(pathexcel,pathfigs);