text:
function organizeExperiments_Radiomics(pathFeatures,pathLR,pathCR,outcomes,timeToEvent,nPatient,nonTextName,fSetNames,paramSEP,baselineSEP,freedomSEP,textType,textName)

startpath = pwd;

% INITIALIZATION
roi = 'GTVtot'; % Using "GTVtot" ROI for all experiments
permTrain = randperm(nPatient.HGJ + nPatient.CHUS); permTest = randperm(nPatient.HMR + nPatient.CHUM); % Random permutation of instances.
nNonText = numel(nonTextName);
nFset = numel(fSetNames);
nameOutcomes = fieldnames(outcomes.HGJ); nOutcomes = numel(nameOutcomes);
nTextType = numel(textType);

% DATA ORGANIZATION
training = struct; testing = struct;
for o = 1:nOutcomes

    % Organizing non-texture features
    cd(pathFeatures)
    train1 = load(['nonText_HGJ_',roi]); train1 = struct2cell(train1); train1 = train1{1};
    train2 = load(['nonText_CHUS_',roi]); train2 = struct2cell(train2); train2 = train2{1};
    test1 = load(['nonText_HMR_',roi]); test1 = struct2cell(test1); test1 = test1{1};
    test2 = load(['nonText_CHUM_',roi]); test2 = struct2cell(test2); test2 = test2{1};
    for i = 1:nNonText
        training.(nameOutcomes{o}).nonText.(nonTextName{i}).Data = [train1.(nonTextName{i}).Data;train2.(nonTextName{i}).Data];
        testing.(nameOutcomes{o}).nonText.(nonTextName{i}).Data = [test1.(nonTextName{i}).Data;test2.(nonTextName{i}).Data];
        training.(nameOutcomes{o}).nonText.(nonTextName{i}).Data = training.(nameOutcomes{o}).nonText.(nonTextName{i}).Data(permTrain);
        testing.(nameOutcomes{o}).nonText.(nonTextName{i}).Data = testing.(nameOutcomes{o}).nonText.(nonTextName{i}).Data(permTest);
    end

    % Organizing texture features
    for f = 1:nFset
        if strcmp(fSetNames{f},'PET') % Using only PET
            train1 = load(['text_HGJ_PT_',roi]); train1 = struct2cell(train1); train1 = train1{1};
            train2 = load(['text_CHUS_PT_',roi]); train2 = struct2cell(train2); train2 = train2{1};
            test1 = load(['text_HMR_PT_',roi]); test1 = struct2cell(test1); test1 = test1{1};
            test2 = load(['text_CHUM_PT_',roi]); test2 = struct2cell(test2); test2 = test2{1};
            train1 = {train1}; train2 = {train2};
            test1 = {test1}; test2 = {test2};
            training.(nameOutcomes{o}).text.(fSetNames{f}).param = paramSEP;
            training.(nameOutcomes{o}).text.(fSetNames{f}).baseline = baselineSEP;
            training.(nameOutcomes{o}).text.(fSetNames{f}).freedom = freedomSEP;
            training.(nameOutcomes{o}).text.(fSetNames{f}).paramName = {'Scale','Quant.algo','Ng'};
            textCellName = {'PT'};
        elseif strcmp(fSetNames{f},'CT') % Using only CT
            train1 = load(['text_HGJ_CT_',roi]); train1 = struct2cell(train1); train1 = train1{1};
            train2 = load(['text_CHUS_CT_',roi]); train2 = struct2cell(train2); train2 = train2{1};
            test1 = load(['text_HMR_CT_',roi]); test1 = struct2cell(test1); test1 = test1{1};
            test2 = load(['text_CHUM_CT_',roi]); test2 = struct2cell(test2); test2 = test2{1};
            train1 = {train1}; train2 = {train2};
            test1 = {test1}; test2 = {test2};
            training.(nameOutcomes{o}).text.(fSetNames{f}).param = paramSEP;
            training.(nameOutcomes{o}).text.(fSetNames{f}).baseline = baselineSEP;
            training.(nameOutcomes{o}).text.(fSetNames{f}).freedom = freedomSEP;
            training.(nameOutcomes{o}).text.(fSetNames{f}).paramName = {'Scale','Quant.algo','Ng'};
            textCellName = {'CT'};
        elseif strcmp(fSetNames{f},'PETCT') % Using both PET and CT
            train1_1 = load(['text_HGJ_PT_',roi]); train1_1 = struct2cell(train1_1); train1_1 = train1_1{1};
            train1_2 = load(['text_HGJ_CT_',roi]); train1_2 = struct2cell(train1_2); train1_2 = train1_2{1};
            train2_1 = load(['text_CHUS_PT_',roi]); train2_1 = struct2cell(train2_1); train2_1 = train2_1{1};
            train2_2 = load(['text_CHUS_CT_',roi]); train2_2 = struct2cell(train2_2); train2_2 = train2_2{1};
            test1_1 = load(['text_HMR_PT_',roi]); test1_1 = struct2cell(test1_1); test1_1 = test1_1{1};
            test1_2 = load(['text_HMR_CT_',roi]); test1_2 = struct2cell(test1_2); test1_2 = test1_2{1};
            test2_1 = load(['text_CHUM_PT_',roi]); test2_1 = struct2cell(test2_1); test2_1 = test2_1{1};
            test2_2 = load(['text_CHUM_CT_',roi]); test2_2 = struct2cell(test2_2); test2_2 = test2_2{1};
            train1 = {train1_1,train1_2}; train2 = {train2_1,train2_2};
            test1 = {test1_1,test1_2}; test2 = {test2_1,test2_2};
            training.(nameOutcomes{o}).text.(fSetNames{f}).param = paramSEP;
            training.(nameOutcomes{o}).text.(fSetNames{f}).baseline = baselineSEP;
            training.(nameOutcomes{o}).text.(fSetNames{f}).freedom = freedomSEP;
            training.(nameOutcomes{o}).text.(fSetNames{f}).paramName = {'Scale','Quant.algo','Ng'};
            textCellName = {'PT','CT'};
        end

        for i = 1:numel(train1)
            sizeParam = size(train1{1});
            training.(nameOutcomes{o}).text.(fSetNames{f}).(textCellName{i}) = cell(sizeParam);
            testing.(nameOutcomes{o}).text.(fSetNames{f}).(textCellName{i}) = cell(sizeParam);
            for n = 1:numel(train1{1})
                for type = 1:nTextType
                    for text = 1:numel(textName{type})
                        training.(nameOutcomes{o}).text.(fSetNames{f}).(textCellName{i}){n}.(textType{type}).(textName{type}{text}).Data = [train1{i}{n}.(textType{type}).(textName{type}{text}).Data;train2{i}{n}.(textType{type}).(textName{type}{text}).Data];
                        testing.(nameOutcomes{o}).text.(fSetNames{f}).(textCellName{i}){n}.(textType{type}).(textName{type}{text}).Data = [test1{i}{n}.(textType{type}).(textName{type}{text}).Data;test2{i}{n}.(textType{type}).(textName{type}{text}).Data];
                        training.(nameOutcomes{o}).text.(fSetNames{f}).(textCellName{i}){n}.(textType{type}).(textName{type}{text}).Data = training.(nameOutcomes{o}).text.(fSetNames{f}).(textCellName{i}){n}.(textType{type}).(textName{type}{text}).Data(permTrain);
                        testing.(nameOutcomes{o}).text.(fSetNames{f}).(textCellName{i}){n}.(textType{type}).(textName{type}{text}).Data = testing.(nameOutcomes{o}).text.(fSetNames{f}).(textCellName{i}){n}.(textType{type}).(textName{type}{text}).Data(permTest);
                    end
                end
            end
        end
    end
end

% Organizing outcomes
for o = 1:nOutcomes
    training.(nameOutcomes{o}).outcome = [outcomes.HGJ.(nameOutcomes{o});outcomes.CHUS.(nameOutcomes{o})];
    testing.(nameOutcomes{o}).outcome = [outcomes.HMR.(nameOutcomes{o});outcomes.CHUM.(nameOutcomes{o})];
    training.(nameOutcomes{o}).timeToEvent = [timeToEvent.HGJ.(nameOutcomes{o});timeToEvent.CHUS.(nameOutcomes{o})];
    testing.(nameOutcomes{o}).timeToEvent = [timeToEvent.HMR.(nameOutcomes{o});timeToEvent.CHUM.(nameOutcomes{o})];
    training.(nameOutcomes{o}).outcome = training.(nameOutcomes{o}).outcome(permTrain);
    testing.(nameOutcomes{o}).outcome = testing.(nameOutcomes{o}).outcome(permTest);
    training.(nameOutcomes{o}).timeToEvent = training.(nameOutcomes{o}).timeToEvent(permTrain);
    testing.(nameOutcomes{o}).timeToEvent = testing.(nameOutcomes{o}).timeToEvent(permTest);
end

cd(pathLR), save('training','training'), save('testing','testing'), save('permTrain','permTrain'), save('permTest','permTest')
cd(pathCR), save('training','training'), save('testing','testing'), save('permTrain','permTrain'), save('permTest','permTest')

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function varargout = edgesEval( model, varargin )
% Run and evaluate structured edge detector on BSDS500.
%
% This function first runs the trained structured edge detector on every
% test or validation image in BSDS then call edgesEvalDir.m to perform the
% actual edge evaluation. edgesEval is specific to the structured edge
% detector and BSDS, edgesEvalDir is general purpose edge evaluation code.
% For example usage of edgesEval see edgesDemo.
%
% USAGE
%  varargout = edgesEval( model, prms )
%
% INPUTS
%  model      - structured edge model trained with edgesTrain
%  prms       - parameters (struct or name/value pairs)
%   .dataType   - ['test'] should be either 'test' or 'val'
%   .name       - [''] name to append to evaluation
%   .opts       - {} list of model opts to overwrite
%   .show       - [0] if true plot results using edgesEvalPlot
%   .pDistr     - [{'type','parfor'}] parameters for fevalDistr
%   .cleanup    - [0] if true delete temporary files
%   .thrs       - [99] number or vector of thresholds for evaluation
%   .maxDist    - [.0075] maximum tolerance for edge match
%
% OUTPUTS
%  varargout  - same outputs as edgesEvalDir
%
% EXAMPLE
%
% See also edgesDemo, edgesDetect, edgesEvalDir, edgesEvalPlot
%
% Structured Edge Detection Toolbox      Version 3.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-microsoft.com]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the MSR-LA Full Rights License [see license.txt]

% get default parameters
dfs={'dataType','test', 'name','', 'opts',{}, 'show',0, ...
  'pDistr',{{'type','parfor'}}, 'cleanup',0, 'thrs',99, 'maxDist',.0075 };
p=getPrmDflt(varargin,dfs,1);

% load model and update model.opts acoording to opts
if( ischar(model) ), model=load(model); model=model.model; end
for i=1:length(p.opts)/2, model.opts.(p.opts{i*2-1})=p.opts{i*2}; end
rgbd=model.opts.rgbd; model.opts.nms=1;

% get list of relevant directories (image, depth, gt, results)
name = [model.opts.modelFnm,p.name];
imgDir = fullfile(model.opts.bsdsDir,'images',p.dataType);
depDir = fullfile(model.opts.bsdsDir,'depth',p.dataType);
gtDir  = fullfile(model.opts.bsdsDir,'groundTruth',p.dataType);
resDir = fullfile(model.opts.modelDir,p.dataType,name);
assert(exist(imgDir,'dir')==7); assert(exist(gtDir,'dir')==7);

% run edgesDetect() on every image in imgDir and store in resDir
if( ~exist(fullfile([resDir '-eval'],'eval_bdry.txt'),'file') )
  if(~exist(resDir,'dir')), mkdir(resDir); end
  ids=dir(imgDir); ids=ids([ids.bytes]>0); ids={ids.name}; n=length(ids);
  ext=ids{1}(end-2:end); for i=1:n, ids{i}=ids{i}(1:end-4); end
  res=cell(1,n); for i=1:n, res{i}=fullfile(resDir,[ids{i} '.png']); end
  do=false(1,n); for i=1:n, do(i)=~exist(res{i},'file'); end
  ids=ids(do); res=res(do); m=length(ids);
  parfor i=1:m, id=ids{i};
    I = imread(fullfile(imgDir,[id '.' ext])); D=[];
    if(rgbd), D=single(imread(fullfile(depDir,[id '.png'])))/1e4; end
    if(rgbd==1), I=D; elseif(rgbd==2), I=cat(3,single(I)/255,D); end
    E=edgesDetect(I,model); imwrite(uint8(E*255),res{i});
  end
end

% perform actual evaluation using edgesEvalDir
varargout=cell(1,max(1,nargout));
[varargout{:}] = edgesEvalDir('resDir',resDir,'gtDir',gtDir,...
  'pDistr',p.pDistr,'cleanup',p.cleanup,'thrs',p.thrs,'maxDist',p.maxDist);
if( p.show ), figure(p.show); edgesEvalPlot(resDir,name); end

end

author: Cloud-CV

repo: object-proposals


