text:
function saveaspdf(h,filename)
%SAVEASPDF Save a figure as a clean pdf file ready for publication.
%   saveaspdf saves the current figure to the path defined in the global
%   SAVEASPDF_PATH, or to your Desktop if no path is specified. The
%   figure's 'Name' property is used as the filename.
%
%   saveaspdf(filename) saves the current figure as filename to the path
%   defined in the global SAVEASPDF_PATH, or to your Desktop if no path is
%   specified.
%
%   saveaspdf(h,filename) saves the figure with handle h as filename to the
%   path defined in the global SAVEASPDF_PATH, or to your Desktop if no
%   path is specified.
%
%   saveaspdf defines several globals which affect its output:
%
%      SAVEASPDF_DPI = 300            - The resolution of the pdf.
%      SAVEASPDF_PATH = ~/Desktop     - Where the pdf is saved.
%      SAVEASPDF_SCALELINEWIDTH = 1.8 - All lines and fonts are scaled by
%      SAVEASPDF_SCALEFONTSIZE = 1.1    these factors before printing. The
%                                       figure is returned to its original
%                                       state after printing.
%      SAVEASPDF_SIZE = []            - Set the output dimensions to [width
%                                       height] in centimeters. If empty,
%                                       the figure's dimensions are used.
%      SAVEASPDF_TRANSPARENT = false  - Set to true change the background
%                                       to be transparent (experimental).
%
%   See also saveas, print.

%   Note: add the following lines to Matlab's ghostscript.m for better
%   quality prints of figures with gradients or transparency. Source:
%   www.histed.org/wiki/index.php?title=Controlling_matlab's_pdf_output.
%
%   % In between "fprintf(rsp_fid, '-sOutputFile="%s"\n', pj.GhostName );"
%   % and "fclose(rsp_fid);":
%   if strcmp(pj.GhostDriver, 'pdfwrite')
%       % Uncomment to embed images into file as ZIP not jpeg compressed
%       fprintf(rsp_fid, ['-dAutoFilterColorImages=false ' ...
%           '-dColorImageFilter=/FlateEncode\n']);
%       % control resolution/downsampling of embedded images
%       fprintf(rsp_fid, '-dDownsampleColorImages=false\n');
%       %fprintf(rsp_fid, '-dColorImageDownsampleThreshold=1.5\n');
%       %fprintf(rsp_fid, '-dColorImageResolution=300\n');
%   end

%   Version: 07/07/11
%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)

% Define globals.
global SAVEASPDF_DPI;
global SAVEASPDF_PATH;
global SAVEASPDF_SCALELINEWIDTH;
global SAVEASPDF_SCALEFONTSIZE;
global SAVEASPDF_SIZE;
global SAVEASPDF_TRANSPARENT;

% Set up default values.
if nargin == 1 && ischar(h)
    filename = h;
    h = gcf;
elseif nargin == 0
    h = gcf;
    filename = get(h,'Name');
    if isempty(filename)
        filename = ['Figure ' int2str(gcf)];
    end
end
if ~isempty(SAVEASPDF_DPI)
    dpi = SAVEASPDF_DPI;
else
    dpi = 300;
end
if ~isempty(SAVEASPDF_PATH)
    path = SAVEASPDF_PATH;
else
    if ispc, path = 'USERPROFILE'; else path = 'HOME'; end
    path = fullfile(getenv(path),'Desktop');
end
filename = fullfile(path,regexprep(filename,'\.[a-zA-Z]+','.pdf'));
if ~isempty(SAVEASPDF_SCALELINEWIDTH)
    scalelw = SAVEASPDF_SCALELINEWIDTH;
else
    scalelw = 1.8;
end
if ~isempty(SAVEASPDF_SCALEFONTSIZE)
    scalefs = SAVEASPDF_SCALEFONTSIZE;
else
    scalefs = 1.1;
end
if ~isempty(SAVEASPDF_SIZE)
    size = SAVEASPDF_SIZE;
else
    size = [];
end
if ~isempty(SAVEASPDF_TRANSPARENT)
    trans = SAVEASPDF_TRANSPARENT;
else
    trans = false;
end

% Backup current settings.
oldUnits = get(h,'Units');
oldPaperUnits = get(h,'PaperUnits');
oldPaperPosition = get(h,'PaperPosition');
oldPaperSize = get(h,'PaperSize');
axes = findall(h,'type','axes');
lines = [findall(h,'type','line'); axes];
text = [findall(h,'type','text'); axes];
oldFigColor = get(h,'Color');
oldAxesColor = get(axes,'Color');
oldLineWidth = get(lines,'LineWidth');
oldFontSize = get(text,'FontSize');
if ~iscell(oldAxesColor), oldAxesColor = {oldAxesColor}; end
if ~iscell(oldLineWidth), oldLineWidth = num2cell(oldLineWidth); end
if ~iscell(oldFontSize), oldFontSize = num2cell(oldFontSize); end

% Scale the line width and font size, reset units and update the positions.
set(lines,{'LineWidth'},num2cell(cellfun(@(x)x*scalelw,oldLineWidth)));
set(text,{'FontSize'},num2cell(cellfun(@(x)round(x*scalefs),oldFontSize)));
set(h,'Units','centimeters');
set(h,'PaperUnits','centimeters');
if isempty(size)
    curPosition = get(h,'Position');
    set(h,'PaperPosition',[0,0,curPosition(3:4)]);
    set(h,'PaperSize',curPosition(3:4));
else
    set(h,'PaperPosition',[0,0,size]);
    set(h,'PaperSize',size);
end
if trans
    set(h,'Color','none');
    set(axes,'Color','none');
end

% Print the pdf.
print(h,'-dpdf',sprintf('-r%d',dpi),filename);

% Restore the previous settings.
if trans
    set(h,'Color',oldFigColor);
    set(axes,{'Color'},oldAxesColor);
end
set(h,'Units',oldUnits);
set(h,'PaperUnits',oldPaperUnits);
set(h,'PaperPosition',oldPaperPosition);
set(h,'PaperSize',oldPaperSize);
set(text,{'FontSize'},oldFontSize);
set(lines,{'LineWidth'},oldLineWidth);

author: Sable

repo: mcbench-benchmarks


text:
function [suspicious_index lof] = LOF(A, k)
%#####################################################################
%# Local Outlier Factor                                              #
%# Authors: Markus M. Breunig, Hans-Peter Kriegel,                   #
%#          Raymond T. Ng, J?rg Sander                               #
%# Oringinal paper :                                                 #
%# LOF: Identifying Density-Based Local Outliers                     #
%# e-mail : { breunig | kriegel | sander }                           #
%#          @dbs.informatik.uni-muenchen.de                          #
%#          rng@cs.ubc.ca                                            #
%# Programmer: Yi-Ren Yeh(yirenyeh@gmail.com)                        #
%# modified by: Zi-Wen Gui(evan176@hotmail.com)                      #
%#                                                                   #
%#                                                                   #
%# Inputs                                                            #
%#   A: the data matrix, each row reprent an instance                #
%#   k: the number of nearest neighbors                              #
%#                                                                   #
%# Outputs                                                           #
%#   lof: the local outlier factor for each instance                 #
%#   suspicious_index: the ranking of instances according to their   #
%#                     suspiciuous score                             #
%#                     For example, suspicious_index(i)=j means the  #
%#                     ith instance is in jth position in the ranking#
%#####################################################################

