text:
function mAP = rfcn_test_vid_multiframe(conf, imdb, roidb, varargin)
% Video-level testing using final D&T model and RPN proposals
% --------------------------------------------------------
% D&T implementation
% Modified from MATLAB  R-FCN (https://github.com/daijifeng001/R-FCN/)
% and Faster R-CNN (https://github.com/shaoqingren/faster_rcnn)
% Copyright (c) 2017, Christoph Feichtenhofer
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------

%% inputs
    ip = inputParser;
    ip.addRequired('conf',                              @isstruct);
    ip.addRequired('imdb',                              @isstruct);
    ip.addRequired('roidb',                             @isstruct);
    ip.addParamValue('net_def_file',    '', 			@isstr);
    ip.addParamValue('net_file',        '', 			@isstr);
    ip.addParamValue('cache_name',      '', 			@isstr);
    ip.addParamValue('suffix',          '',             @isstr);
    ip.addParamValue('ignore_cache',    false,          @islogical);
    ip.addParamValue('ab_fetch_path',    '',          @isstr);
    ip.addParamValue('visualize_every',   Inf,          @isscalar);
    ip.addParamValue('test_classes',    [] 			);
    ip.addParamValue('chain_boxes_forward',    false, 			@islogical);
    ip.addParamValue('time_stride',   1,          @isscalar);
    ip.addParamValue('thresh_tracks',   0.1,          @isscalar);

    ip.parse(conf, imdb, roidb, varargin{:});
    opts = ip.Results;
    visualize = false;

%%  set cache dir
    cache_dir = fullfile(conf.root_path, 'output', 'rfcn_cachedir', opts.cache_name, imdb.name);
    mkdir_if_missing(cache_dir);

%%  init log
    timestamp = datestr(datevec(now()), 'yyyymmdd_HHMMSS');
    mkdir_if_missing(fullfile(cache_dir, 'log'));
    log_file = fullfile(cache_dir, 'log', ['test_', timestamp, '.txt']);
    diary(log_file);

    num_images = length(imdb.image_ids);
    num_classes = imdb.num_classes;

    try
      aboxes = cell(num_classes, 1);
      if opts.ignore_cache
          throw('');
      end
      for i = 1:num_classes
        load(fullfile(cache_dir, [imdb.classes{i} '_boxes_' imdb.name opts.suffix]));
        aboxes{i} = boxes;
      end
         load(fullfile(cache_dir, ['tracks_' imdb.name opts.suffix]));
    catch
%%      testing
        % init caffe net
        caffe_log_file_base = fullfile(cache_dir, 'caffe_log');
        caffe.init_log(caffe_log_file_base);
        caffe_net = caffe.Net(opts.net_def_file, 'test');
        caffe_net.copy_from(opts.net_file);

        % set random seed
        prev_rng = seed_rand(conf.rng_seed);
        caffe.set_random_seed(conf.rng_seed);

        % set gpu/cpu
        if conf.use_gpu
            caffe.set_mode_gpu();
        else
            caffe.set_mode_cpu();
        end

        % determine the maximum number of rois in testing
        max_rois_num_in_gpu = 10000;

        disp('opts:');
        disp(opts);
        disp('conf:');
        disp(conf);

        %heuristic: keep an average of 160 detections per class per images prior to NMS
        max_per_set = 160 * num_images;
        % heuristic: keep at most 400 detection per class per image prior to NMS
        max_per_image = 400;
        % detection thresold for each class (this is adaptively set based on the max_per_set constraint)
        thresh = -inf * ones(num_classes, 1);
        % top_scores will hold one minheap of scores per class (used to enforce the max_per_set constraint)
        top_scores = cell(num_classes, 1);
        % all detections are collected into:
        %    all_boxes[cls][image] = N x 5 array of detections in
        %    (x1, y1, x2, y2, score)
        aboxes = cell(num_classes, 1);
        box_inds = cell(num_classes, 1);
        aboxes_track = cell(num_classes, 1);

        for i = 1:num_classes
            aboxes{i} = cell(length(imdb.image_ids), 1);
            box_inds{i} = cell(length(imdb.image_ids), 1);
        end
        ascores_track = aboxes{1};;
        aboxes_track = cell(length(imdb.image_ids), conf.nFramesPerVid);
        aboxes_fwd = aboxes;
        aboxes_bwd = aboxes;
        track_forward_temp = cell(length(imdb.image_ids), 1);
        if isfield(conf,'sample_vid') && conf.sample_vid
            [inds_val,frames_val] = video_generate_random_minibatch([], [], imdb, conf, false);
            frames_val = cellfun(@(x) padarray(x,conf.nFramesPerVid-mod(length(x),conf.nFramesPerVid),'symmetric','post'), frames_val, 'UniformOutput' , false); % sample first frames
        end
        if ~isfield(conf,'nFramesPerVid')
          conf.nFramesPerVid = 1;
        end
        t_stride = conf.nFramesPerVid;
        t_stride = opts.time_stride ;
        count = 0;
        t_start = tic;
        center_frame = ceil(conf.nFramesPerVid/2);
        for v = 1:numel(inds_val)
          for kk = 1:t_stride:numel(frames_val{v}) - conf.nFramesPerVid + 1
              sub_db_inds = frames_val{v}(kk:1:kk+conf.nFramesPerVid-1);
              d = roidb.rois(sub_db_inds);
              if ~isempty(opts.test_classes),
                testThis = false;
                for ii = 1:length(opts.test_classes)
                  if any(opts.test_classes(ii) == d(1).class), testThis = true; end
                end
                if ~testThis, continue; end
              end
              count = count + 1;
              fprintf('%s: test (%s) vid %d/%d frame %d/%d', procid(), imdb.name, v,numel(inds_val), kk,numel(frames_val{v}));
              th = tic;

              [image_roidb_val, ~, ~] = rfcn_prepare_image_roidb_offline(conf, imdb, roidb.rois(sub_db_inds) , conf.root_path, sub_db_inds);

              im = {}; boxes = {};
              for b = 1:numel(sub_db_inds)
                  im{end+1} = imread(image_roidb_val(b).image_path);
                  boxes{end+1} = image_roidb_val(b).boxes(~image_roidb_val(b).gt,:);
              end

              if opts.chain_boxes_forward && ~isempty(track_forward_temp{sub_db_inds(1)})
                boxes{1} = cat(1, boxes{1}, track_forward_temp{sub_db_inds(1)});
              end

              [boxes_frames, scores_frames, track_boxes] = rfcn_im_detect(conf, caffe_net, im, boxes, max_rois_num_in_gpu);
              fprintf(' time %.3fs\n', toc(th)/numel(im));

              if ~iscell(track_boxes), track_boxes = {track_boxes}; end;
              max_scores = sum(scores_frames{2}, 2);
              tracklets = (max_scores > opts.thresh_tracks) ;
              aboxes_track{sub_db_inds(2),1} = [aboxes_track{sub_db_inds(2),1}; track_boxes{1}(tracklets,1:4)];
              aboxes_track{sub_db_inds(2),2}  = [aboxes_track{sub_db_inds(2),2}; boxes_frames{2}(tracklets,1:4)];
              aboxes_track{sub_db_inds(2),3}  = [aboxes_track{sub_db_inds(2),3}; track_boxes{2}(tracklets,1:4)];
              ascores_track{sub_db_inds(2)} = [ascores_track{sub_db_inds(2)}; scores_frames{2}(tracklets,:)] ;

              track_forward_temp{sub_db_inds(2)}  =  boxes_frames{2}(tracklets,1:4);

              %% resort frame boxes
              count2 = 0;
              for i=sub_db_inds(:)'
                count2 = count2+1;
                boxes = boxes_frames{count2};  scores = scores_frames{count2}; im_f = im{count2};
                for j = 1:num_classes
                    inds = find(scores(:, j) > thresh(j));
                    if ~isempty(inds)
                        [~, ord] = sort(scores(inds, j), 'descend');
                        ord = ord(1:min(length(ord), max_per_image));
                        inds = inds(ord); inds_tmp = inds;
                        cls_boxes = boxes(inds, (1+(j-1)*4):((j)*4));
                        cls_scores = scores(inds, j);
                        aboxes{j}{i} = [aboxes{j}{i}; cat(2, single(cls_boxes), single(cls_scores))];
                        if ~isempty(box_inds{j}{i}), inds=  inds + numel(box_inds{j}{i}); end
                        box_inds{j}{i} = [box_inds{j}{i}; inds ];
                    else
                        aboxes{j}{i} = [aboxes{j}{i}; zeros(0, 5, 'single')];
                        box_inds{j}{i} = box_inds{j}{i};
                    end
                end

                if mod(count,opts.visualize_every) == 0
                 boxes_cell = cell(length(imdb.classes), 1); boxes_track = boxes_cell; boxes_track2 = boxes_cell;

                  vis_thres = 0.9;
                  for k = 1:length(boxes_cell)
                      boxes_cell{k} = [boxes(:, (1+(k-1)*4):(k*4)), scores(:, k)];
                      I = boxes_cell{k}(:, 5) >= vis_thres;
                      boxes_cell{k} = boxes_cell{k}(I, :);
                    if count2==2
                      I =  ascores_track{i,1}(:, k) >= vis_thres;
                      boxes_track{k} = [aboxes_track{i,1}(I,1:4), ascores_track{i,1}(I, k)];
                      boxes_track2{k} = [aboxes_track{i,3}(I,1:4), ascores_track{i,1}(I, k)];
                    end
                  end
                  subplot(2,3,count2);
                  showboxes(uint8(im{count2}(:,:,1:3)), boxes_cell, imdb.classes, 'default'); drawnow;
                  title(sprintf('result for frame: %s', image_roidb_val(count2).image_id)) ;

                  if (count2) == 2
                      subplot(2,3,4);

                      showboxes(uint8(im{1}(:,:,1:3)), boxes_track, imdb.classes, 'default'); drawnow;
                      title(sprintf('tracked frame 2 -> 1 : %s', image_roidb_val(count2).image_id)) ;
                      subplot(2,3,5);
                      showboxes(uint8(im{3}(:,:,1:3)), boxes_track2, imdb.classes, 'default'); drawnow;
                      title(sprintf('tracked frame 2 -> 3 : %s', image_roidb_val(count2).image_id)) ;

                  end
                end
              end
              if mod(count, 1000) == 0
                  for j = 1:num_classes
                  [aboxes{j}, box_inds{j}, thresh(j)] = ...
                      keep_top_k(aboxes{j}, box_inds{j}, i, max_per_set, thresh(j));
                  end
                  disp(thresh);
              end
          end
        end
        for j = 1:num_classes
            [aboxes{j}, box_inds{j}, thresh(j)] = ...
                keep_top_k(aboxes{j}, box_inds{j}, i, max_per_set, thresh(j));
        end
        disp(thresh);

        for i = 1:num_classes
            top_scores{i} = sort(top_scores{i}, 'descend');
            if (length(top_scores{i}) > max_per_set)
                thresh(i) = top_scores{i}(max_per_set);
            end

            % go back through and prune out detections below the found threshold
            for j = 1:length(imdb.image_ids)
                if ~isempty(aboxes{i}{j})
                    I = find(aboxes{i}{j}(:,end) < thresh(i));
                    aboxes{i}{j}(I,:) = [];
                    box_inds{i}{j}(I,:) = [];
                end

            end

            save_file = fullfile(cache_dir, [imdb.classes{i} '_boxes_' imdb.name opts.suffix]);
            boxes = aboxes{i};
            inds = box_inds{i};
            save(save_file, 'boxes', 'inds');
            clear boxes inds;
        end
        save_file = fullfile(cache_dir, ['tracks_' imdb.name opts.suffix]);
        save(save_file, 'aboxes_track', 'ascores_track');

        caffe.reset_all();
        rng(prev_rng);

   %% generate object tracks
    track_boxes = repmat({cell(size(aboxes{1}))},size(aboxes));

    [inds_val,frames_val] = video_generate_random_minibatch([], [], imdb, conf, false);
    paths = cell(numel(inds_val),numel(aboxes));

    for v = 1:numel(inds_val)
      tic_toc_print('generating tracks for video: %d/%d\n', v, length(inds_val));
      for c = 1:numel(aboxes)
        frameIdx = frames_val{v};

        frameBoxes =  aboxes{c}(frameIdx);
        nonempty_frames = (~cellfun(@isempty,frameBoxes));

        if any(nonempty_frames)
            frameBoxes = frameBoxes(nonempty_frames);
            frameIdx = frameIdx(nonempty_frames);
            paths{v,c} = make_tubes(frameBoxes, 25, true,{aboxes_track(frameIdx,:), ascores_track(frameIdx), c} );

            for j=1:numel(paths{v,c})
                for k=1:numel(frameIdx)
                    track_boxes{c}{frameIdx(k)} = [track_boxes{c}{frameIdx(k)}; [paths{v,c}(j).boxes(k,1:4) ...
                     paths{v,c}(j).scores(k) ] ] ;
                end
            end

        end
      end
    end

    end

    %%
    % ------------------------------------------------------------------------
    % Peform AP evaluation
    % ------------------------------------------------------------------------

    if isequal(imdb.eval_func, @imdb_eval_voc)
        new_parpool();
        parfor model_ind = 1:num_classes
          cls = imdb.classes{model_ind};
          res(model_ind) = imdb.eval_func(cls, aboxes{model_ind}, imdb, opts.cache_name, opts.suffix);
        end
    else
%       % sanity test:
%       boxes_gt = cell(num_classes,1);
%       classes = {roidb.rois.class};
%       boxes = {roidb.rois.boxes};
%       boxes = cellfun(@(x) cat(2,x,ones(size(x,1),1)),boxes,'UniformOutput' , false);
%       for c=1:num_classes
%         boxes_gt{c} = cell(num_images,1);
%       end
%       for j=1:num_images
%         if ~isempty(classes{j})
%           for c=unique(classes{j})'
%             boxes_gt{c}{j} = boxes{j}(classes{j}==c,:);
%           end
%         end
%       end
%       res_gt = imdb.eval_func(boxes_gt, imdb, conf, opts.suffix);
    % ilsvrc
    res = imdb.eval_func(aboxes, imdb, conf, opts.suffix);
    res = imdb.eval_func(track_boxes, imdb, conf, opts.suffix, false);
    res = imdb.eval_func(track_boxes, imdb, conf, opts.suffix, false, true);

    for v = 1:numel(inds_val)
        frameIdx = frames_val{v};
        for k=1:numel(frameIdx)
          im = imread(imdb.image_at(frameIdx(k)));
          boxes_cell = cellfun(@(x) x(frameIdx(k)), track_boxes);
          boxes_cell = cellfun(@(x) x(x(:,5) > .2,:), boxes_cell, 'UniformOutput' , false);

          figure(1);
          showboxes(uint8(im(:,:,1:3)), boxes_cell, imdb.classes, 'default'); drawnow;
          title(sprintf('tube result for frame: %s', imdb.image_ids{frameIdx(k)})) ;
          drawnow;
        end
      end
    end

    if ~isempty(res)
        fprintf('\n~~~~~~~~~~~~~~~~~~~~\n');
        fprintf('Results:\n');
        aps = [res(:).ap]' * 100;
        disp(aps);
        disp(mean(aps));
        fprintf('~~~~~~~~~~~~~~~~~~~~\n');
        mAP = mean(aps);
    else
        mAP = nan;
    end

    diary off;
end

% ------------------------------------------------------------------------
function [boxes, box_inds, thresh] = keep_top_k(boxes, box_inds, end_at, top_k, thresh)
% ------------------------------------------------------------------------
    % Keep top K
    X = cat(1, boxes{1:end_at});
    if isempty(X)
        return;
    end
    scores = sort(X(:,end), 'descend');
    thresh = scores(min(length(scores), top_k));
    for image_index = 1:end_at
        if ~isempty(boxes{image_index})
            bbox = boxes{image_index};
            keep = find(bbox(:,end) >= thresh);
            boxes{image_index} = bbox(keep,:);
            box_inds{image_index} = box_inds{image_index}(keep);
        end
    end
end

author: feichtenhofer

repo: Detect-Track


text:
%% Define params
s_cam_definitions;

%% Interpolate all data used for camino masks
mkdir(mask_dir);
files_to_interp = {roi_map_filename};%,ecc_map_filename,var_map_filename,seg_filename,t1_filename};
for ff = 1:length(files_to_interp)
    filename = files_to_interp{ff};
    %[foo1, xform_filename, ext, foo2] = fileparts(filename);
    %xform_filename = ['xform_', xform_filename, ext];
    camAlignImageToDTI( fullfile(dti_dir,b0_filename), ...
            fullfile(ret_dir,filename),...
            fullfile(mask_dir,filename));
end

%% Load all reinterpolated data
ni_seg = niftiRead(fullfile(mask_dir,seg_filename));
ni_var = niftiRead(fullfile(mask_dir,var_map_filename));
ni_ecc = niftiRead(fullfile(mask_dir,ecc_map_filename));
ni_rois = niftiRead(fullfile(mask_dir,roi_map_filename));

%% Gray Matter Mask
% RH
ni_gm = ni_seg;
ni_gm.fname = fullfile(mask_dir,'gm_rh.nii.gz');
% Grab only right hemisphere GM
ni_gm.data(ni_seg.data~=6) = 0;
ni_gm.data(ni_gm.data>0) = 1;
% Only GM within occipital lobe
ni_gm.data(:,occ_exit:end,:)=0;
% Write file
writeFileNifti(ni_gm);

% LH
ni_gm = ni_seg;
ni_gm.fname = fullfile(mask_dir,'gm_lh.nii.gz');
% Grab only right hemisphere GM
ni_gm.data(ni_seg.data~=5) = 0;
ni_gm.data(ni_gm.data>0) = 1;
% Only GM within occipital lobe
ni_gm.data(:,occ_exit:end,:)=0;
% Write file
writeFileNifti(ni_gm);

%% Do the White Matter mask
% RH
ni_wm = ni_seg;
ni_wm.fname = fullfile(mask_dir,'wm_rh.nii.gz');
% Initialize
ni_wm.data=zeros(size(ni_seg.data));
% Get the white matter of RH
ni_wm.data(ni_seg.data==4)=1;
% exit plane
ni_wm.data(:,occ_exit+1:end,:)=0;
%write the file
writeFileNifti(ni_wm);

% RH
ni_wm = ni_seg;
ni_wm.fname = fullfile(mask_dir,'wm_lh.nii.gz');
% Initialize
ni_wm.data=zeros(size(ni_seg.data));
% Get the white matter of LH
ni_wm.data(ni_seg.data==3)=1;
% exit plane
ni_wm.data(:,occ_exit+1:end,:)=0;
%write the file
writeFileNifti(ni_wm);

%% Also create med-res white matter mask for seeds
%camAlignImageToDTI( fullfile(dti_dir,b0_filename), fullfile(mask_dir,'wm_rh.nii.gz'), fullfile(mask_dir,'wm_rh_lowres.nii.gz'), [1 1 1]);
%camAlignImageToDTI( fullfile(dti_dir,b0_filename), fullfile(mask_dir,'wm_lh.nii.gz'), fullfile(mask_dir,'wm_lh_lowres.nii.gz'), [1 1 1]);

%% Do the White Matter exit plane
% RH
ni_exit = ni_seg;
ni_exit.fname = fullfile(mask_dir,'wm_rh_exit.nii.gz');
% Initialize
ni_exit.data=zeros(size(ni_seg.data));
% Get the white matter of RH
ni_exit.data(ni_seg.data==4)=1;
% exit plane
ni_exit.data(:,occ_exit+1:end,:)=0;
ni_exit.data(:,1:occ_exit-1,:)=0;
%write the file
writeFileNifti(ni_exit);

% RH
ni_exit = ni_seg;
ni_exit.fname = fullfile(mask_dir,'wm_lh_exit.nii.gz');
% Initialize
ni_exit.data=zeros(size(ni_seg.data));
% Get the white matter of LH
ni_exit.data(ni_seg.data==3)=1;
% exit plane
ni_exit.data(:,occ_exit+1:end,:)=0;
ni_exit.data(:,1:occ_exit-1,:)=0;
%write the file
writeFileNifti(ni_exit);

%% Occ Lobe Mask
% RH
ni_occ = ni_seg;
ni_occ.fname = fullfile(mask_dir,'occ_rh_mask.nii.gz');
ni_gm = niftiRead(fullfile(mask_dir,'gm_rh.nii.gz'));
% Grab only right hemisphere GM and WM
ni_occ.data(ni_seg.data~=6 & ni_seg.data~=4) = 0;
ni_occ.data(ni_occ.data>0) = 1;
% Make sure all the GM voxels are included
ni_occ.data(ni_gm.data>0) = 1;
% Only GM within occipital lobe
ni_occ.data(:,occ_exit+1:end,:)=0;
% Write file
writeFileNifti(ni_occ);

% LH
ni_occ = ni_seg;
ni_occ.fname = fullfile(mask_dir,'occ_lh_mask.nii.gz');
ni_gm = niftiRead(fullfile(mask_dir,'gm_lh.nii.gz'));
% Grab only left hemisphere GM and WM
ni_occ.data(ni_seg.data~=5 & ni_seg.data~=3) = 0;
ni_occ.data(ni_occ.data>0) = 1;
% Make sure all the GM voxels are included
ni_occ.data(ni_gm.data>0) = 1;
% Only GM within occipital lobe
ni_occ.data(:,occ_exit+1:end,:)=0;
% Write file
writeFileNifti(ni_occ);

%% Create a background mask in native diffusion space
%ni_b0 = niftiRead(fullfile(dti_dir,b0_filename));
%lowres_pixdim = ni_b0.pixdim;
%camAlignImageToDTI( fullfile(dti_dir,b0_filename), fullfile(mask_dir,'occ_rh_mask.nii.gz'), fullfile(mask_dir,'bgmask_rh.nii.gz'), lowres_pixdim);
%camAlignImageToDTI( fullfile(dti_dir,b0_filename), fullfile(mask_dir,'occ_lh_mask.nii.gz'), fullfile(mask_dir,'bgmask_lh.nii.gz'), lowres_pixdim);

%% Define outside occ lobe for exclusion purposes
% RH
ni_occ = niftiRead( fullfile(mask_dir,'occ_rh_mask.nii.gz') );
ni_not_occ = niftiRead(fullfile(mask_dir,'occ_rh_mask.nii.gz'));
ni_not_occ.fname = fullfile(mask_dir,'not_occ_rh_mask.nii.gz');
ni_not_occ.data(ni_occ.data==0) = 1;
ni_not_occ.data(ni_occ.data==1) = 0;
writeFileNifti(ni_not_occ);
% LH
ni_occ = niftiRead( fullfile(mask_dir,'occ_lh_mask.nii.gz') );
ni_not_occ = niftiRead(fullfile(mask_dir,'occ_lh_mask.nii.gz'));
ni_not_occ.fname = fullfile(mask_dir,'not_occ_lh_mask.nii.gz');
ni_not_occ.data(ni_occ.data==0) = 1;
ni_not_occ.data(ni_occ.data==1) = 0;
writeFileNifti(ni_not_occ);

%% Image that defines the occipital exit plane and GM as two ROIs
% RH
ni_occ_exit = niftiRead(fullfile(mask_dir,'wm_rh_exit.nii.gz'));
ni_occ_exit.fname = fullfile(mask_dir,'occ_rh_exit.nii.gz');
ni_gm = niftiRead(fullfile(mask_dir,'gm_rh.nii.gz'));
ni_occ_exit.data(ni_gm.data>0) = 2;
writeFileNifti(ni_occ_exit);
% Now the mask
ni_occ_exit_mask = ni_occ_exit;
ni_occ_exit_mask.fname = fullfile(mask_dir,'occ_rh_exit_mask.nii.gz');
ni_occ_exit_mask.data(ni_occ_exit_mask.data>0) = 1;
writeFileNifti(ni_occ_exit_mask);

% LH
ni_occ_exit = niftiRead(fullfile(mask_dir,'wm_lh_exit.nii.gz'));
ni_occ_exit.fname = fullfile(mask_dir,'occ_lh_exit.nii.gz');
ni_gm = niftiRead(fullfile(mask_dir,'gm_lh.nii.gz'));
ni_occ_exit.data(ni_gm.data>0) = 2;
writeFileNifti(ni_occ_exit);
% Now the mask
ni_occ_exit_mask = ni_occ_exit;
ni_occ_exit_mask.fname = fullfile(mask_dir,'occ_lh_exit_mask.nii.gz');
ni_occ_exit_mask.data(ni_occ_exit_mask.data>0) = 1;
writeFileNifti(ni_occ_exit_mask);

%% Create label map, each voxel in GM map has some meaningful number
% Sweep out various variance explained thresholds
var_thresh = 0.1:0.05:0.3;
for vv = var_thresh
    filename = fullfile(mask_dir,sprintf('atlas_%s_%.2f.nii.gz',var_atlas_root,vv));
    data = ni_var.data;
    data(ni_var.data>vv) = 1;
    data(data<1) = 0;
    dtiWriteNiftiWrapper(data,ni_var.qto_xyz,filename);
    filename = fullfile(mask_dir,sprintf('atlas_%s_%s_%.2f.nii.gz',ecc_atlas_root,var_atlas_root,vv));
    ecc = ni_ecc.data;
    ecc = ecc.*data;
    ecc(ecc>0 & ecc<=5) = 5;
    ecc(ecc>50 & ecc<=10) = 10;
    ecc(ecc>10) = 15;
    dtiWriteNiftiWrapper(ecc,ni_var.qto_xyz,filename);
end

% Ventral ROIs
data = zeros(ni_rois.dim);
data_nov1 = zeros(ni_rois.dim);
for ll = roi_ventral_lbls
    data(ni_rois.data==ll) = 1;
    if ll~=1 && ll~=2
        data_nov1(ni_rois.data==ll) = 1;
    end
end
filename = fullfile(mask_dir,sprintf('atlas_%s.nii.gz',roi_ventral_atlas_root));
dtiWriteNiftiWrapper(data,ni_rois.qto_xyz,filename);
filename = fullfile(mask_dir,sprintf('atlas_%s.nii.gz',roi_ventral_nov1_atlas_root));
dtiWriteNiftiWrapper(data_nov1,ni_rois.qto_xyz,filename);
% Ecc Ventral
var = 0.2;
filename = fullfile(mask_dir,sprintf('atlas_%s_%s_%.2f.nii.gz',ecc_atlas_root,var_atlas_root,var));
ecc = niftiRead(filename);
ecc.data(data==0) = 0;
filename = fullfile(mask_dir,sprintf('atlas_%s_%s_%s_%.2f.nii.gz',roi_ventral_atlas_root,ecc_atlas_root,var_atlas_root,var));
dtiWriteNiftiWrapper(ecc.data,ecc.qto_xyz,filename);

% Dorsal ROIs
data = zeros(ni_rois.dim);
data_nov1 = zeros(ni_rois.dim);
for ll = roi_dorsal_lbls
    data(ni_rois.data==ll) = 1;
    if ll~=1 && ll~=2
        data_nov1(ni_rois.data==ll) = 1;
    end
end
filename = fullfile(mask_dir,sprintf('atlas_%s.nii.gz',roi_dorsal_atlas_root));
dtiWriteNiftiWrapper(data,ni_rois.qto_xyz,filename);
filename = fullfile(mask_dir,sprintf('atlas_%s.nii.gz',roi_dorsal_nov1_atlas_root));
dtiWriteNiftiWrapper(data_nov1,ni_rois.qto_xyz,filename);
% Ecc Dorsal
var = 0.2;
filename = fullfile(mask_dir,sprintf('atlas_%s_%s_%.2f.nii.gz',ecc_atlas_root,var_atlas_root,var));
ecc = niftiRead(filename);
ecc.data(data==0) = 0;
filename = fullfile(mask_dir,sprintf('atlas_%s_%s_%s_%.2f.nii.gz',roi_dorsal_atlas_root,ecc_atlas_root,var_atlas_root,var));
dtiWriteNiftiWrapper(ecc.data,ecc.qto_xyz,filename);

% Visual Field Cluster ROIs

% C1
data = zeros(ni_rois.dim);
for ll = roi_c1_lbls
    data(ni_rois.data==ll) = 1;
end
filename = fullfile(mask_dir,sprintf('atlas_%s.nii.gz',roi_c1_atlas_root));
dtiWriteNiftiWrapper(data,ni_rois.qto_xyz,filename);
filename = fullfile(mask_dir,sprintf('atlas_%s_%s_%.2f.nii.gz',ecc_atlas_root,var_atlas_root,var));
% C1 ECC
var = 0.2;
ecc = niftiRead(filename);
ecc.data(data==0) = 0;
filename = fullfile(mask_dir,sprintf('atlas_%s_%s_%s_%.2f.nii.gz',roi_c1_atlas_root,ecc_atlas_root,var_atlas_root,var));
dtiWriteNiftiWrapper(ecc.data,ecc.qto_xyz,filename);

% CMT
data = zeros(ni_rois.dim);
for ll = roi_cmt_lbls
    data(ni_rois.data==ll) = 1;
end
filename = fullfile(mask_dir,sprintf('atlas_%s.nii.gz',roi_cmt_atlas_root));
dtiWriteNiftiWrapper(data,ni_rois.qto_xyz,filename);
filename = fullfile(mask_dir,sprintf('atlas_%s_%s_%.2f.nii.gz',ecc_atlas_root,var_atlas_root,var));
% CMT ECC
var = 0.2;
ecc = niftiRead(filename);
ecc.data(data==0) = 0;
filename = fullfile(mask_dir,sprintf('atlas_%s_%s_%s_%.2f.nii.gz',roi_cmt_atlas_root,ecc_atlas_root,var_atlas_root,var));
dtiWriteNiftiWrapper(ecc.data,ecc.qto_xyz,filename);

% CVO
data = zeros(ni_rois.dim);
for ll = roi_cvo_lbls
    data(ni_rois.data==ll) = 1;
end
filename = fullfile(mask_dir,sprintf('atlas_%s.nii.gz',roi_cvo_atlas_root));
dtiWriteNiftiWrapper(data,ni_rois.qto_xyz,filename);
filename = fullfile(mask_dir,sprintf('atlas_%s_%s_%.2f.nii.gz',ecc_atlas_root,var_atlas_root,var));
% CVO ECC
var = 0.2;
ecc = niftiRead(filename);
ecc.data(data==0) = 0;
filename = fullfile(mask_dir,sprintf('atlas_%s_%s_%s_%.2f.nii.gz',roi_cvo_atlas_root,ecc_atlas_root,var_atlas_root,var));
dtiWriteNiftiWrapper(ecc.data,ecc.qto_xyz,filename);

% cv3ab
data = zeros(ni_rois.dim);
for ll = roi_cv3ab_lbls
    data(ni_rois.data==ll) = 1;
end
filename = fullfile(mask_dir,sprintf('atlas_%s.nii.gz',roi_cv3ab_atlas_root));
dtiWriteNiftiWrapper(data,ni_rois.qto_xyz,filename);
filename = fullfile(mask_dir,sprintf('atlas_%s_%s_%.2f.nii.gz',ecc_atlas_root,var_atlas_root,var));
% cv3ab ECC
var = 0.2;
ecc = niftiRead(filename);
ecc.data(data==0) = 0;
filename = fullfile(mask_dir,sprintf('atlas_%s_%s_%s_%.2f.nii.gz',roi_cv3ab_atlas_root,ecc_atlas_root,var_atlas_root,var));
dtiWriteNiftiWrapper(ecc.data,ecc.qto_xyz,filename);

%% Plot number of exit voxels vs. variance explained threshold
hem_names = {'lh','rh'};
var_thresh = 0.1:0.05:0.3;
exit_data = zeros(length(hem_names),length(var_thresh));
for hh = 1:2
    hem = hem_names{hh};
    % Load exit wm
    filename = fullfile(mask_dir,sprintf('wm_%s_exit.nii.gz',hem));
    ni = niftiRead(filename);
    num_vox_exit = sum(ni.data(:)>0);
    % Go through each variance explained level
    for vv = 1:length(var_thresh)
        var = var_thresh(vv);
        filename = fullfile(mask_dir,sprintf('%s_%s_%d_%s_%.2f_%s.nii.gz',cbs_root,hem,num_tracks,var_atlas_root,var,lbl_post));
        ni = niftiRead(filename);
        exit_data(hh,vv) = sum(ni.data(:)>0) / num_vox_exit;
    end
end

figure;
bar(var_thresh,exit_data');
legend(hem_names);
xlabel('Variance Explained');
title('Fraction of voxels connected to retinotopic occ. cortex');

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ y, m, d ] = day_borrow_eg_civil ( y, m, d )

%*****************************************************************************80
%
%% DAY_BORROW_EG_CIVIL borrows days from months in an Egyptian Civil date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, M, D, a year, month, and day
%    representing a date.  On input, D might be negative.  On output,
%    M should have decreased by one month, and D gone up by the
%    number of days in the month we "cashed in".  Y may be affected
%    if the input value of M was 1.
%
  while ( d <= 0 )

    m = m - 1;

    [ y, m ] = month_borrow_eg_civil ( y, m );

    days = month_length_eg_civil ( y, m );

    d = d + days;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [s,it]=display(c,tab,it)

  if nargin>1
    [s,it]=display(c.algorithm,tab,it);
    return;
  end

  disp(c.algorithm.name);
  disp(['data dimensions:'])

  global X; global Y;

  if isempty(c.myX)
    s=(['X = ' num2str(size(X(c.index,c.findex),1)) 'x' ...
	  num2str(size(X(c.index,c.findex),2))  ]);
  else
    s=(['X = ' num2str(size(c.myX,1)) 'x' num2str(size(c.myX,2))  ]);
  end

  if isempty(c.myY)
    if ~isempty(Y)
      disp([s '  Y = ' num2str(size(Y(c.index,:),1)) 'x' ...
	    num2str(size(Y(c.index,:),2))  ]);
    else
      disp([s '  Y = 0x0']);
    end
  else
    disp([s '  Y = ' num2str(size(c.myY,1)) 'x' num2str(size(c.myY,2))  ]);
  end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%% Edge Element Discretization of Maxwell Equations
% We test Maxwell solvers in iFEM.

clear all; close all
rowNames ={'h=1/2';'h=1/4';'h=1/8'};
colHeaders = {'H^curl Error','L^2 Error'};

%% The data of the pde
%
% * pde = Maxwelldata1; % zero Neumann boundary condition and curl u = 0
% * pde = Maxwelldata2; % non-homogenous Neumann boundary condition
% * pde = Maxwelldata3; % polynomial data and curl u = 0
% * pde = Maxwelldata4; % zero Dirichlet boundary condition
% * pde = Maxwelldata5; % linear polynomial data
% * pde = planewavedataC; % plane wave with complex coefficients
% * pde = planewavedata1; % plane wave with real coefficients

%% Positive Definite Case
% curl curl E + E = f.

help Maxwelldata2
%%
% Dirichlet boundary condition
[node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],0.5);
pde = Maxwelldata2;
bdFlag = setboundary3(node,elem,'Dirichlet');
option.solver = 'cg';
cubeMaxwell;
%%
makeHtmlTable([energyErr L2Err],[],rowNames,colHeaders);

