diff --git a/LICENSE.txt b/LICENSE.txt
new file mode 100644
index 0000000000000000000000000000000000000000..65a52f5ce43a68a6e2ddda5e3d9f138d6895ade5
--- /dev/null
+++ b/LICENSE.txt
@@ -0,0 +1,7 @@
+Copyright 2018 Vedrana Andersen Dahl, vand@dtu.dk
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
\ No newline at end of file
diff --git a/bonesub.mat b/bonesub.mat
new file mode 100644
index 0000000000000000000000000000000000000000..bd77965c669ceec5496b473d391650f178f99d05
Binary files /dev/null and b/bonesub.mat differ
diff --git a/functions/color_skeleton.m b/functions/color_skeleton.m
new file mode 100644
index 0000000000000000000000000000000000000000..715e594b1deb9db723ec2131270ebdee9c28fcef
--- /dev/null
+++ b/functions/color_skeleton.m
@@ -0,0 +1,37 @@
+function [U,cmap,cbar] = color_skeleton(B,S,D)
+%COLOR_SKELETON Colorize binarization, skeletonization and distances
+% [U,CMAP,CBAR] = COLOR_SKELETON(B,S,D) produces an indexed
+% image and a corresponding colormap for visualization of binarization,
+% skeletonization and distances.
+%
+% B - binary volume (zero outsice, one inside)
+% S - binary volume for skeleton of both inside and outside
+% D - signed distance field (negative for inside, positive for outside)
+% U - indexed uint8 image
+% CMAP - associated colormap
+% CBAR - struct containing inputs for associated colorbar
+%
+% Example:
+% [U,cmap,cbar] = color_skeleton(B,S,D);
+% imagesc(U)
+% colormap(cmap)
+% colorbar(cbar{:})
+%
+% Copyright 2018 Vedrana Andersen Dahl, vand@dtu.dk
+
+
+% Colormap values correspond to
+% 0 - outside (black)
+% 1 - inside (white)
+% 2--128 positive (outside) distance (127 colors from winter colormap)
+% 129--255 negative (inside) distance (127 colors from autumn colormap)
+
+U = uint8(B);
+ds = D(S);
+minmax = [min(ds),max(ds)]; % min should be negative
+ds(ds>0) = 2+126*(minmax(2)-ds(ds>0))/minmax(2); % flipping for symmetry
+ds(ds<0) = 129+126*ds(ds<0)/minmax(1);
+U(S) = ds;
+cmap = [0 0 0; 1 1 1; autumn(127);flip(winter(127))];
+cbtl = round(10*[minmax(2),0,-minmax(1)])/10; % colorbar tick labels
+cbar = {'Limits',[2,255],'Ticks',[2,128.5,255],'TickLabels',cbtl};
diff --git a/functions/compute_area.m b/functions/compute_area.m
new file mode 100644
index 0000000000000000000000000000000000000000..2d0a501b5213bd38c29cf80e4b52ad4de24abbb5
--- /dev/null
+++ b/functions/compute_area.m
@@ -0,0 +1,8 @@
+function area = compute_area(B)
+%COMPUTE_AREA Area of the binarized object
+% AREA = COMPUTE_AREA(B)
+
+Ax = B(2:end,:,:)~=B(1:end-1,:,:);
+Ay = B(:,2:end,:)~=B(:,1:end-1,:);
+Az = B(:,:,2:end)~=B(:,:,1:end-1);
+area = sum(Ax(:))+sum(Ay(:))+sum(Az(:));
\ No newline at end of file
diff --git a/functions/fit_distribution.m b/functions/fit_distribution.m
new file mode 100644
index 0000000000000000000000000000000000000000..de275e52051f102efdba7871542422bcf3d596fa
--- /dev/null
+++ b/functions/fit_distribution.m
@@ -0,0 +1,34 @@
+function measures = fit_distribution(r,show_histogram)
+% FIT_DISTRIBUTION Fits lognormal distribution to data
+% MEASURES = FIT_DISTRIBUTION(R,SHOW_HISTOGRAM)
+% R is a vector of data
+% SHOW_HISTOGRAM is a flag idicating whether plot should be produced
+% MEASURES are valuses extracted from fitted distribution: median
+% (geometric mean), geometric standard deviation, mean, variance and mode
+%
+% Copyright 2018 Vedrana Andersen Dahl, vand@dtu.dk
+
+measures = NaN(1,5);
+if numel(r)>1
+ pd = fitdist(r(:),'lognormal'); % pd.mu and pd.sigma are mean and standard deviation of the associated normal distribution
+ measures(1) = exp(pd.mu); % median of the fitted distribution (geometric mean)
+ measures(2) = exp(pd.sigma); % geometric standard deviation of the fitted distribution
+ measures(3) = exp(pd.mu+pd.sigma.^2/2); % mean of lognormal distribution
+ measures(4) = (exp(pd.sigma.^2)-1)*measures(3).^2; % variance of the lognormal distribution
+ measures(5) = exp(pd.mu-pd.sigma.^2); % mode of the lognormal distribution
+end
+
+if nargin>1 && show_histogram && numel(r)>1
+ h = histfit(r(:),ceil(max(r(:))),'lognormal'); hold on
+ x = h(2).XData;
+ y = h(2).YData;
+ title({['mean = ',num2str(measures(3)),', std = ',num2str(sqrt(measures(4)))],...
+ ['median = ',num2str(measures(1)),', gstd = ',num2str(measures(2))],...
+ ['mode = ',num2str(measures(5))]})
+ l(1) = plot(measures(5)*[1,1],[0,interp1(x,y,measures(5))],'c--');
+ l(2) = plot(measures(1)*[1,1],[0,interp1(x,y,measures(1))],'g--');
+ l(3) = plot(measures(3)*[1,1],[0,interp1(x,y,measures(3))],'m--');
+ xlim([0,max(x)])
+ xlabel('radius (pixels)'), ylabel('count')
+ legend(l,{'mode','median (geometric mean)','mean'},'Location','NE')
+end
diff --git a/functions/show_matvol.m b/functions/show_matvol.m
new file mode 100644
index 0000000000000000000000000000000000000000..4006df02889f9b13449c29d9863cb9ef582e4510
--- /dev/null
+++ b/functions/show_matvol.m
@@ -0,0 +1,47 @@
+function show_matvol(V,clim)
+% SHOW_MATVOL Shows volumetric data
+% show_matvol(V,clim)
+% Inputs: volume, grayscale color limits (optional)
+% Copyright 2018 Vedrana Andersen Dahl, vand@dtu.dk
+
+figure('Units','Normalized','Position',[0.1 0.3 0.5 0.6],...
+ 'KeyPressFcn',@key_press);
+dim = size(V);
+if nargin<2 || isempty(clim)
+ clim = [min(double(V(:))),max(double(V(:)))];
+end
+z = round(0.5*(dim(3)+1));
+%update_drawing
+i = imagesc(V(:,:,z),clim);
+title(['slice ',num2str(z),'/',num2str(dim(3))]), axis image ij
+colormap gray, drawnow
+
+%%%%%%%%%% CALLBACK FUNCTIONS %%%%%%%%%%
+ function key_press(~,object)
+ % keyboard commands
+ switch object.Key
+ case 'uparrow'
+ z = min(z+1,dim(3));
+ case 'downarrow'
+ z = max(z-1,1);
+ case 'rightarrow'
+ z = min(z+10,dim(3));
+ case 'leftarrow'
+ z = max(z-10,1);
+ case 'pagedown'
+ z = min(z+50,dim(3));
+ case 'pageup'
+ z = max(z-50,1);
+ end
+ update_drawing
+ end
+
+%%%%%%%%%% HELPING FUNCTIONS %%%%%%%%%%
+ function update_drawing
+ set(i,'CData',V(:,:,z))
+ title(['slice ',num2str(z),'/',num2str(dim(3))])
+ drawnow
+ end
+
+end
+
diff --git a/porosity_analysis_bone.m b/porosity_analysis_bone.m
new file mode 100644
index 0000000000000000000000000000000000000000..6c93a1dd99b3c65fabf518020995c26471842021
--- /dev/null
+++ b/porosity_analysis_bone.m
@@ -0,0 +1,155 @@
+% Copyright 2018 Vedrana Andersen Dahl, vand@dtu.dk
+clear
+close all
+addpath functions
+
+load bonesub
+gray_th = 3*10^4;
+div = [0.5 30.5,70.5,100.5,130.5,170.5 size(V,3)+0.5]; % layer divisions
+midlayer = floor(0.5*(div(1:end-1) + div(2:end))); % middle slice of a layer
+
+%% volume inspection -- use arrows
+show_matvol(permute(V,[3,2,1]))
+hold on, plot([1;size(V,2)]*ones(1,5),[1;1]*div(2:end-1),'r','LineWidth',2)
+text(size(V,2)/2*ones(1,6),midlayer,num2cell(1:6),'Color','r','FontSize',20)
+show_matvol(V)
+
+%% still figures
+x = 100;
+figure, imagesc(squeeze(V(x,:,:))'); axis image, colormap gray
+hold on, plot([1;size(V,2)]*ones(1,5),[1;1]*div(2:end-1),'r','LineWidth',2)
+text(size(V,2)/2*ones(1,6),midlayer,num2cell(1:6),'Color','r','FontSize',20)
+title(['Longitudinal slice, x=',num2str(x)]), xlabel('y'), ylabel('z')
+
+figure
+for i=1:numel(midlayer)
+ subplot(2,3,i), imagesc(V(:,:,midlayer(i))), axis image, colormap gray
+ title(['Region ' num2str(i), ', z=',num2str(midlayer(i))])
+ xlabel('y'), ylabel('x')
+end
+
+%% binarization
+B = medfilt3(medfilt3(V,[3,3,3])>gray_th,[3,3,3]);
+show_matvol(B)
+
+%% still figures
+figure, imagesc(squeeze(B(x,:,:))'); axis image, colormap gray
+hold on, plot([1;size(B,2)]*ones(1,5),[1;1]*div(2:end-1),'r','LineWidth',2)
+text(size(V,2)/2*ones(1,6),midlayer,num2cell(1:6),'Color','r','FontSize',20)
+title(['Longitudinal slice, x=',num2str(x)]), xlabel('y'), ylabel('z')
+
+figure
+for i=1:numel(midlayer)
+ subplot(2,3,i), imagesc(B(:,:,midlayer(i))), axis image, colormap gray
+ title(['Region ' num2str(i), ', z=',num2str(midlayer(i))])
+ xlabel('y'), ylabel('x')
+end
+
+%% distance transform
+distances = bwdist(~B) - bwdist(B) - B + 0.5; % distance from voxel center to interface
+show_matvol(distances)
+
+%% skeletonization of fibres and air
+% criterium: local maxima in any of the 3 planes for larger number of sceleton points)
+conn = cat(3,zeros(3,3),[0 1 0; 1 1 1; 0 1 0],zeros(3,3));
+skeleton = ((imregionalmax(distances,conn) | imregionalmax(distances,permute(conn,[2 3 1]))...
+ | imregionalmax(distances,permute(conn,[3 1 2]))) & B)...
+ | (imregionalmin(distances,conn) | imregionalmin(distances,permute(conn,[2 3 1]))...
+ | imregionalmin(distances,permute(conn,[3 1 2]))) & ~B;
+
+% alternative with 3D local maxima, which yields less points
+% conn = cat(3,[0 0 0; 0 1 0; 0 0 0],[0 1 0; 1 1 1; 0 1 0],[0 0 0; 0 1 0; 0 0 0]);
+% skeleton = (imregionalmax(distances,conn) & B) | (imregionalmax(distances,conn) & ~B);
+
+% cleaning
+% skeleton = skeleton & abs(distances)>1;
+skeleton([1,end],:,:) = false; skeleton(:,[1,end],:) = false; skeleton(:,:,[1,end]) = false;
+
+%% visualizing single slice with a special colormap
+z = 100; % and only a slice
+[U,cmap,cbar] = color_skeleton(B(:,:,z),skeleton(:,:,z),distances(:,:,z));
+figure, imagesc(U), colormap(cmap), axis image, title(['Slice z=',num2str(z)])
+colorbar(cbar{:})
+
+[U,cmap,cbar] = color_skeleton(squeeze(B(x,:,:))',squeeze(skeleton(x,:,:))',squeeze(distances(x,:,:))');
+figure, imagesc(U), colormap(cmap), axis image,
+title(['Longitudinal slice, x=',num2str(x)]), xlabel('y'), ylabel('z')
+colorbar(cbar{:})
+
+%% visualizing multiple slices using special colormap
+[U,cmap,cbar] = color_skeleton(B,skeleton,distances);
+figure
+for i=1:numel(midlayer)
+ subplot(2,3,i), imagesc(U(:,:,midlayer(i)),[0,255]), axis image, colormap(cmap)
+ title(['region ' num2str(i), ', z=',num2str(midlayer(i))])
+ xlabel('y'), ylabel('x')
+end
+
+%% volume slicing using special colormap
+show_matvol(U,[0,255]), colormap(cmap); % full volume
+colorbar(cbar{:})
+
+%% measuring porosity for a whole volume
+volume_all = sum(B(:))/numel(B); % phase volume per block volume
+area_all = compute_area(B)/numel(B); % phase surface per block volume
+
+% measures from distance fields
+figure
+inside_all = fit_distribution(distances(skeleton & B),true);
+a = gca; title([{'Thickness for a whole sample'};a.Title.String])
+xlim([0 30]), xlabel('Radius (pixels)'), ylabel('Count')
+
+figure
+outside_all = fit_distribution(-distances(skeleton & ~B),true);
+a = gca; title([{'Spacing for a whole sample'};a.Title.String])
+xlim([0 30]), xlabel('Radius (pixels)'), ylabel('Count')
+
+%% performing layered collection of measures
+
+volume_layer = zeros(numel(midlayer),1);
+area_layer = zeros(numel(midlayer),1);
+inside_layer = zeros(numel(midlayer),5);
+outside_layer = zeros(numel(midlayer),5);
+
+fig_in = figure;
+fig_out = figure;
+for i = 1:numel(midlayer)
+ range = ceil(div(i)):floor(div(i+1));
+ B_layer = B(:,:,range);
+ distances_layer = distances(:,:,range);
+ skeleton_layer = skeleton(:,:,range);
+
+ volume_layer(i) = sum(B_layer(:))/numel(B_layer);
+ area_layer(i) = compute_area(B_layer)/numel(B_layer);
+
+ % measures from distance fields
+ figure(fig_in), subplot(2,3,i)
+ inside_layer(i,:) = fit_distribution(distances_layer(skeleton_layer & B_layer),true);
+ legend off, xlim([0 25])
+ title(['Thickness for region ' num2str(i)])
+
+ figure(fig_out), subplot(2,3,i)
+ outside_layer(i,:) = fit_distribution(-distances_layer(skeleton_layer & ~B_layer),true);
+ legend off, xlim([0 25])
+ title(['Separation for region ' num2str(i)])
+end
+
+%% saving results
+% fid = fopen('results.tex','w');
+% fprintf(fid, 'Measure & Sample & Region 1 & Region 2 & Region 3 & Region 4 & Region 5 & Region 6 \\\\\\hline\n');
+% fprintf(fid, 'Bone volume & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f \\\\\n',[volume_all; volume_layer]);
+% fprintf(fid, 'Area & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f \\\\\n',[area_all; area_layer]);
+% fprintf(fid, 'Area per bone volume & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f \\\\\n',[area_all; area_layer]./[volume_all; volume_layer]);
+% fprintf(fid, 'Thickness, mode & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f \\\\\n',[inside_all(5); inside_layer(:,5)]);
+% fprintf(fid, 'Thickness, median & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f \\\\\n',[inside_all(1); inside_layer(:,1)]);
+% fprintf(fid, 'Thickness, mean & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f \\\\\n',[inside_all(3); inside_layer(:,3)]);
+% fprintf(fid, 'Separation, mode & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f \\\\\n',[outside_all(5); outside_layer(:,5)]);
+% fprintf(fid, 'Separation, median & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f \\\\\n',[outside_all(1); outside_layer(:,1)]);
+% fprintf(fid, 'Separation, mean & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f \\\\\n',[outside_all(3); outside_layer(:,3)]);
+% fclose(fid);
+%
+% for savefig = [3,4,6,7,9,10,11,13,14,15,16]
+% set(savefig,'Units','Normalized','Position',[0.1 0.1 0.6 0.8])
+% saveas(savefig,['Figure_',num2str(savefig),'.eps'],'epsc');
+% end
+