% function SmoothLaterals
%Smooth out the lateral root coordinates and put in right order
close all
clear all
%Load dataset lists
if ~exist('mpd')
    mpd=uigetdir(pwd,['Select the main project folder']);
end; load([mpd '\dataset.mat']);

fixattach=input('Press 1 if you ran SegLaterals5 or an earlier version for this data set, 0 otherwise');

for s=1:size(imgpth,1);
    c_data=dataset(s,:);
    c_imgpth=[mpd imgpth(s,end-35:end)];
    c_resfol=[mpd resfol(s,end-19:end)];
    if str2num(c_data(end))>0
        disp(['Processing ' c_data '...']);
        load([c_resfol '/individroots.mat'],'attg','xyg','aom');
        load([c_resfol '/primary.mat'])
        load([c_resfol '/primaryS.mat'])
        im=imread(c_imgpth); figure; imshow(im); hold on; title(c_data)
        xygg=[];
        %plot(prix,priy,'.y','MarkerSize',.5)
        plot(prixyS(:,1),prixyS(:,2),'.y','MarkerSize',.5)
       %         plot(xyg(:,1),xyg(:,2),'.w','MarkerSize',.5)
        %         mask=logical(ones(size(im,2),size(im,2)));
        
        attgg=[];
        for k=1:max(xyg(:,3)); %for each individ root
            %find and plot the root
            totrt=[];
            rtID=k;
            rtInd=find(xyg(:,3)==rtID);
            txyg=xyg(rtInd,1:2);
            taom=aom(rtInd);%aom made by SegLaterals11
            %             plot(txyg(:,1),txyg(:,2),'.g') %root
            %             for j=1:length(txyg)
            %                 text(txyg(j,1),txyg(j,2),num2str(j),'color','w') %root
            %             end
            
            %correcting attachment points
            d=distance(prixyS(:,1:2),attg(k,:));%Finding nearest coordinate on the primary
            [val st]=min(d);
            attgg=[attgg;prixyS(st,1), prixyS(st,2)]; %This is the proper attachment point