%%
% Neumann boundary condition
[node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],0.5);
pde = Maxwelldata2;
bdFlag = setboundary3(node,elem,'Neumann');
option.solver = 'cg';
cubeMaxwell;
%%
makeHtmlTable([energyErr L2Err],[],rowNames,colHeaders);

%%
% Optimal first order of convergence is achieved. HX preconditioned CG
% converges around 20 steps.

%% Indefinite with real coefficients
% curl curl E - E = f.

help planewavedata1
%%
% Dirichlet boundary condition
[node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],0.5);
pde = planewavedata1;
bdFlag = setboundary3(node,elem,'Dirichlet');
option.solver = 'cg';
cubeMaxwell;
%%
makeHtmlTable([energyErr L2Err],[],rowNames,colHeaders);

%%
% Neumann boundary condition
[node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],0.5);
pde = planewavedata1;
bdFlag = setboundary3(node,elem,'Neumann');
option.solver = 'cg';
cubeMaxwell;
%%
makeHtmlTable([energyErr L2Err],[],rowNames,colHeaders);

%%
% Optimal first order of convergence is achieved. HX preconditioned CG
% converges around 40 steps although the system is indefinite.

%% Indefinite: complex coefficients, real solution
% curl curl E - (1-i)E = f.
option.solver = 'gmres';

