text:
% timecrossf() - compute the time frequency content of an array rows and
%                cross-coherence between the rows for. If the rows depict
%                different components activation, these results can then
%                be used to analyse time relation between these
%                components  (using the brainmovie function). You may
%                also input electrode activities into this function to
%                analyse the synchronisation and phase relashionship
%                between the electrodes.
%
% Usage:
%    >> [ allersps, allitcs, allcrossfs, allcrossfphis, times, freqs] ...
%        = timecrossf(data, frames, tlimits, srate, cycles, ...);
%
% Inputs:
%   data        - multi-row data (nb_input,frames*nepochs) for the first
%                 condition. Each row is a different input. The function
%                 will conpute the time frequency content of each row and
%                 the cross-coherence between the rows. If data is a cell
%                 array { data1 data2 }, the function will compute the
%                 time-frequency decomposition for these 2 arrays and the
%                 difference between them.
%   frames      - frames per epoch                        {768}
%   tlimits     - epoch time limits (ms) [mintime maxtime]{[-1000 2000]}
%   srate       - data sampling rate (Hz)                 {256}
%   ...         - see timef or crossf command for details about other
%                 parameters
%
% Outputs:
%   allersps      - cell array (size nb_row,1) of ERSP (event-related
%                   power spectrum content of each row of the data array)
%   allitcs       - cell array (size nb_row,1) of ITC (event-related
%                   inter-trial coherence content of each row of the data
%                   array)
%   allcrossfs    - cell array (size nb_row, nb_row) of cross-coherence
%                   between the row of the data array (only the upper
%                   diagonal part of the cell array is used)
%   allcrossfphis - cell array (size nb_row, nb_row) of cross-coherence
%                   phase between the row of the data array (only the
%                   upper diagonal part of the cell array is used)
%   times         - time array returned by the crossf or timef functions
%   freqs         - frequency array returned by the crossf or timef
%                   functions
%
% Example1: generate a movie, 1 condition, below 10Hz dynamics
%
%   % ICA decomposition, activity in array actICA (31 x (176 x n)),
%   % graph for 6 components, 176 points from -100 ms to 600 ms after
%   % stimulus onset, 250 Hz sampling rate
%   [ allersps, allitcs, allcrossfs, allcrossfphis, times, freqs] ...
%        = timecrossf( actICA(1:6,:), 176, [-100 600], 250, 1, 32, 100);
%   brainmovie( ersps, itcs, crossfs_amp, crossfs_phase, times, [1:2] )
%   %[1:2] indicates the frequency rows to take into account (freqs(1:2))
%
% Example2: generate a movie, 2 condition, below 10Hz dynamics
%
%   % Same as above with actICA2 (31 x (176 x n2)) the activity of the
%   % components for a different condition
%   [ allersps, allitcs, allcrossfs, allcrossfphis, times, freqs] ...
%        = timecrossf( actICA(1:6,:), 176, [-100 600], 250, 1, 32, 100);
%   [ allersps(:,2), allitcs(:,2), allcrossfs(:,:,2), ...
%	                          allcrossfphis(:,:,2)] ...
%        = timecrossf( actICA2(1:6,:), 176, [-100 600], 250, 1, 32, 100);
%   brainmovie( allersps, allitcs, allcrossfs, allcrossfphis, times, ...
%                [1:2] )
%
% See also timef(), crossf(), brainmovie()

% arno@salk.edu, Arnaud Delorme, CNL / Salk Institute, 2001

% This program is free software; you can redistribute it and/or
% modify it.
% 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.

function [ ALLERSP, ALLITC, ALLCROSSF, ALLCROSSFANGLE, times, freqs ] = timecrossf(data, frames, tlimits, srate, cycle, varargin);

if nargin < 3
	help timecrossf;
	return;
end;

SAVE    = 0;	% save graphs in jpeg format
%XDISP   = 'cole:0.0';
nbcompo = size( data, 1);
if iscell(data)
    data1 = data{1};
    data2 = data{2};
    nbcompo = size( data1, 1);
end;

% compute all time frequency maps
% -------------------------------
compter = 1;
for numcompo = 1:nbcompo
	if iscell(data)
        [ersp,itc,powbase,times,freqs,erspboot,itcboot] = newtimef({ data1(numcompo,:) data2(numcompo,:) }, ...
                                                          frames, tlimits, srate, cycle, varargin{:});
        size(ersp{1})
        ALLERSP{numcompo,1}     = applyboot( ersp{1}, erspboot{1});
        ALLERSP{numcompo,2}     = applyboot( ersp{2}, erspboot{2});
        ALLERSP{numcompo,3}     = applyboot( ersp{3}, erspboot{3});
        ALLITC {numcompo,1}     = applyboot( abs(itc{1}) , itcboot{1});
        ALLITC {numcompo,2}     = applyboot( abs(itc{2}) , itcboot{2});
        if ~isreal(itc{3}), itc{3} = abs(itc{3}); end;
        ALLITC {numcompo,3}     = applyboot( itc{3}, itcboot{3});
    else
        [ersp,itc,powbase,times,freqs,erspboot,itcboot] = newtimef( data(numcompo,:), ...
                                                          frames, tlimits, srate, cycle, varargin{:});
        ALLERSP{numcompo,1}    = applyboot( ersp, erspboot);
        ALLITC{numcompo,1}     = applyboot( itc,  itcboot);
    end;
	if SAVE
		%command = sprintf('print -djpeg timef%2.2d', numcompo); eval(command); close;
		command = sprintf('hgsave(''timef%2.2d'');', numcompo); eval(command); close;
	end;
	compter = compter + 1 ;
end;

% compute all cross coherence maps
% --------------------------------
ALLCROSSF = cell(nbcompo, nbcompo);
ALLCROSSFANGLE = cell(nbcompo, nbcompo);
for index1 = 1:nbcompo
	for index2 = 1:nbcompo
		if index2 > index1
            if iscell(data)
                [coh,mcoh,timesout,freqsout,cohboot,cohangles] = newcrossf({ data1(index1,:) data2(index1,:)}, ...
                                                                  { data1(index2,:) data2(index2,:)}, frames,  ...
                                                                  tlimits, srate, cycle, varargin{:});
                size(coh{1})
                ALLCROSSF      { index1, index2, 1 } = applyboot(coh{1}, cohboot{1});
                ALLCROSSF      { index1, index2, 2 } = applyboot(coh{2}, cohboot{2});
                ALLCROSSF      { index1, index2, 3 } = applyboot(coh{3}, cohboot{3});
                ALLCROSSFANGLE { index1, index2, 1 } = cohangles{1};
                ALLCROSSFANGLE { index1, index2, 2 } = cohangles{2};
                ALLCROSSFANGLE { index1, index2, 3 } = cohangles{3};
            else
                [coh,mcoh,timesout,freqsout,cohboot,cohangles] = newcrossf(data(index1,:), data(index2,:), frames,  ...
                                                                  tlimits, srate, cycle, varargin{:});
                ALLCROSSF      { index1, index2 } = applyboot(coh, cohboot);
                ALLCROSSFANGLE { index1, index2 } = cohangles;
                if SAVE
                    %command = sprintf('print -djpeg crossf%2.2d-%2.2d%', index1, index2); eval(command); close;
                    command = sprintf('hgsave(''crossf%2.2d-%2.2d%'');', index1, index2); eval(command); close;
                end;
            end;
		end;
	end;
end;

return;

function array = applyboot(array, arrayboot)
    if isempty(arrayboot), return; end;
    if size(arrayboot,3) == 2
        array(find((array > arrayboot(:,:,1)) & (array < arrayboot(:,:,2)))) = 0;
    elseif size(arrayboot,2) > 2
        array(find(array < arrayboot(:,:))) = 0;
    elseif size(arrayboot,2) == 2
        array(find((array > repmat(arrayboot(:,1),[1 size(array,2)])) & ...
                   (array < repmat(arrayboot(:,2),[1 size(array,2)])))) = 0;
    else
        array(find(array < repmat(arrayboot(:,1),[1 size(array,2)]))) = 0;
    end;

author: PatternRecognition

repo: OpenBMI


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

%%
clear; close all; clc;

%% Syntax
% |[varargout]=sv3(varargin);|

%% Description
% 3D slice viewer

%% Example: Visualizing MRI data

%%
% Example image data
load mri;
M=squeeze(D); %example image data set
v=2./[1,1,.4]; %example voxel size

%%
% Using |sv3|

hf=sv3(M,v);