%             plot(attgg(k,1),attgg(k,2),'or','markerfacecolor','r') %root attchmnt
%             text(attgg(k,1), attgg(k,2),num2str(rtID),'color','w') %root ID
            

 
            
            %aom=sum(mod(txyg,1),2)==0|sum(mod(txyg,1),2)==1;
            %Points from auto are exact, aom=1=auto, aom=0=manual
            
            %AG to Lauren: I kept this and adapted it but did not check why
            %this is needed (making a point as nan)
            if size(taom,1)>1
            if taom(1)==0 & taom(2)==1
                txyg(1,:)=nan;
            end
            else
                disp('check issue')
      
            end
            
            axy=txyg(taom==1,:); %auto points
            mxy=txyg(taom==0,:); %manual points
            
            %             if distance(txyg(2,1:2),txyg(3,1:2))<=sqrt(2) | distance(txyg(3,1:2),txyg(4,1:2))<=sqrt(2) %If they are close to each other they're from the autotrace
            
            %% DEAL WITH AUTO POINTS
            if size(axy,1)>2
                axy=round(axy);
                attgg=round(attgg);
                %Find pos closest to attchmnt
                axy(1,:)=[];
                d=distance(axy(:,1:2),attgg(k,:));%Finding nearest coordinate on the primary
                [val st]=min(d);
                
                % Find the tip (of autotrace) and trace attchmnt to tip (get rid of small stumps etc)
                linInd = sub2ind(size(pripath),axy(:,2),axy(:,1));
                rtmsk=zeros(size(pripath));
                rtmsk(linInd)=1;
                rtmsk=bwmorph(rtmsk,'bridge');
                % rtmsk=bwmorph(rtmsk,'spur');
                %                 rtmsk=imdilate(rtmsk,ones(3));
                %                 rtmsk=bwmorph(rtmsk,'skel');
                rtmsk(axy(st,2), axy(st,1))=1;
                D1=bwdistgeodesic(logical(rtmsk), axy(st,1), axy(st,2));
                D1(D1==Inf)=nan;
                rtlen=max(max(D1));
                [ytip xtip]=find(D1==rtlen);
                ytip=ytip(1); xtip=xtip(1);
                %                 plot(xtip,ytip,'*m')
                
                D2=bwdistgeodesic(logical(rtmsk),xtip,ytip);
                D=D1+D2;
                rtmsk=zeros(size(pripath));
                rtmsk(D==D(ytip,xtip))=1;
                
                rtrc=bwtraceboundary(rtmsk,[axy(st,2), axy(st,1)],'s',8);
                % rty=rtrc(:,1); rtx=rtrc(:,2); plot(rtx,rty,'r')
                
                %% put in right order
                rtrc=unique(rtrc,'rows','first');
                rty=rtrc(:,1);
                rtx=rtrc(:,2);
                xygd=[axy(st,1), axy(st,2)];
                for j=1:length(rtx);
                    [v i]=min(distance([rtx rty],xygd(end,1:2))); %Find nearest point
                    if  rtx(i)==xtip & rty(i)==ytip
                        break
                    end
                    xygd=[xygd; rtx(i) rty(i)];
                    rtx(i,:)=[]; rty(i,:)=[]; %Remove from search
                end
                %plot(xygg(:,1),xygg(:,2),'m')
                
                
                %% Smooth out
                [xys]=smoothline2(xygd,10,'median'); %remove jagged spots
                
                %% Add the attachment point and store every 4th value
                rt=[attgg(k,:); xys(1:4:end,:)];
                plot(rt(:,1),rt(:,2),'-g')
                
                totrt=[totrt;rt];

            end
            
            %% Add on any manual traces
            
            
            %% DEAL WITH MANUAL POINTS
            if size(mxy,1)>1
                mrtx=[];mrty=[];
                %                 %% spline %sometimes does weird things
                %                 npoints=round((max(mxy(:,2))-min(mxy(:,2)))/5);
                %                 t=1:length(mxy);
                %                 ts=1:length(mxy)/npoints:length(mxy);
                %                 mxy=spline(t,mxy(:,1:2)',ts);
                %                 plot(mxy(1,:),mxy(2,:),'c'); %plot(mxy(1,:),mxy(2,:),'c')
                if size(axy,1)<1
                    if fixattach==1
                        %we have to re-extract attachment point here
                        d=distance(prixyS(:,1:2),mxy(1,:));%Finding nearest coordinate on the primary
                        [val2 st2]=min(d);
                        attgg(k,1:2)=[prixyS(st2,1), prixyS(st2,2)]; %This is the proper attachment point
                        mxy(1,:)=attgg(k,1:2);
                    else
                         mxy=cat(1,attgg(k,:),mxy); %AG - adding base of manually extracted lateral root
                    end
                end
                plot(mxy(:,1),mxy(:,2),'.m');
                for m=1:length(mxy)-1
                    d=distance(mxy(m+1,:),mxy(m,:));
                    dx=(mxy(m+1,1)-mxy(m,1));
                    dy=(mxy(m+1,2)-mxy(m,2));
                    %                     b=mxy(m,2)-dy/dx*mxy(m,1);
                    %                     ax=mxy(m,1):dx/(d/5):mxy(m+1,1);
                    %                     ay=ax*dy/dx+b;
                    
                    ax=mxy(m,1)+(0:d/5)*dx/length(0:d/5);
                    ay=mxy(m,2)+(0:d/5)*dy/length(0:d/5);
                    totrt=[totrt;ax' ay'];
                    plot(ax,ay,'.c');
                end
                
                plot(mxy(:,1),mxy(:,2),'.m');
                
            end
            plot(attgg(k,1),attgg(k,2),'or','markerfacecolor','r') %root attchmnt
            text(attgg(k,1), attgg(k,2),num2str(rtID),'color','w') %root ID
            totrt(isnan(totrt(:,1)),:)=[];
            %replace first point of lateral root by exact attachment point
            %to smoothed primary root
            totrt(1,:)=[];
            totrt=cat(1,attgg(k,:),totrt);
            
            plot(totrt(1,1),totrt(1,2),'*b','markerfacecolor','b')
            plot(totrt(end,1),totrt(end,2),'om','markerfacecolor','m')
            xygg=[xygg;totrt rtID*ones(size(totrt,1),1)];
            
            
            
        end
        attg=attgg;
        save([c_resfol '/individrootsG.mat'],'attg','xygg');
    end
end