help planewavedataC
%%
% Dirichlet boundary condition
[node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],0.5);
pde = planewavedataC;
bdFlag = setboundary3(node,elem,'Dirichlet');
cubeMaxwell;
%%
makeHtmlTable([energyErr L2Err],[],rowNames,colHeaders);

%%
% Neumann boundary condition
[node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],0.5);
pde = planewavedataC;
bdFlag = setboundary3(node,elem,'Neumann');
cubeMaxwell;
%%
makeHtmlTable([energyErr L2Err],[],rowNames,colHeaders);

%%
% Optimal first order of convergence is achieved. HX preconditioned GMRES
% converges around 90 steps although the system is indefinite.
%
% Bug: For Neumann boundary condition, the rate of L2 error is not quite right.

%% Indefinite: real coefficents, complex solution
option.solver = 'cg';

help planewavedata
%%
% Dirichlet boundary condition
[node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],0.5);
bdFlag = setboundary3(node,elem,'Dirichlet');
planewaveMaxwell;
%%
makeHtmlTable([energyErr L2Err],[],rowNames,colHeaders);

%%
% Neumann boundary condition
[node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],0.5);
bdFlag = setboundary3(node,elem,'Neumann');
planewaveMaxwell;
%%
makeHtmlTable([energyErr L2Err],[],rowNames,colHeaders);

%%
% The computation of error can't handle complex functions. So in
% planewaveMaxwell the error between uI and uh is computed. Therefore
% slightly better rate of convergence is observed. The rate of L2 error is
% almost second order.

repo: Aero-matlab


text:
function meanShiftImExplore( I, X, sigSpt, sigRng, show )
% Visualization to help choose sigmas for meanShiftIm.
%
% Displays the original image I, and prompts user to select a point on the
% image.  For given point, calculates the distance (both spatial and range)
% to every other point in the image.   It shows the results in a number of
% panes, which include 1) the original image I, 2) Srange - similarity
% based on range only, 3) Seuc - similarity based on Euclidean distance
% only, and 4) overall similarity.  Finally, in each image the green dot
% (possibly occluded) shows the original point, and the blue dot shows the
% new mean of the window after 1 step of meanShift.
%
% USAGE
%  meanShiftImExplore( I, X, sigSpt, sigRng, [show] )
%
% INPUTS
%  I       - MxN image for display
%  X       - MxNxP data array, P may be 1 (X may be same as I)
%  sigSpt  - integer specifying spatial standard deviation
%  sigRng  - value specifying the standard deviation of the range data
%  show    - [1] will display results in figure(show)
%
% OUTPUTS
%
% EXAMPLE
%  I=double(imread('cameraman.tif'))/255;
%  meanShiftImExplore( I, I, 5, .2, 1 );
%
% See also MEANSHIFTIM
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<5 ); show = 1; end
[mrows, ncols, p] = size(X);

%%% get input point
figure(show); clf; im(I);
[c,r] = ginput(1);
r=round(r); c=round(c);

%%% get D and S
[gridRs, gridCs] = ndgrid( 1:mrows, 1:ncols );
Deuc = ((gridRs-r).^2 + (gridCs-c).^2) / sigSpt^2;
x = X(r,c,:); x = x(:)';  Xflat = reshape(X,[],p);
Drange = pdist2( x, Xflat );
Drange = reshape( Drange, mrows, ncols ) / sigRng^2;
D = Drange + Deuc;

S = exp( -D );
Srange = exp( -Drange );
Seuc = exp( -Deuc );

%%% new c and r [stretched for display]
c2 = (gridCs .* S); c2 = sum( c2(:) ) / sum(S(:));
r2 = (gridRs .* S); r2 = sum( r2(:) ) / sum(S(:));
%c2 = c+(c2-c)*2; r2 = r+(r2-r)*2;