%Find the nearest neighbors by "KDTree" for each elements
[k_index, k_dist] = knnsearch(A,A,'k',k+1,'nsmethod','kdtree');
%Ignore first element(itself) at nearest neighbors
k_index = k_index(:,2:end);
%Get k-distance
k_dist1 = k_dist(:,end);
%Get row length of matrix A
n = length(A(:,1));
%Initializa lrd_value vector
lrd_value = zeros(n,1);
%Calculate lrd for each elements
for i = 1:n
    lrd_value(i) = lrd(A, i, k_dist1,k_index, k);
end
%Initializa lof vector
lof = zeros(n,1);
%Calculate LOF
for i = 1:n
    lof(i) = sum(lrd_value(k_index(i,:))/lrd_value(i))/k;
end
%Sort lof factor
[non,suspicious_index]=sort(lof,'descend');

%=========================================================================
function lrd_value = lrd(A, index_p, k_dist,k_index, k)
%Calculate the reachability distance for nearest neighbors
Temp = repmat(A(index_p,:),k,1) - A(k_index(index_p,:),:);
Temp = sqrt(sum(Temp.^2,2));
%max{k-distance(b), d(a, b)}
rd_dsit = max([Temp k_dist(k_index(index_p,:))],[],2);
%Calculate the local reachability density for each elements
lrd_value = k/sum(rd_dsit);

author: dsmi-lab-ntust

repo: AnomalyDetectionToolbox


text:
function [Gamma,scale] = id_Gamma_inference(L,Pi,order)
% inference for independent samples (ignoring time structure)
% scale is the likelihood

L(L<realmin) = realmin;
Gamma = repmat(Pi,size(L,1),1) .* L;
Gamma = Gamma(1+order:end,:);
scale = sum(Gamma,2);
Gamma = rdiv(Gamma,scale);

end

author: OHBA-analysis

repo: HMM-MAR


text:
% batch_export_dvh_to_excel.m
%
% This script exports DVHs for the specified structure names to Excel.
%
% The binwidth to use for DVH calculation is defined in CERROptions.json.
%
% APA, 1/28/2019

% Directory containing CERR files
dirPath = 'path:\to\cerr\files';

% Output Excel file
xlSaveFile = 'path:\to\expoty\excel\dvh.xlsx';

% Names of structures to calculate and export DVHs
strC = {'Lungs_combined', 'Lung_L_combined', 'Lung_R_combined',...
    'Heart_combined', 'Outer_combined', 'Stomach_combined', 'PTV_combined'};

