%ts_histo_tibia_PONE_SAMprocessing.m % histomorphometric analysis from SAM images % % This script contains the SAM image processing described in the paper: % Iori G, Schneider J, Reisinger A, Heyer F, Peralta L, Wyers C, et al. Large cortical bone pores in the tibia are associated with proximal femur strength. PLOS ONE. doi:10.1371/journal.pone.0215405 % % preliminary operations: % - SAM RF signals are already converted to maps of the surface's acoustic impedance % - All images are rotated in a common APML ref. system: ANT looks towards the bottom of the image % - Endosteum masks already drawn % % SAM images, masks and master table for the analysis can be downloaded from: % DOI % % dependencies: % - MATLAB R2018a or later % - simpleCT https://github.com/gianthk/simpleCT % DOI: 10.5281/zenodo.2628321 % - MATLAB and Octave Functions for Computer Vision and Image Processing, Peter Kovesi % https://www.peterkovesi.com/matlabfns/ % - export_fig (if graphics TRUE) https://de.mathworks.com/matlabcentral/fileexchange/23629-export_fig % other dependencies that are included in the simpleCT toolbox % - dir2 https://de.mathworks.com/matlabcentral/fileexchange/21791-search-files-recursively-dir2 % - fwhm https://de.mathworks.com/matlabcentral/fileexchange/10590-fwhm % - max_inscribed_circle https://de.mathworks.com/matlabcentral/fileexchange/30805-maximum-inscribed-circle-using-distance-transform % - inpoly https://github.com/dengwirda/inpoly % % This script does: % 1. load Scanning Acoustic Microscopy (SAM) image of cortical bone % 2. Convert to acoustic impedance % 3. Wellner's adaptive threshold % 4. clean the binary mask % 5. correct the periosteum (PERI), endosteum (ENDO) and cortical bone (CORT) masks % 6. load evaluation (ROI) masks for specific compartments % 7. measure macroscopic morph (shaft cross-section geometry) % 8. calculate Ct.Th distribution9. calculate cortical thickness (Ct.Th) distribution % 10.1 calculate pore properties from binary image of the pores (Ct.Po, Po.D, Po.Dm) % 10.2 calculate pore properties for each eval ROI % 11. store single sample results % 12-13.save results % ______________________________________________________ % % Author: Gianluca Iori % BSRT - Charite Berlin % Created on: 25/04/2018 % Last update: 25/03/2019 % ______________________________________________________ %% init clear all close all clc graphics = true; currentdir = pwd; % master table filename_master = 'tacosound_tibia_left_master.csv'; % the master file is contained in the archive that you downloaded % path definitions homedir = currentdir; resultsdir = [homedir filesep 'results' filesep]; evalmasksdir = [homedir filesep 'masks' filesep '00_eval']; ROIUSmasksdir = [homedir filesep 'masks' filesep '00_eval' filesep 'ROIUS']; % create dir struct for analysis mkdir(resultsdir); mkdir([homedir filesep 'masks' filesep]); mkdir([homedir filesep 'masks' filesep 'tissue' filesep]); mkdir([homedir filesep 'masks' filesep 'endo' filesep]); mkdir([homedir filesep 'masks' filesep 'peri' filesep]); mkdir([homedir filesep 'masks' filesep 'cort' filesep]); mkdir([homedir filesep 'masks' filesep 'pores' filesep]); mkdir(evalmasksdir); mkdir(ROIUSmasksdir); warning('copy the files in the appropriate folders before running the analysis'); resultsdir = [resultsdir filesep datestr(now, 'yymmdd')]; resultsingledir = [resultsdir filesep 'single']; if exist(resultsdir)==0, mkdir(resultsdir); end if exist(resultsingledir)==0, mkdir(resultsingledir); end resultsfile = [resultsdir filesep 'histo.mat']; % SAM images for the analysis can be downloaded here: DOI % you should copy the files in the appropriate folders before running the analysis: % - SAM images -> homedir % - endosteum masks -> homedir/masks/endo % - ROI masks -> homedir/masks/00_eval Each folder contains the specific ROIs for the analysis. % ROIUS (anteromedial tibia) are provided. % You can draw more ROI masks and put them here.. they will be processed) % parameters pixelsize = 0.012; % [mm] 100MHz SAM ctth_BinEdges = 0.00:0.05:10.00; % [mm] ctth_bins = 0.025:0.05:9.975; % [mm] PoDm_BinEdges = 0.000:0.005:0.600; % [mm] PoDm_bins = 0.0025:0.005:0.5975; % [mm] PoDmprops = {'mean' 'median' 'std' 'min' 'max' 'Width' 'Mean' 'Peak' 'Variance' 'Skewness' 'Peakedness' '0.05 percentile' '0.95 percentile'}; PoDm_thresh = [0, 0.060 0.090 0.120]; % [mm] allporesDm = []; allporesDmUS = []; PoAprops = {'mean' 'median' 'std' 'min' 'max' 'Width' 'Peak' 'Variance' 'Skewness' 'Kurtosis' '0.10 quantile' '0.90 quantile'}; PoA_thresh = [0 0.005 0.008 0.012 0.016 0.032 0.05]; % [mm^2] PoDegCirprops = {'mean' 'median' 'std' 'min' 'max'}; sample = 1; %% load master if isempty(filename_master) || ~exist(filename_master) [FILENAME, PATHNAME] = uigetfile('*.csv;*.CSV', 'Select master table'); cd(currentdir); filename_master = [PATHNAME FILENAME]; end master = medtooltable.readtable(filename_master); cd(currentdir); %% eval all samples for sample = 1:length(master.specimen) close all %% checks samplename = master.specimen{sample}; samplename = regexprep(samplename, '_', ''); FILENAME = [homedir filesep samplename '.tif']; TISSUEfilename = [homedir filesep 'masks' filesep 'tissue' filesep samplename '.tif']; ENDOfilename = [homedir filesep 'masks' filesep 'endo' filesep samplename '.tif']; PERIfilename = [homedir filesep 'masks' filesep 'peri' filesep samplename '.tif']; CORTfilename = [homedir filesep 'masks' filesep 'cort' filesep samplename '.tif']; isright = ~isempty(regexp(samplename, 'R')); % LEFT (L) or RIGHT (R) flag % run flag in mastertable if ~master.run(sample), continue; end fprintf('evaluate specimen: \t%s\n', char(master.specimen(sample))); %% 1. load Scanning Acoustic Microscopy (SAM) image of cortical bone % already rotated image in the APSI (Anterior-Posterior-Superior-Inferior) system SAM = imread(FILENAME); %% 2. convert to acoustic impedance IMP = double(SAM)/4369 ; %% 3. Wellner's adaptive threshold [BW] = threshbone(IMP, 'adaptive', 3, 22, 'gaussian', 'relative'); % adaptive threshold %% 4. clean the binary mask % remove non connected particles BW2 = removeparticles(BW, 'connectivity'); % remove single pixel pores BW2 = bwmorph(BW2, 'fill'); % ..show before and after if graphics figure; imshowpair(BW2, BW); axis image; title('tissue mask cleaned'); end %% 5. correct the masks PERI, ENDO and CORT % Read the endosteum mask which was drawn manually using the set of rules proposed in (Malo et. al. 2013) ENDOfilename_in = [homedir filesep 'masks' filesep 'endo' filesep samplename '.tif']; endo_BW = imread(ENDOfilename_in); endo_BW0 = endo_BW; % get boundary of the endosteum mask endo_bound_BW = bwmorph(~endo_BW, 'remove'); % calculate the periosteum mask peri_BW = imfill(BW2 | endo_bound_BW, 'holes'); % get boundary of the periosteum mask peri_bound_BW = bwmorph(peri_BW, 'remove'); % final version of the endosteum mask from a combination with the tissue mask endo_BW = endosteummask(endo_BW0 & BW2 | peri_bound_BW, 'regionprops'); % cortical bone mask as the combination of the two cort_BW = endo_BW & peri_BW; % remove non connected objects with an area < 300 pixels from the binary image cort_BW = bwareaopen(cort_BW, 300); % save the masks imwrite(endo_BW , ENDOfilename); imwrite(peri_BW , PERIfilename); imwrite(cort_BW , CORTfilename); if graphics figure; endo = gcf; imshowpair(IMP, endo_BW); axis image; title('endo'); figure; peri = gcf; imshowpair(IMP, peri_BW); axis image; title('peri'); figure; cort = gcf; imshowpair(IMP, cort_BW); axis image; title('cort'); end %% 6. load evaluation (ROI) masks fprintf('Loading ROI masks...\n'); evalmasks_content = dir2(evalmasksdir); for j=1:length(evalmasks_content) evalmasks(j).name = evalmasks_content(j).name; evalmasks(j).filename = [evalmasksdir filesep evalmasks_content(j).name filesep samplename '.tif']; evalmasks(j).BW = logical(imread(evalmasks(j).filename)); evalmasks(j).cort_BW = evalmasks(j).BW & cort_BW; evalmasks(j).polygon = periosteumcontour(flipud(evalmasks(j).BW)); evalmasks(j).polygon = periosteumcontour(evalmasks(j).BW); end %% 7. macroscopic morph (shaft cross-section geometry): (BA, TA, CA, BV/TV, Ixx, Iyy) tic % 6.2 Tissue Area (only the bone tissue) (T.Ar) TAr = sum(BW2(:)); % [pixels] % 6.3 Cortical Area (area of the cortical compartment) (Ct.Ar) CtAr = sum(cort_BW(:)); % [pixels] % 6.4 Core Area (the cross-sectional area occupied by the bone) (C.Ar) CAr = sum(peri_BW(:)); % [pixels] % 6.5 areal portion of cortex (Ct.Wba) % Ratio between the cortical tissue area and the whole bone area (containing both the trabecular and cortical tissue). [Malo et. al. 2013] CtWba = sum(cort_BW(:) & BW2(:)) / CAr; % [%] % 6.6 moment of inertia MoI = momentofinertia(BW2); % [pixels^2] % 6.7 second moments of area [Ixx, Iyy] = momentofinertia_area(BW2); % [pixels^4] % convert to mm TAr = TAr * pixelsize^2; % [mm^2] CtAr = CtAr * pixelsize^2; % [mm^2] CAr = CAr * pixelsize^2; % [mm^2] MoI = MoI * pixelsize^2; % [mm^2] Ixx = Ixx * pixelsize^4; % [mm^4] Iyy = Iyy * pixelsize^4; % [mm^4] fprintf('calculation of macroscopic morphological properties done!\t'); toc %% 8. Ct.Th distribution tic % get point clouds on the endosteum and periosteum contours peri_points = periosteumcontour(imclose(peri_BW, strel('disk', 30))); endo_points = endosteumcontour(endo_BW); % get only those belonging to a specific ROI for j=1:length(evalmasks_content) evalmasks(j).endo_points = endo_points(inpolygon(endo_points(:,1), endo_points(:,2), evalmasks(j).polygon(:,1), evalmasks(j).polygon(:,2)),:); evalmasks(j).peri_points = peri_points(inpolygon(peri_points(:,1), peri_points(:,2), evalmasks(j).polygon(:,1), evalmasks(j).polygon(:,2)),:); end ctth_filename = [resultsdir filesep samplename '_CtTh']; % output figure if graphics [ctth, ctth_array] = CtTh(peri_points, endo_points, 'pdist2', 1, 'peak'); CtTh_fig = gcf; savefig(ctth_filename); export_fig(ctth_filename, '-dpng'); % you will need export_fig for this.. you can modify this to 'print' if you don't have export_fig close(CtTh_fig); else [ctth, ctth_array] = CtTh(peri_points, endo_points, 'pdist2', 0, 'peak'); end % convert Ct.Th to mm ctth = ctth * pixelsize; ctth_array = ctth_array * pixelsize; % histogram of the minimum Ct.Th ctth_hist = histcounts(ctth_array, ctth_BinEdges, 'Normalization', 'pdf'); % repeat for each eval ROI for j=1:length(evalmasks_content) if graphics [evalmasks(j).ctth, evalmasks(j).ctth_array] = CtTh(evalmasks(j).peri_points, evalmasks(j).endo_points, 'pdist2', 1, 'peak'); CtTh_fig = gcf; ctth_filename = [evalmasks(j).filename(1:end-4) '_CtTh']; savefig(ctth_filename); export_fig(ctth_filename, '-dpng'); close(CtTh_fig); else [evalmasks(j).ctth, evalmasks(j).ctth_array] = CtTh(evalmasks(j).peri_points, evalmasks(j).endo_points, 'pdist2', 0, 'peak'); end evalmasks(j).ctth = evalmasks(j).ctth * pixelsize; evalmasks(j).ctth_array = evalmasks(j).ctth_array * pixelsize; evalmasks(j).ctth_hist = histcounts(evalmasks(j).ctth_array, ctth_BinEdges, 'Normalization', 'pdf'); end fprintf('calculation of cortical thickness done!\t'); toc %% 10.1 pore props from BW of the pores (Ct.Po, Po.D, Po.Dm) tic % pores mask pores_BW = ~BW2 & cort_BW; % remove pores connected to the periosteum pores_BW = pores_BW & ~bwareaopen(pores_BW | ~peri_BW, 5000); % remove pores connected to the endosteum pores_BW = pores_BW & ~bwareaopen(pores_BW | ~endo_BW, 5000); % save it imwrite(pores_BW , [homedir filesep 'masks' filesep 'pores' filesep samplename '.tif']); % Ct.Po CtPo = 100*(sum(pores_BW(:))/sum(cort_BW(:))); % [%] clear PoDmres PoDm_L PoDm_M PoDm_A PoDm_P PoDm_array PoA_array % Calculate pore properties PoSTATS = regionprops(pores_BW, 'Area'); % Calculate Po.Dm distribution for k=1:length(PoSTATS) PoA_array(k) = PoSTATS(k).Area; end PoDm_array = PoDm(pores_BW); PoDm_array = PoDm_array * pixelsize; % [mm] allporesDm = [allporesDm; PoDm_array]; PoA_array = PoA_array * pixelsize^2; % [mm^2] % Degree of circularity PoDegCir_array = (pi*(PoDm_array.^2)/4)./PoA_array'; % Po.Dm and Po.A histogram [PoDm_hist, PoDm_BinEdges] = histcounts(PoDm_array, 'Normalization', 'pdf'); PoDm_bins = PoDm_BinEdges(1:end-1) + (PoDm_BinEdges(2)-PoDm_BinEdges(1))/2; PoDm_binwidth = PoDm_BinEdges(2:end) - PoDm_BinEdges(1:end-1); [PoA_hist, PoA_BinEdges] = histcounts(PoA_array, 'Normalization', 'pdf'); PoA_bins = PoA_BinEdges(1:end-1) + (PoA_BinEdges(2)-PoA_BinEdges(1))/2; PoA_binwidth = PoA_BinEdges(2:end) - PoA_BinEdges(1:end-1); % relative proportion of pores with DIAMETER over a certain threshold and ammount of surface occupied per square mm DmAtot = PoDm_binwidth*PoDm_hist'; % total area under the curve for t = 1:length(PoDm_thresh) idx = find(PoDm_array>PoDm_thresh(t)); idx_hist = find(PoDm_bins>PoDm_thresh(t)); DmAt = PoDm_binwidth(idx_hist)*PoDm_hist(idx_hist)'; % area under the curve after threshold % relative proportion of pores with diameter over a certain threshold PoDm_rel(t) = DmAt / DmAtot; % ammount of surface occupied PoDmA_rel(t) = sum(PoA_array(idx)); end PoDmA_rel = PoDmA_rel / CtAr; % 1/BVTV [0-1] PoDmCtPo_rel = PoDmA_rel*100; % [%] PoDmA_rel = (PoDmA_rel/PoDmA_rel(1))*100; % [%] % relative proportion of pores with AREA over a certain threshold and ammount of surface occupied AAtot = PoA_binwidth*PoA_hist'; % total area under the curve for t = 1:length(PoA_thresh) idx = find(PoA_array>PoA_thresh(t)); idx_hist = find(PoA_bins>PoA_thresh(t)); AAt = PoA_binwidth(idx_hist)*PoA_hist(idx_hist)'; % area under the curve after threshold % relative proportion of pores with diameter over a certain threshold PoA_rel(t) = AAt / AAtot; % ammount of surface occupied PoAA_rel(t) = sum(PoA_array(idx)); end PoAA_rel = PoAA_rel / CtAr; % 1/BVTV [0-1] PoACtPo_rel = PoAA_rel*100; % [%] PoAA_rel = (PoAA_rel/PoAA_rel(1))*100; % [%] % ksdensity of the Po.Dm distribution [prob, idx] = ksdensity(PoDm_array); % Po.Dm distr properties PoDmres = [mean(PoDm_array) median(PoDm_array) std(PoDm_array) min(PoDm_array) max(PoDm_array) fwhm(PoDm_bins, PoDm_hist) idx(find(prob == max(prob))) var(PoDm_array) skewness(PoDm_array) kurtosis(PoDm_array) quantile(PoDm_array, 0.1) quantile(PoDm_array, 0.9)]; % ksdensity of the Po.A distribution [prob, idx] = ksdensity(PoA_array); % Po.A distr properties PoA = [mean(PoA_array) median(PoA_array) std(PoA_array) min(PoA_array) max(PoA_array) fwhm(PoA_bins, PoA_hist) idx(find(prob == max(prob))) var(PoA_array) skewness(PoA_array) kurtosis(PoA_array) quantile(PoA_array, 0.1) quantile(PoA_array, 0.9)]; % Pore degree of circularity results PoDegCir = [mean(PoDegCir_array) median(PoDegCir_array) std(PoDegCir_array) min(PoDegCir_array) max(PoDegCir_array)]; if graphics % Save figure to results PoDmhist_filename = [resultsdir filesep samplename '_PoDmHist']; figure; histogram(PoDm_array*1000, PoDm_BinEdges*1000, 'Normalization', 'pdf'); hold on; % ksdensity plot(PoDm_bins*1000, ksdensity(1000*PoDm_array(:), 1000*PoDm_bins), '-.r','LineWidth',1.5) xlabel('Po.Dm [microns]') ylabel('Probability'); title('TOT'); % axis([0 300 0 1]) axis 'auto y' grid MINOR legend('hist', 'ksdensity'); PoDmhist_fig = gcf; savefig(PoDmhist_filename); export_fig(PoDmhist_filename, '-dpng'); close(PoDmhist_fig); PoAhist_filename = [resultsdir filesep samplename '_PoAHist']; figure; histogram(PoA_array, PoA_BinEdges, 'Normalization', 'pdf'); hold on; % ksdensity plot(PoA_bins, ksdensity(PoA_array(:), PoA_bins), '-.r','LineWidth',1.5) xlabel('Po.A [mm^2]') ylabel('Probability'); title('TOT'); % axis([0 0.04 0 1]) axis 'auto y' grid MINOR legend('hist', 'ksdensity'); PoAhist_fig = gcf; savefig(PoAhist_filename); export_fig(PoAhist_filename, '-dpng'); close(PoAhist_fig); end % Po.D of pores with diameter > threshold for t = 1:length(PoDm_thresh) PoD(t) = sum(PoDm_array>PoDm_thresh(t)) / CtAr; % [1/mm^2] Pon(t) = sum(PoDm_array>PoDm_thresh(t)); end % Po.D of pores with diameter > threshold for t2 = 1:length(PoA_thresh) PoD(t+t2) = sum(PoA_array>PoA_thresh(t2)) / CtAr; % [1/mm^2] Pon(t+t2) = sum(PoA_array>PoA_thresh(t2)); end %% 10.2 pore props for each evalmask for j=1:length(evalmasks_content) CtAr_evalmask = sum(evalmasks(j).cort_BW(:)) * pixelsize^2; % [mm^2] % pores mask % evalmasks(j).pores_BW = ~BW2 & imerode(evalmasks(j).cort_BW, strel('disk', 4)); evalmasks(j).pores_BW = pores_BW & evalmasks(j).cort_BW; imwrite(evalmasks(j).pores_BW , [evalmasks(j).filename(1:end-4) '_pores.tif']); % Ct.Po evalmasks(j).CtPo = 100*(sum(evalmasks(j).pores_BW(:))/sum(evalmasks(j).cort_BW(:))); % Calculate Po.Dm distribution PoSTATS = regionprops(evalmasks(j).pores_BW, 'Area'); evalmasks(j).PoA_array = []; for k=1:length(PoSTATS) evalmasks(j).PoA_array(k) = PoSTATS(k).Area; end evalmasks(j).PoDm_array = PoDm(evalmasks(j).pores_BW); evalmasks(j).PoDm_array = evalmasks(j).PoDm_array * pixelsize; allporesDmUS = [allporesDmUS; evalmasks(j).PoDm_array]; evalmasks(j).PoA_array = evalmasks(j).PoA_array * pixelsize^2; % Degree of circularity evalmasks(j).PoDegCir_array = (pi*(evalmasks(j).PoDm_array.^2)/4)./(evalmasks(j).PoA_array'); % Calculate Po.Dm histogram evalmasks(j).PoDm_hist = histcounts(evalmasks(j).PoDm_array, PoDm_BinEdges, 'Normalization', 'pdf'); evalmasks(j).PoA_hist = histcounts(evalmasks(j).PoA_array, PoA_BinEdges, 'Normalization', 'pdf'); % relative proportion of pores with diameter over a certain threshold % and ammount of surface occupied per square mm DmAtot = PoDm_binwidth*evalmasks(j).PoDm_hist'; % total area under the curve for t = 1:length(PoDm_thresh) idx = find(evalmasks(j).PoDm_array>PoDm_thresh(t)); idx_hist = find(PoDm_bins>PoDm_thresh(t)); DmAt = PoDm_binwidth(idx_hist)*evalmasks(j).PoDm_hist(idx_hist)'; % area under the curve after threshold % relative proportion of pores with diameter over a certain threshold evalmasks(j).PoDm_rel(t) = DmAt / DmAtot; % ammount of surface occupied per square mm evalmasks(j).PoDmA_rel(t) = sum(evalmasks(j).PoA_array(idx)); end evalmasks(j).PoDmA_rel = evalmasks(j).PoDmA_rel / CtAr_evalmask; % [0-1] evalmasks(j).PoDmCtPo_rel = evalmasks(j).PoDmA_rel*100; % [%] evalmasks(j).PoDmA_rel = (evalmasks(j).PoDmA_rel/evalmasks(j).PoDmA_rel(1))*100; % [%] % relative proportion of pores with Area over a certain threshold % and ammount of surface occupied per square mm AAtot = PoA_binwidth*evalmasks(j).PoA_hist'; % total area under the curve for t = 1:length(PoA_thresh) idx = find(evalmasks(j).PoA_array>PoA_thresh(t)); idx_hist = find(PoA_bins>PoA_thresh(t)); AAt = PoA_binwidth(idx_hist)*evalmasks(j).PoA_hist(idx_hist)'; % area under the curve after threshold % relative proportion of pores with diameter over a certain threshold evalmasks(j).PoA_rel(t) = AAt / AAtot; % ammount of surface occupied per square mm evalmasks(j).PoAA_rel(t) = sum(evalmasks(j).PoA_array(idx)); end evalmasks(j).PoAA_rel = evalmasks(j).PoAA_rel / CtAr_evalmask; % [0-1] evalmasks(j).PoACtPo_rel = evalmasks(j).PoAA_rel*100; % [%] evalmasks(j).PoAA_rel = (evalmasks(j).PoAA_rel/evalmasks(j).PoAA_rel(1))*100; % [%] % ksdensity of the Po.Dm distribution [prob, idx] = ksdensity(evalmasks(j).PoDm_array); % Po.Dm distr properties evalmasks(j).PoDm = [mean(evalmasks(j).PoDm_array) median(evalmasks(j).PoDm_array) std(evalmasks(j).PoDm_array) min(evalmasks(j).PoDm_array) max(evalmasks(j).PoDm_array) fwhm(PoDm_bins, evalmasks(j).PoDm_hist) idx(find(prob == max(prob))) var(evalmasks(j).PoDm_array) skewness(evalmasks(j).PoDm_array) kurtosis(evalmasks(j).PoDm_array) quantile(evalmasks(j).PoDm_array, 0.1) quantile(evalmasks(j).PoDm_array, 0.9)]; % ksdensity of the Po.Dm distribution [prob, idx] = ksdensity(evalmasks(j).PoA_array); % Po.A distr properties evalmasks(j).PoA = [mean(evalmasks(j).PoA_array) median(evalmasks(j).PoA_array) std(evalmasks(j).PoA_array) min(evalmasks(j).PoA_array) max(evalmasks(j).PoA_array) fwhm(PoA_bins, evalmasks(j).PoA_hist) idx(find(prob == max(prob))) var(evalmasks(j).PoA_array) skewness(evalmasks(j).PoA_array) kurtosis(evalmasks(j).PoA_array) quantile(evalmasks(j).PoA_array, 0.1) quantile(evalmasks(j).PoA_array, 0.9)]; % Pore degree of circularity results evalmasks(j).PoDegCir = [mean(evalmasks(j).PoDegCir_array) median(evalmasks(j).PoDegCir_array) std(evalmasks(j).PoDegCir_array) min(evalmasks(j).PoDegCir_array) max(evalmasks(j).PoDegCir_array)]; if graphics % Save figure to results PoDmhist_filename = [resultsdir filesep samplename '_PoDmHist_' evalmasks(j).name]; figure; histogram(evalmasks(j).PoDm_array*1000, PoDm_BinEdges*1000, 'Normalization', 'pdf'); hold on; % ksdensity plot(PoDm_bins*1000, ksdensity(1000*evalmasks(j).PoDm_array(:), 1000*PoDm_bins), '-.r','LineWidth',1.5) xlabel('Po.Dm [microns]') ylabel('Probability'); title(evalmasks(j).name); % axis([0 300 0 1]) axis 'auto y' grid MINOR legend('hist', 'ksdensity'); PoDmhist_fig = gcf; savefig(PoDmhist_filename); export_fig(PoDmhist_filename, '-dpng'); close(PoDmhist_fig); PoAhist_filename = [resultsdir filesep samplename '_PoAHist_' evalmasks(j).name]; figure; histogram(evalmasks(j).PoA_array, PoA_BinEdges, 'Normalization', 'pdf'); hold on; % ksdensity plot(PoA_bins, ksdensity(evalmasks(j).PoA_array(:), PoA_bins), '-.r','LineWidth',1.5) xlabel('Po.A [mm^2]') ylabel('Probability'); title(evalmasks(j).name); % axis([0 0.04 0 1]) axis 'auto y' grid MINOR legend('hist', 'ksdensity'); PoAhist_fig = gcf; savefig(PoAhist_filename); export_fig(PoAhist_filename, '-dpng'); close(PoAhist_fig); end % Po.D of pores with diameter > threshold for t = 1:length(PoDm_thresh) evalmasks(j).PoD(t) = sum(evalmasks(j).PoDm_array>PoDm_thresh(t)) / CtAr_evalmask; % [1/mm^2] evalmasks(j).Pon(t) = sum(evalmasks(j).PoDm_array>PoDm_thresh(t)); end % Po.D of pores with Area > threshold for t2 = 1:length(PoA_thresh) evalmasks(j).PoD(t+t2) = sum(evalmasks(j).PoA_array>PoA_thresh(t2)) / CtAr_evalmask; % [1/mm^2] evalmasks(j).Pon(t+t2) = sum(evalmasks(j).PoA_array>PoA_thresh(t2)); end end %% 11. save single sample results resfile = [resultsingledir filesep samplename '.mat']; res.script = mfilename; res.sample = samplename; res.evalmasks = evalmasks; res.pixelsize = pixelsize; res.CtTh_BinEdges = ctth_BinEdges; res.CtTh_bins = ctth_bins; res.PoDm_array = PoDm_array; res.PoDm_BinEdges = PoDm_BinEdges; res.PoDm_bins = PoDm_bins; res.PoDmprops = PoDmprops; res.PoDm_thresh = PoDm_thresh; res.PoDm = PoDmres; res.PoDm_hist = PoDm_hist; res.PoDm_rel = PoDm_rel; res.PoDmCtPo_rel = PoDmCtPo_rel; res.PoDmA_rel = PoDmA_rel; res.PoA_array = PoA_array; res.PoA_BinEdges = PoA_BinEdges; res.PoA_bins = PoA_bins; res.PoAprops = PoAprops; res.PoA_thresh = PoA_thresh; res.PoA = PoA; res.PoA_hist = PoA_hist; res.PoA_rel = PoA_rel; res.PoACtPo_rel = PoACtPo_rel; res.PoAA_rel = PoAA_rel; res.TAr = TAr; res.CtAr = CtAr; res.CAr = CAr; res.CtWba = CtWba; res.MoI = MoI; res.Ixx = Ixx; res.Iyy = Iyy; res.peri_points = peri_points; res.endo_points = endo_points; res.CtTh = ctth; res.CtTh_array = ctth_array; res.CtPo = CtPo; res.PoDegCirprops = PoDegCirprops; res.PoDegCir = PoDegCir; res.PoD = PoD; res.Pon = Pon; save(resfile, 'res', '-v7.3'); %% 12. store results for all samples histo.script = mfilename; histo.samplename{sample} = samplename; histo.pixelsize = pixelsize; histo.PoDmprops = PoDmprops; histo.PoDm_thresh = PoDm_thresh; histo.PoDm_CORT(sample,:) = PoDmres; histo.PoDm_rel_CORT(sample,:) = PoDm_rel; histo.allporesDm_CORT = allporesDm; histo.PoDmCtPo_rel_CORT(sample,:) = PoDmCtPo_rel; histo.PoDmA_rel_CORT(sample,:) = PoDmA_rel; histo.PoDmCtPo_rel_CORT(sample,:) = PoDmCtPo_rel; histo.PoDmA_rel_CORT(sample,:) = PoDmA_rel; histo.PoAprops = PoAprops; histo.PoA_thresh = PoA_thresh; histo.PoA_rel_CORT(sample,:) = PoA_rel; histo.PoACtPo_rel_CORT(sample,:) = PoACtPo_rel; histo.PoAA_rel_CORT(sample,:) = PoAA_rel; histo.PoA_CORT(sample,:) = PoA; histo.TAr(sample) = TAr; histo.CtAr(sample) = CtAr; histo.CAr(sample) = CAr; histo.CtWba(sample) = CtWba; histo.MoI(sample) = MoI; histo.Ixx(sample) = Ixx; histo.Iyy(sample) = Iyy; histo.CtTh_BinEdges = ctth_BinEdges; histo.CtTh_bins = ctth_bins; histo.CtTh_CORT(sample) = ctth; histo.CtPo_CORT(sample) = CtPo; histo.PoD_CORT(sample,:) = PoD; histo.Pon_CORT(sample,:) = Pon; histo.PoDegCir_CORT(sample,:) = PoDegCir; histo.PoDegCirprops = PoDegCirprops; for j=1:length(evalmasks_content) histo.evalmasks{j} = evalmasks(j).name; histo.allporesDm{j} = allporesDmUS; histo.CtTh(sample,j) = evalmasks(j).ctth; histo.CtPo(sample,j) = evalmasks(j).CtPo; histo.PoDm(sample,:,j) = evalmasks(j).PoDm; histo.PoDm_rel(sample,:,j) = evalmasks(j).PoDm_rel; histo.PoA_rel(sample,:,j) = evalmasks(j).PoA_rel; histo.PoDmCtPo_rel(sample,:,j) = evalmasks(j).PoDmCtPo_rel; histo.PoDmA_rel(sample,:,j) = evalmasks(j).PoDmA_rel; histo.PoACtPo_rel(sample,:,j) = evalmasks(j).PoACtPo_rel; histo.PoAA_rel(sample,:,j) = evalmasks(j).PoAA_rel; histo.PoA(sample,:,j) = evalmasks(j).PoA; histo.PoA(sample,:,j) = evalmasks(j).PoA; histo.PoDegCir(sample,:,j) = evalmasks(j).PoDegCir; histo.PoD(sample,:,j) = evalmasks(j).PoD; histo.Pon(sample,:,j) = evalmasks(j).Pon; end end %% 13. save results save(resultsfile, 'histo', '-v7.3');