%%% show
figure(show); clf;
subplot(2,2,1); im(I);
hold('on'); plot( c, r, '.g' ); plot( c2, r2, '.b' ); hold('off');
subplot(2,2,2); im(Srange);
hold('on'); plot( c, r, '.g' ); plot( c2, r2, '.b' ); hold('off');
subplot(2,2,3); im(Seuc);
hold('on'); plot( c, r, '.g' ); plot( c2, r2, '.b' ); hold('off');
subplot(2,2,4); im(S);
hold('on'); plot( c, r, '.g' ); plot( c2, r2, '.b' ); hold('off');

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function [ vO ] = ConvolutionOverlapSave( vS, vK, convShape )
% ----------------------------------------------------------------------------------------------- %
% [ vO ] = ConvolutionDft( vS, vK, convShape )
% Applying 1D Linear Convolution using the Overlap and Save approach. The
% function calculates the optimal DFT Window.
% Input:
%   - vS                -   Input 1D Convolution Signal.
%                           Structure: Vector (signalLength, 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - vK                -   Input 1D Convolution Kernel.
%                           Structure: Vector (kernelLength, 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - convShape         -   Convolution Shape.
%                           The shape of the convolution which the output
%                           convolution matrix should represent. The
%                           options should match MATLAB's conv() function
%                           - Full / Same / Valid.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3}.
% Output:
%   - vS                -   Input 1D Convolution Output Vector.
%                           Structure: Vector (outputLength, 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References:
%   1.  MATLAB's 'conv()' - https://www.mathworks.com/help/matlab/ref/conv.html.
%   2.  Overlap and Save Method (Wikipedia) - https://en.wikipedia.org/wiki/Overlap%E2%80%93save_method.
% Remarks:
%   1.  The signals must be columns signals.
%   2.  It is assumed that the input signal is not shorter than the kernel.
% TODO:
%   1.  C
%   Release Notes:
%   -   1.0.000     27/04/2020  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

CONVOLUTION_SHAPE_FULL         = 1;
CONVOLUTION_SHAPE_SAME         = 2;
CONVOLUTION_SHAPE_VALID        = 3;

signalLength    = size(vS, 1); %<! K
kernelLength    = size(vK, 1); %<! M
dftLength       = CalcOptimalDftLength(signalLength, kernelLength); %<! N
convWinLength   = dftLength - kernelLength + 1; %<! L

paddSignalLength = ceil((signalLength + kernelLength - 1) / convWinLength) * dftLength;

vSS = [zeros(kernelLength - 1, 1); vS; zeros(paddSignalLength - kernelLength - 1 + signalLength, 1)];

numSteps    = ceil((signalLength + kernelLength - 1) / convWinLength);
vKD         = fft(vK, dftLength);
vO          = zeros(numSteps * convWinLength, 1);
vOO         = zeros(dftLength, 1);
idxPos      = 0;
for ii = 1:numSteps
    firstIdx = idxPos + 1;
    lastIdx = idxPos + dftLength;
    vOO(:) = ifft(fft(vSS(firstIdx:lastIdx)) .* vKD, 'symmetric');
    lastIdx = idxPos + convWinLength;
    vO(firstIdx:lastIdx) = vOO(kernelLength:dftLength);
    idxPos = idxPos + convWinLength;
end

switch(convShape)
    case(CONVOLUTION_SHAPE_FULL)
        idxFirst    = 1;
        idxLast     = signalLength + kernelLength - 1;
    case(CONVOLUTION_SHAPE_SAME)
        idxFirst    = 1 + floor(kernelLength / 2);
        idxLast     = idxFirst + signalLength - 1;
    case(CONVOLUTION_SHAPE_VALID)
        idxFirst    = kernelLength;
        idxLast     = (signalLength + kernelLength - 1) - kernelLength + 1;
end

vO          = vO(idxFirst:idxLast);

end

function [ dftLength ] = CalcOptimalDftLength( signelLength, kernelLength )

oConvOs = @(dftLength) (dftLength * log2(dftLength) + dftLength) / (dftLength - kernelLength + 1);

outputLength = signelLength + kernelLength;

firstPow2   = ceil(log2(kernelLength));
lastPow2    = ceil(log2(outputLength));

pow2 = firstPow2;
optNumOps = oConvOs(2 ^ pow2);

for pow2 = (firstPow2 + 1):lastPow2
    currNumOps = oConvOs(2 ^ pow2);
    if(currNumOps > optNumOps)
        break;
    end
    optNumOps = currNumOps;
end

dftLength = 2 ^ pow2;

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function[z]=celldiv(x,y)
%CELLDIV  Division acting on each element in a cell array.
%
%   Z=CELLDIV(X,Y) where X and Y are both cell arrays of N arrays, with the
%   corresponding elements in X and in Y having the same size, returns the
%   cell array Z containing their ratios,
%
%       Z{1}=Y{1}./X{1}, Z{2}=Y{2}./X{2},..., Z{N}=Y{N}./X{N}.
%
%   One of X or Y may also be a scalar or a numeric array of the same
%   length as the other input argument.
%
%   Usage: z=celldiv(x,y);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2015--2019 J.M. Lilly --- type 'help jlab_license' for details

if ~iscell(x)&&iscell(y)
    z=y;
    if length(x)==1
        x=x.*ones(length(y),1);
    end
    if length(x)==length(y)
        for i=1:length(x)
            z{i}=y{i}./x(i);
        end
    else
        error('LENGTH(X) and LENGTH(Y) must be the same.')
    end
elseif iscell(x)&&~iscell(y)
    z=x;
    if length(y)==1
        y=y.*ones(length(x),1);
    end
    if length(x)==length(y)
        for i=1:length(x)
            z{i}=y(i)./x{i};
        end
    else
        error('LENGTH(X) and LENGTH(Y) must be the same.')
    end
elseif iscell(x)&&iscell(y)
    z=y;
    for i=1:length(x)
        z{i}=y{i}./x{i};
    end
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [yd2] = ft22yd2(ft2)
% Convert area from square feet to square yards.
% Chad A. Greene 2012
yd2 = ft2/9;

author: Sable

repo: mcbench-benchmarks


text:
% ct_beta_study
% study how to adjust beta when pixel size changes in X-ray CT recon.
% results:
% changing nx has neglible effect, as expected
% changing pixel size has minimal effect when fwhm is converted to mm!
% values used CT studies:
% l2b=9, fov=500, fwhm = 1.62 mm
% l2b=9, fov=250, fwhm = 1.63 mm
% so we can  use the same beta regardless of pixel size and get nearly same fwhm
% caution: if other factors change, such as the # of views, may need to revisit!

%nx = 256;
%nx = 128;
nx = 64;
ny = nx;

if ~isvar('R')
	mask = true(nx,ny);
	l2b = 9;
	R = Robject(mask, 'edge_type', 'tight', 'beta', 2^l2b);
end

pixel_size_list = [500/512 250/512 0.25];
for ii=1:length(pixel_size_list)
	pixel_size = pixel_size_list(ii);

	sys = ct_sys('nx', nx, 'ny', ny, 'pixel_size', pixel_size);

	G = Gtomo2_dscmex(sys, 'nthread', 2);

%	psf = qpwls_psf(G, C, 2^l2b, mask, 1);
	psf = qpwls_psf(G, R, 1, mask, 1);
	fwhm = fwhm2(psf);
	printf('fwhm in mm = %g', fwhm * pixel_size)
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function obj = iatgetDemoObject
% IATGETDEMOOBJECT Returns the object used by the motion detection demo.
%
%    IATGETDEMOOBJECT finds the VIDEOINPUT object used by the motion detection
%    demo and returns it.
%
%    See also VIDEOINPUT.

% DT 4/2004
% Copyright 2004 The MathWorks, Inc.

% Look for existing objects created by the demo.
obj = imaqfind('Tag', 'Motion Detection Object');

if isempty(obj)
    obj = localCreateObj;
    stop(obj)
    % Clear the callbacks.
    obj.TimerFcn = '';
    obj.FramesAcquiredFcn = '';
    obj.StartFcn = '';
    obj.StopFcn = '';
    obj.TriggerFcn = '';

    % Configure the object for manual trigger so that GETSNAPSHOT returns
    % quicker.
    triggerconfig(obj, 'Manual');

    % Configure the object to return YCbCr data which the demo requires.
    set(obj, 'ReturnedColorSpace', 'YCbCr');

    % Save the video resolution in the base workspace.
    videoRes = obj.ROIPosition;
    assignin('base', 'videoRes', [videoRes(4) videoRes(3)]);
else
    obj = obj{1};

    % Save the video resolution in the base workspace.
    videoRes = obj.ROIPosition;
    assignin('base', 'videoRes', [videoRes(4) videoRes(3)]);
end

function obj = localCreateObj
% localCreateObj searches for the first supported device with a
% default format.  Once found, it creates an object for that device and
% returns the object to the calling function.  Some devices do not have
% default formats and can not be used since the demo doesn't provide a
% method of specifying camera files.

% Determine the available adaptors.
info = imaqhwinfo;

adaptorName = [];
adaptorID = [];
adaptorIndex = 1;

% Search through all of the installed adaptors until an available adaptor
% with hardware that has a default format is found.
while (isempty(adaptorName) && isempty(adaptorID))

    if (adaptorIndex > length(info.InstalledAdaptors))
        error('imaq:motiondetection:nohardware', 'No hardware compatible with the motion detection demo was found.');
    end

    % Get information on the current adaptor.
    curAdaptor = info.InstalledAdaptors{adaptorIndex};
    adaptorInfo = imaqhwinfo(curAdaptor);

    % For the current adaptor, search through each installed device looking
    % for a device that has a default format.
    idIndex = 1;
    while (isempty(adaptorID) && (idIndex <= length(adaptorInfo.DeviceIDs)))
        curID = adaptorInfo.DeviceIDs{idIndex};
        if ~isempty(adaptorInfo.DeviceInfo(idIndex).DefaultFormat)
            adaptorName = curAdaptor;
            adaptorID = curID;
        end

        % If the current device doesn't have a default format check the
        % next one.
        idIndex = idIndex + 1;
    end

    % Check the next adaptor.
    adaptorIndex = adaptorIndex + 1;
end

% Create the object.
obj = videoinput(adaptorName, adaptorID);

% Configure the object's tag so that it can be found later by using
% IMAQFIND.
set(obj, 'Tag', 'Motion Detection Object');

% End of localCreateObj


text:
%This is the Matlab interface to the OMP2 MEX implementation.
%The function is not for independent use, only through omp2.m.

%OMP2MEX Matlab interface to the OMP2 MEX implementation.
%  GAMMA = OMP2MEX(D,X,DtX,XtX,G,EPSILON,SPARSE_G,MSGDELTA,MAXATOMS,PROFILE)
%  invokes the OMP2 MEX function according to the specified parameters. Not
%  all the parameters are required. Those among D, X, DtX, XtX and G which
%  are not specified should be passed as [].
%
%  EPSILON - the target error.
%  SPARSE_G - returns a sparse GAMMA when nonzero, full GAMMA when zero.
%  MSGDELTA - the delay in secs between messages. Zero means no messages.
%  MAXATOMS - the max number of atoms per signal, negative for no max.
%  PROFILE - nonzero means that profiling information should be printed.

%  Ron Rubinstein
%  Computer Science Department
%  Technion, Haifa 32000 Israel
%  ronrubin@cs
%
%  April 2009

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function display(obj)
disp([num2str(obj.value) ' +/- ' num2str(obj.error)])


text:
function varargout = fillPolygon(varargin)
%FILLPOLYGON Fill a polygon specified by a list of points
%
%   fillPolygon(POLY);
%   Fills the interior of the polygon specified by POLY. The boundary of
%   the polygon is not drawn, see 'drawPolygon' to do it.
%   POLY is a single [N*2] array.
%   If POLY contains NaN-couples, each portion between the [NaN;NaN] will
%   be filled separately.
%
%   fillPolygon(PX, PY);
%   Specifies coordinates of the polygon in separate arrays.
%
%
%   H = fillPolygon(...);
%   Also returns a handle to the created patch
%
%
%   See also:
%     polygons2d, drawCurve, drawPolygon

%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 07/04/2005.
%

%   HISTORY
%   2008-05-07 add psb to specify drawing options
%   2008/10/15 add psb to draw polygons with holes

% check input
if isempty(varargin)
    error('need to specify a polygon');
end

% case of a set of polygons stored in a cell array
var = varargin{1};
if iscell(var)
    N = length(var);
    h = zeros(N, 1);
    for i = 1:N
        % check for empty polygons
        if ~isempty(var{i})
            h(i) = fillPolygon(var{i}, varargin{2:end});
        end
    end

    % setup output values
    if nargout > 0
        varargout{1} = h;
    end
    return;
end

% Extract coordinates of polygon vertices
if size(var, 2) > 1
    % first argument is a polygon array
    px = var(:, 1);
    py = var(:, 2);
    varargin(1) = [];
else
    % arguments 1 and 2 correspond to x and y coordinate respectively
    if length(varargin) < 2
        error('should specify either a N*2 array, or 2 N*1 vectors');
    end

    px = varargin{1};
    py = varargin{2};
    varargin(1:2) = [];
end

% Find position of breaks, and copies first point of each loop at the end
inds = find(isnan(px(:)));
i1 = [inds ; length(px)+1];
i0 = [1 ; inds+1];
px(i1, :) = px(i0, :);
py(i1, :) = py(i0, :);

% set default line format
if isempty(varargin)
    varargin = {'b'};
end

% fill the polygon with desired style
h = fill(px, py, varargin{:}, 'lineStyle', 'none');

% output
if nargout > 0
    varargout{1} = h;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [xc,mode] = XCorr1(x,y,varargin)

%XCorr1 - Compute cross-correlograms and their modes for pairs of curves.
%
%  USAGE
%
%    [xc,mode] = XCorr1(x,y,<options>)
%
%    x              first set of curves (MxN matrix: M curves, N bins)
%    y              second set of curves
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'type'        'linear' for linear curves (default), 'circular' otherwise
%    =========================================================================
%

% Copyright (C) 2012 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% Defaults
type = 'linear';

% Check number of parameters
if nargin < 2 | mod(length(varargin),2) ~= 0,
  error('Incorrect number of parameters (type ''help XCorr1'' for details).');
end

% Check parameter sizes
if ~isdmatrix(x) || ~isdmatrix(y),
	error('Both sets of curves should be MxN matrices (type ''help XCorr1'' for details).');
end
if ~all(size(x)==size(y)),
	error('Both sets of curves should have the same sizes (type ''help XCorr1'' for details).');
end

% Parse parameter list
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+2) ' is not a property (type ''help XCorr1'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'type',
			type = varargin{i+1};
			if ~isstring_FMAT(type,'linear','circular'),
				error('Incorrect value for property ''type'' (type ''help XCorr1'' for details).');
			end
		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help XCorr1'' for details).']);
	end
end

[m,n] = size(x);

if strcmp(type,'linear'),
	xc = nan(m,2*n-1);
	for i = 1:m,
		% Cross-correlagram i
		xci = xcorr(x(i,:),y(i,:),'coeff');
		% mode of cross-correlogram i (in [-1..1])
		mi = round(mean(find(xci==max(xci))));
		mi = mi/n-1;
		% Store in matrices
		xc(i,:) = xci;
		mode(i,1) = mi;
	end
else
	xc = nan(m,2*n-1);
	for i = 1:m,
		% Cross-correlagram i
		xci = xcorr(x(i,:),y(i,:),'coeff');
		% mode of cross-correlogram i (in [-1..1])
		mi = round(mean(find(xci==max(xci))));
		mi = mi/n-1;
		% For circular curves, very large modes must be corrected
		if abs(mi) > 0.5,
			if mi > 0,
				xci = [xci(n+1:end) xci(1:n)];
				mi = mi - 1;
			else
				xci = [xci(n:end) xci(1:n-1)];
				mi = mi +1;
			end
		end
		xc(i,:) = xci;
		mode(i,1) = mi;
	end
end

%  	xc = nan(m,2*n-1);
%  	for i = 1:m,
%  		% Cross-correlagram i
%  		X = repmat(x(i,:),n,1);
%  		Y = repmat(y(i,:),1,n+1);
%  		Y(n+1:n+1:end) = [];
%  		Y = reshape(Y,n,[]);
%  		xci = sum(X.*Y,2)';
%  		xci = [xci(2:end) xci];
%  		% mode of cross-correlogram i (in [0..1])
%  		% (we need to compute mean values of x, but x is circular, so we transform it to angles)
%  		a = find(xci==max(xci));
%  		a = exp(j*a*2*pi/n);
%  		mi = round(wrap(angle(mean(a)),2)/(2*pi)*n)/n-1;
%  		if abs(mi) > 0.5,
%  			if mi > 0,
%  				mi = mi - 1;
%  			else
%  				mi = mi +1;
%  			end
%  		end
%  		% Zero outside primary field
%  		m = round((mi+1)*n);
%  		h = ceil(n/2);
%  		xci(m+h:end) = 0;
%  		xci(1:m-h) = 0;
%  		% Store in matrices
%  		xc(i,:) = xci;
%  		mode(i,1) = mi;
%  	end
%  	% Scale so that autocorrelation at shift 0 is 1
%  	x0 = sum(x.^2,2);
%  	y0 = sum(y.^2,2);
%  	scale = sqrt(x0.*y0);
%  	xc = xc./repmat(scale,1,2*n-1);
%
%  	end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% test for level set on meshes

path(path, '../toolbox_graph_data/off/');
if not(exist('name'))
name = 'elephant-50kv';
name = 'bunny';
name = 'david50kf';
name = 'hand';
end

options.name = name;
[vertex,face] = read_off([name '.off']);

func = 'xaxis';
func = 'eigen';

switch func
    case 'xaxis'
        F = rescale(vertex(1,:))';
        tau = linspace(0.05,.95, 10);
    case 'eigen'
        options.symmetrize = 1;
        options.normalize = 1;
        disp('--> Computing Laplacian');
        L = compute_mesh_laplacian(vertex,face,'conformal',options);
        opts.disp = 0;
        p = 150;
        disp('--> Extracting Eigenvectors');
        [U,D] = eigs(L,p, 'SM', opts);
        D = diag(abs(D)); [D,I] = sort(D); U = U(:,I);
        nb = 10; sel = round(linspace(10,p,nb));
        F = U(:, sel);
        tau = 0;
end

rep = 'results/levelsets-meshes/';
if not(exist(rep))
    mkdir(rep);
end

disp('--> Displaying');
for it=1:size(F,2)
    progressbar(it,size(F,2));
    f = F(:,it);
    options.niter_averaging = 3;
    f = perform_mesh_smoothing(face, vertex, f, options);
    % extract level set
    [v1,v2] = compute_levelset_mesh(vertex,face,f,tau,options);

    % display
    lw = 3;
    fw = perform_histogram_matching(f, linspace(0,1,length(f)));
    options.face_vertex_color = fw;
    clf;
    hold on;
    plot_mesh(vertex,face,options);
    for i=1:size(v1,2)
        h = plot3( [v1(1,i) v2(1,i)], [v1(2,i) v2(2,i)], [v1(3,i) v2(3,i)], 'k' );
        set(h, 'LineWidth', lw);
    end
    hold off;
    colormap jet(256);
    saveas(gcf, [rep name '-' func '-' num2string_fixeddigit(it,2), '.png'], 'png');
end

author: gpeyre

repo: matlab-toolboxes


text:
s=figure;
set(s,'DefaultTextFontSize',10)
set(s,'DefaultAxesFontSize',10)
set(s,'DefaultAxesFontName','Times New Roman')
set(s,'DefaultTextFontName','Times New Roman')
set(s,'Units','inches')
pos=get(s,'Position');
set(s,'Position',[pos(1) pos(2) 6 4])
set(s,'PaperType','a4letter')


text:
%%*********************************************************************
%% This is an example of using the interface
%% to solve the following theta problem of a graph
%% with adjacency matrix G:
%% max <E,X>
%% s.t. trace(X)==1 ; X positive semidefinite
%%      X(i,j) == 0 for (i,j) in Edge_set
%% SDPNAL+:
%% Copyright (c) 2017 by
%% Yancheng Yuan, Kim-Chuan Toh, Defeng Sun and Xinyuan Zhao
%%*********************************************************************

 load theta6.mat
 [IE,JE] = find(triu(G,1));
 n = length(G);
 %%
 model = ccp_model('Example_theta');
     X = var_sdp(n, n);
     model.add_variable(X);
     model.maximize(sum(X));
     model.add_affine_constraint(trace(X) == 1);
     options=2;
     if (options==1)
        model.add_affine_constraint(X(IE,JE) == 0);
        model.add_affine_constraint(X >= 0);
     else %% another way to specify the constraints
        U = 1e20*(1-spones(G+G'));
        model.add_affine_constraint(0 <= X <= U);
     end
 model.solve;
%%*********************************************************************

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function prefs = editPreferences;
% GUI to set mrVista preference variables.
%
%  prefs = editPreferences;
%
% Returns a prefs struct with each preference field and its current value.
% This is the same struct that would be returned if you ran the command:
%	prefs = getpref('VISTA');
%
% NOTE ABOUT PREFERENCES: The preferences edited here are part of the
% "preferences group" named 'VISTA'; they are only a subset of
% the preferences which may affect mrVista's behavior. For instance,
% the function mrmPreferences will edit a separate set of preferences under
% the 'mesh' preferences group, which determines behavior for dealing with
% surface meshes.
%
% TODO: sub-group the preferences according to related ideas, such as 'File
% Paths', 'Display', 'ROIs', etc. Build a more general GUI to allow the
% user to switch between these preference groups.
%
% ras, 07/2008.

%% get the current values of the preferences
prefs = getpref('VISTA');

%% build the dialog
% (we only build dialogs for certain items; some, for instance, are only
% relevant to mrVista2 GUI behavior.)
dlg(1).fieldName = 'verbose';
dlg(end).style = 'popup';
dlg(end).string = 'How much feedback when using mrVista?';
dlg(end).list = {'0 -- quiet, minimal feedback' ...
				 '1 -- moderate, most status updates, waitbars shown' ...
				 '2 -- verbose, highest level of feedback'};
dlg(end).value = prefsVerboseCheck + 1;

dlg(end+1).fieldName = 'fileFormat';
dlg(end).style = 'popup';
dlg(end).string = 'Save maps / anatomies in which format?';
dlg(end).list = {'default' 'nifti'};
dlg(end).value = prefsFormatCheck;

% this preference only needed for MATLAB <2007a
% (later versions always use JAVA)
v = ver('Matlab');
if str2num(v.Version) < 7.4
	dlg(end+1).fieldName = 'javaOn';
	dlg(end).style = 'checkbox';
	dlg(end).string = 'Use JAVA to render figures';
	if isfield(prefs, dlg(end).fieldName)
		dlg(end).value = prefs.(dlg(end).fieldName);
	else
		dlg(end).value = 0;
	end
end

dlg(end+1).fieldName = 'defaultAnatomyPath';
dlg(end).style = 'edit';
dlg(end).string = 'Default Anatomy Directory (if ''3DAnatomy'' link not present)';
dlg(end).value = prefsDefaultAnatPath;

dlg(end+1).fieldName = 'defaultROIPath';
dlg(end).style = 'edit';
dlg(end).string = 'Default Shared ROI Directory (relative to anatomy directory)';
if isfield(prefs, dlg(end).fieldName)
	dlg(end).value = prefs.(dlg(end).fieldName);
else
	dlg(end).value = 'ROIs';
end

xHairList = {'Regular Crosshairs' 'BrainVoyager style (gap near intersection)'};
dlg(end+1).fieldName = 'xHairMethod';
dlg(end).style = 'popup';
dlg(end).string = 'Method for rendering crosshairs (3-view)?';
dlg(end).list = xHairList;
if isfield(prefs, dlg(end).fieldName)
	dlg(end).value = prefs.(dlg(end).fieldName);
else
	dlg(end).value = 2;
end

dlg(end+1).fieldName = 'xHairColor';
dlg(end).style = 'number';
dlg(end).string = 'Crosshairs color (3-view)?';
if isfield(prefs, dlg(end).fieldName)
	dlg(end).value = prefs.(dlg(end).fieldName)(:)';
else
	dlg(end).value = [1 .5 .5]; % [.8 .8 .6];
end
dlg(end+1).fieldName = 'extendedGrayFields';
dlg(end).style = 'checkbox';
dlg(end).string = 'Gray views: Extended ''leftNodes''/''rightNodes'' fields';
if isfield(prefs, dlg(end).fieldName)
	dlg(end).value = prefs.(dlg(end).fieldName);
else
	dlg(end).value = 1;
end

dlg(end+1).fieldName = 'autoComputeG2VMap';
dlg(end).style = 'checkbox';
dlg(end).string = 'Auto-compute Gray->Vertex Map for Meshes';
if isfield(prefs, dlg(end).fieldName)
	dlg(end).value = prefs.(dlg(end).fieldName);
else
	dlg(end).value = 1;
end

dlg(end+1).fieldName = 'autoComputeV2GMap';
dlg(end).style = 'checkbox';
dlg(end).string = 'Auto-compute Vertex->Gray Map for Meshes';
if isfield(prefs, dlg(end).fieldName)
	dlg(end).value = prefs.(dlg(end).fieldName);
else
	dlg(end).value = 1;
end

dlg(end+1).fieldName = 'savePrefs';
dlg(end).style = 'checkbox';
dlg(end).string = 'Save GUI settings when you close a view window';
if isfield(prefs, dlg(end).fieldName)
	dlg(end).value = prefs.(dlg(end).fieldName);
else
	dlg(end).value = 1;
end

%% put up the dialog; get response
[resp ok] = generalDialog(dlg, 'Set mrVista Preferences', [.25 .25 .5 .5]);

if ~ok
	if prefs.verbose >= 1
		disp('Edit Preferences Canceled.')
		return
	end
end

% parse response for some fields (most should already be in the final
% format)
resp.verbose = cellfind(dlg(1).list, resp.verbose) - 1;
resp.xHairMethod = cellfind(xHairList, resp.xHairMethod);

%% set the preferences based on the dialog
for f = fieldnames(resp)'
	setpref( 'VISTA', f{1}, resp.(f{1}) );
end

% return the updated prefs
prefs = getpref('VISTA');

return


text:
function varargout=wrsamp(varargin)
%
% wrsamp(tm,data,fileName,Fs,gain,format)
%
%    Wrapper to WFDB WRSAMP:
%         http://www.physionet.org/physiotools/wag/wrsamp-1.htm
%
% Writes data into a WFDB compatible *.dat and *.hea files. The files will be saved at
% the current directory (if the record is in the current directory) or, if a using
% a PhysioNet web record , a subdirectory in the current directory, with
% the relative path determined by recordName. You can use RDSAMP to
% verify that the write was completed successfully (see example below). For
% detailed information on the parameters described below, please see:
% http://www.physionet.org/physiotools/wag/header-5.htm
%
%
%
% LIMITATIONS:
%        The input variables 'tm' and 'data' should be integer values.
%		 Currently the MATLAB/Octave wrappers do take into account Baseline values.
%        Signals that have baseline values will not be converted properly to physical units.
%        You will have to edit the generated header files manually to include baseline variables and update the checksums.
%
%
% Required Parameters:
%
% tm
%       Nx1 vector of integers representing sample index.
%
% data
%       NxM vector of integers (M channels) that will be written to file.
%
% fileName
%
%       String specifying the file name. WRSAMP will write the signal file in
%       the current directory as "fileName.dat", and create a header file in the
%       current directory "fileName.hea" for the specified record.
%
% Optional Parameters:
%
% Fs
%       A 1x1 double specifying the sampling frequency (in samples per second per
%       signal) for the output signals (default: 250). It affects the output header
%       file only. This option has no effect on the output signal file, which
%       contains one sample per signal for each line of input.
%
% gain
%       A 1x1 or Mx1 vector of doubles specifying gain value in A/D units per millivolt
%       for the output signals (default: 200).  This option the output header file only.
%       This option has no effect on the output signal file.
%
% format
%       String specifying the signal's default format (default: "16"). For
%       information on the available WFDB formats please see:
%       http://www.physionet.org/physiotools/wag/signal-5.htm
%
%%Example-
%%Read signal in raw units
%[signal,Fs,tm]=rdsamp('challenge/2013/set-a/a01',[],[],[],4);
%[siginfo,Fs]=wfdbdesc('challenge/2013/set-a/a01');
%%Write a copy to file
%wrsamp(tm,signal(:,1),'a01Copy',Fs(1),200,siginfo(1).Format)
%%Check that the signals match
%[signalCopy,Fs,tm]=rdsamp('a01Copy',[],[],[],4);
%err=sum(signalCopy ~= signal(:,1))
%
%
% Written by Ikaro Silva, 2013
% Last Modified: April 10, 2015
% Since 0.0.1
%
%
% See also rdsamp, wfdbdesc, mat2wfdb
%

%endOfHelp

javaWfdbExec=getWfdbClass('wrsamp');

%Set default pararamter values
inputs={'tm','data','fileName','Fs','gain','format'};
Fs=[];
gain=[];
format=[];
for n=1:nargin
    if(~isempty(varargin{n}))
        eval([inputs{n} '=varargin{n};']);
    end
end

wfdb_argument={'-c','-z','-o',fileName};

if(~isempty(Fs))
    wfdb_argument{end+1}='-F';
    wfdb_argument{end+1}=num2str(Fs);
end

if(~isempty(gain))
    wfdb_argument{end+1}='-G';
    wfdb_argument{end+1}=num2str(gain);
end

if(~isempty(format))
    wfdb_argument{end+1}='-O';
    wfdb_argument{end+1}=format;
end

del=repmat([' '],size(tm));
data=[num2str(tm) del num2str(data)];
javaWfdbExec.setArguments(wfdb_argument);
javaWfdbExec.execWithStandardInput(cellstr(data));

author: ikarosilva

repo: wfdb-app-toolbox


text:
function result = requestJavaAdapter(object)
%REQUESTJAVAADAPTER Support function for GUIDE

%   Copyright 1984-2006 The MathWorks, Inc.
%   $Revision: 1915 $

%%%%%%%  CAUTION                                           %%%%%%%%
%%%%%%%  This file is duplicated in both uitools and guide %%%%%%%%%
%%%%%%%  %TODO - determine if this functionality can be    %%%%%%%%%
%%%%%%%  broken out or replaced so this file doesn't exist %%%%%%%%%
%%%%%%%  in two places                                     %%%%%%%%%

  len = length(object);
  if len == 1
    if ishandle(object) && ~isjava(object)
      result = java(handle(object));
    else
      error('MATLAB:requestJavaAdapter:InvalidInput', '''requestJavaAdapter'' argument must be a handle list.');
    end
  else
    if ~isempty(object) && all(ishandle(object) & ~isjava(object))
      result = cell(len, 1);
      for i = 1:len
        result{i} = java(handle(object(i)));
      end
    else
      error('MATLAB:requestJavaAdapter:InvalidArgument', '''requestJavaAdapter'' argument must be a handle list.');
    end
  end

author: spm


text:
function pklt( noisy_file, outfile)
%
%  Implements a perceptually-motivated subspace algorithm  [1].
%
%
%  Usage:  pklt(noisyFile, outputFile)
%
%         infile - noisy speech file in .wav format
%         outputFile - enhanced output file in .wav format
%
%
%  Example call:  pklt('sp04_babble_sn10.wav','out_pklt.wav');
%
%  References:
%   [1] Jabloun, F. and Champagne, B. (2003). Incorporating the human hearing
%   	 properties in the signal subspace approach for speech enhancement. IEEE
%   	 Trans. on Speech and Audio Processing, 11(6), 700-708.
%
% Authors: Yi Hu and Philipos C. Loizou
%
% Copyright (c) 2006 by Philipos C. Loizou
% $Revision: 0.0 $  $Date: 10/09/2006 $
%-------------------------------------------------------------------------

if nargin<2
   fprintf('Usage: pklt(noisyfile.wav,outFile.wav) \n\n');
   return;
end

vad_thre= 1.2; % 2.378 is the value from Mittal's paper
mu_vad= 0.98; % mu to use in vad

[noisy_speech, Srate, NBITS]= wavread( noisy_file);
subframe_dur= 4;  %subframe length is 4 ms
len= floor( Srate* subframe_dur/ 1000);
P= len; % sub-frame length - for 8k sampling rate, len is 32
frame_dur= 32; % frame length in msecs
N= frame_dur* Srate/ 1000;
Nover2= N/ 2; % window overlap in 50% of frame size
K= N;
frame_window= hamming( N);
subframe_window= hamming( P);
eta_v= .08; % used for gain values calc.

% ====noise covariance matrix estimation
L120=floor( 120* Srate/ 1000);
% assume the first 120ms is noise only
noise= noisy_speech( 1: L120);

noise_autoc= xcorr( noise, len- 1, 'biased');
% from -(len- 1) to (len- 1)
% obtain the autocorrelation functions
Rn= toeplitz( noise_autoc( len: end));
% form a Toeplitz matrix to obtain the noise signal covariance matrix
bartlett_win= bartlett( 2* len- 1);
n_autoc_win= noise_autoc.* bartlett_win;

for k= 0: N- 1
    Phi_w( k+ 1)= n_autoc_win( P: 2*P-1)' * 2* ...
        cos( 2* pi* k* (0: P- 1)'/ N)- n_autoc_win( P);
end
Phi_w= Phi_w';
% make Phi_w column vector

n_start= 1;

Nframes= floor( length( noisy_speech)/ (N/ 2))- 1;  % number of frames
x_overlap= zeros( Nover2, 1);

%===============================  Start Processing =====================

for n=1: Nframes

    noisy= noisy_speech( n_start: n_start+ N- 1);
    noisy_autoc= xcorr( noisy, len- 1, 'biased');
    Ry= toeplitz( noisy_autoc( len: 2* len- 1));

    % Use simple VAD algorithm to update noise cov matrix, Rn
    %
    vad_ratio= Ry(1,1)/ Rn(1,1);
    if (vad_ratio<= vad_thre) % noise dominant
        Rn= mu_vad* Rn+ (1- mu_vad)* Ry;
        noise_autoc_sc= Rn( 1, :)';
        % single sided noise autocorrelation
        noise_autoc= [flipud( noise_autoc_sc( 2: end)); ...
            noise_autoc_sc];
        n_autoc_win= noise_autoc.* bartlett_win;
        % compute Phi_w
        for k= 0: N- 1
            Phi_w( k+ 1)= n_autoc_win( P: 2*P-1)' * 2* ...
                cos( 2* pi* k* (0: P- 1)'/ N)- n_autoc_win( P);
        end
        % Phi_w is column vector
    end
    % =================

    Rx= Ry- Rn;

    [U, D]= eig( Rx);
    dD= diag( D); % retrieving diagonal elements
    dD_Q= find( dD> 0); % index for those eigenvalues greater than 0
    Lambda= dD( dD_Q);
    U1= U( :, dD_Q);
    % eigenvector for those eigenvalues greater than 0

    U1_fft= fft( U1, N);
    V= abs( U1_fft).^ 2;
    Phi_B= V* Lambda/ P;

    %==calculating masking threshold
    Phi_mask= mask( Phi_B( 1: N/ 2+ 1), N, Srate, NBITS);
    Phi_mask= [Phi_mask; flipud( Phi_mask( 2: N/ 2))];

    Theta= V'* Phi_mask/ K;
    Ksi= V'* Phi_w/ K;

    gain_vals= exp( -eta_v* Ksi./ min( Lambda, Theta));
    G= diag( gain_vals);
    H= U1* G* U1';

    % first step of synthesis for subframe
    sub_start= 1;
    sub_overlap= zeros( P/2, 1);
    for m= 1: (2*N/P- 1)
        sub_noisy= noisy( sub_start: sub_start+ P- 1);
        enhanced_sub_tmp= (H* sub_noisy).* subframe_window;
        enhanced_sub( sub_start: sub_start+ P/2- 1)= ...
            enhanced_sub_tmp( 1: P/2)+ sub_overlap;
        sub_overlap= enhanced_sub_tmp( P/2+1: P);
        sub_start= sub_start+ P/2;
    end
    enhanced_sub( sub_start: sub_start+ P/2- 1)= sub_overlap;

    xi= enhanced_sub'.* frame_window;
    xfinal( n_start: n_start+ Nover2- 1)= x_overlap+ xi( 1: Nover2);
    x_overlap= xi( Nover2+ 1: N);

    n_start= n_start+ Nover2;

end

xfinal( n_start: n_start+ Nover2- 1)= x_overlap;

wavwrite(xfinal, Srate, NBITS, outfile);

%=======================E N D ===================================

function M= mask( Sx, dft_length, Fs, nbits)
% Author: Patrick J. Wolfe
%         Signal Processing Group
%         Cambridge University Engineering Department
%         p.wolfe@ieee.org
% Johnston perceptual model initialisation
frame_overlap= dft_length/ 2;
freq_val = (0:Fs/dft_length:Fs/2)';
half_lsb = (1/(2^nbits-1))^2/dft_length;

freq= freq_val;
thresh= half_lsb;
crit_band_ends = [0;100;200;300;400;510;630;770;920;1080;1270;...
        1480;1720;2000;2320;2700;3150;3700;4400;5300;6400;7700;...
        9500;12000;15500;Inf];

% Maximum Bark frequency
%
imax = max(find(crit_band_ends < freq(end)));

% Normalised (to 0 dB) threshold of hearing values (Fletcher, 1929)
% as used  by Johnston.  First and last thresholds are corresponding
% critical band endpoint values, elsewhere means of interpolated
% critical band endpoint threshold values are used.
%
abs_thr = 10.^([38;31;22;18.5;15.5;13;11;9.5;8.75;7.25;4.75;2.75;...
        1.5;0.5;0;0;0;0;2;7;12;15.5;18;24;29]./10);
ABSOLUTE_THRESH = thresh.*abs_thr(1:imax);

% Calculation of tone-masking-noise offset ratio in dB
%
OFFSET_RATIO_DB = 9+ (1:imax)';

% Initialisation of matrices for bark/linear frequency conversion
% (loop increments i to the proper critical band)
%
num_bins = length(freq);
LIN_TO_BARK = zeros(imax,num_bins);
i = 1;
for j = 1:num_bins
    while ~((freq(j) >= crit_band_ends(i)) & ...
            (freq(j) < crit_band_ends(i+1))),
        i = i+1;
    end
    LIN_TO_BARK(i,j) = 1;
end

% Calculation of spreading function (Schroeder et al., 82)

spreading_fcn = zeros(imax);
summ = 0.474:imax;
spread = 10.^((15.81+7.5.*summ-17.5.*sqrt(1+summ.^2))./10);
for i = 1:imax
    for j = 1:imax
        spreading_fcn(i,j) = spread(abs(j-i)+1);
    end
end

% Calculation of excitation pattern function

EX_PAT = spreading_fcn* LIN_TO_BARK;

% Calculation of DC gain due to spreading function

DC_GAIN = spreading_fcn* ones(imax,1);

%Sx = X.* conj(X);

C = EX_PAT* Sx;

% Calculation of spectral flatness measure SFM_dB
%
[num_bins num_frames] = size(Sx);
k = 1/num_bins;
SFM_dB = 10.*log10((prod(Sx).^k)./(k.*sum(Sx)+eps)+ eps);

% Calculation of tonality coefficient and masked threshold offset
%
alpha = min(1,SFM_dB./-60);
O_dB = OFFSET_RATIO_DB(:,ones(1,num_frames)).*...
    alpha(ones(length(OFFSET_RATIO_DB),1),:) + 5.5;

% Threshold calculation and renormalisation, accounting for absolute
% thresholds

T = C./10.^(O_dB./10);
T = T./DC_GAIN(:,ones(1,num_frames));
T = max( T, ABSOLUTE_THRESH(:, ones(1, num_frames)));

% Reconversion to linear frequency scale

%M = 1.* sqrt((LIN_TO_BARK')*T);
M= LIN_TO_BARK'* T;

author: jtkim-kaist

repo: Speech-enhancement


text:
function newSegmentation = grow_region(this, startingSegmentation, threshold, nIter)
% Iterative region growing using starting segmentation and threshold.
%
%   Y = MrImage()
%   Y.grow_region(startingSegmentation, threshold, nIter)
%
% This is a method of class MrImage.
%
% IN
%   startingSegmentation    MrImage-object containing the starting
%                           segmentation, e.g. obtained via binarizing
%                           this
%   threshold               Threshold above which voxel connected to the
%                           object are considered part of the object.
%   nIter                   Number of maximum iterations. Maximum path
%                           length that can be grown.
%
% OUT
%
% EXAMPLE
%   Y.grow_region(Y.binarize(500), 300, 20);
%
%   See also MrImage

% Author:   Saskia Bollmann
% Created:  2020-03-30
% Copyright (C) 2020 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

% initialize
trueLocations = startingSegmentation.copyobj();

for n = 1:nIter

    % create searchLocations
    searchLocations = trueLocations.imdilate(strel('disk', 1));
    % exclude already found locations
    searchLocations = searchLocations - trueLocations;
    % apply in this
    searchVoxel = this .* searchLocations;
    % threshold
    newLocations = searchVoxel.binarize(threshold);
    % check new found locations
    disp([num2str(sum(newLocations.data(:))), ' new voxel(s) found.']);
    if ~any(newLocations.data(:))
        break
    end
    % add to true locations
    trueLocations = trueLocations + newLocations;

end
newSegmentation = trueLocations;

end

author: translationalneuromodeling

repo: tapas


text:
function L = compute_point_laplacian(pts, type, rings, options)

% compute_point_laplacian - compute a laplacian matrix L, L = -(Laplacian operator).
%
%   L = compute_point_laplacian(pts, type, rings, options)
%
%   If options.symmetrize=1 and options.normalize=0 then
%       L = D-W
%   If options.symmetrize=1 and options.normalize=1 then
%       L = eye(n)-D^{-1/2}*W*D^{-1/2}
%   If options.symmetrize=0 and options.normalize=1 then
%       L = eye(n)-D^{-1}*W.
%   where D=diag(sum(W,2)) and W is the unormalized weight matrix
%   (see compute_mesh_weight).
%
%   type can 'combinatorial', 'distance', 'conformal'.
%
%   See also compute_point_weight.
%
%   Changed default value for symmetrize from 1 to 0 (JJCAO)
%   Copyright (c) 2009 jjcao

options.null = 0;
normalize = getoptions(options, 'normalize', 0);
symmetrize = getoptions(options, 'symmetrize', 1);

W = compute_point_weight(pts, type, rings, options);
n = size(W,1);

if symmetrize==1 && normalize==0
    L = diag(sum(W,2)) - W;
elseif symmetrize==1 && normalize==1
    L = speye(n) - diag(sum(W,2).^(-1/2)) * W * diag(sum(W,2).^(-1/2));
elseif symmetrize==0 && normalize==1
    L = speye(n) - diag(sum(W,2).^(-1)) * W;
else
    error('Does not work with symmetrize=0 and normalize=0');
end

author: taiya

repo: cloudcontr

save_path: github-repos/MATLAB/taiya-cloudcontr


text:
function [xStateUpdate,PUpdate,innov,Pzz,W]=BLUEPolarMeasUpdateApprox(xStatePred,PPred,z,R)
%%BLUEPOLARMEASUPDATEAPPROX Perform the measurement update step in the
%                approximate best linear unbiased estimator (BLUE) for a
%                Cartesian state consisting of 2D position and velocity
%                when given a measurement in monostatic polar coordinates.
%
%INPUTS: xStatePred The 4X1 target state consisting of position components
%                   followed by velocity components:
%                   xStatePred=[x;y;xDot;yDot].
%             PPred The 3X3 covariance matrix associated with the predicted
%                   target state estimate.
%                 z A 2X1 one-way (monostatic) polar measurement (no
%                   refraction) with components ordered [range;azimuth]
%                   with azimuth in radians measured in radians,
%                   counterclockwise from the x-axis.
%                 R The 2X2 diagonal covariance matrix associated with the
%                   measurement z. Any cross terms present in R will be
%                   ignored.
%
%OUTPUTS: xStateUpdate The 4 X 1 updated state vector.
%              PUpdate The updated 4 X 4 state covariance matrix.
%           innov, Pzz The 2X1 innovation and the 2X2 innovation covariance
%                      matrix are returned in case one wishes to analyze
%                      the consistency of the estimator or use those values
%                      in gating or likelihood evaluation.
%                    W The gain used in the update. This can be useful
%                      when gating and using the function
%                      calcMissedGateCov.
%
%The algorithm is taken from Section IV of [1]. One equation that is not
%directly written in the above journal article is taken from the conference
%version [2].
%
%REFERENCES:
%[1] Z. Zhao, X. R. Li, and V. P. Jilkov, "Best linear unbiased filtering
%    with nonlinear measurements for target tracking," IEEE Transactions on
%    Aerospace and Electronic Systems, vol. 40, no. 4, pp. 1324-1336,
%    Oct. 2004.
%[2] Z. Zhao, X. R. Li, and V. P. Jilkov, "Optimal linear unbiased
%    filtering with polar measurements for target tracking," in Proceedings
%    of the 5th International Conference on Information Fusion, Annapolis,
%    MD, 8-11 Jul. 2002, pp. 1527-1534.
%
%February 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%Change the ordering of the elements in the state to match that used in the
%paper.
permuteIdx=[1;3;2;4];
undoPermuteIdx=[1;3;2;4];%The inverse permutation
xStateBar=xStatePred(permuteIdx);
PBar=PPred(permuteIdx,permuteIdx);

xBar=xStateBar(1);
yBar=xStateBar(3);

covXTilde=PBar(1,1);
covYTilde=PBar(3,3);
covXTileYTilde=PBar(1,3);

sigmaR2=R(1,1);
sigmaTheta2=R(2,2);

%The Taylor series approximation to the expected value of
%(y^2-x^2)/(x^2+y^2) from Section VI of the journal article.
denom=(xBar^2+yBar^2)^3;
ExyRat1=(yBar^2-xBar^2)/(xBar^2+yBar^2)+...
        2*yBar^2*(yBar^2-3*xBar^2)*covXTilde/denom+...
        4*xBar*yBar*(xBar^2-yBar^2)*covXTileYTilde/denom-...
        2*xBar^2*(xBar^2-3*yBar^2)*covYTilde/denom;

%The Taylor series approximation to the expected value of x*y/(x^2+y^2)
%from Section VI is taken from the form written out explicitly in Equation
%18 of the conference paper.
ExyRat2=xBar*yBar/(xBar^2+yBar^2)+(1/2)*(...
        2*xBar*yBar*(xBar^2-3*yBar^2)*covXTilde/denom+...
        (6*xBar^2*yBar^2-xBar^4-yBar^4)*covXTileYTilde/denom+...
        2*xBar*yBar*(yBar^2-3*xBar^2)*covYTilde/denom);

lambda1=exp(-sigmaTheta2/2);
lambda2=(1/2)*(1+exp(-2*sigmaTheta2));
lambda3=(1/2)*(1-exp(-2*sigmaTheta2));

%The value of S11 is from Section IV of the journal article.
Pzz(1,1)=lambda2*covXTilde+lambda3*covYTilde+(1/2)*sigmaR2+lambda3*yBar^2+...
    (lambda2-lambda1^2)*xBar^2+(1/2)*sigmaR2*exp(-2*sigmaTheta2)*(-ExyRat1);
Pzz(2,2)=lambda2*covYTilde+lambda3*covXTilde+(1/2)*sigmaR2+lambda3*xBar^2+...
    (lambda2-lambda1^2)*yBar^2+(1/2)*sigmaR2*exp(-2*sigmaTheta2)*(ExyRat1);
Pzz(1,2)=exp(-2*sigmaTheta2)*covXTileYTilde+(exp(-2*sigmaTheta2)-...
    lambda1^2)*xBar*yBar-sigmaR2*exp(-2*sigmaTheta2)*ExyRat2;
Pzz(2,1)=Pzz(1,2);

%Equation 17
covXTildeZTilde=lambda1*[PBar(:,1),PBar(:,3)];
W=covXTildeZTilde/Pzz;
zCartPred=lambda1*[xBar;yBar];

zCart=pol2Cart(z);
innov=zCart-zCartPred;
xStateUpdate=xStateBar+W*(zCart-zCartPred);
PUpdate=PBar-W*Pzz*W';

%Make the ordering go back to the way it was.
xStateUpdate=xStateUpdate(undoPermuteIdx);
PUpdate=PUpdate(undoPermuteIdx,undoPermuteIdx);

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function tests = test_VBA_random
% Unit Tests for VBA_random

tests = functiontests (localfunctions);

% -------------------------------------------------------------------------
function test_arbitrary (testCase)
    K = 3;
    p = rand (K, 1);
    p = p / sum (p);
    vals_1 = 'abc';
    vals_k = ['ab'; 'cd'; 'ef'];
    M = size(vals_k, 2);

    % should fail on invalid parameters
    shouldFail = @() VBA_random ('Arbitrary', [0.5; 0.5], vals_1);
    testCase.verifyError(shouldFail, 'VBA:invalidInput');
    shouldFail = @() VBA_random ('Arbitrary', zeros (3, 1), vals_1);
    testCase.verifyError(shouldFail, 'VBA:invalidInput');
    shouldFail = @() VBA_random ('Arbitrary', p, vals_k, 8, 8);
    testCase.verifyError(shouldFail, 'VBA:invalidInput');

    % should return one sample by default
    actual = VBA_random ('Arbitrary', p, vals_1);
    testCase.verifyNumElements (actual, 1);
    actual = VBA_random ('Arbitrary', p, vals_k);
    testCase.verifySize (actual, [M, 1]);

    % should return matrix on scalar N
    N = 8;
    actual = VBA_random ('Arbitrary', p, vals_1, N);
    testCase.verifySize (actual, [N, N]);
    actual = VBA_random ('Arbitrary', p, vals_k, N);
    testCase.verifySize (actual, [M, N]);

    % should return matrix of asked dimension
    N = num2cell (1 + randi (5, 1, 4));
    actual = VBA_random ('Arbitrary', p, vals_1, N{:});
    testCase.verifySize (actual, [N{:}]);

    % should return sample according to distribution
    % + univariate
    actual = VBA_random ('Arbitrary', p, vals_1, 1, 1e6);
    % - support
    testCase.verifyEqual (unique (actual), vals_1);
    % + density
    for i = 1 : numel (p)
        testCase.verifyEqual (mean (actual == vals_1(i)), p(i), 'AbsTol', 1e-2);
    end
    % + multivariate
    actual = VBA_random ('Arbitrary', p, vals_k, 1e6);
    % - support
    testCase.verifyEqual (unique (actual', 'rows'), vals_k);
    % + density
    for i = 1 : numel (p)
        isEq = all(bsxfun(@eq, actual, vals_k(i,:)'));
        testCase.verifyEqual (mean (isEq, 2), p(i), 'AbsTol', 1e-2);
    end

% -------------------------------------------------------------------------
function test_bernoulli (testCase)
    p = rand ();

    % should fail on invalid probability
    shouldFail = @() VBA_random ('Bernoulli', -3);
    testCase.verifyError(shouldFail, 'VBA:invalidInput');
    % should fail on invalid size
    shouldFail = @() VBA_random ('Bernoulli', rand (3), 2);
    testCase.verifyError(shouldFail, 'VBA:invalidInput');

    % should return one sample by default
    actual = VBA_random ('Bernoulli', p);
    testCase.verifyNumElements (actual, 1);

    % should return matrix on scalar N
    N = 3;
    actual = VBA_random ('Bernoulli', p, N);
    testCase.verifySize (actual, [N, N]);

    % should return matrix of asked dimension
    N = num2cell (1 + randi (5, 1, 4));
    actual = VBA_random ('Bernoulli', p,  N{:});
    testCase.verifySize (actual, [N{:}]);

    % should return matrix on matrix p
    actual = VBA_random ('Bernoulli', rand (N{:}));
    testCase.verifySize (actual, [N{:}]);

    % should return sample according to distribution
    actual = VBA_random ('Bernoulli', p, 1, 1e6);
    % + support
    testCase.verifyEqual (unique (actual), [0 1]);
    % + mean
    testCase.verifyEqual (mean (actual), p, 'AbsTol', 1e-2);

    % should deal with nan
    actual = VBA_random ('Bernoulli', [p NaN]);
    testCase.verifyTrue (~ isnan (actual(1)));
    testCase.verifyTrue (isnan (actual(2)));

% -------------------------------------------------------------------------
function test_binomial (testCase)
    n = 2;
    p = rand ();

    % should fail on invalid parameters
    shouldFail = @() VBA_random ('Binomial',0, p);
    testCase.verifyError(shouldFail, 'VBA:invalidInput');
    shouldFail = @() VBA_random ('Binomial',n, -3);
    testCase.verifyError(shouldFail, 'VBA:invalidInput');

    % should return one sample by default
    actual = VBA_random ('Binomial', n, p);
    testCase.verifyNumElements (actual, 1);

    % should return matrix on scalar N
    N = 3;
    actual = VBA_random ('Binomial', n, p, N);
    testCase.verifySize (actual, [N, N]);

    % should return matrix of asked dimension
    N = num2cell (1 + randi (5, 1, 4));
    actual = VBA_random ('Binomial', n, p,  N{:});
    testCase.verifySize (actual, [N{:}]);

    % should return sample according to distribution
    actual = VBA_random ('Binomial', n, p, 1, 1e6);
    % + support
    testCase.verifyEqual (unique (actual), 0 : n);

    % + mean
    testCase.verifyEqual (mean (actual), n * p, 'AbsTol', 1e-2);
    % + variance
    testCase.verifyEqual (var (actual), n * p * (1 - p), 'AbsTol', 1e-2);
    % + density
    for k = 0 : n
        expected = nchoosek (n, k) * (p ^ k) * ((1 - p) ^ (n - k));
        testCase.verifyEqual (mean (actual == k), expected, 'AbsTol', 1e-2);
    end

% -------------------------------------------------------------------------
function test_categorical (testCase)
    K = 3;
    p = rand (K, 1);
    p = p / sum(p);

    % should fail on invalid parameters
    shouldFail = @() VBA_random ('Categorical', zeros (K, 1));
    testCase.verifyError(shouldFail, 'VBA:invalidInput');

    % should return one sample by default
    actual = VBA_random ('Categorical', p);
    testCase.verifyNumElements (actual, 1);

    % should return matrix on scalar N
    N = 3;
    actual = VBA_random ('Categorical', p, N);
    testCase.verifySize (actual, [N, N]);

    % should return matrix of asked dimension
    N = num2cell (1 + randi (5, 1, 4));
    actual = VBA_random ('Categorical', p,  N{:});
    testCase.verifySize (actual, [N{:}]);

    % should return sample according to distribution
    actual = VBA_random ('Categorical', p, 1, 1e6);
    % + support
    testCase.verifyEqual (unique (actual), 1 : K);
    % + mean
    testCase.verifyEqual (mean (actual), (1 : K) * p, 'AbsTol', 1e-2);
   % + density
    for k = 1 : K
        testCase.verifyEqual (mean (actual == k), p(k), 'AbsTol', 1e-2);
    end

    % -------------------------------------------------------------------------
function test_dirichlet (testCase)
    K = 3;
    alpha = 1 + randi (10, K, 1);

    % should fail on invalid parameters
    shouldFail = @() VBA_random ('Dirichlet', rand);
    testCase.verifyError(shouldFail, 'VBA:invalidInput');
    shouldFail = @() VBA_random ('Dirichlet', alpha, 2, 3);
    testCase.verifyError(shouldFail, 'VBA:invalidInput');

    % should return one sample by default
    actual = VBA_random ('Dirichlet', alpha);
    testCase.verifySize (actual, [K, 1]);

    % should return N samples
    N = 8;
    actual = VBA_random ('Dirichlet', alpha, N);
    testCase.verifySize (actual, [K, N]);

    % should return sample according to distribution
    N = 1e6;
    actual = VBA_random ('Dirichlet', alpha, N);
    % + support
    testCase.verifyTrue (VBA_isInRange (actual, [0 1]));
    testCase.verifyEqual (sum (actual), ones (1, N), 'AbsTol', 1e-13);
    % + moments
    [expected.m, expected.v] = VBA_dirichlet_moments(alpha);
    testCase.verifyEqual (mean (actual, 2), expected.m , 'AbsTol', 1e-2);
    testCase.verifyEqual (var (actual, [], 2), expected.v , 'AbsTol', 1e-2);

% -------------------------------------------------------------------------
function test_gamma (testCase)
    a = 10 * rand ();
    b = 10 * rand ();

    % should fail on invalid parameters
    shouldFail = @() VBA_random ('Gamma', - 1, b);
    testCase.verifyError(shouldFail, 'VBA:invalidInput');
    shouldFail = @() VBA_random ('Gamma', a, - 1);
    testCase.verifyError(shouldFail, 'VBA:invalidInput');

    % should return one sample by default
    actual = VBA_random ('Gamma', a, b);
    testCase.verifyNumElements (actual, 1);

    % should return matrix on scalar N
    N = 3;
    actual = VBA_random ('Gamma', a, b, N);
    testCase.verifySize (actual, [N, N]);

    % should return matrix of asked dimension
    N = num2cell (1 + randi (5, 1, 4));
    actual = VBA_random ('Gamma', a, b,  N{:});
    testCase.verifySize (actual, [N{:}]);

    % should return sample according to distribution
    actual = VBA_random ('Gamma', a, b, 1, 1e6);
    % + support
    testCase.verifyTrue (all (actual > 0));
    % + mean
    testCase.verifyEqual (mean (actual), a * b, 'RelTol', 1e-2);
    % + variance
    testCase.verifyEqual (var (actual), a * b ^ 2, 'RelTol', 1e-2);

  % -------------------------------------------------------------------------
function test_gaussian (testCase)

    % + univariate
    mu_1 = randn ();
    Sigma_1 = rand () ^ 2;
    % + multivariate
    k = 2;
    mu_k = randn (k, 1);
    Sigma_k = rand (k);
    Sigma_k = Sigma_k * Sigma_k';

     % should fail on invalid parameters
    shouldFail = @() VBA_random ('Gaussian', mu_1, Sigma_k);
    testCase.verifyError(shouldFail, 'VBA:invalidInput');
    shouldFail = @() VBA_random ('Gaussian', mu_k, Sigma_k, 2, 2);
    testCase.verifyError(shouldFail, 'VBA:invalidInput');

    % should return one sample by default
    % + univariate
    actual = VBA_random ('Gaussian', mu_1, Sigma_1);
    testCase.verifyNumElements (actual, 1);
    % + multivariate
    actual = VBA_random ('Gaussian', mu_k, Sigma_k);
    testCase.verifySize (actual, [k, 1]);

    % should return matrix on scalar N
    N = 5;
    % + univariate
    actual = VBA_random ('Gaussian', mu_1, Sigma_1, N);
    testCase.verifySize (actual, [N, N]);
    % + multivariate
    actual = VBA_random ('Gaussian', mu_k, Sigma_k, N);
    testCase.verifySize (actual, [k, N]);

    % should return matrix of asked dimension
    N = num2cell (1 + randi (5, 1, 4));
    actual = VBA_random ('Gaussian', mu_1, Sigma_1,  N{:});
    testCase.verifySize (actual, [N{:}]);

    % should return sample according to distribution
    % + univariate
    actual = VBA_random ('Gaussian', mu_1, Sigma_1, 1, 1e6);
    % - mean
    testCase.verifyEqual (mean (actual), mu_1, 'AbsTol', 1e-2);
    % - variance
    testCase.verifyEqual (var (actual), Sigma_1, 'AbsTol', 1e-2);
    % + multivariate
    actual = VBA_random ('Gaussian', mu_k, Sigma_k, 1e6);
    % - mean
    testCase.verifyEqual (mean (actual, 2), mu_k, 'AbsTol', 1e-2);
    % - variance
    testCase.verifyEqual (cov (actual'), Sigma_k, 'AbsTol', 1e-2);

% -------------------------------------------------------------------------
function test_multinomial (testCase)
    n = 2;
    K = 3;
    p = rand (K, 1);
    p = p / sum(p);

    % should fail on invalid parameters
    shouldFail = @() VBA_random ('Multinomial',0, p);
    testCase.verifyError(shouldFail, 'VBA:invalidInput');
    shouldFail = @() VBA_random ('Multinomial',n, .5);
    testCase.verifyError(shouldFail, 'VBA:invalidInput');
    shouldFail = @() VBA_random ('Multinomial',n, [1; 1]);
    testCase.verifyError(shouldFail, 'VBA:invalidInput');
    shouldFail = @() VBA_random ('Multinomial',n, p, 8, 8);
    testCase.verifyError(shouldFail, 'VBA:invalidInput');

    % should return one sample by default
    actual = VBA_random ('Multinomial', n, p);
    testCase.verifySize (actual, [K, 1]);

    % should return matrix on scalar N
    N = 8;
    actual = VBA_random ('Multinomial', n, p, N);
    testCase.verifySize (actual, [K, N]);

    % should return sample according to distribution
    N = 1e6;
    actual = VBA_random ('Multinomial', n, p, N);
    % + support
    testCase.verifyEqual (unique (actual)', 0 : n);
    testCase.verifyEqual (sum (actual), n * ones (1, N), 'AbsTol', 1e-13);
    % + density
    for i = 1 : K
        expected = n * p(i);
        testCase.verifyEqual (mean (actual(i, :)), expected, 'AbsTol', 1e-2);
        expected = n * p(i) * (1 - p(i));
        testCase.verifyEqual (var (actual(i, :)), expected, 'AbsTol', 1e-2);
    end

    % should deal with nan
    actual = VBA_random ('Multinomial', n, nan(K, 1));
    testCase.verifySize (actual, [K, 1]);
    testCase.verifyTrue (all (isnan (actual)));

author: MBB-team

repo: VBA-toolbox


text:
function [J_x0,J_xf,J_t] = sys_Dg(neq,t,x0,xf)

global sys_params A B C D

% J_x0 and J_xf are row vectors of length n.
% J_t is not used.

F_NUM = neq(5);

J_x0 = zeros(1, length(x0));
J_xf = 0*C;


text:
function [ xy, line_pointer, line_data ] = xyl_example ( point_num, ...
  line_num, line_data_num )

%*****************************************************************************80
%
%% XYL_EXAMPLE sets data suitable for a pair of XY and XYL files.
%
%  Discussion:
%
%    There are 13 points.
%    There are 3 lines.
%    There are 15 line data items.
%
%         4 12-11
%         /\ | |
%        /  \| |
%       /   13 |
%      /      \10
%     /        \
%    5          3
%    |          |
%    |     9--8 |
%    |     |  | |
%    |     |  | |
%    |     6--7 |
%    |          |
%    1----------2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, integer LINE_NUM, the number of lines.
%
%    Input, integer LINE_DATA_NUM, the number of line items.
%
%    Output, real XY(2,POINT_NUM), the point coordinates.
%
%    Output, integer LINE_POINTER(LINE_NUM+1), pointers to the
%    first line item for each line.
%
%    Output, integer LINE_DATA(LINE_DATA_NUM), indices
%    of points that form lines.
%
  xy(1:2,1:point_num) = [ ...
     0.0,   0.0; ...
     6.0,   0.0; ...
     6.0,   7.0; ...
     3.0,  10.0; ...
     0.0,   7.0; ...
     4.0,   1.0; ...
     5.0,   1.0; ...
     5.0,   4.0; ...
     4.0,   4.0; ...
     5.0,   8.0; ...
     5.0,  11.0; ...
     4.0,  11.0; ...
     4.0,   9.0 ]';

  line_pointer(1:line_num+1) = [ 1, 7, 12, 16 ];

  line_data(1:line_data_num) = [ ...
     1,  2,  3,  4,  5,  1, ...
     6,  7,  8,  9,  6, ...
    10, 11, 12, 13 ];

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function zi = interp(DEM,xi,yi,method)

%INTERP interpolate to query locations
%
% Syntax
%
%     zi = interp(DEM,xi,yi)
%     zi = interp(DEM,xi,yi,method)
%
% Description
%
%     INTERP uses the griddedInterpolant class to interpolate values in the
%     instance of GRIDobj (DEM) to query locations at xi and yi.
%     If DEM.Z is an integer class, interp will convert it to single
%     precision and use linear interpolation as default. If DEM.Z is
%     logical, nearest neigbhor will be used by default.
%
% Input arguments
%
%     DEM     instance of GRIDobj
%     xi,yi   x- and y-coordinates of query locations
%     method  interpolation method (default = 'linear'). See the
%             documentation of the griddedInterpolant class for further
%             methods
%
% Output arguments
%
%     zi      interpolated values at query locations
%
% Example
%
%     DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
%     [x,y] = getoutline(DEM);
%     xy = rand(20,2);
%     xy(:,1) = xy(:,1)*(max(x)-min(x)) + min(x);
%     xy(:,2) = xy(:,2)*(max(y)-min(y)) + min(y);
%     z = interp(DEM,xy(:,1),xy(:,2));
%
%
% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 17. August, 2017

narginchk(3,4)

if nargin == 3;
    method = 'linear';
else
    method = validatestring(method,...
        {'linear','nearest','spline','pchip','cubic'},'interp','method',4);
end

% created griddedInterpolant class
[x,y] = getcoordinates(DEM);
% flip y to get monotonic increasing grid vectors
y     = y(end:-1:1);

if isinteger(DEM.Z)
    cla   = class(DEM.Z);
    DEM.Z = single(DEM.Z);
    convoutput = true;
elseif islogical(DEM.Z)
    cla   = class(DEM.Z);
    DEM.Z = single(DEM.Z);
    if nargin == 3;
        method = 'nearest';
    end
    convoutput = true;
else
    convoutput = false;
end

F     = griddedInterpolant({x,y},flipud(DEM.Z)',method,'none');

% interpolate
zi     = F(xi,yi);

if convoutput
    zi = cast(zi,cla);
end

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
function y = daub10_transform ( n, x )

%*****************************************************************************80
%
%% DAUB10_TRANSFORM computes the DAUB10 transform of a vector.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of the vector.
%    N must be a power of 2 and at least 4.
%
%    Input, real X(N), the vector to be transformed.
%
%    Output, real Y(N), the transformed vector.
%
  c = [ ...
    0.1601023979741929; ...
    0.6038292697971895; ...
    0.7243085284377726; ...
    0.1384281459013203; ...
   -0.2422948870663823; ...
   -0.0322448695846381; ...
    0.0775714938400459; ...
   -0.0062414902127983; ...
   -0.0125807519990820; ...
    0.0033357252854738 ];
  p = 9;
  y(1:n,1) = x(1:n);
  m = n;
  q = floor ( ( p - 1 ) / 2 );

  while ( 4 <= m )

    i = 1;
    z(1:m,1) = 0.0;

    for j = 1 : 2 : m - 1

      mh = floor ( m / 2 );
      for k = 0 : 2 : p - 1
        j0 = i4_wrap ( j + k,     1, m );
        j1 = i4_wrap ( j + k + 1, 1, m );
        z(i,1)    = z(i,1)    + c(  k+1) * y(j0) + c(  k+2) * y(j1);
        z(i+mh,1) = z(i+mh,1) + c(p-k+1) * y(j0) - c(p-k  ) * y(j1);
      end

      i = i + 1;

    end

    y(1:m,1) = z(1:m);

    m = floor ( m / 2 );

  end

  return
end


text:
function g = sigmoid(z)
%SIGMOID Compute sigmoid functoon
%   J = SIGMOID(z) computes the sigmoid of z.

g = 1 ./ (1 + exp(-z));

end

author: atinesh-s

repo: Coursera-Machine-Learning-Stanford


text:
% DEMSILHOUETTEGP2 Model silhouette data with independent MLP GPs.

% FORMAT
% DESC runs a simple regression on the Agawal and Triggs data.
%
% SEEALSO : gpCreate, demInterpolation
%
% COPYRIGHT : Neil D. Lawrence, 2008

% GP

randn('seed', 1e7)
rand('seed', 1e7)

dataSetName = 'silhouette';
experimentNo = 2;

% load data
[X, y, XTest, yTest] = mapLoadData(dataSetName);

% Set up the model
options = gpOptions('ftc');
options.kern{1} = 'mlp';

% Scale outputs to variance 1.
options.scale2var1 = true;

% Use the full Gaussian process model.
q = size(X, 2);
d = size(y, 2);
model = gpCreate(q, d, X, y, options);

display = 1;
iters = 1000;

model = gpOptimise(model, display, iters);
modelDisplay(model)

% Save results
fileBaseName = modelWriteResult(dataSetName, experimentNo);
demSilhouettePlot

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function h=semilogy(varargin)
%PREAL/PLOT  Overloaded semilogy function for class PREAL.

global useUnitsFlag

for k=1:length(varargin)
    if isa(varargin{k},'preal'), varargin{k}=double(varargin{k}); end
end
if nargout==0
    semilogy(varargin{:});
else
    h=semilogy(varargin{:});
end

author: Sable

repo: mcbench-benchmarks


text:
function [E, Es, O] = detectEdge(Iin, Din, id, C, model, sc, outFile, cacheFile)
  if(isstr(Iin)), Iin = imread(Iin); end
  if(isstr(Din)), Din = imread(Din); end
  if(isstr(model)), load(model); end

  opts = model.opts;
  rgbd = opts.rgbd;

  model.opts.multiscale=1;          % for top accuracy set multiscale=1
  model.opts.nTreesEval=4;          % for top speed set nTreesEval=1
  model.opts.nThreads=1;            % max number threads for evaluation
  model.opts.sharpen=0;             % for top speed set sharpen=0

  if(rgbd == 3),
    colorModel = model.colorModel.model;
    cues = getAllCues(Iin, Din, C, colorModel, opts.rgbd3opts.vars, cacheFile);
  end

  for i = 1:length(sc),
    I = imresize(Iin, sc(i), 'lanczos3');
    ng = zeros(size(I,1), size(I,2), 0);

    if(rgbd),
      D = imresize(Din, sc(i), 'nearest');
      D=single(D)/1e4;
    end
    if(rgbd==1),
      I=D;
    elseif(rgbd==2),
      I=cat(3,single(I)/255,D);
    elseif(rgbd==3),
      I=cat(3,single(I)/255,D,1e-3./D);
      ng = cat(3, cues{:});
      ng = imresize(ng, sc(i), 'lanczos3');
    end
    model.opts.nms = 0;
    [Es{i}, O{i}] = edgesDetect(I,model,ng);
    E{i} = edgesNmsMex(Es{i}, O{i}, 1, 5, 1.01, model.opts.nThreads);

    % model.opts.nms = 1;
    % [E{i} O{i}] = edgesDetect(I,model,ng);
  end
  if(~isempty(outFile)), save(outFile, 'E', 'Es', 'O'); end
end

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth