text:
%%test network
cam1_size=size(label_test_cam1);
cam2_size=size(label_test_cam2);
cam3_size=size(label_test_cam3);
cam1_feature1=[];
cam1_feature2=[];
cam1_feature3=[];
cam1_score=[];
test_batch_data=zeros(224,224,3,param.test_batch_size,'single');
for m=0:floor(cam1_size(1)/param.test_batch_size)-1
    fprintf('process cam1:%d/%d\n',m,floor(cam1_size(1)/param.test_batch_size)-1);
    for n=1:param.test_batch_size
        im_data=imread(test_image_name_cam1{m*param.test_batch_size+n});
        im_data=imresize(single(im_data),[224 224]);
        im_data=im_data(:,:,[3,2,1]);
        im_data=permute(im_data,[2,1,3]);
        test_batch_data(:,:,:,n)=im_data;
        test_batch_data(:,:,1,n)=test_batch_data(:,:,1,n)-104;
        test_batch_data(:,:,2,n)=test_batch_data(:,:,2,n)-117;
        test_batch_data(:,:,3,n)=test_batch_data(:,:,3,n)-123;
    end
    net.blobs('data').set_data(test_batch_data);
    net.forward_prefilled;
    cam1_feature1=[cam1_feature1;(squeeze(net.blobs('data1_p_1').get_data))'];
    cam1_feature2=[cam1_feature2;(squeeze(net.blobs('data1_p_2').get_data))'];
    cam1_feature3=[cam1_feature3;(squeeze(net.blobs('data1_p_3').get_data))'];
    cam1_score=[cam1_score;(squeeze(net.blobs('sig').get_data))'];
end
test_batch_data=zeros(224,224,3,param.test_batch_size,'single');
for m=(floor(cam1_size(1)/param.test_batch_size))*param.test_batch_size+1:cam1_size(1)
    fprintf('process cam1:%d/%d\n',m,cam1_size(1));
        count=m-(floor(cam1_size(1)/param.test_batch_size))*param.test_batch_size;
        im_data=imread(test_image_name_cam1{m});
        im_data=imresize(single(im_data),[224 224]);
        im_data=im_data(:,:,[3,2,1]);
        im_data=permute(im_data,[2,1,3]);
        test_batch_data(:,:,:,count)=im_data;
        test_batch_data(:,:,1,count)=test_batch_data(:,:,1,count)-104;
        test_batch_data(:,:,2,count)=test_batch_data(:,:,2,count)-117;
        test_batch_data(:,:,3,count)=test_batch_data(:,:,3,count)-123;
end
net.blobs('data').set_data(test_batch_data);
net.forward_prefilled;
index=cam1_size(1)-(floor(cam1_size(1)/param.test_batch_size))*param.test_batch_size;
result1=(squeeze(net.blobs('data1_p_1').get_data))';
result2=(squeeze(net.blobs('data1_p_2').get_data))';
result3=(squeeze(net.blobs('data1_p_3').get_data))';
score=(squeeze(net.blobs('sig').get_data))';
cam1_feature1=[cam1_feature1;result1(1:index,:)];
cam1_feature2=[cam1_feature2;result2(1:index,:)];
cam1_feature3=[cam1_feature3;result3(1:index,:)];
cam1_score=[cam1_score;score(1:index,:)];
gallery_feature=[];
gallery_label=[];
for m=1:param.test_person_num
    test_label_index=find(label_test_cam1==(m));
    if isempty(test_label_index)
        continue;
    end
    a=bsxfun(@times,cam1_feature1(test_label_index,:),cam1_score(test_label_index,1));
    b=bsxfun(@times,cam1_feature2(test_label_index,:),cam1_score(test_label_index,2));
    c=bsxfun(@times,cam1_feature3(test_label_index,:),cam1_score(test_label_index,3));
    gallery_feature=[gallery_feature;sum(a/sum(cam1_score(test_label_index,1))),sum(b/sum(cam1_score(test_label_index,2))),sum(c/sum(cam1_score(test_label_index,3)))];
%     gallery_feature=[gallery_feature;sum(a/sum(cam1_score(test_label_index,1)))];
%     gallery_feature=[gallery_feature;sum(b/sum(cam1_score(test_label_index,2)))];
%     gallery_feature=[gallery_feature;sum(c/sum(cam1_score(test_label_index,3)))];
    gallery_label=[gallery_label;m];
end
%extract feature cam3
cam3_feature1=[];
cam3_feature2=[];
cam3_feature3=[];
cam3_score=[];
test_batch_data=zeros(224,224,3,param.test_batch_size,'single');
for m=0:floor(cam3_size(1)/param.test_batch_size)-1
        fprintf('process cam3:%d/%d\n',m,floor(cam3_size(1)/param.test_batch_size)-1);
    for n=1:param.test_batch_size
        im_data=imread(test_image_name_cam3{m*param.test_batch_size+n});
        im_data=imresize(single(im_data),[224 224]);
        im_data=im_data(:,:,[3,2,1]);
        im_data=permute(im_data,[2,1,3]);
        test_batch_data(:,:,:,n)=im_data;
        test_batch_data(:,:,1,n)=test_batch_data(:,:,1,n)-104;
        test_batch_data(:,:,2,n)=test_batch_data(:,:,2,n)-117;
        test_batch_data(:,:,3,n)=test_batch_data(:,:,3,n)-123;
    end
    net.blobs('data').set_data(test_batch_data);
    net.forward_prefilled;
    cam3_feature1=[cam3_feature1;(squeeze(net.blobs('data1_p_1').get_data))'];
    cam3_feature2=[cam3_feature2;(squeeze(net.blobs('data1_p_2').get_data))'];
    cam3_feature3=[cam3_feature3;(squeeze(net.blobs('data1_p_3').get_data))'];
    cam3_score=[cam3_score;(squeeze(net.blobs('sig').get_data))'];
end
test_batch_data=zeros(224,224,3,param.test_batch_size,'single');
for m=(floor(cam3_size(1)/param.test_batch_size))*param.test_batch_size+1:cam3_size(1)
        fprintf('process cam3:%d/%d\n',m,cam3_size(1));
        count=m-(floor(cam3_size(1)/param.test_batch_size))*param.test_batch_size;
        im_data=imread(test_image_name_cam3{m});
        im_data=imresize(single(im_data),[224 224]);
        im_data=im_data(:,:,[3,2,1]);
        im_data=permute(im_data,[2,1,3]);
        test_batch_data(:,:,:,count)=im_data;
        test_batch_data(:,:,1,count)=test_batch_data(:,:,1,count)-104;
        test_batch_data(:,:,2,count)=test_batch_data(:,:,2,count)-117;
        test_batch_data(:,:,3,count)=test_batch_data(:,:,3,count)-123;
end
net.blobs('data').set_data(test_batch_data);
net.forward_prefilled;
index=cam3_size(1)-(floor(cam3_size(1)/param.test_batch_size))*param.test_batch_size;
result1=(squeeze(net.blobs('data1_p_1').get_data))';
result2=(squeeze(net.blobs('data1_p_2').get_data))';
result3=(squeeze(net.blobs('data1_p_3').get_data))';
score=(squeeze(net.blobs('sig').get_data))';
cam3_feature1=[cam3_feature1;result1(1:index,:)];
cam3_feature2=[cam3_feature2;result2(1:index,:)];
cam3_feature3=[cam3_feature3;result3(1:index,:)];
cam3_score=[cam3_score;score(1:index,:)];
for m=1:param.test_person_num
    test_label_index=find(label_test_cam3==(m));
    if isempty(test_label_index)
        continue;
    end
    a=bsxfun(@times,cam3_feature1(test_label_index,:),cam3_score(test_label_index,1));
    b=bsxfun(@times,cam3_feature2(test_label_index,:),cam3_score(test_label_index,2));
    c=bsxfun(@times,cam3_feature3(test_label_index,:),cam3_score(test_label_index,3));
    gallery_feature=[gallery_feature;sum(a/sum(cam3_score(test_label_index,1))),sum(b/sum(cam3_score(test_label_index,2))),sum(c/sum(cam3_score(test_label_index,3)))];
%     gallery_feature=[gallery_feature;sum(a/sum(cam3_score(test_label_index,1)))];
%     gallery_feature=[gallery_feature;sum(b/sum(cam3_score(test_label_index,2)))];
%     gallery_feature=[gallery_feature;sum(c/sum(cam3_score(test_label_index,3)))];
    gallery_label=[gallery_label;m];
end
%extract feature cam2
cam2_feature1=[];
cam2_feature2=[];
cam2_feature3=[];
cam2_score=[];
test_batch_data=zeros(224,224,3,param.test_batch_size,'single');
for m=0:floor(cam2_size(1)/param.test_batch_size)-1
        fprintf('process cam2:%d/%d\n',m,floor(cam2_size(1)/param.test_batch_size)-1);
    for n=1:param.test_batch_size
        im_data=imread(test_image_name_cam2{m*param.test_batch_size+n});
        im_data=imresize(single(im_data),[224 224]);
        im_data=im_data(:,:,[3,2,1]);
        im_data=permute(im_data,[2,1,3]);
        test_batch_data(:,:,:,n)=im_data;
        test_batch_data(:,:,1,n)=test_batch_data(:,:,1,n)-104;
        test_batch_data(:,:,2,n)=test_batch_data(:,:,2,n)-117;
        test_batch_data(:,:,3,n)=test_batch_data(:,:,3,n)-123;
    end
    net.blobs('data').set_data(test_batch_data);
    net.forward_prefilled;
    cam2_feature1=[cam2_feature1;(squeeze(net.blobs('data1_p_1').get_data))'];
    cam2_feature2=[cam2_feature2;(squeeze(net.blobs('data1_p_2').get_data))'];
    cam2_feature3=[cam2_feature3;(squeeze(net.blobs('data1_p_3').get_data))'];
    cam2_score=[cam2_score;(squeeze(net.blobs('sig').get_data))'];
end
test_batch_data=zeros(224,224,3,param.test_batch_size,'single');
for m=(floor(cam2_size(1)/param.test_batch_size))*param.test_batch_size+1:cam2_size(1)
        fprintf('process cam2:%d/%d\n',m,cam2_size(1));
        count=m-(floor(cam2_size(1)/param.test_batch_size))*param.test_batch_size;
        im_data=imread(test_image_name_cam2{m});
        im_data=imresize(single(im_data),[224 224]);
        im_data=im_data(:,:,[3,2,1]);
        im_data=permute(im_data,[2,1,3]);
        test_batch_data(:,:,:,count)=im_data;
        test_batch_data(:,:,1,count)=test_batch_data(:,:,1,count)-104;
        test_batch_data(:,:,2,count)=test_batch_data(:,:,2,count)-117;
        test_batch_data(:,:,3,count)=test_batch_data(:,:,3,count)-123;
end
net.blobs('data').set_data(test_batch_data);
net.forward_prefilled;
index=cam2_size(1)-(floor(cam2_size(1)/param.test_batch_size))*param.test_batch_size;
result1=(squeeze(net.blobs('data1_p_1').get_data))';
result2=(squeeze(net.blobs('data1_p_2').get_data))';
result3=(squeeze(net.blobs('data1_p_3').get_data))';
score=(squeeze(net.blobs('sig').get_data))';
cam2_feature1=[cam2_feature1;result1(1:index,:)];
cam2_feature2=[cam2_feature2;result2(1:index,:)];
cam2_feature3=[cam2_feature3;result3(1:index,:)];
cam2_score=[cam2_score;score(1:index,:)];
prob_feature=[];
prob_label=[];
for m=1:param.test_person_num
    test_label_index=find(label_test_cam2==(m));
    if isempty(test_label_index)
        continue;
    end
    a=bsxfun(@times,cam2_feature1(test_label_index,:),cam2_score(test_label_index,1));
    b=bsxfun(@times,cam2_feature2(test_label_index,:),cam2_score(test_label_index,2));
    c=bsxfun(@times,cam2_feature3(test_label_index,:),cam2_score(test_label_index,3));
    prob_feature=[prob_feature;sum(a/sum(cam2_score(test_label_index,1))),sum(b/sum(cam2_score(test_label_index,2))),sum(c/sum(cam2_score(test_label_index,3)))];

%     prob_feature=[prob_feature;sum(a/sum(cam2_score(test_label_index,1)))];
%     prob_feature=[prob_feature;sum(b/sum(cam2_score(test_label_index,2)))];
%     prob_feature=[prob_feature;sum(c/sum(cam2_score(test_label_index,3)))];

    prob_label=[prob_label;m];
end
% cal cmc
prob_norm=bsxfun(@rdivide,prob_feature,sum(abs(prob_feature).^2,2).^(1/2));
gallery_norm=bsxfun(@rdivide,gallery_feature,sum(abs(gallery_feature).^2,2).^(1/2));
score_matrix=prob_norm*gallery_norm';
rank1_hit=0;
rank5_hit=0;
rank10_hit=0;
rank20_hit=0;
problen=size(prob_feature,1);
for m=1:problen
    [~,location_temp]=sort(score_matrix(m,:),2,'descend');
    location=gallery_label(location_temp);
    if find(location(1)==prob_label(m))
        rank1_hit=rank1_hit+1;
    end
    if find(location(1:5)==prob_label(m))
        rank5_hit=rank5_hit+1;
    end
    if find(location(1:10)==prob_label(m))
        rank10_hit=rank10_hit+1;
    end
    if find(location(1:20)==prob_label(m))
        rank20_hit=rank20_hit+1;
    end
end
rank_acc=[rank1_hit/problen,rank5_hit/problen,rank10_hit/problen,rank20_hit/problen];
fin=fopen(param.result_save_file,'a');
fprintf(fin,'up split_index:%d,iter:%d, rank1:%f,rank5:%f,rank10:%f,rank20:%f\n',split_index,iter,rank_acc(1),rank_acc(2),rank_acc(3),rank_acc(4));
fprintf('split_index:%d,iter:%d, rank1:%f,rank5:%f,rank10:%f,rank20:%f\n',split_index,iter,rank_acc(1),rank_acc(2),rank_acc(3),rank_acc(4));
fclose(fin);

author: liuyuisanai

repo: Quality-Aware-Network


text:
%GROUPRECTANGLES_MEANSHIFT  Groups the object candidate rectangles using meanshift
%
%     [rects, weights] = cv.groupRectangles_meanshift(rects, weights, scales)
%     [...] = cv.groupRectangles_meanshift(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __rects__ Input cell array of rectangles, where each rectangle is
%   represented as a 4-element vector `{[x,y,w,h], ...}`, or a numeric
%   Nx4/Nx1x4/1xNx4 array.
% * __weights__ Input vector of associated weights.
% * __scales__ Input vector of corresponding rectangles scales.
%
% ## Output
% * __rects__ Output array of retained and grouped rectangles. Same format as
%   input `rects` (either Nx4 numeric array or cell array).
% * __weights__ Output updated weights.
%
% ## Options
% * __DetectThreshold__ detection threshold (weight) above which resulting
%   modes are kept. default 0.0
% * __WinDetSize__ window size `[w,h]`. default [64,128]
%
% See also: cv.groupRectangles, cv.SimilarRects
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [PSNR, yRGB_est] = CBM3D(yRGB, zRGB, sigma, profile, print_to_screen, colorspace)
%
%  CBM3D is algorithm for attenuation of additive white Gaussian noise from
%  color RGB images. This algorithm reproduces the results from the article:
%
%  [1] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, "Color image
%   denoising via sparse 3D collaborative filtering with grouping constraint in
%   luminance-chrominance space," submitted to IEEE Int. Conf. Image Process.,
%   January 2007, in review, preprint at http://www.cs.tut.fi/~foi/GCF-BM3D.
%
%  FUNCTION INTERFACE:
%
%  [PSNR, yRGB_est] = CBM3D(yRGB, zRGB, sigma, profile, print_to_screen, colorspace)
%
%  ! The function can work without any of the input arguments,
%   in which case, the internal default ones are used !
%
%  BASIC USAGE EXAMPLES:
%
%     Case 1) Using the default parameters (i.e., image name, sigma, etc.)
%
%      [PSNR, yRGB_est] = CBM3D;
%
%     Case 2) Using an external noisy image:
%
%      % Read an RGB image and scale its intensities in range [0,1]
%      yRGB = im2double(imread('image_House256rgb.png'));
%      % Generate the same seed used in the experimental results of [1]
%      randn('seed', 0);
%      % Standard deviation of the noise --- corresponding to intensity
%      %  range [0,255], despite that the input was scaled in [0,1]
%      sigma = 25;
%      % Add the AWGN with zero mean and standard deviation 'sigma'
%      zRGB = yRGB + (sigma/255)*randn(size(yRGB));
%      % Denoise 'zRGB'. The denoised image is 'yRGB_est', and 'NA = 1'
%      %  because the true image was not provided
%      [NA, yRGB_est] = CBM3D(1, zRGB, sigma);
%      % Compute the putput PSNR
%      PSNR = 10*log10(1/mean((yRGB(:)-yRGB_est(:)).^2))
%      % show the noisy image 'zRGB' and the denoised 'yRGB_est'
%      figure; imshow(min(max(zRGB,0),1));
%      figure; imshow(min(max(yRGB_est,0),1));
%
%     Case 3) If the original image yRGB is provided as the first input
%      argument, then some additional information is printed (PSNRs,
%      figures, etc.). That is, "[NA, yRGB_est] = BM3D(1, zRGB, sigma);" in the
%      above code should be replaced with:
%
%      [PSNR, yRGB_est] = CBM3D(yRGB, zRGB, sigma);
%
%
%  INPUT ARGUMENTS (OPTIONAL):
%     1) yRGB (M x N x 3): Noise-free RGB image (needed for computing PSNR),
%                           replace with the scalar 1 if not available.
%     2) zRGB (M x N x 3): Noisy RGBimage (intensities in range [0,1] or [0,255])
%     3) sigma (double)  : Std. dev. of the noise (corresponding to intensities
%                            in range [0,255] even if the range of zRGB is [0,1])
%     4) profile (char)  : 'np' --> Normal Profile
%                          'lc' --> Fast Profile
%     5) print_to_screen : 0 --> do not print output information (and do
%                                not plot figures)
%                          1 --> print information and plot figures
%     6) colorspace (char): 'opp'   --> use opponent colorspace
%                          'yCbCr' --> use yCbCr colorspace
%
%  OUTPUTS:
%     1) PSNR (double)          : Output PSNR (dB), only if the original
%                                 image is available, otherwise PSNR = 0
%     2) yRGB_est (M x N x 3): Final RGB estimate (in the range [0,1])
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright (c) 2007-2011 Tampere University of Technology.
% All rights reserved.
% This work should only be used for nonprofit purposes.
%
% AUTHORS:
%     Kostadin Dabov, email: dabov _at_ cs.tut.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% In case, there is no input image (zRGB or yRGB), then use the filename
%%%%  below to read an original image (might contain path also). Later,
%%%%  artificial AWGN noise is added and this noisy image is processed
%%%%  by the CBM3D.
%%%%
image_name = [
%    'kodim12.png'
   'image_Lena512rgb.png'
%     'image_House256rgb.png'
%    'image_Peppers512rgb.png'
%    'image_Baboon512rgb.png'
%    'image_F16_512rgb.png'
   ];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  Quality/complexity trade-off
%%%%
%%%%  'np' --> Normal Profile (balanced quality)
%%%%  'lc' --> Low Complexity Profile (fast, lower quality)
%%%%
%%%%  'high' --> High Profile (high quality, not documented in [1])
%%%%
%%%%  'vn' --> This profile is automatically enabled for high noise
%%%%           when sigma > 40
%%%%
%%%%  'vn_old' --> This is the old 'vn' profile that was used in [1].
%%%%           It gives inferior results than 'vn' in most cases.
%%%%
if (exist('profile') ~= 1)
    profile         = 'np'; %% default profile
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  Specify the std. dev. of the corrupting noise
%%%%
if (exist('sigma') ~= 1),
   sigma                = 50; %% default standard deviation of the AWGN
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  Colorspace in which we perform denoising. BM is applied to the first
%%%%  component and the matching information is re-used for the other two.
%%%%
if (exist('colorspace') ~= 1),
    colorspace              = 'opp'; %%% (valid colorspaces are: 'yCbCr' and 'opp')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Following are the parameters for the Normal Profile.
%%%%

%%%% Select transforms ('dct', 'dst', 'hadamard', or anything that is listed by 'help wfilters'):
transform_2D_HT_name     = 'bior1.5'; %% transform used for the HT filt. of size N1 x N1
transform_2D_Wiener_name = 'dct';     %% transform used for the Wiener filt. of size N1_wiener x N1_wiener
transform_3rd_dim_name   = 'haar';    %% transform used in the 3-rd dim, the same for HT and Wiener filt.

%%%% Hard-thresholding (HT) parameters:
N1                  = 8;   %% N1 x N1 is the block size used for the hard-thresholding (HT) filtering
Nstep               = 3;   %% sliding step to process every next reference block
N2                  = 16;  %% maximum number of similar blocks (maximum size of the 3rd dimension of a 3D array)
Ns                  = 39;  %% length of the side of the search neighborhood for full-search block-matching (BM), must be odd
tau_match           = 3000;%% threshold for the block-distance (d-distance)
lambda_thr2D        = 0;   %% threshold parameter for the coarse initial denoising used in the d-distance measure
lambda_thr3D        = 2.7; %% threshold parameter for the hard-thresholding in 3D transform domain
beta                = 2.0; %% parameter of the 2D Kaiser window used in the reconstruction

%%%% Wiener filtering parameters:
N1_wiener           = 8;
Nstep_wiener        = 3;
N2_wiener           = 32;
Ns_wiener           = 39;
tau_match_wiener    = 400;
beta_wiener         = 2.0;

%%%% Block-matching parameters:
stepFS              = 1;  %% step that forces to switch to full-search BM, "1" implies always full-search
smallLN             = 'not used in np'; %% if stepFS > 1, then this specifies the size of the small local search neighb.
stepFSW             = 1;
smallLNW            = 'not used in np';
thrToIncStep        = 8;  %% used in the HT filtering to increase the sliding step in uniform regions

if strcmp(profile, 'lc') == 1,

    Nstep               = 6;
    Ns                  = 25;
    Nstep_wiener        = 5;
    N2_wiener           = 16;
    Ns_wiener           = 25;

    thrToIncStep        = 3;
    smallLN             = 3;
    stepFS              = 6*Nstep;
    smallLNW            = 2;
    stepFSW             = 5*Nstep_wiener;

end

% Profile 'vn' was proposed in
%  Y. Hou, C. Zhao, D. Yang, and Y. Cheng, 'Comment on "Image Denoising by Sparse 3D Transform-Domain
%  Collaborative Filtering"', accepted for publication, IEEE Trans. on Image Processing, July, 2010.
% as a better alternative to that initially proposed in [1] (which is currently in profile 'vn_old')
if (strcmp(profile, 'vn') == 1) | (sigma > 40),

    N2                  = 32;
    Nstep               = 4;

    N1_wiener           = 11;
    Nstep_wiener        = 6;

    lambda_thr3D        = 2.8;
    thrToIncStep        = 3;
    tau_match_wiener    = 3500;
    tau_match           = 25000;

    Ns_wiener           = 39;

end

% The 'vn_old' profile corresponds to the original parameters for strong noise proposed in [1].
if (strcmp(profile, 'vn_old') == 1) & (sigma > 40),

    transform_2D_HT_name = 'dct';

    N1                  = 12;
    Nstep               = 4;

    N1_wiener           = 11;
    Nstep_wiener        = 6;

    lambda_thr3D        = 2.8;
    lambda_thr2D        = 2.0;
    thrToIncStep        = 3;
    tau_match_wiener    = 3500;
    tau_match           = 5000;

    Ns_wiener           = 39;

end

decLevel = 0;        %% dec. levels of the dyadic wavelet 2D transform for blocks (0 means full decomposition, higher values decrease the dec. number)
thr_mask = ones(N1); %% N1xN1 mask of threshold scaling coeff. --- by default there is no scaling, however the use of different thresholds for different wavelet decompoistion subbands can be done with this matrix

if strcmp(profile, 'high') == 1,

    decLevel     = 1;
    Nstep        = 2;
    Nstep_wiener = 2;
    lambda_thr3D = 2.5;
    vMask = ones(N1,1); vMask((end/4+1):end/2)= 1.01; vMask((end/2+1):end) = 1.07; %% this allows to have different threhsolds for the finest and next-to-the-finest subbands
    thr_mask = vMask * vMask';
    beta         = 2.5;
    beta_wiener  = 1.5;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Note: touch below this point only if you know what you are doing!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Check whether to dump information to the screen or reamin silent
dump_output_information = 1;
if (exist('print_to_screen') == 1) & (print_to_screen == 0),
    dump_output_information = 0;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Create transform matrices, etc.
%%%%
[Tfor, Tinv]   = getTransfMatrix(N1, transform_2D_HT_name, decLevel);  %% get (normalized) forward and inverse transform matrices
[TforW, TinvW] = getTransfMatrix(N1_wiener, transform_2D_Wiener_name); %% get (normalized) forward and inverse transform matrices

if (strcmp(transform_3rd_dim_name, 'haar') == 1) | (strcmp(transform_3rd_dim_name(end-2:end), '1.1') == 1),
    %%% If Haar is used in the 3-rd dimension, then a fast internal transform is used, thus no need to generate transform
    %%% matrices.
    hadper_trans_single_den         = {};
    inverse_hadper_trans_single_den = {};
else
    %%% Create transform matrices. The transforms are later applied by
    %%% matrix-vector multiplication for the 1D case.
    for hpow = 0:ceil(log2(max(N2,N2_wiener))),
        h = 2^hpow;
        [Tfor3rd, Tinv3rd]   = getTransfMatrix(h, transform_3rd_dim_name, 0);
        hadper_trans_single_den{h}         = single(Tfor3rd);
        inverse_hadper_trans_single_den{h} = single(Tinv3rd');
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 2D Kaiser windows used in the aggregation of block-wise estimates
%%%%
if beta_wiener==2 & beta==2 & N1_wiener==8 & N1==8 % hardcode the window function so that the signal processing toolbox is not needed by default
    Wwin2D = [ 0.1924    0.2989    0.3846    0.4325    0.4325    0.3846    0.2989    0.1924;
        0.2989    0.4642    0.5974    0.6717    0.6717    0.5974    0.4642    0.2989;
        0.3846    0.5974    0.7688    0.8644    0.8644    0.7688    0.5974    0.3846;
        0.4325    0.6717    0.8644    0.9718    0.9718    0.8644    0.6717    0.4325;
        0.4325    0.6717    0.8644    0.9718    0.9718    0.8644    0.6717    0.4325;
        0.3846    0.5974    0.7688    0.8644    0.8644    0.7688    0.5974    0.3846;
        0.2989    0.4642    0.5974    0.6717    0.6717    0.5974    0.4642    0.2989;
        0.1924    0.2989    0.3846    0.4325    0.4325    0.3846    0.2989    0.1924];
    Wwin2D_wiener = Wwin2D;
else
    Wwin2D           = kaiser(N1, beta) * kaiser(N1, beta)'; % Kaiser window used in the aggregation of the HT part
    Wwin2D_wiener    = kaiser(N1_wiener, beta_wiener) * kaiser(N1_wiener, beta_wiener)'; % Kaiser window used in the aggregation of the Wiener filt. part
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% If needed, read images, generate noise, or scale the images to the
%%%% [0,1] interval
%%%%
if (exist('yRGB') ~= 1) | (exist('zRGB') ~= 1)
    yRGB        = im2double(imread(image_name));  %% read a noise-free image
    randn('seed', 0);                          %% generate seed
    zRGB        = yRGB + (sigma/255)*randn(size(yRGB)); %% create a noisy image
else % external images
    image_name = 'External image';

    % convert zRGB to double precision
    zRGB = double(zRGB);

    % convert yRGB to double precision
    yRGB = double(yRGB);

    % if zRGB's range is [0, 255], then convert to [0, 1]
    if (max(zRGB(:)) > 10), % a naive check for intensity range
        zRGB = zRGB / 255;
    end

    % if yRGB's range is [0, 255], then convert to [0, 1]
    if (max(yRGB(:)) > 10), % a naive check for intensity range
        yRGB = yRGB / 255;
    end
end

if (size(zRGB,3) ~= 3) | (size(zRGB,4) ~= 1),
    error('CBM3D accepts only input RGB images (i.e. matrices of size M x N x 3).');
end

% Check if the true image yRGB is a valid one; if not, then we cannot compute PSNR, etc.
yRGB_is_invalid_image = (length(size(zRGB)) ~= length(size(yRGB))) | (size(zRGB,1) ~= size(yRGB,1)) | (size(zRGB,2) ~= size(yRGB,2)) | (size(zRGB,3) ~= size(yRGB,3));
if (yRGB_is_invalid_image),
    dump_output_information = 0;
end

[Xv, Xh, numSlices] = size(zRGB);              %%% obtain image sizes

if numSlices ~= 3
    fprintf('Error, an RGB color image is required!\n');
    return;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Change colorspace, compute the l2-norms of the new color channels
%%%%
[zColSpace l2normLumChrom] = function_rgb2LumChrom(zRGB, colorspace);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Print image information to the screen
%%%%
if dump_output_information == 1,
    fprintf(sprintf('Image: %s (%dx%dx%d), sigma: %.1f\n', image_name, Xv, Xh, numSlices, sigma));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Step 1. Basic estimate by collaborative hard-thresholding and using
%%%% the grouping constraint on the chrominances.
%%%%
tic;
y_hat = bm3d_thr_color(zColSpace, hadper_trans_single_den, Nstep, N1, N2, lambda_thr2D,...
    lambda_thr3D, tau_match*N1*N1/(255*255), (Ns-1)/2, sigma/255, thrToIncStep, single(Tfor), single(Tinv)', inverse_hadper_trans_single_den, single(thr_mask), 'unused arg', 'unused arg', l2normLumChrom, Wwin2D, smallLN, stepFS );
estimate_elapsed_time = toc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Step 2. Final estimate by collaborative Wiener filtering and using
%%%% the grouping constraint on the chrominances.
%%%%
tic;
yRGB_est = bm3d_wiener_color(zColSpace, y_hat, hadper_trans_single_den, Nstep_wiener, N1_wiener, N2_wiener, ...
    'unused_arg', tau_match_wiener*N1_wiener*N1_wiener/(255*255), (Ns_wiener-1)/2, sigma/255, 'unused arg', single(TforW), single(TinvW)', inverse_hadper_trans_single_den, 'unused arg', 'unused arg', l2normLumChrom, Wwin2D_wiener, smallLNW, stepFSW );
wiener_elapsed_time = toc;

yRGB_est = double(yRGB_est);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Convert back to RGB colorspace
%%%%
yRGB_est = function_LumChrom2rgb(yRGB_est, colorspace);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Calculate final estimate's PSNR and ISNR, print them, and show the
%%%% denoised image
%%%%
PSNR = 0; %% Remains 0 if the true image yRGB is not available
if (~yRGB_is_invalid_image), % then we assume yRGB is a valid image
    PSNR = 10*log10(1/mean((yRGB(:)-yRGB_est(:)).^2));
end

if dump_output_information == 1,
    fprintf(sprintf('FINAL ESTIMATE (total time: %.1f sec), PSNR: %.2f dB\n', ...
        wiener_elapsed_time + estimate_elapsed_time, PSNR));

    figure, imshow(min(max(zRGB,0),1)); title(sprintf('Noisy %s, PSNR: %.3f dB (sigma: %d)', ...
        image_name(1:end-4), 10*log10(1/mean((yRGB(:)-zRGB(:)).^2)), sigma));

    figure, imshow(min(max(yRGB_est,0),1)); title(sprintf('Denoised %s, PSNR: %.3f dB', ...
        image_name(1:end-4), PSNR));
end

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Some auxiliary functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Tforward, Tinverse] = getTransfMatrix (N, transform_type, dec_levels)
%
% Create forward and inverse transform matrices, which allow for perfect
% reconstruction. The forward transform matrix is normalized so that the
% l2-norm of each basis element is 1.
%
% [Tforward, Tinverse] = getTransfMatrix (N, transform_type, dec_levels)
%
%  INPUTS:
%
%   N               --> Size of the transform (for wavelets, must be 2^K)
%
%   transform_type  --> 'dct', 'dst', 'hadamard', or anything that is
%                       listed by 'help wfilters' (bi-orthogonal wavelets)
%                       'DCrand' -- an orthonormal transform with a DC and all
%                       the other basis elements of random nature
%
%   dec_levels      --> If a wavelet transform is generated, this is the
%                       desired decomposition level. Must be in the
%                       range [0, log2(N)-1], where "0" implies
%                       full decomposition.
%
%  OUTPUTS:
%
%   Tforward        --> (N x N) Forward transform matrix
%
%   Tinverse        --> (N x N) Inverse transform matrix
%

if exist('dec_levels') ~= 1,
    dec_levels = 0;
end

if N == 1,
    Tforward = 1;
elseif strcmp(transform_type, 'hadamard') == 1,
    Tforward    = hadamard(N);
elseif (N == 8) & strcmp(transform_type, 'bior1.5')==1 % hardcoded transform so that the wavelet toolbox is not needed to generate it
    Tforward =  [ 0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274;
       0.219417649252501   0.449283757993216   0.449283757993216   0.219417649252501  -0.219417649252501  -0.449283757993216  -0.449283757993216  -0.219417649252501;
       0.569359398342846   0.402347308162278  -0.402347308162278  -0.569359398342846  -0.083506045090284   0.083506045090284  -0.083506045090284   0.083506045090284;
      -0.083506045090284   0.083506045090284  -0.083506045090284   0.083506045090284   0.569359398342846   0.402347308162278  -0.402347308162278  -0.569359398342846;
       0.707106781186547  -0.707106781186547                   0                   0                   0                   0                   0                   0;
                       0                   0   0.707106781186547  -0.707106781186547                   0                   0                   0                   0;
                       0                   0                   0                   0   0.707106781186547  -0.707106781186547                   0                   0;
                       0                   0                   0                   0                   0                   0   0.707106781186547  -0.707106781186547];
elseif (N == 8) & strcmp(transform_type, 'dct')==1 % hardcoded transform so that the signal processing toolbox is not needed to generate it
    Tforward = [ 0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274;
       0.490392640201615   0.415734806151273   0.277785116509801   0.097545161008064  -0.097545161008064  -0.277785116509801  -0.415734806151273  -0.490392640201615;
       0.461939766255643   0.191341716182545  -0.191341716182545  -0.461939766255643  -0.461939766255643  -0.191341716182545   0.191341716182545   0.461939766255643;
       0.415734806151273  -0.097545161008064  -0.490392640201615  -0.277785116509801   0.277785116509801   0.490392640201615   0.097545161008064  -0.415734806151273;
       0.353553390593274  -0.353553390593274  -0.353553390593274   0.353553390593274   0.353553390593274  -0.353553390593274  -0.353553390593274   0.353553390593274;
       0.277785116509801  -0.490392640201615   0.097545161008064   0.415734806151273  -0.415734806151273  -0.097545161008064   0.490392640201615  -0.277785116509801;
       0.191341716182545  -0.461939766255643   0.461939766255643  -0.191341716182545  -0.191341716182545   0.461939766255643  -0.461939766255643   0.191341716182545;
       0.097545161008064  -0.277785116509801   0.415734806151273  -0.490392640201615   0.490392640201615  -0.415734806151273   0.277785116509801  -0.097545161008064];
elseif (N == 8) & strcmp(transform_type, 'dst')==1 % hardcoded transform so that the PDE toolbox is not needed to generate it
    Tforward = [ 0.161229841765317   0.303012985114696   0.408248290463863   0.464242826880013   0.464242826880013   0.408248290463863   0.303012985114696   0.161229841765317;
       0.303012985114696   0.464242826880013   0.408248290463863   0.161229841765317  -0.161229841765317  -0.408248290463863  -0.464242826880013  -0.303012985114696;
       0.408248290463863   0.408248290463863                   0  -0.408248290463863  -0.408248290463863                   0   0.408248290463863   0.408248290463863;
       0.464242826880013   0.161229841765317  -0.408248290463863  -0.303012985114696   0.303012985114696   0.408248290463863  -0.161229841765317  -0.464242826880013;
       0.464242826880013  -0.161229841765317  -0.408248290463863   0.303012985114696   0.303012985114696  -0.408248290463863  -0.161229841765317   0.464242826880013;
       0.408248290463863  -0.408248290463863                   0   0.408248290463863  -0.408248290463863                   0   0.408248290463863  -0.408248290463863;
       0.303012985114696  -0.464242826880013   0.408248290463863  -0.161229841765317  -0.161229841765317   0.408248290463863  -0.464242826880013   0.303012985114696;
       0.161229841765317  -0.303012985114696   0.408248290463863  -0.464242826880013   0.464242826880013  -0.408248290463863   0.303012985114696  -0.161229841765317];
elseif strcmp(transform_type, 'dct') == 1,
    Tforward    = dct(eye(N));
elseif strcmp(transform_type, 'dst') == 1,
    Tforward    = dst(eye(N));
elseif strcmp(transform_type, 'DCrand') == 1,
    x = randn(N); x(1:end,1) = 1; [Q,R] = qr(x);
    if (Q(1) < 0),
        Q = -Q;
    end;
    Tforward = Q';
else %% a wavelet decomposition supported by 'wavedec'
    %%% Set periodic boundary conditions, to preserve bi-orthogonality
    dwtmode('per','nodisp');

    Tforward = zeros(N,N);
    for i = 1:N
        Tforward(:,i)=wavedec(circshift([1 zeros(1,N-1)],[dec_levels i-1]), log2(N), transform_type);  %% construct transform matrix
    end
end

%%% Normalize the basis elements
Tforward = (Tforward' * diag(sqrt(1./sum(Tforward.^2,2))))';

%%% Compute the inverse transform matrix
Tinverse = inv(Tforward);

return;

function [y, A, l2normLumChrom]=function_rgb2LumChrom(xRGB, colormode)
% Forward color-space transformation   ( inverse transformation is function_LumChrom2rgb.m )
%
% Alessandro Foi - Tampere University of Technology - 2005 - 2006   Public release v1.03 (March 2006)
% -----------------------------------------------------------------------------------------------------------------------------------------------
%
% SYNTAX:
%
%   [y A l2normLumChrom] = function_rgb2LumChrom(xRGB, colormode);
%
% INPUTS:
%   xRGB  is RGB image with range [0 1]^3
%
%   colormode = 'opp', 'yCbCr', 'pca', or a custom 3x3 matrix
%
%       'opp'     Opponent color space ('opp' is equirange version)
%       'yCbCr'   The standard yCbCr (e.g. for JPEG images)
%       'pca'     Principal components   (note that this transformation is renormalized to be equirange)
%
% OUTPUTS:
%   y  is color-transformed image (with range typically included in or equal to [0 1]^3, depending on the transformation matrix)
%
%   l2normLumChrom (optional) l2-norm of the transformation (useful for noise std calculation)
%   A  transformation matrix  (used necessarily if colormode='pca')
%
%   NOTES:  -  If only two outputs are used, then the second output is l2normLumChrom, unless colormode='pca';
%           -  'opp' is used by default if no colormode is specified.
%
%
% USAGE EXAMPLE FOR PCA TRANSFORMATION:
%  %%%%  -- forward color transformation --
%    if colormode=='pca'
%       [zLumChrom colormode] = function_rgb2LumChrom(zRGB,colormode); % 'colormode' is assigned a 3x3 transform matrix
%    else
%       zLumChrom = function_rgb2LumChrom(zRGB,colormode);
%    end
%
%  %%%% [ ... ]  Some processing  [ ... ]
%
%  %%%%  -- inverse color transformation --
%    zRGB=function_LumChrom2rgb(zLumChrom,colormode);
%

if nargin==1
    colormode='opp';
end
change_output=0;
if size(colormode)==[3 3]
    A=colormode;
    l2normLumChrom=sqrt(sum(A.^2,2));
else
    if strcmp(colormode,'opp')
        A=[1/3 1/3 1/3; 0.5  0  -0.5; 0.25  -0.5  0.25];
    end
    if strcmp(colormode,'yCbCr')
        A=[0.299   0.587   0.114;   -0.16873660714285  -0.33126339285715   0.5;   0.5  -0.4186875  -0.0813125];
    end
    if strcmp(colormode,'pca')
        A=princomp(reshape(xRGB,[size(xRGB,1)*size(xRGB,2) 3]))';
        A=A./repmat(sum(A.*(A>0),2)-sum(A.*(A<0),2),[1 3]);  %% ranges are normalized to unitary length;
    else
        if nargout==2
            change_output=1;
        end
    end
end

%%%% Make sure that each channel's intensity range is [0,1]
maxV = sum(A.*(A>0),2);
minV = sum(A.*(A<0),2);
yNormal = (reshape(xRGB,[size(xRGB,1)*size(xRGB,2) 3]) * A' - repmat(minV, [1 size(xRGB,1)*size(xRGB,2)])') * diag(1./(maxV-minV)); % put in range [0,1]
y = reshape(yNormal, [size(xRGB,1) size(xRGB,2) 3]);

%%%% The l2-norm of each of the 3 transform basis elements
l2normLumChrom = diag(1./(maxV-minV))*sqrt(sum(A.^2,2));

if change_output
    A=l2normLumChrom;
end

return;

function yRGB=function_LumChrom2rgb(x,colormode)
% Inverse color-space transformation   ( forward transformation is function_rgb2LumChrom.m )
%
% Alessandro Foi - Tampere University of Technology - 2005 - 2006   Public release v1.03 (March 2006)
% -----------------------------------------------------------------------------------------------------------------------------------------------
%
% SYNTAX:
%
%   yRGB = function_LumChrom2rgb(x,colormode);
%
% INPUTS:
%  x  is color-transformed image (with range typically included in or equal to [0 1]^3, depending on the transformation matrix)
%
%  colormode = 'opp', 'yCbCr', or a custom 3x3 matrix (e.g. provided by the forward transform when 'pca' is selected)
%
%       'opp'      opponent color space ('opp' is equirange version)
%       'yCbCr'    standard yCbCr (e.g. for JPEG images)
%
% OUTPUTS:
%   x  is RGB image (with range [0 1]^3)
%
%
% NOTE:    'opp' is used by default if no colormode is specified
%

if nargin==1
    colormode='opp';
end
if size(colormode)==[3 3]
    A=colormode;
    B=inv(A);
else
    if strcmp(colormode,'opp')
        A =[1/3 1/3 1/3; 0.5  0  -0.5; 0.25  -0.5  0.25];
        B =[1 1 2/3;1 0 -4/3;1 -1 2/3];
    end
    if strcmp(colormode,'yCbCr')
        A=[0.299   0.587   0.114;   -0.16873660714285  -0.33126339285715   0.5;   0.5  -0.4186875  -0.0813125];
        B=inv(A);
    end
end

%%%% Make sure that each channel's intensity range is [0,1]
maxV = sum(A.*(A>0),2);
minV = sum(A.*(A<0),2);
xNormal = reshape(x,[size(x,1)*size(x,2) 3]) * diag(maxV-minV) +  repmat(minV, [1 size(x,1)*size(x,2)])'; % put in range [0,1]
yRGB = reshape(xNormal * B', [ size(x,1) size(x,2) 3]);

return;

author: xialeiliu

repo: RankIQA

save_path: github-repos/MATLAB/xialeiliu-RankIQA


text:
function [hd, ind1, ind2] = hausdorffDistance(pts1, pts2)
%HAUSDORFFDISTANCE  Hausdorff distance between two point sets
%
%   HD = hausdorffDistance(PTS1, PTS2)
%   Computes the Hausdorff distance between the two point sets PTS1 and
%   PTS2. The Hausdorf distance can be used to compare two shapes.
%
%   The distance between a point x and a set Y is given by:
%     d(x, Y) = inf { d(x,y) | y in Y }
%   The distance between two non empty sets X and Y is given by:
%     d(X, Y) = sup { d(x,Y) | x in X }
%   The Hausdorff distance between sets X and Y distance is defined as the
%   maximum of d(X,Y) and d(Y,X):
%     HD(X,Y) = max { d(X,Y), d(Y,X) }
%
%
%   Example
%   % Compute Hausdorff distance between an ellipse and a rectangle
%     % first define two shapes
%     rect = resamplePolygon(orientedBoxToPolygon([20 30 80 40 30]), 60);
%     poly = ellipseToPolygon([20 30 40 20 30], 500);
%     % display the shapes
%     figure; hold on
%     drawPolygon(poly, 'b');
%     drawPolygon(rect, 'g');
%     axis equal;
%     % compute hausdorff distance
%     [hd ind1 ind2] = hausdorffDistance(poly, rect);
%     p1h = poly(ind1, :);
%     p2h = rect(ind2, :);
%     drawPoint([p1h;p2h], 'mo');
%     drawEdge([p1h p2h], 'm')
%
%   See also
%   minDistancePoints
%
%   References
%   http://en.wikipedia.org/wiki/Hausdorff_distance
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2012-05-04,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2012 INRA - Cepia Software Platform.

% distance from pts1 to pts2
[dists1, ind12] = minDistancePoints(pts1, pts2);
[max1, ind11] = max(dists1);

% distance from pts2 to pts1
[dists2, ind22] = minDistancePoints(pts2, pts1);
[max2, ind21] = max(dists2);

% keep the max of the two distances
hd = max(max1, max2);

% keep the rigt indices
if max1 > max2
    ind1 = ind11;
    ind2 = ind12(ind11);
else
    ind1 = ind22(ind21);
    ind2 = ind21;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function estimateClinicalRF_Test(pathRF,nSplit,nBoot,seed)

startpath = pwd;

% INITIALIZATIONS
cd(pathRF)
training = load('training'); training = struct2cell(training); training = training{1};
testing = load('testing'); testing = struct2cell(testing); testing = testing{1};
nameOutcomes = fieldnames(training.outcomes); nOutcomes = numel(nameOutcomes);
fSetNames = fieldnames(training.textures.(nameOutcomes{1})); nFset = numel(fSetNames);
nClinical = size(training.clinical.table,2);
clinicNames = training.clinical.table.Properties.VariableNames';
testSplit = 1/3; % Proportion of test cases in the splits

% COMPUTATIONS
for o = 1:nOutcomes
    outcome = training.outcomes.(nameOutcomes{o});
    rng(seed), [trainSets,testSets] = getSplits(outcome,nSplit,testSplit);
    for f = 1:nFset
        text = training.textures.(nameOutcomes{o}).(fSetNames{f}); nText = size(text,2);
        indLeft = 1:nClinical;
        indChosen = [];
        for t = 1:nClinical
            maxMetric = 0;
            for i = 1:(nClinical-t+1)
                tableComb = [text,training.clinical.table(:,[indChosen,indLeft(i)])];
                cat = logical([zeros(1,nText),training.clinical.categories([indChosen,indLeft(i)])]);
                meanAUC = 0; meanSens = 0; meanSpec = 0;
                for s = 1:nSplit
                    Ttrain = tableComb(trainSets(:,s),:); Ttest = tableComb(testSets(:,s),:);
                    Ytrain = outcome(trainSets(:,s)); Ytest = outcome(testSets(:,s));
                    rng(seed), [RF] = trainRF_table(Ttrain,Ytrain,cat,nBoot);
                    [prob] = predictRF(Ttest,RF);
                    [aucSplit,sensSplit,specSplit,~] = calcPerformMetrics(prob,Ytest,0.5);
                    meanAUC = meanAUC + aucSplit;
                    meanSens = meanSens + sensSplit;
                    meanSpec = meanSpec + specSplit;
                end
                meanAUC = meanAUC/nSplit; meanSens = meanSens/nSplit; meanSpec = meanSpec/nSplit;
                %metric = 0.5*meanAUC + 0.5*(1-abs(meanSens - meanSpec));
                metric = meanAUC;
                if metric >= maxMetric
                    choice = indLeft(i);
                    indBye = i;
                    maxMetric = metric;
                end
            end
            indChosen = [indChosen,choice];
            var = clinicNames(indChosen);
            fprintf('\n  * Optimal added clinical variables for order %u, "%s", "%s" are: ',t,nameOutcomes{o},fSetNames{f});
            for v = 1:numel(var)-1
                fprintf([var{v},', '])
            end
            fprintf([var{end}])
            testing.clinical.bestAdd.(nameOutcomes{o}).(fSetNames{f}).(['Order',num2str(t)]).Indexes = indChosen;
            testing.clinical.bestAdd.(nameOutcomes{o}).(fSetNames{f}).(['Order',num2str(t)]).AUC = maxMetric;
            indLeft(indBye) = [];
        end
    end
end
save('testing','testing')
fprintf('\n')

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function [ BaseN,Gray ] = GrayNumberConvertion(Value,Base,Digits )

BaseN = zeros( 1,Digits+1 );
Gray  = zeros( 1,Digits+1 );

for i=1:Digits
    BaseN(i) =  mod(( Value / Base^(i-1)),Base );
    Value = Value - BaseN(i);
end

Shift = 0;
for i = (Digits+1):(-1):1
    Gray( i ) = mod( ( BaseN( i ) - Shift ),Base);
    Shift = Shift + Gray(i);
end

end

author: Sable

repo: mcbench-benchmarks


text:
%
% Smooth point-set registration method using neighboring constraints
% -------------------------------------------------------------------
%
% Authors: Gerard Sanromà, René Alquézar and Francesc Serratosa
%
% Contact: gsanorma@gmail.com
% Date: 15/02/2012
%
% Returns discrete {0,1} matrix with the highest coefficients in each row
% and column from matrix M
%
% Input
%   M: matrix of coefficients
%
% Output
%   R: discrete {0,1} output matrix with
%

function R = clean(M)

    M = M - min(M(:)); % translate so that the new min is zero
    fils = size(M,1);
    cols = size(M,2);
    R = zeros(size(M));

    while sum(M(:)) > 0
        [v i] = max(M(:));
        [f,c] = ind2sub(size(M),i);
        R(i) = 1;
        M(f,:) = zeros(1,cols);
        M(:,c) = zeros(fils,1);
    end


text:
function mapcolormap

load climcm_2
colormap(map);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [tree ndnm]=updatetree(tree,ndnm,al)
  apr=find([tree.sym]==al);
  if isempty(apr)
      apr=find([tree.wgt]==0);

      tree(end+1).left=0;
      tree(end).right=0;
      tree(end).sym=al;
      tree(end).parent=apr;
      tree(end).wgt=1;
      tree(end).nn=ndnm;

      ndnm=ndnm-1;

      tree(apr).wgt=tree(apr).wgt+1;
      tree(apr).right=length(tree);

      tree(end+1).left=0;
      tree(end).right=0;
      tree(end).sym=0;
      tree(end).parent=apr;
      tree(end).wgt=0;
      tree(end).nn=ndnm;

      ndnm=ndnm-1;
      tree(apr).left=length(tree);

      apr=tree(apr).parent;
      if apr==0
		return;
      end
  end
      while(1)
      ndwt=tree(apr).wgt;
      ndn=tree(apr).nn;
      wt=find([tree.wgt]==ndwt);
      mx=max([tree(wt).nn]);
      idx=find([tree.nn]==mx);
      if (ndn==mx) | (tree(apr).parent==idx)
          tree(apr).wgt=tree(apr).wgt+1;
      else
          temp=tree(idx);
          tree(idx)=tree(apr);
          tree(apr)=temp;

          temp=idx;
          idx=apr;
          apr=temp;

          temp=tree(idx).nn;
          tree(idx).nn=tree(apr).nn;
          tree(apr).nn=temp;

          temp=tree(idx).parent;
          tree(idx).parent=tree(apr).parent;
          tree(apr).parent=temp;

          id=tree(apr).left;
          if id~=0
              tree(id).parent=apr;
          end
          id=tree(apr).right;
          if id~=0
              tree(id).parent=apr;
          end

          id=tree(idx).left;
          if id~=0
              tree(id).parent=idx;
          end
          id=tree(idx).right;
          if id~=0
              tree(id).parent=idx;
          end

          tree(apr).wgt=tree(apr).wgt+1;
      end
          if tree(apr).parent==0
              break;
          end
          apr=tree(apr).parent;

      end
end


text:
function hess = cstepHess(fun,x,isTril,h)
%CSTEPHESS  Complex Step Hessian of a Scalar Function
%
%  cstepHess uses a complex step approach to finite difference to ensure
%  robust estimate of the Hessian. Note your function MUST be
%  written in Matlab code only and must not contain any if/else statements.
%
%  If it does not work - ensure you are using .' (transpose) and not ' (complex conjugate).
%
%   hess = cstepHess(fun,x) calculates the Hessian of fun at the point x.
%
%   hess = cstepHess(fun,x,isTril) specifies if the returned Hessian
%   should be Symmetric Lower Triangular.
%
%   hess = cstepHess(fun,x,isTril,h) specifies the step-size. This defaults to
%   1e-3.

%   Copyright (C) 2013 Jonathan Currie (IPL)
%
%   This code follows ideas from "New Complex-Step Derivative Approximations
%   with Application to Second-Order Kalman Filtering", by Kok-Lam Lai, John
%   L. Crassidis and Yang Cheng.

if(nargin < 4), h = 1e-3; end
if(nargin < 3 || isempty(isTril)), isTril = false; end

%Constants from paper
I = sqrt(2)/2*(1i + 1);
n = length(x); hess = zeros(n,n);
%Diagonal Elements
for k = 1:n
    xu = x; xu(k) = xu(k) + I*h;
    xl = x; xl(k) = xl(k) - I*h;
    hess(k,k) = imag((fun(xu) + fun(xl))/h^2);
end
%Lower Triangular Elements
lam = 1; k = n-1;
while(k > 0)
    for phi = 1:k
        xu = x; xu(phi:phi+lam) = xu(phi:phi+lam) + I*h;
        xl = x; xl(phi:phi+lam) = xl(phi:phi+lam) - I*h;
        Fsum = sum(sum(hess(phi:phi+lam,phi:phi+lam)));
        hess(phi+lam,phi) = (imag((fun(xu) + fun(xl))/h^2) - Fsum)/2;
        hess(phi,phi+lam) = hess(phi+lam,phi); %copy to upper tri
    end
    k = k - 1;
    lam = lam + 1;
end
%Make tril if requested
if(isTril), hess = tril(hess); end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
report_this_filefun(mfilename('fullpath'));

[ c, ia, ib ] = intersect(a(:,10), aq);

author: CelsoReyes


text:
function [wmProb,faErr,mdErr,eigVal] = dtiFindWhiteMatter(eigVal,b0,xformToAcpc)
%
%  wmProb = dtiFindWhiteMatter(dt6_or_eigVals,b0,[xformToAcpc])
%
% Example usage:
% b0=niftiRead('dti06/bin/b0.nii.gz');
% b0 = double(b0.data);
% dt6=niftiRead('dti06/bin/tensors.nii.gz');
% dt6 = double(squeeze(dt6.data(:,:,:,1,[1 3 6 2 4 5])));
% wmProb = dtiFindWhiteMatter(dt6,b0);
% dtiWriteNiftiWrapper(uint8(round(wmProb*255)),xformToAcpc,'dti06/bin/wmProb.nii.gz',1./255)
%
% HISTORY:
% 2007.07.24 RFD: wrote it.

% This var isn't used any more, but many callers pass it in.
if(~exist('xformToAcpc','var')), xformToAcpc = []; end

finalThresh = 0.5;
histClipB0 = false;

targetMd = [0.6 1.0]; % in micrometers^2/msec
targetFa = [0.4 1.0];

if(~isempty(b0))
   if(~isa(b0,'double'))
      b0 = double(b0);
   end
   % Make sure the b0 is histogram clipped
   if(max(b0(:))>1)
      b0 = mrAnatHistogramClip(b0,0.4,0.99);
   end
end
%templateFile = which('whiteMatter.nii.gz');

% Ensure diffusivity units are in micrometers^2/msec
[curUnitStr,scale] = dtiGuessDiffusivityUnits(eigVal);
eigVal = eigVal.*scale;
if(size(eigVal,4)==6)
    [eigVec,eigVal] = dtiEig(eigVal);
    clear eigVec;
end

%% Compute FA Error term
[fa,md] = dtiComputeFA(eigVal);
faErr = zeros(size(fa));
faErr(fa<targetFa(1)) = (targetFa(1)-fa(fa<targetFa(1)))./targetFa(1);
faErr(fa>targetFa(2)) = (fa(fa>targetFa(2))-targetFa(2))./targetFa(2);
faErr(faErr>1) = 1;
faErr(isnan(faErr)) = 1;

% things with bogus fa>=1 values are always bad
%faErr(fa>=1) = 1;

%% Compute MD Error term
% TODO: ensure md is in micrometers^2/msec!
mdErr = zeros(size(md));
mdErr(md<targetMd(1)) = (targetMd(1)-md(md<targetMd(1)))./targetMd(1);
mdErr(md>targetMd(2)) = (md(md>targetMd(2))-targetMd(2))./targetMd(2);
mdErr(mdErr>1) = 1;
mdErr(isnan(mdErr)) = 1;

%% Compute location prior from a template
%locPriorNi = niftiRead(templateFile);
sz = size(faErr);
%xf = locPriorNi.qto_ijk*xform;
%[sampZ,sampY] = meshgrid([1:sz(2)],[1:sz(1)]);
%sampY = sampY(:); sampZ = sampZ(:);
%coords = mrAnatXformCoords(xf,[ones(size(sampY)) sampY sampZ]);
%locPrior = dtiSMooth3(locPrior,7);
% TODO: get a location prior.
locPrior = ones(sz);

if(~isempty(b0))
    %% Compute the b0 error.
    % Since the b0 intensity scale is arbitrary, we'll use the info we have so
    % far to get an empirical estimate of the desired b0 range. We assume that
    % the b0 within the target region is roughly uniform, a safe assumption
    % for all white matter regions.
    %
    % The b0 adds an important check against artifacts, since certain common
    % anatomical features (such as large sinuses) create very low b0 values,
    % but the FA is high and MD is often within the normal tissue range.
    %
    % Get an initial estimate based on current maps
    img = (1-faErr).*(1-mdErr).*locPrior;
    curGuess = img>finalThresh*0.95;

    % We want to be sure to remove the junk outside the brain. We do
    % that with a hard mask based on normalized b0 intensity.
    brainMask = dtiCleanImageMask(mrAnatHistogramClip(b0,0.4,0.99)>0.4);
    curGuess = curGuess&brainMask;

    mnB0 = mode(b0(curGuess(:)));
    stdB0 = std(b0(curGuess(:)));
    % b0Err is essentially a z-score clipped at ~4SDs and normalized to the 0-1
    % range. We clip and normalize symmetrically (with 'abs').
    b0Err = abs((b0-mnB0)./stdB0);
    if(histClipB0)
       b0Err(b0Err>7) = 7;
       b0Err = b0Err-5;
       b0Err(b0Err<0) = 0;
       b0Err = b0Err./2;
    else
       b0Err(b0Err>5) = 5;
       b0Err = b0Err-1;
       b0Err(b0Err<0) = 0;
       b0Err = b0Err./4;
    end
else
    b0Err = 0;
end

%% Compute the final score
wmProb = (1-faErr).*(1-mdErr).*(1-b0Err).*locPrior;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function lfmaKernDisplay(varargin)

% LFMAKERNDISPLAY Display parameters of the LFMA kernel.
% FORMAT
% DESC displays the parameters of the LFMA
% kernel and the kernel type to the console.
% ARG kern : the kernel to display.
%
% FORMAT does the same as above, but indents the display according
% to the amount specified.
% ARG kern : the kernel to display.
% ARG spacing : how many spaces to indent the display of the kernel by.
%
% SEEALSO : lfmKernParamInit, modelDisplay, kernDisplay
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

lfmKernDisplay(varargin{:});

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [x,y,z] = doseCOM(doseStruct)
%"doseCOM"
%   Find [x,y,z] center of mass of dose contained in doseStruct.  Non
%   uniform slice width is taken into consideration.
%
%   doseStruct is planC{indexS.dose}(doseNum) where doseNum is the number
%   of the desired dose center of mass.
%
%JRA 3/12/04
%
%Usage:
%   function [x,y,z] = doseCOM(doseStruct)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

[xVals, yVals, zVals] = getDoseXYZVals(doseStruct);

dA = getDoseArray(doseStruct);

[xMesh, yMesh] = meshgrid(xVals, yVals);

for i=1:size(dA, 3);
    slice = dA(:,:,i);
    %Weighted mesh.
    wX = xMesh .* slice;
    wY = yMesh .* slice;
    sliceDose(i) = sum(slice(:));
    if sliceDose(i) == 0
        x(i) = 0;
        y(i) = 0;
    else
        x(i) = sum(wX(:)) / sliceDose(i);
        y(i) = sum(wY(:)) / sliceDose(i);
    end
end

%Get voxel thickness at each slice.
zDiff = diff(zVals);
dividers = [zVals(1)-zDiff(1)/2 zDiff/2 + zVals(1:end-1) zVals(end)+zDiff(end)/2];
voxThickness = diff(dividers);
totalHeight = sum(voxThickness);

totalDose = sum(sliceDose);
%Weight x,y,z values by dose on that slice/voxelThickness.
x = sum(x     .* sliceDose / voxThickness) / totalDose * totalHeight;
y = sum(y     .* sliceDose / voxThickness) / totalDose * totalHeight;
z = sum(zVals .* sliceDose / voxThickness) / totalDose * totalHeight;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function K = rbfXnoneKernCompute(rbfKern, noneKern, x1, x2)

% RBFXNONEKERNCOMPUTE Compute a cross kernel between RBF and NONE kernels.
% FORMAT
% DESC computes cross kernel terms between an RBF kernel and a dummy kernel for
% the multiple output kernel.
% ARG rbfKern : the kernel structure associated with the first
% rbf kernel.
% ARG noneKern : the kernel structure associated with the dummy kernel.
% ARG x : inputs for which kernel is to be computed.
% RETURN k : block of values from kernel matrix.
%
% FORMAT
% DESC computes cross kernel terms between two RBF kernels for
% the multiple output kernel.
% ARG rbfKern : the kernel structure associated with the
% rbf  kernel.
% ARG noneKern : the kernel structure associated with the dummy kernel.
% ARG x1 : row inputs for which kernel is to be computed.
% ARG x2 : column inputs for which kernel is to be computed.
% RETURN k : block of values from kernel matrix.
%
% SEEALSO : multiKernParamInit, multiKernCompute, rbfKernParamInit
%
% COPYRIGHT : Neil D. Lawrence, 2008

% KERN

if nargin < 4
  x2 = x1;
end

K = zeros(size(x1, 1), size(x2, 1));


text:
% Test file for trigtech/any.m.

function pass = test_any(pref)

if ( nargin < 1 )
    pref = trigtech.techPref();
end

testclass = trigtech();

% Check behavior for any() down columns.
pass(1) = ~any(testclass);

f = testclass.make(@(x) [sin(pi*x) 0*x cos(pi*x)]);
pass(2) = isequal(any(f), [1 0 1]);

% Check behavior for any() across rows.
g = any(f, 2);
pass(3) = isequal(g.coeffs, 1);

f = testclass.make(@(x) [0*x 0*x]);
g = any(f, 2);
pass(4) = isequal(g.coeffs, 0);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [tt]=tt_elem_reverse(tt)
% Elementwise reverse of a tensor train vector
% NOT tested on a multidimensional case
%
% use tt_qutrtoepl(core(tt_shf(d)'*tt_elem_reverse(tt))) - to build
% upper-toeplitz matrix in "normal", blin, indexing, i.e. second element
% goes to S^{+1}, usw.

d = tt.d;
tt = core2cell(tt);

for i=1:d
    n = size(tt{i}, 2);
    tt{i} = tt{i}(:,(n:-1:1),:);
end;

tt = cell2core(tt_tensor, tt);

end

author: oseledets

repo: TT-Toolbox


text:
% Gtomo_nufft_test.m
% Test the Gtomo_nufft object

% create Gtomo_nufft class object
if ~isvar('An'), printm 'setup'
	ig = image_geom('nx', 64, 'ny', 60, 'fov', 480);
	ig.mask = ig.circ > 0;
	if 1
		sg = sino_geom('par', 'nb', 70, 'na', 80, 'dr', 7, ...
			'orbit_start', -15, 'strip_width', 'd');
		fan_arg = {}; % parallel beam
	else
		sg = sino_geom('ge1', 'orbit_start', -15, 'down', 8, ...
			'strip_width', 'd');
		fan_arg = { % fan beam
			'offset_s', sg.offset_s, ...
			'dis_src_det', sg.dsd, ...
			'dis_iso_det', sg.dod, ...
		};
	end

	% test non-rect image basis function and beam shape (S. Matej)
	if 0
		basis.type  = 'KB';
		basis.diam  = 4;
		basis.shape = 10.4;
		basis.m     = 2;
		basis.dim   = 2;
		basis.kernel=[];

		% KB: FWHM~1.0 : J=6, alpha=40, m=2
		beam.type  = 'KB';
		beam.diam  = 6;
		beam.shape = 40.;
		beam.m     = 2;
	end

	cpu etic

%{
	arg_geom = { % for old version
		fan_arg{:}, ...
		'orbit', sg.orbit, ...
		'orbit_start', sg.orbit_start, ...
		'dx', ig.dx, ...
		'ds', sg.d, ...
		'yscale', -ig.dy/ig.dx, ...
		'strip_width', sg.strip_width, ...
	};
%}

	arg = {
		'interp', {'table', 2^11, 'minmax:kb'}, ...
...%		'basis', basis, ...
...%		'beam', beam, ...
		'is.complex', 0, ...
...%		'kaiser', ...	% use KB interpolator
...%		'uniform', ...
...%		'bilin', ...
...%		'is.test', 1, ...
		'class', 'fatrix2', ...
	};

%	Ao = Gtomo_nufft(ig.mask, [sg.nb sg.na], arg_geom{:}, arg{:});
	An = Gtomo_nufft_new(sg, ig, arg{:});
%	jf_equal(struct(Ao), struct(An)), return

	cpu etoc 'An pre time:'
end

%
% make test image
%
if 0
%	x = ig.zeros;
%	x(12, 8) = 1;
%	x((nx/4+1):(3*nx/4),(ny/4+1):(3*ny/4)) = 1;
	x = ig.unitv(ig.nx/4,ig.ny/3);	% point source
elseif 1
	x = ellipse_im(ig, [0 0 [1 1]*ig.ny*ig.dx*0.9/2 0 1], 'oversample', 2);
	x(end/4:end/2,end/4:end/2) = 1.25;
else
	ix = [-(nx-1)/2:(nx-1)/2]/nx;
	iy = [-(ny-1)/2:(ny-1)/2]/ny;
	[ix iy] = ndgrid(ix, iy);
	x = exp(-(ix.^2 + iy.^2) / 10^2);
end

% check real / imaginary
if 1
	y = An * x(ig.mask);
	y = sg.shape(y);
	im plc 3 3
	im(1, x, 'x')
	im(2, real(y), 'real(An*x)'), cbar
	im(3, imag(y), 'imag(An*x)'), cbar
prompt
end

% DSFT version for comparison
if ~isvar('Ad'), printm 'Ad'
%	Ao = Gtomo_nufft(ig.mask, [sg.nb sg.na], arg_geom{:}, arg{:}, ...
%		'is.dsft2', 1, 'is.dsft1', 1);
	Ad = Gtomo_nufft_new(sg, ig, arg{:}, 'is.dsft2', 1, 'is.dsft1', 1);
%	jf_equal(struct(Ao), struct(Ad)), return
end

% create a "strip integral" system object
if ~isvar('As'), printm 'As'
%	Ao = aspire_pair(sg, ig, 'strip_width', An.arg.strip_width, ...
%		'support', 'all');
%		'support', ig.mask);
%	Ao = Gtomo2_dscmex(Ao, 'nthread', 2); % multiprocessor!
	As = Gtomo2_dscmex(sg, ig, 'nthread', 2); % multiprocessor!
%	jf_equal(struct(Ao), struct(As)), return
end

% compare back projectors
if 0
	y = sg.ones;
	y = sg.unitv(sg.nb/2+7, 9);
%	y = yd / 100;
	xn = An' * y;
	xd = Ad' * y;
	xs = As' * y;
	im(4, xn, 'An''*y'), cbar
	im(5, xd, 'Ad''*y'), cbar
	im(6, xs, 'As''*y'), cbar
	im(8, xn-xd, 'An''*y - Ad''*y'), cbar
	im(9, xn-xs, 'An''*y - As''*y'), cbar
%	printf('back nrms = %g%%', nrms(x2(:), x1(:)) * 100)
prompt
end

% compare projectors
if 1
	yn = An * x;
	yd = Ad * x;
	ys = As * x;
	im(4, yd, 'Ad*x'), cbar
	im(5, yn-yd, 'An*x-Ad*x'), cbar
	im(7, ys, 'As*x'), cbar
	im(8, yn-ys, 'An*x-As*x'), cbar
%	printf('forward nrms = %g%%', nrms(y2(:), y1(:)) * 100)
prompt
end

if 1, printm 'single ray'
	y = sg.unitv(sg.nb/2+9,1);
	im pl 2 3
	im(1, As'*y, 'back ray strip'), cbar h
	im(2, An'*y, 'back ray nufft'), cbar h
	im(3, Ad'*y, 'back ray dsft'), cbar h
	im(4, An'*y - Ad'*y, 'nufft - dsft'), cbar h
	im(5, An'*y - As'*y, 'nufft - strip'), cbar h
	im(6, Ad'*y - As'*y, 'dsft - strip'), cbar h
return
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [data] = bz_shuffleCircular(data)
% circularly shifts each row

for i=1:size(data,1)
    data(i,:) = circshift(data(i,:),randi(size(data,2)));
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [hz,hprimez] = cnnactivate(z,type)
%CNNACTIVATE Activation function for CNNs.
%   [HZ,HPRIMEZ] = CNNACTIVATE(Z,TYPE) outputs activation values, HZ =
%   h(z), and its derivative, HPRIMEZ = h'(z), evaluated at each element
%   of z, where z is a vector or matrix (see Sections 14.5 and 14.6 in
%   DIPUM3E). The type of activation function used is specified in
%   string TYPE, valid values of which are:
%
%   TYPE                            EQUATION IMPLEMENTED
%
%   'sigmoid' (the default)         h(z) = 1./(1 + exp(-z))
%                                   hprime(z) = h.*(1 - h)
%
%   'tanh'                          h(z) = tanh(z)
%                                   hprime(z) = 1 - h.^2
%
%   'ReLU'                          h(z) = max(0,z)
%                                   hprime(z) = 1 for z > 0
%                                             = 0 otherwise
%
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Set default.
if nargin == 1
    type = 'sigmoid';
end

% Set to lower case to reduce possible user case typo errors.
type = lower(type);

% Compute h and hprime as functions of z.
switch type
    case 'sigmoid'
        hz = 1./(1 + exp(-z));
        hprimez = hz.*(1 - hz);
    case 'tanh'
        hz = tanh(z);
        hprimez = 1 - hz.^2;
    case 'relu'
       % Original based on the definition of ReLU. Works well provided
       % that alpha is small enough.
        hz = max(0,z);
        hprimez = zeros(size(z));
        hprimez(z > 0) = 1;
end

author: dipum

repo: dipum-toolbox

save_path: github-repos/MATLAB/dipum-dipum-toolbox


text:
%kImpRaw 'imports a raw data file and bins it '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros ImpRaw.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'First Input data object'
% Integer: sizex 'SizeX', default: 256: 'x dimension of the data'
% Integer: binx 'BinX', default: 1: 'number of pixels to bin along x (if below 0 it is extracted, above 0 summed)'
% Toggle: bt 'Byte Type', default: 0: 'True: Is bytes; False: Is Short'
% Integer: sizey 'SizeY', default: 256: 'y dimension of the data'
% Integer: biny 'BinY', default: 1: 'number of pixels to bin along y (if below 0 it is extracted, above 0 summed)'
% Integer: sizez 'SizeZ', default: 10: 'z dimension of the data'
% Integer: binz 'BinZ', default: 1: 'number of pixels to bin along z (if below 0 it is extracted, above 0 summed)'
% Integer: sizet 'SizeT', default: 1: 'time dimension of the data'
% Integer: sizee 'SizeE', default: 1: 'element dimension of the data'
% Integer: offset 'offset', default: 0: 'Offset in bytes'
% OutputFile: o 'Output', required: 'Resulting output data object'
%
% Example: o = kImpRaw(i, {'i','';'sizex',256;'binx',1;'bt',0;'sizey',256;'biny',1;'sizez',10;'binz',1;'sizet',1;'sizee',1;'offset',0;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% ImpRaw - imports a raw data file and bins it
%
%  DESCRIPTION
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kImpRaw(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kImpRaw(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'sizex', 256;'binx', 1;'bt', 0;'sizey', 256;'biny', 1;'sizez', 10;'binz', 1;'sizet', 1;'sizee', 1;'offset', 0;'o', '__output'};
maxval={0,0,0,0,0,0,0,0,0,0,0,0};
minval={0,0,0,0,0,0,0,0,0,0,0,0};
istoggle=[0,0,0,1,0,0,0,0,0,0,0,0];
was_set=istoggle * 0;
paramtype={'InputFile','Integer','Integer','Toggle','Integer','Integer','Integer','Integer','Integer','Integer','Integer','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'impraw"  -k'],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function a = kahan_inverse ( alpha, n )

%*****************************************************************************80
%
%% KAHAN_INVERSE returns the inverse of the KAHAN matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA, the scalar that defines A.  A typical
%    value is 1.2.  The "interesting" range of ALPHA is 0 < ALPHA < PI.
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  ci = cos ( alpha );

  for i = 1 : n
    for j = 1 : n

      if ( i == j )
        a(i,j) = 1.0;
      elseif ( i == j - 1 )
        a(i,j) = ci;
      elseif ( i < j )
        a(i,j) = ci * ( 1.0 + ci )^(j-i-1);
      else
        a(i,j) = 0.0;
      end

    end
  end
%
%  Scale the columns.
%
  for j = 1 : n
    si = sin ( alpha )^j;
    a(1:n,j) = a(1:n,j) / si;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Question No:2

% Detect the edges in an image
% a)Using a vertical Sobel mask with a threshold determined automatically.
% b)Using a vertical Sobel mask with a specified threshold
% c)Both vertical and horizontal edges with a specified threshold
% d)Edges at 45 degree with image filter using a specified mask and a
%   specified threshold
% e)Edges at -45 degree with image filter using a specified mask and a
%   specified threshold

function sobeledge(x)
f=imread(x);
f=im2double(f);
choice=0;
while (choice~=6)
choice=input('1: Vertical Sobel Mask - Automatic Threshold\n2: Vertical Sobel Mask - Specified Threshold\n3: Horizontal and Vertical Edges - Specified Threshold\n4: Edges at 45 degree - Specified Threshold\n5: Edges at -45 degree - Specified Threshold\n6: Exit\n Enter your choice : ');
switch choice
    case 1
        [VSFAT Threshold]=edge(f,'sobel','vertical');
        figure, imshow(f),title('Original Image'),figure,imshow(VSFAT),title('Sobel Filter - Automatic Threshold');
        Threshold
    case 2
        Threshold=input('Enter a Threshold value : ');
        VSFST=edge(f,'sobel',Threshold,'vertical');
        figure, imshow(f),title('Original Image'),figure,imshow(VSFST),title('Sobel Filter - Specified Threshold');
        Threshold
    case 3
        Threshold=input('Enter a Threshold value : ');
        SFST=edge(f,'sobel',Threshold);
        figure, imshow(f),title('Original Image'),figure,imshow(SFST),title('Sobel Filter (Horizontal and Vertical) - Specified Threshold');
        Threshold
    case 4
        Threshold=input('Enter a Threshold value : ');
        s45=[-2 -1 0;-1 0 1;0 1 2];
        SFST45=imfilter(f,s45,'replicate');
        SFST45=SFST45>=Threshold;
        figure, imshow(f),title('Original Image'),figure,imshow(SFST45),title('Sobel Filter (45 Degree) - Specified Threshold');
        Threshold
    case 5
        Threshold=input('Enter a Threshold value : ');
        sm45=[0 1 2;-1 0 1;-2 -1 0];
        SFSTM45=imfilter(f,sm45,'replicate');
        SFSTM45=SFSTM45>=Threshold;
        figure, imshow(f),title('Original Image'),figure,imshow(SFSTM45),title('Sobel Filter (45 Degree) - Specified Threshold');
        Threshold
    case 6
        display('Program Exited');
    otherwise
        display('\nWrong Choice\n');
    end
end


text:
function output = callipqp(varargin)

% Retrieve needed data
interfacedata = varargin{1};
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
K       = interfacedata.K;
x0      = interfacedata.x0;
Q       = 2*interfacedata.Q;
ub      = interfacedata.ub;
lb      = interfacedata.lb;

showprogress('Calling IPQP',options.showprogress);

if K.f>0
    A = -F_struc(1:K.f,2:end);
    b = F_struc(1:K.f,1);
else
    A = [];
    b = [];
end

if K.l>0
C = -F_struc(1+K.f:K.f+K.l,2:end);
d = F_struc(1+K.f:K.f+K.l,1);
else
    C = [];
    d = [];
end

if ~isempty(ub)
    C = [C;eye(length(c))];
    d = [d;ub];
end
if ~isempty(lb)
    C = [C;-eye(length(c))];
    d = [d;-lb];
end

solvertime = tic;
[x,problem] = ipqp(2*Q,c,C,d,A,b);
solvertime = toc(solvertime);

% Internal format for duals
D_struc = [];

switch problem
case 0
    problem = 0;
    infostr = yalmiperror(problem,'IPQP');
case 1
    problem = 1;
    infostr = yalmiperror(problem,'IPQP');
otherwise

    problem = 2;
    infostr = yalmiperror(problem,'IPQP');
end

% Save all data sent to solver?
if options.savesolverinput
%     solverinput.A = A;
%     solverinput.b = b;
%     solverinput.C = Aq;
%     solverinput.q = beq;
%     solverinput.c = c;
%     solverinput.H = Q;
%     solverinput.options = options.nag;
else
    solverinput = [];
end

% Save all data from the solver?
if options.savesolveroutput
%     solveroutput.x = x;
%     solveroutput.fmin = fmin;
%     solveroutput.flag = flag;
%     solveroutput.output=output;
%     solveroutput.lambda=lambda;
else
    solveroutput = [];
end

% Standard interface
output = createOutputStructure(x,D_struc,[],problem,infostr,solverinput,solveroutput,solvertime);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function sixdof = hom2six(M)
% function sixdof = hom2six(M)
%
% Convert homogenous transformation (4x4 matrix) to 6-dof representation
% as used in SPM8. That is, this function assumes that
% M = Trans([x;y;z]) * RotX(a) * RotY(b) * RotZ(c)
% and returns the 6 parameters [x,y,z,a,b,c] that generate M.
%
% Note that this representation is not unique in case b = +/- pi/2.

% (C) 2010 S. Klanke

%Rotation part is Rx(a)*Ry(b)*Rz(c)
%[           cb*cc,           cb*sc,              sb]
%[ -sa*sb*cc-ca*sc, -sa*sb*sc+ca*cc,           sa*cb]
%[ -ca*sb*cc+sa*sc, -ca*sb*sc-sa*cc,           ca*cb]

if M(1,3) == 1   % sb = 1, cb = 0
%[  0, 0, 1]
%[ -sa*cc-ca*sc, -sa*sc+ca*cc, 0]
%[ -ca*cc+sa*sc, -ca*sc-sa*cc, 0]
% not unique -> set a = 0 - > ca=1, sa=0
%[  0   0  1]
%[-sc  cc  0]
%[-cc -sc  0]
	a = 0;
	b = pi/2;
	c = atan2(-M(3,2),M(2,2));
elseif M(1,3) == -1  % sb = -1, cb = 0
%[  0, 0, -1]
%[ +sa*cc-ca*sc, +sa*sc+ca*cc, 0]
%[ +ca*cc+sa*sc, +ca*sc-sa*cc, 0]
% not unique -> set a = 0 - > ca=1, sa=0
%[  0   0  -1]
%[-sc  cc   0]
%[ cc  sc   0]
    a = 0;
	b = -pi/2;
	c = atan2(M(3,2),M(2,2));
else
   % cos(b)~=0
   c = atan2(M(1,2),M(1,1));
   b = atan2(M(1,3),norm(M(1,1:2)));
   a = atan2(M(2,3),M(3,3));
end

sixdof = [M(1:3,4)' a b c];

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ a, b ] = p57_lim ( )

%*****************************************************************************80
%
%% P57_LIM returns the integration limits for problem 57.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, the limits of integration.
%
  a = 0.0;
  b = 1.0;

  return
end


text:
function msg = mhe_init()

s = 'http://byu.apmonitor.com';
b = 'mhe';

apm(s,b,'clear all');

% load model and data
apm_load(s,b,'model.apm');
csv_load(s,b,'mhe.csv');

% configure MV / CV
apm_info(s,b,'FV','U');
apm_info(s,b,'FV','tau');
apm_info(s,b,'FV','a1');
apm_info(s,b,'FV','a2');
apm_info(s,b,'MV','Q1');
apm_info(s,b,'MV','Q2');
apm_info(s,b,'SV','TH1');
apm_info(s,b,'SV','TH2');
apm_info(s,b,'CV','TC1');
apm_info(s,b,'CV','TC2');

% tune FVs
apm_option(s,b,'U.dmax',1.0);
apm_option(s,b,'U.lower',5);
apm_option(s,b,'U.upper',15);

apm_option(s,b,'tau.dmax',1.0);
apm_option(s,b,'tau.lower',4.0);
apm_option(s,b,'tau.upper',8.0);

apm_option(s,b,'a1.dmax',0.001);
apm_option(s,b,'a1.lower',0.003);
apm_option(s,b,'a1.upper',0.03);

apm_option(s,b,'a2.dmax',0.001);
apm_option(s,b,'a2.lower',0.002);
apm_option(s,b,'a2.upper',0.02);

% turn on FVs as degrees of freedom, later
apm_option(s,b,'U.status',0);
apm_option(s,b,'tau.status',0);
apm_option(s,b,'a1.status',0);
apm_option(s,b,'a2.status',0);

apm_option(s,b,'U.fstatus',0);
apm_option(s,b,'tau.fstatus',0);
apm_option(s,b,'a1.fstatus',0);
apm_option(s,b,'a2.fstatus',0);

% read Q, don't let optimize use MV
apm_option(s,b,'Q1.status',0);
apm_option(s,b,'Q1.fstatus',1);

apm_option(s,b,'Q2.status',0);
apm_option(s,b,'Q2.fstatus',1);

% include CV in objective function
apm_option(s,b,'TC1.status',1);
apm_option(s,b,'TC1.fstatus',1);
%apm_option(s,b,'TC1.meas_gap',0.5);

apm_option(s,b,'TC2.status',1);
apm_option(s,b,'TC2.fstatus',1);
%apm_option(s,b,'TC2.meas_gap',0.5);

% web-viewer option, update every second
apm_option(s,b,'apm.web_plot_freq',3);

% dynamic estimation
apm_option(s,b,'apm.imode',5);
apm_option(s,b,'apm.ev_type',2);
apm_option(s,b,'apm.nodes',3);
apm_option(s,b,'apm.solver',3);
apm_option(s,b,'apm.coldstart',1);

msg = 'initialization complete';

end

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function result=cpcr(x,y,varargin)

%CPCR performs a classical principal components regression.
% First, classical PCA is applied to the predictor variables x (see cpca.m) and
% k components are retained. Then a multiple linear regression method (see mlr.m)
% is performed of the response variable y on the k principal components.
%
% I/O: result=cpcr(x,y,'k',2);
%
% Required input arguments:
%      x : Data matrix of the explanatory variables
%          (n observations in rows, p variables in columns)
%      y : Data matrix of the response variables
%          (n observations in rows, q variables in columns)
%
% Optional input argument:
%      k : Number of principal components to compute. If k is missing,
%          a scree plot is drawn which allows to select
%          the number of principal components.
%  plots : If equal to one (default), a menu is shown which allows to draw several plots,
%          such as a score outlier map and a regression outlier map.
%          If 'plots' is equal to zero, all plots are suppressed.
%          See also makeplot.m
%
% The output is a structure containing
%
%   result.slope      : Classical slope
%   result.int        : Classical intercept
%   result.fitted     : Classcial prediction vector
%   result.res        : Classical residuals
%   result.sigma      : Estimated variance-covariance matrix of the residuals
%   result.rsquared   : R-squared value
%   result.k          : Number of components used in the regression
%   result.sd         : Classical score distances
%   result.od         : Classical orthogonal distances
%   result.resd       : Residual distances (when there are several response variables).
%                       If univariate regression is performed, it contains the standardized residuals.
%   result.cutoff     : Cutoff values for the score, orthogonal and residual distances.
%   result.flag      : The observations whose orthogonal distance is larger than result.cutoff.od
%                      (orthogonal outliers => result.flag.od) and/or whose residual distance is
%                      larger than result.cutoff.resd (bad leverage points/vertical outliers => result.flag.resd)
%                      can be considered as outliers and receive a flag equal to zero (=> result.flag.all).
%                      The regular observations, including the good leverage points, receive a flag 1.
%   result.class      : 'CPCR'
%   result.cpca       : Full output of the classical PCA part (see cpca.m)
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
%              http://wis.kuleuven.be/stat/robust.html
%
% Written by S. Verboven
% Last Update: 05/04/2003

default=struct('plots',1,'k',0);
list=fieldnames(default);
options=default;
IN=length(list);
i=1;
counter=1;
if nargin>2
    %
    % placing inputfields in array of strings
    %
    for j=1:nargin-2
        if rem(j,2)~=0
            chklist{i}=varargin{j};
            i=i+1;
        end
    end
    % Checking which default parameters have to be changed
    % and keep them in the structure 'options'.
    while counter<=IN
        index=strmatch(list(counter,:),chklist,'exact');
        if ~isempty(index) % in case of similarity
            for j=1:nargin-2 % searching the index of the accompanying field
                if rem(j,2)~=0 % fieldnames are placed on odd index
                    if strcmp(chklist{index},varargin{j})
                        I=j;
                    end
                end
            end
            options=setfield(options,chklist{index},varargin{I+1});
            index=[];
        end
        counter=counter+1;
    end
end

%classical PCA
if options.k==0
    out.cpca=cpca(x,'plots',0);
else
    out.cpca=cpca(x,'plots',0,'k',options.k);
end
%model with intercept
T=[out.cpca.T(:,1:out.cpca.k) ones(size(out.cpca.T,1),1)];
%Multivariate linear regression
out.a=inv(T'*T)*T'*y;
out.fitted=T*out.a;
len=size(out.a,1);
p=size(T,2);
q=size(y,2);
geg=[T,y];
[n,m]=size(geg);
S=cov(geg);
Sx=S(1:p-1,1:p-1);
Sxy=S(1:p-1,p+1:m);
Syx=Sxy';
Sy=S(p+1:m,p+1:m);
Se=Sy-out.a(1:p-1,1:q)'*Sx*out.a(1:p-1,1:q); %variance of errors
%regression coefficients slope and intercept ([\beta \alpha])
out.coeffs=[out.cpca.P(:,1:out.cpca.k)*out.a(1:len-1,:); out.a(len,:)-out.cpca.M*out.cpca.P(:,1:out.cpca.k)*out.a(1:len-1,:)];    %%coefficients in the original space;
out.slope=out.cpca.P(:,1:out.cpca.k)*out.a(1:len-1,:);
out.int=out.a(len,:)-out.cpca.M*out.cpca.P(:,1:out.cpca.k)*out.a(1:len-1,:);
out.res=y-out.fitted;
out.sigma=Se;
if q==1
    out.stdres=out.res./sqrt(diag(out.sigma));
end
out.k=out.cpca.k;
STTm=sum((y-repmat(mean(y),length(y),1)).^2);
SSE=sum(out.res.^2);
out.rsquared=1-SSE/STTm;
out.class='CPCR';

%calculating residual distances
if q>1
    if (-log(det(Se)/(m-1)))>50
        out.resd='singularity';
    else
       cen=zeros(q,1);
       out.resd=sqrt(libra_mahalanobis(out.res,cen','cov',Se))';
   end
else % q==1
    out.resd=out.stdres; %standardized residuals
end
out.cutoff=out.cpca.cutoff;
out.cutoff.resd=sqrt(chi2inv(0.975,size(y,2)));
%computing flags
out.flag.od=out.cpca.flag.od;
out.flag.sd=out.cpca.flag.sd;
out.flag.resd=abs(out.resd)<=out.cutoff.resd;
out.flag.all=(out.flag.od & out.flag.resd);

result=struct( 'slope',{out.slope}, 'int',{out.int},'fitted',{out.fitted},'res',{out.res},...
    'cov',{out.sigma},'rsquared',{out.rsquared},...
    'k',{out.k},'sd', {out.cpca.sd},'od',{out.cpca.od},'resd',{out.resd},...
    'cutoff',{out.cutoff},'flag',{out.flag},'class',{out.class},...
    'cpca',{out.cpca});

try
    if options.plots
        makeplot(result)
    end
catch %output must be given even if plots are interrupted
    %> delete(gcf) to get rid of the menu
    end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function cosmo_wizard_set_config()
% GUI-based 'wizard' to set CoSMoMVPA configuration file
%
% cosmo_wizard_set_config()
%
% This 'wizard' can be used to set the CoSMoMVPA configuration paths
% for exercises and demonstration outputs (in CoSMoMVPA's 'examples/'
% directory).
%
% It asks to select the directories for the 'tutorial_data_path'
% and 'output_data_path'. The former requires tutorial data, which can be
% downloaded from cosmomvpa.org/download.html
%
% It will store a '.cosmomvpa.cfg' file in either the user's home directory
% (on Linux / OSX), or in any available userpath() directory.
%
% Selecting cancel at any time will close this wizard.
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

% simulate go-to
while true
    if isempty(which('cosmo_config'))
        % store original path.
        % (should we ask to user whether the path should be set)
        orig_path=path();
        path_resetter=onCleanup(@()path(orig_path));

        addpath(get_cosmo_mvpa_dir());
        cosmo_set_path();
    end

    % read existing configuration, if it exists
    [config,fn]=get_cosmo_config();

    if isempty(fn)
        % no existing configuration, find location where to store it
        [fn,error_msg]=get_output_filename();

        if ~isempty(error_msg)
            break
        end
    end

    % ask for paths and update configuration
    [config,error_msg]=update_config(config);
    if ~isempty(error_msg)
        break;
    end

    % store configuration
    error_msg=write_config(fn,config);

    if ~isempty(error_msg)
        break;
    end

    % all done without error, quit this function
    msgbox('Done! Happy CoSMoMVPA-ing!');
    return
end

error('%s did not complete, reason: %s',mfilename(),error_msg);

function mvpa_dir=get_cosmo_mvpa_dir()
    cur_dir=fileparts(mfilename('fullpath'));
    mvpa_dir=fullfile(fileparts(cur_dir),'mvpa');

function error_msg=write_config(fn,config)
    config_str=config2str(config);

    title='Confirm writing CoSMoMVPA configuration';
    body=sprintf(['The following CoSMoMVPA configuration\n\n%s\n'...
                    'will be written in this file:\n\n%s\n\nContinue?'],...
                    config_str,fn);

    if show_ok_cancel_dialog(title,body)
        error_msg='';
        cosmo_config(fn,config);
        return;
    end

    error_msg='Writing configuration file was cancelled';

function str=config2str(config)
    rows=cellfun(@(key)sprintf('%s=%s\n',key,config.(key)),...
                        fieldnames(config),...
                        'UniformOutput',false);
    str=sprintf('%s',rows{:});

function response_is_ok=show_ok_cancel_dialog(title,body)
    response=questdlg(body,...
                        title,...
                        'Cancel','OK','OK');

    response_is_ok=strcmp(response,'OK');

function [config,error_msg]=update_config(config)
    % get keys
    path_keys=get_path_keys();
    path_key_fields=fieldnames(path_keys);
    n_paths=numel(path_key_fields);

    start_path='';
    for k=1:n_paths
        key=path_key_fields{k};
        desc=path_keys.(key);

        title=sprintf('Set ''%s''',key);
        body=sprintf(['In the next dialog, select '...
                                'the directory ' desc]);

        if isfield(config,key)
            current_config_path=config.(key);
            if isdir(current_config_path)
                % if it already has a value, use that as starting value
                start_path=current_config_path;
                body=sprintf(['%s\n\nThe current setting '...
                                'for ''%s'' is:\n%s'],...
                                body,key,start_path);
            end
        end

        response_is_ok=show_ok_cancel_dialog(title,body);

        if response_is_ok

            result_dir=uigetdir(start_path,desc);

            if ischar(result_dir)
                start_path=result_dir;
                config.(key)=result_dir;
                continue;
            end
        end

        % no valid response
        error_msg=sprintf(['Selecting directory for key ''%s'' '...
                                'was cancelled'],key);
        return;
    end

    error_msg='';

function path_keys=get_path_keys()
    path_keys=struct();
    path_keys.tutorial_data_path=[...
                        'where CoSMoMVPA tutorial and example'...
                        'data is stored. This data is required '...
                        'for running the exercises and '...
                        'demonstrations that come with '...
                        'CoSMoMVPA (in the examples/ '...
                        'directory)\n'...
                        'If you have not '...
                        'downloaded this data, select ''Cancel'' '...
                        'and see the '...
                        'website for download instructions:\n\n'...
                        '    http://cosmomvpa.org/download.html'...
                        ''];
    path_keys.output_data_path=[...
                        'where output from running the exercises and '...
                        'demonstrations (such as fMRI volumes with '...
                        'searchlight result maps) is be stored. This '...
                        'can be '...
                        'any directory of your liking that you have '...
                        'write permissions to'];

function [path_fn,error_msg]=get_output_filename()
    fn='.cosmomvpa.cfg';
    error_msg='';

    if isunix()
        home_directory=getenv('HOME');
        candidate_paths={home_directory};
    else
        candidate_paths={};
    end

    if ~isempty(which('userpath','builtin'))
        % this does not work on GNU Octave
        candidate_paths{end+1}=userpath();
    end

    % not preferred, but no other options available
    % (typically, GNU Octave on Windows)
    candidate_paths{end+1}=get_cosmo_mvpa_dir();

    sep=pathsep();
    candidates_strjoined=cosmo_strjoin(candidate_paths,sep);
    candidates_paths=cosmo_strsplit(candidates_strjoined,sep);

    n_candidates=numel(candidates_paths);
    for k=1:n_candidates
        p=candidates_paths{k};
        if isempty(p)
            continue;
        end

        path_fn=fullfile(p,fn);
        if can_write(path_fn);
            return
        end
    end

    path_fn=[];
    error_msg=sprintf(['No suitable directory found where %s '...
                'could be written. The following directories were '...
                'tried but failed: ''%s''\n\n'],...
                fn,cosmo_strjoin(candidates_paths,''', '''));

function tf=can_write(fn)
    file_existed=exist(fn,'file');

    % try to open in append mode
    fid=fopen(fn,'a');
    file_closer=onCleanup(@()fclose(fid));
    tf=fid>0;

    clear file_closer;

    if ~file_existed
        delete(fn);
    end

function [config,fn]=get_cosmo_config()
    warning_state=cosmo_warning();
    cleaner=onCleanup(@()cosmo_warning(warning_state));
    cosmo_warning('off');

    [config,fn]=cosmo_config();

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
clear, clc;

% This is an example for running the function tree_LeastR
%
%  Problem:
%
%  min  1/2 || A x - y||^2 + z * sum_j w_j ||x_{G_j}||
%
%  G_j's are nodes with tree structure
%
%  The tree structured group information is contained in
%  opts.ind, which is a 3 x nodes matrix, where nodes denotes the number of
%  nodes of the tree.
%
%  opts.ind(1,:) contains the starting index
%  opts.ind(2,:) contains the ending index
%  opts.ind(3,:) contains the corresponding weight (w_j)
%
%  Note:
%  1) If each element of x is a leaf node of the tree and the weight for
%  this leaf node are the same, we provide an alternative "efficient" input
%  for this kind of node, by creating a "super node" with
%  opts.ind(1,1)=-1; opts.ind(2,1)=-1; and opts.ind(3,1)=the common weight.
%
%  2) If the features are well ordered in that, the features of the left
%  tree is always less than those of the right tree, opts.ind(1,:) and
%  opts.ind(2,:) contain the "real" starting and ending indices. That is to
%  say, x( opts.ind(1,j):opts.ind(2,j) ) denotes x_{G_j}. In this case,
%  the entries in opts.ind(1:2,:) are within 1 and n.
%
%
%  If the features are not well ordered, please use the input opts.G for
%  specifying the index so that
%   x( opts.G ( opts.ind(1,j):opts.ind(2,j) ) ) denotes x_{G_j}.
%  In this case, the entries of opts.G are within 1 and n, and the entries of
%  opts.ind(1:2,:) are within 1 and length(opts.G).
%
%% Related papers
%
% [1] Jun Liu and Jieping Ye, Moreau-Yosida Regularization for
%     Grouped Tree Structure Learning, NIPS 2010
%
%%

cd ..
cd ..

root=cd;
addpath(genpath([root '/SLEP']));
                     % add the functions in the folder SLEP to the path

% change to the original folder
cd Examples/tree;

m=50;  n=100;       % The data matrix is of size m x n

%randNum=10;           % a random number

% ---------------------- generate random data ----------------------
%randn('state',(randNum-1)*3+1);
A=randn(m,n);        % the data matrix

%randn('state',(randNum-1)*3+2);
xOrin=full(sprandn(n, 1,1));
xOrin(1:50,:)=0;

%randn('state',(randNum-1)*3+3);
noise=randn(m,1);
y=A*xOrin +...
    noise*0.01;      % the response

%% In this example, the tree is set as:
%
% root, 1:100, with weight 0
% its children nodes, 1:50, and 51:100
%
% For 1:50, its children are 1:20, 21:40, and 41:50
%
% For 51:100, its children are 51:70, and 71:100
%
% These nodes in addition have each individual features (they contain) as
% children nodes.
%
%%

%% One efficient way
% We make use of the fact that the indices of the left nodes of the tree
% are smaller than the right nodes.

%----------------------- Set optional items -----------------------
opts=[];

% Starting point
opts.init=2;        % starting from a zero point

% Termination
opts.tFlag=5;       % run .maxIter iterations
opts.maxIter=100;   % maximum number of iterations

% regularization
opts.rFlag=1;       % use ratio

% Normalization
opts.nFlag=0;       % without normalization

% Group Property
opts.ind=[[-1, -1, 1]',... % leave nodes (each node contains one feature)
    [1, 20, sqrt(20)]', [21, 40, sqrt(20)]',... % the layer above the leaf
    [41, 50, sqrt(10)]', [51, 70, sqrt(20)]', [71,100, sqrt(30)]',...
    [1, 50, sqrt(50)]', [51, 100, sqrt(50)]']; % the higher layer

%----------------------- Run the code tree_LeastR -----------------------
z=0.1;
tic;
[x1, funVal1, ValueL1]= tree_LeastR(A, y, z, opts);
toc;

%% An alternative way
% We make use of the fact that the indices of the left nodes of the tree
% are smaller than the right nodes.
%%

%----------------------- Set optional items -----------------------
opts=[];

% Starting point
opts.init=2;        % starting from a zero point

% Termination
opts.tFlag=5;       % run .maxIter iterations
opts.maxIter=100;   % maximum number of iterations

% regularization
opts.rFlag=1;       % use ratio

% Normalization
opts.nFlag=0;       % without normalization

% Group Property
opts.ind=[[-1, -1, 1]',... % leave nodes (each node contains one feature)
    [1, 20, sqrt(20)]', [21, 40, sqrt(20)]',... % the layer above the leaf
    [41, 50, sqrt(10)]', [51, 70, sqrt(20)]', [71,100, sqrt(30)]',...
    [101, 150, sqrt(50)]', [151, 200, sqrt(50)]']; % the higher layer
opts.G=[1:100, 1:100];

%----------------------- Run the code tree_LeastR -----------------------
z=0.1;
tic;
[x2, funVal2, ValueL2]= tree_LeastR(A, y, z, opts);
toc;

author: zzstefan

repo: BrainNetClass