diff --git a/code/InSegtFibre_2D.m b/code/InSegtFibre_2D.m
new file mode 100755
index 0000000000000000000000000000000000000000..bf5ee958f04809ae6f9060c1573004c6c7088ac5
--- /dev/null
+++ b/code/InSegtFibre_2D.m
@@ -0,0 +1,221 @@
+% 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')
diff --git a/code/InsegtFibre_3D.m b/code/InsegtFibre_3D.m
new file mode 100755
index 0000000000000000000000000000000000000000..33dc83721940e987d868bbceb14785f6dc4ebf7a
--- /dev/null
+++ b/code/InsegtFibre_3D.m
@@ -0,0 +1,254 @@
+% 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
diff --git a/code/InsegtFibre_orientations.m b/code/InsegtFibre_orientations.m
new file mode 100755
index 0000000000000000000000000000000000000000..86cbcc38e3bff275064376c3bd25ff779e37d166
--- /dev/null
+++ b/code/InsegtFibre_orientations.m
@@ -0,0 +1,138 @@
+% 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
+
diff --git a/code/InsegtFibre_validate3DTracks.m b/code/InsegtFibre_validate3DTracks.m
new file mode 100755
index 0000000000000000000000000000000000000000..29ae8d1a6fc20c3296fe748d3a632463903c58b0
--- /dev/null
+++ b/code/InsegtFibre_validate3DTracks.m
@@ -0,0 +1,68 @@
+% 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)
diff --git a/code/functions/calc_eltsdict.m b/code/functions/calc_eltsdict.m
new file mode 100755
index 0000000000000000000000000000000000000000..6edad49199094c0b1ae84ce3bc2d86506867368d
--- /dev/null
+++ b/code/functions/calc_eltsdict.m
@@ -0,0 +1,8 @@
+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
+
diff --git a/code/functions/calcu_psize.m b/code/functions/calcu_psize.m
new file mode 100755
index 0000000000000000000000000000000000000000..76ee312d02214630eb8e97e8273ffb72aa251b07
--- /dev/null
+++ b/code/functions/calcu_psize.m
@@ -0,0 +1,15 @@
+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
+
diff --git a/code/functions/imshow3D.m b/code/functions/imshow3D.m
new file mode 100755
index 0000000000000000000000000000000000000000..15f49c4f93cf5a8c5857581b427ffff6297fa39d
--- /dev/null
+++ b/code/functions/imshow3D.m
@@ -0,0 +1,291 @@
+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
diff --git a/code/functions/loadVolumeRoI.m b/code/functions/loadVolumeRoI.m
new file mode 100755
index 0000000000000000000000000000000000000000..42631336b15b4ad85b71c9a6d82c583cd093227b
--- /dev/null
+++ b/code/functions/loadVolumeRoI.m
@@ -0,0 +1,42 @@
+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
+
diff --git a/code/functions/num_dictAtoms.m b/code/functions/num_dictAtoms.m
new file mode 100755
index 0000000000000000000000000000000000000000..968a080312d92111256017531e6462ebbadb0d48
--- /dev/null
+++ b/code/functions/num_dictAtoms.m
@@ -0,0 +1,8 @@
+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
+
diff --git a/code/functions/selectFoV.m b/code/functions/selectFoV.m
new file mode 100755
index 0000000000000000000000000000000000000000..8b3af9a06b1e4c2079cb104bf76b09d029f0622f
--- /dev/null
+++ b/code/functions/selectFoV.m
@@ -0,0 +1,29 @@
+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
+
diff --git a/code/functions/track_slicer.m b/code/functions/track_slicer.m
new file mode 100755
index 0000000000000000000000000000000000000000..7500058f54d9c5e3093fd574641ab6c4fdeab5f5
--- /dev/null
+++ b/code/functions/track_slicer.m
@@ -0,0 +1,109 @@
+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
diff --git a/code/scripts/addpaths_inSegtFibre.m b/code/scripts/addpaths_inSegtFibre.m
new file mode 100755
index 0000000000000000000000000000000000000000..9e96bfd90a2a552977cb960a5d67971f5d83316f
--- /dev/null
+++ b/code/scripts/addpaths_inSegtFibre.m
@@ -0,0 +1,12 @@
+%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
diff --git a/code/scripts/create_dict.m b/code/scripts/create_dict.m
new file mode 100755
index 0000000000000000000000000000000000000000..5ef7e62937405a54169181e9d0f3675242a4daf7
--- /dev/null
+++ b/code/scripts/create_dict.m
@@ -0,0 +1,31 @@
+%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')
diff --git a/code/scripts/indicate_dataFolder.m b/code/scripts/indicate_dataFolder.m
new file mode 100755
index 0000000000000000000000000000000000000000..1899d0f51838d081c8b284dfcf92bcb903228c93
--- /dev/null
+++ b/code/scripts/indicate_dataFolder.m
@@ -0,0 +1,19 @@
+%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
diff --git a/code/scripts/input_depthRoI.m b/code/scripts/input_depthRoI.m
new file mode 100644
index 0000000000000000000000000000000000000000..a6c8c5d5ac6f521f4b301372569c7c756ba58732
--- /dev/null
+++ b/code/scripts/input_depthRoI.m
@@ -0,0 +1,19 @@
+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
+
diff --git a/code/scripts/load_multiple_trackSets.m b/code/scripts/load_multiple_trackSets.m
new file mode 100755
index 0000000000000000000000000000000000000000..0a25847e5002a7d7752b130fd765720e28b411fb
--- /dev/null
+++ b/code/scripts/load_multiple_trackSets.m
@@ -0,0 +1,15 @@
+%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);
+
diff --git a/code/scripts/setup_dict.m b/code/scripts/setup_dict.m
new file mode 100755
index 0000000000000000000000000000000000000000..b6e8f81b6736bb5b963faa5330c3ae395f500c28
--- /dev/null
+++ b/code/scripts/setup_dict.m
@@ -0,0 +1,15 @@
+%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
+
diff --git a/code/texture_gui/code/batch_processing_script.m b/code/texture_gui/code/batch_processing_script.m
new file mode 100755
index 0000000000000000000000000000000000000000..9d87ea3884a90d90352f9ae69dd3b186bd652b83
--- /dev/null
+++ b/code/texture_gui/code/batch_processing_script.m
@@ -0,0 +1,39 @@
+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
diff --git a/code/texture_gui/code/functions/biadjacency_matrix.cpp b/code/texture_gui/code/functions/biadjacency_matrix.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..efb16ea08cc771d7e59ff26f4b2c0434d870ab15
--- /dev/null
+++ b/code/texture_gui/code/functions/biadjacency_matrix.cpp
@@ -0,0 +1,103 @@
+/*=================================================================
+* 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++;
+ }
+ }
+}
+
diff --git a/code/texture_gui/code/functions/biadjacency_matrix.mexa64 b/code/texture_gui/code/functions/biadjacency_matrix.mexa64
new file mode 100755
index 0000000000000000000000000000000000000000..703f486d9d5d8524e9e27c48c259071a086f90a1
Binary files /dev/null and b/code/texture_gui/code/functions/biadjacency_matrix.mexa64 differ
diff --git a/code/texture_gui/code/functions/biadjacency_matrix.mexmaci64 b/code/texture_gui/code/functions/biadjacency_matrix.mexmaci64
new file mode 100755
index 0000000000000000000000000000000000000000..f40f6e00dc3e5bf9fbbabf158abccc084bfc6a93
Binary files /dev/null and b/code/texture_gui/code/functions/biadjacency_matrix.mexmaci64 differ
diff --git a/code/texture_gui/code/functions/biadjacency_matrix.mexw64 b/code/texture_gui/code/functions/biadjacency_matrix.mexw64
new file mode 100755
index 0000000000000000000000000000000000000000..6985266719edbacda108ffd1a788c7c78c26a019
Binary files /dev/null and b/code/texture_gui/code/functions/biadjacency_matrix.mexw64 differ
diff --git a/code/texture_gui/code/functions/build_dictionary.m b/code/texture_gui/code/functions/build_dictionary.m
new file mode 100755
index 0000000000000000000000000000000000000000..6750331e3cb3640e4c36a65d336d4b938b765326
--- /dev/null
+++ b/code/texture_gui/code/functions/build_dictionary.m
@@ -0,0 +1,23 @@
+function dictionary = build_dictionary(image,dictionary_options)
+
+image = normalize_image(image);
+dictionary.options = dictionary_options;
+
+switch dictionary.options.method%dictionary_options.method
+ case 'euclidean'
+ dictionary.tree = build_km_tree(image,...
+ dictionary_options.patch_size,...
+ dictionary_options.branching_factor,...
+ dictionary_options.number_training_patches,...
+ dictionary_options.number_layers,...
+ dictionary_options.normalization);
+ case 'nxcorr'
+ dictionary.tree = build_km_tree_xcorr(image,...
+ dictionary_options.patch_size,...
+ dictionary_options.branching_factor,...
+ dictionary_options.number_training_patches,...
+ dictionary_options.number_layers);
+ otherwise
+ error('Unknown dictionary method.')
+
+end
diff --git a/code/texture_gui/code/functions/build_km_tree.cpp b/code/texture_gui/code/functions/build_km_tree.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..a1e1116738f18af5569fc9b55bc14ee31b694ebb
--- /dev/null
+++ b/code/texture_gui/code/functions/build_km_tree.cpp
@@ -0,0 +1,421 @@
+/*=================================================================
+* syntax: T = build_km_tree(I, M, b, t, L, n); OR T = build_km_tree(I, M, b, t, L);
+*
+* build_km_tree - build km-tree matrix from image
+*
+* Input: - I: X-by-Y intensity image
+* - M: patch size (length of edge)
+* - L: number of dictionary layers. This parameter is limited
+* such that the average number of patches in a leafnode is
+* greater than five
+* - b: branching factor
+* - t: number of training patches
+* - n: normalization (true or false), defaults to false
+*
+* Output: - T: MMl-by-K matrix where l is the number of layers
+* in the image (1 for grayscale and 3 for RGB)
+* and K is the number of nodes in the tree.
+*
+* Author: Anders Dahl, abda@dtu.dk, december 2015.
+*=================================================================*/
+
+#include "mex.h"
+#include <stdio.h>
+#include <math.h>
+#include "matrix.h"
+#include <vector>
+#include <algorithm>
+
+#include <iostream>
+using namespace std;
+
+// struct for image
+struct im_st
+{
+ double *im_data; // pointer to image data
+ int rows, cols, layers, n_pix; // rows, cols and layers are the image dimensions and n_pix = rows*cols
+};
+
+// struct for the tree
+struct tree_st
+{
+ double *tree_data;
+ int n_dim, n_nodes, branch_fac, M, Mh;
+};
+
+// struct for image patches
+struct im_patch
+{
+ double *patch_data; // pointer to the data
+ int idx; // id used for sorting the patches according to the tree
+ bool operator<(const im_patch& rhs) const {return idx < rhs.idx;} // operator needed for sorting
+};
+
+
+// Function for sampling patches from the image into the patch arrays
+// inputs reference to the image struct, tree struct, patch struct and position of the sampling coordinate.
+// There is no check if the sampling is outside the image
+void sample_patch(im_st& im, tree_st& tree, im_patch& patch, int r_im, int c_im, bool normalize)
+{
+ int id_l, id_r, id_i; // iterators for looking up image data
+ int id_p = 0; // iterator for looking up patch data
+ double sum_sq = 0, pix_val; // variables for normalization
+
+ for ( int l = 0; l < im.layers; l++ ){ // image is sampled by three nested loops
+ id_l = im.n_pix*l;
+ for ( int i = c_im-tree.Mh; i <= c_im+tree.Mh; i++ ){
+ id_r = id_l + i*im.rows;
+ for ( int j = r_im-tree.Mh; j <= r_im+tree.Mh; j++ ){
+ id_i = id_r + j;
+ pix_val = *(im.im_data + id_i);
+ *(patch.patch_data + id_p) = pix_val;
+ sum_sq += pix_val*pix_val;
+ id_p++;
+ }
+ }
+ }
+
+ if ( normalize ){ // normalization to unit length
+ double inv_sq = 1;
+ if ( sum_sq > 0 ){
+ inv_sq = 1/sqrt(sum_sq); // divide by sum of squares
+ }
+ for ( int i = 0; i < tree.n_dim; i++ ){
+ *(patch.patch_data + i) = (*(patch.patch_data + i))*inv_sq;
+ }
+ }
+
+
+}
+
+
+// function for randomly permuted set of indices
+// n is the numbers to choose from, n_set is the number of random indices
+// that is returned. Returns a vector of indices
+vector<int> randperm( int n, int n_set ) {
+ if ( n_set > n ){ // check if n_set exceeds n
+ n_set = n;
+ }
+
+ vector<int> rid;
+ rid.reserve(n); // vector of indices
+ for ( int i = 0; i < n; i++ ){ // set all indices in order
+ rid.push_back(i);
+ }
+
+ int t, id; // place holders for id and temporary number
+ int r_max = RAND_MAX; // maximum random number
+ for ( int i = 0; i < n_set; i++ ){
+ // choose a random number between i and n-1 and swap place between i and id
+ if ( LONG_MAX > RAND_MAX && n-i-1>RAND_MAX ){ // not enough with a random integer up til RAND_MAX
+ id = ((rand()*(r_max+1)+rand()) % (n-i)) + i;
+ }
+ else{
+ id = (rand() % (n-i)) + i;
+ }
+ t = rid[id];
+ rid[id] = rid[i];
+ rid[i] = t;
+ }
+ rid.resize(n_set); // set size to n_set
+ return rid;
+}
+
+
+// copy values from a patch array into the tree array at node
+void set_values(tree_st& tree, im_patch& patch, int node){
+ int idx = tree.n_dim*node;
+ for ( int i = 0; i < tree.n_dim; i++ ){
+ *(tree.tree_data + idx) = *(patch.patch_data + i);
+ idx++;
+ }
+}
+
+// add values to vector of cluster center points
+void add_values(vector<double>& center_sum, im_patch& patch, int id, int n_dim){
+ int idx = n_dim*id;
+ for ( int i = 0; i < n_dim; i++ ){
+ center_sum[idx] += *(patch.patch_data + i);
+ idx++;
+ }
+}
+
+// estimates the squared Euclidian distance between an image patch and a tree node
+double get_dist(tree_st& tree, im_patch& patch, int node)
+{
+ double d = 0, tmp;
+ int id = tree.n_dim*node;
+
+ for ( int i = 0; i < tree.n_dim; i++ ){
+ tmp = *(tree.tree_data + id) - *(patch.patch_data + i);
+ d += tmp*tmp;
+ id++;
+
+ }
+
+ return d;
+}
+
+// k-means-function taking a reference to the vector of image patches and a
+// tree struct as input f and t gives the image patches that should be clustered.
+// node is the first node in the tree included in the clustering
+void k_means( vector<im_patch>& patches, tree_st& tree, int f, int t, int node )
+{
+ // vectors holding the sum of values in the cluster and a vector containing the change
+ vector<double> centSum(tree.branch_fac*tree.n_dim), centChange(tree.branch_fac);
+ vector<int> centCount(tree.branch_fac); // vector for counting the number of points in a cluster
+ double min_d, d, tmp;//, diff; // variables for clustering
+ // variables for cluster id and index of previous cluseter, which will be overwritten by new cluster id
+ int id = 0, id_in = patches[f].idx;
+
+ if ( t-f > tree.branch_fac ){ // check if there are enough point to carry out the clustering
+ // initialize the center positions
+ for ( int i = 0; i < tree.branch_fac; i++ ){
+ set_values(tree, patches[f+i], node+i);
+ }
+
+ // run clutering for 20 iterations - only little change happens after 10 iterations
+ for ( int n_run = 0; n_run < 30; n_run++){
+
+ for ( int i = f; i < t; i++ ){ // go through the patches from f to t
+ min_d = get_dist(tree, patches[i], node); // initially set min distance and id to the first
+ id = 0;
+ for ( int j = 1; j < tree.branch_fac; j++ ){ // run throgh the other points
+ d = get_dist(tree, patches[i], node + j); // get the distance
+ if ( d < min_d ){ // if the this cluster is closer set this as min distance
+ min_d = d;
+ id = j;
+ }
+ }
+ add_values(centSum, patches[i], id, tree.n_dim); // add the patch to the closest cluster
+ centCount[id]++; // count the extra patch
+ // update the idx to the child idx - note that all layers start with idx = 0
+ patches[i].idx = (id + id_in*tree.branch_fac);
+ }
+
+ // update the clusters in the tree and calculate the center change (not used for anything)
+ id = node*tree.n_dim;
+ int id_c = 0;
+
+ for ( int i = 0; i < tree.branch_fac; i++ ){ // go through all new clusters
+ if ( centCount[i] == 0 ){
+ centCount[i] = 1;
+ }
+ for ( int j = 0; j < tree.n_dim; j++ ){ // go through cluster pixels
+ tmp = centSum[id_c]/centCount[i];
+ //diff = (tmp - *(tree.tree_data + id)); // for calculating center change
+ //centChange[i] += diff*diff;
+ *(tree.tree_data + id) = tmp;
+ id_c++;
+ id++;
+ }
+ }
+
+ // set counter and sum to zero
+ fill(centSum.begin(), centSum.end(), 0);
+ fill(centCount.begin(), centCount.end(),0);
+ fill(centChange.begin(), centChange.end(), 0);
+ }
+ }
+}
+
+// runs through the patches vector to find the last element with id
+int get_to( vector<im_patch>& patches, int id )
+{
+ int to = 0;
+ for ( int i = 0; i < patches.size(); i++ ){
+ if ( patches[i].idx == id ){
+ to = i;
+ }
+ }
+ return to+1;
+}
+
+// Main function for building the km-tree. Takes the image and tree struct
+// and the number of training patches as argument
+void build_km_tree ( im_st& im, tree_st& tree, int n_train, bool normalize ) {
+ // allocate memory for the image patches
+ double* im_patch_data = new double[n_train*tree.M*tree.M*im.layers];
+
+ int rows_c = im.rows-tree.M+1, cols_c = im.cols-tree.M+1; // number of rows and cols within sampling area
+ int n_im_patches = rows_c*cols_c; // number of pixels in the image for sampling - inside boundary
+
+ // checks that the number of training patches is not exceeding the number of pixels in the sample area
+ if ( n_im_patches < n_train ){
+ n_train = n_im_patches;
+ }
+
+ vector<int> r_id = randperm(n_im_patches, n_train); // indices of random patches
+
+ vector<im_patch> patches; // vector of image patches
+ patches.resize(n_train); // allocate memory
+
+ int r, c, idx = 0; // variables used for sampling the image
+ // sample image patches
+ for (int i = 0; i < n_train; i++ )
+ {
+ c = r_id[i]/rows_c; // column is floored because of int
+ r = r_id[i]-c*rows_c; // row is rest after column
+ patches[i].idx = 0; // inital id is 0
+ patches[i].patch_data = im_patch_data + idx; // pointer to patch memory
+ sample_patch(im, tree, patches[i], r + tree.Mh, c + tree.Mh, normalize); // sampel in image with added boundary
+ idx += tree.n_dim; // step number of patch pixels forward
+ }
+
+ // k-means tree
+ int n_layer = (int)ceil(log((double)tree.n_nodes)/log((double)tree.branch_fac)); // number of layers in the tree
+ int n_in_layer; // number of nodes in layer
+ int t, f; // keeps track of patches that belong to a certain cluster
+ int node = 0; // node number that will be updated
+
+ // go through the layers in the tree
+ for (int i = 0; i < n_layer; i++ )
+ {
+ t = 0; // start at 0
+ n_in_layer = pow((double)tree.branch_fac,i); // number of nodes in current layer of the tree
+ sort(patches.begin(), patches.end()); // sort the patches according to their current id
+ for ( int j = 0; j < n_in_layer; j++ ) // go through the nodes in the layer and cluster that node
+ {
+ f = t; // patch j from
+ t = get_to(patches,j); // patch j to
+ // check that the node does not exceed the size of the tree
+ if ( node + tree.branch_fac <= tree.n_nodes ){
+ k_means( patches, tree, f, t, node );
+ }
+ else {
+ break;
+ }
+ node += tree.branch_fac; // next node
+ }
+ }
+
+ delete[] im_patch_data; // free up patch memory
+}
+
+
+// The gateway routine
+void mexFunction( int nlhs, mxArray *plhs[],
+ int nrhs, const mxArray *prhs[])
+{
+ // input image (I), patch size (M*M), number of nodes in the tree (n), branching
+ // factor (b), and number of training patches (n_train). Outputs the km-tree (tree)
+ double *I, *tree; // pointers to image and tree
+ int b, M, L, n, ndim, n_train; // variables
+ const int *dim; // image dimensinos
+ int dtree[2]; // tree dimensions
+ bool normalize = false; // decide if vectors of image patches should be normalized to unit length
+
+ /* Check for proper number of arguments. */
+ /* NOTE: You do not need an else statement when using
+ mexErrMsgTxt within an if statement. It will never
+ get to the else statement if mexErrMsgTxt is executed.
+ (mexErrMsgTxt breaks you out of the MEX-file.)
+ */
+ if(nrhs < 5 || nrhs > 6)
+ mexErrMsgTxt("Five or six inputs required.");
+ if(nlhs != 1)
+ mexErrMsgTxt("One output required.");
+
+ if ( !mxIsDouble(prhs[0]) )
+ mexErrMsgTxt("Image should be double!");
+
+ // Create a pointer to the input matrix.
+ I = mxGetPr(prhs[0]);
+
+ // input passing
+ double *Md, *bd, *Ld, *n_train_d;
+ Md = mxGetPr(prhs[1]);
+ M = (int)Md[0];
+
+ bd = mxGetPr(prhs[2]);
+ b = (int)bd[0];
+
+ // check if number of clusters is smaller than branching factor
+ if ( n < b ){
+ n = b;
+ }
+
+ n_train_d = mxGetPr(prhs[3]);
+ n_train = (int)n_train_d[0];
+
+ // determine number of tree nodes
+ Ld = mxGetPr(prhs[4]);
+ L = (int)Ld[0]; // layers in tree
+ n = 0;
+ int n_tmp = 0;
+ int max_n = (double)n_train/1.0;
+// int max_n = (double)n_train/5.0;
+ for ( int i = 0; i < L; i++ ){
+ n_tmp += pow((double)b,(i+1));
+ if ( n_tmp > max_n ){
+ L = i+1;
+ break;
+ }
+ n = n_tmp;
+ }
+ printf("Number of nodes in resulting tree: %d in %d layers.\n", n, L);
+
+ if ( nrhs == 6 ){
+ bool *normalize_d;
+ normalize_d = (bool *)mxGetData(prhs[5]);
+ normalize = normalize_d[0];
+ }
+
+ // check input properties
+ if ( 1 - (M % 2) || M < 1)
+ mexErrMsgTxt("M must be odd and positive.");
+
+ if ( n < 1 )
+ mexErrMsgTxt("n must be positive.");
+
+ if ( b < 1 )
+ mexErrMsgTxt("b must be positive.");
+
+ // Get the dimensions of the matrix input.
+ ndim = mxGetNumberOfDimensions(prhs[0]);
+ if (ndim != 2 && ndim != 3)
+ mexErrMsgTxt("search_km_tree only works for 2-dimensional or 3-dimensional images.");
+
+ // image dimensions
+ dim = mxGetDimensions(prhs[0]);
+
+
+ // image struct
+ im_st Im;
+ Im.im_data = I;
+ Im.rows = dim[0];
+ Im.cols = dim[1];
+ if ( ndim == 3 )
+ {
+ Im.layers = dim[2];
+ }
+ else
+ {
+ Im.layers = 1;
+ }
+ Im.n_pix = Im.rows*Im.cols;
+
+ dtree[0] = Im.layers*M*M;
+ dtree[1] = n;
+
+ // Set the output pointer to the output matrix. Array initialized to zero.
+ plhs[0] = mxCreateNumericArray(2, dtree, mxDOUBLE_CLASS, mxREAL);
+
+ // Create a C pointer to a copy of the output matrix.
+ tree = mxGetPr(plhs[0]);
+ for (int i = 0; i < dtree[0]*dtree[1]; i++ )
+ *(tree + i) = -1;
+
+ // tree struct
+ tree_st Tree;
+ Tree.tree_data = tree;
+ Tree.n_dim = dtree[0];
+ Tree.n_nodes = dtree[1];
+ Tree.branch_fac = b;
+ Tree.M = M;
+ Tree.Mh = (int)(0.5*((double)M-1.0));
+
+ // build the km-tree using C++ subroutine
+ build_km_tree ( Im, Tree, n_train, normalize );
+
+}
diff --git a/code/texture_gui/code/functions/build_km_tree.mexa64 b/code/texture_gui/code/functions/build_km_tree.mexa64
new file mode 100755
index 0000000000000000000000000000000000000000..64f090ccc0f775fd55f0c576b8e29df61ad8ff29
Binary files /dev/null and b/code/texture_gui/code/functions/build_km_tree.mexa64 differ
diff --git a/code/texture_gui/code/functions/build_km_tree.mexmaci64 b/code/texture_gui/code/functions/build_km_tree.mexmaci64
new file mode 100755
index 0000000000000000000000000000000000000000..af834dcdddc581cf19474646a2383c938c2441ca
Binary files /dev/null and b/code/texture_gui/code/functions/build_km_tree.mexmaci64 differ
diff --git a/code/texture_gui/code/functions/build_km_tree.mexw64 b/code/texture_gui/code/functions/build_km_tree.mexw64
new file mode 100755
index 0000000000000000000000000000000000000000..fc2e3cefa7e292a1cc3fd9baf7d06414b15fa888
Binary files /dev/null and b/code/texture_gui/code/functions/build_km_tree.mexw64 differ
diff --git a/code/texture_gui/code/functions/build_km_tree_xcorr.cpp b/code/texture_gui/code/functions/build_km_tree_xcorr.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..0e15bbf7e5a66e4f0102ecc7e8c4692a933097f5
--- /dev/null
+++ b/code/texture_gui/code/functions/build_km_tree_xcorr.cpp
@@ -0,0 +1,427 @@
+/*=================================================================
+* syntax: T = build_km_tree_xcorr(I, M, L, b, t);
+*
+* build_km_tre:xcorr - build km-tree matrix from image based on cross correlation
+*
+* Input: - I: X-by-Y intensity image
+* - M: patch size (length of edge)
+* - L: number of dictionary layers. This parameter is limited
+* such that the average number of patches in a leafnode is
+* greater than five
+* - b: brancing factor
+* - t: number of training patches
+*
+* Output: - T: MMl-by-K matrix where l is the number of layers
+* in the image (1 for grayscale and 3 for RGB)
+* and K is the number of nodes in the tree.
+*
+* Author: Anders Dahl, abda@dtu.dk, december 2015.
+*=================================================================*/
+
+#include "mex.h"
+#include <stdio.h>
+#include <math.h>
+#include "matrix.h"
+#include <vector>
+#include <algorithm>
+
+#include <iostream>
+using namespace std;
+
+// struct for image
+struct im_st
+{
+ double *im_data; // pointer to image data
+ int rows, cols, layers, n_pix; // rows, cols and layers are the image dimensions and n_pix = rows*cols
+};
+
+// struct for the tree
+struct tree_st
+{
+ double *tree_data;
+ int n_dim, n_nodes, branch_fac, M, Mh;
+ double inv_ndim;
+};
+
+// struct for image patches
+struct im_patch
+{
+ double *patch_data; // pointer to the data
+ int idx; // id used for sorting the patches according to the tree
+ bool operator<(const im_patch& rhs) const {return idx < rhs.idx;} // operator needed for sorting
+};
+
+// nomalize patch vectors to zero mean and unit Euclidean length
+void norm_patch(double* patch_data, double& mu, int n_dim){
+
+ double sum_sq = 0;
+
+ for ( int i = 0; i < n_dim; i++ ){
+ *(patch_data + i) = *(patch_data + i) - mu;
+ sum_sq += (*(patch_data + i))*(*(patch_data + i));
+ }
+
+ double inv_sq = 1;
+ if ( sum_sq > 0 ){
+ inv_sq = 1.0/sqrt(sum_sq); // divide by sum of squares
+ }
+
+ for ( int i = 0; i < n_dim; i++ ){
+ *(patch_data + i) = (*(patch_data + i))*inv_sq;
+ }
+}
+
+
+// Function for sampling patches from the image into the patch arrays
+// inputs reference to the image struct, tree struct, patch struct and position of the sampling coordinate.
+// There is no check if the sampling is outside the image
+void sample_patch(im_st& im, tree_st& tree, im_patch& patch, int r_im, int c_im)
+{
+ int id_l, id_r, id_i; // iterators for looking up image data
+ int id_p = 0; // iterator for looking up patch data
+ double sum_sq = 0, sum = 0, pix_val, mu, tmp; // variables for normalization
+
+ for ( int l = 0; l < im.layers; l++ ){ // image is sampled by three nested loops
+ id_l = im.n_pix*l;
+ for ( int i = c_im-tree.Mh; i <= c_im+tree.Mh; i++ ){
+ id_r = id_l + i*im.rows;
+ for ( int j = r_im-tree.Mh; j <= r_im+tree.Mh; j++ ){
+ id_i = id_r + j;
+ pix_val = *(im.im_data + id_i);
+ *(patch.patch_data + id_p) = pix_val;
+ sum += pix_val;
+ id_p++;
+ }
+ }
+ }
+
+ mu = sum*tree.inv_ndim;
+
+ norm_patch(patch.patch_data, mu, tree.n_dim);
+}
+
+
+// function for randomly permuted set of indices
+// n is the numbers to choose from, n_set is the number of random indices
+// that is returned. Returns a vector of indices
+vector<int> randperm( int n, int n_set ) {
+ if ( n_set > n ){ // check if n_set exceeds n
+ n_set = n;
+ }
+
+ vector<int> rid;
+ rid.reserve(n); // vector of indices
+ for ( int i = 0; i < n; i++ ){ // set all indices in order
+ rid.push_back(i);
+ }
+
+ int t, id; // place holders for id and temporary number
+ int r_max = RAND_MAX; // maximum random number
+ for ( int i = 0; i < n_set; i++ ){
+ // choose a random number between i and n-1 and swap place between i and id
+ if ( LONG_MAX > RAND_MAX && n-i-1>RAND_MAX ){ // not enough with a random integer up til RAND_MAX
+ id = ((rand()*(r_max+1)+rand()) % (n-i)) + i;
+ }
+ else{
+ id = (rand() % (n-i)) + i;
+ }
+ t = rid[id];
+ rid[id] = rid[i];
+ rid[i] = t;
+ }
+ rid.resize(n_set); // set size to n_set
+ return rid;
+}
+
+
+// copy values from a patch array into the tree array at node
+void set_values(tree_st& tree, im_patch& patch, int node){
+ int idx = tree.n_dim*node;
+ for ( int i = 0; i < tree.n_dim; i++ ){
+ *(tree.tree_data + idx) = *(patch.patch_data + i);
+ idx++;
+ }
+}
+
+// add values to vector of cluster center points
+void add_values(vector<double>& center_sum, im_patch& patch, int id, int n_dim){
+ int idx = n_dim*id;
+ for ( int i = 0; i < n_dim; i++ ){
+ center_sum[idx] += *(patch.patch_data + i);
+ idx++;
+ }
+}
+
+// estimates the normalized cross correlation between an image patch and a tree node
+double get_corr(tree_st& tree, im_patch& patch, int node)
+{
+ double corr = 0, tmp_t, tmp_p;
+ int id = tree.n_dim*node;
+
+ for ( int i = 0; i < tree.n_dim; i++ ){
+ tmp_t = *(tree.tree_data + id);
+ tmp_p = *(patch.patch_data + i);
+ corr += tmp_t*tmp_p;
+ id++;
+ }
+
+ return corr;
+}
+
+// k-means-function taking a reference to the vector of image patches and a
+// tree struct as input f and t gives the image patches that should be clustered.
+// node is the first node in the tree included in the clustering
+void k_means( vector<im_patch>& patches, tree_st& tree, int f, int t, int node )
+{
+ // vectors holding the sum of values in the cluster and a vector containing the change
+ vector<double> centSum(tree.branch_fac*tree.n_dim), centChange(tree.branch_fac);
+ vector<int> centCount(tree.branch_fac); // vector for counting the number of points in a cluster
+ double max_corr, corr, tmp, mu, sum;//, diff; // variables for clustering
+ // variables for cluster id and index of previous cluseter, which will be overwritten by new cluster id
+ int id = 0, id_in = patches[f].idx;
+
+ if ( t-f > tree.branch_fac ){ // check if there are enough point to carry out the clustering
+ // initialize the center positions
+ for ( int i = 0; i < tree.branch_fac; i++ ){
+ set_values(tree, patches[f+i], node+i);
+ }
+
+ // run clutering for 10 iterations - only little change happens after 10 iterations
+ for ( int n_run = 0; n_run < 10; n_run++){
+
+ for ( int i = f; i < t; i++ ){ // go through the patches from f to t
+ max_corr = get_corr(tree, patches[i], node); // initially set min distance and id to the first
+ id = 0;
+ for ( int j = 1; j < tree.branch_fac; j++ ){ // run throgh the other points
+ corr = get_corr(tree, patches[i], node + j); // get the cross correlation
+ if ( corr > max_corr ){ // if the this cluster is closer set this as max correlation
+ max_corr = corr;
+ id = j;
+ }
+ }
+ add_values(centSum, patches[i], id, tree.n_dim); // add the patch to the closest cluster
+ centCount[id]++; // count the extra patch
+ // update the idx to the child idx - note that all layers start with idx = 0
+ patches[i].idx = (id + id_in*tree.branch_fac);
+ }
+
+ // update the clusters in the tree and calculate the center change (not used for anything)
+ id = node*tree.n_dim;
+ int id_c = 0;
+
+
+ for ( int i = 0; i < tree.branch_fac; i++ ){ // go through all new clusters
+ sum = 0;
+ if ( centCount[i] == 0 ){
+ centCount[i] = 1;
+ }
+ for ( int j = 0; j < tree.n_dim; j++ ){ // go through cluster pixels
+ tmp = centSum[id_c]/centCount[i];
+ //diff = (tmp - *(tree.tree_data + id)); // for calculating center change
+ //centChange[i] += diff*diff;
+ *(tree.tree_data + id) = tmp;
+ sum += tmp;
+ id_c++;
+ id++;
+ }
+ mu = sum*tree.inv_ndim;
+ norm_patch(tree.tree_data+id-tree.n_dim, mu, tree.n_dim);
+ }
+
+ // set counter and sum to zero
+ fill(centSum.begin(), centSum.end(), 0);
+ fill(centCount.begin(), centCount.end(),0);
+ fill(centChange.begin(), centChange.end(), 0);
+ }
+ }
+}
+
+// runs through the patches vector to find the last element with id
+int get_to( vector<im_patch>& patches, int id )
+{
+ int to = 0;
+ for ( int i = 0; i < patches.size(); i++ ){
+ if ( patches[i].idx == id ){
+ to = i;
+ }
+ }
+ return to+1;
+}
+
+// Main function for building the km-tree. Takes the image and tree struct
+// and the number of training patches as argument
+void build_km_tree ( im_st& im, tree_st& tree, int n_train ) {
+ // allocate memory for the image patches
+ double* im_patch_data = new double[n_train*tree.M*tree.M*im.layers];
+
+ int rows_c = im.rows-tree.M+1, cols_c = im.cols-tree.M+1; // number of rows and cols within sampling area
+ int n_im_patches = rows_c*cols_c; // number of pixels in the image for sampling - inside boundary
+
+ // checks that the number of training patches is not exceeding the number of pixels in the sample area
+ if ( n_im_patches < n_train ){
+ n_train = n_im_patches;
+ }
+
+ vector<int> r_id = randperm(n_im_patches, n_train); // indices of random patches
+
+ vector<im_patch> patches; // vector of image patches
+ patches.resize(n_train); // allocate memory
+
+ int r, c, idx = 0; // variables used for sampling the image
+ // sample image patches
+ for (int i = 0; i < n_train; i++ )
+ {
+ c = r_id[i]/rows_c; // column is floored because of int
+ r = r_id[i]-c*rows_c; // row is rest after column
+ patches[i].idx = 0; // inital id is 0
+ patches[i].patch_data = im_patch_data + idx; // pointer to patch memory
+ sample_patch(im, tree, patches[i], r + tree.Mh, c + tree.Mh); // sampel in image with added boundary
+ idx += tree.n_dim; // step number of patch pixels forward
+ }
+
+ // k-means tree
+ int n_layer = (int)ceil(log((double)tree.n_nodes)/log((double)tree.branch_fac)); // number of layers in the tree
+ int n_in_layer; // number of nodes in layer
+ int t, f; // keeps track of patches that belong to a certain cluster
+ int node = 0; // node number that will be updated
+
+ // go through the layers in the tree
+ for (int i = 0; i < n_layer; i++ )
+ {
+ t = 0; // start at 0
+ n_in_layer = pow((double)tree.branch_fac,i); // number of nodes in current layer of the tree
+ sort(patches.begin(), patches.end()); // sort the patches according to their current id
+ for ( int j = 0; j < n_in_layer; j++ ) // go through the nodes in the layer and cluster that node
+ {
+ f = t; // patch j from
+ t = get_to(patches,j); // patch j to
+ // check that the node does not exceed the size of the tree
+ if ( node + tree.branch_fac <= tree.n_nodes ){
+ k_means( patches, tree, f, t, node );
+ }
+ else {
+ break;
+ }
+ node += tree.branch_fac; // next node
+ }
+ }
+
+ delete[] im_patch_data; // free up patch memory
+}
+
+
+// The gateway routine
+void mexFunction( int nlhs, mxArray *plhs[],
+ int nrhs, const mxArray *prhs[])
+{
+ // input image (I), patch size (M*M), number of nodes in the tree (n), branching
+ // factor (b), and number of training patches (n_train). Outputs the km-tree (tree)
+ double *I, *tree; // pointers to image and tree
+ int b, M, L, n, ndim, n_train; // variables
+ const int *dim; // image dimensinos
+ int dtree[2]; // tree dimensions
+
+ /* Check for proper number of arguments. */
+ /* NOTE: You do not need an else statement when using
+ mexErrMsgTxt within an if statement. It will never
+ get to the else statement if mexErrMsgTxt is executed.
+ (mexErrMsgTxt breaks you out of the MEX-file.)
+ */
+ if(nrhs != 5)
+ mexErrMsgTxt("Five inputs required.");
+ if(nlhs != 1)
+ mexErrMsgTxt("One output required.");
+
+ // Create a pointer to the input matrix.
+ I = mxGetPr(prhs[0]);
+
+ // input passing
+ double *Md, *bd, *Ld, *n_train_d;
+ Md = mxGetPr(prhs[1]);
+ M = (int)Md[0];
+
+ bd = mxGetPr(prhs[2]);
+ b = (int)bd[0];
+
+ // check if number of clusters is smaller than branching factor
+ if ( n < b ){
+ n = b;
+ }
+
+ n_train_d = mxGetPr(prhs[3]);
+ n_train = (int)n_train_d[0];
+
+ // determine number of tree nodes
+ Ld = mxGetPr(prhs[4]);
+ L = (int)Ld[0]; // layers in tree
+ n = 0;
+ int n_tmp = 0;
+ int max_n = (double)n_train/5.0;
+ for ( int i = 0; i < L; i++ ){
+ n_tmp += pow((double)b,(i+1));
+ if ( n_tmp > max_n ){
+ L = i+1;
+ break;
+ }
+ n = n_tmp;
+ }
+ printf("Number of nodes in resulting tree: %d in %d layers.\n", n, L);
+
+ // check input properties
+ if ( 1 - (M % 2) || M < 1)
+ mexErrMsgTxt("M must be odd and positive.");
+
+ if ( n < 1 )
+ mexErrMsgTxt("n must be positive.");
+
+ if ( b < 1 )
+ mexErrMsgTxt("b must be positive.");
+
+ // Get the dimensions of the matrix input.
+ ndim = mxGetNumberOfDimensions(prhs[0]);
+ if (ndim != 2 && ndim != 3)
+ mexErrMsgTxt("search_km_tree only works for 2-dimensional or 3-dimensional images.");
+
+ // image dimensions
+ dim = mxGetDimensions(prhs[0]);
+
+
+ // image struct
+ im_st Im;
+ Im.im_data = I;
+ Im.rows = dim[0];
+ Im.cols = dim[1];
+ if ( ndim == 3 )
+ {
+ Im.layers = dim[2];
+ }
+ else
+ {
+ Im.layers = 1;
+ }
+ Im.n_pix = Im.rows*Im.cols;
+
+ dtree[0] = Im.layers*M*M;
+ dtree[1] = n;
+
+ // Set the output pointer to the output matrix. Array initialized to zero.
+ plhs[0] = mxCreateNumericArray(2, dtree, mxDOUBLE_CLASS, mxREAL);
+
+ // Create a C pointer to a copy of the output matrix.
+ tree = mxGetPr(plhs[0]);
+ for (int i = 0; i < dtree[0]*dtree[1]; i++ )
+ *(tree + i) = -1;
+
+ // tree struct
+ tree_st Tree;
+ Tree.tree_data = tree;
+ Tree.n_dim = dtree[0];
+ Tree.inv_ndim = 1/(double)Tree.n_dim;
+ Tree.n_nodes = dtree[1];
+ Tree.branch_fac = b;
+ Tree.M = M;
+ Tree.Mh = (int)(0.5*((double)M-1.0));
+ // build the km-tree using C++ subroutine
+ build_km_tree ( Im, Tree, n_train );
+
+}
diff --git a/code/texture_gui/code/functions/build_km_tree_xcorr.mexa64 b/code/texture_gui/code/functions/build_km_tree_xcorr.mexa64
new file mode 100755
index 0000000000000000000000000000000000000000..4f23c01adc36f93d1ea6e2b3a2df0915bac3247f
Binary files /dev/null and b/code/texture_gui/code/functions/build_km_tree_xcorr.mexa64 differ
diff --git a/code/texture_gui/code/functions/build_km_tree_xcorr.mexmaci64 b/code/texture_gui/code/functions/build_km_tree_xcorr.mexmaci64
new file mode 100755
index 0000000000000000000000000000000000000000..592f56e947462978fc07a59c067768e67277fcc6
Binary files /dev/null and b/code/texture_gui/code/functions/build_km_tree_xcorr.mexmaci64 differ
diff --git a/code/texture_gui/code/functions/build_km_tree_xcorr.mexw64 b/code/texture_gui/code/functions/build_km_tree_xcorr.mexw64
new file mode 100755
index 0000000000000000000000000000000000000000..a0b030b73702276a175b4c59a684ad0e0d2e5ff6
Binary files /dev/null and b/code/texture_gui/code/functions/build_km_tree_xcorr.mexw64 differ
diff --git a/code/texture_gui/code/functions/compile_mex_functions.m b/code/texture_gui/code/functions/compile_mex_functions.m
new file mode 100755
index 0000000000000000000000000000000000000000..bab817b10e16d154f102e50dbd6c3311f70f473c
--- /dev/null
+++ b/code/texture_gui/code/functions/compile_mex_functions.m
@@ -0,0 +1,8 @@
+% compile mex files
+mex -largeArrayDims biadjacency_matrix.cpp
+mex build_km_tree.cpp % based on Euclidean distance
+mex search_km_tree.cpp % based on Euclidean distance
+mex build_km_tree_xcorr.cpp % based on normalized cross correlation
+mex search_km_tree_xcorr.cpp % based on normalized cross correlation
+mex probability_search_km_tree.cpp % based on normalized cross correlation
+
diff --git a/code/texture_gui/code/functions/compute_mappings.m b/code/texture_gui/code/functions/compute_mappings.m
new file mode 100755
index 0000000000000000000000000000000000000000..9abe4f73c4c25df190e78ea63377b3ffd8cf7335
--- /dev/null
+++ b/code/texture_gui/code/functions/compute_mappings.m
@@ -0,0 +1,21 @@
+function [mappings,A] = compute_mappings(image,dictionary)
+
+switch dictionary.options.method
+ case 'euclidean'
+ A = search_km_tree(image,...
+ dictionary.tree,...
+ dictionary.options.branching_factor,...
+ dictionary.options.normalization);
+ case 'nxcorr'
+ A = search_km_tree_xcorr(image,...
+ dictionary.tree,...
+ dictionary.options.branching_factor);
+ otherwise
+ error('Unknown dictionary method.')
+end
+B = biadjacency_matrix(A,dictionary.options.patch_size);
+
+[rc,nm] = size(B);
+mappings.T1 = sparse(1:nm,1:nm,1./(sum(B,1)+eps),nm,nm)*B';
+mappings.T2 = sparse(1:rc,1:rc,1./(sum(B,2)+eps),rc,rc)*B;
+
diff --git a/code/texture_gui/code/functions/image_texture_gui.m b/code/texture_gui/code/functions/image_texture_gui.m
new file mode 100755
index 0000000000000000000000000000000000000000..a3145772eb60534de8814ff7da8dc12b3465869c
--- /dev/null
+++ b/code/texture_gui/code/functions/image_texture_gui.m
@@ -0,0 +1,805 @@
+function image_texture_gui(image,dict,nr_labels,labeling)
+%IMAGE_TEXTURE_GUI Interactive texture segmentation of the image
+% IMAGE_TEXTURE_GUI(IMAGE,DICTOP,NR_LABELS,LABELING)
+% Called without imput arguments starts in a demo mode. Input:
+% IMAGE, may be rgb or grayscale, 0-to-1 or uint8
+% DICT, may be either:
+% - dictopt struct containing dictinary parameters,
+% - dictionary struct containing the dictionary,
+% - mappings struct containing the transformation matrices,
+% - not given or empty, so default options are loaded
+% NR_LABELS, number of labels, defaults to 2.
+% LABELING, optional initial labeling, given as either
+% - binary matrix of the size nr_pixels-by-nr_labels,
+% - labeling image of the same dimensions as IMAGE
+%
+% Keyboard controls (TODO update this for newest version):
+% L, shift+L and numerical values change label
+% T, shift+T, uparrow and downarros change pen thickness
+% W and shift+W change show option
+% M and shift+M change method
+% O and shift+O change overwrite option
+% R and shift+R change regularize option
+% U and shift+U change live update option
+% C and shift+C change colormap
+% A and shift+A change opacity (alpha)
+% S saves a project
+% E export
+% F freeze
+%
+% Author: vand@dtu.dk, 2015
+% char(100+[18 -3 10 0 -36 0 16 17 -54 0 7])
+%
+% TODO:
+% - when LABELING given, update segmentation upon startup - DONE, TESTING
+% - when zooming, make sure mouse overlay circle not displayed (frozen)
+% - option for changing opacity - DONE, TESTING
+
+%%%%%%%%%% DEFAULTS %%%%%%%%%%
+if nargin<1 % default example image
+ image = imread('bag.png');
+ dict.method = 'euclidean';
+ dict.patch_size = 15;
+ dict.branching_factor = 2;
+ dict.number_layers = 5;
+ dict.number_training_patches = 1000;
+ dict.normalization = false;
+ nr_labels = 2;
+else
+ if nargin<2 || isempty(dict)% default dictionary
+ dict.method = 'euclidean';
+ dict.patch_size = 3;
+ dict.branching_factor = 2;
+ dict.number_layers = 5;
+ dict.number_training_patches = 1000;
+ dict.normalization = false;
+ end
+ if nargin<3 || isempty(nr_labels) % defalult number of labels
+ nr_labels = 2;
+ else
+ nr_labels = double(nr_labels); % for colormaps to work properly
+ end
+end
+
+[r,c,~] = size(image);
+if nargin<4 || isempty(labeling) % default, unlabeled initial labeling
+ LABELING = zeros(r*c,nr_labels);
+else
+ parse_labeling_input(labeling)
+end
+
+%%%%%%%%%% PROBABILITY COMPUTING METHODS %%%%%%%%%%
+% to add a new method, add a function to the method_options list
+method_options = {...
+ @(labeling)distributed(labeling),...
+ @(labeling)two_max(labeling),...
+ @(labeling)two_max_over(labeling),...
+ @(labeling)two_cont(labeling),...
+ @(labeling)two_cont_over(labeling),...
+ };
+
+%%%%%%%%%% SETTINGS %%%%%%%%%%
+
+% method options
+METHOD_INDEX = 1; % initial method is the first on the list
+METHOD = method_options{METHOD_INDEX};
+METHOD_NAME = method_name_str;
+
+% brush label options
+LABEL = 1; % initial label is 1
+
+% brush thickness options
+thickness_options = [1 2 3 4 5 10 20 30 40 50 100 -1]; % last option (value -1) is 'fill'
+thickness_options_string = num2cell(thickness_options);
+thickness_options_string{end} = 'fill';
+THICKNESS_INDEX = 5; % initial pencil thickness is the fifth option
+RADIUS = thickness2radius(thickness_options(THICKNESS_INDEX)); % pencil radius
+
+% results visualization options
+show_options(1:2) = {'segmentation','overlay'};
+show_options((1:nr_labels)+2) = num2cell([repmat('probability ',[nr_labels,1]),num2str((1:nr_labels)')],2);
+SHOW_INDEX = 1;
+
+% overwrite option
+overwrite_options = {'no','yes'};
+OVERWRITE = false; % initialy no overwrite
+
+% regularization (smoothness) options
+regularize_options = [0 1 2 3 4 5 10 15 20 25 30 50];
+REGULARIZE_INDEX = 1; % initialy no regularization
+
+% visualization colormap options (for labelings and results)
+colormap_options = {@(x)[0.5,0.5,0.5;cool(x)], @(x)[0.5,0.5,0.5;spring(x)],...
+ @(x)[0.5,0.5,0.5;parula(x)]}; % gray color for unlabeled
+COLORMAP_INDEX = 1; % current colormap
+COLORS = colormap_options{COLORMAP_INDEX}(nr_labels); % visualization colors for label overlay
+
+% visualization opacity options (for labelings and results)
+color_weight_options = 0.1:0.1:0.9;
+COLOR_WEIGHT_INDEX = 2;
+COLOR_WEIGHT = color_weight_options(COLOR_WEIGHT_INDEX);
+
+% live update option
+live_update_options = {'off','on'};
+LIVE_UPDATE = true; % initially on
+
+% other settings
+nr_circle_pts = 16; % number of points defining a circular pencil
+
+%%%%%%%%%% INITIALIZATION AND LAYOUT %%%%%%%%%%
+
+[image,image_rgb,image_gray] = normalize_image(image); % impose 0-to-1 rgb
+LABELING_OVERLAY = image_rgb; % image overlaid labeling
+SEGMENTATION_OVERLAY = 0.5+zeros(r,c,3); % segmentation, optionally overlay
+
+fmar = [0.2 0.2]; % discance from screen edge to figure (x and y)
+amar = [0.02 0.02]; % margin around axes, relative to figure
+my = 0.85:0.04:0.95; % menu items y position
+mw = 0.15; % menu items width
+mx = 0.05:0.15:0.8;
+mh = 0.03; % menu items height
+cw = (0.25-0.15)/(nr_labels+1); % colorcube width
+cx = 0.15:cw:0.25; % colorcubes x position
+pointer_char = 'X';
+
+fig = figure('Units','Normalized','Position',[fmar,1-2*fmar],...
+ 'Pointer','watch','KeyPressFcn',@key_press,'InvertHardCopy', 'off',...
+ 'Name','Texture segmentation GUI');
+
+labeling_axes = axes('Units','Normalized','Position',[0,0,0.5,0.85]+[amar,-2*amar]);
+imagesc(LABELING_OVERLAY), axis image off, hold on
+segmentation_axes = axes('Units','Normalized','Position',[0.5,0,0.5,0.85]+[amar,-2*amar]);
+imagesc(SEGMENTATION_OVERLAY,[0,nr_labels]), axis image off, hold on
+
+clean_toolbar % also defines linkaxes tool
+
+uicontrol('String','Label [L] : ','Style','text','HorizontalAlignment','right',...
+ 'Units','Normalized','Position',[mx(1),my(3),mw,mh]);
+labels_text = uicontrol('String',num2str(LABEL),...
+ 'BackgroundColor',COLORS(LABEL+1,:),...
+ 'Style','text','HorizontalAlignment','left',...
+ 'Units','Normalized','Position',[mx(2),my(3),0.03,mh]);
+label_pointer = cell(nr_labels+1,1);
+for k = 1:nr_labels+1
+ label_pointer{k} = uicontrol('String',' ','Style','text',...
+ 'HorizontalAlignment','center','BackgroundColor',COLORS(k,:),...
+ 'Units','Normalized','Position',[cx(k),my(2),cw,mh],...
+ 'Enable','Inactive','ButtonDownFcn',@label_click,'UserData',k-1);
+end
+set(label_pointer{LABEL+1},'String',pointer_char);
+
+uicontrol('String','Thickness [T] : ','Style','text','HorizontalAlignment','right',...
+ 'Units','Normalized','Position',[mx(1),my(1),mw,mh]);
+thickness_text = uicontrol('String',thickness_options_string(THICKNESS_INDEX),...
+ 'Style','text','HorizontalAlignment','left',...
+ 'Units','Normalized','Position',[mx(2),my(1),mw,mh]);
+
+uicontrol('String','Show [W] : ','Style','text','HorizontalAlignment','right',...
+ 'Units','Normalized','Position',[mx(5),my(2),mw,mh]);
+show_text = uicontrol('String',show_options{SHOW_INDEX},...
+ 'Style','text','HorizontalAlignment','left',...
+ 'Units','Normalized','Position',[mx(6),my(2),mw,mh]);
+
+uicontrol('String','Live update [U] : ','Style','text','HorizontalAlignment','right',...
+ 'Units','Normalized','Position',[mx(5),my(3),mw,mh]);
+update_text = uicontrol('String',live_update_options(LIVE_UPDATE+1),...
+ 'Style','text','HorizontalAlignment','left',...
+ 'Units','Normalized','Position',[mx(6),my(3),mw,mh]);
+
+uicontrol('String','Method [M] : ','Style','text','HorizontalAlignment','right',...
+ 'Units','Normalized','Position',[mx(3),my(3),mw,mh]);
+method_text = uicontrol('String',METHOD_NAME,...
+ 'Style','text','HorizontalAlignment','left',...
+ 'Units','Normalized','Position',[mx(4),my(3),mw,mh]);
+
+uicontrol('String','Overwrite [O] : ','Style','text','HorizontalAlignment','right',...
+ 'Units','Normalized','Position',[mx(3),my(2),mw,mh]);
+overlay_text = uicontrol('String',overwrite_options(OVERWRITE+1),...
+ 'Style','text','HorizontalAlignment','left',...
+ 'Units','Normalized','Position',[mx(4),my(2),mw,mh]);
+
+uicontrol('String','Regularize [R] : ','Style','text','HorizontalAlignment','right',...
+ 'Units','Normalized','Position',[mx(3),my(1),mw,mh]);
+regularize_text = uicontrol('String',regularize_options(REGULARIZE_INDEX),...
+ 'Style','text','HorizontalAlignment','left',...
+ 'Units','Normalized','Position',[mx(4),my(1),mw,mh]);
+
+drawnow % pointer shows busy system
+
+LIMITS = [1,c-0.5,1,r+0.5]; % to capture zoom
+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);
+
+%%%%%%%%%% TEXTURE REPRESENTATION %%%%%%%%%%
+% initiolization: building a texture representation of the image
+T1 = 0;
+T2 = 0;
+parsing_dictionary_input; % T1 and T2 are assigned here
+
+if nargin>3
+ PROBABILITY = METHOD(LABELING);
+else
+ PROBABILITY = 1/nr_labels*ones(size(LABELING));
+ %labeling_overwrite % relevant only if we allow loading all settings
+ %regularize % relevant only if we allow loading all settings
+end
+compute_overlays
+
+% ready to draw
+set(fig,'Pointer','arrow','WindowButtonDownFcn',@start_draw,...
+ 'WindowButtonMotionFcn',@pointer_motion);
+XO = []; % current drawing point
+
+show_overlays(get_pixel);
+uiwait % waits with assigning output until a figure is closed
+
+%%%%%%%%%% CALLBACK FUNCTIONS %%%%%%%%%%
+ function key_press(~,object)
+ % keyboard commands
+ key = object.Key;
+ numkey = str2double(key);
+ if ~isempty(numkey) && numkey<=nr_labels;
+ change_label(numkey)
+ else
+ switch key
+ case 'l'
+ change_label(move_once(LABEL+1,nr_labels+1,...
+ any(strcmp(object.Modifier,'shift')))-1)
+ case 'uparrow'
+ THICKNESS_INDEX = move_once(THICKNESS_INDEX,...
+ numel(thickness_options),false);
+ RADIUS = thickness2radius(thickness_options(THICKNESS_INDEX));
+ set(thickness_text,'String',...
+ thickness_options_string(THICKNESS_INDEX));
+ show_overlays(get_pixel);
+ case 'downarrow'
+ THICKNESS_INDEX = move_once(THICKNESS_INDEX,...
+ numel(thickness_options),true);
+ RADIUS = thickness2radius(thickness_options(THICKNESS_INDEX));
+ set(thickness_text,'String',...
+ thickness_options_string(THICKNESS_INDEX));
+ show_overlays(get_pixel);
+ case 't'
+ THICKNESS_INDEX = move_once(THICKNESS_INDEX,...
+ numel(thickness_options),any(strcmp(object.Modifier,'shift')));
+ RADIUS = thickness2radius(thickness_options(THICKNESS_INDEX));
+ set(thickness_text,'String',...
+ thickness_options_string(THICKNESS_INDEX));
+ show_overlays(get_pixel);
+ case 'w'
+ SHOW_INDEX = move_once(SHOW_INDEX,numel(show_options),...
+ any(strcmp(object.Modifier,'shift')));
+ set(show_text,'String',show_options{SHOW_INDEX})
+ compute_overlays
+ show_overlays(get_pixel);
+ case 'c'
+ COLORMAP_INDEX = move_once(COLORMAP_INDEX,...
+ length(colormap_options),...
+ any(strcmp(object.Modifier,'shift')));
+ COLORS = colormap_options{COLORMAP_INDEX}(nr_labels);
+ set(labels_text,'BackgroundColor',COLORS(LABEL+1,:));
+ for kk = 1:nr_labels+1
+ set(label_pointer{kk},'BackgroundColor',COLORS(kk,:))
+ end
+ compute_overlays
+ show_overlays(get_pixel);
+ case 'a'
+ COLOR_WEIGHT_INDEX = move_once(COLOR_WEIGHT_INDEX,...
+ length(color_weight_options),...
+ any(strcmp(object.Modifier,'shift')));
+ COLOR_WEIGHT = color_weight_options(COLOR_WEIGHT_INDEX);
+ compute_overlays
+ show_overlays(get_pixel);
+ case 'm'
+ METHOD_INDEX = move_once(METHOD_INDEX,...
+ length(method_options),...
+ any(strcmp(object.Modifier,'shift')));
+ METHOD = method_options{METHOD_INDEX};
+ METHOD_NAME = method_name_str;
+ set(method_text,'String',METHOD_NAME)
+ PROBABILITY = METHOD(LABELING);
+ labeling_overwrite
+ regularize
+ compute_overlays
+ show_overlays(get_pixel);
+ case 'o'
+ OVERWRITE = ~OVERWRITE;
+ set(overlay_text,'String',overwrite_options{OVERWRITE+1})
+ PROBABILITY = METHOD(LABELING);
+ labeling_overwrite
+ regularize
+ compute_overlays
+ show_overlays(get_pixel);
+ case 'u'
+ LIVE_UPDATE = ~LIVE_UPDATE;
+ set(update_text,'String',live_update_options{LIVE_UPDATE+1})
+ if LIVE_UPDATE
+ PROBABILITY = METHOD(LABELING);
+ labeling_overwrite
+ regularize
+ compute_overlays
+ show_overlays(get_pixel);
+ end
+ case 'r'
+ REGULARIZE_INDEX = move_once(REGULARIZE_INDEX,...
+ numel(regularize_options),...
+ any(strcmp(object.Modifier,'shift')));
+ set(regularize_text,'String',...
+ regularize_options(REGULARIZE_INDEX))
+ PROBABILITY = METHOD(LABELING);
+ labeling_overwrite
+ regularize
+ compute_overlays
+ show_overlays(get_pixel);
+ case 's'
+ save_displayed
+ case 'e'
+ export_LS
+ case 'f'
+ freeze_LS
+ case 'q'
+ close(fig)
+ end
+ end
+ end
+
+ function label_click(source,~)
+ change_label(source.UserData)
+ end
+
+ function change_label(new_label)
+ set(label_pointer{LABEL+1},'String',' ');
+ LABEL = new_label;
+ set(label_pointer{LABEL+1},'String',pointer_char);
+ set(labels_text,'String',num2str(LABEL),...
+ 'BackgroundColor',COLORS(LABEL+1,:));
+ end
+
+ function adjust_limits(~,~)
+ % response to zooming and panning
+ LIMITS([1,2]) = get(labeling_axes,'XLim');
+ LIMITS([3,4]) = get(labeling_axes,'YLim');
+ 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
+
+ function start_draw(~,~)
+ % click in the image
+ x = get_pixel;
+ if is_inside(x)
+ if RADIUS>0 % thickness>0
+ XO = x;
+ M = disc(XO,RADIUS,nr_circle_pts,[r,c]);
+ set(fig,'WindowButtonMotionFcn',@drag_and_draw,...
+ 'WindowButtonUpFcn',@end_draw)
+ else % fill
+ M = fill(x);
+ end
+ update(M);
+ end
+ end
+
+ function drag_and_draw(~,~)
+ % drag after clicking in the image
+ x = get_pixel;
+ M = stadium(XO,x,RADIUS,nr_circle_pts,[r,c]);
+ update(M);
+ XO = x;
+ end
+
+ function end_draw(~,~)
+ % release click after clicking in the image
+ M = stadium(XO,get_pixel,RADIUS,nr_circle_pts,[r,c]);
+ update(M);
+ XO = [];
+ set(fig,'WindowButtonMotionFcn',@pointer_motion,...
+ 'WindowButtonUpFcn','')
+ end
+
+ function pointer_motion(~,~)
+ % move around without clicking
+ if strcmp(zoom_handle.Enable,'off') && ...
+ strcmp(pan_handle.Enable,'off') % not zooming or panning
+ x = get_pixel;
+ if is_inside(x)
+ set(fig,'Pointer','crosshair')
+ else
+ set(fig,'Pointer','arrow')
+ end
+ show_overlays(x);
+ end
+ end
+
+%%%%%%%%%% HELPING FUNCTIONS %%%%%%%%%%
+
+ function [L,S] = membership2indexed
+ % computes labeling and segmentation as indexed r-by-c images
+ % from membership r*c-by-nr_labels representation
+ [maxlab,L] = max(LABELING,[],2);
+ L(maxlab==0) = 0;
+ L = uint8(reshape(L,[r,c]));
+ [maxprob,S] = max(PROBABILITY,[],2);
+ S(sum(PROBABILITY==maxprob(:,ones(nr_labels,1)),2)>1) = 0;
+ S = uint8(reshape(S,[r,c]));
+ end
+
+ function save_displayed
+ % saves mat file with user settings and images as separate files
+ [file,path] = uiputfile('settings.mat','Save settings as');
+ if ~isequal(file,0) && ~isequal(path,0)
+ matfile = fullfile(path,file);
+ roothname = matfile(1:find(matfile=='.',1,'last')-1);
+ current_settings.method = METHOD_NAME;
+ current_settings.show = show_options{SHOW_INDEX};
+ current_settings.overwrite = OVERWRITE;
+ current_settings.regularize = regularize_options(REGULARIZE_INDEX);
+ save(matfile,'current_settings')
+ % saving displayed images
+ imwrite(LABELING_OVERLAY,[roothname,'_labels_displayed.png'])
+ imwrite(SEGMENTATION_OVERLAY,[roothname,'_results_displayed.png'])
+ [L,S] = membership2indexed;
+ imwrite(L,COLORS,[roothname,'_labels_indexed.png'])
+ imwrite(S,COLORS,[roothname,'_segmentation_indexed.png'])
+ end
+ end
+
+ function export_LS
+ % saves variables to workspace
+ % TODO, consider using:
+ % export2wsdlg({'Labeling','Segmentation'},{'gui_L','gui_S'},{L,S})
+ button = questdlg({'Exporting variables to the base workspace',...
+ 'will close texture segmentation gui and',...
+ 'might overwrite existing variables'},...
+ 'Exporting variables','OK','Cancel','OK');
+ if strcmp(button,'OK')
+ [L,S] = membership2indexed;
+ [~,dictprob] = METHOD(LABELING);
+ assignin('base','gui_L',L)
+ assignin('base','gui_S',S)
+ assignin('base','gui_dictprob',dictprob)
+ % assignin('base','gui_dictprob',T1*LABELING)
+ end
+ close(fig)
+ end
+
+ function freeze_LS
+ % feezes, closes and opens segmentation_correction_gui
+ button = questdlg({'Freezing segmentation will close texture segmentation gui',...
+ 'and open segmentation correction gui.'},...
+ 'Freezing segmentation','OK','Cancel','OK');
+ if strcmp(button,'OK')
+ [L,S] = membership2indexed;
+ close(fig)
+ segmentation_correction_gui(image_rgb,S,nr_labels,L);
+ end
+ end
+
+ function a = is_inside(x)
+ % check if x is inside image limits
+ a = inpolygon(x(1),x(2),LIMITS([1,2,2,1]),LIMITS([3,3,4,4]));
+ end
+
+ function p = get_pixel
+ % get cursor position
+ p = get(labeling_axes,'CurrentPoint');
+ p = round(p(1,[1,2]));
+ end
+
+ function show_overlays(x)
+ % overlay a circular region where the pointer is and show
+ shown_left = LABELING_OVERLAY;
+ shown_right = SEGMENTATION_OVERLAY;
+ if RADIUS>0 % thickness>0
+ P = repmat(disc(x,RADIUS,nr_circle_pts,[r,c]),[1,1,3]);
+ shown_left(P(:)) = 0.5+0.5*LABELING_OVERLAY(P(:));
+ shown_right(P(:)) = 0.5+0.5*SEGMENTATION_OVERLAY(P(:));
+ end
+ % we have to imagesc(shown) to remove overlay if needed
+ axes(labeling_axes), cla, imagesc(shown_left)
+ axes(segmentation_axes), cla, imagesc(shown_right)
+ end
+
+ function update(M)
+ % change the state of the segmentation by updating LABELING with a
+ % mask M, and updating PROBABILITY
+ labcol = zeros(1,nr_labels);
+ if LABEL>0 % not unlabeling
+ labcol(LABEL) = 1;
+ end
+ LABELING(M(:),:) = repmat(labcol,[sum(M(:)),1]);
+ if LIVE_UPDATE
+ PROBABILITY = METHOD(LABELING); % PROBABILITY computed
+ labeling_overwrite
+ regularize
+ compute_overlays(true) % compute both overlays
+ else
+ compute_overlays(false) % comput only left overlay
+ end
+ % computing overlay images
+ show_overlays(get_pixel); % showing overlay and pointer
+ end
+
+ function compute_overlays(compute_both)
+ if nargin<1
+ compute_both = true; % default computes overlay for both images
+ end
+ % computes overlays but not pointer overalay
+ % TODO: a lot of this does not need to be recalculated
+ LABELING_OVERLAY = reshape(LABELING*COLORS(2:end,:),size(image_rgb)).*...
+ (COLOR_WEIGHT+(1-COLOR_WEIGHT)*image_gray);
+ unlabeled = repmat(~any(LABELING,2),[3,1]); % pixels not labeled
+ LABELING_OVERLAY(unlabeled) = image_rgb(unlabeled);
+ if compute_both
+ if SHOW_INDEX<3 % showing segmentation or overlay
+ maxprob = max(PROBABILITY,[],2);
+ maxprobloc = PROBABILITY == maxprob(:,ones(nr_labels,1));
+ uncertain = sum(maxprobloc,2)>1; % pixels with max probability at two or more classes
+ maxprobloc(uncertain,:) = 0;
+ if SHOW_INDEX==1 % segmentation
+ SEGMENTATION_OVERLAY = reshape([uncertain,maxprobloc]*COLORS,size(image_rgb));
+ else % SHOW_INDEX==2 overlay
+ SEGMENTATION_OVERLAY = reshape([uncertain,maxprobloc]*COLORS,...
+ size(image_rgb)).*(COLOR_WEIGHT+(1-COLOR_WEIGHT)*image_gray);
+ end
+ else % 'probability x'
+ pw = SHOW_INDEX-2; % probability to show
+ minpw = min(PROBABILITY(:,pw));
+ maxpw = max(PROBABILITY(:,pw));
+ % TODO scaling should be better, relative to 1/nr_labels
+ SEGMENTATION_OVERLAY = reshape((PROBABILITY(:,pw)-minpw)/(maxpw-minpw)*[1,1,1],size(image_rgb));
+ end
+ end
+ end
+
+% TODO disc shold be saved as a list of index shifts with respect to
+% the central pixel, and change only when thickness changes
+ function M = disc(x,r,N,dim)
+ % disc shaped mask in the image
+ angles = (0:2*pi/N:2*pi*(1-1/N));
+ X = x(1)+r*cos(angles);
+ Y = x(2)+r*sin(angles);
+ M = poly2mask(X,Y,dim(1),dim(2));
+ end
+
+ function M = stadium(x1,x2,r,N,dim)
+ % stadium shaped mask in the image
+ angles = (0:2*pi/N:pi)-atan2(x1(1)-x2(1),x1(2)-x2(2));
+ X = [x1(1)+r*cos(angles), x2(1)+r*cos(angles+pi)];
+ Y = [x1(2)+r*sin(angles), x2(2)+r*sin(angles+pi)];
+ M = poly2mask(X,Y,dim(1),dim(2));
+ end
+
+ function M = fill(x)
+ [maxL,labL] = max(LABELING,[],2);
+ label_image = reshape(maxL.*labL,[r,c]);
+ M = bwselect(label_image==label_image(x(2),x(1)),x(1),x(2),4);
+ end
+
+ function [I,I_rgb,I_gray] = normalize_image(I)
+ % initialization: normalize image
+ if isa(I,'uint8')
+ I = double(I)/255;
+ end
+ if isa(I,'uint16')
+ I = double(I)/65535;
+ end
+ if size(I,3)==3 % rgb image
+ I_gray = repmat(rgb2gray(I),[1 1 3]);
+ I_rgb = I;
+ else % assuming grayscale image
+ I_gray = repmat(I,[1,1,3]);
+ I_rgb = I_gray;
+ end
+ end
+
+ function name = method_name_str
+ method_str = func2str(METHOD);
+ s1 = find(method_str==')',1,'first')+1;
+ s2 = find(method_str=='(',1,'last')-1;
+ name = method_str(s1:s2);
+ end
+
+ function n = move_once(n,total,reverse)
+ % moves option index once, respecting total number of options
+ if ~reverse
+ n = mod(n,total)+1;
+ else
+ n = mod(n+total-2,total)+1;
+ end
+ end
+
+ function labeling_overwrite
+ % labeled areas get assigned probability 1
+ if OVERWRITE
+ labeled = any(LABELING,2);
+ PROBABILITY(labeled,:) = LABELING(labeled,:); % overwritting labeled
+ end
+ end
+
+ function regularize
+ sigma = regularize_options(REGULARIZE_INDEX);
+ if sigma>0
+ filter = fspecial('gaussian',[2*ceil(sigma)+1,1],sigma);
+ PROBABILITY = reshape(PROBABILITY,[r,c,nr_labels]);
+ PROBABILITY = imfilter(PROBABILITY,filter,'replicate');
+ PROBABILITY = imfilter(PROBABILITY,filter','replicate');
+ PROBABILITY = reshape(PROBABILITY,[r*c,nr_labels]);
+ end
+ end
+
+ function r = thickness2radius(t)
+ r = t/2+0.4;
+ end
+
+ function clean_toolbar
+ set(fig,'MenuBar','none','Toolbar','figure');
+ toolbar = findall(fig,'Type','uitoolbar');
+ all_tools = allchild(toolbar);
+ % removing tools
+ for i=1:numel(all_tools)
+ if isempty(strfind(all_tools(i).Tag,'Pan'))&&...
+ isempty(strfind(all_tools(i).Tag,'Zoom'))&&...
+ isempty(strfind(all_tools(i).Tag,'SaveFigure'))&&...
+ isempty(strfind(all_tools(i).Tag,'PrintFigure'))&&...
+ isempty(strfind(all_tools(i).Tag,'DataCursor'))
+ delete(all_tools(i)) % keeping only Pan, Zoom, Save and Print
+ end
+ end
+ % adding a tool
+ [icon,~,alpha] = imread('linkaxesicon.png');
+ icon = double(icon)/255;
+ icon(alpha==0)=NaN;
+ uitoggletool(toolbar,'CData',icon,...
+ 'TooltipString','Link Axes','Tag','LinkAxes',...
+ 'OnCallback',{@link_axes,'xy'},...
+ 'OffCallback',{@link_axes,'off'});
+ % changing the order of tools
+ all_tools = allchild(toolbar);
+ set(toolbar,'children',all_tools([2,1,3:end]));
+ end
+
+ function link_axes(~,~,flag)
+ linkaxes([labeling_axes,segmentation_axes],flag)
+ end
+
+ function parsing_dictionary_input
+ % either dictionary_options, dictionary or mappings should be given
+ if isfield(dict,'patch_size') % dictionary options given
+ dictionary = build_dictionary(image,dict);
+ mappings = compute_mappings(image,dictionary);
+ T1 = mappings.T1;
+ T2 = mappings.T2;
+ elseif isfield(dict,'tree')% dictinary given
+ mappings = compute_mappings(image,dict);
+ T1 = mappings.T1;
+ T2 = mappings.T2;
+ elseif isfield(dict,'T1')% mapping given
+ T1 = dict.T1;
+ T2 = dict.T2;
+ else
+ error('Could not parse dictionary input.')
+ end
+ end
+
+ function parse_labeling_input(labeling)
+ % parsing labeling input, either labeling or labeling image
+ dim_l = size(labeling);
+ labeling = double(labeling);
+ if numel(dim_l)==3 %% rgb labeling image, to be turned to 2D image
+ max_l = max(labeling(:))+1;
+ labeling = labeling(:,:,1)+max_l*labeling(:,:,2)+max_l^2*labeling(:,:,3);
+ dim_l = size(labeling);
+ end
+ if numel(dim_l)==2 % either LABELING or 2D labeling image
+ if all(dim_l==[r*c,nr_labels]) % LABELING
+ LABELING = labeling;
+ elseif all(dim_l==[r,c]) % labeling image
+ l = unique(labeling(:)); % present labels
+ if numel(union(l,0:nr_labels))~=nr_labels+1
+ labeling_old = labeling;
+ for li = 1:numel(l)
+ labeling(labeling_old==l(li)) = li-1;
+ end
+ end
+ % now we have labels from 0 to nr_labels
+ i = (labeling(:)-1)*r*c + (1:r*c)'; % indices in LABELINGS
+ i = i(labeling(:)>0); % only labeled parts
+ LABELING = zeros(r*c,nr_labels);
+ LABELING(i) = 1;
+ else
+ error('Could not parse labeling input.')
+ end
+ else
+ error('Could not parse labeling input.')
+ end
+ end
+
+%%%%%%%%%% LABELINGS TO PROBABILITIES METHODS %%%%%%%%%%
+
+ function [probabilities,dictprob] = distributed(labelings)
+ % unlabeled pixels have label weights DISTRIBUTED equally
+ labelings(~any(labelings,2),:) = 1/nr_labels; % distributing
+ dictprob = T1*labelings;
+ probabilities = T2*dictprob; % computing probabilities
+ end
+
+ function [probabilities,dictprob] = two_max(labelings)
+ % DISTRIBUTED tresholded and repeated
+ labelings(~any(labelings,2),:) = 1/nr_labels; % distributing
+ probabilities = T2*(T1*labelings); % probabilities
+ maxprob = max(probabilities,[],2); % finding max probabilities
+ labelings_new = probabilities == maxprob(:,ones(nr_labels,1)); % new labeling is where max prob was
+ uncertain = sum(labelings_new,2)>1; % pixels with max probability at two or more classes
+ labelings_new(uncertain,:) = 1/nr_labels; % distributing at uncertain
+ dictprob = T1*labelings_new;
+ probabilities = T2*dictprob;
+ end
+
+ function [probabilities,dictprob] = two_max_over(labelings)
+ % DISTRIBUTED tresholded, overwriten and repeated
+ known_labels = any(labelings,2);
+ labelings(~known_labels,:) = 1/nr_labels; % distributing
+ probabilities = T2*(T1*labelings); % probabilities
+ maxprob = max(probabilities,[],2); % finding max probabilities
+ labelings_new = probabilities == maxprob(:,ones(nr_labels,1)); % new labeling is where max prob was
+ uncertain = sum(labelings_new,2)>1; % pixels with max probability at two or more classes
+ labelings_new(uncertain,:) = 1/nr_labels; % distributing at uncertain
+ labelings_new(known_labels,:) = labelings(known_labels,:); % OVERWRITING
+ dictprob = T1*labelings_new;
+ probabilities = T2*dictprob;
+ end
+
+ function [probabilities,dictprob] = two_cont(labelings)
+ % DISTRIBUTED repeated
+ labelings(~any(labelings,2),:) = 1/nr_labels; % distributing
+ labelings_new = T2*(T1*labelings); % probabilities
+ dictprob = T1*labelings_new;
+ probabilities = T2*dictprob;
+ end
+
+ function [probabilities,dictprob] = two_cont_over(labelings)
+ % DISTRIBUTED overwriten and repeated
+ known_labels = any(labelings,2);
+ labelings(~known_labels,:) = 1/nr_labels; % distributing
+ labelings_new = T2*(T1*labelings); % probabilities
+ labelings_new(known_labels,:) = labelings(known_labels,:); % OVERWRITING
+ dictprob = T1*labelings_new;
+ probabilities = T2*dictprob;
+ end
+
+% function probabilities = normalized(labelings)
+% probabilities = T2*(T1*labelings); % computing probabilities
+% % normalizing probabilities so that they sum to 1
+% probabilities = probabilities./repmat(sum(probabilities,2),[1,size(probabilities,2)]);
+% probabilities(isnan(probabilities)) = 0;
+% end
+%
+% function probabilities = normalized_2_max(labelings)
+% probabilities = T2*(T1*labelings); % probabilities
+% maxprob = max(probabilities,[],2); % finding max probabilities
+% labelings_new = probabilities == maxprob(:,ones(nr_labels,1)); % new labeling is where max prob was
+% uncertain = sum(labelings_new,2)>1; % pixels with max probability at two or more classes
+% labelings_new(uncertain,:) = 0; % zero at uncertain
+% probabilities = T2*(T1*labelings_new);
+% probabilities = probabilities./repmat(sum(probabilities,2),[1,size(probabilities,2)]);
+% probabilities(isnan(probabilities)) = 0;
+% end
+%
+% function probabilities = normalized_2_cont(labelings)
+% labelings_new = T2*(T1*labelings); % probabilities
+% probabilities = T2*(T1*labelings_new);
+% probabilities = probabilities./repmat(sum(probabilities,2),[1,size(probabilities,2)]);
+% probabilities(isnan(probabilities)) = 0;
+% end
+
+end
diff --git a/code/texture_gui/code/functions/linkaxesicon.png b/code/texture_gui/code/functions/linkaxesicon.png
new file mode 100755
index 0000000000000000000000000000000000000000..f1c1759a401c8c3cc0dbbb158860bcb9e18ab4ae
Binary files /dev/null and b/code/texture_gui/code/functions/linkaxesicon.png differ
diff --git a/code/texture_gui/code/functions/normalize_image.m b/code/texture_gui/code/functions/normalize_image.m
new file mode 100755
index 0000000000000000000000000000000000000000..8dd729f89488dfc3714d2fb55c1fc88dcdde0360
--- /dev/null
+++ b/code/texture_gui/code/functions/normalize_image.m
@@ -0,0 +1,9 @@
+function I = normalize_image(I)
+% we work with [0 1] images
+if isa(I,'uint8')
+ I = double(I)/255;
+end
+if isa(I,'uint16')
+ I = double(I)/65535;
+end
+end
\ No newline at end of file
diff --git a/code/texture_gui/code/functions/probability_search_km_tree.cpp b/code/texture_gui/code/functions/probability_search_km_tree.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..489281bd1148ea5b59977876e996ea9e822744a9
--- /dev/null
+++ b/code/texture_gui/code/functions/probability_search_km_tree.cpp
@@ -0,0 +1,305 @@
+/*=================================================================
+* syntax: Q = probability_search_km_tree(I, T, P, b, n); OR Q = search_km_tree(I, T, P, b);
+*
+* serach_km_tree - build probability image from dictionary and intensity image
+*
+* Input: - I: X-by-Y image
+* - T: MMl-by-K tree matrix where l is the number of layers
+* in the image (1 for grayscale and 3 for RGB)
+ - P: MMc-by-K probability matrix where c is the number of
+ classes
+* - b: branching factor
+* - n: normalization (true or false), defaults to false
+*
+* Output: - Q: X-by-Y-by-c probability image
+*
+* Author: Anders Dahl, abda@dtu.dk, February 2016.
+*=================================================================*/
+
+#include "mex.h"
+#include <stdio.h>
+#include <math.h>
+#include "matrix.h"
+#include <vector>
+
+#include <iostream>
+using namespace std;
+
+// struct for the tree
+struct tree_st
+{
+ double *tree_data;
+ double *probability_data;
+ int n_dim, n_dprob, n_nodes, branch_fac, M, Mh, n_cls;
+ vector<int> idx_disp;
+};
+
+// struct for image
+struct im_st
+{
+ double *im_data;
+ int rows, cols, layers, n_pix;
+};
+
+// estimate the distance between a vector in tree given by the node and a
+// patch in the image given by the row and column
+double get_dist(vector<double>& patch, tree_st& tree, int& node)
+{
+ double d = 0, tmp;
+ int id_t = tree.n_dim*node;
+
+ for ( int i = 0; i < tree.n_dim; i++ ){
+ tmp = patch[i] - *(tree.tree_data + id_t);
+ id_t++;
+ d += tmp*tmp;
+ }
+
+ return d;
+}
+
+
+// Function for sampling patches from the image into the patch arrays
+// inputs reference to the image struct, tree struct, patch struct and position of the sampling coordinate.
+// There is no check if the sampling is outside the image
+// vector<double> sample_patch(im_st& im, int& M, int r_im, int c_im, bool& normalize)
+vector<double> sample_patch(im_st& im, int& M, int r_im, int c_im, bool& normalize)
+{
+ int id_l, id_r, id_i; // iterators for looking up image data
+ int id_p = 0; // iterator for looking up patch data
+ double sum_sq = 0, pix_val; // variables for normalization
+ int n_dim = M*M*im.layers; // number of dimensions computed here, becasue tree is not included
+ vector<double> patch(n_dim);
+ int Mh = (M-1)/2;
+
+ for ( int l = 0; l < im.layers; l++ ){ // image is sampled by three nested loops (layers, columns, rows)
+ id_l = im.n_pix*l;
+ for ( int i = c_im-Mh; i <= c_im+Mh; i++ ){
+ id_r = id_l + i*im.rows;
+ for ( int j = r_im-Mh; j <= r_im+Mh; j++ ){
+ id_i = id_r + j;
+ pix_val = *(im.im_data + id_i);
+ patch[id_p] = pix_val;
+ sum_sq += pix_val*pix_val; // sum of squares for normalization
+ id_p++;
+ }
+ }
+ }
+ if ( normalize ){ // if the patch should be normalized to unit length
+ double inv_sq = 1;
+ if ( sum_sq > 0 ){
+ inv_sq = 1/sqrt(sum_sq); // inverse sum of squares
+ }
+ for ( int i = 0; i < n_dim; i++ ){
+ patch[i] = patch[i]*inv_sq;
+ }
+ }
+ return patch;
+}
+
+
+// The tree search function
+int search_tree(im_st& im, tree_st& tree, int& r, int& c, bool& normalize)
+{
+ int node = 0, node_min = -1, node_min_level, next_node; // variables for searching the tree
+ double d_min = 10e100, d, d_min_level;
+
+ vector<double> patch = sample_patch(im, tree.M, c, r, normalize); // get the pixel values in a patch
+ while ( node < tree.n_nodes ){ // go through the tree
+ if ( *(tree.tree_data + node*tree.n_dim) == -1 ){ // check if node is a leaf-node
+ return node_min;
+ }
+
+ d_min_level = 10e100; // set minimum distance to high value
+ for ( int i = 0; i < tree.branch_fac; i++ ){ // go through nodes at level
+ next_node = node + i;
+ d = get_dist(patch, tree, next_node);
+
+ if ( d < d_min_level ){ // set current node to the minimum distance
+ d_min_level = d;
+ node_min_level = next_node;
+ }
+ }
+ if ( d_min_level < d_min ){ // set overall minimum distance and minimum node
+ d_min = d_min_level;
+ node_min = node_min_level;
+ }
+ node = (node_min_level+1)*tree.branch_fac; // go to the child node
+ }
+ return node_min;
+}
+
+void add_probability(tree_st& tree, double *Q, int& idx, int& node){
+ int id_node = node*tree.n_dprob;
+ for ( int i = 0; i < tree.idx_disp.size(); i++ ){
+ *(Q + idx + tree.idx_disp[i]) += *(tree.probability_data + id_node + i);
+ }
+}
+
+
+void add_count(vector<int>& ct_im, int& col, int& row, int& rows, tree_st& tree){
+ for ( int i = col-tree.Mh; i <= col+tree.Mh; i++ ){
+ for ( int j = row-tree.Mh; j <= row+tree.Mh; j++ ){
+ ct_im[i*rows+j]++;
+ }
+ }
+}
+
+// The tree search function applied to the entire image - border is zero and interior is in 1,...,n
+void probability_search_image(im_st& im, tree_st& tree, double *Q, bool& normalize)
+{
+ vector<int> ct_im;
+ ct_im.resize(im.rows*im.cols);
+ int node;
+ int idx = tree.Mh*im.rows; // increase with empty rows at border
+ for ( int i = tree.Mh; i < im.cols-tree.Mh; i++ ){
+ idx += tree.Mh; // first Mh pixels are border
+ for ( int j = tree.Mh; j < im.rows-tree.Mh; j++ ){
+ node = search_tree(im, tree, i, j, normalize); // find assignment
+ add_probability(tree, Q, idx, node);
+ idx++;
+ add_count(ct_im, i, j, im.rows, tree);
+ }
+ idx += tree.Mh; // last Mh pixels are border
+ }
+
+ int l_id;
+ for ( int j = 0; j < tree.n_cls; j++) {
+ l_id = j*im.rows*im.cols;
+ for ( int i = 0; i < ct_im.size(); i++ ){
+ *(Q+i+l_id) = *(Q+i+l_id)/(double)ct_im[i];
+ }
+ }
+}
+
+
+// The gateway routine
+void mexFunction( int nlhs, mxArray *plhs[],
+ int nrhs, const mxArray *prhs[])
+{
+ // input image (I), tree (tree), probabilities (P) and output probabilities (Q)
+ double *I, *tree, *P, *Q;
+ int b, M, ndim, ndtree, ndprob;
+ const int *dim, *dtree, *dprob;
+ bool normalize = false;
+ /* Check for proper number of arguments. */
+ /* NOTE: You do not need an else statement when using
+ mexErrMsgTxt within an if statement. It will never
+ get to the else statement if mexErrMsgTxt is executed.
+ (mexErrMsgTxt breaks you out of the MEX-file.)
+ */
+ if(nrhs < 4 || nrhs > 5)
+ mexErrMsgTxt("Four or five inputs required.");
+ if(nlhs != 1)
+ mexErrMsgTxt("One output required.");
+
+ // Create a pointer to the input matrix.
+ I = mxGetPr(prhs[0]);
+ tree = mxGetPr(prhs[1]);
+ P = mxGetPr(prhs[2]);
+
+ double *bd;
+ bd = mxGetPr(prhs[3]);
+ b = (int)bd[0];
+
+ if ( nrhs == 5 ){
+ bool *normalize_d;
+ normalize_d = (bool *)mxGetData(prhs[4]);
+ normalize = normalize_d[0];
+ }
+
+ if ( b < 1 )
+ mexErrMsgTxt("b must be positive.");
+
+ // Get the dimensions of the matrix input.
+ ndim = mxGetNumberOfDimensions(prhs[0]);
+ if (ndim != 2 && ndim != 3)
+ mexErrMsgTxt("probability_search_km_tree only works for 2-dimensional or 3-dimensional images.");
+
+ ndtree = mxGetNumberOfDimensions(prhs[1]);
+ if (ndtree != 2)
+ mexErrMsgTxt("probability_search_km_tree requires tree to be a matrix (2-dimensional).");
+
+ ndprob = mxGetNumberOfDimensions(prhs[2]);
+ if (ndprob != 2)
+ mexErrMsgTxt("probability_search_km_tree requires probability matrix to be a matrix (2-dimensional).");
+
+ dim = mxGetDimensions(prhs[0]);
+ dtree = mxGetDimensions(prhs[1]);
+ dprob = mxGetDimensions(prhs[2]);
+
+ if ( dtree[1] != dprob[1] )
+ mexErrMsgTxt("Tree matrix and probability matrix must have the same number of columns.");
+
+ if ( ndim == 3 )
+ {
+ M = (int)sqrt((double)dtree[0]/(double)dim[2]);
+ }
+ else
+ {
+ M = (int)sqrt((double)dtree[0]);
+ }
+
+ if ( 1 - (M % 2) || M < 1)
+ mexErrMsgTxt("M must be odd and positive.");
+
+
+ // image struct
+ im_st Im;
+ Im.im_data = I;
+ Im.rows = dim[0];
+ Im.cols = dim[1];
+ if ( ndim == 3 )
+ {
+ Im.layers = dim[2];
+ }
+ else
+ {
+ Im.layers = 1;
+ }
+ Im.n_pix = Im.rows*Im.cols;
+
+ // tree struct
+ tree_st Tree;
+ Tree.tree_data = tree;
+ Tree.probability_data = P;
+ Tree.n_dim = dtree[0];
+ Tree.n_nodes = dtree[1];
+ Tree.n_dprob = dprob[0];
+ Tree.branch_fac = b;
+ Tree.M = M;
+ Tree.Mh = (int)(0.5*(double)(M-1.0));
+ // Number of classes
+ Tree.n_cls = dprob[0]/(M*M);
+ // Compute relative displacement vector
+ Tree.idx_disp.reserve(dprob[0]);
+ int id;
+ for ( int k = 0; k < Tree.n_cls; k++ ){
+ id = -Tree.Mh*Im.rows - Tree.Mh + k*Im.rows*Im.cols;
+ for ( int i = 0; i < M; i++ ){
+ for ( int j = 0; j < M; j++ ){
+ Tree.idx_disp.push_back(id);
+ id++;
+ }
+ id+=Im.rows-M;
+ }
+ }
+// for ( int i = 0; i < Tree.idx_disp.size(); i++ ){
+// cout << Tree.idx_disp[i] << " ";
+// }
+// cout << endl;
+
+
+ if ( M*M*Im.layers != Tree.n_dim )
+ mexErrMsgTxt("Dimensions of the tree and the image does not fit.");
+
+ int ndpout[3] = {Im.rows, Im.cols, Tree.n_cls};
+
+ // Set the output pointer to the output matrix. Array initialized to zero.
+ plhs[0] = mxCreateNumericArray(3, ndpout, mxDOUBLE_CLASS, mxREAL);
+
+ // Create a C pointer to a copy of the output matrix.
+ Q = mxGetPr(plhs[0]);
+ // Search the tree using the C++ subroutine
+ probability_search_image(Im, Tree, Q, normalize);
+}
+
diff --git a/code/texture_gui/code/functions/probability_search_km_tree.mexa64 b/code/texture_gui/code/functions/probability_search_km_tree.mexa64
new file mode 100755
index 0000000000000000000000000000000000000000..80a9ec4cb9d02af3b31b2d1b6f0df54720cd8d27
Binary files /dev/null and b/code/texture_gui/code/functions/probability_search_km_tree.mexa64 differ
diff --git a/code/texture_gui/code/functions/probability_search_km_tree.mexmaci64 b/code/texture_gui/code/functions/probability_search_km_tree.mexmaci64
new file mode 100755
index 0000000000000000000000000000000000000000..86ce9e57e1823c0b786e2bedbfb5a8fbdd0aada4
Binary files /dev/null and b/code/texture_gui/code/functions/probability_search_km_tree.mexmaci64 differ
diff --git a/code/texture_gui/code/functions/probability_search_km_tree.mexw64 b/code/texture_gui/code/functions/probability_search_km_tree.mexw64
new file mode 100755
index 0000000000000000000000000000000000000000..14ed3e6c6fea7689c9c02ffc4a86221931b35546
Binary files /dev/null and b/code/texture_gui/code/functions/probability_search_km_tree.mexw64 differ
diff --git a/code/texture_gui/code/functions/process_image.m b/code/texture_gui/code/functions/process_image.m
new file mode 100755
index 0000000000000000000000000000000000000000..6dce460669eb429803cf9132fe20e5957ae30045
--- /dev/null
+++ b/code/texture_gui/code/functions/process_image.m
@@ -0,0 +1,16 @@
+function [S,P] = process_image(im,dictionary)
+% segmentation and probabilities from and image and a dictionary
+
+im = normalize_image(im);
+P = probability_search_km_tree(im, ...
+ dictionary.tree, ...
+ dictionary.dictprob, ...
+ dictionary.options.branching_factor, ...
+ dictionary.options.normalization);
+
+[maxprob,S] = max(P,[],3);
+uncertain = sum(P == maxprob(:,:,ones(size(P,3),1)),3)>1;
+S(uncertain) = 0;
+
+end
+
diff --git a/code/texture_gui/code/functions/search_km_tree.cpp b/code/texture_gui/code/functions/search_km_tree.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..3b9fcbcb4adc6a2c885724f29f1dcdc24fcbd425
--- /dev/null
+++ b/code/texture_gui/code/functions/search_km_tree.cpp
@@ -0,0 +1,239 @@
+/*=================================================================
+* syntax: A = search_km_tree(I, T, b, n); OR A = search_km_tree(I, T, b);
+*
+* serach_km_tree - build assignment image from intensity image
+*
+* Input: - I: X-by-Y image
+* - T: MMl-by-K tree matrix where l is the number of layers
+* in the image (1 for grayscale and 3 for RGB)
+* - b: brancing factor
+* - n: normalization (true or false), defaults to false
+*
+* Output: - A: X-by-Y assignment matrix
+*
+* Author: Anders Dahl, abda@dtu.dk, december 2015.
+*=================================================================*/
+
+#include "mex.h"
+#include <stdio.h>
+#include <math.h>
+#include "matrix.h"
+#include <vector>
+
+#include <iostream>
+using namespace std;
+
+// struct for the tree
+struct tree_st
+{
+ double *tree_data;
+ int n_dim, n_nodes, branch_fac, M, Mh;
+};
+
+// struct for image
+struct im_st
+{
+ double *im_data;
+ int rows, cols, layers, n_pix;
+};
+
+// estimate the distance between a vector in tree given by the node and a
+// patch in the image given by the row and column
+double get_dist(vector<double>& patch, tree_st& tree, int& node)
+{
+ double d = 0, tmp;
+ int id_t = tree.n_dim*node;
+
+ for ( int i = 0; i < tree.n_dim; i++ ){
+ tmp = patch[i] - *(tree.tree_data + id_t);
+ id_t++;
+ d += tmp*tmp;
+ }
+
+ return d;
+}
+
+
+// Function for sampling patches from the image into the patch arrays
+// inputs reference to the image struct, tree struct, patch struct and position of the sampling coordinate.
+// There is no check if the sampling is outside the image
+// vector<double> sample_patch(im_st& im, int& M, int r_im, int c_im, bool& normalize)
+vector<double> sample_patch(im_st& im, int& M, int r_im, int c_im, bool& normalize)
+{
+ int id_l, id_r, id_i; // iterators for looking up image data
+ int id_p = 0; // iterator for looking up patch data
+ double sum_sq = 0, pix_val; // variables for normalization
+ int n_dim = M*M*im.layers; // number of dimensions computed here, becasue tree is not included
+ vector<double> patch(n_dim);
+ int Mh = (M-1)/2;
+
+ for ( int l = 0; l < im.layers; l++ ){ // image is sampled by three nested loops (layers, columns, rows)
+ id_l = im.n_pix*l;
+ for ( int i = c_im-Mh; i <= c_im+Mh; i++ ){
+ id_r = id_l + i*im.rows;
+ for ( int j = r_im-Mh; j <= r_im+Mh; j++ ){
+ id_i = id_r + j;
+ pix_val = *(im.im_data + id_i);
+ patch[id_p] = pix_val;
+ sum_sq += pix_val*pix_val; // sum of squares for normalization
+ id_p++;
+ }
+ }
+ }
+ if ( normalize ){ // if the patch should be normalized to unit length
+ double inv_sq = 1;
+ if ( sum_sq > 0 ){
+ inv_sq = 1/sqrt(sum_sq); // inverse sum of squares
+ }
+ for ( int i = 0; i < n_dim; i++ ){
+ patch[i] = patch[i]*inv_sq;
+ }
+ }
+ return patch;
+}
+
+
+// The tree search function
+int search_tree(im_st& im, tree_st& tree, int& r, int& c, bool& normalize)
+{
+ int node = 0, node_min = -1, node_min_level, next_node; // variables for searching the tree
+ double d_min = 10e100, d, d_min_level;
+
+ vector<double> patch = sample_patch(im, tree.M, c, r, normalize); // get the pixel values in a patch
+ while ( node < tree.n_nodes ){ // go through the tree
+ if ( *(tree.tree_data + node*tree.n_dim) == -1 ){ // check if node is a leaf-node
+ return node_min;
+ }
+
+ d_min_level = 10e100; // set minimum distance to high value
+ for ( int i = 0; i < tree.branch_fac; i++ ){ // go through nodes at level
+ next_node = node + i;
+ d = get_dist(patch, tree, next_node);
+
+ if ( d < d_min_level ){ // set current node to the minimum distance
+ d_min_level = d;
+ node_min_level = next_node;
+ }
+ }
+ if ( d_min_level < d_min ){ // set overall minimum distance and minimum node
+ d_min = d_min_level;
+ node_min = node_min_level;
+ }
+ node = (node_min_level+1)*tree.branch_fac; // go to the child node
+ }
+ return node_min;
+}
+
+// The tree search function applied to the entire image - border is zero and interior is in 1,...,n
+void search_image(im_st& im, tree_st& tree, double *A, bool& normalize)
+{
+ int idx = tree.Mh*im.rows; // increase with empty rows at border
+ for ( int i = tree.Mh; i < im.cols-tree.Mh; i++ ){
+ idx += tree.Mh; // first Mh pixels are border
+ for ( int j = tree.Mh; j < im.rows-tree.Mh; j++ ){
+ *(A + idx) = search_tree(im, tree, i, j, normalize) + 1; // find assignment
+ idx++;
+ }
+ idx += tree.Mh; // last Mh pixels are border
+ }
+}
+
+
+// The gateway routine
+void mexFunction( int nlhs, mxArray *plhs[],
+ int nrhs, const mxArray *prhs[])
+{
+ // input image (I), tree (tree) and output assignment (A)
+ double *I, *A, *tree;
+ int b, M, ndim, ndtree;
+ const int *dim, *dtree;
+ bool normalize = false;
+ /* Check for proper number of arguments. */
+ /* NOTE: You do not need an else statement when using
+ mexErrMsgTxt within an if statement. It will never
+ get to the else statement if mexErrMsgTxt is executed.
+ (mexErrMsgTxt breaks you out of the MEX-file.)
+ */
+ if(nrhs < 3 || nrhs > 4)
+ mexErrMsgTxt("Three or four inputs required.");
+ if(nlhs != 1)
+ mexErrMsgTxt("One output required.");
+
+ // Create a pointer to the input matrix.
+ I = mxGetPr(prhs[0]);
+ tree = mxGetPr(prhs[1]);
+
+ double *bd;
+ bd = mxGetPr(prhs[2]);
+ b = (int)bd[0];
+
+ if ( nrhs == 4 ){
+ bool *normalize_d;
+ normalize_d = (bool *)mxGetData(prhs[3]);
+ normalize = normalize_d[0];
+ }
+
+ if ( b < 1 )
+ mexErrMsgTxt("b must be positive.");
+
+ // Get the dimensions of the matrix input.
+ ndim = mxGetNumberOfDimensions(prhs[0]);
+ if (ndim != 2 && ndim != 3)
+ mexErrMsgTxt("search_km_tree only works for 2-dimensional or 3-dimensional images.");
+
+ ndtree = mxGetNumberOfDimensions(prhs[1]);
+ if (ndtree != 2)
+ mexErrMsgTxt("search_km_tree only works for 2-dimensional tree.");
+
+ dim = mxGetDimensions(prhs[0]);
+ dtree = mxGetDimensions(prhs[1]);
+
+ if ( ndim == 3 )
+ {
+ M = (int)sqrt((double)dtree[0]/(double)dim[2]);
+ }
+ else
+ {
+ M = (int)sqrt((double)dtree[0]);
+ }
+
+ if ( 1 - (M % 2) || M < 1)
+ mexErrMsgTxt("M must be odd and positive.");
+
+
+ // tree struct
+ tree_st Tree;
+ Tree.tree_data = tree;
+ Tree.n_dim = dtree[0];
+ Tree.n_nodes = dtree[1];
+ Tree.branch_fac = b;
+ Tree.M = M;
+ Tree.Mh = (int)(0.5*(double)(M-1.0));
+
+ // image struct
+ im_st Im;
+ Im.im_data = I;
+ Im.rows = dim[0];
+ Im.cols = dim[1];
+ if ( ndim == 3 )
+ {
+ Im.layers = dim[2];
+ }
+ else
+ {
+ Im.layers = 1;
+ }
+ Im.n_pix = Im.rows*Im.cols;
+
+ if ( M*M*Im.layers != Tree.n_dim )
+ mexErrMsgTxt("Dimensions of the tree and the image does not fit.");
+
+ // Set the output pointer to the output matrix. Array initialized to zero.
+ plhs[0] = mxCreateNumericArray(ndtree, dim, mxDOUBLE_CLASS, mxREAL);
+
+ // Create a C pointer to a copy of the output matrix.
+ A = mxGetPr(plhs[0]);
+ // Search the tree using the C++ subroutine
+ search_image(Im, Tree, A, normalize);
+}
+
diff --git a/code/texture_gui/code/functions/search_km_tree.mexa64 b/code/texture_gui/code/functions/search_km_tree.mexa64
new file mode 100755
index 0000000000000000000000000000000000000000..1c413daea71625531c3473517f3c177df8cadac0
Binary files /dev/null and b/code/texture_gui/code/functions/search_km_tree.mexa64 differ
diff --git a/code/texture_gui/code/functions/search_km_tree.mexmaci64 b/code/texture_gui/code/functions/search_km_tree.mexmaci64
new file mode 100755
index 0000000000000000000000000000000000000000..f6786b110627679ad91fbfd91cb9022bdf2b1a71
Binary files /dev/null and b/code/texture_gui/code/functions/search_km_tree.mexmaci64 differ
diff --git a/code/texture_gui/code/functions/search_km_tree.mexw64 b/code/texture_gui/code/functions/search_km_tree.mexw64
new file mode 100755
index 0000000000000000000000000000000000000000..cbe6d41e8ffdb3c5057845234410e9daeeebe734
Binary files /dev/null and b/code/texture_gui/code/functions/search_km_tree.mexw64 differ
diff --git a/code/texture_gui/code/functions/search_km_tree_xcorr.cpp b/code/texture_gui/code/functions/search_km_tree_xcorr.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..20c722b6018fe6907bb0ed9d4c0dbe87a133bc0e
--- /dev/null
+++ b/code/texture_gui/code/functions/search_km_tree_xcorr.cpp
@@ -0,0 +1,239 @@
+#include "mex.h"
+#include <stdio.h>
+#include <math.h>
+#include "matrix.h"
+#include <vector>
+
+#include <iostream>
+using namespace std;
+
+/*
+ * This is a MEX-file for MATLAB.
+ *
+ * Anders Bjorholm Dahl, 8/12-2015
+ *
+ *
+ */
+
+// struct for the tree
+struct tree_st
+{
+ double *tree_data;
+ int n_dim, inv_ndim, n_nodes, branch_fac, M, Mh;
+};
+
+// struct for image
+struct im_st
+{
+ double *im_data;
+ int rows, cols, layers, n_pix;
+};
+
+// estimate the distance between a vector in tree given by the node and a
+// patch in the image given by the row and column
+double get_corr(vector<double>& patch, tree_st& tree, int& node)
+{
+ double corr = 0, tmp;
+ int id_t = tree.n_dim*node;
+
+ for ( int i = 0; i < tree.n_dim; i++ ){
+ corr += patch[i]*(*(tree.tree_data + id_t));
+ id_t++;
+ }
+
+ return corr;
+}
+
+// nomalize patch vectors to zero mean and unit Euclidean length
+void norm_patch(vector<double>& patch, double& mu, int n_dim){
+
+ double sum_sq = 0;
+
+ for ( int i = 0; i < n_dim; i++ ){
+ patch[i] = patch[i] - mu;
+ sum_sq += patch[i]*patch[i];
+ }
+
+ double inv_sq = 1;
+ if ( sum_sq > 0 ){
+ inv_sq = 1.0/sqrt(sum_sq); // divide by sum of squares
+ }
+
+ for ( int i = 0; i < n_dim; i++ ){
+ patch[i] = patch[i]*inv_sq;
+ }
+}
+
+
+// Function for sampling patches from the image into the patch arrays
+// inputs reference to the image struct, tree struct, patch struct and position of the sampling coordinate.
+// There is no check if the sampling is outside the image
+// vector<double> sample_patch(im_st& im, int& M, int r_im, int c_im, bool& normalize)
+vector<double> sample_patch(im_st& im, int& M, int r_im, int c_im)
+{
+ int id_l, id_r, id_i; // iterators for looking up image data
+ int id_p = 0; // iterator for looking up patch data
+ double sum = 0, pix_val, mu; // variables for normalization
+ int n_dim = M*M*im.layers; // number of dimensions computed here, becasue tree is not included
+ vector<double> patch(n_dim);
+ int Mh = (M-1)/2;
+
+ for ( int l = 0; l < im.layers; l++ ){ // image is sampled by three nested loops (layers, columns, rows)
+ id_l = im.n_pix*l;
+ for ( int i = c_im-Mh; i <= c_im+Mh; i++ ){
+ id_r = id_l + i*im.rows;
+ for ( int j = r_im-Mh; j <= r_im+Mh; j++ ){
+ id_i = id_r + j;
+ pix_val = *(im.im_data + id_i);
+ patch[id_p] = pix_val;
+ sum += pix_val; // sum for estimating the mean
+ id_p++;
+ }
+ }
+ }
+
+ mu = sum/((double)n_dim); // mean value
+ norm_patch(patch, mu, n_dim); // nomalization
+
+ return patch;
+}
+
+
+// The tree search function
+int search_tree(im_st& im, tree_st& tree, int& r, int& c)
+{
+ int node = 0, node_max = -1, node_max_level, next_node; // variables for searching the tree
+ double corr_max = -1, corr, corr_max_level;
+
+ vector<double> patch = sample_patch(im, tree.M, c, r); // get the pixel values in a patch
+ while ( node < tree.n_nodes ){ // go through the tree
+ if ( *(tree.tree_data + node*tree.n_dim) == -1 ){ // check if node is a leaf-node
+ return node_max;
+ }
+
+ corr_max_level = -1; // set correlation to low value
+ for ( int i = 0; i < tree.branch_fac; i++ ){ // go through nodes at level
+ next_node = node + i;
+ corr = get_corr(patch, tree, next_node);
+
+ if ( corr > corr_max_level ){ // set current node to the maximum correlation
+ corr_max_level = corr;
+ node_max_level = next_node;
+ }
+ }
+ if ( corr_max_level > corr_max ){ // set overall maximum correlation and corresponding node
+ corr_max = corr_max_level;
+ node_max = node_max_level;
+ }
+ node = (node_max_level+1)*tree.branch_fac; // go to the child node
+ }
+ return node_max;
+}
+
+// The tree search function applied to the entire image - border is zero and interior is in 1,...,n
+void search_image(im_st& im, tree_st& tree, double *A)
+{
+ int idx = tree.Mh*im.rows; // increase with empty rows at border
+ for ( int i = tree.Mh; i < im.cols-tree.Mh; i++ ){
+ idx += tree.Mh; // first Mh pixels are border
+ for ( int j = tree.Mh; j < im.rows-tree.Mh; j++ ){
+ *(A + idx) = search_tree(im, tree, i, j) + 1; // find assignment
+ idx++;
+ }
+ idx += tree.Mh; // last Mh pixels are border
+ }
+}
+
+
+// The gateway routine
+void mexFunction( int nlhs, mxArray *plhs[],
+ int nrhs, const mxArray *prhs[])
+{
+ // input image (I), tree (tree) and output assignment (A)
+ double *I, *A, *tree;
+ int b, M, ndim, ndtree;
+ const int *dim, *dtree;
+ /* Check for proper number of arguments. */
+ /* NOTE: You do not need an else statement when using
+ mexErrMsgTxt within an if statement. It will never
+ get to the else statement if mexErrMsgTxt is executed.
+ (mexErrMsgTxt breaks you out of the MEX-file.)
+ */
+ if(nrhs != 3 )
+ mexErrMsgTxt("Three inputs required.");
+ if(nlhs != 1)
+ mexErrMsgTxt("One output required.");
+
+ // Create a pointer to the input matrix.
+ I = mxGetPr(prhs[0]);
+ tree = mxGetPr(prhs[1]);
+
+ double *bd;
+ bd = mxGetPr(prhs[2]);
+ b = (int)bd[0];
+
+ if ( b < 1 )
+ mexErrMsgTxt("b must be positive.");
+
+ // Get the dimensions of the matrix input.
+ ndim = mxGetNumberOfDimensions(prhs[0]);
+ if (ndim != 2 && ndim != 3)
+ mexErrMsgTxt("search_km_tree only works for 2-dimensional or 3-dimensional images.");
+
+ ndtree = mxGetNumberOfDimensions(prhs[1]);
+ if (ndtree != 2)
+ mexErrMsgTxt("search_km_tree only works for 2-dimensional tree.");
+
+ dim = mxGetDimensions(prhs[0]);
+ dtree = mxGetDimensions(prhs[1]);
+
+ if ( ndim == 3 )
+ {
+ M = (int)sqrt((double)dtree[0]/(double)dim[2]);
+ }
+ else
+ {
+ M = (int)sqrt((double)dtree[0]);
+ }
+
+ if ( 1 - (M % 2) || M < 1)
+ mexErrMsgTxt("M must be odd and positive.");
+
+
+ // tree struct
+ tree_st Tree;
+ Tree.tree_data = tree;
+ Tree.n_dim = dtree[0];
+ Tree.inv_ndim = 1/((double)Tree.n_dim);
+ Tree.n_nodes = dtree[1];
+ Tree.branch_fac = b;
+ Tree.M = M;
+ Tree.Mh = (int)(0.5*(double)(M-1.0));
+
+ // image struct
+ im_st Im;
+ Im.im_data = I;
+ Im.rows = dim[0];
+ Im.cols = dim[1];
+ if ( ndim == 3 )
+ {
+ Im.layers = dim[2];
+ }
+ else
+ {
+ Im.layers = 1;
+ }
+ Im.n_pix = Im.rows*Im.cols;
+
+ if ( M*M*Im.layers != Tree.n_dim )
+ mexErrMsgTxt("Dimensions of the tree and the image does not fit.");
+
+ // Set the output pointer to the output matrix. Array initialized to zero.
+ plhs[0] = mxCreateNumericArray(ndtree, dim, mxDOUBLE_CLASS, mxREAL);
+
+ // Create a C pointer to a copy of the output matrix.
+ A = mxGetPr(plhs[0]);
+ // Search the tree using the C++ subroutine
+ search_image(Im, Tree, A);
+}
+
diff --git a/code/texture_gui/code/functions/search_km_tree_xcorr.mexa64 b/code/texture_gui/code/functions/search_km_tree_xcorr.mexa64
new file mode 100755
index 0000000000000000000000000000000000000000..0ac82cdef7ccd8a4ca18fcd4b65ed8df7fd7ec13
Binary files /dev/null and b/code/texture_gui/code/functions/search_km_tree_xcorr.mexa64 differ
diff --git a/code/texture_gui/code/functions/search_km_tree_xcorr.mexmaci64 b/code/texture_gui/code/functions/search_km_tree_xcorr.mexmaci64
new file mode 100755
index 0000000000000000000000000000000000000000..5f224caefe0987623c73d07353e150a8b390ada6
Binary files /dev/null and b/code/texture_gui/code/functions/search_km_tree_xcorr.mexmaci64 differ
diff --git a/code/texture_gui/code/functions/search_km_tree_xcorr.mexw64 b/code/texture_gui/code/functions/search_km_tree_xcorr.mexw64
new file mode 100755
index 0000000000000000000000000000000000000000..8d771be06e99b9fa7ee3d64aa69f4f7669f16de2
Binary files /dev/null and b/code/texture_gui/code/functions/search_km_tree_xcorr.mexw64 differ
diff --git a/code/texture_gui/code/functions/segmentation_correction_gui.m b/code/texture_gui/code/functions/segmentation_correction_gui.m
new file mode 100755
index 0000000000000000000000000000000000000000..6df8c878f3ef343279513f8cf5f7cd0009c45ea1
--- /dev/null
+++ b/code/texture_gui/code/functions/segmentation_correction_gui.m
@@ -0,0 +1,406 @@
+function segmentation_correction_gui(image,segmentation,nr_labels,labeling)
+%IMAGE_TEXTURE_GUI Interactive segmentation correction
+
+%%%%%%%%%% DEFAULTS %%%%%%%%%%
+
+if nargin<1 % default example image
+ %image = imread('bag.png');
+ image = imread('football.jpg');
+end
+[r,c,~] = size(image);
+if nargin<2 || isempty(segmentation)% default segmentation
+ segmentation = zeros(r,c,'uint8');
+else
+ segmentation = uint8(segmentation);
+end
+if nargin<3 || isempty(nr_labels) % defalult number of labels
+ nr_labels = double(max(2,max(segmentation(:))));
+end
+if nargin<4 || isempty(labeling) % default, unlabeled initial correction
+ % TODO allow corrections to be inferred from labelings
+ DRAWING = zeros(r,c,'uint8');
+else
+ DRAWING = uint8(labeling);
+end
+
+%%%%%%%%%% SETTINGS %%%%%%%%%%
+
+% labels
+LABEL = 1; % initial label is 1
+
+% thickness
+thickness_options = [1 2 3 4 5 10 20 30 40 50 100 -1]; % last option (value -1) is 'fill'
+thickness_options_string = num2cell(thickness_options);
+thickness_options_string{end} = 'fill';
+THICKNESS_INDEX = 5; % initial pencil thickness is the fifth option
+RADIUS = thickness2radius(thickness_options(THICKNESS_INDEX)); % pencil radius
+
+% show
+show_options(1:2) = {'segmentation','overlay'};
+SHOW_INDEX = 1;
+
+% colormap
+colormap_options = {@(x)[0.5,0.5,0.5;cool(x)], @(x)[0.5,0.5,0.5;spring(x)],...
+ @(x)[0.5,0.5,0.5;parula(x)]}; % gray color for unlabeled
+COLORMAP_INDEX = 1; % current colormap
+COLORS = colormap_options{COLORMAP_INDEX}(nr_labels); % visualization colors for label overlay
+% other settings
+color_weight = 0.4; % color weight for label overlay
+nr_circle_pts = 16; % number of points defining a circular pencil
+
+%%%%%%%%%% INITIALIZATION AND LAYOUT %%%%%%%%%%
+[image,gray] = normalize_image(image); % impose 0-to-1 rgb
+CORRECTED = uint8(segmentation);
+CORRECTED(DRAWING~=0) = DRAWING(DRAWING~=0); % if initial corrections given
+DRAWING_OVERLAY = image; % image overlaid labeling
+CORRECTED_OVERLAY = ind2rgb(CORRECTED,COLORS);
+
+fmar = [0.2 0.2]; % discance from screen edge to figure (x and y)
+amar = [0.05 0.05]; % margin around axes, relative to figure
+my = 0.8:0.05:0.9; % menu items y position
+mh = 0.03; % menu items height
+cw = (0.4-0.3)/(nr_labels+1); % colorcube width
+cx = 0.3:cw:0.4; % colorcubes x position
+
+fig = figure('Units','Normalized','Position',[fmar,1-2*fmar],...
+ 'Pointer','watch','KeyPressFcn',@key_press,'InvertHardCopy', 'off',...
+ 'Name','Segmentation correction GUI');
+clean_toolbar
+
+labeling_axes = axes('Units','Normalized','Position',[0,0,0.5,0.8]+[amar,-2*amar]);
+imagesc(DRAWING_OVERLAY), axis image off, hold on
+segmentation_axes = axes('Units','Normalized','Position',[0.5,0,0.5,0.8]+[amar,-2*amar]);
+imagesc(CORRECTED_OVERLAY,[0,nr_labels]), axis image off, hold on
+
+uicontrol('String','Label [L] : ','Style','text','HorizontalAlignment','right',...
+ 'Units','Normalized','Position',[0,my(3),0.25,mh]);
+labels_text = uicontrol('String',num2str(LABEL),...
+ 'BackgroundColor',COLORS(LABEL+1,:),...
+ 'Style','text','HorizontalAlignment','left',...
+ 'Units','Normalized','Position',[0.25,my(3),0.03,mh]);
+label_pointer = cell(nr_labels+1,1);
+label_colorcubes = cell(nr_labels+1,1);
+label_char = repmat(' ',[1,nr_labels+1]);
+label_char(LABEL+1)='|';
+for k = 1:nr_labels+1
+ label_colorcubes{k} = uicontrol('String',' ','Style','text',...
+ 'BackgroundColor',COLORS(k,:),...
+ 'Units','Normalized','Position',[cx(k),my(3),cw,mh/2]);
+ label_pointer{k} = uicontrol('String',label_char(k),'Style','text',...
+ 'HorizontalAlignment','center',...
+ 'Units','Normalized','Position',[cx(k),my(3)+mh/2,cw,mh/2]);
+end
+
+uicontrol('String','Thickness [T] : ','Style','text','HorizontalAlignment','right',...
+ 'Units','Normalized','Position',[0,my(2),0.25,mh]);
+thickness_text = uicontrol('String',thickness_options_string(THICKNESS_INDEX),...
+ 'Style','text','HorizontalAlignment','left',...
+ 'Units','Normalized','Position',[0.25,my(2),0.25,mh]);
+
+uicontrol('String','Show [W] :','Style','text','HorizontalAlignment','right',...
+ 'Units','Normalized','Position',[0,my(1),0.25,mh]);
+show_text = uicontrol('String',show_options{SHOW_INDEX},...
+ 'Style','text','HorizontalAlignment','left',...
+ 'Units','Normalized','Position',[0.25,my(1),0.25,mh]);
+
+drawnow % pointer shows busy system
+
+LIMITS = [1,c-0.5,1,r+0.5]; % to capture zoom
+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);
+
+compute_overlays
+
+% ready to draw
+set(fig,'Pointer','arrow','WindowButtonDownFcn',@start_draw,...
+ 'WindowButtonMotionFcn',@pointer_motion);
+XO = []; % current drawing point
+uiwait % waits with assigning output until a figure is closed
+
+%%%%%%%%%% CALLBACK FUNCTIONS %%%%%%%%%%
+ function key_press(~,object)
+ % keyboard commands
+ key = object.Key;
+ numkey = str2double(key);
+ if ~isempty(numkey) && numkey<=nr_labels;
+ label_char(LABEL+1)=' ';
+ LABEL = numkey;
+ set(labels_text,'String',num2str(LABEL),...
+ 'BackgroundColor',COLORS(LABEL+1,:));
+ label_char(LABEL+1)='|';
+ for kk = 1:nr_labels+1
+ set(label_pointer{kk},'String',label_char(kk));
+ end
+ else
+ switch key
+ case 'l'
+ label_char(LABEL+1)=' ';
+ LABEL = move_once(LABEL+1,nr_labels+1,...
+ any(strcmp(object.Modifier,'shift')))-1;
+ set(labels_text,'String',num2str(LABEL),...
+ 'BackgroundColor',COLORS(LABEL+1,:));
+ label_char(LABEL+1)='|';
+ for kk = 1:nr_labels+1
+ set(label_pointer{kk},'String',label_char(kk));
+ end
+ case 'uparrow'
+ THICKNESS_INDEX = move_once(THICKNESS_INDEX,...
+ numel(thickness_options),false);
+ RADIUS = thickness2radius(thickness_options(THICKNESS_INDEX));
+ set(thickness_text,'String',...
+ thickness_options_string(THICKNESS_INDEX));
+ show_overlays(get_pixel);
+ case 'downarrow'
+ THICKNESS_INDEX = move_once(THICKNESS_INDEX,...
+ numel(thickness_options),true);
+ RADIUS = thickness2radius(thickness_options(THICKNESS_INDEX));
+ set(thickness_text,'String',...
+ thickness_options_string(THICKNESS_INDEX));
+ show_overlays(get_pixel);
+ case 't'
+ THICKNESS_INDEX = move_once(THICKNESS_INDEX,...
+ numel(thickness_options),any(strcmp(object.Modifier,'shift')));
+ RADIUS = thickness2radius(thickness_options(THICKNESS_INDEX));
+ set(thickness_text,'String',...
+ thickness_options_string(THICKNESS_INDEX));
+ show_overlays(get_pixel);
+ case 'w'
+ SHOW_INDEX = move_once(SHOW_INDEX,numel(show_options),...
+ any(strcmp(object.Modifier,'shift')));
+ set(show_text,'String',show_options{SHOW_INDEX})
+ compute_overlays
+ show_overlays(get_pixel);
+ case 'c'
+ COLORMAP_INDEX = move_once(COLORMAP_INDEX,...
+ length(colormap_options),...
+ any(strcmp(object.Modifier,'shift')));
+ COLORS = colormap_options{COLORMAP_INDEX}(nr_labels);
+ set(labels_text,'BackgroundColor',COLORS(LABEL+1,:));
+ for kk = 1:nr_labels+1
+ set(label_colorcubes{kk},'BackgroundColor',COLORS(kk,:))
+ end
+ compute_overlays
+ show_overlays(get_pixel);
+ case 's'
+ save_project
+ case 'e'
+ export_DC
+ case 'q'
+ close(fig)
+ end
+ end
+ end
+
+ function adjust_limits(~,~)
+ % response to zooming and panning
+ LIMITS([1,2]) = get(labeling_axes,'XLim');
+ LIMITS([3,4]) = get(labeling_axes,'YLim');
+ 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
+
+ function start_draw(~,~)
+ % click in the image
+ x = get_pixel;
+ if is_inside(x)
+ if RADIUS>0 % thickness>0
+ XO = x;
+ M = disc(XO,RADIUS,nr_circle_pts,[r,c]);
+ set(fig,'WindowButtonMotionFcn',@drag_and_draw,...
+ 'WindowButtonUpFcn',@end_draw)
+ else % fill
+ M = fill(x);
+ end
+ update(M);
+ end
+ end
+
+ function drag_and_draw(~,~)
+ % drag after clicking in the image
+ x = get_pixel;
+ M = stadium(XO,x,RADIUS,nr_circle_pts,[r,c]);
+ update(M);
+ XO = x;
+ end
+
+ function end_draw(~,~)
+ % release click after clicking in the image
+ M = stadium(XO,get_pixel,RADIUS,nr_circle_pts,[r,c]);
+ update(M);
+ XO = [];
+ set(fig,'WindowButtonMotionFcn',@pointer_motion,...
+ 'WindowButtonUpFcn','')
+ end
+
+ function pointer_motion(~,~)
+ % move around without clicking
+ if strcmp(zoom_handle.Enable,'off') && ...
+ strcmp(pan_handle.Enable,'off') % not zooming or panning
+ x = get_pixel;
+ if is_inside(x)
+ set(fig,'Pointer','crosshair')
+ else
+ set(fig,'Pointer','arrow')
+ end
+ show_overlays(x);
+ end
+ end
+
+%%%%%%%%%% HELPING FUNCTIONS %%%%%%%%%%
+
+ function save_project
+ % save mat file with user input and setting
+ % save images as separate files
+ [file,path] = uiputfile('correction.mat','Save project as');
+ if ~isequal(file,0) && ~isequal(path,0)
+ matfile = fullfile(path,file);
+ roothname = matfile(1:find(matfile=='.',1,'last')-1);
+ current_settings.show = show_options{SHOW_INDEX};
+ save(matfile,'DRAWINGS','CORRECTED',...
+ 'dictopt','nr_labels','current_settings')
+ imwrite(DRAWING_OVERLAY,[roothname,'_drawing.png'])
+ imwrite(CORRECTED_OVERLAY,[roothname,'_corrected.png'])
+ end
+ end
+
+ function export_DC
+ button = questdlg({'Exporting variables to the base workspace',...
+ 'might overwrite existing variables'},...
+ 'Exporting variables','OK','Cancel','OK');
+ if strcmp(button,'OK')
+ assignin('base','gui_D',DRAWING)
+ assignin('base','gui_C',CORRECTED)
+ end
+
+ end
+
+ function a = is_inside(x)
+ % check if x is inside image limits
+ a = inpolygon(x(1),x(2),LIMITS([1,2,2,1]),LIMITS([3,3,4,4]));
+ end
+
+ function p = get_pixel
+ % get cursor position
+ p = get(labeling_axes,'CurrentPoint');
+ p = round(p(1,[1,2]));
+ end
+
+ function show_overlays(x)
+ % overlay a circular region where the pointer and show
+ shown_left = DRAWING_OVERLAY;
+ shown_right = CORRECTED_OVERLAY;
+ if RADIUS>0 % thickness>0
+ P = repmat(disc(x,RADIUS,nr_circle_pts,[r,c]),[1,1,3]);
+ shown_left(P(:)) = 0.5+0.5*DRAWING_OVERLAY(P(:));
+ shown_right(P(:)) = 0.5+0.5*CORRECTED_OVERLAY(P(:));
+ end
+ % we have to imagesc(shown) to remove overlay if needed
+ axes(labeling_axes), cla, imagesc(shown_left)
+ axes(segmentation_axes), cla, imagesc(shown_right)
+ end
+
+ function update(M)
+ % change the state of the segmentation by updating LABELINGS with a
+ % mask M, and updating PROBABILITIES
+ DRAWING(M(:)) = LABEL;
+ if LABEL>0
+ CORRECTED(M(:)) = LABEL;
+ else
+ CORRECTED(M(:)) = segmentation(M(:));
+ end
+ compute_overlays % computing overlay images
+ show_overlays(get_pixel); % showing overlay and pointer
+ end
+
+ function compute_overlays
+ % computes overlays but not pointer overalay
+ uncorrected = repmat(DRAWING==0,[1 1 3]);
+ DRAWING_OVERLAY = uncorrected.*image + ~uncorrected.*...
+ (color_weight.*ind2rgb(DRAWING,COLORS) + (1-color_weight).*gray);
+ switch show_options{SHOW_INDEX}
+ case 'segmentation'
+ CORRECTED_OVERLAY = ind2rgb(CORRECTED,COLORS);
+ case 'overlay'
+ CORRECTED_OVERLAY = color_weight*...
+ ind2rgb(CORRECTED,COLORS) + (1-color_weight)*gray;
+ end
+ end
+
+% TODO disc shold be saved as a list of index shifts with respect to
+% the central pixel, and change only when thickness changes
+ function M = disc(x,r,N,dim)
+ % disc shaped mask in the image
+ angles = (0:2*pi/N:2*pi*(1-1/N));
+ X = x(1)+r*cos(angles);
+ Y = x(2)+r*sin(angles);
+ M = poly2mask(X,Y,dim(1),dim(2));
+ end
+
+ function M = stadium(x1,x2,r,N,dim)
+ % stadium shaped mask in the image
+ angles = (0:2*pi/N:pi)-atan2(x1(1)-x2(1),x1(2)-x2(2));
+ X = [x1(1)+r*cos(angles), x2(1)+r*cos(angles+pi)];
+ Y = [x1(2)+r*sin(angles), x2(2)+r*sin(angles+pi)];
+ M = poly2mask(X,Y,dim(1),dim(2));
+ end
+
+ function M = fill(x)
+ M = bwselect(DRAWING==DRAWING(x(2),x(1)),x(1),x(2),4);
+ end
+
+ function [I,G] = normalize_image(I)
+ % initialization: normalize image
+ if isa(I,'uint8')
+ I = double(I)/255;
+ end
+ if isa(I,'uint16')
+ I = double(I)/65535;
+ end
+ if size(I,3)==3 % rgb image
+ G = repmat(rgb2gray(I),[1 1 3]);
+ else % assuming grayscale image
+ I = repmat(I,[1,1,3]);
+ G = I;
+ end
+ end
+
+ function n = move_once(n,total,reverse)
+ % moves option index once, respecting total number of options
+ if ~reverse
+ n = mod(n,total)+1;
+ else
+ n = mod(n+total-2,total)+1;
+ end
+ end
+
+ function r = thickness2radius(t)
+ r = t/2+0.4;
+ end
+
+ function clean_toolbar
+ set(fig,'MenuBar','none','Toolbar','figure');
+ all_tools = allchild(findall(fig,'Type','uitoolbar'));
+ %my_tool = uipushtool(toolbar,'CData',rand(16,16,3),'Separator','on',...
+ % 'TooltipString','my tool','Tag','my tool');
+ for i=1:numel(all_tools)
+ if isempty(strfind(all_tools(i).Tag,'Pan'))&&...
+ isempty(strfind(all_tools(i).Tag,'Zoom'))&&...
+ isempty(strfind(all_tools(i).Tag,'SaveFigure'))&&...
+ isempty(strfind(all_tools(i).Tag,'PrintFigure'))&&...
+ isempty(strfind(all_tools(i).Tag,'DataCursor'))
+ delete(all_tools(i)) % keeping only Pan, Zoom, Save and Print
+ end
+ end
+ end
+
+end
diff --git a/code/texture_gui/code/functions/update_dictionary.m b/code/texture_gui/code/functions/update_dictionary.m
new file mode 100755
index 0000000000000000000000000000000000000000..11627b2d4d2026bcde2b7f23fb82a432405c7c48
--- /dev/null
+++ b/code/texture_gui/code/functions/update_dictionary.m
@@ -0,0 +1,9 @@
+function dictionary = update_dictionary(dictionary,gui_dictprob)
+
+M = dictionary.options.patch_size;
+nr_dict_patches = size(dictionary.tree,2);
+L = size(gui_dictprob,2); % nr labels
+gui_dictprob = reshape(gui_dictprob,[M^2,nr_dict_patches,L]);
+gui_dictprob = permute(gui_dictprob,[1,3,2]);
+gui_dictprob = reshape(gui_dictprob,[M^2*L,nr_dict_patches]);
+dictionary.dictprob = gui_dictprob;
\ No newline at end of file
diff --git a/code/texture_gui/code/reusing_labels_script.m b/code/texture_gui/code/reusing_labels_script.m
new file mode 100755
index 0000000000000000000000000000000000000000..f25eb922bd6ce18a3b51f9658c4596ed6e318df2
--- /dev/null
+++ b/code/texture_gui/code/reusing_labels_script.m
@@ -0,0 +1,23 @@
+clear
+close all
+addpath functions
+
+% example of re-using externaly made labeling image
+im = imread('../data/randen15B.png'); % randen 5 textures
+labeling = double(imread('../data/randen_labels_rgb.png'));
+dictopt.patch_size = 15;
+dictopt.branching_factor = 4;
+dictopt.number_layers = 4;
+dictopt.number_training_patches = 2000;
+dictopt.normalization = false;
+dictopt.method = 'euclidean';
+image_texture_gui(im,dictopt,5,labeling)
+
+%% an image exported from gui can be re-used
+image_texture_gui(im,dictopt,5) %% <- EXPORT (E) labeling here
+image_texture_gui(im,dictopt,5,gui_L) %% <- use labeling here
+
+%% an image saved from gui can be re-used
+image_texture_gui(im,dictopt,5) %% <- SAVE (S) labeling here using some filename
+image_texture_gui(im,dictopt,5,imread('filename_labels_indexed.png')) %% <- use labeling here
+
diff --git a/code/texture_gui/code/test.m b/code/texture_gui/code/test.m
new file mode 100755
index 0000000000000000000000000000000000000000..b05b7fdaef8a69b9cb8668204e2d70dbbfe29803
--- /dev/null
+++ b/code/texture_gui/code/test.m
@@ -0,0 +1,31 @@
+addpath('functions')
+%%
+fig = figure;
+imagesc(rand(200,200));
+set(fig,'MenuBar','none','Toolbar','figure');
+
+%%
+toolbar = findall(fig,'Type','uitoolbar');
+all_tools = allchild(toolbar);
+% removing tools
+for i=1:numel(all_tools)
+ t = get(all_tools(i),'Tag');
+ if isempty(strfind(t,'Pan'))&&...
+ isempty(strfind(t,'Zoom'))&&...
+ isempty(strfind(t,'SaveFigure'))&&...
+ isempty(strfind(t,'PrintFigure'))&&...
+ isempty(strfind(t,'DataCursor'))
+ delete(all_tools(i)) % keeping only Pan, Zoom, Save and Print
+ end
+end
+%% adding a tool
+[icon,~,alpha] = imread('linkaxesicon.png');
+icon = double(icon)/255;
+icon(alpha==0)=NaN;
+uitoggletool(toolbar,'CData',icon,...
+ 'TooltipString','Link Axes','Tag','LinkAxes',...
+ 'OnCallback',{@link_axes,'xy'},...
+ 'OffCallback',{@link_axes,'off'});
+%% changing the order of tools
+all_tools = allchild(toolbar);
+set(toolbar,'children',all_tools([2,1,3:end]));
\ No newline at end of file
diff --git a/code/texture_gui/code/texture_gui_uses_script.m b/code/texture_gui/code/texture_gui_uses_script.m
new file mode 100755
index 0000000000000000000000000000000000000000..aa2c156a988022be9cfc72a49b46b3728a8d7483
--- /dev/null
+++ b/code/texture_gui/code/texture_gui_uses_script.m
@@ -0,0 +1,29 @@
+clear
+close all
+
+addpath functions
+
+%% demo: no input
+image_texture_gui
+
+%% usage 1: only image (default dictionary options )
+im = double(imread('bag.png'))/255;
+image_texture_gui(im)
+
+%% usage 2: image and dictionary options
+dictopt.method = 'euclidean';
+dictopt.patch_size = 11;
+dictopt.branching_factor = 2;
+dictopt.number_layers = 5;
+dictopt.number_training_patches = 30000;
+dictopt.normalization = false;
+image_texture_gui(im,dictopt,5)
+
+%% usage 3: image and dictionary
+dictionary = build_dictionary(im,dictopt);
+image_texture_gui(im,dictionary,2)
+
+%% usage 4: image and mappings
+mappings = compute_mappings(im,dictionary);
+image_texture_gui(im,mappings,2)
+