% function PotentialNodules(mpd)
%To extract potential nodule coordinates through automatic image processing

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

for s=1:size(imgpth,1);
    c_data=dataset(s,:);
    c_imgpth=imgpth(s,:);
    c_resfol=resfol(s,:);
    if str2num(c_data(end))>0
        if ~exist([c_resfol '/nodules_auto.mat'])
            disp(['Processing ' c_data '...']);
            load([c_resfol '/skel.mat'],'skel','mask');
            load([c_resfol '/primary.mat'],'pripath','prirc');
            load([c_resfol '/latroots.mat'],'pskel','attm')
            im=imread(c_imgpth);
            
            %% Convert to hsv
            i2=rgb2hsv(im);
            h=i2(:,:,1);
            s=i2(:,:,2);
            v=i2(:,:,3);
            
            %% Find potential nodules based on protrusions
            %     se = strel('disk',3); %Get rid of fine roots
            %     mask3=imerode(mask2,ones(3));
            %     skel2=bwmorph(mask3,'thin','Inf');
            
            eptsm=bwmorph(pskel,'endpoints'); %Map of end points (i.e root tips, or nodules)
            bpts=bwmorph(pskel-attm,'branchpoints'); %Map of branching points
            %     [yb xb]=find(pskel==1);
            %     figure(fo); hold on; plot(xb,yb,'.g')
            %     [yb xb]=find(apts==1);
            %     figure(fo); hold on; plot(xb,yb,'oy')
            %     [yb xb]=find(bpts==1);
            %     figure(fo); hold on; plot(xb,yb,'.r')
            %     [yb xb]=find(epts==1);
            %     figure(fo); hold on; plot(xb,yb,'ob')
            
            skelseg=pskel-bpts; %segment the roots and nodules
            eptsm2=bwmorph(skelseg,'endpoints'); %Need to segment more for 8-neighbour conncomp
            skelseg=skelseg-eptsm2+eptsm;
            cc=bwconncomp(skelseg,8);
            nodmap=zeros(size(skel));
            nx=[]; ny=[];
            for k=1:cc.NumObjects
                
                if numel(cc.PixelIdxList{k})<size(skel,1)*.01 & numel(cc.PixelIdxList{k})>1%If it's a short segment
                    %                         [pny pnx]=find(epts(cc.PixelIdxList{k})==1);
                    %                         if size(yn,1)>1
                    if sum(sum(eptsm(cc.PixelIdxList{k})))==1 %If it has an end point
                        %                         tempmask=zeros(size(skelseg));
                        %                         tempmask(cc.PixelIdxList{k})=1;
                        %                         se = strel('disk', 8);
                        %                         tempmask=imdilate(tempmask,se);%+nodmask;
                        
                        if median(s(cc.PixelIdxList{k}))>.3 & median(h(cc.PixelIdxList{k}))<.2 & median(v(cc.PixelIdxList{k}))>.3 %If it has the h/s/v of a nodule
                            %                         if median(double(r(cc.PixelIdxList{k}))./double(b(cc.PixelIdxList{k})))>1.25
                            %                         if mean(s(tempmask==1))>.3 & mean(h(tempmask==1))<.2 & mean(v(tempmask==1))>.3 & mean(mask(tempmask==1))>0.8 %If it has the h/s/v of a nodule
                            tmask=zeros(size(nodmap));
                            tmask(cc.PixelIdxList{k})=1;
                            
                            %% Check in rounded area if it still has nodule colour (a thin root will be surrounded by black)
                            se = strel('disk', 3);
                            tmask=imdilate(tmask,se);%+nodmask;
                            [r c]=find(tmask==1);
                            linInd=sub2ind(size(skel),r,c);
                            if mean(s(linInd))>.3 & mean(h(linInd))<.2 & mean(v(linInd))>.3 & mean(mask(linInd))>.8%If it has the h/s/v of a nodule
                                %                         if median(double(r(cc.PixelIdxList{k}))./double(b(cc.PixelIdxList{k})))>1.25
                                %                         if mean(s(tempmask==1))>.3 & mean(h(tempmask==1))<.2 & mean(v(tempmask==1))>.3 & mean(mask(tempmask==1))>0.8 %If it has the h/s/v of a nodule
                                nx=cat(1,nx,mean(c));
                                ny=cat(1,ny,mean(r));
                                nodmap(cc.PixelIdxList{k})=1;
                            end
                            
                            %                             [row col]=ind2sub(size(skelseg),cc.PixelIdxList{k});
                            %                             nx=cat(1,nx,mean(col));
                            %                             ny=cat(1,ny,mean(row));
                            %                             nodh=cat(1,nodh, mean(h(cc.PixelIdxList{k})));
                            %                             nods=cat(1,nods, mean(s(cc.PixelIdxList{k})));
                            %                             nodv=cat(1,nodv, mean(v(cc.PixelIdxList{k})));
                            %                         end
                        end
                    end
                end
            end
            %     [yb xb]=find(nodmap==1);
            %     figure(fo); hold on; plot(xb,yb,'.m')
            
            
            %% Check other regions for potential nodule that was not protruding enough to be detected
            %     se = strel('disk', 5);
            %     nodmap=imdilate(nodmap,se);%+nodmask;
            %             cc=bwconncomp(nodmap2,4);
            %     ny=[]; ny=[];
            %%Get colour channels
            r=double(im(:,:,1));
            b=double(im(:,:,2));
            nodcol=zeros(size(nodmap));
            nodcol(r./b>1.25&v>.65&s>.2)=1;
            for j=1:length(prirc)   %%Remove area above primary
                nodcol(1:prirc(j,1),prirc(j,2))=0;
            end
            %         figure;imshow(nodcol)
            %     se = strel('disk',1);
            %     nodcol=imdilate(nodcol,se);
            %      figure;imshow(nodcol)
            %     nodcol=imerode(nodcol,se);
            %      figure;imshow(nodcol)
            
            %%Smooth it out
            nodcol=bwmorph(nodcol,'majority',15);
            %         figure;imshow(nodcol)
            %                     nodcol(s>mean(nods)&h<mean(nodh)&v>mean(nodv))=1;
            cc=bwconncomp(nodcol,8);
            for k=1:cc.NumObjects
                tmask=zeros(size(nodmap));
                tmask(cc.PixelIdxList{k})=1;
                stats=regionprops(tmask,'Area','Centroid','MajorAxisLength','MinorAxisLength');
                
                if stats.Area>prod(size(mask))*0.000025 &  stats.Area<prod(size(mask))*0.0001 %If big enough, but not too big
                    %                         figure; imshow(tmask); title([num2str(stats.MajorAxisLength) ' ' num2str(stats.MinorAxisLength) ' '  num2str(stats.MajorAxisLength/stats.MinorAxisLength)])
                    
                    if stats.MajorAxisLength/stats.MinorAxisLength<2.5 %If relatively round
                        %                 if sum(nodmap(cc.PixelIdxList{k}))==0 %If a nodule hasn't already been detected there
                        %                 [r c]=ind2sub(size(nodmap),cc.PixelIdxList{k});
                        nx=cat(1,nx,stats.Centroid(1));
                        ny=cat(1,ny,stats.Centroid(2));
                        nodmap(cc.PixelIdxList{k})=1;
                        %                 end
                        %             end
                    end
                end
            end
            figure; imshow(im); hold on;
            [yb xb]=find(skelseg==1);
            plot(nx,ny,'.b')
            [yb xb]=find(nodmap==1);
            plot(nx,ny,'.m')
            
            save([c_resfol '/nodules_auto.mat'],'nx','ny','nodmap');
        end
    end
end