%Find all CERR files
fileC = {};
if strcmpi(dirPath,'\') || strcmpi(dirPath,'/')
    filesTmp = getCERRfiles(dirPath(1:end-1));
else
    filesTmp = getCERRfiles(dirPath);
end
fileC = [fileC filesTmp];

filesNotConvertedC = {};

pathStr = getCERRPath;
optName = [pathStr 'CERROptions.json'];
optS = opts4Exe(optName);

dvhC = {};

%Loop over CERR plans
for iFile = 1:length(fileC)

    planC = loadPlanC(fileC{iFile},tempdir);
    planC = quality_assure_planC(fileC{iFile},planC);
    planC = updatePlanFields(planC);
    indexS = planC{end};

    strNameC = {planC{indexS.structures}.structureName};

    for indStr = 1:length(strC)
                %DVH
                strNum = find(strcmpi(strC{indStr},strNameC));
                doseNum = length(planC{indexS.dose});
                if ~isempty(strNum)
                    [dosesV, volsV] = getDVH(strNum, doseNum, planC);
                    [doseBinsV, volsHistV] = doseHist(dosesV, volsV, optS.DVHBinWidth);
                    doseBinsV = doseBinsV * 100; % Gy to cGy
                    fVol = volsHistV;
                    dvhC{iFile,indStr} = [doseBinsV; fVol];
                else
                    dvhC{iFile,indStr} = [];
                end
    end

end

for iDvh = 1:size(dvhC,2)
    dvhM = [dvhC{:,iDvh}];
    doseHistV = unique(dvhM(1,:));
    dvhM = zeros(length(fileC),length(doseHistV));
    doseHistC{iDvh} = doseHistV;
    for iFile = 1:size(dvhC,1)
        if ~isempty(dvhC{iFile,iDvh})
            indMin = findnearest(doseHistV,min(dvhC{iFile,iDvh}(1,:)));
            indMax = findnearest(doseHistV,max(dvhC{iFile,iDvh}(1,:)));
            dvhM(iFile,indMin:indMax) = dvhC{iFile,iDvh}(2,:);
        end
    end
    dvhToWriteC{iDvh} = dvhM;
end

[~,namC] = cellfun(@fileparts,fileC,'UniformOutput',false);
for indStr = 1:length(strC)
    xlswrite(xlSaveFile,namC',strC{indStr},['A2:A',num2str(length(fileC)+1)])
    colEnd = xlsColNum2Str(length(doseHistC{indStr})+1);
    xlswrite(xlSaveFile,doseHistC{indStr},strC{indStr},...
        ['B1:',colEnd{1},'1'])
    colEnd = xlsColNum2Str(size(dvhToWriteC{indStr},2)+1);
    xlswrite(xlSaveFile,dvhToWriteC{indStr},strC{indStr},...
        ['B2:',colEnd{1},num2str(size(dvhToWriteC{indStr},1)+1)])
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%GRIDBDRX Reconstruct interior boundaries.
%
%   [ GRID ] = GRIDBDRX( GRID, ORDER ) Reconstructs interior
%   boundaries and updates the GRID.B field. ORDER (default true)
%   changes the order of boundaries so external boundaries are first,
%   and internal last.
%
%   See also GRIDBDR

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function results_image = PTKGetAirwayImageFromCentreline(label_bronchi, airway_root, template, colour_by_segment_index)
    % PTKGetAirwayImageFromCentreline. Generates an image of the airways
    %     coloured according to a number of subtrees. The start branch of the
    %     subtrees are specified in label_bronchi. Any bronchi serving more than
    %     one subtree is coloured grey.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    results_image = template.BlankCopy;

    % Create a vector containing every centreline voxel in the tree, and a cell
    % array, containing a vector of centreline points for each label segment
    all_centreline_voxels = [];
    segment_centreline_voxels_extended = [];
    for label_bronchus_index = 1 : length(label_bronchi)
        centreline_voxels = PTKTreeUtilities.GetCentrelineVoxelsForTheseBranchesExtended(label_bronchi(label_bronchus_index), template);
        segment_centreline_voxels_extended{label_bronchus_index} = centreline_voxels;
        all_centreline_voxels = [all_centreline_voxels; centreline_voxels];
    end
    all_centreline_voxels = unique(all_centreline_voxels);

    results_image_raw = zeros(template.ImageSize, 'uint8');

    bronchi_to_do = airway_root;

    % Find the voxel index of the endpoint of each label bronchus. We want to ensure
    % that the bronchus
    endpoint_voxels = [];
    for label_bronchus_index = 1 : length(label_bronchi)
        end_centreline_point = PTKTreeUtilities.GetEndmostPoint(label_bronchi(label_bronchus_index));
        endpoint_voxels(label_bronchus_index) = PTKTreeUtilities.CentrelinePointsToLocalIndices(end_centreline_point, template);
    end

    while ~isempty(bronchi_to_do)
        airway_bronchus = bronchi_to_do(end);
        bronchi_to_do(end) = [];
        voxel_indices_in_segment = airway_bronchus.GetAllAirwayPoints;
        voxel_indices_in_segment = template.GlobalToLocalIndices(voxel_indices_in_segment);

        label_bronchus_indices = [];

        for label_bronchus_index = 1 : length(label_bronchi)
            if ismember(endpoint_voxels(label_bronchus_index), voxel_indices_in_segment)
                label_bronchus_indices = label_bronchus_index;
            end
        end

        if isempty(label_bronchus_indices)
            number_of_centreline_points_in_bronchus = sum(ismember(all_centreline_voxels, voxel_indices_in_segment));

            % Search through all segments to see if the voxels match a bronchus
            % Only match bronchi if more than 50% of the centreline voxels in the
            % bronchus are part of this bronchus
            for label_bronchus_index = 1 : length(label_bronchi)
                number_of_centreline_points_in_bronchus_for_this_centreline = sum(ismember(segment_centreline_voxels_extended{label_bronchus_index}, voxel_indices_in_segment));
                if number_of_centreline_points_in_bronchus_for_this_centreline > (number_of_centreline_points_in_bronchus/2)
                    label_bronchus_indices = [label_bronchus_indices, label_bronchus_index];
                end
            end
        end
        % Colour the bronchus if it matches only one segment; otherwise colour it grey
        if length(label_bronchus_indices) == 1
            if colour_by_segment_index
                colour = label_bronchus_indices;
            else
                colour = 1 + mod(label_bronchus_indices - 1, 6);
            end
        else
            if colour_by_segment_index
                colour = 0;
            else
                colour = 7;
            end
        end
        results_image_raw(voxel_indices_in_segment) = colour;

        bronchi_to_do = [bronchi_to_do airway_bronchus.Children];
    end

    results_image.ChangeRawImage(results_image_raw);
    results_image.ImageType = PTKImageType.Colormap;
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [varargout]=quiver3Dpatch(x,y,z,ux,uy,uz,c,a)

% function [F,V,C]=quiver3Dpatch(x,y,z,ux,uy,uz,c,a)
% ------------------------------------------------------------------------
%
% This function allows plotting of colored 3D arrows by generating patch
% data (faces F, vertices V and color data C). The patch data which
% allows plotting of 3D quiver arrows with specified (e.g. colormap driven)
% color. To save memory n arrows are created using only n*6 faces and n*7
% vertices. The vector "a" defines arrow length scaling where a(1) is the
% smallest arrow length and a(2) the largest. Use the PATCH command to plot
% the arrows:
%
% [F,V,C]=quiver3Dpatch(x,y,z,ux,uy,uz,a)
% patch('Faces',F,'Vertices',V,'CData',C,'FaceColor','flat');
%
% Below is an example illustrating color specifications for (combined)
% patch data.
%
%%% EXAMPLE
% % Simulating 3D volume and vector data
% n=25;
% [X,Y,Z]=meshgrid(linspace(-4.77,4.77,n));
% phi=(1+sqrt(5))/2;
% M=2 - (cos(X + phi*Y) + cos(X - phi*Y) + cos(Y + phi*Z) + cos(Y - phi*Z) + cos(Z - phi*X) + cos(Z + phi*X));
%
% % Simulating vector data
% %Vector data here based on the gradient of the image
% [u,v,w] = gradient(M);
% G=hypot(hypot(u,v),w); %Vector lenghts
%
% a=[min(G(:)) max(G(:))]; %Arrow length scaling
% L=G>0.7; %Logic indices for arrows to make a selection if required
% [Fv,Vv,Cv]=quiver3Dpatch(X(L),Y(L),Z(L),u(L),v(L),w(L),G(L),a);
%
% figure;
% xlabel('X','FontSize',10);ylabel('Y','FontSize',10);zlabel('Z','FontSize',10);
% title('Colormap driven vectors','FontSize',15);
% patch('Faces',Fv,'Vertices',Vv,'EdgeColor','none', 'CData',Cv,'FaceColor','flat','FaceAlpha',1);
%
% colormap jet; colorbar;
% % caxis([min(Cv(:)) max(Cv(:))]);
% view(3); grid on; axis equal; axis vis3d;
% set(gca,'FontSize',10);
% camlight headlight; lighting flat;
% drawnow;
%
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2014/01/13 Updated example in help
% 2014/11/11 Updated to allow for RGB color input data (e.g. C is an nx3 array)
% 2016/06/30 Added handling of empty "a" parameter. Reverts to default
% normal length range.
%------------------------------------------------------------------------

%%

%Convert to columns
x=x(:); y=y(:); z=z(:);
ux=ux(:); uy=uy(:); uz=uz(:);

%Spherical coordinates
[THETA_vec,PHI_vec,R_vec] = cart2sph(ux,uy,uz);

%% Setting arrow size properties

if isempty(a)
    a=[min(R_vec(:)) max(R_vec(:))];
end

%Arrow length
if  min(R_vec(:))==max(R_vec(:)) %If all radii are equal, or if just 1 vector is used
    arrow_length=a(2)*ones(size(R_vec)); %All arrow lengths become a(2)
else
    %Scale arrow lengths between a(1) and a(2)
    arrow_length=R_vec-min(R_vec(:));
    arrow_length=a(1)+((arrow_length./max(arrow_length(:))).*(a(2)-a(1)));
end

%Other arrow dimensions as functions of arrow length and phi (golden ratio)
phi=(1+sqrt(5))/2;
head_size=arrow_length./(phi);
head_width=head_size./(2.*phi);
stick_width=head_width./(2.*phi);
if numel(a)==3
    head_width=a(3)*head_width;
    stick_width=a(3)*stick_width;
end
h=sin((2/3).*pi).*stick_width;
ha=sin((2/3).*pi).*head_width;

%% Creating arrow triangle vertices coordinates

X_tri=[zeros(size(x))  zeros(size(x)) zeros(size(x))...
    head_size.*ones(size(x)) head_size.*ones(size(x)) head_size.*ones(size(x))...
    arrow_length];
Y_tri=[-(0.5.*stick_width).*ones(size(x)) (0.5.*stick_width).*ones(size(x))  zeros(size(x))...
    -(0.5.*head_width).*ones(size(x))  (0.5.*head_width).*ones(size(x))  zeros(size(x))...
    zeros(size(x))];
Z_tri=[-(0.5.*stick_width.*tan(pi/6)).*ones(size(x))...
    -(0.5.*stick_width.*tan(pi/6)).*ones(size(x))...
    (h-(0.5.*stick_width.*tan(pi/6))).*ones(size(x))...
    -(0.5.*head_width.*tan(pi/6)).*ones(size(x))...
    -(0.5.*head_width.*tan(pi/6)).*ones(size(x))...
    (ha-(0.5.*head_width.*tan(pi/6))).*ones(size(x))...
    zeros(size(x))];

% Rotating vertices
[THETA_ar,PHI_ar,R_vec_ar] = cart2sph(X_tri,zeros(size(Y_tri)),Z_tri);
PHI_ar=PHI_ar+PHI_vec*ones(1,size(THETA_ar,2));
[X_arg,Y_arg,Z_arg] = sph2cart(THETA_ar,PHI_ar,R_vec_ar);
Y_arg=Y_arg+Y_tri;
[THETA_ar,PHI_ar,R_vec_ar] = cart2sph(X_arg,Y_arg,Z_arg);
THETA_ar=THETA_ar+THETA_vec*ones(1,size(THETA_ar,2));
[X_arg,Y_arg,Z_arg] = sph2cart(THETA_ar,PHI_ar,R_vec_ar);

X_arg=X_arg+x*ones(1,size(THETA_ar,2)); X_arg=X_arg';
Y_arg=Y_arg+y*ones(1,size(THETA_ar,2)); Y_arg=Y_arg';
Z_arg=Z_arg+z*ones(1,size(THETA_ar,2)); Z_arg=Z_arg';

V=[X_arg(:) Y_arg(:) Z_arg(:)];

%% Creating faces matrix

%Standard vertex order for 6 face arrow style
F_order=[1 2 7; 2 3 7; 3 1 7; 4 5 7; 5 6 7; 6 4 7;];
no_nodes=size(X_tri,2);
b=(no_nodes.*((1:1:numel(x))'-1))*ones(1,3);

%Loops size(F_order,1) times
F=zeros(numel(x)*size(F_order,1),3); %Allocating faces matrix
for f=1:1:size(F_order,1)
    Fi=ones(size(x))*F_order(f,:)+b;
    F(1+(f-1)*numel(x):f*numel(x),:)=Fi;
end

%     %Alternative without for loop, not faster for some tested problems
%     F_order=(ones(numel(x),1)*[1 2 7 2 3 7 3 1 7 4 5 7 5 6 7 6 4 7])';
%     F_order1=ones(numel(x),1)*(1:6);
%     F_order2=ones(numel(x),1)*[2 3 1 5 6 4];
%     F_order=[F_order1(:) F_order2(:)];
%     F_order(:,3)=7;
%     b=repmat(((no_nodes.*(0:1:numel(x)-1)')*ones(1,3)),[6,1]);
%     F=F_order+b;

%% Collect face and vertex output
varargout{1}=F;
varargout{2}=V;

%% Add color specification if requested

if nargout==3
    if isempty(c) %If empty specify vector magnitudes as colormap driven color
        C=repmat(R_vec,[size(F_order,1),1]);
    else %If user specified color replicate to match # of added faces for arrow
        %c my be an nx3 array to allow for RGB type color data
        C=repmat(c,[size(F_order,1),1]);
    end
    varargout{3}=C;
end

end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [g]=gsp_design_expwin(G, bmax, a)
%GSP_DESIGN_EXPWIN create an expwin window of lenth N with parameter a
%   Usage: g = gsp_design_expwin(G);
%          g = gsp_design_expwin(G, bmax);
%          g = gsp_design_expwin(G, bmax, a);
%
%   Input parameters:
%       G       : Graph structure
%       bmax    : Maximum relative band (default 0.2)
%       a       : Slope parameter (default 1).
%
%   Output parameters
%       g       : filter
%
%   This function design the following filter:
%
%   .. g(x) = s( (1-x) / bmax /lmax )
%
%   .. math:: g(x) = s\left(\frac{1-x}{\lambda_{\text{max}} b_{\text{max}} }\right)
%
%   where  $s(x)$ is the step function
%
%   ..          /   0                                      if x < -1
%   ..  s(x) = | exp(-a/x) / ( exp(-a/x) + exp(-a/(1-x)) ) if x in [-1, 1]
%   ..          \   1                                      if x > 1
%
%   .. math:: s(x)=\begin{cases} 0 & \mbox{if }x<-1 \\ \frac{e^{-\frac{a}{x}}}{e^{-\frac{a}{x}}+e^{-\frac{a}{1-x}}} & \mbox{if }x\in[-1,1]\\ 1 & \mbox{if }x>1 \end{cases}
%
%   It uses a clever exponential construction to obtain an infinitely
%   differentiable function that is band limited!
%
%   This function will compute the maximum eigenvalue of the laplacian. To
%   be more efficient, you can precompute it using::
%
%       G = gsp_estimate_lmax(G);
%
%   Example:::
%
%         Nf = 4;
%         G = gsp_sensor(100);
%         G = gsp_estimate_lmax(G);
%         g = gsp_design_expwin(G);
%         gsp_plot_filter(G,g);
%

% Author: Nathanael Perraudin
% Date  : 18 December 2014

if nargin<3
    a = 1;
end

if nargin<2
    bmax = 0.2;
end

if numel(bmax)>1
    g = cell(numel(bmax),1);
    for ii = 1:numel(bmax);
        g{ii} = gsp_design_expwin(G,bmax(ii),a);
    end
    return
end

if isstruct(G)
    if ~isfield(G,'lmax')
%         if param.verbose
            fprintf('GSP_DESIGN_EXPWIN has to compute lmax \n')
%         end
        G = gsp_estimate_lmax(G);
    end
   lmax = G.lmax;
else
   lmax = G;
end

g = @(x) gsp_smooth_downstep(x/bmax/lmax, a , 1);

end

repo: gspbox


text:
function finetuned_model_path = train_LocNet_model(...
        image_db_train, image_db_val, model, conf)
% train_LocNet_model: trains the LocNet network
%
% INPUTS:
% 1) image_db_train: struct with the specifiers of the training dataset
% 2) image_db_val: struct with the specifiers of the testing dataset
% 3) model: struct with the localization model
% 4) conf: struct with the configuration options of the training process
%
% OUTPUT:
% 1) finetuned_model_path: string with the path to the trained caffe model
% produced from the training procedure.
%
% This file is part of the code that implements the following paper:
% Title      : "LocNet: Improving Localization Accuracy for Object Detection"
% Authors    : Spyros Gidaris, Nikos Komodakis
% Institution: Universite Paris Est, Ecole des Ponts ParisTech
% ArXiv link : http://arxiv.org/abs/1511.07763
% code       : https://github.com/gidariss/LocNet
%
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2016 Spyros Gidaris
%
% Title     : "LocNet: Improving Localization Accuracy for Object Detection"
% ArXiv link: http://arxiv.org/abs/1511.07763
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------

mkdir_if_missing(conf.finetune_rst_dir);
timestamp = datestr(datevec(now()), 'yyyymmdd_HHMMSS');
log_file  = fullfile(conf.finetune_rst_dir, 'output', ['LocNet_' timestamp, '.txt']);
mkdir_if_missing(fileparts(log_file));
finetuned_model_path = '';
diary(log_file);

image_db_train_file = fullfile(conf.finetune_rst_dir,'train_file.mat');
image_db_val_file   = fullfile(conf.finetune_rst_dir,'val_file.mat');

% load training dataset
try
    load(image_db_train_file,'image_db_train');
catch
    image_db_train = load_image_dataset(...
        'dataset',          image_db_train.dataset, ...
        'image_set',        image_db_train.image_set, ...
        'voc_year',         image_db_train.voc_year, ...
        'proposals_method', image_db_train.proposals_method,...
        'use_flips',        image_db_train.use_flips);
    image_db_train.cache_regions_dir = ...
        fullfile(pwd,'cached_regions', ...
        ['LocNet_model_', image_db_train.image_set_name, '_', image_db_train.proposals_suffix,filesep]);
    mkdir_if_missing(image_db_train.cache_regions_dir);

    fprintf('Prepare training data:\n')
    image_db_train = setup_image_db(image_db_train, conf);
    fprintf('Done\n')
    save(image_db_train_file,'image_db_train','-v7.3');
end

% load testing dataset
try
    load(image_db_val_file,'image_db_val');
catch
    image_db_val = load_image_dataset(...
        'dataset',          image_db_val.dataset, ...
        'image_set',        image_db_val.image_set, ...
        'voc_year',         image_db_val.voc_year, ...
        'proposals_method', image_db_val.proposals_method,...
        'use_flips',        image_db_val.use_flips);

    image_db_val.cache_regions_dir = ...
        fullfile(pwd,'cached_regions', ...
        ['LocNet_model_', image_db_val.image_set_name, '_', image_db_val.proposals_suffix,filesep]);
    mkdir_if_missing(image_db_val.cache_regions_dir);

    fprintf('Prepare validation data:\n')
    image_db_val = setup_image_db(image_db_val, conf);
    fprintf('Done\n')
    save(image_db_val_file,'image_db_val','-v7.3');
end

fprintf('conf:\n');
disp(conf);
fprintf('model:\n');
disp(model);
fprintf('pooler:\n');
disp(model.pooler_loc);
fprintf('model.loc_params:\n');
disp(model.loc_params);

current_dir = pwd;
cd(fileparts(conf.finetune_net_def_file))
% initialize caffe solver
[solver, iter_] = InitializeSolver(conf);

fprintf('Starting iteration %d\n', iter_);

prev_rng = seed_rand();
last_finetuned_model_prev_path = '';

rng(prev_rng);

iter_ = solver.iter();
mean_error_test = TestNet(conf, model, solver, image_db_val);
ShowState(iter_, NaN, mean_error_test);
if conf.test_only
    mean_error_test = TestNet(conf, model, solver, image_db_val);
    ShowState(iter_, NaN, mean_error_test);
    return;
end

% start training
shuffled_inds   = {};
while(iter_ < conf.max_iter)
    % train the localization model (for conf.test_interval number of iterations)
    [mean_error_train, shuffled_inds] = TrainNet(conf, model, solver, image_db_train, shuffled_inds);
    % test the localization model (for conf.test_iter number of iterations)
    mean_error_test = TestNet( conf, model, solver, image_db_val);

    iter_ = solver.iter();
    ShowState(iter_, mean_error_train, mean_error_test);
    diary; diary; % flush diary

    last_finetuned_model_path = fullfile(conf.finetune_rst_dir, ...
        sprintf('%s_iter_%d.caffemodel', conf.snapshot_prefix, iter_));
    assert(exist(last_finetuned_model_path,'file')>0);
    fprintf('last save as %s\n', last_finetuned_model_path);

    % delete the previously saved model
    deletePreviousModel(last_finetuned_model_prev_path);
    last_finetuned_model_prev_path = last_finetuned_model_path;
end
cd(current_dir);
finetuned_model_path = last_finetuned_model_path;
fprintf('Optimization is finished\n')
end

function [solver, iter_] = InitializeSolver(conf)
solver = caffe.Solver(conf.finetune_net_def_file);
iter_   = 1;

if isfield(conf,'solver_state_file') && ~isempty(conf.solver_state_file)
    solver.restore(conf.solver_state_file);
    [directory,filename] = fileparts(conf.solver_state_file);
    model_file = [directory,filesep, filename, '.caffemodel'];
    solver.net.copy_from(model_file);
    solver.test_nets(1).copy_from(model_file);
    iter_   = solver.iter();
    rng('shuffle')
else
    solver.net.copy_from(conf.net_file);
    solver.test_nets(1).copy_from(conf.net_file);
end
end

function deletePreviousModel(finetuned_model_prev_path)
if exist(finetuned_model_prev_path, 'file')
	delete(finetuned_model_prev_path);
end
[directory,filename] = fileparts(finetuned_model_prev_path);
finetuned_model_prev_path_solver_state = [directory,filesep, filename,'.solverstate'];
if exist(finetuned_model_prev_path_solver_state, 'file')
    delete(finetuned_model_prev_path_solver_state);
end
end

function ShowState(iter, train_error, test_error)
fprintf('\n------------------------- Iteration %d -------------------------\n', iter);
fprintf('Error - Training : %.4f - Testing : %.4f\n', train_error, test_error);
end

function [mean_loc_error] = TestNet(conf, model, solver, image_db_val)
fprintf('Testing: \n'); th = tic;

% sample testing images and testing regions
conf.save_data_test_file = fullfile(conf.finetune_rst_dir, ...
    sprintf('test_file_ims_per_batch%d_batch_size%d_test_iter%d.mat', ...
    conf.ims_per_batch, conf.batch_size, conf.test_iter));
try
    ld = load(conf.save_data_test_file);
    shuffled_inds_val    = ld.shuffled_inds_val;
    windows_bbox_loc_val = ld.windows_bbox_loc_val;
catch
    % sample test images
    shuffled_inds_val    = GenerateRandomMinibatch([], image_db_val, conf.ims_per_batch);
    shuffled_inds_val    = shuffled_inds_val(randperm(length(shuffled_inds_val), conf.test_iter));
    windows_bbox_loc_val = cell(length(shuffled_inds_val),1);
    % sample test regions from the test images
    for i = 1:length(shuffled_inds_val)
        windows_bbox_loc_val{i} = SampleImageRegionsForLoc(conf, image_db_val.images(shuffled_inds_val{i}));
    end
    save(conf.save_data_test_file, 'shuffled_inds_val', 'windows_bbox_loc_val');
end

mean_loc_error = 0;
mean_loc_loss  = 0;

bbox_inits_loc  = cell(length(shuffled_inds_val),1);
predictions_loc = cell([1,1,1,length(shuffled_inds_val)]);
bbox_target_loc = cell(length(shuffled_inds_val),1);
for i = 1:length(shuffled_inds_val)
    image_roidb_this_minibath = image_db_val.images(shuffled_inds_val{i}); % validation image
    windows_bbox_loc_this = windows_bbox_loc_val{i}; % sampled regions in the validation image

    % get image caffe blobs
    [im_blob, im_scales, image_sizes] = GetImageBlob(model, image_roidb_this_minibath);

    % get regions caffe blobs and target blobs
    [loc_blobs, bbox_target_loc{i}, bbox_inits_loc{i}] = ...
        GetLocNetMinibatchData(model.loc_params, model, image_roidb_this_minibath, ...
        windows_bbox_loc_this, im_scales, image_sizes);
    % set network input
    net_inputs = [{im_blob}, loc_blobs];

    caffe_reshape_net_as_input(solver.test_nets(1), net_inputs);
    caffe_set_input(solver.test_nets(1), net_inputs);
    solver.test_nets(1).forward_prefilled();

    % get output
    mean_loc_loss        = mean_loc_loss + solver.test_nets(1).blobs('loss_loc').get_data();
    predictions_loc{i}   = solver.test_nets(1).blobs('preds_loc').get_data();
    loc_accuracy         = get_loc_accuracy(model.loc_params, predictions_loc{i}, loc_blobs);
    mean_loc_error       = mean_loc_error + (1-loc_accuracy);
end

[mAR_before, mAR_after, mAR_unbound] = ...
    get_mean_Average_Recall_stats(bbox_inits_loc, bbox_target_loc, predictions_loc, model.loc_params);
mean_loc_error = mean_loc_error / length(shuffled_inds_val);
mean_loc_loss  = mean_loc_loss  / length(shuffled_inds_val);

fprintf('elapsed time %.2f sec | localization: loss %.4f error %.4f mAR %.2f --> mAR %.2f (upper bound %.2f)\n', ...
    toc(th), mean_loc_loss, mean_loc_error, ...
    mAR_before*100, mAR_after*100, mAR_unbound*100);
end

function [loc_accuracy] = get_loc_accuracy(conf, predictions_loc, loc_blobs)
% compute the accuracy of the predicted probability vectors w.r.t. the
% target probability vectors (only for the LocNet models).

switch conf.loc_type
    case {'inout','borders','combined'}
        target_vectors = loc_blobs{2};
        target_weights = loc_blobs{3};
        mask           = target_weights > 0;
        preds_vectors  = single(predictions_loc > 0.5);
        loc_accuracy   = sum(target_vectors(mask) == preds_vectors(mask)) / sum(mask(:));
    otherwise
        loc_accuracy = 0;
end

end

function [mAR_before, mAR_after, mAR_unbound] = ...
    get_mean_Average_Recall_stats(bbox_inits_loc, bbox_target_loc, predictions_loc, conf)

bbox_target_loc = cell2mat(bbox_target_loc);
bbox_inits_loc  = cell2mat(bbox_inits_loc);
class_indices   = bbox_inits_loc(:,5);
bbox_inits_loc  = bbox_inits_loc(:,1:4);

switch conf.loc_type
    case 'bboxreg'
        predictions_loc = squeeze(cell2mat(predictions_loc(:)'));

        % permute from [(4*num_classes) x NumBoxes] -> [NumBoxes x (4*num_classes)]
        predictions_loc = single(permute(predictions_loc, [2, 1]));
        % reshape from [NumBoxes x (4*num_classes)] -> [NumBoxes x 4 x num_classes]
        predictions_loc = reshape(predictions_loc, ...
            [size(predictions_loc,1), 4, conf.num_classes]);

        pred_bbox_loc = decode_reg_vals_to_bbox_targets(...
            bbox_inits_loc, predictions_loc, class_indices);

        if nargout >= 3
            % find the ideal predictions...
            target_reg_values = encode_bbox_targets_to_reg_vals(...
                bbox_inits_loc, bbox_target_loc, conf.num_classes);

            pred_bbox_loc_ideal = decode_reg_vals_to_bbox_targets(...
                bbox_inits_loc, target_reg_values, class_indices);

            % upper bound on the mAR with the CNN-based bounding box
            % regression model. Normaly mAR = 1.0.
            mAR_unbound = compute_mAR_of_bboxes(bbox_inits_loc, pred_bbox_loc_ideal(:,1:4), ...
                bbox_target_loc(:,1:4), class_indices, conf.num_classes);
        end
    case {'inout','borders','combined'}
        predictions_loc = squeeze(cell2mat(predictions_loc));

        pred_bbox_loc = decode_loc_probs_to_bbox_targets(...
            bbox_inits_loc, class_indices, predictions_loc, conf);

        if nargout >= 3
            % find the ideal predictions...
            target_loc_prob_vectors = encode_bbox_target_to_loc_probs(...
                bbox_inits_loc, bbox_target_loc, conf);

            pred_bbox_loc_ideal = decode_loc_probs_to_bbox_targets(...
                bbox_inits_loc, class_indices, target_loc_prob_vectors, conf);

            % upper bound on the mAR with LocNet model. Because of the
            % quantization of the search region in M horizontal stripes and
            % M vertical stripes the mAR < 1.0
            mAR_unbound = compute_mAR_of_bboxes(bbox_inits_loc, pred_bbox_loc_ideal(:,1:4), ...
                bbox_target_loc(:,1:4), class_indices, conf.num_classes);
        end
    otherwise
        error('Invalid localization type %s',conf.loc_type)
end

mAR_before = compute_mAR_of_bboxes(bbox_inits_loc, bbox_inits_loc, ...
            bbox_target_loc(:,1:4), class_indices, conf.num_classes);

mAR_after = compute_mAR_of_bboxes(bbox_inits_loc, pred_bbox_loc(:,1:4), ...
            bbox_target_loc(:,1:4), class_indices, conf.num_classes);
end

function [mean_loc_error, shuffled_inds] = TrainNet(conf, model, solver, image_db_train, shuffled_inds)
mean_loc_error = 0;
mean_loc_loss = 0;
fprintf('Training:\n'); th = tic;
for i = 1:conf.test_interval
    % sample training images
    [shuffled_inds, sub_db_inds] = GenerateRandomMinibatch(...
        shuffled_inds, image_db_train, conf.ims_per_batch);
    image_roidb_this_minibath = image_db_train.images(sub_db_inds);

    % get the training image blobs that will be fed to caffe
    [im_blob, im_scales, image_sizes] = GetImageBlob(model, image_roidb_this_minibath);

    % sample training regions for localization
    windows_bbox_loc = SampleImageRegionsForLoc(conf, image_roidb_this_minibath);

    % get training regions caffe blobs and the target blobs
    loc_blobs = GetLocNetMinibatchData(model.loc_params, model, ...
        image_roidb_this_minibath, windows_bbox_loc, im_scales, image_sizes);
    % set network input
    net_inputs = [{im_blob}, loc_blobs];

    caffe_reshape_net_as_input(solver.net, net_inputs);
    caffe_set_input(solver.net, net_inputs);
    solver.step(1);

    % get output
    mean_loc_loss   = mean_loc_loss + solver.net.blobs('loss_loc').get_data();
    predictions_loc = solver.net.blobs('preds_loc').get_data();
    loc_accuracy    = get_loc_accuracy(model.loc_params, predictions_loc, loc_blobs);
    mean_loc_error  = mean_loc_error + (1-loc_accuracy);

    if (mod(i, conf.display_step) == 0)
        fprintf('%2d/%d:\telapsed time %.2f sec --- error %.4f loss %.4f\n', ...
            i, conf.test_interval, toc(th), mean_loc_error/i, mean_loc_loss/i);
        th = tic;
        diary; diary; % flush diary
    end
end
mean_loc_error = mean_loc_error/i;
end

function [mAR, AR] = compute_mAR_of_bboxes(all_bbox_init, all_bbox_pred, all_bbox_targets, class_ids, num_classes)
% compute the mAR of the predicted bounding boxes w.r.t. the target
% bounding boxes.

AR = zeros(1,num_classes);
fprintf('AveRecall: ')
for c = 1:num_classes
    is_this_cls  = class_ids==c;
    bbox_pred    = all_bbox_pred(is_this_cls,:);
    bbox_targets = all_bbox_targets(is_this_cls,:);
    is_not_gt    = ~all(all_bbox_init(is_this_cls,:) == bbox_targets, 2);
    bbox_pred    = bbox_pred(is_not_gt,:);
    bbox_targets = bbox_targets(is_not_gt,:);
    AR(c)        = compute_ave_recall_of_bbox(bbox_pred, bbox_targets);
    fprintf('[%d]=%.2f(%d) ',c, AR(c),size(bbox_targets,1));
end

fprintf('\n')
mAR = mean(AR);
end

function [shuffled_inds, sub_inds] = GenerateRandomMinibatch(...
    shuffled_inds, image_db, ims_per_batch)
% This function comes from the Faster R-CNN code and the matlab
% re-implementation of Fast-RCNN (https://github.com/ShaoqingRen/faster_rcnn)

% shuffle training data per batch
if isempty(shuffled_inds)
    fprintf('!!!!!!!!!!!!!!!!!!!!!!!!New Epoch!!!!!!!!!!!!!!!!!!!!!!!!!!!\n');
    % make sure each minibatch, only has horizontal images or vertical
    % images, to save gpu memory
    hori_image_inds = image_db.image_sizes(:,2) >= image_db.image_sizes(:,1);
    vert_image_inds = ~hori_image_inds;
    hori_image_inds = single(find(hori_image_inds));
    vert_image_inds = single(find(vert_image_inds));

    % random perm
    lim = floor(length(hori_image_inds) / ims_per_batch) * ims_per_batch;
    hori_image_inds = hori_image_inds(randperm(length(hori_image_inds), lim));
    lim = floor(length(vert_image_inds) / ims_per_batch) * ims_per_batch;
    vert_image_inds = vert_image_inds(randperm(length(vert_image_inds), lim));

    % combine sample for each ims_per_batch
    hori_image_inds = reshape(hori_image_inds, ims_per_batch, []);
    vert_image_inds = reshape(vert_image_inds, ims_per_batch, []);

    shuffled_inds = [hori_image_inds, vert_image_inds];
    shuffled_inds = shuffled_inds(:, randperm(size(shuffled_inds, 2)));
    shuffled_inds = num2cell(shuffled_inds, 1);
end

if nargout > 1
    % generate minibatch training data
    sub_inds = shuffled_inds{1};
    assert(length(sub_inds) == ims_per_batch);
    shuffled_inds(1) = [];
end

end

function sampled_regions_loc = SampleImageRegionsForLoc(conf, image_roidb)

num_images = length(image_roidb);
assert(mod(conf.batch_size, num_images) == 0, ...
    sprintf('num_images %d must divide BATCH_SIZE %d', num_images, conf.batch_size));

rois_per_image  = conf.batch_size / num_images; % regions per image
num_images  = length(image_roidb);
sampled_regions_loc = cell(num_images, 1);
for i = 1:num_images
    num_regions = size(image_roidb(i).windows_bbox_loc,1);
    % sample rois_per_image number of regions from this image
    indices = randperm(num_regions, min(rois_per_image, num_regions));
    windows_bbox_loc_this = image_roidb(i).windows_bbox_loc(indices,:);
    region_inds_loc = windows_bbox_loc_this(:,1);
    % load the entire set of training regions of this image (cached in disk)
    regions = load_img_regions(image_roidb(i).regions);
    % keep only the regions of this minibatch
    bboxes = regions(region_inds_loc,1:4);
    sampled_regions_loc{i} = single([bboxes,windows_bbox_loc_this(:,2:3)]);
end
end

function [im_blob, im_scales, image_sizes, im_scaled_size] = GetImageBlob(model, image_roidb)
% This function comes from the Faster R-CNN code and the matlab
% re-implementation of Fast-RCNN (https://github.com/ShaoqingRen/faster_rcnn)

num_images     = length(image_roidb);
random_scale_inds = randi(length(model.scales), num_images, 1);
processed_imgs = cell(num_images, 1);
im_scales      = nan(num_images, 1);
image_sizes    = nan(num_images, 2);
im_scaled_size = nan(num_images, 2);
for i = 1:num_images
    img = get_image(image_roidb(i).image_path);
    image_sizes(i,1) = size(img,1);
    image_sizes(i,2) = size(img,2);
    target_size = model.scales(random_scale_inds(i));
    [processed_imgs{i}, im_scales(i), im_scaled_size(i,:)] = prepare_img_blob(img, ...
        target_size, model.mean_pix, model.max_size);
end

im_blob = img_list_to_blob(processed_imgs);
end

function image_db = setup_image_db(image_db, conf)
num_imgs = length(image_db.image_paths);
image_db.all_windows_bbox_loc = cell(num_imgs, 1);
num_regions_fg_loc = zeros(num_imgs, 1);

cache_regions_dir = image_db.cache_regions_dir;
all_region_files  = strcat(cache_regions_dir, getImageIdsFromImagePaths( image_db.image_paths ),'.mat');

for i = 1:num_imgs
    % find for each image its training bounding boxes and their
    % corresponding target bounding boxes
    [image_db.all_regions{i}, image_db.all_windows_bbox_loc{i}, num_regions_fg_loc(i)] = ...
        setup_img_data(image_db.all_bbox_gt{i}, image_db.all_regions{i}, conf);

    % cache the training bounding boxes in the hard disk (for large
    % datasets, such as MSCOCO, is in-efficient memory wise to keep the
    % training boxes in the memory)
    save_img_regions(image_db.all_regions{i}, all_region_files{i});
    image_db.all_regions{i} = [];
    tic_toc_print('setup img data %d / % d (%.2f)\n', i, num_imgs, 100*i/num_imgs);
end

image_db.all_regions = [];
image_db.all_regions = all_region_files;

is_valid = (num_regions_fg_loc > 0);

fprintf('valid: %d/%d (%.3f) \n', sum(is_valid),length(is_valid),sum(is_valid)/length(is_valid));

image_db.image_paths          = image_db.image_paths(is_valid);
image_db.all_regions          = image_db.all_regions(is_valid);
image_db.all_windows_bbox_loc = image_db.all_windows_bbox_loc(is_valid);
image_db.all_bbox_gt          = image_db.all_bbox_gt(is_valid);
num_regions_fg_loc            = num_regions_fg_loc(is_valid);

image_db.images = cell2struct(...
    [image_db.image_paths, image_db.all_regions, ...
     image_db.all_windows_bbox_loc, image_db.all_bbox_gt], ...
    {'image_path', 'regions', 'windows_bbox_loc', 'bbox_gt'}, 2);

image_db = rmfield(image_db, 'image_paths');
image_db = rmfield(image_db, 'all_regions');
image_db = rmfield(image_db, 'all_windows_bbox_loc');
image_db = rmfield(image_db, 'all_bbox_gt');

fprintf('Total number of fg regions for localization: %d\n', sum(num_regions_fg_loc))

num_classes = length(image_db.classes);
fg_per_class = zeros(num_classes,1);
gt_per_class = zeros(num_classes,1);
labels_per_img = -ones(num_imgs,num_classes,'single');
for i = 1:num_imgs
    class_ids         = image_db.images(i).windows_bbox_loc(:,end);
    fg_per_class_this = accumarray(class_ids, 1, [num_classes, 1]);
    fg_per_class      = fg_per_class+fg_per_class_this;
    gt_per_class_this = accumarray(image_db.images(i).bbox_gt(:,5),1,[num_classes, 1]);
    gt_per_class      = gt_per_class + gt_per_class_this;

    unique_class_ids = unique(class_ids);
    labels_per_img(i,unique_class_ids) = 1;
end

% print the number of training boxes and grount truth bounding boxes per
% category
fprintf('Regions / GTs per class:\n');
for i = 1:num_classes
    fprintf('Category %d: FG %d GT %d R %.2f\n',i, fg_per_class(i), gt_per_class(i), fg_per_class(i)/gt_per_class(i));
end

image_db.gt_per_class   = gt_per_class;
image_db.fg_per_class   = fg_per_class;
image_db.labels_per_img = labels_per_img;
end

function [bbox_proposals, windows_bbox_loc, num_regions_fg_loc] = ...
    setup_img_data(bbox_gt, bbox_proposals, conf)

num_classes           = conf.num_classes;
fg_threshold_bbox_loc = conf.fg_threshold_bbox_loc;

bbox_gt_label = bbox_gt(:,5);
bbox_all = [bbox_gt(:,1:4); bbox_proposals(:,1:4)];
overlaps = zeros(size(bbox_all,1), num_classes, 'single');

windows_bbox_loc = cell(num_classes, 1);
for class_idx = 1:num_classes
    this_class_ind = find(bbox_gt_label == class_idx);

    windows_bbox_loc{class_idx} = zeros(0,3,'single');
    if ~isempty(this_class_ind)
        overlap_class = single(boxoverlap(bbox_all(:,1:4), bbox_gt(this_class_ind,1:4)));
        [overlaps(:, class_idx), gt_index] = max(overlap_class, [], 2);

        this_class_pr_index = find(overlaps(:, class_idx) >= fg_threshold_bbox_loc);
        if ~isempty(this_class_pr_index)
            this_class_gt_index = this_class_ind(gt_index);
            this_class_gt_index = this_class_gt_index(this_class_pr_index);
            windows_bbox_loc{class_idx} = ...
                single([this_class_pr_index(:), this_class_gt_index(:), ...
                ones(length(this_class_pr_index),1)*class_idx]);
        end
    end
end
windows_bbox_loc = single(cell2mat(windows_bbox_loc)); % label, overlap, bbox_init, bbox_target
num_regions_fg_loc = size(windows_bbox_loc,1);

% remove unecesary bbox proposals...
pr_ind = windows_bbox_loc(:,1);
[pr_ind_u, ~, pr_ind_new] = unique(pr_ind);
bbox_proposals = bbox_all(pr_ind_u,1:4);
assert(all(all(bbox_proposals(pr_ind_new,1:4)==bbox_all(pr_ind,1:4))))
windows_bbox_loc(:,1) = pr_ind_new;
end

function save_img_regions(regions, regions_file)
if ~(exist(regions_file,'file')>0)
    save(regions_file, 'regions', '-v7.3');
end
end

function regions = load_img_regions(regions_file)
load(regions_file, 'regions');
end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
function par = parametrize(varargin)
%PARAMETRIZE Parametrization of a polyline, based on edges lengths.
%
%   PAR = parametrize(POLY);
%   Returns a parametrization of the curve defined by the serie of points,
%   based on euclidean distance between two consecutive points.
%   POLY is a N-by-2 array, representing coordinates of vertices. The
%   result PAR is N-by-1, and contains the cumulative length of edges until
%   corresponding vertex.
%
%   PAR = parametrize(PX, PY);
%   is the same, but specify points coordinates in separate column vectors.
%
%   PAR = parametrize(..., 'normalize', 1);
%   PAR = parametrize(..., 'normalize', true);
%   Rescales the result such that the last element of PAR is 1.
%
%   Example
%     % Parametrize a circle approximation
%     poly = circleToPolygon([0 0 1], 200);
%     p = parametrize(poly);
%     p(end)
%     ans =
%         6.2829
%
%   See also
%   polygons2d, polylineLength

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2003-04-06
% Copyright 2003-2022 INRA - TPV URPOI - BIA IMASTE

%% Process inputs

% extract vertex coordinates
if size(varargin{1}, 2) > 1
    % vertices in a single array
    pts = varargin{1};
    varargin(1) = [];

elseif length(varargin) == 2
    % points as separate arrays
    pts = [varargin{1} varargin{2}];
    varargin(1:2) = [];

end

% by default, do not normalize
normalize = false;

% extract options
while length(varargin) > 1
    param = varargin{1};
    switch lower(param)
        case 'normalize'
            normalize = varargin{2};
        otherwise
            error('Unknown parameter name: %s', param);
    end
    varargin(1:2) = [];
end

%% Parametrize polyline

% compute cumulative sum of euclidean distances between consecutive
% vertices, setting distance of first vertex to 0.
if size(pts, 2) == 2
    % process points in 2D
    par = [0 ; cumsum(hypot(diff(pts(:,1)), diff(pts(:,2))))];
else
    % process points in arbitrary dimension
    par = [0 ; cumsum(sqrt(sum(diff(pts).^2, 2)))];
end

% eventually rescale between 0 and 1
if normalize
    par = par / par(end);
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom