function status = prostate_pathomics

% Pathomic Feature Calculator

%% ------------------------------------------------------------------------

show_results = 0;

fid = fopen(sprintf('tile_seg_output1-%s.csv',datestr(now(),'mm_dd_yyyy')),'a');

fid2 = fopen(sprintf('tile_seg_output2-%s.csv',datestr(now(),'mm_dd_yyyy')),'a');

anot = {'Seminal_Vesicles','Atrophy','HGPIN','G3','G4FG','G4CG','G5','Tissue'}; % annotation classes

scanner = {'Huron','Olympus'}; % whole slide image scanners

for s = 1:numel(scanner)

    d = dir(sprintf('%s/*/*.tiff',scanner{s}));

    date = {d.date}';

    input_images = {d.name}';

    folder = {d.folder}';

    for i = 1:numel(d)

        grade_dir = strsplit(folder{i},'/');

        grade{i} = grade_dir{end};

        grade_number(i) = find(contains(anot,grade{i}));

    end

    for i = 1:numel(d)

        grade_dir = strsplit(folder{i},'/');

        grade{i} = grade_dir{end};

        grade_number(i) = find(contains(anot,grade{i}));

    end

    for x = 1:numel(input_images)

        fprintf('Processing: %s\n',input_images{x});

        im = imread(input_images{x});

        im_resized = imresize(im,0.4);

        % Image color deconvolution

        im_resized = im_resized(:,:,1:3);

        I = rgb2gray(im_resized);

        tissue_mask = I<200;

        tissue_mask = ~bwareaopen(~tissue_mask,6000); % removes rip noise

        tissue_mask = imgaussfilt(uint8(tissue_mask),4);

        tissue_mask = uint8(bwareaopen(tissue_mask,50000)); % help flood lumen

        tissue_mask_label = bwlabel(tissue_mask);

        [lumen_mask,epithelium_cells,epithelium_mask,stroma_mask] = gimmeSegs_prostate2(im_resized,0,0,0);

        disp('Cleaning up segmentations');

        epithelium_mask = epithelium_mask .* tissue_mask;

        stroma_mask = stroma_mask .* tissue_mask;

        lumen_mask_noise = lumen_mask .* tissue_mask; % masks small lumen and some epithelium

        lumen_mask_noise = bwareafilt(logical(lumen_mask_noise),[0 150]);

        % Clean up lumen

        if strcmp(scanner{s},'Huron')

            lumen_mask = im_resized(:,:,2)>240; % huron lumen more white than oly

        end

        se = strel('disk',1);

        lumen_mask = lumen_mask & ~lumen_mask_noise;

        lumen_mask = imclose(lumen_mask,se);

        lumen_mask = imfill(lumen_mask,'holes');

        lumen_mask = bwareaopen(lumen_mask,100); %100

        lumen_mask = imgaussfilt(uint8(lumen_mask),1);

% -------------------------------------------------------------------------

        % Clean up stroma

        stroma_mask = stroma_mask & ~lumen_mask & ~epithelium_cells;

        stroma_mask = bwmorph(stroma_mask,'bridge');

        stroma_mask = ~bwareaopen(~stroma_mask,200);

        stroma_mask = bwareaopen(stroma_mask,20);

        stroma_mask = imgaussfilt(uint8(stroma_mask),1);

% -------------------------------------------------------------------------

        % Clean up epithelium

        epithelium_mask = logical(epithelium_mask) + epithelium_cells;

        if grade_number(x) < 4 | grade_number(x) == 8

            se = strel('disk',4);

            epithelium_mask = imclose(epithelium_mask,se);

        elseif grade_number(x) == 6

            se = strel('disk',1);

            epithelium_mask = imclose(epithelium_mask,se);

        end

        epithelium_mask = logical(epithelium_mask) & ~lumen_mask & ~stroma_mask;

        epithelium_mask = bwmorph(epithelium_mask,'thick',2);

        epithelium_mask = bwareaopen(epithelium_mask,250);

        epithelium_mask = ~bwareaopen(~epithelium_mask,500);

        epithelium_mask = bwareaopen(epithelium_mask,500);

        epithelium_mask = imgaussfilt(uint8(epithelium_mask),1);

        epithelium_mask_filled = imfill(epithelium_mask,'holes') & ~lumen_mask & ~stroma_mask;

        epithelium_mask_filled = ~bwareaopen(~epithelium_mask_filled,100);

        epithelium_mask_filled = bwareaopen(epithelium_mask_filled,500);

        epithelium_mask_filled = imgaussfilt(uint8(epithelium_mask_filled),1);