%%
%
% <<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 [idx, dst] = WrapperForKsearch(dataset,testset,k)
% This wrapper is used because we cannot pass an ANN class object to
% parfeval in MATLAB R2017A.
% See reference (https://www.cs.umd.edu/~mount/ANN/Files/1.1.2/ANNmanual_1.1.pdf) returns the squared distances
% ANN: A Library for
% Approximate Nearest Neighbor Searching
% David M. Mount and Sunil Arya
% Version 1.1.2
% Release Date: Jan 27, 2010

% Modified by William to automatically do a projection and get the meters
% distance on a sphere
[dataset(:,1),dataset(:,2)] = m_ll2xy(dataset(:,1),dataset(:,2));
[testset(:,1),testset(:,2)] = m_ll2xy(testset(:,1),testset(:,2));
dataset(isnan(dataset(:,1)),:) = [];
% Find the nearest neighbors
anno = ann(dataset');
idx = ksearch(anno, testset',k,0);
idx = idx';
close(anno);

% The vector of long lat pairs
dst = zeros(length(testset),k);
for kk = 1:k
    long = zeros(length(testset)*2,1);
    lat  = zeros(length(testset)*2,1);
    long(1:2:end) = testset(:,1); long(2:2:end) = dataset(idx(:,kk),1);
    lat(1:2:end) = testset(:,2); lat(2:2:end) = dataset(idx(:,kk),2);
    % Get spherical earth distances
    dt = m_xydist(long,lat);
    dst(:,kk) = dt(1:2:end)*1e3;
end

end

author: CHLNDDEV


text:
% IAIF Glottal Inverse Filtering
%   [g,dg,a,ag] = iaif(x,fs,p_vt,p_gl,d,hpfilt)
%
% Description
%   This function estimates vocal tract linear prediction coefficients and
%   the glottal volume velocity waveform from a speech signal frame using
%   Iterative Adaptive Inverse Filtering (IAIF) method.
%
% Inputs
%   x      : Speech signal frame [samples]
%   fs     : Sampling frequency [Hz]
%   p_vt   : Order of LPC analysis for vocal tract
%   p_gl   : Order of LPC analysis for glottal source
%   d      : Leaky integration coefficient (e.g. 0.99)
%   hpfilt : High-pass filter flag (0: do not apply, 1...N: apply N times)
%
% Outputs
%   g      : Glottal volume velocity waveform
%   dg     : Glottal volume velocity derivative waveform
%   a      : LPC coefficients of vocal tract
%   ag     : LPC coefficients of source spectrum
%
% Notes
%   This function does not perform pitch synchronous analysis. This ensures
%   the robustness of the method regardless of the GCI estimation
%   performance.
%
% Example
%   Simplest, type g = iaif(x,fs) to estimate the glottal flow of a speech
%   frame.
%
% References
%  [1] P. Alku, "Glottal wave analysis with pitch synchronous iterative
%      adaptive inverse filtering", Speech Communication, vol. 11, no. 2-3,
%      pp. 109–118, 1992.
%
% Copyright (c) 2013 Aalto University
%
% License
%  This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file 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 Lesser General Public License for more
%  details.
%
% This function is part of the Common Speech Processing Repository
% http://covarep.github.io/covarep/
%
% Octave Compatible
%
% Author
%  Tuomo Raitio tuomo.raitio@aalto.fi
%
% $Id <info set by the versioning system> $

function [g,dg,a,ag,hpfilter_out] = iaif(x,fs,p_vt,p_gl,d,hpfilt,hpfilter_in)

% Set default parameters
if nargin < 7
    hpfilter_in = [];
    if nargin < 6
        hpfilt = 1;
        if nargin < 5
            d = 0.99;
            if nargin < 4
                p_gl = 2*round(fs/4000);
                if nargin < 3
                    p_vt = 2*round(fs/2000)+4;
                    if nargin < 2
                        disp('Error: Not enough input arguments.');
                    end
                end
            end
        end
    end
end
preflt = p_vt+1;

% Ensure column vector form
x = x(:);

% High-pass filter speech in order to remove possible low frequency
% fluctuations (Linear-phase FIR, Fc = 70 Hz)
hpfilter_out = [];
if hpfilt > 0
    Fstop = 40;                 % Stopband Frequency
    Fpass = 70;                 % Passband Frequency
    Nfir = round(300/16000*fs); % FIR numerator order
    if mod(Nfir,2) == 1
        Nfir = Nfir + 1;
    end

    % it is very very expensive to calculate the firls filter! However, as
    % long as the fs does not change, the firls filter does not change.
    % Therefore, the computed filter is returned and can be passed to this
    % function later on to avoid the calculated of the (same) filter.
    if ~isempty(hpfilter_in)
        B = hpfilter_in;
    else
        B = hpfilter_fir(Fstop,Fpass,fs,Nfir);
    end
    hpfilter_out = B;

    for i = 1:hpfilt
        x = filter(B,1,[x ; zeros(round(length(B)/2)-1,1)]);
        x = x(round(length(B)/2):end);
    end
end

% Estimate the combined effect of the glottal flow and the lip radiation
% (Hg1) and cancel it out through inverse filtering. Note that before
% filtering, a mean-normalized pre-frame ramp is appended in order to
% diminish ripple in the beginning of the frame. The ramp is removed after
% filtering.
if length(x)>p_vt
    win = hanning(length(x));
    signal = [linspace(-x(1),x(1),preflt)' ; x];
    idx = preflt+1:numel(signal);

    Hg1 = lpc(x.*win,1);
    y = filter(Hg1,1,signal);
    y = y(idx);

    % Estimate the effect of the vocal tract (Hvt1) and cancel it out through
    % inverse filtering. The effect of the lip radiation is canceled through
    % intergration. Signal g1 is the first estimate of the glottal flow.
    Hvt1 = lpc(y.*win,p_vt);
    g1 = filter(Hvt1,1,signal);
    g1 = filter(1,[1 -d],g1);
    g1 = g1(idx);

    % Re-estimate the effect of the glottal flow (Hg2). Cancel the contribution
    % of the glottis and the lip radiation through inverse filtering and
    % integration, respectively.
    Hg2 = lpc(g1.*win,p_gl);
    y = filter(Hg2,1,signal);
    y = filter(1,[1 -d],y);
    y = y(idx);

    % Estimate the model for the vocal tract (Hvt2) and cancel it out through
    % inverse filtering. The final estimate of the glottal flow is obtained
    % through canceling the effect of the lip radiation.
    Hvt2 = lpc(y.*win,p_vt);
    dg = filter(Hvt2,1,signal);
    g = filter(1,[1 -d],dg);
    g = g(preflt+1:end);
    dg = dg(idx);

    % Set vocal tract model to 'a' and glottal source spectral model to 'ag'
    a = Hvt2;
    ag = Hg2;
else g=[];
    dg=[];
    a=[];
    ag=[];
    disp('IAIF - frame not analysed!!')
end

function B = hpfilter_fir(Fstop,Fpass,fs,N)
% FIR least-squares Highpass filter design using the FIRLS function
%
% Tuomo Raitio
% 20.8.2013

% Calculate the coefficients using the FIRLS function.
B  = firls(N, [0 Fstop Fpass fs/2]/(fs/2), [0 0 1 1], [1 1]);

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function x7 = testNet(x, net, Params)

    N=size(x,1);
    W=net.W;B=net.B;
    dropP=Params.dropP;
    %% forward
    [x1,~]=dropout(x,dropP(1));%N*V
    x2=x1*W{1}+repmat(B{1},N,1);%N*neuronN
    x3=max(0,x2);%ReLU
    [x4, ~]=dropout(x3,dropP(2));
    x5=x4*W{2}+repmat(B{2},N,1);
    x6=(exp(x5)-exp(-x5))./(exp(x5)+exp(-x5));%Sigmoid%N*neuronN
    x7=x6*W{3}+repmat(B{3},N,1);%N*M
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function circmo()
    %   This subroutine "circle"  selects the Ni closest earthquakes
    %   around a interactively selected point.  Resets ZG.newcat and ZG.newt2
    %   Operates on "primeCatalog".
    %
    % axis: hmo
    % plots to: plos as xk
    % inCatalog: a
    % outCatalog: newt2, newcat
    % mouse controlled
    % closest events
    %
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    %  Input Ni:
    %
    report_this_filefun();
    ZG=ZmapGlobal.Data;

    delete(findobj('Tag','plos1'))

    axes(hmo)

    % interactively get the circle of interest
    shape=ShapeCircle();

    [ZG.newt2, max_km] = selectCircle(newa, shape.toStruct());

    disp(['Radius of selected Circle: ' num2str(max_km)  ' km' ])

    set(gca,'NextPlot','add')
    plot(ZG.newt2.Longitude,ZG.newt2.Latitude,'xk','Tag','plos1');

    ZG.newcat = ZG.newt2;                   % resets ZG.newcat and ZG.newt2

    % plot cumulative number
    ctp=CumTimePlot(ZG.newt2);
    ctp.plot();

end

author: CelsoReyes


text:
classdef EAGMOEAD < ALGORITHM
% <multi> <real/integer/label/binary/permutation>
% External archive guided MOEA/D
% LGs --- 8 --- The number of learning generations

%------------------------------- Reference --------------------------------
% X. Cai, Y. Li, Z. Fan, and Q. Zhang, An external archive guided
% multiobjective evolutionary algorithm based on decomposition for
% combinatorial optimization, IEEE Transactions on Evolutionary
% Computation, 2015, 19(4): 508-523.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            LGs = Algorithm.ParameterSet(8);

            %% Generate the weight vectors
            [W,Problem.N] = UniformPoint(Problem.N,Problem.M);
            T = ceil(Problem.N/10);

            %% Detect the neighbours of each solution
            B = pdist2(W,W);
            [~,B] = sort(B,2);
            B = B(:,1:T);

            %% Generate random population
            Population = Problem.Initialization();
            Archive    = Population;            % External archive
            s          = zeros(Problem.N,LGs);	% Number of successful solutions in last several generations

            %% Optimization
            while Algorithm.NotTerminated(Archive)
                [MatingPool,offspringLoc] = MatingSelection(B,s);
                Offspring  = OperatorGAhalf(Problem,Population(MatingPool));
                Population = UpdatePopulation(Population,Offspring,offspringLoc,W,B);
                [Archive,sucessful] = UpdateArchive(Archive,Offspring);
                % Update the number of successful solutions generated by
                % each subproblem in the last LGs generations
                if any(sucessful)
                    s(:,mod(ceil(Problem.FE/Problem.N),LGs)+1) = hist(offspringLoc(sucessful),1:Problem.N)';
                end
            end
        end
    end
end


text:
function dt6 = dtiMakeDt6FromFsl(b0FileName, t1FileName)
%
% dt6 = dtiMakeDt6FromFsl([b0FileName], [t1FileName])
%
% Takes a set of images from FSL's DTI image files and saves them out as
% in a dt6 matlab file. You should specify the FSL S0 (usually b=0) file as
% input (a file dialog will pop-up if you don't).
%
% The t1 file is optional- if it is omitted, then the b=0 will be used as
% the anatomical reference.
%
% HISTORY:
% 2005.08.04 RFD (bob@white.stanford.edu) wrote it (based on dtiMakeDt6)

slice = [0,0,0]; % slice to display for checking coregistration
bb = [-80,80; -120,90; -60,90]';
mmDt = [2 2 2];
mmAnat = [1 1 1];
autoAlign = false;
showFigs = true;

if ~exist('b0FileName','var') | isempty(b0FileName)
    [f, p] = uigetfile({'*.nii.gz','NIFTI gz';'*.*','All files'}, 'Select the FSL S0 file...');
    if(isnumeric(f)) error('Need an S0 file to continue...'); end
    b0FileName = fullfile(p, f);
    disp(b0FileName);
end
b0 = niftiRead(b0FileName);
% Fix ill-specified quaternion xforms so that they do something reasonable.
if(all(b0.qto_ijk(1:3,4) == [0 0 0]'))
    % *** HACK! [38 69 28] is the ac for one normalized data set- maybe it
    % works for many?
    sz = size(b0.data);
    origin = round(sz(1:3)./2.*[0.95 1.1 0.9]);
    b0.qto_ijk(1:3,4) = origin'+0.5;
    b0.qto_xyz = inv(b0.qto_ijk);
end
b0.data = double(b0.data);
% Clip vals <0 (these should be junk for a b=0 image)
b0.data(b0.data<0) = 0;

if ~exist('t1FileName','var')
    [f, p] = uigetfile({'*.nii.gz','NIFTI gz';'*.hdr','Analyze files (*.hdr)';'*.*','All files'},...
        'Select the T1 file...', b0FileName);
    if(isnumeric(f))
        t1FileName = [];
    else
        t1FileName = fullfile(p, f);
        disp(t1FileName);
    end
end
if(~isempty(t1FileName))
    t1 = niftiRead(t1FileName);
    % Fix ill-specified quaternion xforms so that they do something reasonable.
    if(all(t1.qto_ijk(1:3,4) == [0 0 0]'))
        sz = size(t1.data);
        origin = round(sz(1:3)./2.*[0.95 1.1 0.9]);
        t1.qto_ijk(1:3,4) = origin'+0.5;
        t1.qto_xyz = inv(t1.qto_ijk);
    end
    t1.data = double(t1.data);
    t1.data(t1.data<0) = 0;
else
    t1 = b0;
end

[datapath,basename] = fileparts(b0FileName);
us = strfind(basename,'_');
basename = basename(1:us(end));

if(~exist('outPathName','var') | isempty(outPathName))
    outPathName = fullfile(datapath, [basename 'dt6.mat']);
end
if(exist(outPathName,'file') | exist([outPathName '.mat'],'file'))
  disp('output file exists- please rename...');
  [f,p] = uiputfile('*.mat', 'Select output file...', outPathName);
  if(isnumeric(f)) error('User cancelled.'); end
  outPathName = fullfile(p,f);
end

disp(['Data will be saved in ' outPathName '...']);

% *** TO DO: fix this- here we assume that the data use a NIFTI-style
% header and this it's quaternion xform is specified correctly.
t1.acpcXform = t1.qto_xyz;
b0.acpcXform = b0.qto_xyz;
% reslicing introduces a 1/2 voxel (1mm shift. We correct for that here.
b0.acpcXform(1:3,4) = b0.acpcXform(1:3,4)-mmDt'./2;

h = figure;
showFigure(h, t1, b0, bb, slice);
%set(h,'Position', [10, 50, 1000, 1000]);
set(h, 'PaperPositionMode', 'auto');
print(h, '-dpng', '-r90', [outPathName,'_headerAlign.png']);
if(~showFigs) close(h); end

if(autoAlign)
    disp('Coregistering (using spm2 tools)...');
    img = t1.data;
    img = mrAnatHistogramClip(img, 0.50, 0.995);
    VG.uint8 = uint8(round(img*255));
    VG.mat = t1.acpcXform;

    img = b0.data;
    img = mrAnatHistogramClip(img,0.40, 0.99);
    VF.uint8 = uint8(round(img*255));
    VF.mat = b0.acpcXform;
    p = defaults.coreg.estimate;
    % NOTE: there seems to be a consistent 1mm (1/2 of B0 voxel) translation
    % between the t1 and B0.
    transRot = spm_coreg(VG,VF,p);
    transRot(1:3) = transRot(1:3)+mmDt/2;
    b0.acpcXform = inv(VF.mat\spm_matrix(transRot(:)'));
    fig = figure;
    showFigure(fig, t1, b0, bb, [0,0,0], ['auto aligned']);
    set(fig, 'PaperPositionMode', 'auto');
    print(fig, '-dpng', '-r90', [outPathName,'_autoAlign.png']);
end

% Now build and save the dt6 data
eigVal = zeros([size(b0.data) 3]);
eigVec = zeros([size(b0.data) 3 3]);
for(ii=1:3)
    tmp = niftiRead(fullfile(datapath, [basename 'L' num2str(ii) '.nii.gz']));
    eigVal(:,:,:,ii) = double(tmp.data);
end
for(ii=1:3)
    tmp = niftiRead(fullfile(datapath, [basename 'V' num2str(ii) '.nii.gz']));
    eigVec(:,:,:,:,ii) = double(tmp.data);
    eigVec(:,:,:,2,ii) = -eigVec(:,:,:,2,ii);
end

dt6 = dtiRebuildTensor(eigVec, eigVal);
clear eigVal eigVec;
%dt6_tmp = dtiLoadTensorElements(fullfile(fileparts(b0.baseFname), 'TensorElements.float.'));
%%dt6_tmp = permute(dt6_tmp,[2 1 3 4]);
%for(ii=1:6)
%    dt6(:,:,:,ii) = applyCannonicalXform(dt6_tmp(:,:,:,ii), b0.cannonical, b0.mmPerVox);
%end
%dt6 = permute(dt6,[2 1 3 4]);
%clear dt6_tmp;

% Reslice everything
disp('Interpolating tensors...');
%dt6A = dtiResliceTensorAffine(dt6, inv(b0.acpcXform), b0.cannonical_mmPerVox, bb, mmDt);
%dt6A = permute(dt6A,[2 1 3 4]);
dt6 = mrAnatResliceSpm(dt6, inv(b0.acpcXform), bb, mmDt, [7 7 7 0 0 0], showFigs);
dt6(isnan(dt6)) = 0;
% NOTE: we want to apply b0.acpcXform to the tensors, even though we
% resliced them with inv(b0.acpcXform). inv(b0.acpcXform) maps from
% the new space to the old- the correct mapping for the interpolation,
% since we interpolate by creating a grid in the new space and fill it by
% pulling data from the old space. But the tensor reorientation should be
% done using the old-to-new space mapping (b0.acpcXform).
rigidXform = dtiFiniteStrainDecompose(b0.acpcXform);
[t,r] = affineDecompose([rigidXform,[0 0 0]';[0 0 0 1]]);
if(all(all(rigidXform==eye(3))))
    disp('No PPD correction needed- rotation is zero.');
else
    fprintf('Applying PPD rotation [%0.4f %0.4f %0.4f]...\n',r);
    dt6 = dtiXformTensors(dt6, rigidXform);
end
%figure; imagesc(makeMontage(dt6_new(:,:,:,1))); axis equal; colormap gray

disp('Interpolating B0...');
[b0_img,newXform] = mrAnatResliceSpm(b0.data, inv(b0.acpcXform), bb, mmDt, [7 7 7 0 0 0], showFigs);
b0_img(b0_img<0) = 0;
b0_img(isnan(b0_img)) = 0;

% create and apply a brain mask to the dt6 data
brainMask = mrAnatHistogramClip(b0_img, 0.4, 0.99);
brainMask = brainMask > 0.1;
% [c,v] = mrAnatHistogramSmooth(b0_img,256,0.05);
% % derivative of the thresholded derivative gives us the locations where a
% % rising or falling trend begins. In general, we can just take the first
% % rising trend (ie. just after the background noise) to find the brain.
% % We then back off from that by 25% (an empirically determined heuristic).
% peakStart = find(diff(diff(c)>0)>0);
% thresh = v(peakStart(1)+2);
% thresh = thresh.*0.75;
% brainMask = b0_img > thresh;
brainMask = dtiCleanImageMask(brainMask, 9);
figure; imagesc(makeMontage(brainMask)); axis equal tight; colormap(gray);
title('Brain Mask (based on b=0 image)');
dt6(repmat(~brainMask, [1,1,1,6])) = 0;

disp('Interpolating T1...');
anat.img = mrAnatResliceSpm(t1.data, inv(t1.acpcXform), bb, mmAnat, [], showFigs);
anat.img(anat.img<0) = 0;
anat.img(isnan(anat.img)) = 0;
anat.img = int16(round(anat.img));

mmPerVox = mmDt;
notes = 'Built from FSL data';
b0 = int16(round(b0_img));
anat.mmPerVox = mmAnat;

origin = bb(1,[1,2,3])-mmDt./2;
xformToAnat = [diag(mmDt./mmAnat) [0 0 0]'; [0 0 0 1]];
xformToAcPc = [diag(mmDt) origin'; [0 0 0 1]]; %swapXY*[diag(mm) origin'; [0 0 0 1]];
anat.xformToAcPc = [diag(mmAnat) origin'; [0 0 0 1]];
% rp = t1.talairach.refPoints;
% % Why swap coords like this? Well, the real reason has to do with the
% % history of ComputeTalairach. We could figure this out in a more general
% % way using t1.talairach.refPoints.mat, which tells us how the coords were
% % transformed in ComputeTalairach. In the end, the following reordering
% % works for our data, since the pre-processing is consistent.
% mm = t1.cannonical_mmPerVox([2 3 1]);
anat.talScale.notes = ['Scale factors to go from subject mm to Talairach mm.' ...
                       'Eg. talCoord = anat.talScale.sac * imgCoord.' ...
                       'The scales are (in order) superior, inferior, left, right ' ...
                       'and anterior of AC, ac-to-pc and posterior of pc.'];
% anat.talScale.sac = 72./sqrt(sum(((rp.acXYZ-rp.sacXYZ).*mm).^2));
% anat.talScale.iac = 42./sqrt(sum(((rp.acXYZ-rp.iacXYZ).*mm).^2));
% anat.talScale.lac = 62./sqrt(sum(((rp.acXYZ-rp.lacXYZ).*mm).^2));
% anat.talScale.rac = 62./sqrt(sum(((rp.acXYZ-rp.racXYZ).*mm).^2));
% anat.talScale.aac = 68./sqrt(sum(((rp.acXYZ-rp.aacXYZ).*mm).^2));
% anat.talScale.acpc = 24./sqrt(sum(((rp.acXYZ-rp.pcXYZ).*mm).^2));
% % The PPC is referenced to the PC. It is at Talairach (0,-102,0) and the PC
% % is at (0,-24,0), so it is 78 mm beyond the PC.
% anat.talScale.ppc = 78./sqrt(sum(((rp.pcXYZ-rp.ppcXYZ).*mm).^2));
anat.talScale.sac = 1.0;
anat.talScale.iac = 1.0;
anat.talScale.lac = 1.0;
anat.talScale.rac = 1.0;
anat.talScale.aac = 1.0;
anat.talScale.acpc = 1.0;
anat.talScale.ppc = 1.0;

disp(['Saving to ' outPathName '...']);
l = license('inuse'); h = hostid;
created = ['Created at ' datestr(now,31) ' by ' l(1).user ' on ' h{1} '.'];
save(outPathName, 'dt6', 'mmPerVox', 'notes', 'xformToAnat', 'xformToAcPc', 'b0', 'anat', 'created');
if(nargout==0)
    clear dt6;
end
return;

function showFigure(fig, t1, b0, bb, slice, figName)
if(~exist('figName','var')) figName = 'Interpolated slices'; end
% Get X,Y and Z (L-R, A-P, S-I) slices from T1 and b0 volumes
[t1Xsl] = dtiGetSlice(t1.acpcXform,t1.data,3,slice(3),bb);
[t1Ysl] = dtiGetSlice(t1.acpcXform,t1.data,2,slice(2),bb);
[t1Zsl] = dtiGetSlice(t1.acpcXform,t1.data,1,slice(1),bb);
[b0Xsl] = dtiGetSlice(b0.acpcXform,b0.data,3,slice(3),bb);
[b0Ysl] = dtiGetSlice(b0.acpcXform,b0.data,2,slice(2),bb);
[b0Zsl] = dtiGetSlice(b0.acpcXform,b0.data,1,slice(1),bb);
% Max values for image scaling
t1mv = max([t1Xsl(:); t1Ysl(:); t1Zsl(:)])+0.000001;
b0mv = max([b0Xsl(:); b0Ysl(:); b0Zsl(:)])+0.000001;
% Create XxYx3 RGB images for each of the axis slices. The green and
% blue channels are from the T1, the red channel is an average of T1
% and b=0.
Xsl(:,:,1) = t1Xsl./t1mv.*.5 + b0Xsl./b0mv.*.5;
Xsl(:,:,2) = t1Xsl./t1mv.*.5; Xsl(:,:,3) = t1Xsl./t1mv.*.5;
Ysl(:,:,1) = t1Ysl./t1mv.*.5 + b0Ysl./b0mv.*.5;
Ysl(:,:,2) = t1Ysl./t1mv.*.5; Ysl(:,:,3) = t1Ysl./t1mv.*.5;
Zsl(:,:,1) = t1Zsl./t1mv.*.5 + b0Zsl./b0mv.*.5;
Zsl(:,:,2) = t1Zsl./t1mv.*.5; Zsl(:,:,3) = t1Zsl./t1mv.*.5;

% Show T1 slices
figure(fig); set(fig, 'NumberTitle', 'off', 'Name', figName);
figure(fig); subplot(3,3,1); imagesc(bb(:,1), bb(:,2), t1Xsl);
colormap(gray); axis equal tight xy;
figure(fig); subplot(3,3,2); imagesc(bb(:,3), bb(:,1), t1Ysl);
colormap(gray); axis equal tight xy;
figure(fig); subplot(3,3,3); imagesc(bb(:,3), bb(:,2), t1Zsl);
colormap(gray); axis equal tight xy;
axis equal tight;

% Show b=0 slices
figure(fig); subplot(3,3,4); imagesc(bb(:,1), bb(:,2), b0Xsl);
colormap(gray); axis equal tight xy;
figure(fig); subplot(3,3,5); imagesc(bb(:,3), bb(:,1), b0Ysl);
colormap(gray); axis equal tight xy;
figure(fig); subplot(3,3,6); imagesc(bb(:,3), bb(:,2), b0Zsl);
colormap(gray); axis equal tight xy;
axis equal tight;

% Show combined slices
figure(fig); subplot(3,3,7); imagesc(bb(:,1), bb(:,2), Xsl);
axis equal tight xy;
figure(fig); subplot(3,3,8); imagesc(bb(:,3), bb(:,1), Ysl);
axis equal tight xy;
figure(fig); subplot(3,3,9); imagesc(bb(:,3), bb(:,2), Zsl);
axis equal tight xy;
axis equal tight;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function test_tutorial_networkanalysis

% MEM 3gb
% WALLTIME 00:30:00
% DEPENDENCY ft_networkanalysis

%% read the continuous data and segment into 2 seconds epochs
cfg = [];
cfg.dataset              = dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/SubjectRest.ds');
cfg.trialdef.triallength = 2;
cfg.trialdef.ntrials     = Inf;

cfg = ft_definetrial(cfg);
cfg.continuous  = 'yes';
cfg.channel     = {'MEG'};
data = ft_preprocessing(cfg);

datadir = dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/networkanalysis');
cd(datadir);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BYPASS THE INTERACTIVE PART, DECLARE TRIALS 38 91 153 AS BAD
%% make a visual inspection and reject bad trials/sensors
%cfg = [];
%cfg.method  = 'summary';
%cfg.channel = 'MEG';
%cfg.layout  = 'CTF275.lay';
%dataclean = ft_rejectvisual(cfg, data);
%
%% you can check the rejected trial numbers by typing
%trlind = [];
%for i=1:length(dataclean.cfg.artfctdef.summary.artifact)
%  trlind(i) = find(data.sampleinfo(:,1)==dataclean.cfg.artfctdef.summary.artifact(i));
%end;
%disp(trlind);

cfg = [];
cfg.trials = setdiff(1:numel(data.trial), [38 91 153]);
dataclean  = ft_selectdata(cfg, data);

%% downsample the data to speed up component analysis
cfg = [];
cfg.resamplefs = 60;
cfg.detrend    = 'yes';
datads = ft_resampledata(cfg, dataclean);

%% use ICA in order to identify cardiac and blink components
cfg = [];
cfg.method          = 'runica';
cfg.runica.maxsteps = 50;
cfg.randomseed      = 0;
comp = ft_componentanalysis(cfg, datads);
%load(fullfile(datadir,'comp.mat'));

%% visualize components

% these were the indices of the bad components that were identified
% they may be different if you re-run the ICA decomposition
badcomp = [1 3 8 23];

cfg = [];
cfg.channel    = badcomp;
cfg.layout     = 'CTF275.lay';
cfg.compscale  = 'local';
cfg.continuous = 'yes';
ft_databrowser(cfg, comp);

cfg = [];
cfg.component = badcomp;
dataica = ft_rejectcomponent(cfg, comp);

%% compute the power spectrum
cfg              = [];
cfg.output       = 'pow';
cfg.method       = 'mtmfft';
cfg.taper        = 'dpss';
cfg.foilim       = [1 20];
cfg.tapsmofrq    = 2;
cfg.keeptrials   = 'no';
fft_data = ft_freqanalysis(cfg, dataica);

%% compute the planar transformation
load ctf275_neighb; % loads the neighbourhood structure for the channels

dataicatmp = dataica;
dataicatmp.grad = data.grad;

cfg = [];
cfg.neighbours    = neighbours;
cfg.planarmethod  = 'sincos';
planar = ft_megplanar(cfg, dataicatmp);
clear dataicatmp;

%% compute the power spectrum
cfg              = [];
cfg.output       = 'pow';
cfg.method       = 'mtmfft';
cfg.taper        = 'dpss';
cfg.foilim       = [1 20];
cfg.tapsmofrq    = 2;
cfg.keeptrials   = 'no';
fft_data_planar = ft_freqanalysis(cfg, planar);

%% plot the topography and the spectrum
figure;

cfg = [];
cfg.layout = 'CTF275.lay';
cfg.xlim   = [9 11];
subplot(2,1,1); ft_topoplotER(cfg, fft_data);

cfg = [];
cfg.channel = {'MRO22', 'MRO32', 'MRO33'};
subplot(2,1,2); ft_singleplotER(cfg, fft_data);

%% load the required geometrical information

%% load the template source model, which is in MNI coordinates
%template = load('standard_sourcemodel3d2cm.mat');
% load mri % individual mri
% load hdm % individual volume model
%
% %% compute the source model
% cfg = [];
% cfg.sourcemodel.warpmni   = 'yes';
% cfg.sourcemodel.template  = template.sourcemodel;
% cfg.sourcemodel.nonlinear = 'yes'; % use non-linear normalization
% cfg.mri            = mri;
% sourcemodel        = ft_prepare_sourcemodel(cfg);

load hdm
load sourcemodel_4k

%% check for the correct alignment of sensors (green) headmodel(transparent) and sourcemodel(blue)
figure;

% make the headmodel surface transparent
ft_plot_headmodel(hdm, 'edgecolor', 'none');
alpha 0.4

% add the source model positions and sensors
ft_plot_mesh(ft_convert_units(sourcemodel, 'cm'),'vertexcolor',sourcemodel.sulc);
ft_plot_sens(dataclean.grad);

view([0 -90 0])

%% compute sensor level Fourier spectra
cfg            = [];
cfg.method     = 'mtmfft';
cfg.output     = 'fourier';
cfg.keeptrials = 'yes';
cfg.tapsmofrq  = 2;
cfg.foi        = 10;
freq           = ft_freqanalysis(cfg, dataica);

%% compute the leadfield
cfg             = [];
cfg.sourcemodel        = sourcemodel;
cfg.headmodel   = hdm;
cfg.channel     = {'MEG'};
lf = ft_prepare_leadfield(cfg, freq);

%% compute the actual source reconstruction
cfg                   = [];
cfg.frequency         = freq.freq;
cfg.grad              = freq.grad;
cfg.method            = 'pcc';
cfg.sourcemodel              = lf;
cfg.headmodel         = hdm;
cfg.keeptrials        = 'yes';
cfg.pcc.lambda        = '10%';
cfg.pcc.projectnoise  = 'yes';
cfg.pcc.fixedori     = 'yes';
source = ft_sourceanalysis(cfg, freq);

%% reduce the source reconstructed data to the dominant orientation
cfg = [];
cfg.projectmom = 'yes';
source_proj = ft_sourcedescriptives(cfg,source);

% and provide the dimension and grid positions of the MRI template positions again
source_proj.dim = template.sourcemodel.dim;
source_proj.pos = template.sourcemodel.pos;

[ftver, ftdir] = ft_version;
if isunix
   templatefile = [ftdir '/template/anatomy/single_subj_T1.nii'];
elseif ispc
  templatefile =  [ftdir '\template\anatomy\single_subj_T1.nii'];
end

template_mri = ft_read_mri(templatefile);

cfg              = [];
cfg.parameter    = 'nai';
cfg.interpmethod = 'nearest';
source_int  = ft_sourceinterpolate(cfg, source_proj, template_mri);

%% plot the neural activity index (power/noise)
cfg               = [];
cfg.method        = 'ortho';
cfg.funparameter  = 'nai';
cfg.maskparameter = cfg.funparameter;
cfg.funcolorlim   = [0.0 10];
cfg.opacitylim    = [3 10];
cfg.opacitymap    = 'rampup';
cfg.funcolormap   = 'jet';
ft_sourceplot(cfg, source_int);

cfg               = [];
cfg.method        = 'ortho';
cfg.funparameter  = 'nai';
cfg.location      = [16 -80 -2];
cfg.maskparameter = cfg.funparameter;
cfg.funcolorlim   = [0.0 10];
cfg.opacitylim    = [3 10];
cfg.opacitymap    = 'rampup';
cfg.funcolormap   = 'jet';
ft_sourceplot(cfg,source_int);

%% add the scalp topographies to the figure

cfg = [];
cfg.layout = 'CTF275.lay';
cfg.xlim   = [9.777650 11.309908];
subplot(2,2,2); ft_topoplotER(cfg,fft_data);
title('axial gradient')

cfg = [];
fft_data_planar_cmb = ft_combineplanar(cfg, fft_data_planar);

cfg = [];
cfg.layout = 'CTF275.lay';
cfg.xlim   = [9.777650 11.309908];
subplot(2,2,3); ft_topoplotER(cfg, fft_data_planar_cmb);
title('planar gradient')

%% compute the power spectrum again but keep the individual trials
cfg              = [];
cfg.output       = 'pow';
cfg.method       = 'mtmfft';
cfg.taper        = 'dpss';
cfg.foilim       = [8 12];
cfg.tapsmofrq    = 2;
cfg.keeptrials   = 'yes';
fft_data = ft_freqanalysis(cfg, dataica);

%% identify the indices of trials with high and low alpha power
tmp     = mean(fft_data.powspctrm,3);           % mean over frequencies between 8-12Hz
ind     = find(mean(tmp,1)==max(mean(tmp,1)));  % find the sensor where power is max
indlow  = find(tmp(:,ind)<=median(tmp(:,ind)));
indhigh = find(tmp(:,ind)>=median(tmp(:,ind)));

%% compute the planar transformation
cfgneigh          = [];
cfgneigh.feedback = 'no';
cfgneigh.method   = 'template';

cfg = [];
cfg.neighbours    = ft_prepare_neighbours(cfgneigh, dataica);
cfg.planarmethod  = 'sincos';
planar = ft_megplanar(cfg, dataica);

%% compute the power spectrum
cfg              = [];
cfg.output       = 'pow';
cfg.method       = 'mtmfft';
cfg.taper        = 'dpss';
cfg.foilim       = [1 20];
cfg.tapsmofrq    = 2;
cfg.keeptrials   = 'no';

cfg.trials       = indlow;
tmp = ft_freqanalysis(cfg, planar);
fft_data_planar_low = ft_combineplanar([], tmp);

cfg.trials       = indhigh;
tmp = ft_freqanalysis(cfg, planar);
fft_data_planar_high = ft_combineplanar([], tmp);

%% and also do the axial representation for comparison purposes
cfg              = [];
cfg.output       = 'pow';
cfg.method       = 'mtmfft';
cfg.taper        = 'dpss';
cfg.foilim       = [1 20];
cfg.tapsmofrq    = 2;
cfg.keeptrials   = 'no';
cfg.trials       = indlow;
fft_data_low = ft_freqanalysis(cfg,dataica);

cfg.trials       = indhigh;
fft_data_high = ft_freqanalysis(cfg,dataica);

%% compute the difference between high and low
cfg = [];
cfg.parameter = 'powspctrm';
cfg.operation = 'subtract';
diff        = ft_math(cfg, fft_data_high,        fft_data_low);
diff_planar = ft_math(cfg, fft_data_planar_high, fft_data_planar_low);
diff_axial  = ft_math(cfg, fft_data_high,        fft_data_low);

%% plot the topography of the difference along with the spectra
figure;

cfg = [];
cfg.layout = 'CTF275.lay';
cfg.xlim   = [9.777650 11.309908];

subplot(1,3,1); ft_topoplotER(cfg, diff_planar); title ('planar gradient')

subplot(1,3,2); ft_topoplotER(cfg, diff_axial); title ('axial gradient')

cfg = [];
cfg.channel = {'MLO21', 'MRO31'};
subplot(1,3,3); ft_singleplotER(cfg,fft_data_high, fft_data_low);

obj = subplot(1,3,3);
set(obj, 'Position', [.7 .37 .2 .3])
title('')
legend('high alpha','low alpha','Location','northoutside','Orientation','horizontal');

%% compute fourier spectra for frequency of interest
cfg            = [];
cfg.method     = 'mtmfft';
cfg.output     = 'fourier';
cfg.keeptrials = 'yes';
cfg.tapsmofrq  = 2;
cfg.foi        = 10;

cfg.trials     = indlow;
freq_low = ft_freqanalysis(cfg, dataica);

cfg.trials     = indhigh;
freq_high = ft_freqanalysis(cfg, dataica);

%% compute the beamformer filters based on the entire data
cfg                   = [];
cfg.frequency         = freq.freq;
cfg.grad              = freq.grad;
cfg.method            = 'pcc';
cfg.sourcemodel              = lf;
cfg.headmodel         = hdm;
cfg.keeptrials        = 'yes';
cfg.pcc.lambda        = '5%';
cfg.pcc.projectnoise  = 'yes';
cfg.pcc.keepfilter    = 'yes';
source = ft_sourceanalysis(cfg, freq);

% use the precomputed filters
cfg                   = [];
cfg.frequency         = freq.freq;
cfg.grad              = freq.grad;
cfg.method            = 'pcc';
cfg.sourcemodel              = lf;
cfg.sourcemodel.filter       = source.avg.filter;
cfg.headmodel         = hdm;
cfg.keeptrials        = 'yes';
cfg.pcc.lambda        = '5%';
cfg.pcc.projectnoise  = 'yes';
source_low  = ft_sourceanalysis(cfg, freq_low);
source_high = ft_sourceanalysis(cfg, freq_high);

%% project dipole moments along the dominant orientation
cfg = [];
cfg.projectmom = 'yes';
source_proj_high = ft_sourcedescriptives(cfg,source_high);
source_proj_low  = ft_sourcedescriptives(cfg,source_low);
source_proj_high.dim = template.sourcemodel.dim;
source_proj_high.pos = template.sourcemodel.pos;
source_proj_low.dim  = template.sourcemodel.dim;
source_proj_low.pos  = template.sourcemodel.pos;

%% interpolate the results onto an anatomical brain template
[ftver, ftdir] = ft_version;
if isunix
   templatefile = [ftdir '/template/anatomy/single_subj_T1.nii'];
elseif ispc
  templatefile =  [ftdir '\template\anatomy\single_subj_T1.nii'];
end

template_mri = ft_read_mri(templatefile);

cfg              = [];
cfg.parameter    = 'pow';
cfg.interpmethod = 'nearest';
source_int_high  = ft_sourceinterpolate(cfg, source_proj_high, template_mri);
source_int_low   = ft_sourceinterpolate(cfg, source_proj_low,  template_mri);

cfg  = [];
cfg.operation = '(x1-x2)/x2';
cfg.parameter = 'pow';
source_int = ft_math(cfg, source_int_high, source_int_low);

cfg  = [];
cfg.operation = '(x1-x2)/x2';
cfg.parameter = 'pow';
source_int = ft_math(cfg, source_int_high, source_int_low);

% up to 50 percent of maximum
source_int.mask = source_int.pow > max(source_int.pow(:))*.5;

% copy the anatomy
source_int.anatomy = source_int_high.anatomy;

cfg = [];
cfg.method        = 'ortho';
cfg.funparameter  = 'pow';
cfg.maskparameter = 'mask';
cfg.funcolorlim   = [-1 1];
cfg.location      = [16 -78 38];
cfg.funcolormap   = 'jet';
ft_sourceplot(cfg, source_int);

cfg = [];
cfg.layout = 'CTF275.lay';
cfg.xlim   = [9.777650 11.309908];
subplot(2,2,2); ft_topoplotER(cfg, diff_planar);

%% reduce memory demands and compute connectivity

% compute the sparse representation
source_sparse = ft_source2sparse(source_proj);

% then compute connectivity
cfg=[];
cfg.method  ='coh';
cfg.complex = 'absimag';
source_conn = ft_connectivityanalysis(cfg, source_sparse);

%% reassign the connectivity
cohspctrm_full = nan(size(source_proj.pos,1));
for i=1:size(source_conn.pos)
  pos1 = source_conn.pos(i,1:3);
  pos2 = source_conn.pos(i,4:6);
  ind1 = find(source_proj.pos(:,1)==pos1(1) & source_proj.pos(:,2)==pos1(2) & source_proj.pos(:,3)==pos1(3));
  ind2 = find(source_proj.pos(:,1)==pos2(1) & source_proj.pos(:,2)==pos2(2) & source_proj.pos(:,3)==pos2(3));
  cohspctrm_full(ind1,ind2) = source_conn.cohspctrm(i);
end

% then go to the 'full' representation again
source_conn_full = [];
source_conn_full.pos       = source_proj.pos;
source_conn_full.dim       = source_proj.dim;
source_conn_full.inside    = source_proj.inside;
source_conn_full.cohspctrm = cohspctrm_full;
source_conn_full.dimord    = 'pos_pos';

%% compute graph metric
cfg = [];
cfg.method    = 'degrees';
cfg.parameter = 'cohspctrm';
cfg.threshold = .1;
network = ft_networkanalysis(cfg,source_conn_full);

%% interpolate on anatomical mri and plot the result
network.pos     = template.sourcemodel.pos;
network.dim     = template.sourcemodel.dim;
network.inside  = template.sourcemodel.inside;

cfg              = [];
cfg.parameter    = 'degrees';
cfg.interpmethod = 'nearest';
network_int  = ft_sourceinterpolate(cfg, network, template_mri);

%%
cfg               = [];
cfg.method        = 'ortho';
cfg.funparameter  = 'avg.degrees';
cfg.funcolormap   = 'jet';
cfg.location      = 'max';
cfg.funparameter  = 'degrees';
ft_sourceplot(cfg, network_int);

% we increase the threshold here to highlight the dominant links
edge = source_conn.cohspctrm >.15;
dlmwrite('edge.edge',edge,'\t');

load standard_sourcemodel3d2cm.mat

node = zeros(372,6);
node(:,1:3) = sourcemodel.pos(sourcemodel.inside,:);
node(:,4)   = 4;
node(:,5)   = network.degrees(network.inside);
node(:,6)   = 0;
node(:,1:3) = node(:,1:3)*10;
dlmwrite('node.node',node,' ');

BrainNet_MapCfg('/yourpath/mesh.nv','/yourpath/node.node','/yourpath/edge.edge');
view([0 -90 0])

cfg = [];
cfg.method = 'plv';
source_conn = ft_connectivityanalysis(cfg, source_sparse);

%% then compute graph metric
source_conn.dim     = source_sparse.dim;
source_conn.outside = source_proj.outside;

% then go to the 'full' representation again
source_conn_full = ft_source2full(source_conn);
source_conn_full.dimord='pos_pos';

%%
fn=fieldnames(source_conn_full);
cfg = [];
cfg.method = 'degrees';
cfg.parameter = fn{4};
cfg.threshold = .5;
deg = ft_networkanalysis(cfg,source_conn_full);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [Q, trace] = Update(Q , V , V2, knn , p , r , a, ap , trace, alpha, gamma, lambda, f )
% Update update de Qtable
% V: previous values from state before taking action (a)
% V2: current values after action (a)
% r: reward received from the environment after taking action (a) in state
%                                             s1 and reaching the state s2
% knn: the features in previous states
% p:  the features coeficients (probabilities a.k.a. unit activations)
% a:  the last executed action
% Q: the current Q-table
% alpha: learning rate
% gamma: discount factor
% lambda: eligibility trace decay factor
% trace : eligibility trace vector
% f: true if episode ends

trace(knn,:) = 0.0; %optional trace reset
trace(knn,a) = p;

if f==true
     delta  =  r - V(a);
else
    delta  =  ( r + gamma .* V2(ap) ) - V(a);
end

Q      =  Q  + alpha .* delta .* trace;
trace  =  gamma * lambda .* trace;
end

author: epfl-lasa


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

%*****************************************************************************80
%
%% BESSEL_J1_VALUES returns some values of the J1 Bessel function.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      BesselJ[1,x]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 August 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  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 = 21;

  fx_vec = [ ...
      0.3275791375914652E+00, ...
      0.6604332802354914E-01, ...
     -0.3390589585259365E+00, ...
     -0.5767248077568734E+00, ...
     -0.4400505857449335E+00, ...
      0.0000000000000000E+00, ...
      0.4400505857449335E+00, ...
      0.5767248077568734E+00, ...
      0.3390589585259365E+00, ...
     -0.6604332802354914E-01, ...
     -0.3275791375914652E+00, ...
     -0.2766838581275656E+00, ...
     -0.4682823482345833E-02, ...
      0.2346363468539146E+00, ...
      0.2453117865733253E+00, ...
      0.4347274616886144E-01, ...
     -0.1767852989567215E+00, ...
     -0.2234471044906276E+00, ...
     -0.7031805212177837E-01, ...
      0.1333751546987933E+00, ...
      0.2051040386135228E+00 ];

  x_vec = [ ...
     -5.0E+00, ...
     -4.0E+00, ...
     -3.0E+00, ...
     -2.0E+00, ...
     -1.0E+00, ...
      0.0E+00, ...
      1.0E+00, ...
      2.0E+00, ...
      3.0E+00, ...
      4.0E+00, ...
      5.0E+00, ...
      6.0E+00, ...
      7.0E+00, ...
      8.0E+00, ...
      9.0E+00, ...
     10.0E+00, ...
     11.0E+00, ...
     12.0E+00, ...
     13.0E+00, ...
     14.0E+00, ...
     15.0E+00 ];

  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

author: johannesgerer

repo: jburkardt-m


text:
classdef Gene < handle
    %indices of cellular constituents
    properties (SetAccess = private)
        wholeCellModelIDs     %Whole Cell model IDs
        names                 %names
        startCoordinates      %start coordinates
        lengths               %lengths
        strands               %strands
        essential             %is gene essential (Y=>yes, N=>no, M=>maybe)

        mRNAIndexs            %mRNA genes within geneWholeCellModelIDs
        rRNAIndexs            %rRNA genes within geneWholeCellModelIDs
        sRNAIndexs            %sRNA genes within geneWholeCellModelIDs
        tRNAIndexs            %tRNA genes within geneWholeCellModelIDs
        ribosomalRRNAIndexs   %ribosomal RNA genes within geneWholeCellModelIDs
    end

    methods
        function this = Gene()
        end

        function initializeConstants(this, knowledgeBase, ~)
            this.wholeCellModelIDs   = {knowledgeBase.genes.wholeCellModelID}';
            this.names               = {knowledgeBase.genes.name}';
            this.startCoordinates    = [knowledgeBase.genes.startCoordinate]';
            this.lengths             = [knowledgeBase.genes.sequenceLength]';
            this.strands             = 2-[knowledgeBase.genes.direction]';
            this.essential           = {knowledgeBase.genes.essential}';

            this.mRNAIndexs          = double([knowledgeBase.mRNAGenes.idx])';
            this.rRNAIndexs          = double([knowledgeBase.rRNAGenes.idx])';
            this.sRNAIndexs          = double([knowledgeBase.sRNAGenes.idx])';
            this.tRNAIndexs          = double([knowledgeBase.tRNAGenes.idx])';
            this.ribosomalRRNAIndexs = double(knowledgeBase.ribosomalRRNAIndexs);
        end
    end

    %helper functions
    methods
        function value = getIndexs(this, wholeCellModelIDs)
            [~, value] = ismember(wholeCellModelIDs, this.wholeCellModelIDs);
        end

        function value = getIndexsByPosition(this, pos)
            value = find(...
                pos >= this.startCoordinates & ...
                pos <= this.startCoordinates + this.lengths - 1 ...
                , 1);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
% function to overwrite ops1 fields with ops0 fields if settings have
% changed for new runs
function ops = opsChanges(ops1, ops0)

ops = ops1;
f = fieldnames(ops0);

for j = 1:length(f)
    ops.(f{j}) = ops0.(f{j});
end

author: cortex-lab


text:
function  SubPopulation = SelectPop(SubPopulation,Problem)

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    for i = 1 : length(SubPopulation)
        pop = SubPopulation{i};
        [~,rank] = sort(pop.objs);
        SubPopulation{i} = pop(rank(1:Problem.N/2));
    end
end


text:
function  [value, isterminal, direction] = events_slidePos(~,Z,P)

[~, contacts] = dynamics_slidePos([],Z,P);
th = Z(1,:);
dx = Z(4,:);
H = contacts(1,:);
V = contacts(2,:);

%%%% CONSTRAINTS %%%%
% 1) -pi/2 < th
% 2) th < pi/2
% 3) dx >= 0
% 4) V >= 0
% 5) H <= 0

% As a rule, each constraint will be satisfied if it's event function value
% is positive. This makes things easier at the FSM level.
n = length(th);
value = zeros(5,n);
isterminal = true(size(value));
direction = -ones(size(value));

%%%% HACK %%%%
% avoid the singularity in the dynamics at th = pi/2
th_crit = pi/2 - 1e-10;
%%%% DONE %%%%

value(1,:) = th + th_crit ;
value(2,:) = th_crit - th;
value(3,:) = dx;
value(4,:) = V;
value(5,:) = -H;

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [alpha, b0] = qp_learn(K,Y,C)

  a=qpc(K,Y,-C*(Y==-1),C*(Y==1),zeros(length(K),1),1);
  a=a.*Y;
  alpha=a;

  H = K.*(Y*Y');

  class1 = find (a > max(a)/1e3 & a < 0.99*C & Y == 1);
  class2 = find (a > max(a)/1e3 & a < 0.99*C & Y == -1);
  nAsv = length (class1);
  nBsv = length (class2);
  if (nAsv == 0) & (nBsv==0)
  b0=(min(H(find(a < 0.99*C & Y== 1),:)*a)-      min(H(find(a < 0.99*C & Y==-1),:)*a))/2;
    warning ('Threshold might be inacurate');
 else

  A = sum(H(class1,:)*a)-nAsv;
  B = sum(H(class2,:)*a)-nBsv;
  b0 = -(A-B)/(nAsv+nBsv);
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function LSHparam = trainLSH(LSHparam)

% Input:
%          LSHparam
%              LSHparam.nbits---number of bits (nbits do not need to be a multiple of 8)
% Output:
%             LSHparam:
%                 LSHparam.w---random projection

dim = LSHparam.dim;
nbits = LSHparam.nbits;

W = randn(dim, nbits);

LSHparam.w = W;

author: willard-yuan

repo: hashing-baseline-for-image-retrieval


text:
function  [Z] =  MCWNNM_ADMM( Y, NSig, Par )
% This routine solves the following weighted nuclear norm optimization problem with column weights,
%
% min_{X, Z} ||W(Y-X)||_F^2 + ||Z||_w,*  s.t.  X = Z
%
% Inputs:
%        Y      -- 3p^2 x M dimensional noisy matrix, D is the data dimension, and N is the number of image patches.
%        NSig -- 3p^2 x 1 dimensional vector of weights
%        Par   -- structure of parameters
% Output:
%        Z      -- 3p^2 x M dimensional denoised matrix

% tol = 1e-8;
if ~isfield(Par, 'maxIter')
    Par.maxIter = 10;
end
if ~isfield(Par, 'rho')
    Par.rho = 1;
end
if ~isfield(Par, 'mu')
    Par.mu = 1;
end
if ~isfield(Par, 'display')
    Par.display = true;
end
%% Initializing optimization variables
% intialize
mNSig = min(NSig);
W = (mNSig+eps) ./ (NSig+eps);

Z = zeros(size(Y));
A = zeros(size(Y));
%% Start main loop
iter = 0;
PatNum       = size(Y,2);
TempC  = Par.Constant * sqrt(PatNum) * mNSig^2;
while iter < Par.maxIter
    iter = iter + 1;

    % update X, fix Z and A
    % min_{X} ||W * Y - W * X||_F^2 + 0.5 * rho * ||X - Z + 1/rho * A||_F^2
    X = diag(1 ./ (W.^2 + 0.5 * Par.rho)) * (diag(W.^2) * Y + 0.5 * Par.rho * Z - 0.5 * A);

    % update Z, fix X and A
    % min_{Z} ||Z||_*,w + 0.5 * rho * ||Z - (X + 1/rho * A)||_F^2
    Temp = X + A/Par.rho;
    [U, SigmaTemp, V] =   svd(full(Temp), 'econ');
    [SigmaZ, svp] = ClosedWNNM(diag(SigmaTemp), 2/Par.rho*TempC, eps);
    Z =  U(:, 1:svp) * diag(SigmaZ) * V(:, 1:svp)';
    %     % check the convergence conditions
    %     stopC = max(max(abs(X - Z)));
    %     if Par.display && (iter==1 || mod(iter,10)==0 || stopC<tol)
    %         disp(['iter ' num2str(iter) ',mu=' num2str(Par.mu,'%2.1e') ...
    %             ',rank=' num2str(rank(Z,1e-4*norm(Z,2))) ',stopALM=' num2str(stopC,'%2.3e')]);
    %     end
    %     if stopC < tol
    %         break;
    %     else
    % update the multiplier A, fix Z and X
    A = A + Par.rho * (X - Z);
    Par.rho = min(1e4, Par.mu * Par.rho);
    %     end
end
return;

author: csjunxu


text:
  function [xs, info] = l1_regress_admm1(yi, A, varargin)
%|function [xs, info] = l1_regress_admm1(yi, A, varargin)
%|
%| l1 regression, minimizing cost function
%| cost(x) = pot( y - A x ) where pot(r) = |r|_1 by default.
%| minimized via ADMM algorithm with v = A x split.
%|
%| in
%|	yi	[M 1]		noisy data
%|	A	[M N]		matrix
%|
%| options
%|	x0	[N 1]		initial estimate (default: A \ yi)
%|	shrink			shrink function for pot: shrink(x, reg)
%|				(default: l1, soft thresholding)
%|
%|	niter			# total iterations (default: 1)
%|					(max # if tol used)
%|	isave	[]		list of iterations to archive (default: 'last')
%|	userfun	@		user defined function handle (see default below)
%|					taking arguments (x, iter, userarg{:})
%|	userarg {}		user arguments to userfun (default {})
%|	stop_diff_tol		stop iterations if norm(xnew-xold)/norm(xnew)
%|				is less than this unitless value.  default: 0
%|	stop_diff_norm		use norm(., type) for stop rule
%|				choices: 1 | 2 (default) | inf
%|	chat	0|1		verbosity (default 0)
%|
%| out
%|	xs	[N niter]	estimates each iteration
%|	info	[niter 1]	time each iteration (for default userfun)
%|
%| Copyright 2013-03-21, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end
if nargin == 1 && streq(yi, 'test'), l1_regress_example, return, end

% defaults
arg.x0 = [];
arg.C = [];
arg.beta = 1; % regularization parameter
arg.shrink = [];
arg.rho = 1; % AL penalty parameter

arg.niter = 1;
arg.isave = [];
arg.userfun = @userfun_default;
arg.userarg = {};
arg.stop_diff_tol = 0;
arg.stop_diff_norm = 2;
arg.chat = 0;

arg = vararg_pair(arg, varargin);

if isempty(arg.x0)
	arg.x0 = A \ yi;
end

x = arg.x0;

if isempty(arg.shrink)
	tmp = potential_fun('l1', 1);
	shrink = @(z, reg) tmp.shrink(z, reg);
else
	shrink = arg.shrink;
end

arg.isave = iter_saver(arg.isave, arg.niter);
if arg.stop_diff_tol
	norm_diff = @(x) norm(x(:), arg.stop_diff_norm);
end

cpu etic

np = numel(x);
xs = zeros(np, length(arg.isave));
if any(arg.isave == 0)
	xs(:, arg.isave == 0) = x(:);
end

eta = 0; % dual variable

rho = arg.rho;

% iterate
for iter = 1:arg.niter
	ticker(mfilename, iter, arg.niter)

	xp = x; % previous

	% update primal
	x = A \ (yi + eta);

	% update auxiliary
	tmp = A * x - eta - yi;
	vv = yi + shrink(tmp, 1 / rho);

	% update multiplier
	tmp = A * x - vv;
	eta = eta - tmp;

	if any(arg.isave == iter)
		xs(:, arg.isave == iter) = x;
	end
	info(iter,:) = arg.userfun(x, iter, arg.userarg{:});

	% check norm(xnew-xold) / norm(xnew) vs threshold
	if iter > 1 && arg.stop_diff_tol
		ratio = norm_diff(x - xp) / norm_diff(x);
		if ratio < arg.stop_diff_tol
			if 1 || arg.chat
				printm('stop at iteration %d, diff %g < %g', ...
					iter, ratio, arg.stop_diff_tol)
			end
			if isequal(arg.isave, arg.niter) % saving only last?
				xs = x; % save the 'final' iterate
			else % saving many iterates?
				xs(:, arg.isave > iter) = []; % clear out unused
			end
		break
		end
	end
end

% default user function.
% using this evalin('caller', ...) trick, one can compute anything of interest
function out = userfun_default(x, iter, varargin)
%pr minmax(x)
out = [cpu('etoc')];

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function model = Expectation_MixHP(Seqs, model, alg)

Nk = sum(model.R,1); % 10.51
alpha = model.alpha+Nk; % Dirichlet
Elogpi = psi(0,alpha)-psi(0,sum(alpha)); % 10.66

% calculate responsibility
tic;
EX = zeros(length(Seqs), model.K);

% E-step: evaluate the responsibility using the current parameters
for c = 1:length(Seqs)

    Time = Seqs(c).Time;
    Event = Seqs(c).Mark;
    Tstart = Seqs(c).Start;

    if isempty(alg.Tmax)
        Tstop = Seqs(c).Stop;
    else
        Tstop = alg.Tmax;
        indt = Time < alg.Tmax;
        Time = Time(indt);
        Event = Event(indt);
    end

    N = length(Time);
    % calculate the integral decay function in the log-likelihood function
    G = Kernel_Integration(Tstop - Time, model);

    LL = Elogpi;

    for i = 1:N

        ui = Event(i);
        ti = Time(i);

        Elambdai = sqrt(pi/2).*model.b(ui,:)+eps;
        Vlambdai = (2-pi/2).*(model.b(ui,:)).^2;
        if i>1
            tj = Time(1:i-1);
            uj = Event(1:i-1);

            gij = Kernel(ti-tj, model);
            auiuj = model.beta(uj, :, :, ui);
            pij = repmat(gij, [1,1,model.K,1]).* auiuj;

            tmp = sum(sum(pij,1),2);
            Elambdai = Elambdai + tmp(:)';
            tmp = sum(sum(pij.^2, 1), 2);
            Vlambdai = Vlambdai + tmp(:)';
        end

        LL = LL+log(Elambdai) - Vlambdai./(2*Elambdai.^2);

    end
    LL = LL - (Tstop-Tstart).*sqrt(pi/2).*sum(model.b);
    tmp = sum(sum(repmat(G, [1,1,model.K]).*...
        sum(model.beta(Event,:,:,:),4),1),2);
    LL = LL- tmp(:)';

    XX = (LL - max(LL));
    EX(c,:)=(exp(XX))./sum(exp(XX));

end

model.R = EX;

end

author: HongtengXu

repo: Hawkes-Process-Toolkit


text:
clear all;
addpath('../lib');

%% lets plot 3 cycles of 50Hz AC voltage
f = 50;
Vm = 10;
phi = 0;

% generate the signal
t = [0:0.0001:3/f];
th = 2*pi*f*t;
v = Vm*sin(th+phi);

%% plot now
plt = Plot(t*1E3, v);

plt.XLabel = 'Time, t (ms)'; % xlabel
plt.YLabel = 'Voltage, V (V)'; %ylabel
plt.YTick = [-10, 0, 10]; %[tick1, tick2, .. ]
plt.YLim = [-11, 11]; % [min, max]
plt.XMinorTick = 'off'; % 'on' or 'off'
plt.Colors = {[1, 0, 0], [0, 0, 1]}; %[red, green, blue]

plt.FontName = 'Times';

% Save? comment the following line if you do not want to save
plt.export('plotVoltSomeOther.tiff');

author: masumhabib

repo: PlotPub

save_path: github-repos/MATLAB/masumhabib-PlotPub


text:
function m = gpuDistribute(m)
%  m = gpuDistribute(m)
%---------------------------------------------------------------------
% Distibutes the fields of the model <m> onto the GPU card.
%---------------------------------------------------------------------
% DES

author: dustinstansbury

repo: medal


text:
function [Zica, W, T, mu] = kICA(Z,r)
%
% Syntax:       Zica = kICA(Z,r);
%               [Zica, W, T, mu] = kICA(Z,r);
%
% Inputs:       Z is an d x n matrix containing n samples of d-dimensional
%               data
%
%               r is the number of independent components to compute
%
% Outputs:      Zica is an r x n matrix containing the r independent
%               components - scaled to variance 1 - of the input samples
%
%               W and T are the ICA transformation matrices such that
%               Zr = T \ W' * Zica + repmat(mu,1,n);
%               is the r-dimensional ICA approximation of Z
%
%               mu is the d x 1 sample mean of Z
%
% Description:  Performs independent component analysis (ICA) on the input
%               data using the max-kurtosis ICA algorithm
%
% Reference:    http://www.cs.nyu.edu/~roweis/kica.html
%
% Author:       Brian Moore
%               brimoor@umich.edu
%
% Date:         November 12, 2016
%

% Center and whiten data
mu = mean(Z,2);
T = sqrtm(inv(cov(Z')));
Zcw = T * bsxfun(@minus,Z,mu);

% Max-kurtosis ICA
[W, ~, ~] = svd(bsxfun(@times,sum(Zcw.^2,1),Zcw) * Zcw');
Zica = W(1:r,:) * Zcw;

author: xiuyechen

repo: FishExplorer


text:
% HKS_matching
clear all; close all; clc;

gptoolbox_path = ''; % complete
addpath(genpath(gptoolbox_path));

data_folder = '..\data\';
meshes = dir([data_folder,'*.o*']); % a list of all the meshes in the data folder
meshes = {meshes.name}; % access each mesh name by meshes{i}

cam_folder = [data_folder, 'cams\'];
cams = dir([cam_folder,'*.mat']); % a list of all the cam in the cams folder
cams = {cams.name};

addpath(data_folder); addpath(cam_folder);

%%
% % For dragon model:
% [V, F] = readOBJ(meshes{8});
% cam = cams{8};
% landmarks = [6921; 6827; 9683; 9710]; % right front leg, left front leg, right back leg, left back leg,
% figure; MESH_VIS.mesh(F,V,'cams',cam,'landmarks',landmarks);

% % For cat model:
[V, F] = readOFF(meshes{4});
cam = cams{4};
landmarks = [3048; 1994; 6621; 5560]; % right front leg, left front leg, right back leg, left back leg,
figure; MESH_VIS.mesh(F,V,'cams',cam,'landmarks',landmarks);

% % For xyzrgb_dragon1 model:
% [V, F] = readOBJ(meshes{14});
% cam = cams{8};
% landmarks = [91921; 58423; 51903; 52775]; % right front leg, left front leg, right back leg, left back leg,

%% Matching:
% to compute the HKS, use gptoolbox's: [K,MK] = hks(V,F)
% and to mormalize the hks use: K./MK

author: odedstein

repo: sgi-introduction-course