text:
function CFWCR_VOT()

	% *************************************************************
	% VOT: Always call exit command at the end to terminate Matlab!
	% *************************************************************
	cleanup = onCleanup(@() exit() );

	% *************************************************************
	% VOT: Set random seed to a different value every time.
	% *************************************************************
	RandStream.setGlobalStream(RandStream('mt19937ar', 'Seed', sum(clock)));

	% *************************************************************
	% VOT: init the resources
	% *************************************************************
	[wrapper_path, name, ext] = fileparts(mfilename('fullpath'));
	addpath(wrapper_path);
	cd_ind = strfind(wrapper_path, filesep());
	repo_path = wrapper_path(1:cd_ind(end)-1);
	addpath(repo_path);
	setup_paths();
	vl_setupnn();
	% **********************************
	% VOT: Get initialization data
	% **********************************
	[handle, image, region] = vot('polygon');
	% Initialize the tracker
	disp(image);
	bb_scale = 1;

	% If the provided region is a polygon ...
	if numel(region) > 4
	    % Init with an axis aligned bounding box with correct area and center
	    % coordinate
	    cx = mean(region(1:2:end));
	    cy = mean(region(2:2:end));
	    x1 = min(region(1:2:end));
	    x2 = max(region(1:2:end));
	    y1 = min(region(2:2:end));
	    y2 = max(region(2:2:end));
	    A1 = norm(region(1:2) - region(3:4)) * norm(region(3:4) - region(5:6));
	    A2 = (x2 - x1) * (y2 - y1);
	    s = sqrt(A1/A2);
	    w = s * (x2 - x1) + 1;
	    h = s * (y2 - y1) + 1;
	else
	    cx = region(1) + (region(3) - 1)/2;
	    cy = region(2) + (region(4) - 1)/2;
	    w = region(3);
	    h = region(4);
	end
	init_c = [cx cy];
	init_sz = bb_scale * [w h];

	im_size = size(imread(image));
	im_size = im_size([2 1]);

	init_pos = min(max(round(init_c - (init_sz - 1)/2), [1 1]), im_size);
	init_sz = min(max(round(init_sz), [1 1]), im_size - init_pos + 1);
	region = [init_pos, init_sz];
	frame = 1

	% **********************************
	% VOT: ECO init
	% **********************************
	params = init_param(region);

	max_train_samples = params.nSamples;
	features = params.t_features;

	% Set some default parameters
	params = init_default_params(params);
	if isfield(params, 't_global')
	    global_fparams = params.t_global;
	else
	    global_fparams = [];
	end

	% Init sequence data
	pos = params.init_pos(:)';
	target_sz = params.init_sz(:)';

	init_target_sz = target_sz;

	% Check if color image
	im = imread(image);
	if size(im,3) == 3
	    if all(all(im(:,:,1) == im(:,:,2)))
	        is_color_image = false;
	    else
	        is_color_image = true;
	    end
	else
	    is_color_image = false;
	end

	if size(im,3) > 1 && is_color_image == false
	    im = im(:,:,1);
	end

    params.use_mexResize = false;
    global_fparams.use_mexResize = false;
	% Calculate search area and initial scale factor
	search_area = prod(init_target_sz * params.search_area_scale);
	if search_area > params.max_image_sample_size
	    currentScaleFactor = sqrt(search_area / params.max_image_sample_size);
	elseif search_area < params.min_image_sample_size
	    currentScaleFactor = sqrt(search_area / params.min_image_sample_size);
	else
	    currentScaleFactor = 1.0;
	end

	% target size at the initial scale
	base_target_sz = target_sz / currentScaleFactor;
	% window size, taking padding into account
	switch params.search_area_shape
	    case 'proportional'
	        img_sample_sz = floor( base_target_sz * params.search_area_scale);     % proportional area, same aspect ratio as the target
	    case 'square'
	        img_sample_sz = repmat(sqrt(prod(base_target_sz * params.search_area_scale)), 1, 2); % square area, ignores the target aspect ratio
	    case 'fix_padding'
	        img_sample_sz = base_target_sz + sqrt(prod(base_target_sz * params.search_area_scale) + (base_target_sz(1) - base_target_sz(2))/4) - sum(base_target_sz)/2; % const padding
	    case 'custom'
	        img_sample_sz = [base_target_sz(1)*2 base_target_sz(2)*2]; % for testing
	end

	[features, global_fparams, feature_info] = init_features(features, global_fparams, is_color_image, img_sample_sz, 'odd_cells');
	% Set feature info
	img_support_sz = feature_info.img_support_sz;
	feature_sz = feature_info.data_sz;
	feature_dim = feature_info.dim;
	num_feature_blocks = length(feature_dim);
	feature_reg = permute(num2cell(feature_info.penalty), [2 3 1]);

	% Get feature specific parameters
	feature_params = init_feature_params(features, feature_info);
	feature_extract_info = get_feature_extract_info(features);
	if params.use_projection_matrix
	    compressed_dim = feature_params.compressed_dim;
	else
	    compressed_dim = feature_dim;
	end
	compressed_dim_cell = permute(num2cell(compressed_dim), [2 3 1]);
	% Size of the extracted feature maps
	feature_sz_cell = permute(mat2cell(feature_sz, ones(1,num_feature_blocks), 2), [2 3 1]);

	% Number of Fourier coefficients to save for each filter layer. This will
	% be an odd number.
	filter_sz = feature_sz + mod(feature_sz+1, 2);
	filter_sz_cell = permute(mat2cell(filter_sz, ones(1,num_feature_blocks), 2), [2 3 1]);

	% The size of the label function DFT. Equal to the maximum filter size.
	output_sz = max(filter_sz, [], 1);

	% How much each feature block has to be padded to the obtain output_sz
	pad_sz = cellfun(@(filter_sz) (output_sz - filter_sz) / 2, filter_sz_cell, 'uniformoutput', false);

	% Compute the Fourier series indices and their transposes
	ky = cellfun(@(sz) (-ceil((sz(1) - 1)/2) : floor((sz(1) - 1)/2))', filter_sz_cell, 'uniformoutput', false);
	kx = cellfun(@(sz) -ceil((sz(2) - 1)/2) : 0, filter_sz_cell, 'uniformoutput', false);

	% construct the Gaussian label function using Poisson formula
	sig_y = sqrt(prod(floor(base_target_sz))) * params.output_sigma_factor * (output_sz ./ img_support_sz);
	yf_y = cellfun(@(ky) single(sqrt(2*pi) * sig_y(1) / output_sz(1) * exp(-2 * (pi * sig_y(1) * ky / output_sz(1)).^2)), ky, 'uniformoutput', false);
	yf_x = cellfun(@(kx) single(sqrt(2*pi) * sig_y(2) / output_sz(2) * exp(-2 * (pi * sig_y(2) * kx / output_sz(2)).^2)), kx, 'uniformoutput', false);
	yf = cellfun(@(yf_y, yf_x) yf_y * yf_x, yf_y, yf_x, 'uniformoutput', false);

	% construct cosine window
	cos_window = cellfun(@(sz) single(hann(sz(1)+2)*hann(sz(2)+2)'), feature_sz_cell, 'uniformoutput', false);
	cos_window = cellfun(@(cos_window) cos_window(2:end-1,2:end-1), cos_window, 'uniformoutput', false);

	% Compute Fourier series of interpolation function
	[interp1_fs, interp2_fs] = cellfun(@(sz) get_interp_fourier(sz, params), filter_sz_cell, 'uniformoutput', false);
	% Get the reg_window_edge parameter
	reg_window_edge = {};
	for k = 1:length(features)
	    if isfield(features{k}.fparams, 'reg_window_edge')
	        reg_window_edge = cat(3, reg_window_edge, permute(num2cell(features{k}.fparams.reg_window_edge(:)), [2 3 1]));
	    else
	        reg_window_edge = cat(3, reg_window_edge, cell(1, 1, length(features{k}.fparams.nDim)));
	    end
	end

	% Construct spatial regularization filter
	reg_filter = cellfun(@(reg_window_edge) get_reg_filter(img_support_sz, base_target_sz, params, reg_window_edge), reg_window_edge, 'uniformoutput', false);

	% Compute the energy of the filter (used for preconditioner)
	reg_energy = cellfun(@(reg_filter) real(reg_filter(:)' * reg_filter(:)), reg_filter, 'uniformoutput', false);

	if params.use_scale_filter
	    [nScales, scale_step, scaleFactors, scale_filter, params] = init_scale_filter(params);
	else
	    % Use the translation filter to estimate the scale.
	    nScales = params.number_of_scales;
	    scale_step = params.scale_step;
	    scale_exp = (-floor((nScales-1)/2):ceil((nScales-1)/2));
	    scaleFactors = scale_step .^ scale_exp;
	end

	if nScales > 0
	    %force reasonable scale changes
	    min_scale_factor = scale_step ^ ceil(log(max(5 ./ img_support_sz)) / log(scale_step));
	    max_scale_factor = scale_step ^ floor(log(min([size(im,1) size(im,2)] ./ base_target_sz)) / log(scale_step));
	end

	% Set conjugate gradient uptions
	init_CG_opts.CG_use_FR = true;
	init_CG_opts.tol = 1e-6;
	init_CG_opts.CG_standard_alpha = true;
	init_CG_opts.debug = 0;
	CG_opts.CG_use_FR = params.CG_use_FR;
	CG_opts.tol = 1e-6;
	CG_opts.CG_standard_alpha = params.CG_standard_alpha;
	CG_opts.debug = 0;

	time = 0;

	% Initialize and allocate
	prior_weights = zeros(max_train_samples,1, 'single');
	sample_weights = prior_weights;
	samplesf = cell(1, 1, num_feature_blocks);
	for k = 1:num_feature_blocks
	    samplesf{k} = complex(zeros(max_train_samples,compressed_dim(k),filter_sz(k,1),(filter_sz(k,2)+1)/2,'single'));
	end
	score_matrix = inf(max_train_samples, 'single');

	latest_ind = [];
	frames_since_last_train = inf;
	num_training_samples = 0;
	minimum_sample_weight = params.learning_rate*(1-params.learning_rate)^(2*max_train_samples);

	res_norms = [];
	residuals_pcg = [];
	if frame == 1
        % Extract image region for training sample
        sample_pos = round(pos);
        sample_scale = currentScaleFactor;
        xl = extract_features(im, sample_pos, currentScaleFactor, features, global_fparams, feature_extract_info);

        % Do windowing of features
        xlw = cellfun(@(feat_map, cos_window) bsxfun(@times, feat_map, cos_window), xl, cos_window, 'uniformoutput', false);

        % Compute the fourier series
        xlf = cellfun(@cfft2, xlw, 'uniformoutput', false);

        % Interpolate features to the continuous domain
        xlf = interpolate_dft(xlf, interp1_fs, interp2_fs);

        % New sample to be added
        xlf = compact_fourier_coeff(xlf);

        % Initialize projection matrix
        xl1 = cellfun(@(x) reshape(x, [], size(x,3)), xl, 'uniformoutput', false);
        xl1 = cellfun(@(x) bsxfun(@minus, x, mean(x, 1)), xl1, 'uniformoutput', false);

        if strcmpi(params.proj_init_method, 'pca')
            [projection_matrix, ~, ~] = cellfun(@(x) svd(x' * x), xl1, 'uniformoutput', false);
            projection_matrix = cellfun(@(P, dim) single(P(:,1:dim)), projection_matrix, compressed_dim_cell, 'uniformoutput', false);
        elseif strcmpi(params.proj_init_method, 'rand_uni')
            projection_matrix = cellfun(@(x, dim) single(randn(size(x,2), dim)), xl1, compressed_dim_cell, 'uniformoutput', false);
            projection_matrix = cellfun(@(P) bsxfun(@rdivide, P, sqrt(sum(P.^2,1))), projection_matrix, 'uniformoutput', false);
        elseif strcmpi(params.proj_init_method, 'none')
            projection_matrix = [];
        else
            error('Unknown initialization method for the projection matrix: %s', params.proj_init_method);
        end
        clear xl1 xlw

        % Shift sample
        shift_samp = 2*pi * (pos - sample_pos) ./ (sample_scale * img_support_sz);
        xlf = shift_sample(xlf, shift_samp, kx, ky);

        % Project sample
        xlf_proj = project_sample(xlf, projection_matrix);
    elseif params.learning_rate > 0
        if ~params.use_detection_sample
            % Extract image region for training sample
            sample_pos = round(pos);
            sample_scale = currentScaleFactor;
            xl = extract_features(im, sample_pos, currentScaleFactor, features, global_fparams, feature_extract_info);

            % Project sample
            xl_proj = project_sample(xl, projection_matrix);

            % Do windowing of features
            xl_proj = cellfun(@(feat_map, cos_window) bsxfun(@times, feat_map, cos_window), xl_proj, cos_window, 'uniformoutput', false);

            % Compute the fourier series
            xlf1_proj = cellfun(@cfft2, xl_proj, 'uniformoutput', false);

            % Interpolate features to the continuous domain
            xlf1_proj = interpolate_dft(xlf1_proj, interp1_fs, interp2_fs);

            % New sample to be added
            xlf_proj = compact_fourier_coeff(xlf1_proj);
        else
            % Use the sample that was used for detection
            sample_scale = sample_scale(scale_ind);
            xlf_proj = cellfun(@(xf) xf(:,1:(size(xf,2)+1)/2,:,scale_ind), xtf_proj, 'uniformoutput', false);
        end

        % Shift the sample so that the target is centered
        shift_samp = 2*pi * (pos - sample_pos) ./ (sample_scale * img_support_sz);
        xlf_proj = shift_sample(xlf_proj, shift_samp, kx, ky);
    end
    xlf_proj_perm = cellfun(@(xf) permute(xf, [4 3 1 2]), xlf_proj, 'uniformoutput', false);

	if params.use_sample_merge
        % Find the distances with existing samples
        dist_vector = find_cluster_distances(samplesf, xlf_proj_perm, num_feature_blocks, num_training_samples, max_train_samples, params);

        [merged_sample, new_cluster, merged_cluster_id, new_cluster_id, score_matrix, prior_weights,num_training_samples] = ...
            merge_clusters(samplesf, xlf_proj_perm, dist_vector, score_matrix, prior_weights,...
                           num_training_samples,num_feature_blocks,max_train_samples,minimum_sample_weight,params);
    else
        % Do the traditional adding of a training sample and weight update
        % of C-COT
        [prior_weights, replace_ind] = update_prior_weights(prior_weights, sample_weights, latest_ind, frame, params);
        latest_ind = replace_ind;

        merged_cluster_id = 0;
        new_cluster = xlf_proj_perm;
        new_cluster_id = replace_ind;
    end

    if frame > 1 && params.learning_rate > 0 || frame == 1 && ~params.update_projection_matrix
        % Insert the new training sample
        for k = 1:num_feature_blocks
            if merged_cluster_id > 0
                samplesf{k}(merged_cluster_id,:,:,:) = merged_sample{k};
            end

            if new_cluster_id > 0
                samplesf{k}(new_cluster_id,:,:,:) = new_cluster{k};
            end
        end
    end

    sample_weights = prior_weights;

    train_tracker = (frame < params.skip_after_frame) || (frames_since_last_train >= params.train_gap);
    if train_tracker
        % Used for preconditioning
        new_sample_energy = cellfun(@(xlf) abs(xlf .* conj(xlf)), xlf_proj, 'uniformoutput', false);

        if frame == 1
            if params.update_projection_matrix
                hf = cell(2,1,num_feature_blocks);
                lf_ind = cellfun(@(sz) sz(1) * (sz(2)-1)/2 + 1, filter_sz_cell, 'uniformoutput', false);
                proj_energy = cellfun(@(P, yf) 2*sum(abs(yf(:)).^2) / sum(feature_dim) * ones(size(P), 'single'), projection_matrix, yf, 'uniformoutput', false);
            else
                hf = cell(1,1,num_feature_blocks);
            end
            % Initialize the filter
            for k = 1:num_feature_blocks
                hf{1,1,k} = complex(zeros([filter_sz(k,1) (filter_sz(k,2)+1)/2 compressed_dim(k)], 'single'));
            end

            % Initialize Conjugate Gradient parameters
            CG_opts.maxit = params.init_CG_iter; % Number of initial iterations if projection matrix is not updated
            init_CG_opts.maxit = ceil(params.init_CG_iter / params.init_GN_iter);
            sample_energy = new_sample_energy;
            rhs_samplef = cell(size(hf));
            diag_M = cell(size(hf));
            p = []; rho = []; r_old = [];
        else
            CG_opts.maxit = params.CG_iter;

            if params.CG_forgetting_rate == inf || params.learning_rate >= 1
                % CG will be reset
                p = []; rho = []; r_old = [];
            else
                rho = rho / (1-params.learning_rate)^params.CG_forgetting_rate;
            end
            % Update the approximate average sample energy using the learning
            % rate. This is only used to construct the preconditioner.
            sample_energy = cellfun(@(se, nse) (1 - params.learning_rate) * se + params.learning_rate * nse, sample_energy, new_sample_energy, 'uniformoutput', false);
        end
        % Do training
        if frame == 1 && params.update_projection_matrix
            % Initial Gauss-Newton optimization of the filter and
            % projection matrix.

            % Construct stuff for the proj matrix part
            init_samplef = cellfun(@(x) permute(x, [4 3 1 2]), xlf, 'uniformoutput', false);
            init_samplef_H = cellfun(@(X) conj(reshape(X, size(X,2), [])), init_samplef, 'uniformoutput', false);

            % Construct preconditioner
            diag_M(1,1,:) = cellfun(@(m, reg_energy) (1-params.precond_reg_param) * bsxfun(@plus, params.precond_data_param * m, (1-params.precond_data_param) * mean(m,3)) + params.precond_reg_param*reg_energy, sample_energy, reg_energy, 'uniformoutput',false);
            diag_M(2,1,:) = cellfun(@(m) params.precond_proj_param * (m + params.projection_reg), proj_energy, 'uniformoutput',false);

            projection_matrix_init = projection_matrix;
            for iter = 1:params.init_GN_iter
                % Project sample with new matrix
                init_samplef_proj = cellfun(@(x,P) mtimesx(x, P, 'speed'), init_samplef, projection_matrix, 'uniformoutput', false);
                init_hf = cellfun(@(x) permute(x, [3 4 1 2]), hf(1,1,:), 'uniformoutput', false);

                % Construct the right hand side vector for the filter part
                rhs_samplef(1,1,:) = cellfun(@(xf, yf) bsxfun(@times, conj(permute(xf, [3 4 2 1])), yf), init_samplef_proj, yf, 'uniformoutput', false);

                % Construct the right hand side vector for the projection matrix part
                fyf = cellfun(@(f, yf) reshape(bsxfun(@times, conj(f), yf), [], size(f,3)), hf(1,1,:), yf, 'uniformoutput', false);
                rhs_samplef(2,1,:) = cellfun(@(P, XH, fyf, fi) (2*real(XH * fyf - XH(:,fi:end) * fyf(fi:end,:)) - params.projection_reg * P), ...
                    projection_matrix, init_samplef_H, fyf, lf_ind, 'uniformoutput', false);

                % Initialize the projection matrix increment to zero
                hf(2,1,:) = cellfun(@(P) zeros(size(P), 'single'), projection_matrix, 'uniformoutput', false);

                % do conjugate gradient
                [hf, ~, ~, ~, res_norms_temp] = pcg_ccot(...
                    @(x) lhs_operation_joint(x, init_samplef_proj, reg_filter, feature_reg, init_samplef, init_samplef_H, init_hf, params.projection_reg),...
                    rhs_samplef, init_CG_opts, ...
                    @(x) diag_precond(x, diag_M), ...
                    [], hf);

                % Make the filter symmetric (avoid roundoff errors)
                hf(1,1,:) = symmetrize_filter(hf(1,1,:));

                % Add to the projection matrix
                projection_matrix = cellfun(@plus, projection_matrix, hf(2,1,:), 'uniformoutput', false);

                res_norms = [res_norms; res_norms_temp];
            end
            % Extract filter
            hf = hf(1,1,:);
            % Re-project and insert training sample
            xlf_proj = project_sample(xlf, projection_matrix);
            for k = 1:num_feature_blocks
                samplesf{k}(1,:,:,:) = permute(xlf_proj{k}, [4 3 1 2]);
            end
        else
            % Construct the right hand side vector
            rhs_samplef = cellfun(@(xf) permute(mtimesx(sample_weights, 'T', xf, 'speed'), [3 4 2 1]), samplesf, 'uniformoutput', false);
            rhs_samplef = cellfun(@(xf, yf) bsxfun(@times, conj(xf), yf), rhs_samplef, yf, 'uniformoutput', false);

            % Construct preconditioner
            diag_M = cellfun(@(m, reg_energy) (1-params.precond_reg_param) * bsxfun(@plus, params.precond_data_param * m, (1-params.precond_data_param) * mean(m,3)) + params.precond_reg_param*reg_energy, sample_energy, reg_energy, 'uniformoutput',false);

            % do conjugate gradient
            [hf, ~, ~, ~, res_norms, p, rho, r_old] = pcg_ccot(...
                @(x) lhs_operation(x, samplesf, reg_filter, sample_weights, feature_reg),...
                rhs_samplef, CG_opts, ...
                @(x) diag_precond(x, diag_M), ...
                [], hf, p, rho, r_old);
        end

        % Reconstruct the full Fourier series
        hf_full = full_fourier_coeff(hf);
        frames_since_last_train = 0;
    else
        frames_since_last_train = frames_since_last_train+1;
    end

    % Update the scale filter
    if nScales > 0 && params.use_scale_filter
        scale_filter = scale_filter_update(im, pos, base_target_sz, currentScaleFactor, scale_filter, params);
    end
    % Update the target size (only used for computing output box)
    target_sz = base_target_sz * currentScaleFactor;
	while true
	    % **********************************
	    % VOT: Get next frame
	    % **********************************
	    [handle, image] = handle.frame(handle);
	    if isempty(image)
	        break;
	    end;
	    disp(image);
	    frame = frame + 1;

		% **********************************
		% VOT: ECO update
		% **********************************
		im = imread(image);
    	if size(im,3) > 1 && is_color_image == false
        	im = im(:,:,1);
        end
		if frame > 1
	        old_pos = inf(size(pos));
	        iter = 1;
	        %translation search
	        while iter <= params.refinement_iterations && any(old_pos ~= pos)
	            % Extract features at multiple resolutions
	            sample_pos = round(pos);
	            det_sample_pos = sample_pos;
	            sample_scale = currentScaleFactor*scaleFactors;
	            xt = extract_features(im, sample_pos, sample_scale, features, global_fparams, feature_extract_info);
	            % Project sample
	            xt_proj = project_sample(xt, projection_matrix);

	            % Do windowing of features
	            xt_proj = cellfun(@(feat_map, cos_window) bsxfun(@times, feat_map, cos_window), xt_proj, cos_window, 'uniformoutput', false);

	            % Compute the fourier series
	            xtf_proj = cellfun(@cfft2, xt_proj, 'uniformoutput', false);

	            % Interpolate features to the continuous domain
	            xtf_proj = interpolate_dft(xtf_proj, interp1_fs, interp2_fs);
	            % Compute convolution for each feature block in the Fourier domain
	            scores_fs_feat = cellfun(@(hf, xf, pad_sz) padarray(sum(bsxfun(@times, hf, xf), 3), pad_sz), hf_full, xtf_proj, pad_sz, 'uniformoutput', false);
	            switch params.weights_type
						    case 'constant'
						    	scores_fs_feat{1,1,1} = param.weights(1) *scores_fs_feat{1,1,1};
						    	scores_fs_feat{1,1,2} = param.weights(2) *scores_fs_feat{1,1,1};
						    case 'sigmoid'
						    	coe = params.initial - params.factor./ (1 + exp(-double(frame)/params.divide_denominator));
						    	scores_fs_feat{1,1,1} = 1*scores_fs_feat{1,1,1};
	            		scores_fs_feat{1,1,2} = coe*scores_fs_feat{1,1,2};
							end
	            % Also sum over all feature blocks.
	            % Gives the fourier coefficients of the convolution response.
	            scores_fs = permute(sum(cell2mat(scores_fs_feat), 3), [1 2 4 3]);

	            % Optimize the continuous score function with Newton's method.
	            [trans_row, trans_col, scale_ind] = optimize_scores(scores_fs, params.newton_iterations);
	            % Compute the translation vector in pixel-coordinates and round
	            % to the closest integer pixel.
	            translation_vec = [trans_row, trans_col] .* (img_support_sz./output_sz) * currentScaleFactor * scaleFactors(scale_ind);
	            scale_change_factor = scaleFactors(scale_ind);

	            % update position
	            old_pos = pos;
	            pos = sample_pos + translation_vec;

	            if params.clamp_position
	                pos = max([1 1], min([size(im,1) size(im,2)], pos));
	            end
	            % Do scale tracking with the scale filter
	            if nScales > 0 && params.use_scale_filter
	                scale_change_factor = scale_filter_track(im, pos, base_target_sz, currentScaleFactor, scale_filter, params);
	            end
	            % Update the scale
	            currentScaleFactor = currentScaleFactor * scale_change_factor;

	            % Adjust to make sure we are not to large or to small
	            if currentScaleFactor < min_scale_factor
	                currentScaleFactor = min_scale_factor;
	            elseif currentScaleFactor > max_scale_factor
	                currentScaleFactor = max_scale_factor;
	            end

	            iter = iter + 1;
	        end
		end
	    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	    %% Model update step
	    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	    % Extract sample and init projection matrix
	    if frame == 1
	        % Extract image region for training sample
	        sample_pos = round(pos);
	        sample_scale = currentScaleFactor;
	        xl = extract_features(im, sample_pos, currentScaleFactor, features, global_fparams, feature_extract_info);

	        % Do windowing of features
	        xlw = cellfun(@(feat_map, cos_window) bsxfun(@times, feat_map, cos_window), xl, cos_window, 'uniformoutput', false);

	        % Compute the fourier series
	        xlf = cellfun(@cfft2, xlw, 'uniformoutput', false);

	        % Interpolate features to the continuous domain
	        xlf = interpolate_dft(xlf, interp1_fs, interp2_fs);

	        % New sample to be added
	        xlf = compact_fourier_coeff(xlf);

	        % Initialize projection matrix
	        xl1 = cellfun(@(x) reshape(x, [], size(x,3)), xl, 'uniformoutput', false);
	        xl1 = cellfun(@(x) bsxfun(@minus, x, mean(x, 1)), xl1, 'uniformoutput', false);

	        if strcmpi(params.proj_init_method, 'pca')
	            [projection_matrix, ~, ~] = cellfun(@(x) svd(x' * x), xl1, 'uniformoutput', false);
	            projection_matrix = cellfun(@(P, dim) single(P(:,1:dim)), projection_matrix, compressed_dim_cell, 'uniformoutput', false);
	        elseif strcmpi(params.proj_init_method, 'rand_uni')
	            projection_matrix = cellfun(@(x, dim) single(randn(size(x,2), dim)), xl1, compressed_dim_cell, 'uniformoutput', false);
	            projection_matrix = cellfun(@(P) bsxfun(@rdivide, P, sqrt(sum(P.^2,1))), projection_matrix, 'uniformoutput', false);
	        elseif strcmpi(params.proj_init_method, 'none')
	            projection_matrix = [];
	        else
	            error('Unknown initialization method for the projection matrix: %s', params.proj_init_method);
	        end
	        clear xl1 xlw

	        % Shift sample
	        shift_samp = 2*pi * (pos - sample_pos) ./ (sample_scale * img_support_sz);
	        xlf = shift_sample(xlf, shift_samp, kx, ky);

	        % Project sample
	        xlf_proj = project_sample(xlf, projection_matrix);
	    elseif params.learning_rate > 0
	        if ~params.use_detection_sample
	            % Extract image region for training sample
	            sample_pos = round(pos);
	            sample_scale = currentScaleFactor;
	            xl = extract_features(im, sample_pos, currentScaleFactor, features, global_fparams, feature_extract_info);

	            % Project sample
	            xl_proj = project_sample(xl, projection_matrix);

	            % Do windowing of features
	            xl_proj = cellfun(@(feat_map, cos_window) bsxfun(@times, feat_map, cos_window), xl_proj, cos_window, 'uniformoutput', false);

	            % Compute the fourier series
	            xlf1_proj = cellfun(@cfft2, xl_proj, 'uniformoutput', false);

	            % Interpolate features to the continuous domain
	            xlf1_proj = interpolate_dft(xlf1_proj, interp1_fs, interp2_fs);

	            % New sample to be added
	            xlf_proj = compact_fourier_coeff(xlf1_proj);
	        else
	            % Use the sample that was used for detection
	            sample_scale = sample_scale(scale_ind);
	            xlf_proj = cellfun(@(xf) xf(:,1:(size(xf,2)+1)/2,:,scale_ind), xtf_proj, 'uniformoutput', false);
	        end

	        % Shift the sample so that the target is centered
	        shift_samp = 2*pi * (pos - sample_pos) ./ (sample_scale * img_support_sz);
	        xlf_proj = shift_sample(xlf_proj, shift_samp, kx, ky);
	    end
	    xlf_proj_perm = cellfun(@(xf) permute(xf, [4 3 1 2]), xlf_proj, 'uniformoutput', false);
	    if params.use_sample_merge
	        % Find the distances with existing samples
	        dist_vector = find_cluster_distances(samplesf, xlf_proj_perm, num_feature_blocks, num_training_samples, max_train_samples, params);

	        [merged_sample, new_cluster, merged_cluster_id, new_cluster_id, score_matrix, prior_weights,num_training_samples] = ...
	            merge_clusters(samplesf, xlf_proj_perm, dist_vector, score_matrix, prior_weights,...
	                           num_training_samples,num_feature_blocks,max_train_samples,minimum_sample_weight,params);
	    else
	        % Do the traditional adding of a training sample and weight update
	        % of C-COT
	        [prior_weights, replace_ind] = update_prior_weights(prior_weights, sample_weights, latest_ind, frame, params);
	        latest_ind = replace_ind;

	        merged_cluster_id = 0;
	        new_cluster = xlf_proj_perm;
	        new_cluster_id = replace_ind;
	    end
	    if frame > 1 && params.learning_rate > 0 || frame == 1 && ~params.update_projection_matrix
	        % Insert the new training sample
	        for k = 1:num_feature_blocks
	            if merged_cluster_id > 0
	                samplesf{k}(merged_cluster_id,:,:,:) = merged_sample{k};
	            end

	            if new_cluster_id > 0
	                samplesf{k}(new_cluster_id,:,:,:) = new_cluster{k};
	            end
	        end
	    end
	    sample_weights = prior_weights;

	    train_tracker = (frame < params.skip_after_frame) || (frames_since_last_train >= params.train_gap);
	    if train_tracker
	        % Used for preconditioning
	        new_sample_energy = cellfun(@(xlf) abs(xlf .* conj(xlf)), xlf_proj, 'uniformoutput', false);
	        if frame == 1
	            if params.update_projection_matrix
	                hf = cell(2,1,num_feature_blocks);
	                lf_ind = cellfun(@(sz) sz(1) * (sz(2)-1)/2 + 1, filter_sz_cell, 'uniformoutput', false);
	                proj_energy = cellfun(@(P, yf) 2*sum(abs(yf(:)).^2) / sum(feature_dim) * ones(size(P), 'single'), projection_matrix, yf, 'uniformoutput', false);
	            else
	                hf = cell(1,1,num_feature_blocks);
	            end
	            % Initialize the filter
	            for k = 1:num_feature_blocks
	                hf{1,1,k} = complex(zeros([filter_sz(k,1) (filter_sz(k,2)+1)/2 compressed_dim(k)], 'single'));
	            end

	            % Initialize Conjugate Gradient parameters
	            CG_opts.maxit = params.init_CG_iter; % Number of initial iterations if projection matrix is not updated
	            init_CG_opts.maxit = ceil(params.init_CG_iter / params.init_GN_iter);
	            sample_energy = new_sample_energy;
	            rhs_samplef = cell(size(hf));
	            diag_M = cell(size(hf));
	            p = []; rho = []; r_old = [];
	        else
	            CG_opts.maxit = params.CG_iter;

	            if params.CG_forgetting_rate == inf || params.learning_rate >= 1
	                % CG will be reset
	                p = []; rho = []; r_old = [];
	            else
	                rho = rho / (1-params.learning_rate)^params.CG_forgetting_rate;
	            end
	            % Update the approximate average sample energy using the learning
	            % rate. This is only used to construct the preconditioner.
	            sample_energy = cellfun(@(se, nse) (1 - params.learning_rate) * se + params.learning_rate * nse, sample_energy, new_sample_energy, 'uniformoutput', false);
	        end
	        % Do training
	        if frame == 1 && params.update_projection_matrix
	            % Initial Gauss-Newton optimization of the filter and
	            % projection matrix.

	            % Construct stuff for the proj matrix part
	            init_samplef = cellfun(@(x) permute(x, [4 3 1 2]), xlf, 'uniformoutput', false);
	            init_samplef_H = cellfun(@(X) conj(reshape(X, size(X,2), [])), init_samplef, 'uniformoutput', false);

	            % Construct preconditioner
	            diag_M(1,1,:) = cellfun(@(m, reg_energy) (1-params.precond_reg_param) * bsxfun(@plus, params.precond_data_param * m, (1-params.precond_data_param) * mean(m,3)) + params.precond_reg_param*reg_energy, sample_energy, reg_energy, 'uniformoutput',false);
	            diag_M(2,1,:) = cellfun(@(m) params.precond_proj_param * (m + params.projection_reg), proj_energy, 'uniformoutput',false);

	            projection_matrix_init = projection_matrix;

	            for iter = 1:params.init_GN_iter
	                % Project sample with new matrix
	                init_samplef_proj = cellfun(@(x,P) mtimesx(x, P, 'speed'), init_samplef, projection_matrix, 'uniformoutput', false);
	                init_hf = cellfun(@(x) permute(x, [3 4 1 2]), hf(1,1,:), 'uniformoutput', false);

	                % Construct the right hand side vector for the filter part
	                rhs_samplef(1,1,:) = cellfun(@(xf, yf) bsxfun(@times, conj(permute(xf, [3 4 2 1])), yf), init_samplef_proj, yf, 'uniformoutput', false);

	                % Construct the right hand side vector for the projection matrix part
	                fyf = cellfun(@(f, yf) reshape(bsxfun(@times, conj(f), yf), [], size(f,3)), hf(1,1,:), yf, 'uniformoutput', false);
	                rhs_samplef(2,1,:) = cellfun(@(P, XH, fyf, fi) (2*real(XH * fyf - XH(:,fi:end) * fyf(fi:end,:)) - params.projection_reg * P), ...
	                    projection_matrix, init_samplef_H, fyf, lf_ind, 'uniformoutput', false);

	                % Initialize the projection matrix increment to zero
	                hf(2,1,:) = cellfun(@(P) zeros(size(P), 'single'), projection_matrix, 'uniformoutput', false);

	                % do conjugate gradient
	                [hf, ~, ~, ~, res_norms_temp] = pcg_ccot(...
	                    @(x) lhs_operation_joint(x, init_samplef_proj, reg_filter, feature_reg, init_samplef, init_samplef_H, init_hf, params.projection_reg),...
	                    rhs_samplef, init_CG_opts, ...
	                    @(x) diag_precond(x, diag_M), ...
	                    [], hf);

	                % Make the filter symmetric (avoid roundoff errors)
	                hf(1,1,:) = symmetrize_filter(hf(1,1,:));

	                % Add to the projection matrix
	                projection_matrix = cellfun(@plus, projection_matrix, hf(2,1,:), 'uniformoutput', false);

	                res_norms = [res_norms; res_norms_temp];
	            end
	            % Extract filter
	            hf = hf(1,1,:);

	            % Re-project and insert training sample
	            xlf_proj = project_sample(xlf, projection_matrix);
	            for k = 1:num_feature_blocks
	                samplesf{k}(1,:,:,:) = permute(xlf_proj{k}, [4 3 1 2]);
	            end
	        else
	            % Construct the right hand side vector
	            rhs_samplef = cellfun(@(xf) permute(mtimesx(sample_weights, 'T', xf, 'speed'), [3 4 2 1]), samplesf, 'uniformoutput', false);
	            rhs_samplef = cellfun(@(xf, yf) bsxfun(@times, conj(xf), yf), rhs_samplef, yf, 'uniformoutput', false);

	            % Construct preconditioner
	            diag_M = cellfun(@(m, reg_energy) (1-params.precond_reg_param) * bsxfun(@plus, params.precond_data_param * m, (1-params.precond_data_param) * mean(m,3)) + params.precond_reg_param*reg_energy, sample_energy, reg_energy, 'uniformoutput',false);

	            % do conjugate gradient
	            [hf, ~, ~, ~, res_norms, p, rho, r_old] = pcg_ccot(...
	                @(x) lhs_operation(x, samplesf, reg_filter, sample_weights, feature_reg),...
	                rhs_samplef, CG_opts, ...
	                @(x) diag_precond(x, diag_M), ...
	                [], hf, p, rho, r_old);
	        end

	        % Reconstruct the full Fourier series
	        hf_full = full_fourier_coeff(hf);

	        frames_since_last_train = 0;
	    else
	        frames_since_last_train = frames_since_last_train+1;
	    end
	    % Update the scale filter
	    if nScales > 0 && params.use_scale_filter
	        scale_filter = scale_filter_update(im, pos, base_target_sz, currentScaleFactor, scale_filter, params);
	    end
	    % Update the target size (only used for computing output box)
	    target_sz = base_target_sz * currentScaleFactor;

	    %save position and calculate FPS
	    region = round([pos([2,1]) - (target_sz([2,1]) - 1)/2, target_sz([2,1])]);
	    region = double(region);
	    disp(region);

	    % **********************************
	    % VOT: Report position for frame
	    % **********************************
	    handle = handle.report(handle, region);

	end;
	% **********************************
	% VOT: Output the results
	% **********************************
	handle.quit(handle);

end

function params = init_param(region)

	cnn_params.nn_name = 'imagenet-vgg-m-2048-cut.mat'; % Name of the network
	cnn_params.output_layer = [3 14];
	cnn_params.downsample_factor = [2 1];           % How much to downsample each output layer
	cnn_params.input_size_mode = 'adaptive';        % How to choose the sample size
	cnn_params.input_size_scale = 1;                % Extra scale factor of the input samples to the network (1 is no scaling)
	cnn_params.use_gpu = true;
	cnn_params.gpu_id = [3];

	% Which features to include
	params.t_features = {
	    struct('getFeature',@get_cnn_layers, 'fparams',cnn_params),...
	};

	% Global feature parameters1s
	params.t_global.normalize_power = 2;    % Lp normalization with this p
	params.t_global.normalize_size = true;  % Also normalize with respect to the spatial size of the feature
	params.t_global.normalize_dim = true;   % Also normalize with respect to the dimensionality of the feature

	% Image sample parameters
	params.search_area_shape = 'square';    % The shape of the samples
	params.search_area_scale = 4.0;         % The scaling of the target size to get the search area
	params.min_image_sample_size = 200^2;   % Minimum area of image samples
	params.max_image_sample_size = 250^2;   % Maximum area of image samples

	% Detection parameters
	params.refinement_iterations = 1;       % Number of iterations used to refine the resulting position in a frame
	params.newton_iterations = 5;           % The number of Newton iterations used for optimizing the detection rere
	params.clamp_position = false;          % Clamp the target position to be inside the image

	% Learning parameters
	params.output_sigma_factor = 1/12;			% Label function sigma
	params.learning_rate = 0.012;	 	 				% Learning rate
	params.nSamples = 100;                  % Maximum number of stored training samples
	params.sample_replace_strategy = 'lowest_prior';    % Which sample to replace when the memory is full
	params.lt_size = 0;                     % The size of the long-term memory (where all samples have equal weight)
	params.train_gap = 5;                   % The number of intermediate frames with no training (0 corresponds to training every frame)
	params.skip_after_frame = 1;            % After which frame number the sparse update scheme should start (1 is directly)
	params.use_detection_sample = true;     % Use the sample that was extracted at the detection stage also for learning

	% Factorized convolution parameters
	params.use_projection_matrix = true;    % Use projection matrix, i.e. use the factorized convolution formulation
	params.update_projection_matrix = true; % Whether the projection matrix should be optimized or not
	params.proj_init_method = 'pca';        % Method for initializing the projection matrix
	params.projection_reg = 2e-7;	 	 				% Regularization paremeter of the projection matrix

	% Generative sample space model parameters
	params.use_sample_merge = true;                 % Use the generative sample space model to merge samples
	params.sample_update_criteria = 'Merge';        % Strategy for updating the samples
	params.weight_update_criteria = 'WeightedAdd';  % Strategy for updating the distance matrix
	params.neglect_higher_frequency = false;        % Neglect hiigher frequency components in the distance comparison for speed

	% Conjugate Gradient parameters
	params.CG_iter = 5;                     % The number of Conjugate Gradient iterations in each update after the first frame
	params.init_CG_iter = 10*20;            % The total number of Conjugate Gradient iterations used in the first frame
	params.init_GN_iter = 10;               % The number of Gauss-Newton iterations used in the first frame (only if the projection matrix is updated)
	params.CG_use_FR = false;               % Use the Fletcher-Reeves (true) or Polak-Ribiere (false) formula in the Conjugate Gradient
	params.CG_standard_alpha = true;        % Use the standard formula for computing the step length in Conjugate Gradient
	params.CG_forgetting_rate = 75;	 	 	% Forgetting rate of the last conjugate direction
	params.precond_data_param = 0.7;	 	% Weight of the data term in the preconditioner
	params.precond_reg_param = 0.1;	 	% Weight of the regularization term in the preconditioner
	params.precond_proj_param = 30;	 	 	% Weight of the projection matrix part in the preconditioner

	% Regularization window parameters
	params.use_reg_window = true;           % Use spatial regularization or not
	params.reg_window_min = 1e-4;			% The minimum value of the regularization window
	params.reg_window_edge = 10e-3;         % The impact of the spatial regularization
	params.reg_window_power = 2;            % The degree of the polynomial to use (e.g. 2 is a quadratic window)
	params.reg_sparsity_threshold = 0.12;   % A relative threshold of which DFT coefficients that should be set to zero

	% Interpolation parameters
	params.interpolation_method = 'bicubic';    % The kind of interpolation kernel
	params.interpolation_bicubic_a = -0.75;     % The parameter for the bicubic interpolation kernel
	params.interpolation_centering = true;      % Center the kernel at the feature sample
	params.interpolation_windowing = false;     % Do additional windowing on the Fourier coefficients of the kernel

	% Scale parameters for the translation model
	% Only used if: params.use_scale_filter = false
	params.use_scale_filter = false
	params.number_of_scales = 10;           % Number of scales to run the detector
	params.scale_step = 1.03;               % The scale factor

	params.weights = [1, 2];               	% The weights factor
	params.weights_type    = 'sigmoid'; 		% type: constant, sigmoid
	params.divide_denominator = 100;
	params.initial = 1;
	params.factor = 0;

	% Initialize
	params.init_sz = [region(4), region(3)];
	params.init_pos = [region(2), region(1)] + (params.init_sz - 1)/2;
end

repo: CFWCR


text:
function posV = xscale(h,scale,offset)
%function xscale(h,scale,offset)
%This script recenters and then scales the x-axis of the
%axis with handle h.  The offset value is in the unscaled
%(input) coordinates.
%JOD. LM: 23 Oct 03.
%
% 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/>.

ticklabelS = get(h,'xticklabel');
tmpV = str2num(ticklabelS);
posV = (tmpV-offset)*scale;
newlabelsS = num2str(posV);
set(h,'xticklabel',newlabelsS);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ out ] = VideoToScoreVideoSkip( video, train, s, w, thresh, skip)
%  VideoToScoreVideoSkip - get component video of video/frame passed in
%--------------------------------------------------------------------------
%   Params: video - should be a height X width X 3 X frames --- array
%           train - histograms of the training images
%           s - the window size that each frame will be split up in to form
%               histograms
%           w - the width of the bins for the RGB color
%               histograms
%           thresh - the cutoff distance threshold used to measure whether
%               or not window histograms are close enough to the training
%               histograms.
%           skip - The frames to skip if vidAllCopy is a video entered
%
%   Returns: out - component video
%--------------------------------------------------------------------------

height = length(video(:,1,1,1));
width = length(video(1,:,1,1));
frames = length(video(1,1,1,:));
H = floor(height/s);
W = floor(width/s);
out = false(H,W,frames);

for i = 1:ceil(frames/skip)
    %display(strcat('frame number: ',num2str((i-1)*skip+1)));
    out(:,:,(i-1)*skip+1) = ImageToScoreArray(video(:,:,:,(i-1)*skip+1), train, s, w, thresh);
    for j = (i-1)*skip+2:min(i*skip, frames)
        out(:,:,j) = out(:,:,(i-1)*skip+1);
    end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [A B X E] = slr_iteration(D, labels, lambda, tau, eta, delta, epsilon, maxIter, maxOut)

% Dec 2014
%
% by Jian Lai, jlai1@ntu.edu.sg
% solve the equation (17) the supervised low rand decomposition in our PAMI paper
% min ||A||_* + \lambda ||B||_* + \tau ||X||_F^2 + \eta ||E||_1 s.t. D = A + BX + E.

% input
% D: training data, each column is a atom, with the size m*n
% labels: the training labels of D, with the size n*1
% lambda, tau and eta: the balancing parameters in equation (17)
% delta: the balancing parameter of T in equation (20) and (23)
% epsilon: the tolerance of the convergent checking
% maxIter: max iteration within each subproblem
% maxOut: max iteration of outer loop (defaul value is 4)

% output
% A: the class-specific dictionary
% B: the nonclass-specific dictionary
% X: the coefficient matrix of B
% E: the error matrix

% input checking
if nargin < 6
    error('Too few arguments') ;
end

if nargin < 7
    epsilon = 1e-4;
end

if nargin < 8
    maxIter = 1000;
end

if nargin < 9
    maxIter = 4;
end

% initializate A_i as (19) in our PAMI paper
unilabel = unique(labels);
numlabel = size(unilabel,2);
for classindex = 1:numlabel
    classid = unilabel(classindex);
    Di = D(:,classid==labels);
    [U S V] = svd(Di, 'econ');
	Ai{classindex} = U(:, 1) * diag(S(1,1)) * V(:, 1)';
end

% initialize of A
A = [];
for classindex = 1:numlabel
    A = [A Ai{classindex}];
end

% initialize X to an identity matrix
X = eye(size(A,2));

% initialize B to a zero matrix
B = zeros(size(D));

for outiter = 1:maxOut

    fprintf('The %d iteration of SLR\n', outiter);

    % solve the subproblem (20) in our PAMI paper
    [B, T] = slr_BT_ialm(D, A, X, lambda, delta, epsilon, maxIter);

    % solve the subproblem (23) in our PAMI paper
    [X, T] = slr_XT_ialm(D, A, B, tau, delta, epsilon, maxIter);

    % solve the subproblem (25) in our PAMI paper
    [A, E] = slr_AE_ialm(D, B, X, eta, epsilon, maxIter);

end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function [X, info] = IRmrnsd(A, b, varargin)
%IRmrnsd Modified Residual Norm Steepest Descent method
%
% options  = IRmrnsd('defaults')
% [X,info] = IRmrnsd(A,b)
% [X,info] = IRmrnsd(A,b,K)
% [X,info] = IRmrnsd(A,b,options)
% [X,info] = IRmrnsd(A,b,K,options)
%
% This function implements the Modified Residual Norm Steepest Descent
% method for computing a nonnegatiely constrained least squares solution.
%
% With 'defaults' as input returns the default options.  Otherwise outputs
% the iterates specified in K, using max(K) as MaxIter, and using all other
% default options.  With options as input: uses the user-specified options
% and all the other default options.
%
% Inputs:
%  A : either (a) a full or sparse matrix
%             (b) a matrix object that performs the matrix*vector operation
%             (c) user-defined function handle
%  b : right-hand side vector
%  K : (optional) integer vector that specifies which iterates are returned
%      in X; the maximum number of iterations is assumed to be max(K)
%      [ positive integer | vector of positive components ]
%  options : structure with the following fields (optional)
%      x0         - Initial guess for the iterations
%                   [ array | {'none'} ]
%      MaxIter    - maximum allowed number of cgls iterations
%                   [ {100 } | positive integer ]
%                   NOTE: K overrules MaxIter if both are assigned
%      x_true     - true solution; allows us to returns error norms with
%                   respect to x_true at each iteration
%                   [ array | {'none'} ]
%      NoiseLevel - norm of noise in rhs divided by norm of rhs
%                   [ {'none'} | nonnegative scalar]
%      eta        - safety factor for the discrepancy principle
%                   [ {1.01} | scalar greater than (and close to) 1 ]
%      NE_Rtol    - relative tolerance on the normal equation residual norm
%                   [ {1e-12} | positive integer ]
%      NoStop     - specifies whether the iterations should proceed
%                   after a stopping criterion has been satisfied
%                   [ 'on' | {'off'} ]
%      IterBar    - shows the progress of the iterations
%                   [ {'on'} | 'off' ]
% Note: the options structure can be created using the function IRset.
%
% Outputs:
%   X : computed solutions, stored column-wise (at the iterations listed in K)
%   info: structure with the following fields:
%      its      - number of the last computed iteration
%      saved_iterations - iteration numbers of iterates stored in X
%      StopFlag - string that describes the inner stopping condition:
%                   * Residual tolerance satisfied (discrepancy principle)
%                   * Normal equations residual tolerance satisfied
%                   * Reached maximum number of iterations
%      Rnrm     - relative residual norms at each iteration
%      NE_Rnrm  - normal eqs relative residual norms at each iteration
%      Xnrm     - solution norms at each iteration
%      Enrm     - relative error norms (requires x_true) at each iteration
%      StopReg  - struct containing information about the solution that
%                 satisfies the stopping criterion.  Fields:
%                   It   : iteration where the stopping criterion is satisfied
%                   X    : solution satisfying the stopping criterion
%                   Enrm : the corresponding relative error (requires x_true)
%      BestReg  - struct containing information about the solution that
%                 minimizes Enrm (requires x_true), with the fields:
%                   It   : iteration where the minimum is attained
%                   X    : best solution
%                   Enrm : best relative error
%
% See also: IRconstr_ls, IRfista, IRnnfcgls, IRget, IRset

% Silvia Gazzola, University of Bath
% Per Christian Hansen, Technical University of Denmark
% James G. Nagy, Emory University
% April, 2018.

% This file is part of the IR Tools package and is distributed under the
% 3-Clause BSD License. A separate license file should be provided as part
% of the package.

% Set default values for options.
defaultopt = struct('x0','none', 'MaxIter',100, 'x_true','none', ...
    'NoiseLevel','none', 'eta',1.01, 'NE_Rtol',1e-12, 'IterBar','on', ...
    'NoStop', 'off');

% If input is 'defaults,' return the default options in X.
if nargin==1 && nargout <= 1 && isequal(A,'defaults')
    X = defaultopt;
    return;
end

defaultopt.verbosity = 'on';

% Check for acceptable number of optional input arguments.
switch length(varargin)
    case 0
        K = []; options = [];
    case 1
        if isa(varargin{1}, 'double')
            K = varargin{1}; options = [];
        else
            K = []; options = varargin{1};
        end
    case 2
        if isa(varargin{1}, 'double')
            K = varargin{1}; options = varargin{2};
        else
            K = varargin{2}; options = varargin{1};
        end
        if isfield(options, 'MaxIter') && ~isempty(options.MaxIter) && (~isempty(K) && options.MaxIter ~= max(K))
            warning('The value of MaxIter is discarded; the maximum value in K is taken as MaxIter')
        end
    otherwise
        error('Too many input parameters')
end

if isempty(options)
    options = defaultopt;
end

options = IRset(defaultopt, options);

MaxIter    = IRget(options, 'MaxIter',    [], 'fast');
x_true     = IRget(options, 'x_true',     [], 'fast');
NoiseLevel = IRget(options, 'NoiseLevel', [], 'fast');
eta        = IRget(options, 'eta',        [], 'fast');
NE_Rtol    = IRget(options, 'NE_Rtol',    [], 'fast');
IterBar    = IRget(options, 'IterBar',    [], 'fast');
NoStop     = IRget(options, 'NoStop',     [], 'fast');
verbose    = IRget(options, 'verbosity',  [], 'fast');

verbose = strcmp(verbose, 'on');

% Setting K.
if isempty(K)
    K = MaxIter;
end
% Sorting the iterations (in case they are shuffled in input).
K = K(:); K = sort(K,'ascend'); K = unique(K);
if ~((isreal(K) && (all(K > 0)) && all(K == floor(K))))
    error('K must be a vector of positive real integers')
end
if K(end) ~= MaxIter
    MaxIter = K(end);
end

StopIt = MaxIter;

if isempty(NoiseLevel) || strcmp(NoiseLevel,'none')
    Rtol = 0;
else
    Rtol = eta*NoiseLevel;
end

% We need to find the number of columns in matrix A, but if A is not given
% as a matrix, and no initial guess is given, then we can find it by
% computing A'*b.  Since we need this anyway, it doesn't cost any
% additional work.
trAb = Atransp_times_vec(A, b);
n = length(trAb);

nrmb = norm(b(:));
nrmAtb = norm(trAb(:));

% See if an initial guess is given.  If not, then use 0 as initial guess.
x = IRget(options, 'x0', [], 'fast');

if strcmp(x,'none')
    % the default initial guess for the iterations is defined as the
    % minimizer of || b - A*(alpha*ones(n,1)) ||_2
    coeffx0 = A_times_vec(A, ones(n,1)); alpha = (coeffx0'*b)/norm(coeffx0)^2;
    if alpha <= 0
        alpha = sqrt(eps);
    end
    x = alpha*ones(n,1);
end

if norm(x) == 0
    error(['IRmrnsd cannot handle a zero initial guess. ',...
    'Please consider assigning a different initial guess, or avoid specifying the initial guess in order to have a default value'])
end

% If initial guess has negative values, compensate.
minx = min(x(:));
if minx < 0
     x = x - min(0,minx) + sqrt(eps);
end

noIterBar = strcmp(IterBar,{'off'});

% Declare matrices.
X = zeros(n,length(K));
Xnrm    = zeros(MaxIter,1);
Rnrm    = zeros(MaxIter,1);
NE_Rnrm = zeros(MaxIter,1);
if strcmp(x_true,'none')
    errornorms = false;
else
    errornorms = true;
    Enrm = zeros(MaxIter,1);
    nrmtrue = norm(x_true(:));
    BestReg.It = [];
    BestReg.X = [];
    BestReg.Enrm = [];
    BestEnrm = 1e10;
    BestReg.Xnrm = [];
    BestReg.Rnrm = [];
    BestReg.NE_Rnrm = [];
end

NoStop = strcmp(NoStop,'on');
saved_iterations = zeros(1, length(K));

% Initializing some variables.
r = b - A_times_vec(A,x);
g = -Atransp_times_vec(A, r);
xg = x .* g;
gamma = g(:)' * xg(:);

if ~noIterBar
  h_wait = waitbar(0, 'Running iterations, please wait ...');
end

j = 0;
for k = 1:MaxIter
    if ~noIterBar
        waitbar(k/MaxIter, h_wait)
    end

    s = - x .* g;

    u = A_times_vec(A, s);

    theta = gamma / (u(:)'*u(:));
    neg_ind = s < 0;

    alpha = min( theta, min( -x(neg_ind) ./ s(neg_ind) ) );
    if isempty(alpha)
        alpha = theta;
    end

    x = x + alpha*s;

    r = r - alpha*u;
    AlreadySaved = 0;
    if any(K == k)
        j = j+1;
        X(:,j) = x;
        saved_iterations(j) = k;
        AlreadySaved = 1;
    end

    z = Atransp_times_vec(A, u);

    g = g + alpha*z;
    xg = x .* g;
    gamma = g(:)' * xg(:);

    % Compute norms.
    Xnrm(k)    = norm(x(:));
    Rnrm(k)    = norm(r(:))/nrmb;
    NE_Rnrm(k) = sqrt(gamma)/nrmAtb;
    if errornorms
        Enrm(k) = norm(x_true-x)/nrmtrue;
        if Enrm(k)<BestEnrm
            BestReg.It = k;
            BestReg.X = x;
            BestEnrm = Enrm(k);
            BestReg.Enrm = BestEnrm;
            BestReg.Xnrm = Xnrm(k);
            BestReg.Rnrm = Rnrm(k);
            BestReg.NE_Rnrm = NE_Rnrm(k);
        end
    end
    if (Rnrm(k) <= Rtol)  && (StopIt == MaxIter)
        if verbose
            disp('Residual tolerance satisfied')
        end
        StopFlag = 'Residual tolerance satisfied';
        StopIt = k;
        StopReg.It = k;
        StopReg.X = x;
        if errornorms, StopReg.Enrm = Enrm(k); end
        if ~ NoStop
            if ~AlreadySaved
                j = j+1;
                X(:,j) = x;
                saved_iterations(j) = k;
                AlreadySaved = 1;
            end
            Xnrm    = Xnrm(1:k);
            Rnrm    = Rnrm(1:k);
            NE_Rnrm = NE_Rnrm(1:k);
            if errornorms, Enrm = Enrm(1:k); end
            X = X(:,1:j);
            saved_iterations = saved_iterations(1:j);
            break
        end
    end
    if NE_Rnrm(k) <= NE_Rtol
        if verbose
            disp('Normal equations residual tolerance satisfied')
        end
        StopFlag = 'Normal equations residual tolerance satisfied';
        StopIt = k;
        StopReg.It = k;
        StopReg.X = x;
        if errornorms, StopReg.Enrm = Enrm(k); end
        if ~ NoStop
            if ~AlreadySaved
                j = j+1;
                X(:,j) = x;
                saved_iterations(j) = k;
                AlreadySaved = 1;
            end
            Xnrm    = Xnrm(1:k);
            Rnrm    = Rnrm(1:k);
            NE_Rnrm = NE_Rnrm(1:k);
            if errornorms, Enrm = Enrm(1:k); end
            X = X(:,1:j);
            saved_iterations = saved_iterations(1:j);
            break
        end
    end
end
if k == MaxIter
  if StopIt == MaxIter
    % Stop because max number of iterations reached
    if verbose
        disp('Reached maximum number of iterations')
    end
    StopFlag = 'Reached maximum number of iterations';
    StopReg.It = k;
    StopReg.X = x;
    if errornorms, StopReg.Enrm = Enrm(k); end
    if ~AlreadySaved
        j = j+1;
        X(:,j) = x;
        saved_iterations(j) = k;
    end
    Xnrm    = Xnrm(1:k);
    Rnrm    = Rnrm(1:k);
    NE_Rnrm = NE_Rnrm(1:k);
    if errornorms, Enrm = Enrm(1:k); end
    X = X(:,1:j);
    saved_iterations = saved_iterations(1:j);
  end
end
if ~noIterBar, close(h_wait), end
if nargout==2
  info.its = k;
  info.saved_iterations = saved_iterations(1:j);
  info.StopFlag = StopFlag;
  info.StopReg = StopReg;
  info.Rnrm = Rnrm(1:k);
  info.NE_Rnrm = NE_Rnrm(1:k);
  info.Xnrm = Xnrm(1:k);
  if errornorms
    info.Enrm = Enrm(1:k);
    info.BestReg = BestReg;
  end
end

repo: IRtools


text:
function varargout  =   vect2var(varargin)
%VECT2VAR Unpack vector to variables
%
%     [x1, x2, x3,...] =   vect2var(X) create variables x1, x2, x3... from
%     the vector X in the current workspace such as x1 = X(1), x2 = X(2),
%     x3 = X(3)...
%     vect2var is usefull with the optimization toolbox in which unknown
%     variables can be grouped into a single row vector. vect2var allows
%     the use of common names for variables with a minimum waste of time.
%
%     Example :
%     X     =   [3E5    20   5   0.898];
%     [b f jh ws] =   vect2var(X)
%     Returns :
%       b =
%             300000
%       f =
%           20
%       jh =
%            5
%       ws =
%           0.8980
%
%     See also DEAL, MMV2STRUCT, LISTS

%   vect2var is complementary of mmv2struct developped by D.C. Hanselman
%   Author : Meille Christophe
%   Pharmacokinetics Department, UFR Pharmacy, Marseille France
%   Version 1.07    Date: 2005/11/18

X   =   varargin{1};% Extract vector X
lgt =   nargout;    % Extract number of outputs
if ndims(X)>2       % Error handling
    error('Input vector required.')
elseif length(X)<lgt
    error('Input vector shorter than number of outputs.')
else
    varargout = cell(1,lgt);% preallocating cell array
    for i=1:lgt;% index of outputs
        varargout{i}    =   X(i); % unpack
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function [ crange ] = SpecColorRange( spec,range )
%[ crange ] = SpecColorRange( spec,range ) calculates a good color range for a
%spectrogram and sets the plot to the good range.
%
%
%DLevenstein 2017
%%
if ~exist('range','var')
    range = [2.5 2.5];
end

[~,mu,sig] = zscore(spec);
crange = [min(mu)-range(1)*max(sig) max(mu)+range(2)*max(sig)];
caxis(crange)

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni

% Discrete Fourier Transform  properties

% Parseval's Theorem

n=0:10;
x=1./(n+1);
En=sum(abs(x).^2)

N=length(x);
X=dft(x);
Edft=(1/N)*sum(abs(X).^2)


text:
function [net, info] = cnn_imagenet_se_mcn(varargin)
% CNN_IMAGENET_SE_MCN Evaluate imported PyTorch models on ImageNet val set
% (closely based on the mcn cnn_imagenet.m example)
%
% Copyright (C) 2017 Samuel Albanie
% Licensed under The MIT License [see LICENSE.md for details]

  opts.gpus = 3 ;
  opts.continue = 1 ;
  opts.batchSize = 256 ;
  opts.model = 'SE-ResNet-50-mcn' ;
  opts.modelDir = fullfile(vl_rootnn, 'data/models-import') ;
  opts.dataDir = fullfile(vl_rootnn, 'data/datasets/ILSVRC2012') ;
  opts.labelMap = fullfile(vl_rootnn, 'contrib/mcnSENets/misc/label_map.txt') ;
  [opts, varargin] = vl_argparse(opts, varargin) ;

  opts.expDir = fullfile(vl_rootnn, 'data', ['imagenet12-' opts.model]) ;
  [opts, varargin] = vl_argparse(opts, varargin) ;

  opts.lite = false ;
  opts.numFetchThreads = 12 ;
  opts.imdbPath = fullfile(vl_rootnn, 'data', 'imagenet12', 'imdb.mat');
  opts = vl_argparse(opts, varargin) ;

  opts.train.gpus = opts.gpus ;

% -------------------------------------------------------------------------
%                                                              Prepare data
% -------------------------------------------------------------------------

  if exist(opts.imdbPath, 'file')
    imdb = load(opts.imdbPath) ;
    imdb.imageDir = fullfile(opts.dataDir, 'images');
  else
    imdb = cnn_imagenet_setup_data('dataDir', opts.dataDir, 'lite', opts.lite) ;
    mkdir(opts.expDir) ;
    save(opts.imdbPath, '-struct', 'imdb') ;
  end

  % remap labels to match the order used in training
  imdb = updateLabelMap(imdb, opts) ;

% -------------------------------------------------------------------------
%                                                             Prepare model
% -------------------------------------------------------------------------
  net = load_model(opts.modelDir, opts.model) ;

  % modify the imdb to skip training images
  imdb.images.set(imdb.images.set == 1) = 4 ;
  [net, info] = cnn_train_dag(net, imdb, getBatchFn(opts, net.meta), ...
                              'expDir', opts.expDir, 'gpus', opts.train.gpus, ...
                              'numEpochs', 1, 'continue', opts.continue, ...
                              'batchSize', opts.batchSize) ;

  % hack the checkpoint file to save an empty dag
  modelPath = fullfile(opts.expDir, 'net-epoch-1.mat') ; tmp = load(modelPath) ;
  net_ = net ; net = dagnn.DagNN().saveobj() ; stats = tmp.stats ; state = {} ; %#ok
  save(modelPath, 'net', 'stats', 'state') ; net = net_ ;

% -------------------------------------------------------------------------
function dag = load_model(modelDir, name)
% -------------------------------------------------------------------------
  modelPath = fullfile(modelDir, sprintf('%s.mat', name)) ;
  if ~exist(modelDir, 'dir') , mkdir(modelDir) ; end

  if ~exist(modelPath, 'file')
    fprintf('Downloading the %s model ... this may take a while\n', name) ;
    base = 'http://www.robots.ox.ac.uk/~albanie' ;
    url = sprintf('%s/models/se-nets/%s.mat', base, name) ;
    urlwrite(url, modelPath) ;
  end

  dag = dagnn.DagNN.loadobj(load(modelPath)) ;
  dag.addLayer('softmax', dagnn.SoftMax(), dag.layers(end).outputs, 'prediction', {}) ;
  dag.addLayer('top1err', dagnn.Loss('loss', 'classerror'), ...
               {'prediction','label'}, 'top1err') ;
  dag.addLayer('top5err', dagnn.Loss('loss', 'topkerror', 'opts', {'topK',5}), ...
               {'prediction','label'}, 'top5err') ;

% -------------------------------------------------------------------------
function fn = getBatchFn(opts, meta)
% -------------------------------------------------------------------------
  bopts = struct('useGpu', numel(opts.train.gpus) > 0, ...
                 'imageSize', meta.normalization.imageSize(1:2)) ;
  fn = @(x,y) eval_get_batch(bopts,x,y) ;

% -------------------------------------------------------------------------
function varargout = eval_get_batch(opts, imdb, batch)
% -------------------------------------------------------------------------
  images = strcat([imdb.imageDir filesep], imdb.images.name(batch)) ;
  data = getImageBatch(images, opts, 'prefetch', nargout == 0) ;
  if nargout > 0
    labels = imdb.images.label(batch) ;
    varargout{1} = {'data', data, 'label', labels} ;
  end

% ----------------------------------------
function imdb = updateLabelMap(imdb, opts)
% ----------------------------------------
  labelMap = importdata(opts.labelMap) ;
  keep = imdb.images.label ~= 0 ;
  newLabels = labelMap(imdb.images.label(keep)) ;
  imdb.images.label(keep) = newLabels ;

% -------------------------------------------------
function data = getImageBatch(imagePaths, varargin)
% -------------------------------------------------
% GETIMAGEBATCH  Load and jitter a batch of images
opts.useGpu = false ;
opts.prefetch = false ;
opts.numThreads = 10 ;

% Options that were used during PyTorch training
% Note: that normalisation must occur after the pixel
% values have been rescaled to [0,1]
opts.cropSize = 224 / 256 ;
opts.imageSize = [224, 224] ;
opts.meanImg = [123, 117, 104] ;
%opts.std = [0.229, 0.224, 0.225] ;
opts = vl_argparse(opts, varargin);

args{1} = {imagePaths, ...
           'NumThreads', opts.numThreads, ...
           'Pack', ...
           'Interpolation', 'bilinear', ... % use bilinear to reproduce trainig resize
           'Resize', opts.imageSize(1:2), ...
           'CropSize', opts.cropSize, ...
           'CropAnisotropy', [1 1], ... % preserve aspect ratio
           'CropLocation', 'center'} ; % centre crop for testing

if opts.useGpu, args{end+1} = {'Gpu'} ; end
args = horzcat(args{:}) ;

if opts.prefetch
  vl_imreadjpeg(args{:}, 'prefetch') ;
  data = [] ;
else
  data = vl_imreadjpeg(args{:}) ;
  data = bsxfun(@minus, data{1}, permute(opts.meanImg, [1 3 2])) ;
  %data = data{1} / 255 ; % scale to (almost) [0,1]
  %data = bsxfun(@rdivide, data, permute(opts.std, [1 3 2])) ;
end

repo: MFT


text:
function K = covNNone(hyp, x, z, i)

% Neural network covariance function with a single parameter for the distance
% measure. The covariance function is parameterized as:
%
% k(x^p,x^q) = sf2 * asin(x^p'*P*x^q / sqrt[(1+x^p'*P*x^p)*(1+x^q'*P*x^q)])
%
% where the x^p and x^q vectors on the right hand side have an added extra bias
% entry with unit value. P is ell^-2 times the unit matrix and sf2 controls the
% signal variance. The hyperparameters are:
%
% hyp = [ log(ell)
%         log(sqrt(sf2) ]
%
% Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch, 2010-09-10.
%
% See also COVFUNCTIONS.M.

if nargin<2, K = '2'; return; end                  % report number of parameters
if nargin<3, z = []; end                                   % make sure, z exists
xeqz = isempty(z); dg = strcmp(z,'diag');                       % determine mode

n = size(x,1);
ell2 = exp(2*hyp(1));
sf2 = exp(2*hyp(2));

sx = 1 + sum(x.*x,2);
if dg                                                               % vector kxx
  K = sx./(sx+ell2);
else
  if xeqz                                                 % symmetric matrix Kxx
    S = 1 + x*x';
    K = S./(sqrt(ell2+sx)*sqrt(ell2+sx)');
  else                                                   % cross covariances Kxz
    S = 1 + x*z'; sz = 1 + sum(z.*z,2);
    K = S./(sqrt(ell2+sx)*sqrt(ell2+sz)');
  end
end

if nargin<4                                                        % covariances
  K = sf2*asin(K);
else                                                               % derivatives
  if i==1                                                          % lengthscale
    if dg
      V = K;
    else
      vx = sx./(ell2+sx);
      if xeqz
        V = repmat(vx/2,1,n) + repmat(vx'/2,n,1);
      else
        vz = sz./(ell2+sz); nz = size(z,1);
        V = repmat(vx/2,1,nz) + repmat(vz'/2,n,1);
      end
    end
    K = -2*sf2*(K-K.*V)./sqrt(1-K.*K);
  elseif i==2                                                        % magnitude
    K = 2*sf2*asin(K);
  else
    error('Unknown hyperparameter')
  end
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function temp = tempdim(temp,from,to)
%TEMPDIM  Convert temperature units
%
%   tempOut = TEMPDIM(tempIn, FROM, TO) converts tempIn from the units
%   specified by the string FROM to the units specified by the string
%   TO.  FROM and TO are case-insensitive, and may equal any of the
%   following:
%
%        'farenheit', 'far', or 'f'
%        'rankine', 'rank', or 'r'
%        'celsius', 'cel' or 'c'
%        'kelvin', 'kel', or 'k'
%
%   See also  FAR2CEL,  FAR2KEL,  FAR2RANK,
%             CEL2FAR,  CEL2KEL,  CEL2RANK,
%             KEL2FAR,  KEL2CEL,  KEL2RANK,
%             RANK2FAR, RANK2CEL, RANK2KEL

error(nargchk(3, 4, nargin, 'struct'))

% Warn and convert to real if TEMP is complex.
temp = ignoreComplex(temp, mfilename, 'TEMP');

% Convert units only if there's something to change.
if ~strcmp(from, to)
    temp = applyconversion(temp, from, to);
end

function temp = applyconversion(temp, from, to)

from = lower(from);
to = lower(to);

toIsSupported   = true;
fromIsSupported = true;

switch from
    case {'farenheit','far','f'}
        switch to
            case {'rankine','rank','r'}
                temp = far2rank(temp);
            case {'celsius','cel','c'}
                temp = far2cel(temp);
            case {'kelvin', 'kel','k'}
                temp = far2kel(temp);
            otherwise
                toIsSupported = false;
        end
    case {'rankine','rank','r'}
        switch to
            case {'farenheit','far','f'}
                temp = rank2far(temp);
            case {'celsius','cel','c'}
                temp = rank2cel(temp);
            case {'kelvin', 'kel','k'}
                temp = rank2kel(temp);
            otherwise
                toIsSupported = false;
        end
    case {'celsius','cel','c'}
        switch to
            case {'rankine','rank','r'}
                temp = cel2rank(temp);
            case {'farenheit','far','f'}
                temp = cel2far(temp);
            case {'kelvin', 'kel','k'}
                temp = cel2kel(temp);
            otherwise
                toIsSupported = false;
        end
    case {'kelvin', 'kel','k'}
        switch to
            case {'farenheit','far','f'}
                temp = kel2far(temp);
            case {'celsius','cel','c'}
                temp = kel2cel(temp);
            case {'rankine','rank','r'}
                temp = kel2rank(temp);
            otherwise
                toIsSupported = false;
        end
    otherwise
        fromIsSupported = false;
end

assert(toIsSupported, 'map:distdim:UnsupportedToUnits', ...
    'Unsupported ''TO'' units: %s.', to)

assert(fromIsSupported, 'map:distdim:UnsupportedFromUnits', ...
    'Unsupported ''FROM'' units: %s.', from)


text:
% FIDPATH gets the path of Varian FID data

% Note that the filename one chooses in the path doesn't really
% matter since the Varian *.fid directory contains 4 files,
% and they are always named "fid" (FID data), "procpar" (MR
% parameters), "log" (acquisition on/off), and text (PSD name).

function fpath = fidpath(fpath, checkfid, title)

if nargin < 3
   title = 'Select Varian FID path';
   if nargin < 2
      checkfid = false;
      if nargin < 1
         fpath = cd;
      end
   end
end

if isempty(fpath)
   fpath = cd;
end
if isempty(checkfid)
   checkfid = false;
end
if isempty(title)
   title = 'Select Varian FID path';
end

% Make sure fpath ends with "/"
if (fpath(end) ~= '/')
   fpath = [fpath, '/'];
end

if fileexist(fpath, checkfid)
   return
end

% Check if the user forgot to add .fid extension
fpathnew = [fpath(1:end-1) '.fid/'];
if fileexist(fpathnew, checkfid)
   fpath = fpathnew;
   return
end

% Check if the user forgot to add .par extension
fpathnew = [fpath(1:end-1) '.par/'];
if fileexist(fpathnew, false)
   fpath = fpathnew;
   return
end

curdir = cd;
if exist(fpath, 'dir')
   cd(fpath);
end

% Get FID file path/name
contflag = 1;
while (contflag == 1),
    [fname, fpath] = uigetfile('*', title);
    % Make sure "fid" and "procpar" files exist
    if isequal(fpath, 0) || isequal(fname, 0)
        error('No directory selected.')
    elseif ~fileexist(fpath, checkfid)
        disp('GETPATH unable to locate fid and/or procpar file in selected directory. Try again.')
    else
        contflag = 0;
    end
end

cd(curdir);
drawnow;
return

% LOCAL FUNCTION
function res = fileexist(fpath, checkfid);

if checkfid
   if (exist([fpath 'fid']) ~= 2)
      res = 0;
      return
   end
end

if exist([fpath 'procpar']) == 2
   res = 1;
else
   res = 0;
end

return

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function results = vl_test_imsmooth(varargin)
% VL_TEST_IMSMOOTH
vl_test_init ;

function s = setup()
I = im2double(imread(fullfile(vl_root,'data','spots.jpg'))) ;
I = max(min(vl_imdown(I),1),0) ;
s.I = single(I) ;

function test_pad_by_continuity(s)
% Convolving a constant signal padded with continuity does not change
% the signal.
I = ones(3) ;
for ker = {'triangular', 'gaussian'}
  ker = char(ker) ;
  J  = vl_imsmooth(I, 2, ...
                   'kernel', ker, ...
                   'padding', 'continuity') ;
  vl_assert_almost_equal(J, I, 1e-4, ...
                         'padding by continutiy with kernel = %s', ker) ;
end

function test_kernels(s)
for ker = {'triangular', 'gaussian'}
  ker = char(ker) ;
  for type = {@single, @double}
    for simd = [0 1]
      for sigma = [1 2 7]
        for step = [1 2 3]
          vl_simdctrl(simd) ;
          conv = type{1} ;
          g = equivalent_kernel(ker, sigma) ;
          J  = vl_imsmooth(conv(s.I), sigma, ...
                           'kernel', ker, ...
                           'padding', 'zero', ...
                           'subsample', step) ;
          J_ = conv(convolve(s.I, g, step)) ;
          vl_assert_almost_equal(J, J_, 1e-4, ...
                                 'kernel=%s sigma=%f step=%d simd=%d', ...
                                 ker, sigma, step, simd) ;
        end
      end
    end
  end
end

function g = equivalent_kernel(ker, sigma)
switch ker
  case 'gaussian'
    W = ceil(4*sigma) ;
    g = exp(-.5*((-W:W)/(sigma+eps)).^2) ;
  case 'triangular'
    W = max(round(sigma),1) ;
    g = W - abs(-W+1:W-1) ;
end
g = g / sum(g) ;

function I = convolve(I, g, step)
if strcmp(class(I),'single')
  g = single(g) ;
else
  g = double(g) ;
end
for k=1:size(I,3)
  I(:,:,k) = conv2(g,g,I(:,:,k),'same');
end
I = I(1:step:end,1:step:end,:) ;

author: Cloud-CV

repo: object-proposals


text:
function [  ] = plotAxialSagittalCoronal( vol, scale_fig, title_fig )
%PLOTAXIALSAGITTALCORONAL Summary of this function goes here
%
%   Code refractored from Berkin Bilgic's scripts: "script_Laplacian_unwrap_Sharp_Fast_TV_gre3D.m"
%   and "script_Laplacian_unwrap_Sharp_Fast_TV_gre3D.m"
%   Original source: https://martinos.org/~berkin/software.html
%
%   Original reference:
%   Bilgic et al. (2014), Fast quantitative susceptibility mapping with
%   L1-regularization and automatic parameter selection. Magn. Reson. Med.,
%   72: 1444-1459. doi:10.1002/mrm.25029
im_size = round(size(vol)/2);

figure(), subplot(1,3,1), imagesc( vol(:,:,im_size(3)), scale_fig  ), axis square off, colormap gray, axis image
figure(get(gcf,'Number')), subplot(1,3,3), imagesc( imrotate(squeeze(vol(:,im_size(2),:)),90), scale_fig  ), axis square off, colormap gray, axis image
figure(get(gcf,'Number')), subplot(1,3,2), imagesc( imrotate(squeeze(vol(im_size(1),:,:)),90), scale_fig  ), axis square off, colormap gray, axis image

if nargin == 3
    title(title_fig)
end

drawnow

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function yval = r8poly2_val2 ( dim_num, ndata, tdata, ydata, left, tval )

%*****************************************************************************80
%
%% R8POLY2_VAL2 evaluates a parabolic interpolant through tabular data.
%
%  Discussion:
%
%    This routine is a utility routine used by OVERHAUSER_SPLINE_VAL.
%    It constructs the parabolic interpolant through the data in
%    3 consecutive entries of a table and evaluates this interpolant
%    at a given abscissa value.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of a single data point.
%    DIM_NUM must be at least 1.
%
%    Input, integer NDATA, the number of data points.
%    NDATA must be at least 3.
%
%    Input, real TDATA(NDATA), the abscissas of the data points.
%    The values in TDATA must be in strictly ascending order.
%
%    Input, real YDATA(DIM_NUM,NDATA), the data points
%    corresponding to the abscissas.
%
%    Input, integer LEFT, the location of the first of the three
%    consecutive data points through which the parabolic interpolant
%    must pass.  1 <= LEFT <= NDATA - 2.
%
%    Input, real TVAL, the value of T at which the parabolic
%    interpolant is to be evaluated.  Normally, TDATA(1) <= TVAL <= T(NDATA),
%    and the data will be interpolated.  For TVAL outside this range,
%    extrapolation will be used.
%
%    Output, real YVAL(DIM_NUM), the value of the parabolic
%    interpolant at TVAL.
%

%
%  Check.
%
  if ( left < 1 || ndata-2 < left )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8POLY2_VAL2 - Fatal error!\n' );
    fprintf ( 1, '  LEFT < 1 or NDATA-2 < LEFT.\n' );
    fprintf ( 1, '  LEFT = %d\n', left );
    error ( 'R8POLY2_VAL2 - Fatal error!' );
  end

  if ( dim_num < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8POLY2_VAL2 - Fatal error!\n' );
    fprintf ( 1, '  DIM_NUM < 1.\n' );
    fprintf ( 1, '  DIM_NUM = %d\n', dim_num );
    error ( 'R8POLY2_VAL2 - Fatal error!' );
  end
%
%  Copy out the three abscissas.
%
  t1 = tdata(left);
  t2 = tdata(left+1);
  t3 = tdata(left+2);

  if ( t2 <= t1 || t3 <= t2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8POLY2_VAL2 - Fatal error!\n' );
    fprintf ( 1, '  T2 <= T1 or T3 <= T2.\n' );
    fprintf ( 1, '  T1 = %f\n', t1 );
    fprintf ( 1, '  T2 = %f\n', t2 );
    fprintf ( 1, '  T3 = %f\n', t3 );
    error ( 'R8POLY2_VAL2 - Fatal error!' );
  end
%
%  Construct and evaluate a parabolic interpolant for the data
%  in each dimension.
%
  for i = 1 : dim_num

    y1 = ydata(i,left);
    y2 = ydata(i,left+1);
    y3 = ydata(i,left+2);

    dif1 = ( y2 - y1 ) / ( t2 - t1 );
    dif2 = ( ( y3 - y1 ) / ( t3 - t1 ) ...
           - ( y2 - y1 ) / ( t2 - t1 ) ) / ( t3 - t2 );

    yval(i) = y1 + ( tval - t1 ) * ( dif1 + ( tval - t2 ) * dif2 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function bas = bas3DP1(p,t)

%% USAGE: coefficient of Lagrange P1 basis for tetrahedral mesh
%
% INPUTS:
% p --- np-by-3 vector
% t --- nt-by-4 vector
%
% OUTPUTS:
% bas --- nt-4-4 matrix stores the coefficient of basis value(s)
%         The k-i-j th entry denotes the j-th basis function on k-th cell
%         with the i-th coefficient where c_i (i=1,2,3,4) stands for
%                 c_1 + c_2*x + c_3*y + c_4*z

% Last Modified: 08/07/2020 by Xu Zhang

%%
x1 = p(t(:,1),1); y1 = p(t(:,1),2); z1 = p(t(:,1),3);
x2 = p(t(:,2),1); y2 = p(t(:,2),2); z2 = p(t(:,2),3);
x3 = p(t(:,3),1); y3 = p(t(:,3),2); z3 = p(t(:,3),3);
x4 = p(t(:,4),1); y4 = p(t(:,4),2); z4 = p(t(:,4),3);

A = -x2.*y3.*z1 + x2.*y4.*z1 + x1.*y3.*z2 - x1.*y4.*z2 + x2.*y1.*z3 ...
    -x1.*y2.*z3 + x1.*y4.*z3 - x2.*y4.*z3 + x4.*(-y2.*z1 + y3.*z1 + ...
    y1.*z2 - y3.*z2 - y1.*z3 + y2.*z3) + (-x2.*y1 + x1.*y2 - x1.*y3 + ...
    x2.*y3).*z4 + x3.*(-y4.*z1 - y1.*z2 + y4.*z2 + y2.*(z1 - z4) + y1.*z4);

nt = size(t,1);  bas = zeros(nt,4,4);
bas(:,1,1) = (-x4.*y3.*z2 + x3.*y4.*z2 + x4.*y2.*z3 - x2.*y4.*z3 - x3.*y2.*z4 + x2.*y3.*z4);
bas(:,2,1) =-(-y3.*z2 + y4.*z2 + y2.*z3 - y4.*z3 - y2.*z4 + y3.*z4);
bas(:,3,1) = (-x3.*z2 + x4.*z2 + x2.*z3 - x4.*z3 - x2.*z4 + x3.*z4);
bas(:,4,1) =-(-x3.*y2 + x4.*y2 + x2.*y3 - x4.*y3 - x2.*y4 + x3.*y4);
bas(:,1,2) =-(-x4.*y3.*z1 + x3.*y4.*z1 + x4.*y1.*z3 - x1.*y4.*z3 - x3.*y1.*z4 + x1.*y3.*z4);
bas(:,2,2) = (-y3.*z1 + y4.*z1 + y1.*z3 - y4.*z3 - y1.*z4 + y3.*z4);
bas(:,3,2) =-(-x3.*z1 + x4.*z1 + x1.*z3 - x4.*z3 - x1.*z4 + x3.*z4);
bas(:,4,2) = (-x3.*y1 + x4.*y1 + x1.*y3 - x4.*y3 - x1.*y4 + x3.*y4);
bas(:,1,3) = (-x4.*y2.*z1 + x2.*y4.*z1 + x4.*y1.*z2 - x1.*y4.*z2 - x2.*y1.*z4 + x1.*y2.*z4);
bas(:,2,3) =-(-y2.*z1 + y4.*z1 + y1.*z2 - y4.*z2 - y1.*z4 + y2.*z4);
bas(:,3,3) = (-x2.*z1 + x4.*z1 + x1.*z2 - x4.*z2 - x1.*z4 + x2.*z4);
bas(:,4,3) =-(-x2.*y1 + x4.*y1 + x1.*y2 - x4.*y2 - x1.*y4 + x2.*y4);
bas(:,1,4) =-(-x3.*y2.*z1 + x2.*y3.*z1 + x3.*y1.*z2 - x1.*y3.*z2 - x2.*y1.*z3 + x1.*y2.*z3);
bas(:,2,4) = (-y2.*z1 + y3.*z1 + y1.*z2 - y3.*z2 - y1.*z3 + y2.*z3);
bas(:,3,4) =-(-x2.*z1 + x3.*z1 + x1.*z2 - x3.*z2 - x1.*z3 + x2.*z3);
bas(:,4,4) = (-x2.*y1 + x3.*y1 + x1.*y2 - x3.*y2 - x1.*y3 + x2.*y3);
bas = bas./A;

repo: ifem


text:
function gotcha_pr2cphd( directory_name, output_file, pass, pol, azimuth_range )
%GOTCHA_PR2CPHD Converts a dataset from the AFRL 2D/3D SAR Volumetric
%public release dataset to CPHD file format.
%
% Usage: gotcha_pr2cphd( directory_name, output_file, pass, pol, azimuth_range)
%
% Inputs: directory_name: Data root directory for the GOTCHA public release
%                         dataset.  Assumes all files are kept in the same
%                         directory structure that the dataset is
%                         distributed with.
%         output_file:    Name of CPHD file to save.
%         pass:           What pass to image (1-8).  Default = 1.
%         pol:            What polarization to image ('HH','HV','VH','VV').
%                         Default = 'HH'.
%         azimuth_range:  Two-element vector of the first and last azimuth
%                         angle (in degrees) to convert to CPHD format.
%                         Default = [1 360] (all available angles).
%
% This function is just a wrapper to CONVERT_TO_CPHD, but allows for
% simpler selection of the pass, polarization, and azimuth range desired.
%
% Author: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

% Setup default parameters
if nargin<3
    pass = 1;
end
if nargin<4
    pol = 'HH';
end
if nargin<5
    azimuth_range = [1 360];
end

% Calculate which pulse numbers go with which files in directory for
% desired azimuth range.
last_pulse = 0;
pulse_ranges = zeros(azimuth_range(2),2); % First and last pulse number in every file
for az = 1:azimuth_range(end)
    fname = fullfile(directory_name, sprintf('pass%d',pass),pol,...
        sprintf('data_3dsar_pass%d_az%03d_%s.mat',pass,az,pol));
    load(fname); % Loads data into a structure named 'data'
    num_pulses = size(data.fp,2);
    pulse_ranges(az,:) = last_pulse + [1 num_pulses];
    last_pulse = last_pulse + num_pulses;
end
first_pulse = pulse_ranges(azimuth_range(1),1);
last_pulse = pulse_ranges(azimuth_range(2),2);

convert_to_cphdx(fname,output_file,'pulse_indices',first_pulse:last_pulse);
% convert_to_cphd30(fname,output_file,first_pulse:last_pulse); % Old way

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
%% string get_file_extension(string)
% filePath - string
% ext - string
%
function [ext] = get_file_extension(fileName)
  %ext = fileName(length(fileName)-2:length(fileName));
  [~,~,ext] = fileparts(fileName);
  ext = ext(2:end);
end

author: andrewssobral

repo: lrslibrary


text:
function volvisGUI(x,y,z,v)
% volvisGUI Interactive volume visualization.
%
%   Example:
%   [x,y,z,v] = flow;
%   volvisGUI(x,y,z,v)

%   Copyright 2007 The MathWorks, Inc.

%% Initalize visualization
figure;
s = volumeVisualization(x,y,z,v);
s.addSlicePlane(s.xMin);

%% Add uicontrol
hSlider = uicontrol(...
    'Units','normalized', ...
    'Position',[.75 .05 .2 .05], ...
    'Style','slider', ...
    'Min',s.xMin, ...
    'Max',s.xMax, ...
    'Value',s.xMin, ...
    'Callback',@updateSliderPosition);

%%
    function updateSliderPosition(varargin)
        s.deleteLastSlicePlane();
        x = get(hSlider,'Value');
        s.addSlicePlane(x);
    end

end