% -------------------------------------------------------------------------

        lumen_mask = bwmorph(lumen_mask,'thick',2);

        lumen_mask = uint8(imclearborder(lumen_mask,4)); % removes any lumen touching tile edges

        disp('Segmenting Gland');

        SEL_image = stroma_mask + 2*lumen_mask + 3*epithelium_mask_filled;

        %         figure('Position',[100 2000, 1500, 900]);

        %         subplot(131); imagesc(im_resized); axis image

        %         title(input_images{x});

        %         subplot(132); imagesc(SEL_image); axis image

        %         title('SEL_image');

        %         subplot(133); imagesc(SEL_image_epi); axis image

        %         title('SEL_image separate epi');

        %         saveas(gcf,sprintf('SEL_%s',input_images{x}),'tif');

        % Begin epith_wall_thickness_calculation, so re-map the variables into the legacy code

        bw_lumen = uint8(lumen_mask);

        bw_lumen_label = bwlabel(bw_lumen);

        bw_epith_label = bwlabel(epithelium_mask);

        bw_stroma = stroma_mask;

        bw_epith = uint8(epithelium_mask);

        bw_epith_filled = epithelium_mask_filled;

        bw_epith_filled_label = bwlabel(bw_epith_filled);

        core_mask = tissue_mask;

        bw_cells = uint8(epithelium_cells);

        %         figure('Position',[100 2000, 1500, 900]);

        %         subplot(2,3,1);

        %         imagesc(im_resized); axis image;

        %         title('Raw Image')

        %         subplot(2,3,2);

        %         imagesc(tissue_mask_label); axis image;

        %         title('Tissue Mask')

        %         subplot(2,3,3);

        %         imagesc(bw_epith);axis image;

        %         title('Epithelium Mask');

        %         subplot(2,3,4);

        %         imagesc(bw_epith_label);axis image;

        %         title('Epithelium labels');

        %         subplot(2,3,5);

        %         imagesc(bw_lumen);axis image;

        %         title('Lumen Mask')

        %         subplot(2,3,6);

        %         imagesc(bw_stroma);axis image;

        %         title('Stroma Mask')

        %         % saveas(gcf,sprintf('seg_%s',input_images{x}),'tif');

% -------------------------------------------------------------------------

        % Stroma v epithelium area calculations

        stroma_area = bwarea(bw_stroma);

        epith_area = bwarea(bw_epith);

% -------------------------------------------------------------------------

        [Bl,Ll,Nl,Al] = bwboundaries(bw_lumen_label,'noholes');

        statsl = regionprops(Ll,'Area','Centroid','PixelList');

        [Be,Le,Ne,Ae] = bwboundaries(bw_epith_filled_label,'noholes');

        statse = regionprops(Le,'Area','Centroid','PixelList');

        colors=['b' 'g' 'r' 'c' 'm' 'y'];

        %         colors = ['c' 'c' 'c' 'c' 'c' 'c'];

        ep_lb = [];

%% ------------------------------------------------------------------------

        bw_lumen_roundness = bw_lumen_label;

        bw_lumen_area = bw_lumen_label;

        bw_epith_thickness = bw_epith_filled_label;

        bw_cell_frac = bw_epith_filled_label;

        bw_epith_size = bw_epith_filled_label;

        bw_epith_roundness = bw_epith_filled_label;

        lumen_tort = [];

        epith_size = [];

        cell_frac = [];

        wall_thickness = [];

        min_lum_index_thick = [];

        wall_thick = [];

        epith_tort = [];

%         figure('Position',[100 100 1500 1500])

%         imagesc(im_resized);axis image; hold on;

        for i = 1:length(Bl)

            ep_lb(i) = bw_epith_filled_label(statsl(i).PixelList(1,2),statsl(i).PixelList(1,1)); % this maybe was backwards? 2,1?

            if ep_lb(i) == 0 %This indicates that an epithelium didn't get closed, and these are then excluded

                disp(sprintf('lumen number %i is connected with an epithelium edge and will be excluded',i));

                j = i; % in next step j is used for loop

                boundaryl = Bl{j};

                area = statsl(j).Area;

                delta_sq = diff(boundaryl).^2;

                perimeter = sum(sqrt(sum(delta_sq,2)));

                lumen_tort(j) = 4*pi*area/(perimeter^2);

                bw_lumen_roundness(bw_lumen_roundness == j) = lumen_tort(j); % replaces labeled lumen with tortuosity value

                bw_lumen_area(bw_lumen_area == j) = area;

                try

                    ep_lb(i) = bw_epith_filled_label(statsl(i).PixelList(2,2),statsl(i).PixelList(2,1));

                catch

                    disp('hey-you');

                end

            end

        end

        j = 0;

        inde = 0;

        for i = 1:length(Be) % Loop over the epithelium

            num_skipped = 0;

            if statse(i).Area<100000000 && statse(i).Area>0

                % plot the epithelium

                inde = inde+1;

                output_data.data(x).gland(inde).boundarye = Be{i};

                %             output_data.data(x).gland(inde).epith_size = statse(i).Area;

                ep_masked_cells_im = uint8(bw_epith_filled_label == i).* bw_cells;

                ep_masked_lumen_im = uint8(bw_epith_filled_label == i).* bw_lumen;

                output_data.data(x).gland(inde).epith_size = statse(i).Area - numel(find(ep_masked_lumen_im));

                output_data.data(x).gland(inde).cell_frac = sum(sum(ep_masked_cells_im)) / (statse(i).Area - sum(sum(ep_masked_lumen_im)));

                boundarye = Be{i};

                ep_com = [round(mean(boundarye(:,1))) round(mean(boundarye(:,2)))];

                %             current_core = tissue_mask_label(ep_com(1),ep_com(2));

                skip_pic = 0;

                try

                    glnd_im_out = im_resized(ep_com(1)-111:ep_com(1)+112,ep_com(2)-111:ep_com(2)+112,:);

                catch

                    skip_pic = 1;

                end

                cidx = mod(i,length(colors))+1;

                plot(boundarye(:,2), boundarye(:,1),...

                    colors(cidx),'LineWidth',2);

                areae = statse(i).Area;

                delta_sqe = diff(boundarye).^2;

                perimetere = sum(sqrt(sum(delta_sqe,2)));

                output_data.data(x).gland(inde).epith_tort = 4*pi*areae/perimetere^2;

                % For plots -----------------------------------------------

                epith_size(i) = statse(i).Area - numel(find(ep_masked_lumen_im));

                cell_frac(i) = sum(sum(ep_masked_cells_im)) / (statse(i).Area - sum(sum(ep_masked_lumen_im)));

                epith_tort(i) = 4*pi*areae/perimetere^2;

                %----------------------------------------------------------

                %randomize text position for better visibility

                rndRow = ceil(length(boundarye)/(mod(rand*i,7)+1));

                col = boundarye(rndRow,2); row = boundarye(rndRow,1);

                % Figure out the lumen within

                indl = 0;

                if numel(find(ep_lb==i))>0

                    % loop over the lumen within:

                    for j = find(ep_lb==i)

                        try

                            if statsl(j).Area>1

                                indl = indl + 1;

                                % plot the lumen boundaries

                                boundaryl = Bl{j};

                                cidx = mod(i,length(colors))+1;

                                plot(boundaryl(:,2), boundaryl(:,1),...

                                    colors(cidx),'LineWidth',2);

                                %randomize text position for better visibility

                                rndRow = ceil(length(boundaryl)/(mod(rand*i,7)+1));

                                col = boundaryl(rndRow,2); row = boundaryl(rndRow,1);

                                area = statsl(j).Area;

                                delta_sq = diff(boundaryl).^2;

                                perimeter = sum(sqrt(sum(delta_sq,2)));

                                output_data.data(x).gland(inde).lumen(indl).lumen_tort = 4*pi*area/perimeter^2;

                                output_data.data(x).gland(inde).lumen(indl).wall_thick = 0;

                                output_data.data(x).gland(inde).lumen(indl).area = area;

                                % For plots -------------------------------

                                lumen_tort(j) = 4*pi*area/perimeter^2;

                                bw_lumen_roundness(bw_lumen_roundness == j) = lumen_tort(j); % replaces labeled lumen with tortuosity value

                                bw_lumen_area(bw_lumen_area == j) = area;

                                % -----------------------------------------

                                wall_thickness = [];

                                min_lum_index_thick = [];

                                for jj = 1:size(Bl{j},1)

                                    ind = 1;

                                    for k = 1:5:size(Be{i},1)

                                        wall_thickness(ind) = sqrt((Bl{j}(jj,1)-Be{i}(k,1))^2 + (Bl{j}(jj,2)-Be{i}(k,2))^2);

                                        ind = ind+1;

                                    end

                                    min_lum_index_thick(jj) = min(wall_thickness);

                                    if show_results == 1

                                        if rem(jj,30)==0

                                            h = text(Bl{j}(jj,2), Bl{j}(jj,1), sprintf('%0.2f',min_lum_index_thick(jj)));

                                            set(h,'Color',colors(cidx),...

                                                'FontSize',8);

                                        end

                                    end

                                end

                                output_data.data(x).gland(inde).lumen(indl).wall_thick = mean(min_lum_index_thick);

                                wall_thick(i) = mean(min_lum_index_thick);

                                if show_results == 1

                                    h = text(Bl{j}(1,2)+20, Bl{j}(1,1)+10, sprintf('%i\n%i\n%0.2f\n%0.2f',statse(i).Area,statsl(j).Area,output_data.data(x).gland(inde).lumen(indl).lumen_tort,output_data.data(x).gland(inde).lumen(indl).wall_thick))  ;

                                    set(h,'Color',colors(cidx),...

                                        'FontSize',8);

                                end

                                % output_data.data(x).gland(inde)

                                % disp('results shown');

                                %

%                                 try

                                fprintf(fid,'%s,%s,%i,%s,%i,%i,%i,%f,%f,%f,%f,%f,%f,%f,%f\n',...

                                    scanner{s},input_images{x},x,grade{x},...

                                    grade_number(x),inde,indl,stroma_area,epith_area,...

                                    output_data.data(x).gland(inde).epith_size,...

                                    output_data.data(x).gland(inde).epith_tort,...

                                    output_data.data(x).gland(inde).cell_frac,...

                                    output_data.data(x).gland(inde).lumen(indl).area,...

                                    output_data.data(x).gland(inde).lumen(indl).wall_thick,...

                                    output_data.data(x).gland(inde).lumen(indl).lumen_tort);

                            end

                        catch

                            disp('Debug point')

                        end

                    end

                    try

                        tot_lum_area = [];

                        tot_wall_thick = [];

                        tot_lum_tort = [];

                        for mmm = 1:numel(output_data.data(x).gland(inde).lumen(:))

                            tot_lum_area = [tot_lum_area output_data.data(x).gland(inde).lumen(mmm).area];

                            tot_wall_thick = [tot_wall_thick output_data.data(x).gland(inde).lumen(mmm).wall_thick];

                            tot_lum_tort = [tot_lum_tort output_data.data(x).gland(inde).lumen(mmm).lumen_tort];

                        end

                        fprintf(fid2,'%s,%s,%i,%s,%i,%i,%i,%f,%f,%f,%f,%f,%f,%f,%f\n',...

                            scanner{s},input_images{x},x,grade{x},grade_number(x),inde,...

                            output_data.data(x).gland(inde).epith_size,...

                            output_data.data(x).gland(inde).epith_tort,...

                            output_data.data(x).gland(inde).cell_frac,...

                            mean(tot_lum_area),...

                            mean(tot_wall_thick),...

                            mean(tot_lum_tort),...

                            sum(sum(bw_lumen))/sum(sum(core_mask)),...

                            sum(sum(bw_epith))/sum(sum(core_mask)),...

                            sum(sum(bw_stroma))/sum(sum(core_mask)));

                    catch

                        disp('debug point');

                    end

                else

                    disp('skipping epith, no lumen');

                    ep_masked_cells_im = uint8(bw_epith_filled_label == i).* bw_cells;

                    ep_masked_lumen_im = uint8(bw_epith_filled_label == i).* bw_lumen;

                    epith_size(i) = statse(i).Area - numel(find(ep_masked_lumen_im));

                    cell_frac(i) = sum(sum(ep_masked_cells_im)) / (statse(i).Area - sum(sum(ep_masked_lumen_im)));

                    boundarye = Be{i};

                    areae = statse(i).Area;

                    delta_sqe = diff(boundarye).^2;

                    perimetere = sum(sqrt(sum(delta_sqe,2)));

                    epith_tort(i) = 4*pi*areae/perimetere^2;

                end

            else

                disp('skipping epith, too big');

                boundaryl = Bl{j};

                area = statsl(j).Area;

                delta_sq = diff(boundaryl).^2;

                perimeter = sum(sqrt(sum(delta_sq,2)));

                lumen_tort(j) = 4*pi*area/(perimeter^2);

                bw_lumen_roundness(bw_lumen_roundness == j) = lumen_tort(j); % replaces labeled lumen with tortuosity value

                bw_lumen_area(bw_lumen_area == j) = area;

                ep_masked_cells_im = uint8(bw_epith_filled_label == i).* bw_cells;

                ep_masked_lumen_im = uint8(bw_epith_filled_label == i).* bw_lumen;

                epith_size(i) = statse(i).Area - numel(find(ep_masked_lumen_im));

                cell_frac(i) = sum(sum(ep_masked_cells_im)) / (statse(i).Area - sum(sum(ep_masked_lumen_im)));

                boundarye = Be{i};

                areae = statse(i).Area;

                delta_sqe = diff(boundarye).^2;

                perimetere = sum(sqrt(sum(delta_sqe,2)));

                epith_tort(i) = 4*pi*areae/perimetere^2;

            end

%% -------------------------------------------------------------

            % display epithelium tort, size, and cell frac

            bw_epith_roundness(bw_epith_roundness == i) = epith_tort(i);

            bw_epith_size(bw_epith_size == i) = epith_size(i);

            bw_cell_frac(bw_cell_frac == i) = cell_frac(i);

            try

                bw_epith_thickness(bw_epith_thickness == i) = wall_thick(i);

            catch

                bw_epith_thickness(bw_epith_thickness == i) = 0;

            end

%% -------------------------------------------------------------

        end

%% ------------------------------------------------------------------------

%         % display features

%         figure('Position',[100 100 1500 1500])

%         imshow(im_resized); axis image; hold on;

%         h = imagesc(bw_lumen_roundness); caxis([0 0.7]); % change caxis to see colormap better

%         hold off

%         

%         [M,N] = size(bw_lumen_roundness);

%         alpha_data = bw_lumen_roundness > 0;

%         alpha_data = alpha_data(1:M, 1:N);

%         set(h, 'AlphaData', alpha_data);

%         title('lumen tort 2');

%         

%         figure('Position',[100 100 1500 1500])

%         imshow(im_resized); axis image; hold on;

%         h = imagesc(bw_lumen_area); caxis([0 3*10^4]); % change caxis to see colormap better

%         hold off

%         

%         [M,N] = size(bw_lumen_area);

%         alpha_data = bw_lumen_area > 0;

%         alpha_data = alpha_data(1:M, 1:N);

%         set(h, 'AlphaData', alpha_data);

%         title('lumen area 2');

%         

%         figure('Position',[100 100 1500 1500])

%         imshow(im_resized); axis image; hold on;

%         h = imagesc(bw_epith_thickness);

%         hold off

%         

%         [M,N] = size(bw_epith_thickness);

%         alpha_data = bw_epith_thickness > 0;

%         alpha_data = alpha_data(1:M, 1:N);

%         set(h, 'AlphaData', alpha_data);

%         title('epith thick 2')

%         

%         figure('Position',[100 100 1500 1500])

%         imshow(im_resized); axis image; hold on;

%         h = imagesc(bw_epith_roundness); caxis([0 1]); % change caxis to see colormap better

%         hold off

%         

%         [M,N] = size(bw_epith_roundness);

%         alpha_data = bw_epith_roundness > 0;

%         alpha_data = alpha_data(1:M, 1:N);

%         set(h, 'AlphaData', alpha_data);

%         title('epith tort 2')

%         

%         figure('Position',[100 100 1500 1500])

%         imshow(im_resized); axis image; hold on;

%         h = imagesc(bw_epith_size); %caxis([0 250]);

%         hold off

%         

%         [M,N] = size(bw_epith_size);

%         alpha_data = bw_epith_size > 0;

%         alpha_data = alpha_data(1:M, 1:N);

%         set(h, 'AlphaData', alpha_data);

%         title('epith_size 2');

%         

%         figure('Position',[100 100 1500 1500])

%         imshow(im_resized); axis image; hold on;

%         h = imagesc(bw_cell_frac); caxis([0 1]);

%         hold off

%         

%         [M,N] = size(bw_cell_frac);

%         alpha_data = bw_cell_frac > 0;

%         alpha_data = alpha_data(1:M, 1:N);

%         set(h, 'AlphaData', alpha_data);

%         title('cell frac 2')

%% ------------------------------------------------------------------------

        fprintf('%s done\n',input_images{x});

        % saveas(gcf,sprintf('outline_seg_%s',input_images{x}),'tif');

        close all

    end

end

fclose(fid);

status = 'done';

end