text:
function [Cext,val] = extendLSAPEinstance(C,val)
% Extend a cost matrix for error-correcting matching to its square version
%
%  [Cext,val] = extendCostMtx(C,val)
%
%  inputs:
%  C (n+1)x(m+1) cost matrix of an error-correcting bipartite graph
%  val the value used to fill 2nd and 3rd block of Cext (optional)
%  val = max(max(C)) + 1 if not given
%
%  outputs:
%  Cext (n+m)x(m+n) extended cost matrix
%  valcpt (optional) is the value val used to fill 2nd and 3rd block of
%  Cext
%
%         | c(1,1)   ...  c(1,m-1)  | c(1,m) val   ...      val   |
%         |   .              .      |  val  c(2,m) val ...   .    |
%         |   .              .      |   .                   val   |
%         | c(n-1,1) ... c(n-1,m-1) |  val     ...   val c(n-1,m) |
%  Cext =  -------------------------------------------------------
%         | c(n,1) val  ...    val  |                             |
%         |  val c(n,2) val ..  .   |           0_{m,n}           |
%         |   .                val  |                             |
%         |  val  ...  val c(n,m-1) |                             |
%
%
%   author: Sebastien Bougleux
%   institution: Normandie Univ, UNICAEN, ENSICAEN, CNRS, GREYC
% -----------------------------------------------------------
% This file is part of LSAPE.
% LSAPE is free software: you can redistribute it and/or modify
% it under the terms of the CeCILL-C License. See README file
% for more details.
% -----------------------------------------------------------

    if ~exist('val','var') || nargin == 1
        val = max(max(C)) + 1;
    end

    D = val .* cast(ones(size(C,1)-1),class(C));
    E = val .* cast(ones(size(C,2)-1),class(C));
    D(logical(eye(size(D)))) = C(1:end-1,end);
    E(logical(eye(size(E)))) = C(end,1:end-1);
    Cext = [[C(1:end-1,1:end-1),D];[E,zeros(size(C,2)-1,size(C,1)-1)]];

end

author: optimaltransport


text:
function circle_arc_grid_test ( )

%*****************************************************************************80
%
%% CIRCLE_ARC_GRID_TEST tests the CIRCLE_ARC_GRID library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 November 2011
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CIRCLE_ARC_GRID_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test CIRCLE_ARC_GRID.\n' );

  circle_arc_grid_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CIRCLE_ARC_GRID_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function polygon_grid_display ( n, nv, v, ng, xg )

%*****************************************************************************80
%
%% POLYGON_GRID_DISPLAY displays grid points inside a polygon.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 April 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of subintervals.
%
%    Input, integer NV, the number of vertices in the polygon.
%
%    Input, real V[NV,2], the coordinates of the vertices.
%
%    Input, integer NG, the number of grid points.
%
%    Input, real XG(NG,2), the grid points.
%
  clf

  hold on

  plot ( [v(nv,1),v(1,1)], [v(nv,2),v(1,2)], 'r-', 'LineWidth', 2 );
  for i = 1 : nv - 1
    plot ( v(i:i+1,1), v(i:i+1,2), 'r-', 'LineWidth', 2 );
  end
  plot ( v(:,1), v(:,2), 'r.', 'MarkerSize', 25 );

  vc(1,1) = sum ( v(1:nv,1) ) / nv;
  vc(1,2) = sum ( v(1:nv,2) ) / nv;
  for i = 1 : nv
    plot ( [v(i,1),vc(1,1)], [v(i,2),vc(1,2)], 'b-', 'LineWidth', 2 );
  end
  plot ( xg(:,1), xg(:,2), 'b.', 'MarkerSize', 20 );
  axis equal
  title ( sprintf ( 'Polygonal grid, N = %d, NG = %d', n, ng ) )
  grid on
  hold off

  return
end


text:
function result = normalizekm(km)
% Copyright 2012 Nino Shervashidze, Karsten Borgwardt
% normalizes kernelmatrix km
% such that diag(result) = 1, i.e. K(x,y) / sqrt(K(x,x) * K(y,y))
% @author Karsten Borgwardt
% @date June 3rd 2005
% all rights reserved

nv = sqrt(diag(km));
nm =  nv * nv';
knm = nm .^ -1;
for i = 1:size(knm,1)
for j = 1:size(knm,2)
if (knm(i,j) == Inf)
  knm(i,j) = 0;
end
end
end
result = km .* knm;

author: muhanzhang

repo: DGCNN

save_path: github-repos/MATLAB/muhanzhang-DGCNN


text:
function test_suite = testFunctionHandlesC
%testFunctionHandlesC Test file used by TestFunctionHandlesTest
%   Contains two passing tests that use a test fixture containing an intentional
%   error.

%   Steven L. Eddins
%   Copyright 2008 The MathWorks, Inc.

initTestSuite;

function testData = setUpFcn
testData = 5;

function testA(testData)
assertEqual(testData, 5);

function testB(testData)
assertEqual(testData, 5);

function tearDownFcn(testData)
% This assertion is expected to error.
assertEqual(testData, 20);

author: ricedsp


