Skip to content
Snippets Groups Projects
Commit efcbcf19 authored by monj's avatar monj
Browse files

Added code

parent 5294bfda
Branches
No related tags found
No related merge requests found
Showing
with 1436 additions and 0 deletions
% InSegtFibre_2D
% Written by Monica Jane Emerson, April 2018.
% MIT license
%% Step 0: Add paths of GUI and other functions and scripts used in Insegt Fibre
close all,
clear all,
addpath(genpath('./texture_gui'))
addpath('./scripts')
addpath('./functions')
disp('Step 0 completed')
%% Step 1: Load and visualise the four 2D images
%the path of the 2Ddata folder in your computer
str_datafolder = '../data/2Ddata/'; %default value
%Visualise all four datasets
contents_datafolder = dir(str_datafolder);
num_datasets = length(contents_datafolder)-2;
figure
for k= 1:num_datasets
%read image
im = imread([str_datafolder,contents_datafolder(k+2).name]);
%create subtitle for image
subtitle = [num2str(k),'. ',contents_datafolder(k+2).name(1:end-4)];
ind = find(subtitle=='_');
subtitle(ind) = ' ';
%plot image
subplot(2,2,k), imagesc(im), axis image, colormap gray, title(subtitle)
end
h = msgbox('Inspect the quality of the four different scans');
waitfor(h)
disp('Step 1 completed')
%% Step 2: Select data-set and region of interest (RoI)
close all
%USER INPUT: select dataset
x = inputdlg('Choose a data-set (1, 2, 3 or 4?): ','User input',[1,20],{num2str(4)});
%read the 2D image chosen by the user
dataset = str2double(x{1});
im = imread([str_datafolder,contents_datafolder(dataset+2).name]);
%print the name of the selected data-set
disp(['Selected data-set: ', contents_datafolder(dataset+2).name(1:end-4)]);
%USER INPUT: select RoI
h = msgbox('Crop a region of interest containing 100-1000 fibres.');
waitfor(h)
im_crop = imcrop(im);
%visualise the selected RoI
if isempty(im_crop)
h = msgbox('There was an error in the cropping. Run the section again.');
waitfor(h)
else
figure, imagesc(im_crop), axis image, colormap gray,
title('Region of Interest for training the dictionary')
end
disp('Step 2 completed')
%% Step 3: Set default parameters and the patch size at the scale of the fibres
close all
%USER INPUT: Parameters that give the scale of the fibres in pixels
prompt_px = ['Pixel size [', char(181),'m]: '];
prompt_diam = ['Diameter size [', char(181),'m]: '];
x = str2double(inputdlg({prompt_px, prompt_diam},'User input',[1,30],{num2str(0.96),num2str(7)}));
%compute patch size as the closest odd integer to the fibre diameter
%measured in pixels diam/pixel_size
patch_size = calcu_psize(x(2)/x(1));
%check that the patch size is in the optimal range to ensure precision and
%speed
if ((x(2)/x(1)) < 7)
factor = 9/(x(2)/x(1));
patch_size = calcu_psize(factor*x(2)/x(1));
msg = ['We will upscale your images with a factor ',...
num2str(factor,'%0.2f'),' to obtain better precision. The patch size is now ',...
num2str(patch_size),'.'];
h = msgbox(msg);
waitfor(h)
elseif ((x(2)/x(1)) >= 17)
factor = 11/(x(2)/x(1));
patch_size = calcu_psize(factor*x(2)/x(1));
msg = ['We will downscale your images with a factor ',...
num2str(factor),' to obtain reduce the computational time. The patch size is now ',...
num2str(patch_size),'.'];
h = msgbox(msg);
waitfor(h)
else
factor = 1;
end
%set patch size
dictopt.patch_size = patch_size;
%set the default ditionary parameters
dictopt.method = 'euclidean';
dictopt.branching_factor = 3; %it should at least be 3
dictopt.number_layers = 5; %at least 5. The higher the more accurate, but also slower and more computationally heavy
num_dict_elts = calc_eltsdict(dictopt.branching_factor,dictopt.number_layers);
dictopt.number_training_patches = 5000; %at least an order of magnitude more than the number of dictionary elements
dictopt.normalization = false; %set to true if the global intensity of the slices varies along the depth of the volume
disp('Step 3 completed')
%% Step 4: Run InSegt
%display the important dictionary parameters
fprintf('\nImportant Dictionary Parameters:\n');
disp(['patch size: ', num2str(dictopt.patch_size)]);
disp(['number dictionary elements: ', num2str(calc_eltsdict(dictopt.branching_factor,dictopt.number_layers))]);
%resize cropped region of interest with factor
im_train = imresize(im_crop, factor);
%open the GUI with 2 labels for pixel annotation (fibre centre region or not)
image_texture_gui(im_train,dictopt,2)
disp('Step 4 completed')
%% Step 5: Obtain the centre coordinates from the hard segmentation
close all
%USER INPUT: label used for the centre regions
label = 2; %1 (blue) or 2 (pink)
%calculate the centroids of the connected centre regions
c = regionprops(gui_S==label,'centroid');
centrePoints = cat(1,c.Centroid);
%show fibre centre coordinates over slice
figure(1), subplot(1,2,1), imagesc(im_crop), axis image, colormap gray,
hold on, h = imagesc(ind2rgb(imresize(gui_S==label,1/factor), cool(2))); set(h, 'AlphaData',0.5);
figure(1), subplot(1,2,2), imagesc(im_crop), axis image, colormap gray
hold on, plot(centrePoints(:,1)*1/factor,centrePoints(:,2)*1/factor,'r*');
disp('Step 5 completed')
%% Step 6: Save the dictionary to process other regions of the scan
close all
%train the dictionary with the small RoI
dictionary = build_dictionary(im_train,dictopt); %create the dictionary of intensities
image_texture_gui(im_train,dictionary,2) %learn the dictionary of probabilities by annotating in the GUI
dictionary = update_dictionary(dictionary,gui_dictprob); %update dictionary to include the learnt probalilities
%USER INPUT: select a RoI to process with the learnt dictionary
h = msgbox('Crop a region of interest, it can be as big as you like.');
waitfor(h)
im_p = imcrop(im);
if isempty(im_crop)
h = msgbox('There was an error in the cropping. Run the section again from the cropping.');
waitfor(h)
else
figure, imagesc(im_p), axis image, colormap gray,
title('Region of Interest for finding fibre centres')
end
%apply the dictionary to process the larger RoI
im_process = imresize(im_p,factor);
[S,allP] = process_image(im_process,dictionary);
%calculate the centroids of the connected centre regions
label = 2; %1 (blue) or 2 (pink)
c = regionprops(S==label,'centroid');
centrePoints = cat(1,c.Centroid);
%Centre regions and coordinates over the image to evaluate the accuracy of the segmentation
figure, imagesc(im_process), axis image, colormap gray,
hold on, h = imagesc(ind2rgb(S, cool(2))); set(h, 'AlphaData',0.5);
hold on, plot(centrePoints(:,1),centrePoints(:,2),'y.');
title('Segmented centre regions and fibre detections over the processed image');
if factor ~=1
figure, imagesc(im_p), axis image, colormap gray,
hold on, plot(centrePoints(:,1)*1/factor,centrePoints(:,2)*1/factor,'y.');
title('Centre detections over the original resolution')
end
disp('Step 6 completed')
%% Step 7: Tune the threshold applied to the probability map
close all
%visualise the propabilistic segmentation of fibre centre regions
figure(1), imagesc(allP(:,:,2)), axis image off, colormap gray, title('Centre class probability map'), colorbar
y = 'yes';
while strcmp(y,'yes')
%USER INPUT: threshold value for the probability map
x = inputdlg('Choose threshold: ','User input',[1,20],{num2str(0.5)});
%Centre regions over the image to evaluate the accuracy of the segmentation
figure, imagesc(im_process), axis image, colormap gray,
hold on, h = imagesc(ind2rgb(allP(:,:,2)>str2double(x{1}), cool(2))); set(h, 'AlphaData',0.5);
%USER INPUT: do you want to try another threshold??
y = questdlg('Do you want to try another thresold?: ','User input','yes','no','yes');
end
disp('Step 7 completed')
%% Step 8: Save your manual input to add additional markings later on
close all
%annotate a little in the GUI and save your manual input
h = msgbox(['Annotate a little in the GUI and press [S] before closing the GUI (then close by pressing export [E]). Choose to save',...
' your annotations in the folder "InsegtFibre_workshopDTU/data/saved_data/trainingData".']);
waitfor(h)
image_texture_gui(im_train,dictopt,2)
%save the training region, over which you placed the markings
save('../data/saved_data/trainingData/training_image.mat','im_train');
%Call the GUI with a saved labelling
load('../data/saved_data/trainingData/training_image.mat');
[filename,pathname] = uigetfile('*_labels_indexed.png', 'Choose label image', '../data/saved_data/trainingData/');
labelling = imread([pathname,filename]); %read labelling
image_texture_gui(im_train,dictopt,2,labelling) %call the GUI with the saved labelling
disp('Step 8 completed')
% InSegtFibre_3D
% Written by Monica Jane Emerson, May 2018.
% MIT license
%% Step 0: Add paths of Insegt GUI, other functions and scripts, and prepare workspace
close all,
clear all,
addpath(genpath('./texture_gui'))
addpath('./scripts')
addpath('./functions')
disp('Step 0 completed')
%% Step 1: Locate folder with the CT data
close all
indicate_dataFolder
disp('Step 1 completed')
%% Step 2: Load and visualise the data
close all
[V,depth] = loadVolumeRoI(path_volumeFolder);
figure, imshow3D(V) %visualise the loaded RoI
disp('Step 2 completed')
%% Step 3: Compute a dictionary model from the CT scan
close all
setup_dict %set the parameters of the dictionary model
create_dict %learn the dictionary model from training data
disp('Step 3 completed')
%% Step 4: Obtain the centre-class probabilistic segmentation for a slice
close all
% USER INPUT: Load a saved dictionary
[filename,pathname] = uigetfile('*.mat','Choose dictionary','../data/saved_data/dictionaries/dictionary');
load([pathname,filename]);
% USER INPUT: Choose slice
prompt = ['Choose a slice to segment: [1,', num2str(size(V,3)),']:'];
x = inputdlg(prompt,'User input',[1,30],{num2str(round(size(V,3)/4))});
sliceAnalys = str2double(x{1});
h = msgbox('Segmenting image...');
[~,allP] = process_image(V(:,:,sliceAnalys),dictionary);
delete(h)
Pcentre = allP(:,:,2);
%Visualise for each pixel its probabily of belonging to a fibre centre region
h = figure(1), %h.Units = 'normalized'; h.Position = [0 0 1 1];
subplot(1,2,1), imagesc(Pcentre), axis image, colormap gray, colorbar,
title('Probability of belonging to the central region of a fibre')
figure(1), subplot(1,3,3), histogram(Pcentre(:)),
ylabel('Count of probabilities'), xlabel('Range of probabilities')
disp('Step 4 completed')
%% Step 5: Obtain the fibre centres by thresholding the probabilistic segmentation
% USER INPUT: threshold value for the probability map
thresh = inputdlg('Choose threshold: ','User input',[1,20],{num2str(0.50)});
%Threshold probability map
areas_th = Pcentre > str2double(thresh{1});
%Calculate centroids
concomp = logical(areas_th);
c = regionprops(concomp,'centroid');
centrePoints_2D = cat(1,c.Centroid);
%Visualisations to evaluate the chosen threshold
figure(2), subplot(1,2,1), imagesc(V(:,:,sliceAnalys)), axis image, colormap gray,
hold on, h = imagesc(ind2rgb(areas_th, cool(2))), set(h, 'AlphaData',0.5);
figure(2), subplot(1,2,2), imagesc(V(:,:,sliceAnalys)), axis image, colormap gray
hold on, plot(centrePoints_2D(:,1),centrePoints_2D(:,2),'r*');
disp('Step 5 completed')
%% Step 6: Detect fibre centres in the batch of slices forming the volume
close all
clear centrePoints, clear allP, clear c
% USER INPUT: Load a saved dictionary
[filename,pathname] = uigetfile('*.mat','Choose dictionary','../data/saved_data/dictionaries/dictionary_');
load([pathname,filename]);
% USER INPUT: threshold value for the probability map
thresh = inputdlg('Choose threshold: ','User input',[1,20],{num2str(0.5)});
%Initialise variables
step = 1;
slices = 1:step:size(V,3);
%Start timer
tstart = tic;
h = waitbar(0,'Processing volume...');
for l = 1:numel(slices)
%Compute the probabilistic segmentation
[~,allP] = process_image(V(:,:,slices(l)),dictionary);
%Save probability map
P(:,:,l) = allP(:,:,2);
%Compute the connected components of the thresholded probabilities
c = regionprops(bwlabel(allP(:,:,2) > str2num(thresh{1})),'centroid');
%Store the coordinates of the centre points
centrePoints{l} = cat(1,c.Centroid);
%Update the data processing wait bar
waitbar(l/numel(slices),h)
end
%disp processing time
t_load = toc(tstart);
f = msgbox(['time to process volume: ' , num2str(t_load), 's']);
close(h)
waitfor(f)
%Compute the number of detected fibres per slice
clear len
close all
for l = 1:length(centrePoints)
len(l) = size(centrePoints{l},1);
end
aux = str2double(depth);
figure(1), plot([aux(1):aux(3):aux(2)],len,'.'), axis tight, title('Nr. fibres per slice'),
xlabel('slice number (1 is the top of the sample)')
ylabel('number of detected fibres')
%Visualise the detections in 3D
if ~exist('vx_size')
vx_size = str2double(inputdlg(['Pixel size [', char(181),'m]: '],'User input',[1,30],{num2str(1.1)}));
end
for l = 1:step:length(centrePoints)
figure(2), hold on, plot3(centrePoints{l}(:,1)*vx_size,centrePoints{l}(:,2)*vx_size,l*ones(size(centrePoints{l}(:,1)))*vx_size,'r.','MarkerSize',1)
end
h = figure(2), axis tight, view(3), axis equal
xlabel('x axis [\mum]'),ylabel('y axis [\mum]'), zlabel('z axis [\mum]')
aux = cellstr(num2str(vx_size*linspace(str2double(depth{1}),str2double(depth{2}),length(h.Children.ZTickLabel)+1)','%.f'));
h.Children.ZTickLabel = aux(2:end);
% USER INPUT: Save 3D detections
uisave('centrePoints','../data/saved_data/centreCoordinates/centreCoords_?.mat')
disp('Step 6 completed')
%% Step 7: Compute the 3D fibre trajectories via tracking
% USER INPUT: Load the saved 3D detections
[filename,pathname] = uigetfile('*.mat','Choose the set of 3D coordinates to track','../data/saved_data/centreCoordinates');
load([pathname,filename]);
% USER INPUT: How many pixels is a fibre allowed to move from one slice to the
% next? (function of the step, the diameter and the pixel size).
x = inputdlg('How many pixels is a fibre allowed to move from one slice to the next? (a function of the step, the diameter and the pixel size): ','User input',[1,65],{num2str(5)});
dist = str2num(x{1});
% USER INPUT: Should the tracking start with the detections of the bottom
% or the top slice?
if ~exist('depth')
depth = input_depthRoI(num_slices);
end
top = 'Top of the volume (slice 1)';
bottom = ['Bottom of the volume (slice ', depth{2},')'];
direct = questdlg('Where shall we start the fibres?','Tracking direction',top, bottom,bottom);
if strcmp(direct,bottom)
m = length(centrePoints);
else
m = 1;
end
%Initialise variables
thisCt = centrePoints{m}; %we start tracking from the bottom of the volume (end slice)
f_x = nan(size(thisCt,1),length(centrePoints));
f_y = nan(size(thisCt,1),length(centrePoints));
f_x(:,1) = centrePoints{m}(:,1);
f_y(:,1) = centrePoints{m}(:,2);
%Starting time
tstart = tic;
%Track the coordinates in the z direction
for l = 2:length(centrePoints)
% coordinates of next slice
if strcmp(direct,bottom)
ctNext = centrePoints{end-l+1};
else
ctNext = centrePoints{l};
end
% build kd-tree
tree_down = KDTreeSearcher(ctNext);
% search the kd-tree
[id_down, d_down] = knnsearch(tree_down,thisCt,'K',1);
% build kd-tree
tree_up = KDTreeSearcher(thisCt);
% search the kd-tree
[id_up, d_up] = knnsearch(tree_up,ctNext,'K',1);
% keep matches that are bidirectional
id_top = unique(id_up(id_down));
id_bottom = id_down(unique(id_up(id_down)));
% keep matches that are closer than a distance
idMatch = find(d_down(id_top) < dist);
id_top_th = id_top(idMatch);
id_bottom_th = id_bottom(idMatch);
if strcmp(direct,bottom)
f_x(id_top_th,l) = centrePoints{end-l+1}(id_bottom_th,1);
f_y(id_top_th,l) = centrePoints{end-l+1}(id_bottom_th,2);
else
f_x(id_top_th,l) = centrePoints{l}(id_bottom_th,1);
f_y(id_top_th,l) = centrePoints{l}(id_bottom_th,2);
end
%points to be matched with the next slice
thisCt(id_top_th,:) = ctNext(id_bottom_th,:);
end
t_load = toc(tstart);
disp(['time to track coordinates along the volume: ' , num2str(t_load), 's']);
%Save track coordinates in the structure 'fibres'
fibres.x = f_x;
fibres.y = f_y;
%Visualise fibre tracks
% USER INPUT: How many fibres would you like to visualise (in percentage)?
questions = {'What percentage of the tracked fibres would you like to visualise (%)?: '};
x = inputdlg(questions,'User input',[1,40],{num2str(10)}); %ask for patch size
perc = str2double(x{1})/100;
fib_ind = randperm(size(fibres.x,1),round(perc*size(fibres.x,1))); %indices of fibre to plot
z = [1:size(fibres.x,2)]; %z-values
if ~exist('vx_size')
vx_size = str2double(inputdlg(['Pixel size [', char(181),'m]: '],'User input',[1,30],{num2str(1.1)}));
end
figure(3),clf
for l= 1:numel(fib_ind)
h = figure(3);
hold on, plot3(fibres.x(fib_ind(l),:)*vx_size,fibres.y(fib_ind(l),:)*vx_size,z*vx_size,'LineWidth',1),
end
axis tight, axis equal, view(3),
xlabel('x axis [\mum]'),ylabel('y axis [\mum]'), zlabel('z axis [\mum]')
aux = cellstr(num2str(vx_size*linspace(str2double(depth{1}),str2double(depth{2}),length(h.Children.ZTickLabel)+1)','%.f'));
h.Children.ZTickLabel = aux(2:end);
% USER INPUT: Save fibre tracks
if strcmp(direct,bottom)
direction = [depth{2},'to',depth{1}];
else
direction = [depth{1},'to',depth{2}];
end
uisave('fibres',['../data/saved_data/fibreTracks/fibreTracks_dataSetX_',direction,'step',depth{3},'.mat'])
\ No newline at end of file
% InSegtFibre_orientations
% Written by Monica Jane Emerson, May 2019.
% MIT license
%% Prepare workspace
close all
clear all
addpath('./scripts')
addpath('./functions')
%% Read tracks, z-values (analysed slices) and voxel size
[filename_tracks,pathname_tracks] = uigetfile('*.mat','Choose tracks','../data/saved_data/fibreTracks');
fibres = getfield(load([pathname_tracks,filename_tracks]),'fibres');
depth = input_depthRoI; %get z coordinates (start, end and step)
z_axis = [str2double(depth{1}):str2double(depth{3}):str2double(depth{2})];
vx_size = str2double(inputdlg(['What is the pixel size of your scan [', char(181),'m]: '],'Pixel/voxel size',[1,45],{num2str(1.1)}));
%% Display tracks
h1 = figure(1);
for l= 1:size(fibres.x,1)
figure(1), hold on,
plot3(fibres.x(l,:)*vx_size,fibres.y(l,:)*vx_size,z_axis*vx_size,'k','LineWidth',1),
end
figure(1), axis tight, axis equal, view(2), view(3)
xlabel('x axis [\mum]'),ylabel('y axis [\mum]'), zlabel('z axis [\mum]')
%% Display locations of missed detections
missed_det = isnan(fibres.x');
figure(2), imagesc(missed_det), axis image, axis xy
title('Missed fibre centres')
ylabel('z axis [slices]'),xlabel('fibre ID'),
%% Select fibres for the orientation plot
%Indices of fibres that miss start and/or end detection
%(will be excluded from the orientation plot)
exc_fibres = find((missed_det(1,:)|missed_det(end,:))==1);
%Indices of the fibre kept for the orientation plot
keep_fibres = 1:size(fibres.x,1);
keep_fibres(exc_fibres) = [];
display(['The orientation can been computed for ',num2str(length(keep_fibres)/size(fibres.x,1)*100,'%.2f') ,'% of the tracked fibres']);
%% Compute orientation (inclination and it's direction, azimuth)
%displacement between end and start points
disp_x = fibres.x(keep_fibres,end)-fibres.x(keep_fibres,1);
disp_y = fibres.y(keep_fibres,end)-fibres.y(keep_fibres,1);
disp_r = sqrt(disp_x.^2 + disp_y.^2);
%inclination and azimuth angles
inclination = atan2(disp_r,(size(fibres.x,2)-1)*ones(length(keep_fibres),1));
azimuth = atan2(disp_y,disp_x);
display('The orientations (inclination and azimuth angles) have been computed');
%% Histograms of orientation angles
close all
%Inclination
max_inc = str2double(inputdlg('What is the maximum misalignment angle you expect [0, 90]?','Max. Inclination',[1,35],{num2str(15)}));
bins_inc = [0:max_inc];
int_inc = [0:0.1:max_inc];
h_inc = hist(inclination*180/pi,bins_inc)/length(inclination);
h_inc_int = interp1(bins_inc,h_inc,int_inc);
h1 = figure(1),hold on, plot(int_inc,smooth(h_inc_int))
ylim([0,max(smooth(h_inc_int))]), xlim([0,max_inc]),h1.Children.XTick = linspace(0,max_inc,10);
xlabel('Inclination \phi [\circ] of the fibres'), ylabel('Relative frequency')
title('Amount of fibre misalignment with respect to the scanning z-axis'),
%Azimuth
max_az = 360;
bins_az = [-max_az/2:24:max_az/2];
int_az = [-max_az/2:2.4:max_az/2];
h_az = hist(azimuth*180/pi,bins_az)/length(azimuth);
h_az_int = interp1(bins_az,h_az,int_az);
h2 = figure(2),hold on, plot(int_az,smooth(h_az_int))
xlim([-max_az/2,max_az/2]), h2.Children.XTick = [-180,-90,0,90,180];
ylim([0,max(smooth(h_az_int))])
xlabel('Azimuth \theta [\circ] of the fibres'), ylabel('Relative frequency')
title('Direction of misalignment with respect to scanning x-axis'),
%% Colour-code fibre tracks according to orientation angles
%Azimuth (direction of misalingment)
cmap_az = circshift(colormap(hsv),32);
h3 = figure(3);
for l= 1:length(keep_fibres)
colouring = cmap_az(round(azimuth(l)/(2*pi)*63+1+63/2),:);
figure(3), hold on, plot3(fibres.x(keep_fibres(l),:)*vx_size,fibres.y(keep_fibres(l),:)*vx_size,z_axis*vx_size,'LineWidth',1,'Color',colouring),
end
figure(3), axis tight, axis equal, view(2),
xlabel('x axis [\mum]'),ylabel('y axis [\mum]'), zlabel('z axis [\mum]')
colormap(cmap_az), colorbar, caxis([-180,180]), h3.Children(1).Ticks = -180:30:180;
for k = 1:length(h3.Children(1).Ticks)
h3.Children(1).TickLabels(k) = {[h3.Children(1).TickLabels{k},'']};
end
%Inclination (amount of misalignment)
cmap_inc = linspace(0,1,32)'*[1,1,1];
sat_angle = max_inc/3;
sat = sat_angle*pi/180;
h4 = figure(4);
for l= 1:length(keep_fibres)
if inclination(l) < sat
colouring = [1,1,1]*inclination(l)/sat;
else
colouring = [1,1,1];
end
figure(4), hold on, plot3(fibres.x(keep_fibres(l),:)*vx_size,fibres.y(keep_fibres(l),:)*vx_size,z_axis*vx_size,'LineWidth',1,'Color',colouring),
end
figure(4), axis tight, axis equal, view(2),
xlabel('x axis [\mum]'),ylabel('y axis [\mum]'), zlabel('z axis [\mum]')
colormap(cmap_inc),
colorbar, caxis([0,sat_angle]), h.Children(1).Ticks = linspace(0,sat_angle,5);
h4.Children(1).TickLabels(end) = {['>',num2str(sat_angle)]};
for k = 1:length(h.Children(1).Ticks)
h4.Children(1).TickLabels(k) = {[h4.Children(1).TickLabels{k},'']};
end
%Azimuth weighed by inclination (amount of misalignment)
h5 = figure(5);
for l= 1:length(keep_fibres)
if inclination(l) < sat
colouring = cmap_az(round(azimuth(l)/(2*pi)*63+1+63/2),:)*inclination(l)/sat;
else
colouring = cmap_az(round(azimuth(l)/(2*pi)*63+1+63/2),:);
end
figure(5), hold on, plot3(fibres.x(keep_fibres(l),:)*vx_size,fibres.y(keep_fibres(l),:)*vx_size,z_axis*vx_size,'LineWidth',1,'Color',colouring),
end
figure(5), axis tight, axis equal, view(2),
xlabel('x axis [\mum]'),ylabel('y axis [\mum]'), zlabel('z axis [\mum]')
colormap(cmap_az), colorbar, caxis([-180,180]), h5.Children(1).Ticks = -180:30:180;
for k = 1:length(h.Children(1).Ticks)
h5.Children(1).TickLabels(k) = {[h5.Children(1).TickLabels{k},'']};
end
% VALIDATION_3DFIBREDETECTION
% Written by Monica Jane Emerson, April 2019.
% MIT license
%% Add paths of Insegt, functions and scripts, and prepare workspace
close all, clear all,
addpath('./scripts')
addpath('./functions')
%% Get the measured fibre tracks in the format required by the validation GUI
% USER INPUT: Load the fibre tracks to validate
[filename_tracks,pathname_tracks] = uigetfile('*.mat','Choose the fibre tracks to validate','saved_data/fibreTracks');
load([pathname_tracks,filename_tracks]);
%USER INPUT: Deregister the tracks if necessary
x = questdlg('Have the tracks undergone a transformation (rotation and/or translation) since they were measured from the CT scan?','2D de-registration','Yes','No','No');
if strcmp(x,'Yes')
[filename,pathname] = uigetfile('*.mat','Load the transformation matrix','saved_data/2Dregistration/');
t = load([pathname,filename],'tform');
%USER INPUT: Offset between registration FoV and analysis FoV
offset = str2double(inputdlg('Offset between the FoV used for the analysis and the FoV used for the 2D registration','Offset registration-analysis',[1,50],{num2str(100)}));
aux = [fibres.x(:) + offset*ones(size(fibres.x(:))),fibres.y(:) + offset*ones(size(fibres.x(:))),ones(size(fibres.y(:)))]*t.tform.tdata.Tinv;
fibres.x = reshape(aux(:,1), size(fibres.x)) - offset*ones(size(fibres.x));
fibres.y = reshape(aux(:,2), size(fibres.x)) - offset*ones(size(fibres.x));
end
for l = 1:size(fibres.x,2)
centers{l} = [fibres.x(:,l),fibres.y(:,l),[1:size(fibres.x,1)]'];
end;
%% Get the analysed volume of interest (VoI) in the format read by the validaiton GUI
close all
% USER INPUT: Option to open a volume of interest saved as a tif file
x = questdlg('Is the analysed volume of interest saved as one .tif file, with the slices ordered in the direction of tracking?','Load analysed volume','Yes','No','Yes');
if strcmp(x,'No')
indicate_dataFolder;
[V,depth] = loadVolumeRoI(path_volumeFolder);
% USER INPUT: Tracking direction
top = ['Top to bottom (from slice ',depth{1},')'];
bottom = ['Bottom to top (from slice ',depth{2},')'];
direct = questdlg('In which direction did you track the coordinates?','Tracking direction',top, bottom, bottom);
% USER INPUT: Give a name to the .tif Volume of Interest
if strcmp(direct,bottom)
direction = [depth{2},'to',depth{1},'step',depth{3}];
else
direction = [depth{1},'to',depth{2},'step',depth{3}];
end
name_tif = inputdlg('Name for the .tif file of your volume of interest (VoI): ','Name your VOI',[1,55],{['V_?_',direction,'.tif']});
filename_voi= name_tif{1};
pathname_voi = 'saved_data/volumeOfInterest/';
if strcmp(direct,bottom)
imwrite(V(:,:,end),[pathname_voi,filename_voi]);
for l= 2:size(V,3)
imwrite(V(:,:,end-l+1),[pathname_voi,filename_voi], 'WriteMode','append');
end
else
imwrite(V(:,:,1),[pathname_voi,filename_voi]);
for l= 2:size(V,3)
imwrite(V(:,:,l),[pathname_voi,filename_voi], 'WriteMode','append');
end
end
else
[filename_voi,pathname_voi] = uigetfile('*.tif','Choose volume of interest','saved_data/volumeOfInterest/VOI_v');
end
%% Validate all fibres
track_slicer([pathname_voi,filename_voi],centers)
function [num] = calc_eltsdict(k,n)
%CALC_ELTSDICT Summary of this function goes here
% Written by Monica Jane Emerson, April 2018, MIT License
vect = [1:n];
series = k.^(vect-1);
num = k*sum(series);
end
function [patch_size] = calcu_psize(pixPerDiam)
%CALCU_PSIZE Loads a selected region of interest from a volume.
% Details...
%by Monica Jane Emerson, April 2019, MIT License
if mod(round(pixPerDiam),2)==0
if (round(pixPerDiam)-pixPerDiam)<0
patch_size = round(pixPerDiam)+1;
else
patch_size = round(pixPerDiam)-1;
end
else
patch_size = round(pixPerDiam);
end
end
function imshow3D( Img, disprange )
%IMSHOW3D displays 3D grayscale images in slice by slice fashion with mouse
%based slice browsing and window and level adjustment control.
%
% Usage:
% imshow3D ( Image )
% imshow3D ( Image , [] )
% imshow3D ( Image , [LOW HIGH] )
%
% Image: 3D image MxNxK (K slices of MxN images)
% [LOW HIGH]: display range that controls the display intensity range of
% a grayscale image (default: the widest available range)
%
% Use the scroll bar or mouse scroll wheel to switch between slices. To
% adjust window and level values keep the mouse right button pressed and
% drag the mouse up and down (for level adjustment) or right and left (for
% window adjustment).
%
% "Auto W/L" button adjust the window and level automatically
%
% While "Fine Tune" check box is checked the window/level adjustment gets
% 16 times less sensitive to mouse movement, to make it easier to control
% display intensity rang.
%
% Note: The sensitivity of mouse based window and level adjustment is set
% based on the user defined display intensity range; the wider the range
% the more sensitivity to mouse drag.
%
%
% Example
% --------
% % Display an image (MRI example)
% load mri
% Image = squeeze(D);
% figure,
% imshow3D(Image)
%
% % Display the image, adjust the display range
% figure,
% imshow3D(Image,[20 100]);
%
% See also IMSHOW.
%
% - Maysam Shahedi (mshahedi@gmail.com)
% - Released: 1.0.0 Date: 2013/04/15
% - Revision: 1.1.0 Date: 2013/04/19
%
sno = size(Img,3); % number of slices
S = round(sno/2);
global InitialCoord;
MinV = 0;
MaxV = max(Img(:));
LevV = (double( MaxV) + double(MinV)) / 2;
Win = double(MaxV) - double(MinV);
WLAdjCoe = (Win + 1)/1024;
FineTuneC = [1 1/16]; % Regular/Fine-tune mode coefficients
if isa(Img,'uint8')
MaxV = uint8(Inf);
MinV = uint8(-Inf);
LevV = (double( MaxV) + double(MinV)) / 2;
Win = double(MaxV) - double(MinV);
WLAdjCoe = (Win + 1)/1024;
elseif isa(Img,'uint16')
MaxV = uint16(Inf);
MinV = uint16(-Inf);
LevV = (double( MaxV) + double(MinV)) / 2;
Win = double(MaxV) - double(MinV);
WLAdjCoe = (Win + 1)/1024;
elseif isa(Img,'uint32')
MaxV = uint32(Inf);
MinV = uint32(-Inf);
LevV = (double( MaxV) + double(MinV)) / 2;
Win = double(MaxV) - double(MinV);
WLAdjCoe = (Win + 1)/1024;
elseif isa(Img,'uint64')
MaxV = uint64(Inf);
MinV = uint64(-Inf);
LevV = (double( MaxV) + double(MinV)) / 2;
Win = double(MaxV) - double(MinV);
WLAdjCoe = (Win + 1)/1024;
elseif isa(Img,'int8')
MaxV = int8(Inf);
MinV = int8(-Inf);
LevV = (double( MaxV) + double(MinV)) / 2;
Win = double(MaxV) - double(MinV);
WLAdjCoe = (Win + 1)/1024;
elseif isa(Img,'int16')
MaxV = int16(Inf);
MinV = int16(-Inf);
LevV = (double( MaxV) + double(MinV)) / 2;
Win = double(MaxV) - double(MinV);
WLAdjCoe = (Win + 1)/1024;
elseif isa(Img,'int32')
MaxV = int32(Inf);
MinV = int32(-Inf);
LevV = (double( MaxV) + double(MinV)) / 2;
Win = double(MaxV) - double(MinV);
WLAdjCoe = (Win + 1)/1024;
elseif isa(Img,'int64')
MaxV = int64(Inf);
MinV = int64(-Inf);
LevV = (double( MaxV) + double(MinV)) / 2;
Win = double(MaxV) - double(MinV);
WLAdjCoe = (Win + 1)/1024;
elseif isa(Img,'logical')
MaxV = 0;
MinV = 1;
LevV =0.5;
Win = 1;
WLAdjCoe = 0.1;
end
SFntSz = 9;
LFntSz = 10;
WFntSz = 10;
LVFntSz = 9;
WVFntSz = 9;
BtnSz = 10;
ChBxSz = 10;
if (nargin < 2)
[Rmin Rmax] = WL2R(Win, LevV);
elseif numel(disprange) == 0
[Rmin Rmax] = WL2R(Win, LevV);
else
LevV = (double(disprange(2)) + double(disprange(1))) / 2;
Win = double(disprange(2)) - double(disprange(1));
WLAdjCoe = (Win + 1)/1024;
[Rmin Rmax] = WL2R(Win, LevV);
end
axes('position',[0,0.2,1,0.8]), imshow(Img(:,:,S), [Rmin Rmax])
FigPos = get(gcf,'Position');
S_Pos = [50 45 uint16(FigPos(3)-100)+1 20];
Stxt_Pos = [50 65 uint16(FigPos(3)-100)+1 15];
Wtxt_Pos = [50 20 60 20];
Wval_Pos = [110 20 60 20];
Ltxt_Pos = [175 20 45 20];
Lval_Pos = [220 20 60 20];
BtnStPnt = uint16(FigPos(3)-250)+1;
if BtnStPnt < 300
BtnStPnt = 300;
end
Btn_Pos = [BtnStPnt 20 100 20];
ChBx_Pos = [BtnStPnt+110 20 100 20];
if sno > 1
shand = uicontrol('Style', 'slider','Min',1,'Max',sno,'Value',S,'SliderStep',[1/(sno-1) 10/(sno-1)],'Position', S_Pos,'Callback', {@SliceSlider, Img});
stxthand = uicontrol('Style', 'text','Position', Stxt_Pos,'String',sprintf('Slice# %d / %d',S, sno), 'BackgroundColor', [0.8 0.8 0.8], 'FontSize', SFntSz);
else
stxthand = uicontrol('Style', 'text','Position', Stxt_Pos,'String','2D image', 'BackgroundColor', [0.8 0.8 0.8], 'FontSize', SFntSz);
end
ltxthand = uicontrol('Style', 'text','Position', Ltxt_Pos,'String','Level: ', 'BackgroundColor', [0.8 0.8 0.8], 'FontSize', LFntSz);
wtxthand = uicontrol('Style', 'text','Position', Wtxt_Pos,'String','Window: ', 'BackgroundColor', [0.8 0.8 0.8], 'FontSize', WFntSz);
lvalhand = uicontrol('Style', 'edit','Position', Lval_Pos,'String',sprintf('%6.0f',LevV), 'BackgroundColor', [1 1 1], 'FontSize', LVFntSz,'Callback', @WinLevChanged);
wvalhand = uicontrol('Style', 'edit','Position', Wval_Pos,'String',sprintf('%6.0f',Win), 'BackgroundColor', [1 1 1], 'FontSize', WVFntSz,'Callback', @WinLevChanged);
Btnhand = uicontrol('Style', 'pushbutton','Position', Btn_Pos,'String','Auto W/L', 'FontSize', BtnSz, 'Callback' , @AutoAdjust);
ChBxhand = uicontrol('Style', 'checkbox','Position', ChBx_Pos,'String','Fine Tune', 'BackgroundColor', [0.8 0.8 0.8], 'FontSize', ChBxSz);
set (gcf, 'WindowScrollWheelFcn', @mouseScroll);
set (gcf, 'ButtonDownFcn', @mouseClick);
set(get(gca,'Children'),'ButtonDownFcn', @mouseClick);
set(gcf,'WindowButtonUpFcn', @mouseRelease)
set(gcf,'ResizeFcn', @figureResized)
% -=< Figure resize callback function >=-
function figureResized(object, eventdata)
FigPos = get(gcf,'Position');
S_Pos = [50 45 uint16(FigPos(3)-100)+1 20];
Stxt_Pos = [50 65 uint16(FigPos(3)-100)+1 15];
BtnStPnt = uint16(FigPos(3)-250)+1;
if BtnStPnt < 300
BtnStPnt = 300;
end
Btn_Pos = [BtnStPnt 20 100 20];
ChBx_Pos = [BtnStPnt+110 20 100 20];
if sno > 1
set(shand,'Position', S_Pos);
end
set(stxthand,'Position', Stxt_Pos);
set(ltxthand,'Position', Ltxt_Pos);
set(wtxthand,'Position', Wtxt_Pos);
set(lvalhand,'Position', Lval_Pos);
set(wvalhand,'Position', Wval_Pos);
set(Btnhand,'Position', Btn_Pos);
set(ChBxhand,'Position', ChBx_Pos);
end
% -=< Slice slider callback function >=-
function SliceSlider (hObj,event, Img)
S = round(get(hObj,'Value'));
set(get(gca,'children'),'cdata',Img(:,:,S))
caxis([Rmin Rmax])
if sno > 1
set(stxthand, 'String', sprintf('Slice# %d / %d',S, sno));
else
set(stxthand, 'String', '2D image');
end
end
% -=< Mouse scroll wheel callback function >=-
function mouseScroll (object, eventdata)
UPDN = eventdata.VerticalScrollCount;
S = S - UPDN;
if (S < 1)
S = 1;
elseif (S > sno)
S = sno;
end
if sno > 1
set(shand,'Value',S);
set(stxthand, 'String', sprintf('Slice# %d / %d',S, sno));
else
set(stxthand, 'String', '2D image');
end
set(get(gca,'children'),'cdata',Img(:,:,S))
end
% -=< Mouse button released callback function >=-
function mouseRelease (object,eventdata)
set(gcf, 'WindowButtonMotionFcn', '')
end
% -=< Mouse click callback function >=-
function mouseClick (object, eventdata)
MouseStat = get(gcbf, 'SelectionType');
if (MouseStat(1) == 'a') % RIGHT CLICK
InitialCoord = get(0,'PointerLocation');
set(gcf, 'WindowButtonMotionFcn', @WinLevAdj);
end
end
% -=< Window and level mouse adjustment >=-
function WinLevAdj(varargin)
PosDiff = get(0,'PointerLocation') - InitialCoord;
Win = Win + PosDiff(1) * WLAdjCoe * FineTuneC(get(ChBxhand,'Value')+1);
LevV = LevV - PosDiff(2) * WLAdjCoe * FineTuneC(get(ChBxhand,'Value')+1);
if (Win < 1)
Win = 1;
end
[Rmin, Rmax] = WL2R(Win,LevV);
caxis([Rmin, Rmax])
set(lvalhand, 'String', sprintf('%6.0f',LevV));
set(wvalhand, 'String', sprintf('%6.0f',Win));
InitialCoord = get(0,'PointerLocation');
end
% -=< Window and level text adjustment >=-
function WinLevChanged(varargin)
LevV = str2double(get(lvalhand, 'string'));
Win = str2double(get(wvalhand, 'string'));
if (Win < 1)
Win = 1;
end
[Rmin, Rmax] = WL2R(Win,LevV);
caxis([Rmin, Rmax])
end
% -=< Window and level to range conversion >=-
function [Rmn Rmx] = WL2R(W,L)
Rmn = L - (W/2);
Rmx = L + (W/2);
if (Rmn >= Rmx)
Rmx = Rmn + 1;
end
end
% -=< Window and level auto adjustment callback function >=-
function AutoAdjust(object,eventdata)
Win = double(max(Img(:))-min(Img(:)));
Win (Win < 1) = 1;
LevV = double(min(Img(:)) + (Win/2));
[Rmin, Rmax] = WL2R(Win,LevV);
caxis([Rmin, Rmax])
set(lvalhand, 'String', sprintf('%6.0f',LevV));
set(wvalhand, 'String', sprintf('%6.0f',Win));
end
end
% -=< Maysam Shahedi (mshahedi@gmail.com), April 19, 2013>=-
\ No newline at end of file
function [V,depth] = loadVolumeRoI(path)
%LOADVOLUMEROI Loads a selected region of interest from a volume.
% Details...
%by Monica Jane Emerson, March 2019, MIT License
contents_datafolder = dir(path);
path_firstSlice = [path,contents_datafolder(3).name];
%select FoV in the cross-sectional direction
[rect,~] = selectFoV(path_firstSlice);
%select slices
num_slices = length(contents_datafolder)-2; %depth of the volume
zmin = ['Lower limit in the depth direction [1, ',num2str(num_slices),']:'];
zmax = ['Upper limit in the depth direction [1, ',num2str(num_slices),']:'];
zstep = ['Step (1 to take every slice):'];
prompt = {zmin,zmax,zstep};
depth = inputdlg(prompt,'RoI in the depth direction', [1,50], {'1',num2str(num_slices),'1'});
slices = [str2double(depth{1}):str2double(depth{3}):str2double(depth{2})];
V = zeros([rect(4)-rect(2)+1,rect(3)-rect(1)+1,numel(slices)]); %volume preallocation (for speed)
tstart = tic;
h = waitbar(0,'Loading volume...', 'CreateCancelBtn','setappdata(gcbf,''canceling'',1)');
for k= 1:numel(slices)
im = imread([path,contents_datafolder(slices(k)+2).name]);
V(:,:,k) = double(im(rect(2):rect(4),rect(1):rect(3)))/2^16;
%Processing data wait bar
waitbar(k/numel(slices),h)
if getappdata(h,'canceling')
break
end
if k == numel(slices)
delete(h)
t_load = toc(tstart);
f = msgbox(['time to load volume: ' , num2str(t_load), 's']);
end
end
end
function [num] = num_dictAtoms(k,n)
%NUM_DICTATOMS Summary of this function goes here
%by Monica Jane Emerson, April 2018, MIT License
vect = [1:n];
series = k.^(vect-1);
num = k*sum(series);
end
function [rect,im_rect] = selectFoV(slice_path)
%SELECTFOV is a user-friendly function to select a region of interest from a 2D image.
% The user can choose to select interactively the region of interest (RoI)
% interactively with imcrop, or specify the coordinates of the RoI.
% by Monica Jane Emerson, March 2019, MIT License
%USER INPUT: SLICE_PATH, The location of the slice to crop the RoI from.
% Other user input given through interaction with windows.
autoManual = questdlg('How would you like to select the region of interest (RoI)?','Select analysis region','Interactively','I know the coordinates','Interactively');
im = imread(slice_path);
if strcmp(autoManual,'Interactively')
[~,rect] = imcrop(im);
rect = round(rect);
rect(3:4) = rect(1:2) + rect(3:4) - [1,1];
else
xmin = ['Lower limit for the x value [1, ',num2str(size(im,2)),']:'];
xmax = ['Upper limit for the x value [1, ',num2str(size(im,2)),']:'];
ymin = ['Lower limit for the y value [1, ',num2str(size(im,1)),']:'];
ymax = ['Upper limit for the y value [1, ',num2str(size(im,1)),']:'];
prompt = {xmin,ymin,xmax,ymax};
dimRoI = inputdlg(prompt,'RoI in the cross-sectional direction', [1,60], {'200','200','1799','1799'});
rect = str2double(dimRoI);
end
im_rect = double(im(rect(2):rect(4),rect(1):rect(3)))/2^16;
end
function track_slicer(filename,centers,clim)
%TRACK_SLICER Slice trough data, centers and tracks.
% TRACK_SLICER(FILENAME,CENTERS,CLIM)
% FILENAME, name of the tiff stack with image data.
% CENTERS, detected fibre centers. If CENTERS contains track ids (as
% returned by fibre_tracking function) those will be shown as
% numbers.
% CLIM, optional color limits.
% Author: vand@dtu.dk 2018
if ~iscell(centers) % then tracks are given, and need to be converted
tracks = centers;
tracks = permute(tracks,[2,3,1]);
tracks(:,3,:) = repmat((1:size(tracks,1))',[1 1 size(tracks,3)]); % using slice number not z coordinate
centers = permute(mat2cell(tracks,size(tracks,1),size(tracks,2),ones(size(tracks,3),1)),[3 2 1]);
end
Z = length(imfinfo(filename));
z = round(0.5*(Z+1));
%update_drawing
I = imread(filename,'Index',z);
dim = [size(I),Z];
if nargin<3 || isempty(clim)
clim = [min(I(:)),max(I(:))];
end
show_numbers = true;
fig = figure('Units','Normalized','Position',[0.1 0.3 0.5 0.6],...
'KeyPressFcn',@key_press);
imagesc(I,clim)
hold on, colormap gray, axis image ij
plot(centers{z}(:,1),centers{z}(:,2),'r.')
title(['slice ',num2str(z),'/',num2str(Z)])
drawnow
% to capture zoom
LIMITS = [0.5,dim(2)+0.5,0.5,dim(1)+0.5];
zoom_handle = zoom(fig);
pan_handle = pan(fig);
set(zoom_handle,'ActionPostCallback',@adjust_limits,...
'ActionPreCallback',@force_keep_key_press);
set(pan_handle,'ActionPostCallback',@adjust_limits,...
'ActionPreCallback',@force_keep_key_press);
% if nargout>1
% uiwait(fig) % waits with assigning output until a figure is closed
% end
%%%%%%%%%% CALLBACK FUNCTIONS %%%%%%%%%%
function key_press(~,object)
% keyboard commands
key = object.Key;
switch 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 'pageup'
z = min(z+50,dim(3));
case 'pagedown'
z = max(z-50,1);
case 'n'
show_numbers = ~show_numbers;
end
update_drawing
end
%%%%%%%%%% HELPING FUNCTIONS %%%%%%%%%%
function update_drawing
I = imread(filename,'Index',z);
cla, imagesc(I,clim), hold on
if ~show_numbers || size(centers{z},2)<3 || (LIMITS(2)-LIMITS(1))*(LIMITS(4)-LIMITS(3))>300^2
plot(centers{z}(:,1),centers{z}(:,2),'r.')
else
present = (centers{z}(:,1)<LIMITS(2))&(centers{z}(:,1)>LIMITS(1))&...
(centers{z}(:,2)<LIMITS(4))&(centers{z}(:,2)>LIMITS(3));
text(centers{z}(present,1),centers{z}(present,2),num2cell(centers{z}(present,3)'),...
'HorizontalAlignment','center','VerticalAlignment','middle',...
'Color','r')
end
title(['slice ',num2str(z),'/',num2str(Z)])
drawnow
end
function adjust_limits(~,~)
% response to zooming and panning
LIMITS([1,2]) = get(gca,'XLim');
LIMITS([3,4]) = get(gca,'YLim');
if show_numbers
update_drawing
end
end
function force_keep_key_press(~,~)
% a hack to maintain my key_press while in zoom and pan mode
% http://undocumentedmatlab.com/blog/enabling-user-callbacks-during-zoom-pan
hManager = uigetmodemanager(fig);
[hManager.WindowListenerHandles.Enabled] = deal(false);
set(fig, 'KeyPressFcn', @key_press);
end
end
%ADDPATHS_INSEGTFIBRE Makes all the code of Insegt Fibre accesible
% Written by Monica Jane Emerson, March 2019, MIT License
%USER INPUT: The path of the textureGUI folder in your computer
x = inputdlg('Confirm the default path of the texture_GUI folder: ',...
'User input',[1,30],{'./texture_gui'});
path_textureGUI = x{1};
%Add the functions and Insegt software path
addpath(genpath(path_textureGUI))
addpath('./scripts')
addpath('./functions')
\ No newline at end of file
%CREATE_DICT Updates the dictionary of intensities built by SETUP_DICT with
%a training data-set that can be loaded or created from scratch.
% Written by Monica Jane Emerson, March 2019, MIT License
%USER INPUT: Option to load a labelling that was saved
x = questdlg('Would you like to load a labelling or create a new one?','Annotation of training data','New one','Load','New one');
if strcmp(x,'Load')
[filename,pathname] = uigetfile('*.mat','Choose training image','../data/saved_data/trainingData/imtrain.mat');
load([pathname,filename]);
[filename,pathname] = uigetfile('*labels_indexed.png','Choose label image','../data/saved_data/trainingData/');
labelling = imread([pathname,filename]);
dictionary = build_dictionary(im_train,dictopt); %create the dictionary of intensities
image_texture_gui(im_train,dictionary,2,labelling)
else
%select a small RoI for training, containing 100 fibres approx
prompt = ['Choose a slice for training: [1,', num2str(size(V,3)),']:'];
x = inputdlg(prompt,...
'User input',[1,30],{num2str(round(size(V,3)/2))}); %select slice
path_trainSlice = [path_volumeFolder,contents_datafolder(str2double(x{1})+2).name];
[~,im_train] = selectFoV(path_trainSlice);
uisave('im_train','../data/saved_data/trainingData/imtrain.mat')
dictionary = build_dictionary(im_train,dictopt); %create the dictionary of intensities
image_texture_gui(im_train,dictionary,2) %learn the dictionary of probabilities by annotating in the GUI
end
dictionary = update_dictionary(dictionary,gui_dictprob); %update dictionary to include the learnt probalilities
%USER INPUT: Save dictionary to process the complete scan
uisave('dictionary','../data/saved_data/dictionaries/dictionary_v.mat')
%INDICATE_DATAFOLDER Allows the user to specify the location of their data,
%and displays a slice to ensure that the right data is contained at the
%specified location.
% Written by Monica Jane Emerson, March 2019, MIT License
%USER INPUT: Path to CT volume?
path_volumeFolder = uigetdir('../data/','Open the folder with your CT data');
if ~strcmp(path_volumeFolder(end),'/')
path_volumeFolder = [path_volumeFolder,'/'];
end
%add the path and check the contents of the data folder
addpath(genpath(path_volumeFolder))
contents_datafolder = dir(path_volumeFolder);
%display cross-section, first slice
path_firstSlice = [path_volumeFolder,contents_datafolder(3).name];
im = imread(path_firstSlice);
figure, imagesc(im), axis image, colormap gray, title('First slice of the loaded data-set')
\ No newline at end of file
function [ depth ] = input_depthRoI(varargin)
%INPUT_DEPTHROI asks the user to specify the depth range to analyse
% Details...
%by Monica Jane Emerson, May 2019, MIT License
if nargin == 1
num_slices = varargin{1};
else
num_slices = [];
end
zmin = ['Lower limit in the depth direction [1, ',num2str(num_slices),']:'];
zmax = ['Upper limit in the depth direction [1, ',num2str(num_slices),']:'];
zstep = ['Step (1 to take every slice):'];
prompt = {zmin,zmax,zstep};
depth = inputdlg(prompt,'RoI in the depth direction', [1,50], {'1',num2str(num_slices),'10'});
end
%LOAD_MULTIPLE_TRACKSETS is a user-friendly script to load tracks measured
%from multiple data-sets, e.g. acquired at progressive load steps.
% by Monica Jane Emerson, March 2019, MIT License
%% USER INPUT: Load several sets of tracks corresponding to different steps of a time-lapse data-set
x = inputdlg('How many sets of tracks would you like to load [>=2]?','Nr. analysed matching volumes',[1,60],{num2str(2)});
for l = 1:str2double(x{1})
[filename_tracks,pathname_tracks] = uigetfile('*.mat',['Choose tracks for data-set ',num2str(l)],'saved_data/fibreTracks');
aux = load([pathname_tracks,filename_tracks]);
tracks{l} = aux.fibres;
end
num_datasets = length(tracks);
num_slices = size(tracks{1}.x,2);
%SETUP_DICT Set-ups the default parameters for the dictionary, andasks for the
%data-dependent parameter (patch size)
% Written by Monica Jane Emerson, March 2019, MIT License
%USER INPUT: select patch size, the closest odd integer to factor*diam/pixel_size
x = inputdlg('Choose patch size: ','User input',[1,20],{num2str(11)}); %ask for patch size
dictopt.patch_size = str2double(x{1}); %save patch size
%Set the default parameters and initialise dictionary
dictopt.method = 'euclidean';
dictopt.branching_factor = 3; %>=3
dictopt.number_layers = 5; %>=5. The higher the more accurate, but also slower and more computationally heavy
dictopt.number_training_patches = 15000; %at least 10*num_dictAtoms(branching_factor,number_layers)
dictopt.normalization = false; %set to true if the global intensity of the slices varies along the depth of the volume
clear
close all
addpath functions
% VERSION 1: WITHOUT CORRECTION
% for batch processing dictionary is reused, and manual labeling might need
% correction, so easiest is to build dictionary outside the gui to be able
% to use it later
im = double(imread('../data/slice1.png'))/255;
dictopt.method = 'euclidean';
dictopt.patch_size = 9;
dictopt.branching_factor = 5;
dictopt.number_layers = 4;
dictopt.number_training_patches = 30000;
dictopt.normalization = false;
dictionary = build_dictionary(im,dictopt);
% IMPORTANT:
% once inside gui export dict_labels to workspace (E)
image_texture_gui(im,dictionary,2)
%%
dictionary = update_dictionary(dictionary,gui_dictprob);
%%
figure, imagesc(gui_S), axis image, title('results from gui')
figure
for i=1:5
I = imread(['../data/slice',num2str(i),'.png']);
S = process_image(I,dictionary);
subplot(1,5,i)
imagesc(S), axis image, title(i), drawnow
end
%%
% VERSION 2: WITH (OR WITHOUT) CORRECTION
% TODO: freezeing (F) and correcting
/*=================================================================
* syntax: B = biadjacency_matrix(A,M,K) OR B = biadjacency_matrix(A,M)
*
* BIADJACENCY_MATRIX - build biadjacancy matrix from assignment image
*
* Input: - A: X-by-Y assignment image
* - M: patch size (length of edge)
* - K: number of dictionary patches, defaults to max(A(:))
*
* Output: - B: XY-by-MMK sparse biadjacency matrix
*
* Author: Vedrana Andersen Dahl, vand@dtu.dk, december 2015.
*=================================================================*/
#include <math.h>
#include <stdio.h>
#include <vector>
#include <algorithm>
#include "mex.h"
// struct containing i and j indices of a sparse matrix element
struct ij
{
int i,j;
ij(int i_, int j_) : i(i_), j(j_){}; // constructor
bool operator<(const ij second) const{ // leq operator needed for sorting
return (j<second.j) || ((j==second.j) && (i<second.i));
}
};
// The gateway mex routine
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
/* Check for proper number of input and output arguments */
if ((nlhs != 1) || (nrhs != 2))
mexErrMsgTxt("Usage: B = BIADJANCENCY_MATRIX(A,M) OR "
"B = BIADJANCENCY_MATRIX(A,M,K).\n");
/* Read input */
double *A = mxGetPr(prhs[0]); // assignment image
int X = mxGetM(prhs[0]); // image size X
int Y = mxGetN(prhs[0]); // image size Y
int M = (int)mxGetPr(prhs[1])[0]; // patch size
int K; // number dict patches
if (nrhs==3)
K = (int)mxGetPr(prhs[2])[0];
else{ // assumes number of dict patches is max(A)
K = 0;
for (int a=0; a<X*Y; a++)
if (A[a]>K)
K = A[a];
}
/* Compute some useful sizes */
int c = (M-1)/2; // width of boundary having no assignment
int n = X*Y; // number of image pixels
int m = M*M*K; // number of dict pixels
int s = (X-M+1)*(Y-M+1)*M*M; // number image-dict links (elements in B)
/* Finding elements of B as row-column indices */
std::vector<ij> bij;
bij.reserve(s);
int ic,i,j;
for (int y=0+c; y<Y-c; y++){ // visiting patches centered around pixels
for (int x=0+c; x<X-c; x++){
ic = x+y*X; // central patch pixel
for (int dy=-c; dy<=c; dy++){ // visiting pixels around central
for (int dx=-c; dx<=c; dx++){
i = (x+dx)+(y+dy)*X;
j = (c+dx)+(c+dy)*M+(A[ic]-1)*M*M;
bij.push_back(ij(i,j));
}
}
}
}
/* Sorting elements in bij columnwise */
std::sort (bij.begin(), bij.end());
/* Placeholder for output */
plhs[0] = mxCreateSparseLogicalMatrix(n,m,s); // output mxArray, sparse logical matrix B
if (plhs[0]==NULL)
mexErrMsgTxt("Could not allocate enough memory!\n");
/* Access fields of output mxArray via pointers */
mwIndex *ir = mxGetIr(plhs[0]); // row index (0 indexed)
mwIndex *jc = mxGetJc(plhs[0]); // cumulative number of elements per column
mxLogical *pr = mxGetLogicals(plhs[0]); // element values (will be all true)
/* Converting row-column indices into row-cumulative column */
int k = 0; // for visiting elements of bij
jc[0] = 0; // first element of cumulative sum is 0
for (int bc=0; bc<m; bc++){ // all columns of B
jc[bc+1] = jc[bc];
while (k<bij.size() && bij[k].j==bc){
jc[bc+1]++;
ir[k] = bij[k].i;
pr[k] = true;
k++;
}
}
}
File added
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment