Datasets
Models
Downloads: 1
Diff of
/prostate_pathomics.m
[000000]
..
[be1edc]
Switch to side-by-side view
--- a +++ b/prostate_pathomics.m @@ -0,0 +1,479 @@ + 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 \ No newline at end of file