text:
classdef LaunchVehiclePowerSinkEnum < matlab.mixin.SetGet
    %LaunchVehiclePowerSinkEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        Simple('Simple Power Sink');
    end

    properties
        name(1,:) char
    end

    methods
        function obj = LaunchVehiclePowerSinkEnum(name)
            obj.name = name;
        end
    end

    methods(Static)
        function listBoxStr = getListBoxStr()
            m = enumeration('LaunchVehiclePowerSinkEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('LaunchVehiclePowerSinkEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('LaunchVehiclePowerSinkEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function GBPD = calcGBPD(gB,ebsd,varargin)
% compute the grain boundary plane distribution
%
% Syntax
%
%   GBPD = calcGBPD(gB,ebsd)
%
%   % use a specific halfwidth
%   GBPD = calcGBPD(gB,ebsd,'halfwidth',10*degree)
%
% Input
%  gB   - @grainBoundary
%  ebsd - @EBSD
%
% Output
%
%  GBPD - @S2FunHarmonic
%
% See also
%

%% step 1: extract data

% grain boundary directions
d = gB.direction;

% rotations that rotate the x vector towards the trace normals
omega = 90*degree + angle(d,vector3d.X,vector3d.Z);
rot = rotation.byAxisAngle(zvector,omega);

% the orientations that align the crystallographic x-axis with the trace
% normals
ori = [rot .* ebsd('id',gB.ebsdId(:,1)).orientations;...
  rot .* ebsd('id',gB.ebsdId(:,2)).orientations];

%% step 2: define kernel function

% define a kernel function that is a fibre through the crystallograhic
% z-axis and the crystallographic x-axis
psi = S2FunHarmonic(S2DeLaValleePoussinKernel('halfwidth',5*degree,varargin{:}));

psi = psi.radon;

bw = min(getMTEXpref('maxS2Bandwidth'),psi.bandwidth);
rot = rotation.byAxisAngle(xvector,90*degree);

% multiply this kernel function with the sin of the polar angle
fun = @(v) pi/2*psi.eval(rot*v) .* sin(angle(v,zvector));

% the final kernel function as S2Harmonic
psi = S2FunHarmonicSym.quadrature(fun, 'bandwidth', bw, ori.CS);

%% testing only

%n = 100000;
%cs = crystalSymmetry('3','x||b');

%d = vector3d.byPolar(90*degree,-30*degree)
%omega = 90*degree+angle(d,vector3d.X,vector3d.Z);
%ori = orientation.rand(n,cs);
%omega = rand(n,1);

%rot = rotation.byAxisAngle(zvector,omega);

%ori = rot .* orientation.id(cs);

%% step 3: compute orientation density

% compute the orientation density of the modified boundary orientations
odf = calcDensity(ori,'kernel',SO3DirichletKernel(bw),'harmonic');

%% step 4: convolution
GBPD = conv(odf,psi);

%plot(GBPD)

end

author: mtex-toolbox

repo: mtex

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


text:
classdef nme_bus_nld_acp_node_test < mp.nme_bus_acp

%   MATPOWER
%   Copyright (c) 2021-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end

    methods
        function name = name(obj)
            name = 'bus_nld';
        end

        function obj = add_vvars(obj, nm, dm, idx)
            dme = obj.data_model_element(dm);
            nb = obj.nk;

            %% prepare angle bounds for ref buses
            va_lb = -Inf(nb, 1);
            va_ub =  Inf(nb, 1);
            k = find(dme.type == mp.NODE_TYPE.REF);
            va_lb(k) = dme.va_start(k);
            va_ub(k) = dme.va_start(k);

            nm.add_var('va', ['va_' obj.name], nb, dme.va_start, va_lb, va_ub);
            nm.add_var('vm', ['vm_' obj.name], nb, dme.vm_start, dme.vm_lb, dme.vm_ub);
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
classdef LaunchVehicleSensorReport < matlab.mixin.SetGet
    %LaunchVehicleSensorReport Summary of this class goes here
    %   Detailed explanation goes here

    properties
        sensor AbstractSensor
        lvdData LvdData
    end

    methods
        function obj = LaunchVehicleSensorReport(sensor, lvdData)
            obj.sensor = sensor;
            obj.lvdData = lvdData;
        end

        function [filepath, coverageT, azT, elT, rngT, angleToBoresightT] = generateReport(obj, targets, progressBar, busyStatusLabel, sensorNum, totalSensorNum, reportFolder)
            arguments
                obj(1,1) LaunchVehicleSensorReport
                targets(1,:) AbstractSensorTarget
                progressBar(1,1) wt.ProgressBar
                busyStatusLabel(1,1) matlab.ui.control.Label
                sensorNum(1,1) double
                totalSensorNum(1,1) double
                reportFolder(1,1) string
            end

            progressBar.Indeterminate = false;
            progressBar.startProgress('');
            busyStatusLabel.Text = sprintf('Generating Data for Sensor "%s" [%u of %u]', obj.sensor.getListboxStr(), sensorNum, totalSensorNum);

            bodyInfos = obj.lvdData.celBodyData.getAllBodyInfo();
            stateLog = obj.lvdData.stateLog;

            entries = stateLog.getAllEntries();
            times = [entries.time]';
            coverages = cell(1,numel(targets));
            sensorAzs = cell(1,numel(targets));
            sensorEls = cell(1,numel(targets));
            sensorRngs = cell(1,numel(targets));
            sensorAnglesToBoresight = cell(1,numel(targets));
            for(i=1:length(entries))
                entry = entries(i);

                sensorState = entry.getSensorStateForSensor(obj.sensor);
                scElem = entry.getCartesianElementSetRepresentation(false);
                dcm = entry.attitude.dcm;
                frame = entry.centralBody.getBodyCenteredInertialFrame();

                results = obj.sensor.evaluateSensorTargets(sensorState, targets, scElem, dcm, bodyInfos, frame);

                for(j=1:length(results))
                    result = results(j);

                    coverages{j} = vertcat(coverages{j}, result.resultsBool(:)');

                    [az, el, rng, angle] = result.getBoresightRelativeAngles(sensorState, scElem, dcm);

                    sensorAzs{j} = vertcat(sensorAzs{j}, az);
                    sensorEls{j} = vertcat(sensorEls{j}, el);
                    sensorRngs{j} = vertcat(sensorRngs{j}, rng);
                    sensorAnglesToBoresight{j} = vertcat(sensorAnglesToBoresight{j}, angle);
                end

                progressBar.setProgress(i/length(entries), sprintf('%0.3f %%', 100*i/length(entries)));
            end

            progressBar.Indeterminate = true;
            busyStatusLabel.Text = sprintf('Writing Report File for Sensor "%s" [%u of %u]', obj.sensor.getListboxStr(), sensorNum, totalSensorNum);
            progressBar.setProgress(1, '');
            drawnow;

            coverageT = {};
            filename = sprintf('SensorReport_%s_%s.xls', obj.sensor.getListboxStr(),datestr(now, 'YYYYmmDD_HHMMss'));
            filepath = fullfile(reportFolder, filename);
            for(i=1:length(coverages))
                coverage = coverages{i};

                %compute instant and cumulative coverage
                instantCoverage = NaN(height(coverage), 1);
                cumCoverage = NaN(height(coverage), 1);
                for(j=1:height(coverage))
                    subCoverage = coverage(j,:);
                    instantCoverage(j) = sum(subCoverage) / numel(subCoverage);

                    subCoverage = coverage(1:j,:);
                    sumSubCoverage = sum(subCoverage,1);
                    sumSubCoverage(sumSubCoverage >= 1) = 1;
                    cumCoverage(j) = sum(sumSubCoverage) / numel(sumSubCoverage);
                end

                target = targets(i);
                labelStrs = target.getTargetPtLabelStrs();

                %Info worksheet
                sensorNums = sensorNum * ones(numel(targets),1);
                sensorNames = repmat(string(obj.sensor.getListboxStr()), numel(targets), 1);
                targetNums = [1:numel(targets)]'; %#ok<NBRAK>
                for(j=1:length(targets))
                    targetNames(j,1) = string(targets(j).getListboxStr()); %#ok<AGROW>
                end

                varNames = {'Sensor Number', 'Sensor Name', 'Target Number', 'Target Name'};
                data = table(sensorNums, sensorNames, targetNums, targetNames, 'VariableNames',varNames);
                writetable(data, filepath, 'FileType','spreadsheet', 'Sheet','Information');

                %Sensor Coverage
                worksheetName = sprintf('Coverage, Sensor %u Target %u', sensorNum, i);

                header = horzcat('Time [UT sec]', labelStrs, "Instantaneous Coverage Fraction", "Cumulative Coverage Fraction");
                data = horzcat(times, coverage, instantCoverage, cumCoverage);

                coverageT{i} = writeDataTableToXLSFile(data, header, filepath, worksheetName); %#ok<AGROW>

                %Sensor to Target Az
                worksheetName = sprintf('Azimuth, Sensor %u Target %u', sensorNum, i);

                header = horzcat('Time [UT sec]', labelStrs);
                data = horzcat(times, rad2deg(sensorAzs{i}));

                azT{i} = writeDataTableToXLSFile(data, header, filepath, worksheetName); %#ok<AGROW>

                %Sensor to Target El
                worksheetName = sprintf('Elevation, Sensor %u Target %u', sensorNum, i);

                header = horzcat('Time [UT sec]', labelStrs);
                data = horzcat(times, rad2deg(sensorEls{i}));

                elT{i} = writeDataTableToXLSFile(data, header, filepath, worksheetName); %#ok<AGROW>

                %Sensor to Target Range
                worksheetName = sprintf('Range, Sensor %u Target %u', sensorNum, i);

                header = horzcat('Time [UT sec]', labelStrs);
                data = horzcat(times, sensorRngs{i});

                rngT{i} = writeDataTableToXLSFile(data, header, filepath, worksheetName); %#ok<AGROW>

                %Sensor to Target Range
                worksheetName = sprintf('Angle, Sensor %u Target %u', sensorNum, i);

                header = horzcat('Time [UT sec]', labelStrs);
                data = horzcat(times, rad2deg(sensorAnglesToBoresight{i}));

                angleToBoresightT{i} = writeDataTableToXLSFile(data, header, filepath, worksheetName); %#ok<AGROW>
            end
        end
    end
end

function T = writeDataTableToXLSFile(data, header, filename, worksheetName)
    T = array2table(data);
    T.Properties.VariableNames = header;

    userdata.WorksheetName = worksheetName;
    T.Properties.UserData = userdata;

    writetable(T, filename, 'FileType','spreadsheet', 'Sheet',worksheetName);
end


text:
function [FX,dFdX,dFdP] = f_embedAR(X,P,ut,in)
% AR(1) embedding evolution function

n = size(X,1);
x = X(1:n/2);
z = X(n/2+1:n);
% call native evolution function
[fx,J,dfdP] = VBA_evalFun('f',x,P,ut,in.opt,in.dim,0);
% add AR(1) noise and form derivatives
FX = [fx+z;z];
dFdX = [ J          zeros(n/2)
         eye(n/2)   eye(n/2) ];
dFdP = [dfdP,zeros(size(P,1),n/2)];

author: MBB-team

repo: VBA-toolbox


text:
function [inplane ok] = roiLoadVol2Inplane(inplane, roiName, saveFlag)
%
% [inplane ok] = roiLoadVol2Inplane(inplane, <roiName>, <saveFlag=1>);
%
% A convenient accessor function: for an inplane view,
% load an ROI defined on the volume, and convert it over into
% the inplanes. Saves loading a separate Volume view, if this is
% all you want to do.
%
% If the ROI name is omitted, will prompt the user.
%
% If saveFlag=1, will save the xformed inplane ROI (provided it
% is nonempty).
%
% ras, 05/06.
if notDefined('inplane'), inplane = getSelectedInplane; end
if notDefined('saveFlag'), saveFlag = 1;                end

if ~isequal(inplane.viewType, 'Inplane')
    error('Works only on inplane views.')
end

hV = getSelectedGray;

if notDefined('roiName')
    % dialog
    roiName = getROIfilename(hV, 0);
end

% ensure cell ROI name
if ~iscell(roiName), roiName = {roiName};  end

for i = 1:length(roiName)
    [hV, ok] = loadROI(hV, roiName{i}, 1, [], 1);

    if ~ok
        fprintf('Couldn''t load Volume ROI %s \n', roiName{i});
        continue
    end

    ipROI = vol2ipROI(hV.ROIs(end), hV, inplane);

	% is the ROI empty? Go ahead and xform it, but warn the user and set
	% the 'ok' status to zero:
	if isempty(ipROI.coords)
		warning('[%s]: Xformed empty ROI %s.', mfilename, ipROI.name);
		ok = 0;
	end

    inplane = addROI(inplane, ipROI, 0);

    if saveFlag & ~isempty(ipROI.coords)
        inplane = saveROI(inplane, inplane.ROIs(end), [], saveFlag);
    end
end

inplane = selectROI(inplane, length(inplane.ROIs));

clear hV

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function value = p05_i4 ( action, name, value )

%*****************************************************************************80
%
%% P05_I4 sets or gets I4 parameters for problem 05.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ACTION,
%    'D' to set a parameter to its default value.
%    'S' to set a parameter,
%    'G' to get a parameter,
%    'I' to increment a parameter.
%
%    Input, string NAME, the name of the variable.
%    '#' is the number of calls to the integrand routine.
%
%    Input/output, integer VALUE.
%    * If ACTION = 'I', then VALUE is an input quantity, and is the
%      new value to be added to NAME.
%    * If ACTION = 'S', then VALUE is an input quantity, and is the
%      new value to be assigned to NAME.
%    * If ACTION = 'G', then VALUE is an output quantity, and is the
%      current value of NAME.
%
  persistent calls;

  if ( action == 'D' | action == 'd' )

    if ( name == '#' | name == '*' )
      calls = 0;
    end

  elseif ( action == 'G' | action == 'g' )

    if ( name == '#' )
      value = calls;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P05_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P05_I4 - Fatal error!' );
    end

  elseif ( action == 'I' | action == 'i' )

    if ( name == '#' )
      calls = calls + value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P05_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P05_I4 - Fatal error!' );
    end

  elseif ( action == 'S' | action == 's' )

    if ( name == '#' )
      calls = value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P05_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P05_I4 - Fatal error!' );
    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P05_I4 - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized action = "%s".\n', action );
    error ( 'P05_I4 - Fatal error!' );

  end

  return
end


text:
function sdd = settings_detection_default()
% Default settings for detection.
sdd.thr_method = 'median'; % median, wavelet, plexon
sdd.side = 'negative'; % two, positive, negative, mean
sdd.spike_filter_type = 'butter';
sdd.spike_filter_order = 4;
sdd.spike_Rp = 300; % Hz
sdd.spike_Rs = 3000; % Hz
sdd.pre_thresh = 40; % samples
sdd.post_thresh = 59; % samples
sdd.dead_time = 20; % samples
sdd.sampling_rate = 40000; % sampling rate

write_arguments_to_json(sdd,'config_detection');

author: ramintoosi

repo: ROSS

save_path: github-repos/MATLAB/ramintoosi-ROSS


text:
function download_finviz_data(destination_dir)
clearvars -except destination_dir
source_dir = '/Users/ansonwong/Downloads';
source_filename=strcat(source_dir,'/finviz.csv');
destination_filename=strcat(destination_dir,'/finviz.csv');

% Download finviz data using my default browser
finviz_URL_full=...
  strcat('http://finviz.com/export.ashx?v=152&ft=4',...
            '&c=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72');
web(finviz_URL_full,'-browser');

% Wait for the download to finish
fprintf('Downloaded ''%s'' (success=%d)\n',source_filename,wait_for_existence(source_filename));

% Now that the download has finished, move the file to our destination directory
[status,message,messageid]=movefile(source_filename,destination_filename);
if(status==0), error('Error in moving finviz file!\n'); end
end

function success = wait_for_existence(filename)
success = 1;

seconds_elapsed=0;
wait_seconds=0.5;
while( exist(filename)==0 ),
  pause(wait_seconds);
  seconds_elapsed=seconds_elapsed+wait_seconds;
  if(seconds_elapsed > 120),
    error('Waited too long for %s to exist!\n',filename)
    success=0; return
  end
end

fprintf('Waited %.1f seconds.\n',seconds_elapsed);
end

author: gudbrandtandberg


text:
function [h,p] = gppredcheck(gpstats,alpha)
%GPPREDCHECK Check calibration of Gaussian process prediction.

n = gpstats.last;
idx = 1:n;
p = NaN;

% No predictions available, test failed
if n == 0; h = 1; return; end

zscores = (gpstats.fval(idx) - gpstats.ymu(:,idx))./(gpstats.ys(:,idx));

% A NaN means that some prediction went wrong
if any(isnan(zscores)); h = 1; return ; end

zscores(isnan(zscores)) = [];
n = numel(zscores);

if n < 3
    mychi2inv = @(x,v) 2*gammaincinv(x,v/2);
    plo = mychi2inv(alpha/2,n);
    phi = mychi2inv(1-alpha/2,n);
    t = sum(zscores.^2);
    if t < plo || t > phi || any(isnan([plo,phi]))
        h = 1;
    else
        h = 0;
    end
else
    [h,p] = swtest(zscores, alpha);
end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function y=modulo(x,N);
%MODULO	Congruence of a vector.
%	Y=MODULO(X,N) gives the congruence of each element of the
%	vector X modulo N. These values are strictly positive and
%	lower equal than N.
%
%	See also REM.

%	O. Lemoine - February 1996.
%	Copyright (c) 1996 by CNRS (France).
%
%  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 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if isreal(x)
  y=mod(x,N);
  idx=find(y==0);
  y(idx)=N;
else
  y=mod(real(x),N)+i*mod(imag(x),N);
end;

author: HeLiangHIT


text:
function Hv = lbfgsHvFunc2(v,Hdiag,N,M)
Hv = v/Hdiag - N*(M\(N'*v));

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function [ imgOut ] = hyperOrthorectify( imgIn, altitude, hpbw )
%HYPERORTHORECTIFY Orthorectifies areal observed data.
%   Orthorectifies areal observed data using nearest neighbor interpolation.
%
% Inputs
%   imgIn       Input image (m x n) or (m x n x p)
%   altitude    Sensor altitude (meters)
%   hpbw        Half power beam width (radians).
% Outputs
%   imgOut      Orthorectified image.

% Input parameters
if (ndims(imgIn) == 2)
    [h, w] = size(imgIn);
    p = 1;
elseif (ndims(imgIn) == 3)
    [h, w, p] = size(imgIn);
end

radPerPix = hpbw/w;
x = tan(hpbw/2)*altitude;  % m
gsd = altitude*radPerPix;  % m
n = x/gsd;

outImg = zeros(h, floor(n)*2, p);
for k=1:p
    for j=1:h
        for i=-floor(n):1:floor(n)-1
            boresiteDistance = gsd*i;
            theta = atan(boresiteDistance/altitude);
            imagePix = round(theta / radPerPix);
            imgOut(j, floor(n)+i+1, k) = imgIn(j, (w/2)+imagePix+1, k);
        end
    end
end

author: isaacgerg

repo: matlabHyperspectralToolbox


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR:  rotating an US image
%
%   - load US data
%   - set image viewer and interpolator and trafo=rigid2D
%   - run a loop over different angles
%==============================================================================

clear, close all, help(mfilename)

%% load data
Tdata = double(imread('US.jpg'));
omega = [0,size(Tdata,1),0,size(Tdata,2)];
m     = [192,128] ;
xc    = getCellCenteredGrid(omega,m);

% setup interpolation scheme and image viewer
imgModel('reset','imgModel','linearInter');
viewImage('reset','viewImage','viewImage2D','colormap',gray(256),'axis','off');
fprintf('%20s : %s\n','viewImage',viewImage);
fprintf('%20s : %s\n','image model',imgModel);

% shortcuts for plotting stuff
Grid   = @(X)   plotGrid(X,omega,m,'spacing',8,'color','w');
dimstr = @(m)   sprintf('%s=[%s]',inputname(1),sprintf(' %d',m));
Title  = @(s,t) title([s,sprintf(', t=%s',num2str(t))],'fontsize',30);

%% display initial image
FAIRfigure(1,'figname',mfilename); clf; subplot(1,2,1);
viewImage(imgModel(Tdata,omega,xc),omega,m); hold on; gh = Grid(xc);
title(sprintf('%s, %s','data',dimstr(m)),'fontsize',30);

%% initialize transformation
trafo('set','trafo','rigid2D');
S = @(t) t/6;

% y(t) = R(t)*x+(I-R(t))*c,
%
% R(t) = [cos(t) -sin(t) ], c = [omega(1)]/2
%        [sin(w)  cos(t) ]      [omega(2)]
%
% or y(t) = rigid2D(w,x), where w = [t;(I-R(t))omega'/2].
c  = (omega(2:2:end)+omega(1:2:end))'/2
wc = @(t) [t;(eye(2)-[ cos(t),-sin(t);sin(t),cos(t)])*c];

%=======================================================================================
% the loop over time t
%=======================================================================================
t = 0.35*pi*sin(linspace(0,2.5*pi,101));

for j=1:length(t),
  yc = trafo(wc(t(j)),xc);      % compute the transformed points
  Tc = imgModel(Tdata,omega,yc);   % compute the transformed image

  if j == 1,              % initialize visualization
    figure(1); subplot(1,2,1); %colordef(gcf,'black');
    set(gh,'visible','off'); gh = Grid(yc);
    title(sprintf('%s, %s','data',dimstr(m)),'fontsize',30);
    subplot(1,2,2); vh = viewImage(Tc,omega,m); Title(trafo,t(j));
    FAIRpause;
  else,                   % continue plots
    subplot(1,2,1); set(gh,'visible','off');        gh = Grid(yc);
    axis(omega)
    subplot(1,2,2); set(vh,'cdata',reshape(Tc,m)'); Title(trafo,t(j));
    FAIRpause(1/2000)
  end;
  drawnow; fprintf('.'); if rem(j,50) == 0, fprintf('\n'); end;
end;
fprintf('\n')
%=======================================================================================


text:
function [stringTokens, numberOfTokens] = stringTokenizer( stringTarget, delimiter )
%STRINGTOKENIZER
% stringTarget = ';;ab;cd de=0==';
% delimiter = '; =';

% regexpression
pattern = ['\w*[^' delimiter ']\w*'] ;
[startIndices, endIndices] = regexp(stringTarget, pattern);

stringTokens = [];
index = 1;

for i = 1:length(startIndices)
    stringTokens = [stringTokens {stringTarget(startIndices(i):endIndices(i))}];
end

numberOfTokens = length(stringTokens);

end

author: Sable

repo: mcbench-benchmarks


text:
function [ A ] = calcu_similarity_matrix( X );

% X, data matrix, each ROW contains one signal

[n, p] = size(X);
nn = 7; kk = floor(log2(n)) + 1; %default value
D = dist2(X, X); %distance matrix

A = scale_dist3_knn(D, nn, kk, true);
dd = 1 ./ sum(A);
dd = sqrt(dd);
A = bsxfun(@times, A, dd);
A = A';
A = bsxfun(@times, A, dd);

A = (A + A') / 2;

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function NoiseCovFiles = bst_noisecov(iTargetStudies, iDataStudies, iDatas, Options, isDataCov)
% BST_NOISECOV: Compute noise covariance matrix for a set of studies.
%
% USAGE:  NoiseCovFiles = bst_noisecov(iTargetStudies, iDataStudies, iDatas, Options=[ask], isDataCov=0)
%               Options = bst_noisecov()
%
% INPUT:
%     - iTargetStudies : List of studies indices for which the noise covariance matrix is produced
%     - iDataStudies   : [1,nData] int, List of data files to use for computation (studies indices)
%                        If not defined or [], uses all the recordings from all the studies (iTargetStudies)
%     - iDatas         : [1,nData] int, List of data files to use for computation (data indices)
%     - isDataCov      : If 1, saves the result as the data covariance, if 0 saves it as the noise covariance
%     - Options        : Structure with the following fields (if not defined: asked to the user)
%           |- Baseline        : [tStart, tStop]; range of time values considered as baseline
%           |- RemoveDcOffset  : {'file', 'all'}; 'all' removes the baseline avg file by file; 'all' computes the baseline avg from all the files
%           |- ReplaceFile     : If 1 replaces automatically the previous noisecov file without asking the user for a confirmation
%           |- ChannelTypes    : Cell array with the list of selected channel types

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2009-2022

%% ===== RETURN DEFAULT OPTIONS =====
% Options structure
if (nargin == 0)
    NoiseCovFiles = struct(...
        'Baseline',        [-.1, 0], ...
        'DataTimeWindow',  [0, 1], ...
        'RemoveDcOffset',  'file', ...
        'ReplaceFile',     [], ...
        'ChannelTypes',    []);
    return;
end

%% ===== PARSE INPUTS =====
if (nargin < 5) || isempty(isDataCov)
    isDataCov = [];
end
if (nargin < 4) || isempty(Options)
    Options = [];
end
NoiseCovFiles = {};
isRaw = 0;
% Get source files
if (nargin < 3) || isempty(iDataStudies) || (length(iDataStudies) ~= length(iDatas))
    % Get all the datafiles depending on the target studies
    sStudies = bst_get('Study', iTargetStudies);
    sDatas = [sStudies.Data];
    DataFiles = {sDatas.FileName};
else
    % Unique studies
    uniqueStudies = unique(iDataStudies);
    DataFiles = {};
    % For each study
    for i = 1:length(uniqueStudies)
        % Get study
        iStudy = uniqueStudies(i);
        sStudy = bst_get('Study', iStudy);
        iDataLocal = iDatas(iDataStudies == iStudy);
        % Add files to list
        DataFiles = cat(2, DataFiles, {sStudy.Data(iDataLocal).FileName});
        % Raw data
        isRaw = any(strcmpi({sStudy.Data(iDataLocal).DataType}, 'raw'));
    end
end
% Default file tag
if isDataCov
    strComment = 'Data covariance';
else
    strComment = 'Noise covariance';
end

%% ===== GET DATA CHANNELS =====
% Get channel studies
sTargetStudies = bst_get('Study', iTargetStudies);
% Find a study with channel file
iWithChan = find(~cellfun(@isempty, {sTargetStudies.Channel}), 1);
if isempty(iWithChan)
    error('No channel file.');
end
% Get one channel study
sStudy = sTargetStudies(iWithChan);
% Load channel file
ChannelMat = in_bst_channel(sStudy.Channel.FileName);
% Get all the valid channels
iChan = good_channel(ChannelMat.Channel, [], {'MEG','EEG','ECOG','SEEG'});
nChanAll = length(ChannelMat.Channel);
% Get the possible channe types
ChannelTypes = unique({ChannelMat.Channel(iChan).Type});

%% ===== READ ALL TIME VECTORS =====
% Progress bar
bst_progress('start', 'Read recordings information', 'Analysing input files...', 0, length(DataFiles));
% Regular list of imported data files
if ~isRaw
    nFiles = length(DataFiles);
    DataMats = repmat(struct('Time', [], 'SamplingRate', [], 'nAvg', [], 'Leff', [], 'iEpoch', [], 'iBadTime', []), nFiles);
    % Loop on all input files
    for iFile = 1:length(DataFiles)
        % Load file metadata
        DataMat = in_bst_data(DataFiles{iFile}, 'Time', 'nAvg', 'Leff');
        % Check file time
        if (length(DataMat.Time) < 3)
            error(['File has no time dimension: ' DataFiles{iFile}]);
        end
        % Check sampling rate
        if (iFile == 1)
            SamplingRate = DataMat.Time(2) - DataMat.Time(1);
        elseif (abs(SamplingRate - (DataMat.Time(2) - DataMat.Time(1))) > 1e-6)
            error(['The files you selected have different sampling frequencies. They should not be processed together.' 10 ...
                   'Please only select recordings with the same sampling frequency.']);
        end
        % Save values
        DataMats(iFile).Time = DataMat.Time;
        DataMats(iFile).nAvg = DataMat.nAvg;
        DataMats(iFile).Leff = DataMat.Leff;
        bst_progress('inc', 1);
    end
% Raw file
else
    % Only one raw file allowed
    if (length(DataFiles) > 1)
        error('Only one raw file allowed for covariance computation');
    end
    nFiles = 1;
    % Read the description of the raw file
    RawMat = in_bst_data(DataFiles{1});
    sFile = RawMat.F;
    clear RawMat;
    % Get bad segments/epochs
    [badSeg, badEpochs] = panel_record('GetBadSegments', sFile);

    % Initialize DataMats structure
    DataMats = repmat(struct('Time', [], 'nAvg', [], 'Leff', [], 'iEpoch', [], 'iBadTime', []), 0);
    % Define size of blocks to read
    MAX_BLOCK_SIZE = 10000;
    % Loop on epochs
    for iEpoch = 1:max(length(sFile.epochs), 1)
        % Bad epoch
        if ~isempty(sFile.epochs) && sFile.epochs(iEpoch).bad
            disp(sprintf('NOISECOV> Ignoring epoch #%d (tagged as bad)', iEpoch));
        end
        % Get the bad segments for this epoch
        iBadEpoch = find(badEpochs == iEpoch);
        if ~isempty(iBadEpoch)
            badSegEpoc = badSeg(:, iBadEpoch);
        else
            badSegEpoc = [];
        end
        % Get total number of samples
        if ~isempty(sFile.epochs)
            samples = round(sFile.epochs(iEpoch).times .* sFile.prop.sfreq);
            nAvg = double(sFile.epochs(iEpoch).nAvg);
        else
            samples = round(sFile.prop.times .* sFile.prop.sfreq);
            nAvg = 1;
        end
        totalSmpLength = double(samples(2) - samples(1)) + 1;
        % Number of blocks to split this epoch in
        nbBlocks = ceil(totalSmpLength / MAX_BLOCK_SIZE);
        % For each block
        for iBlock = 1:nbBlocks
            % Get samples indices for this block (start ind = 0)
            smpBlock = samples(1) + [(iBlock - 1) * MAX_BLOCK_SIZE, min(iBlock * MAX_BLOCK_SIZE - 1, totalSmpLength - 1)];
            smpList = smpBlock(1):smpBlock(2);
            % Create a data block
            iNew = length(DataMats) + 1;
            DataMats(iNew).iEpoch = iEpoch;
            DataMats(iNew).Time   = smpList ./ sFile.prop.sfreq;
            DataMats(iNew).nAvg   = nAvg;
            DataMats(iNew).Leff   = nAvg;
            % Remove the portions that have bad segments in them
            iBadTime = [];
            for ix = 1:size(badSegEpoc, 2)
                iBadTime = [iBadTime, find((smpList >= badSegEpoc(1,ix)) & (smpList <= badSegEpoc(2,ix)))];
            end
            if ~isempty(iBadTime)
                DataMats(iNew).iBadTime = iBadTime;
            end
        end
    end
    SamplingRate = 1 ./ sFile.prop.sfreq;
end
% Get number of samples and sampling rates
nSamples = length([DataMats.Time]) - length([DataMats.iBadTime]);
% Close progress bar
bst_progress('stop');
% Check number of actual samples available for computation
if (nSamples == 0)
    error('This selection does not contain any file that can be used for computing the covariance matrix.');
end

%% ===== GET NOISECOV OPTIONS ======
% Get number of files
nBlocks = length(DataMats);
% If the options were not passed in argument
if isempty(Options)
    % Loop to get all the valid times
    allTimes = [];
    for iFile = 1:length(DataMats)
        tmpTime = DataMats(iFile).Time;
        tmpTime(DataMats(iFile).iBadTime) = [];
        allTimes = [allTimes, tmpTime];
    end
    % Prepare GUI options
    guiOptions.timeSamples  = sort(allTimes);
    guiOptions.nFiles       = nFiles;
    guiOptions.nBlocks      = nBlocks;
    guiOptions.nChannels    = length(iChan);
    guiOptions.freq         = 1 ./ SamplingRate;
    guiOptions.isDataCov    = isDataCov;
    guiOptions.ChannelTypes = ChannelTypes;
    % Display dialog window
    Options = gui_show_dialog(strComment, @panel_noisecov, 1, [], guiOptions);
    if isempty(Options)
        return
    end
end
% Get only the channels of the selected modality
if ~isempty(Options.ChannelTypes)
    iChan = good_channel(ChannelMat.Channel, [], Options.ChannelTypes);
end
if isempty(iChan)
    error('No channels are selected from this file.');
end

%% ===== COMPUTE AVERAGE/TIME =====
if strcmpi(Options.RemoveDcOffset, 'all')
    % Compute the average across ALL the time samples of ALL the files
    Favg = zeros(nChanAll, 1);
    Ntotal = zeros(nChanAll, 1);
    % Progress bar
    bst_progress('start', 'Average across time', 'Computing average across time...', 0, nBlocks);
    % Loop on all the files
    for iFile = 1:nBlocks
        bst_progress('inc', 1);
        % Load recordings
        [DataMat, iTimeBaseline] = ReadRecordings();
        if isempty(iTimeBaseline)
            continue
        end
        % Get good channels
        iGoodChan = intersect(find(DataMat.ChannelFlag == 1), iChan);

        % === Compute average ===
        % Favg(iGoodChan)   = Favg(iGoodChan)   + double(DataMat.nAvg) .* sum(DataMat.F(iGoodChan,iTimeBaseline),2);
        % Ntotal(iGoodChan) = Ntotal(iGoodChan) + double(DataMat.nAvg) .* length(iTimeBaseline);
        Favg(iGoodChan)   = Favg(iGoodChan)   + double(DataMat.Leff) .* sum(DataMat.F(iGoodChan,iTimeBaseline),2);
        Ntotal(iGoodChan) = Ntotal(iGoodChan) + double(DataMat.Leff) .* length(iTimeBaseline);
    end
    % Remove zero-values in Ntotal
    Ntotal(Ntotal == 0) = 1;
    % Divide each channel by total number of time samples
    Favg = Favg ./ Ntotal;
end

%% ===== COMPUTE NOISE COVARIANCE =====
Ntotal   = zeros(nChanAll);
NoiseCov = zeros(nChanAll);
FourthMoment = zeros(nChanAll);
% Progress bar
bst_progress('start', strComment, ['Computing ' lower(strComment) '...'], 0, nBlocks);
drawnow
% Loop on all the files
for iFile = 1:nBlocks
    bst_progress('inc', 1);
    % Load recordings
    [DataMat, iTimeBaseline, iTimeCov] = ReadRecordings();
    if isempty(iTimeBaseline)
        continue
    end
    N = length(iTimeCov);

    % === Compute average ===
    if strcmpi(Options.RemoveDcOffset, 'file')
        % Get good channels
        iGoodChan = intersect(find(DataMat.ChannelFlag == 1), iChan);
        % Average baseline values
        Favg = mean(DataMat.F(:,iTimeBaseline), 2);
    end

    % === Compute covariance ===
    % Remove average
    DataMat.F(iGoodChan,:) = bst_bsxfun(@minus, DataMat.F(iGoodChan,:), Favg(iGoodChan,1));
    % Compute covariance for this file
    % fileCov  = DataMat.nAvg .* (DataMat.F(iGoodChan,iTimeCov)    * DataMat.F(iGoodChan,iTimeCov)'   );
    % fileCov2 = DataMat.nAvg .* (DataMat.F(iGoodChan,iTimeCov).^2 * DataMat.F(iGoodChan,iTimeCov)'.^2);
    fileCov  = DataMat.Leff .* (DataMat.F(iGoodChan,iTimeCov)    * DataMat.F(iGoodChan,iTimeCov)'   );
    fileCov2 = DataMat.Leff .* (DataMat.F(iGoodChan,iTimeCov).^2 * DataMat.F(iGoodChan,iTimeCov)'.^2);
    % Add file covariance to accumulator
    NoiseCov(iGoodChan,iGoodChan)     = NoiseCov(iGoodChan,iGoodChan)     + fileCov;
    FourthMoment(iGoodChan,iGoodChan) = FourthMoment(iGoodChan,iGoodChan) + fileCov2;
    Ntotal(iGoodChan,iGoodChan) = Ntotal(iGoodChan,iGoodChan) + N;
end
% Remove zeros from N matrix
Ntotal(Ntotal <= 1) = 2;
% Divide final matrix by number of samples
NoiseCov     = NoiseCov     ./ (Ntotal - 1);
FourthMoment = FourthMoment ./ (Ntotal - 1);
% Display result in the command window
nSamplesTotal = max(Ntotal(:));
% disp(['Number of time samples used for the noise covariance: ' num2str(nSamplesTotal)]);
% Check for NaN values
if (nnz(isnan(NoiseCov)) > 0)
    error('The output covariance contains NaN values. Please check your recordings and tag the bad channels correctly.');
end

%% ===== IMPORTING IN DATABASE =====
% Build file structure
NoiseCovMat = db_template('noisecovmat');
NoiseCovMat.NoiseCov     = NoiseCov;
NoiseCovMat.Comment      = [strComment ': '];
NoiseCovMat.nSamples     = Ntotal;
NoiseCovMat.FourthMoment = FourthMoment;
% Add names of sensors
allTypes = unique({ChannelMat.Channel(iChan).Type});
if all(ismember({'MEG MAG', 'MEG GRAD'}, allTypes))
    allTypes = setdiff(allTypes, {'MEG MAG', 'MEG GRAD'});
    allTypes = union(allTypes, {'MEG'});
end
for i = 1:length(allTypes)
    NoiseCovMat.Comment = [NoiseCovMat.Comment, allTypes{i}];
    if (i < length(allTypes))
        NoiseCovMat.Comment = [NoiseCovMat.Comment, ', '];
    end
end
% Add history entry
if isDataCov
    if isempty(Options.DataTimeWindow)
        strTime = 'Data=[All file], ';
    elseif (max(abs(Options.DataTimeWindow)) > 2)
        strTime = sprintf('Data=[%1.3f, %1.3f]s, ', Options.DataTimeWindow);
    else
        strTime = sprintf('Data=[%d, %d]ms, ', round(Options.DataTimeWindow * 1000));
    end
    if isempty(Options.Baseline)
        strTime = [strTime, 'Baseline=[All file]'];
    elseif (max(abs(Options.Baseline)) > 2)
        strTime = [strTime, sprintf('Baseline=[%1.3f, %1.3f]s', Options.Baseline)];
    else
        strTime = [strTime, sprintf('Baseline=[%d, %d]ms', round(Options.Baseline * 1000))];
    end
else
    if isempty(Options.Baseline)
        strTime = '[All file]';
    elseif (max(abs(Options.Baseline)) > 2)
        strTime = sprintf('[%1.3f, %1.3f]s', Options.Baseline);
    else
        strTime = sprintf('[%d, %d]ms', round(Options.Baseline * 1000));
    end
end
NoiseCovMat = bst_history('add', NoiseCovMat, 'compute', sprintf('Computed based on %d files (%d blocks, %d samples): %s, %s', nFiles, nBlocks, nSamplesTotal, strTime, Options.RemoveDcOffset));
% Add list of input files in history
for iSrc = 1:length(DataFiles)
    NoiseCovMat = bst_history('add', NoiseCovMat, 'src', DataFiles{iSrc});
end
% Save in database
NoiseCovFiles = import_noisecov(iTargetStudies, NoiseCovMat, Options.ReplaceFile, isDataCov);
% Close progress bar
bst_progress('stop');

%% ========================================================================
%  ======= SUPPORT FUNCTIONS ==============================================
%  ========================================================================

%% ===== READ RECORDINGS BLOCK =====
    function [DataMat, iTimeBaseline, iTimeCov] = ReadRecordings()
        iTimeBaseline = [];
        iTimeCov = [];
        if ~isRaw
            bst_progress('text', ['File: ' DataFiles{iFile}]);
            DataMat = in_bst_data(DataFiles{iFile}, 'F', 'ChannelFlag', 'Time', 'nAvg', 'Leff');
        else
            DataMat = DataMats(iFile);
            % If file is block does not contain any baseline segment: skip it
            if (length(DataMats(iFile).Time) < 2) || (~isempty(Options.Baseline) && ((DataMats(iFile).Time(end) < Options.Baseline(1)) || (DataMats(iFile).Time(1) > Options.Baseline(2))))
                return;
            end
            % Read raw block
            UseSsp = 1;
            [DataMat.F, DataMat.Time] = panel_record('ReadRawBlock', sFile, ChannelMat, DataMats(iFile).iEpoch, DataMats(iFile).Time([1,end]), 0, 1, 'no', UseSsp);
            DataMat.ChannelFlag = sFile.channelflag;
        end
        % WE DON'T WANT TO APPLY THE AVERAGE REFERENCE AT THIS STAGE: NOISE COV IS CALCULATED ON ORIGINAL MONTAGE, AND AVG REF IS CALCULATED IN PROCESS_INVERSE.
%         % Apply average reference: separately SEEG, ECOG, EEG
%         if any(ismember(unique({ChannelMat.Channel.Type}), {'EEG','ECOG','SEEG'}))
%             sMontage = panel_montage('GetMontageAvgRef', ChannelMat.Channel, DataMat.ChannelFlag);
%             DataMat.F = sMontage.Matrix * DataMat.F;
%         end
        % If not enough time frames (ie. if data files)
        if (length(DataMat.Time) <= 2) || (size(DataMat.F,2) <= 2)
            return;
        end
        % Check size
        if (size(DataMat.F,1) ~= nChanAll)
            error('Number of channels is not constant.');
        end
        % Get the times required in the block (only when there are bad segments in raw files)
        if isRaw && ~isempty(DataMat.iBadTime)
            DataMat.Time(DataMat.iBadTime) = [];
            DataMat.F(:, DataMat.iBadTime) = [];
        end
        % Get times that are considered as baseline
        if isempty(DataMat.Time)
            iTimeBaseline = [];
        elseif ~isempty(Options.Baseline)
            iTimeBaseline = panel_time('GetTimeIndices', DataMat.Time, Options.Baseline);
        else
            iTimeBaseline = 1:length(DataMat.Time);
        end
        % Get the time indices on which to compute the covariance
        if isDataCov
            if isempty(DataMat.Time)
                iTimeCov = [];
            elseif ~isempty(Options.DataTimeWindow)
                iTimeCov = panel_time('GetTimeIndices', DataMat.Time, Options.DataTimeWindow);
            else
                iTimeCov = 1:length(DataMat.Time);
            end
        else
            iTimeCov = iTimeBaseline;
        end
    end
end

author: brainstorm-tools


text:
function clebsch_gordan_values_test ( )

%*****************************************************************************80
%
%% CLEBSCH_GORDAN_VALUES_TEST demonstrates the use of CLEBSCH_GORDAN_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CLEBSCH_GORDAN_VALUES_TEST:\n' );
  fprintf ( 1, '  CLEBSCH_GORDAN_VALUES returns values of \n' );
  fprintf ( 1, '  the Clebsch-Gordan coefficient.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    '      J1      J2      J3      M1      M2      M3        CG\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, j1, j2, j3, m1, m2, m3, fx ] = clebsch_gordan_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %6.2f  %6.2f  %6.2f  %6.2f  %6.2f  %6.2f  %24.16f\n', ...
      j1, j2, j3, m1, m2, m3, fx );

  end

  return
end

author: johannesgerer


text:
function inside = p08_inside ( m, n, point )

%*****************************************************************************80
%
%% P08_INSIDE reports if a point is inside the region in problem 08.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real POINT(M,N), the coordinates
%    of the points.
%
%    Output, logical INSIDE(N), is TRUE if the point is in the region.
%
  nv = 6;

  c1 = [ 0.0, 0.0 ];
  c2 = [ 0.6, 0.0 ];

  r1 = 1.0;
  r2 = 0.1;
  theta1 = - pi / 12.0;
  theta2 =   pi / 12.0;
%
%  These coordinates are for a square with one corner at (0.9, 0.0),
%  and sides of length 0.2.  It's somewhat bigger than the "bite"
%  it takes out of the piece of pie.
%
  x_poly = [ ...
    0.9000000,  0.0000000,
    1.0414213, -0.1414213,
    1.1828427,  0.0000000,
    1.0414213,  0.1414213 ]';

  inside(1:n) = 0;

  for j = 1 : n
%
%  Is the point inside the circular sector?
%
    if ( circle_sector_contains_point_2d ( r1, c1, theta1, theta2, ...
       point(1:2,j) ) )
%
%  Is the point NOT inside the box?
%
      if ( ~polygon_contains_point_2d ( 4, x_poly, point(1:2,j) ) )
%
%  Is the point NOT inside the circle?
%
        if ( ~circle_imp_contains_point_2d ( r2, c2, point(1:2,j) ) )
          inside(j) = 1;
        end
      end
    end

  end

  return
end


text:
function [ ] = FlipDimRunTimeAnalysis(  )

numIterations   = 100;
maxDim          = 1000;

runTime = 0;

for ii = 1:maxDim
    vX = randn([ii, 1]);
    hRunTime = tic();
    for jj = 1:numIterations
        vY = flip(vX);
    end
    runTime = runTime + toc(hRunTime);
end

disp(['Run Time - ', num2str(runTime)]);

runTime = 0;

for ii = 1:maxDim
    vX = randn([ii, 1]);
    hRunTime = tic();
    for jj = 1:numIterations
        vY = flipdim(vX, 1);
    end
    runTime = runTime + toc(hRunTime);
end

disp(['Run Time - ', num2str(runTime)]);

end

author: RoyiAvital

repo: StackExchangeCodes


text:
%RGB2HSV Dataset overload
%
%  B = RGB2HSV(A)

author: marianux

repo: ecg-kit

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


text:
function OUT = L(X,nlag)
% =======================================================================
% Creates a matrix (or vector) of lagged values (the first obs are NaN)
% =======================================================================
% OUT = L(X,nlag)
% -----------------------------------------------------------------------
% INPUT
%	- X: input matrix (nobs X nvars)
%	- nlag: order of lags (-1 is lag ---- +1 is lead)
% -----------------------------------------------------------------------
% OUTPUT
%	- OUT: lagged matrix (nobs X nvars), the first obs are NaN
% -----------------------------------------------------------------------
% EXAMPLE
%	x = 1:10';
%   L1x = L(x,-1);
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% March 2012. Updated November 2020
% -----------------------------------------------------------------------

init = NaN;
switch(nargin)
    case 1
        error('Missing input: number of lags');
    case 2
       if nlag > 0 % this is lead
           zt = ones(nlag,cols(X))*init;
           OUT = [trimr(X,nlag,0); zt];
       elseif nlag < 0 % this is lag
           zt = ones(abs(nlag),cols(X))*init;
           OUT = [ zt; trimr(X,0,abs(nlag))];
       else
           OUT = X;
       end
    otherwise
        error('Too many inputs for function L');
end

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
% evaluate simulated data - ICA and NMF
% V ~ WH
% V -> N_pix x N_t
% W -> N_pix x N_comp - xth pixel of the ith components
% H -> N_copm x N_t - contribution of the i-th component in the time t

function separNMFICA(sep0, offset0, path_data, path_res, prename, niter, avg_num, savethis, sep_how)

if ~exist('initval', 'var')
    initval = 0;
end

if ~exist('sep_how', 'var')
    sep_how = 'in';
end
ncomp = 3; %number of components to be separated + background
randh=rand(2,1,4); %random coeeficients for hinit

for rr = 1: length(offset0)
    for ll=1 : length(sep0)
        for qq=1 : length(avg_num)
            namedir = [prename num2str(100*sep0(ll)) 'offset_' num2str(offset0(rr))];
            cd ([path_data namedir])
            for mm=1:niter
                %reads the first...
                namefile = [namedir '-iter_' num2str(mm)];
                load ([namefile '.mat'])
                %cat the folowing,,,
                p.catitervec=[1];
                [dpixc, dveccr, blinkmat, p] = catsimul(namedir, p.catitervec);
                if sum(sep_how == 'i')>0 %ICA
                    [icasig{mm}, A{mm}, W{mm}] = fastica (dveccr, 'numOfIC', ncomp, 'g', 'tanh');
                    icapixICA{mm} = reshape(A{mm},p.nx, p.ny, ncomp);
                end

                if sum(sep_how == 'n')>0 %NMF
                    % background estimation:
                    %[out, bg(mm), bg_im]=backgroundoffset(dpixc);
                    [out, bg(mm), bg_im]=backgroundoffset(dpixc, 'no', 5, 20, 8); %empirical values...

                    dvec_bg = ones(p.nx*p.ny, 1);
                    %                     dvec_bg = p.offset*ones(1, p.nx*p.ny); %changed for
                    %                     offset 10...

                    dvec_ind = squeeze(reshape(double(array2im(dpixc_ind)), p.nx*p.ny, 1, 2)); % vectors of resized images
                    %                     sum_dvec_ind = sum(dvec_ind, 1);
                    %                     dvec_ind = dvec_ind./repmat(sum_dvec_ind, p.nx*p.ny,1); %normlaized

                    %                     winit = [f*dvec_ind'; dvec_bg];       %original 'true' points + background
                    % % %                     winittmp = [dvec_ind, dvec_bg];
                    winittmp = [rand(size(dvec_ind)), dvec_bg];
                    sumw = sum(winittmp,1);
                    winit = winittmp./repmat(sumw, p.nx*p.ny, 1); %normalized to 1
                    f = mean(dveccr(:)-bg(mm))/mean(mean(winit(:, 1:2))); %ration of the data/psf
                    %                     winit = [rand(ncomp, p.nx*p.ny); dvec_bg];
                    blinkmatrand = rand(ncomp-1, p.Nt); %uniform random;
                    hinit = [f*blinkmatrand; bg(mm)*sumw(ncomp)*ones(1, p.Nt)];             %random weights will be assigned to firts two and bg fixed
                    % % %                     hinit = [f*blinkmat./repmat(mean(blinkmat,2),1,size(blinkmat,2)); bg(mm)*sumw(ncomp)*ones(1, p.Nt)];             %random weights will be assigned to firts two and bg fixed

                    %                    [w{mm},h{mm}, wtrace{mm},wtrace{mm}]=nmf_test(double(dveccr'),ncomp+1,1,hinit,winit, [3], [3]);
                    %                     [w{mm},h{mm}, wtrace,htrace,ddiv{mm}]=nmf_testconvD(double(dveccr'),ncomp+1,1,hinit,winit, [3], [3]);

                    dvec_bg = ones(p.nx*p.ny, 1);
                    winittmp = [dvec_ind, dvec_bg];
                    sumw = sum(winittmp,1);
                    winit = winittmp./repmat(sumw, p.nx*p.ny, 1); %normalized to 1
                    f = mean(dveccr(:)-bg(mm))/mean(mean(winit(:, 1:2))); %ration of the data/psf

                    p.Nt=1;
                    p.meanblinkmat=mean(mean(blinkmat,2));
                    p.htrue=[0.3; 0.7];
                    htruef=p.meanblinkmat*p.htrue;

                    %hinit = [p.meanblinkmat*rand(2,1); p.offset*sumw(ncomp)*ones(1, p.Nt)];
                    %hinit = [p.meanblinkmat*rand(2,1); p.offset*sumw(ncomp)*ones(1, p.Nt)];
                    for kk=1:4
                        hinit = [p.meanblinkmat*randh(:,:,kk); p.offset*sumw(ncomp)*ones(1, p.Nt)];
                        %                     hinit = [htruef; p.offset*sumw(ncomp)*ones(1, p.Nt)];
                        dveccr=dvec_ind*htruef+p.offset;
                        dpixcd = dip_image(reshape(dveccr,p.nx,p.ny));
                        dpixcdn_tmp = newim(dpixcd);
                        for oo=1:avg_num(qq)
                            dpixcdn_tmp = dpixcdn_tmp + noise(dpixcd,'poisson');
                        end
                        if avg_num(qq)==0
                            dpixcdn_tmp=dpixcd; %noise free
                        end
                        dpixcdn=dpixcdn_tmp/max(avg_num(qq),1); %average of (avg_num(qq))X realization...
                        dveccr=reshape(double(dpixcdn),p.nx*p.ny,1);
                        %[c,w{mm},h{mm}, X1,X2, dhr, minXr, miXvalr, mhdr, htrace, p]
                        [c,w{mm},h{mm}, X1,X2, dh, minX, miXval, mhd, htr,htrace, p]=nmf_S48(double(dveccr),ncomp,1,winit,hinit, [3], [3],p);
                        res{kk}=struct('c',c,'w',w,'h',h,'X1',X1,'X2',X2,'dh',dh,'minX',minX,'miXval',miXval, 'mhd',mhd, 'htr',htr, 'htrace',htrace,'p',p);
                    end

                    %hinit = [f*blinkmat./repmat(mean(blinkmat,2),1,size(blinkmat,2)); bg(mm)*sumw(ncomp)*ones(1, p.Nt)];             %random weights will be assigned to firts two and bg fixed
                    %[c,w{mm},h{mm}, X1,X2, dht, minXt, miXvalt, mhdt, htrace,p]=nmf_S41(double(dveccr),ncomp,1,winit,hinit, [3], [3],p);
                    qqq=[];
                    icapixNMF{mm} = reshape(w{mm},p.nx,p.ny,ncomp);
                end

            end
            p.avg_num = avg_num(qq);
            p.path_data = path_data;
            p.path_res = path_res;

            if savethis == 1
                fprintf('saving data \n');
                if ~(strcmp(path_res, path_data)) %not identical
                    mkdir ([path_res namedir]);
                    cd ([path_res namedir]);
                end
                qqq=[];
                %             save (p.namedir)
                %save ([namedir '_separMulti'],save ([namedir '_separMulti'],'p', 'X1','X2','dhr','minXr', 'miXvalr','dht'))
                save ([namedir '_DdviMap_avg' num2str(avg_num(qq))],'res')
                writedata([],[],p,[namedir '_avg' num2str(avg_num(qq)) '_param'])
            end
        end
    end
end

fprintf('\n')

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [ n_data, x, fx ] = sin_values ( n_data )

%*****************************************************************************80
%
%% SIN_VALUES returns some values of the sine function.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      Sin[x]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Cambridge University Press, 1999,
%    ISBN: 0-521-64314-7,
%    LC: QA76.95.W65.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 13;

  fx_vec = [ ...
     0.00000000000000000000, ...
     0.25881904510252076235, ...
     0.47942553860420300027, ...
     0.50000000000000000000, ...
     0.70710678118654752440, ...
     0.84147098480789650665, ...
     0.86602540378443864676, ...
     1.00000000000000000000, ...
     0.90929742682568169540, ...
     0.14112000805986722210, ...
     0.00000000000000000000, ...
    -0.75680249530792825137, ...
    -0.95892427466313846889 ];

  x_vec = [ ...
    0.0000000000000000000, ...
    0.26179938779914943654, ...
    0.50000000000000000000, ...
    0.52359877559829887308, ...
    0.78539816339744830962, ...
    1.0000000000000000000, ...
    1.0471975511965977462, ...
    1.5707963267948966192, ...
    2.0000000000000000000, ...
    3.0000000000000000000, ...
    3.1415926535897932385, ...
    4.0000000000000000000, ...
    5.0000000000000000000  ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end

repo: jburkardt-m


text:
function test_example_simulateddata_beamformer

% MEM 4gb
% WALLTIME 00:10:00

%
%% Compute forward simulated data and apply a beamformer scan
%
% This example script shows you how to create some simulated channel-level MEG data with a single dipole at a specified location in the head. Subsequently it does a beamformer source reconstruction to localize that source.
%
% create an array with some magnetometers at 12cm distance from the origin
[X, Y, Z] = sphere(10);
pos = unique([X(:) Y(:) Z(:)], 'rows');
pos = pos(pos(:,3)>=0,:);
grad = [];
grad.coilpos = 12*pos;
grad.coilori = pos; % in the outward direction
% grad.tra = eye(length(pos)); % each coils contributes exactly to one channel
for i=1:length(pos)
  grad.label{i} = sprintf('chan%03d', i);
end

% create a spherical volume conductor with 10cm radius
vol.r = 10;
vol.o = [0 0 0];

% note that beamformer scanning will be done with a 1cm grid, so you should
% not put the dipole on a position that will not be covered by a grid
% location later
cfg = [];
cfg.headmodel = vol;
cfg.grad = grad;
cfg.dip.pos = [0 0 4];    % you can vary the location, here the dipole is along the z-axis
cfg.dip.mom = [1 0 0]';   % the dipole points along the x-axis
cfg.relnoise = 10;
cfg.ntrials = 20;
data = ft_dipolesimulation(cfg);

% compute the data covariance matrix, which will capture the activity of
% the simulated dipole
cfg = [];
cfg.covariance = 'yes';
timelock = ft_timelockanalysis(cfg, data);

% do the beamformer source reconstuction on a 1 cm grid
cfg = [];
cfg.headmodel = vol;
cfg.grad = grad;
cfg.resolution = 1;
cfg.method = 'lcmv';
cfg.lcmv.projectnoise = 'yes'; % needed for neural activity index
source = ft_sourceanalysis(cfg, timelock);

% compute the neural activity index, i.e. projected power divided by
% projected noise
cfg = [];
cfg.powmethod = 'none'; % keep the power as estimated from the data covariance, i.e. the induced power
source_nai = ft_sourcedescriptives(cfg, source);

cfg = [];
cfg.method = 'ortho';
cfg.funparameter = 'nai';
cfg.funcolorlim = [1.4 1.5];  % the voxel in the center of the volume conductor messes up the autoscaling
ft_sourceplot(cfg, source_nai);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%PLOTSOM Plot the Self-Organizing Map in 2D
%
%    PLOTSOM(W)
%
% Plot the Self-Organizing Map W, trained by som.m. This is only
% possible if the map is 2D.
%
% SEE ALSO (PRTools Guide)
% SOM


text:
function vifp=vifp_mscale(ref,dist)

% -----------COPYRIGHT NOTICE STARTS WITH THIS LINE------------
% Copyright (c) 2005 The University of Texas at Austin
% All rights reserved.
%
% Permission is hereby granted, without written agreement and without license or royalty fees, to use, copy,
% modify, and distribute this code (the source files) and its documentation for
% any purpose, provided that the copyright notice in its entirety appear in all copies of this code, and the
% original source of this code, Laboratory for Image and Video Engineering (LIVE, http://live.ece.utexas.edu)
% at the University of Texas at Austin (UT Austin,
% http://www.utexas.edu), is acknowledged in any publication that reports research using this code. The research
% is to be cited in the bibliography as:
%
% H. R. Sheikh and A. C. Bovik, "Image Information and Visual Quality", IEEE Transactions on
% Image Processing, (to appear).
%
% IN NO EVENT SHALL THE UNIVERSITY OF TEXAS AT AUSTIN BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL,
% OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OF THIS DATABASE AND ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF TEXAS
% AT AUSTIN HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
% THE UNIVERSITY OF TEXAS AT AUSTIN SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE DATABASE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS,
% AND THE UNIVERSITY OF TEXAS AT AUSTIN HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
%
% -----------COPYRIGHT NOTICE ENDS WITH THIS LINE------------
%
% This software release consists of a MULTISCALE PIXEL DOMAIN, SCALAR GSM implementation of the algorithm described in the paper:
%
% H. R. Sheikh and A. C. Bovik, "Image Information and Visual Quality"., IEEE Transactions on Image Processing, (to appear).
% Download manuscript draft from http://live.ece.utexas.edu in the Publications link.
%
% THE PIXEL DOMAIN ALGORITHM IS NOT DESCRIBED IN THE PAPER. THIS IS A COMPUTATIONALLY SIMPLER
% DERIVATIVE OF THE ALGORITHM PRESENTED IN THE PAPER
%
% Input : (1) img1: The reference image as a matrix
%         (2) img2: The distorted image (order is important)
%
% Output: (1) VIF the visual information fidelity measure between the two images
%
% Default Usage:
%    Given 2 test images img1 and img2, whose dynamic range is 0-255
%
%    vif = vifvec(img1, img2);
%
% Advanced Usage:
%    Users may want to modify the parameters in the code.
%    (1) Modify sigma_nsq to find tune for your image dataset.
% Email comments and bug reports to hamid.sheikh@ieee.org

sigma_nsq=2;

num=0;
den=0;
for scale=1:4

    N=2^(4-scale+1)+1;
    win=fspecial('gaussian',N,N/5);

    if (scale >1)
        ref=filter2(win,ref,'valid');
        dist=filter2(win,dist,'valid');
        ref=ref(1:2:end,1:2:end);
        dist=dist(1:2:end,1:2:end);
    end

    mu1   = filter2(win, ref, 'valid');
    mu2   = filter2(win, dist, 'valid');
    mu1_sq = mu1.*mu1;
    mu2_sq = mu2.*mu2;
    mu1_mu2 = mu1.*mu2;
    sigma1_sq = filter2(win, ref.*ref, 'valid') - mu1_sq;
    sigma2_sq = filter2(win, dist.*dist, 'valid') - mu2_sq;
    sigma12 = filter2(win, ref.*dist, 'valid') - mu1_mu2;

    sigma1_sq(sigma1_sq<0)=0;
    sigma2_sq(sigma2_sq<0)=0;

    g=sigma12./(sigma1_sq+1e-10);
    sv_sq=sigma2_sq-g.*sigma12;

    g(sigma1_sq<1e-10)=0;
    sv_sq(sigma1_sq<1e-10)=sigma2_sq(sigma1_sq<1e-10);
    sigma1_sq(sigma1_sq<1e-10)=0;

    g(sigma2_sq<1e-10)=0;
    sv_sq(sigma2_sq<1e-10)=0;

    sv_sq(g<0)=sigma2_sq(g<0);
    g(g<0)=0;
    sv_sq(sv_sq<=1e-10)=1e-10;

     num=num+sum(sum(log10(1+g.^2.*sigma1_sq./(sv_sq+sigma_nsq))));
     den=den+sum(sum(log10(1+sigma1_sq./sigma_nsq)));

end
vifp=num/den;

author: Linfeng-Tang

repo: Image-Fusion


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MRiLab auto generated file: DO NOT EDIT!     %
% Generated by MRiLab "DoWriteXML2m" Generator %
% MRiLab Version 1.3                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [rfAmp,rfPhase,rfFreq,rfCoil,GzAmp,GyAmp,GxAmp,ADC,Ext,uts,ts,flags]=PSD_SPGR3DME
global VCtl
global VVar
CV1=2e-3;
CV10=0;
CV11=0;
CV12=0;
CV13=0;
CV14=0;
CV2=3e-3;
CV3=abs(rem(VVar.TRCount,2)*2-1);
CV4=0e-3;
CV5=0;
CV6=0;
CV7=0;
CV8=0;
CV9=0;
rfAmpAll=[];
rfPhaseAll=[];
rfFreqAll=[];
rfCoilAll=[];
GzAmpAll=[];
GyAmpAll=[];
GxAmpAll=[];
ADCAll=[];
ExtAll=[];
rfTimeAll=[];
GzTimeAll=[];
GyTimeAll=[];
GxTimeAll=[];
ADCTimeAll=[];
ExtTimeAll=[];
SEtAll=[];
uts=[];
ts=[];
flags=[];
if VCtl.PlotSeq == 1
rfAmp=[];
rfPhase=[];
rfFreq=[];
rfCoil=[];
GzAmp=[];
GyAmp=[];
GxAmp=[];
ADC=[];
Ext=[];
Freq=1;
Notes='regular TR section';
AttributeOpt={'on','off'};
Switch=AttributeOpt{1};
TREnd=Inf;
TRStart=1;
tE=VCtl.TR;
tS=0;
if VVar.TRCount<TRStart | VVar.TRCount>TREnd | mod(VVar.TRCount-TRStart,Freq)~=0 | strcmp(Switch,'off')
% do nothing
else
ts = [ts tS tE];
end
ts = [0 max(ts)-min(ts)];
return;
end
%==============Pulses 1==============
rfAmp=[];
rfPhase=[];
rfFreq=[];
rfCoil=[];
GzAmp=[];
GyAmp=[];
GxAmp=[];
ADC=[];
Ext=[];
rfTime=[];
GzTime=[];
GyTime=[];
GxTime=[];
ADCTime=[];
ExtTime=[];
Freq=1;
Notes='regular TR section';
AttributeOpt={'on','off'};
Switch=AttributeOpt{1};
TREnd=Inf;
TRStart=1;
tE=VCtl.TR;
tS=0;
if isempty(tS) | isempty(tE) | (tS>=tE)
error('SE setting is incorrect for Pulses 1!');
end
if VVar.TRCount<TRStart | VVar.TRCount>TREnd | mod(VVar.TRCount-TRStart,Freq)~=0 | strcmp(Switch,'off')
% do nothing
else
%--------------------
AttributeOpt={'on','off'};
p.AnchorTE=AttributeOpt{1};
AttributeOpt={'Non','Hamming','Hanning'};
p.Apod=AttributeOpt{2};
p.CoilID=1;
p.DupSpacing=0;
p.Duplicates=1;
p.FA=VCtl.FlipAng;
p.Notes='sinc rf pulse';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.TBP=4;
p.dt=20e-6;
p.rfFreq=0;
p.rfPhase=0;
p.tEnd=CV4+0.5e-3;
p.tStart=CV4;
if strcmp(p.Switch,'on')
if strcmp(p.AnchorTE,'on')
switch VCtl.TEAnchor
case 'Start'
VCtl.TEAnchorTime=p.tStart;
case 'Middle'
VCtl.TEAnchorTime=(p.tStart+p.tEnd)/2;
case 'End'
VCtl.TEAnchorTime=p.tEnd;
end
end
[rfAmp1,rfPhase1,rfFreq1,rfCoil1,rfTime1]=rfSinc(p);
if strcmp(VCtl.MultiTransmit,'off')
if VCtl.MasterTxCoil==rfCoil1(1)
rfAmp=[rfAmp rfAmp1];
rfPhase=[rfPhase rfPhase1];
rfFreq=[rfFreq rfFreq1];
rfCoil=[rfCoil rfCoil1];
rfTime=[rfTime rfTime1];
end
else
rfAmp=[rfAmp rfAmp1];
rfPhase=[rfPhase rfPhase1];
rfFreq=[rfFreq rfFreq1];
rfCoil=[rfCoil rfCoil1];
rfTime=[rfTime rfTime1];
end
end
p=[];
%--------------------
AttributeOpt={'on','off'};
p.AnchorTE=AttributeOpt{2};
AttributeOpt={'Non','Hamming','Hanning'};
p.Apod=AttributeOpt{2};
p.CoilID=1;
p.DupSpacing=0;
p.Duplicates=1;
p.FA=0;
p.Notes='sinc rf pulse';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.TBP=4;
p.dt=0.05e-3;
p.rfFreq=0;
p.rfPhase=0;
p.tEnd=VCtl.TR;
p.tStart=CV4+0.6e-3;
if strcmp(p.Switch,'on')
if strcmp(p.AnchorTE,'on')
switch VCtl.TEAnchor
case 'Start'
VCtl.TEAnchorTime=p.tStart;
case 'Middle'
VCtl.TEAnchorTime=(p.tStart+p.tEnd)/2;
case 'End'
VCtl.TEAnchorTime=p.tEnd;
end
end
[rfAmp2,rfPhase2,rfFreq2,rfCoil2,rfTime2]=rfSinc(p);
if strcmp(VCtl.MultiTransmit,'off')
if VCtl.MasterTxCoil==rfCoil2(1)
rfAmp=[rfAmp rfAmp2];
rfPhase=[rfPhase rfPhase2];
rfFreq=[rfFreq rfFreq2];
rfCoil=[rfCoil rfCoil2];
rfTime=[rfTime rfTime2];
end
else
rfAmp=[rfAmp rfAmp2];
rfPhase=[rfPhase rfPhase2];
rfFreq=[rfFreq rfFreq2];
rfCoil=[rfCoil rfCoil2];
rfTime=[rfTime rfTime2];
end
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Gz1Sign=1;
p.Gz2Sign=0;
p.Notes='cartesian phase';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.t1End=CV4+CV2;
p.t1Start=CV4+CV1;
p.t2End=VCtl.TE+CV2;
p.t2Start=VCtl.TE+CV1;
p.tRamp=100e-6;
if strcmp(p.Switch,'on')
[GzAmp1,GzTime1]=GzCartesian(p);
GzAmp=[GzAmp GzAmp1];
GzTime=[GzTime GzTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Gy1Sign=1;
p.Gy2Sign=0;
p.Notes='cartesian phase';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.t1End=CV4+CV2;
p.t1Start=CV4+CV1;
p.t2End=VCtl.TE+CV2;
p.t2Start=VCtl.TE+CV1;
p.tRamp=100e-6;
if strcmp(p.Switch,'on')
[GyAmp1,GyTime1]=GyCartesian(p);
GyAmp=[GyAmp GyAmp1];
GyTime=[GyTime GyTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Gx1Sign=-1;
p.Gx2Sign=1;
p.Gx3Sign=0;
p.Notes='cartesian frequency';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.t1Start=CV4+CV1;
p.t2Middle=VCtl.TE;
p.t3Start=VCtl.TE+CV1;
p.tRamp=100e-6;
if strcmp(p.Switch,'on')
[GxAmp1,GxTime1]=GxCartesian(p);
GxAmp=[GxAmp GxAmp1];
GxTime=[GxTime GxTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Notes='cartesian readout';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tMiddle=VCtl.TE;
if strcmp(p.Switch,'on')
[ADC1,ADCTime1]=ADCCartesian(p);
ADC=[ADC ADC1];
ADCTime=[ADCTime ADCTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Ext=5;
p.Notes='calculate remaining scan time';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tStart=0;
if strcmp(p.Switch,'on')
[Ext1,ExtTime1]=ExtBit(p);
Ext=[Ext Ext1];
ExtTime=[ExtTime ExtTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Ext=1;
p.Notes='reset K space location';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tStart=CV4+0.5e-3;
if strcmp(p.Switch,'on')
[Ext2,ExtTime2]=ExtBit(p);
Ext=[Ext Ext2];
ExtTime=[ExtTime ExtTime2];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Ext=6;
p.Notes='dephase Mxy';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tStart=VCtl.TR*(99/100);
if strcmp(p.Switch,'on')
[Ext3,ExtTime3]=ExtBit(p);
Ext=[Ext Ext3];
ExtTime=[ExtTime ExtTime3];
end
p=[];
%--------------------
SEt=[tS tE];
rfAmp(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
rfPhase(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
rfFreq(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
rfCoil(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
GzAmp(GzTime<0 | GzTime>SEt(2)-SEt(1)) = [];
GyAmp(GyTime<0 | GyTime>SEt(2)-SEt(1)) = [];
GxAmp(GxTime<0 | GxTime>SEt(2)-SEt(1)) = [];
ADC(ADCTime<0 | ADCTime>SEt(2)-SEt(1)) = [];
Ext(ExtTime<0 | ExtTime>SEt(2)-SEt(1)) = [];
rfTime(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
GzTime(GzTime<0 | GzTime>SEt(2)-SEt(1)) = [];
GyTime(GyTime<0 | GyTime>SEt(2)-SEt(1)) = [];
GxTime(GxTime<0 | GxTime>SEt(2)-SEt(1)) = [];
ADCTime(ADCTime<0 | ADCTime>SEt(2)-SEt(1)) = [];
ExtTime(ExtTime<0 | ExtTime>SEt(2)-SEt(1)) = [];
rfAmp(abs(rfAmp)<eps) = 0;
rfTime = rfTime + SEt(1);
GzTime = GzTime + SEt(1);
GyTime = GyTime + SEt(1);
GxTime = GxTime + SEt(1);
ADCTime = ADCTime + SEt(1);
ExtTime = ExtTime + SEt(1);
rfAmpAll=[rfAmpAll rfAmp];
rfPhaseAll=[rfPhaseAll rfPhase];
rfFreqAll=[rfFreqAll rfFreq];
rfCoilAll=[rfCoilAll rfCoil];
GzAmpAll=[GzAmpAll GzAmp];
GyAmpAll=[GyAmpAll GyAmp];
GxAmpAll=[GxAmpAll GxAmp];
ADCAll=[ADCAll ADC];
ExtAll=[ExtAll Ext];
rfTimeAll=[rfTimeAll rfTime];
GzTimeAll=[GzTimeAll GzTime];
GyTimeAll=[GyTimeAll GyTime];
GxTimeAll=[GxTimeAll GxTime];
ADCTimeAll=[ADCTimeAll ADCTime];
ExtTimeAll=[ExtTimeAll ExtTime];
SEtAll=[SEtAll SEt];
end
%====================================
if isempty(rfTimeAll)
error('rf sequence line can not be empty! Master Tx coil element must be used.');
end
if isempty(GzTimeAll)
error('GzSS sequence line can not be empty!');
end
if isempty(GyTimeAll)
error('GyPE sequence line can not be empty!');
end
if isempty(GxTimeAll)
error('GxR sequence line can not be empty!');
end
if isempty(ADCTimeAll)
error('ADC sequence line can not be empty!');
end
if isempty(ExtTimeAll)
error('Ext sequence line can not be empty!');
end
SEflag=repmat([0 0 0 0 0 0]',[1 2]);
rfflag=repmat([1 0 0 0 0 0]',[1 max(size(rfTimeAll))]);
Gzflag=repmat([0 1 0 0 0 0]',[1 max(size(GzTimeAll))]);
Gyflag=repmat([0 0 1 0 0 0]',[1 max(size(GyTimeAll))]);
Gxflag=repmat([0 0 0 1 0 0]',[1 max(size(GxTimeAll))]);
ADCflag=repmat([0 0 0 0 1 0]',[1 max(size(ADCTimeAll))]);
Extflag=repmat([0 0 0 0 0 1]',[1 max(size(ExtTimeAll))]);
ts=[[min(SEtAll) max(SEtAll)] rfTimeAll GzTimeAll GyTimeAll GxTimeAll ADCTimeAll ExtTimeAll]-min(SEtAll);
flags=[SEflag rfflag Gzflag Gyflag Gxflag ADCflag Extflag];
[ts,ind]=sort(ts);
uts=unique(ts);
flags=flags(:,ind);
[rfTime,ind]=sort(rfTimeAll-min(SEtAll));
rfAmp=rfAmpAll(:,ind);
rfPhase=rfPhaseAll(:,ind);
rfFreq=rfFreqAll(:,ind);
rfCoil=rfCoilAll(:,ind);
[GzTime,ind]=sort(GzTimeAll-min(SEtAll));
GzAmp=GzAmpAll(:,ind);
[GyTime,ind]=sort(GyTimeAll-min(SEtAll));
GyAmp=GyAmpAll(:,ind);
[GxTime,ind]=sort(GxTimeAll-min(SEtAll));
GxAmp=GxAmpAll(:,ind);
[ADCTime,ind]=sort(ADCTimeAll-min(SEtAll));
ADC=ADCAll(:,ind);
[ExtTime,ind]=sort(ExtTimeAll-min(SEtAll));
Ext=ExtAll(:,ind);
rfAmp(1) = 0;
rfPhase(1) = 0;
rfFreq(1) = 0;
GzAmp(1) = 0;
GyAmp(1) = 0;
GxAmp(1) = 0;
ADC(1) = 0;
Ext(1) = 0;
rfAmp(end) = 0;
rfPhase(end) = 0;
rfFreq(end) = 0;
GzAmp(end) = 0;
GyAmp(end) = 0;
GxAmp(end) = 0;
ADC(end) = 0;
Ext(end) = 0;
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [mask,numClusters,XYZ] = clusterSizeMask(sizeThresh,height_mask)
% :Usage:
% ::
%
%     function [mask,numClusters,XYZ] = clusterSizeMask(sizeThresh,height_mask)
%
% ..
%    Tor Wager, 10/27/01
% ..

mask = []; numClusters = 0;, XYZ = [];

% Get point list of activated voxels
% ------------------------------------------------------
voxels = mask2voxel(height_mask);	% returns n x 3
voxels = voxels';					% put xyz in a single column

if isempty(voxels)
	disp('No voxels meet height threshold.')
	mask = zeros(size(height_mask));
	return
end

% Get cluster indices of voxels
% ------------------------------------------------------
[cl_index] = spm_clusters(voxels);

% Find index of voxels of sufficient size
% ------------------------------------------------------
if ~isempty(sizeThresh) & sizeThresh > 0
    for i = 1:max(cl_index)
	    a(cl_index == i) = sum(cl_index == i);
    end
else
    sizeThresh = 0;
    a = ones(size(cl_index));
end

which_vox = (a >= sizeThresh);
numClusters = sum(length(unique(cl_index(find(a >= sizeThresh)))));

if numClusters == 0
	disp('No clusters meet extent threshold.')
	mask = zeros(size(height_mask));
    XYZ = [];
	return
end

voxels = voxels(:,which_vox);

% Voxel point list output - 3 vector
% ------------------------------------------------------
XYZ = voxels;

% Convert back to mask for mask output
% ------------------------------------------------------
voxels = voxels';					% convert back to row vectors
mask = voxel2mask(voxels,size(height_mask));

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [ Psi ] = buildPsi2( options, psi, configuration )
%BUILDPSI functions creates vector \Psi(I, S) for given configuration
%   Detailed explanation goes here
%
% INPUT
%   options         ...     structure holds info about components (size, count) and S (statistics of component position)
%   psi             ...     cell Mx1 contains matrices \Psi^q(I, S_i)
%   configuration   ...     matrix 2xM [x0, ..., x_{M-1}; y_0, ..., y_{M-1}] with positions of components
%
% OUTPUT
%   Psi     ...     vector [\Psi_0^q(I, S_0); ... ; \Psi_{M-1}^q(I, S_{M-1}); g(S_0xS_1); ... ; g(S_0xS_{M-1})];
%
% 10-08-10 Michal Uricar
% 11-07-11 Michal Uricar, corners dataset modificatoin

    indices = zeros(1, options.M);

    Psi = [];

    for i = 1 : options.M
        psi_si = psi{i};
        P = configuration(:, i) - [options.S(1, i) - 1; options.S(2, i) - 1];
        siz = [options.S(4, i) - options.S(2, i) + 1; options.S(3, i) - options.S(1, i) + 1];
        indices(i) = (P(1) - 1)*siz(1) + P(2);
        Psi = [Psi; psi_si(:, indices(i))];
    end;

    % g
    %for i = 2 : options.M - 2
    for i = 2 : options.M - 3
        dxdy = configuration(:, i) - configuration(:, 1);
        g = [dxdy; dxdy.^2];
        Psi = [Psi; g];
    end;
    % g5
    dxdy = configuration(:, 6) - configuration(:, 2);
    g = [dxdy; dxdy.^2];
    Psi = [Psi; g];
    % g6
    dxdy = configuration(:, 7) - configuration(:, 3);
    g = [dxdy; dxdy.^2];
    Psi = [Psi; g];
    % g7
    dxdy = configuration(:, 8) - configuration(:, 1);
    g = [dxdy; dxdy.^2];
    Psi = [Psi; g];

end

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
function res = fsample(this, value)
% Method for getting and setting the sampling rate
% FORMAT res = fsample(this)
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: fsample.m 3200 2009-06-12 17:29:40Z vladimir $

if nargin == 1
    res = this.Fsample;
else
    this.Fsample = value;
    res = this;
end


text:
function [modelPlanePyr, modelRegPyr] = sc_planar_structure_pyramid(scaleImgPyr, modelPlane, modelReg)

%
% SC_PLANAR_STRUCTURE_PYRAMID:
%
% Pre-compute the models of planes and regularity according to the image
% pyramid, rescale parameters according to the image sizes.
%
% Output:
%   - modelPlanePyr: plane model
%   - modelRegPyr: regularity model

numLevel = length(scaleImgPyr);

modelPlanePyr = cell(numLevel, 1);
modelRegPyr = cell(numLevel, 1);

planePostProb = modelPlane.postProb;

for iLvl = 1: numLevel
    % === Update plane model ===
    scaleImgCur = scaleImgPyr{iLvl}.imgScale;
    sizeImgCur = scaleImgPyr{iLvl}.imgSize;

    % Update posterior plane probability
    planePostProbCur = imresize(planePostProb, sizeImgCur, 'bicubic');
    planePostProbCurSum = sum(planePostProbCur, 3);
    planePostProbCur = bsxfun(@rdivide, planePostProbCur, planePostProbCurSum);
    modelPlanePyr{iLvl}.numPlane = modelPlane.numPlane;
    modelPlanePyr{iLvl}.planeProb = planePostProbCur;
    modelPlanePyr{iLvl}.mLogLPlaneProb = -log(planePostProbCur);

    % Update rectification matrix and rotation parameters
    for iPlane = 1: modelPlane.numPlane
        H = eye(3);
        vLine = modelPlane.plane{iPlane}.vLine;
        vLine(1:2) = vLine(1:2)/scaleImgCur;
        H(3,:) = vLine;
        modelPlanePyr{iLvl}.rectMat{iPlane} = H;
        modelPlanePyr{iLvl}.rotPar{iPlane} = modelPlane.plane{iPlane}.rotPar;

        % Construct the rotation matrices
        for iTheta = 1:2
            t = modelPlanePyr{iLvl}.rotPar{iPlane}(iTheta);
            Hr = eye(3);
            Hr(1,1) = cos(t);   Hr(1,2) = - sin(t);
            Hr(2,1) = sin(t);   Hr(2,2) =   cos(t);
            modelPlanePyr{iLvl}.rotMat{iPlane, iTheta} = Hr;
        end
    end

    % === Update reguarlity model ===
    for iPlane = 1: modelPlane.numPlane
        modelRegPyr{iLvl}.dispVec{iPlane} = scaleImgCur*modelReg.plane{iPlane}.dispVec;
        modelRegPyr{iLvl}.numDispVec(iPlane) = modelReg.plane{iPlane}.numDispVec;
    end
end

end

repo: StructCompletion


text:
% CONVERT SPHERICAL TO CARTESIAN
function [p] = GetCartesianFromSpherical(r,psi,theta)
% This function calculates the relative position from a spherical
% coordinate system that assumes angles are measured from the
% agents heading vector Xref.
% INPUTS:
% range     - The objects radial seperation (m)
% azimuth   - The objects relative azimuth angle (rad)
% elevation - The objects relative elevation (rad)
% OUTPUTS:
% cartesianPosition  - The new 3D position in local coordinates (m)

% DEFINE THE POSITION AS A VECTOR INTERVAL
p = [cos(psi)*cos(theta);...
     sin(psi)*cos(theta);...
              sin(theta)]*r;
end

repo: OpenMAS


text:
function [ num_int, pint ] = halfspace_imp_triangle_int_3d ( a, b, c, d, t )

%*****************************************************************************80
%
%% HALFSPACE_IMP_TRIANGLE_INT_3D: intersection ( implicit halfspace, triangle ) in 3D.
%
%  Discussion:
%
%    The implicit form of a half-space in 3D may be described as the set
%    of points (X,Y,Z) on or "above" an implicit plane:
%
%      0 <= A * X + B * Y + C * Z + D
%
%    The triangle is specified by listing its three vertices.
%
%    The intersection may be described by the number of vertices of the
%    triangle that are included in the halfspace, and by the location of
%    points between vertices that separate a side of the triangle into
%    an included part and an unincluded part.
%
%    0 vertices, 0 separators    (no intersection)
%    1 vertex, 0 separators      (point intersection)
%    2 vertices, 0 separators    (line intersection)
%    3 vertices, 0 separators    (triangle intersection)
%
%    1 vertex, 2 separators,     (intersection is a triangle)
%    2 vertices, 2 separators,   (intersection is a quadrilateral).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, C, D, the parameters that define the
%    implicit plane, which in turn define the implicit halfspace.
%
%    Input, real T(3,3), the vertices of the triangle.
%
%    Output, integer NUM_INT, the number of intersection points returned,
%    which will always be between 0 and 4.
%
%    Output, real PINT(3,4), the coordinates of the NUM_INT
%    intersection points.  The points will lie in sequence on the triangle.
%    Some points will be vertices, and some may be separators.
%

%
%  Compute the signed distances between the vertices and the plane.
%
  dist1 = a * t(1,1) + b * t(2,1) + c * t(3,1) + d;
  dist2 = a * t(1,2) + b * t(2,2) + c * t(3,2) + d;
  dist3 = a * t(1,3) + b * t(2,2) + c * t(3,3) + d;
%
%  Now we can find the intersections.
%
  [ num_int, pint ] = halfspace_triangle_int_3d ( t, dist1, dist2, dist3 );

  return
end


text:
function L = gaussianLogLikelihood(noise, mu, varsigma, y)

% GAUSSIANLOGLIKELIHOOD Log-likelihood of data under Gaussian noise model.

% NOISE

N = size(mu, 1);
D = size(mu, 2);
varsigma = varsigma + noise.sigma2;
for i = 1:D
  mu(:, i) = mu(:, i) + noise.bias(i);
end
arg = (y - mu);
arg = arg.*arg./varsigma;

L = - 0.5*sum(sum(log(varsigma))) ...
    - 0.5*sum(sum(arg)) ...
    - 0.5*N*D*log(2*pi);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%% plane_intersect
% Below is a demonstration of the features of the |plane_intersect| function

%%
clear; close all; clc;

%% Syntax
% |X=plane_intersect(V1,V2,V3,N1,N2,N3);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [yd3] = mi32yd3(mi3)
% Convert volume from cubic miles to cubic yards.
% Chad Greene 2012
yd3 = mi3*5451776000;