text:
% std_plottf() - plot ERSP/ITC images a component
%              or channel cluster in a STUDY. Also allows plotting scalp
%              maps.
% Usage:
%          >> std_plottf( times, freqs, data, 'key', 'val', ...)
% Inputs:
%  times - [vector] latencies in ms of the data points.
%  freqs - [vector] frequencies in Hz of the data points.
%  data  -  [cell array] mean data for each subject group and/or data
%           condition. For example, to plot mean ERPs from a STUDY
%           for epochs of 800 frames in two conditions from three groups
%           of 12 subjects:
%
%           >> data = { [800x12] [800x12] [800x12];... % 3 groups, cond 1
%                       [800x12] [800x12] [800x12] };  % 3 groups, cond 2
%           >> std_plottf(erp_ms,data);
%
%           By default, parametric statistics are computed across subjects
%           in the three groups. (group,condition) ERP averages are plotted.
%           See below and >> help statcond
%           for more information about the statistical computations.
%
% Optional display parameters:
%  'datatype'    - ['ersp'|'itc'] data type {default: 'ersp'}
%  'titles'      - [cell array of string] titles for each of the subplots.
%                  { default: none}
%
% Statistics options:
%  'groupstats'  - ['on'|'off'] Compute (or not) statistics across groups.
%                  {default: 'off'}
%  'condstats'   - ['on'|'off'] Compute (or not) statistics across groups.
%                  {default: 'off'}

%  'threshold'   - [NaN|real<<1] Significance threshold. NaN -> plot the
%                  p-values themselves on a different figure. When possible,
%                  significance regions are indicated below the data.
%                  {default: NaN}
%  'maskdata'    - ['on'|'off'] when threshold is non-NaN and not both
%                  condition and group statistics are computed, the user
%                  has the option to mask the data for significance.
%                  {defualt: 'off'}
%
% Other plotting options:
%  'plotmode'    - ['normal'|'condensed'] statistics plotting mode:
%                  'condensed' -> plot statistics under the curves
%                  (when possible); 'normal' -> plot them in separate
%                  axes {default: 'normal'}
%  'freqscale'   - ['log'|'linear'|'auto'] frequency plotting scale.
%                  {default: 'auto'}
%  'ylim'        - [min max] ordinate limits for ERP and spectrum plots
%                  {default: all available data}
%
% ITC/ERSP image plotting options:
%  'tftopoopt'   - [cell array] tftopo() plotting options (ERSP and ITC)
%  'caxis'       - [min max] color axis (ERSP, ITC, scalp maps)
%
% Scalp map plotting options:
%  'chanlocs'    - [struct] channel location structure
%
% Author: Arnaud Delorme, CERCO, CNRS, 2006-
%
% See also: pop_erspparams(), pop_erpparams(), pop_specparams(), statcond()

% Copyright (C) 2006 Arnaud Delorme
%
% 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [pgroup, pcond, pinter] = std_plottf(timevals, freqs, data, varargin)

pgroup = [];
pcond  = [];
pinter = [];
if nargin < 2
    help std_plottf;
    return;
end;

opt = finputcheck( varargin, { 'titles'         'cell'   []              cell(20,20);
                               'caxis'          'real'   []              [];
                               'ersplim'        'real'   []              []; % same as above
                               'itclim'         'real'   []              []; % same as above
                               'ylim'           'real'   []              [];
                               'tftopoopt'      'cell'   []              {};
                               'threshold'      'real'   []              NaN;
                               'unitx'          'string' []              'ms'; % just for titles
                               'chanlocs'       'struct' []              struct('labels', {});
                               'freqscale'      'string' { 'log' 'linear' 'auto' }  'auto';
                               'groupstats'     'cell'   []              {};
                               'condstats'      'cell'   []              {};
                               'interstats'     'cell'   []              {};
                               'maskdata'       'string' { 'on' 'off' }   'off';
                               'datatype'       'string' { 'ersp' 'itc' }    'ersp';
                               'plotmode'       'string' { 'normal' 'condensed' }  'normal' }, 'std_plottf');
if isstr(opt), error(opt); end;
if all(all(cellfun('size', data, 3)==1))               opt.singlesubject = 'on'; end;

% remove empty entries
datapresent = ~cellfun(@isempty, data);
for c = size(data,1):-1:1, if sum(datapresent(c,:)) == 0, data(c,:) = []; opt.titles(c,:) = []; if ~isempty(opt.groupstats), opt.groupstats(c) = []; end; end; end;
for g = size(data,2):-1:1, if sum(datapresent(:,g)) == 0, data(:,g) = []; opt.titles(:,g) = []; if ~isempty(opt.condstats ), opt.condstats( g) = []; end; end; end;

if ~isempty(opt.groupstats) & ~isempty(opt.condstats) & strcmpi(opt.maskdata, 'on')
    disp('Cannot use ''maskdata'' option with both condition stat. and group stat. on');
    disp('Disabling statistics');
    opt.groupstats = {}; opt.condstats = {}; opt.maskdata = 'off';
end;
if ~isempty(opt.ersplim), opt.caxis = opt.ersplim; end;
if ~isempty(opt.itclim), opt.caxis = opt.itclim; end;
onecol  = { 'b' 'b' 'b' 'b' 'b' 'b' 'b' 'b' 'b' 'b' };
manycol = { 'b' 'r' 'g' 'k' 'c' 'y' };

nc = size(data,1);
ng = size(data,2);
if nc >= ng, opt.transpose = 'on';
else         opt.transpose = 'off';
end;

% test log frequencies
% --------------------
if length(freqs) > 2 & strcmpi(opt.freqscale, 'auto')
    midfreq = (freqs(3)+freqs(1))/2;
    if midfreq*.9999 < freqs(2) & midfreq*1.0001 > freqs(2), opt.freqscale = 'linear';
    else                                                     opt.freqscale = 'log';
    end;
end;

% condensed plot
% --------------
if strcmpi(opt.plotmode, 'condensed')
    meanplot = zeros(size(data{1},1), size(data{1},2));
    count = 0;
    for c = 1:nc
        for g = 1:ng
            if ~isempty(data{c,g})
                meanplot = meanplot + mean(data{c,g},3);
                count = count+1;
            end;
        end;
    end;
    meanplot = meanplot/count;
    options = { 'chanlocs', opt.chanlocs, 'electrodes', 'off', 'cbar', 'on', ...
            'cmode', 'separate', opt.tftopoopt{:} };

    if strcmpi(opt.freqscale, 'log'), options = { options{:} 'logfreq', 'native' }; end;
    tftopo( meanplot', timevals, freqs, 'title', opt.titles{1}, options{:});
    currentHangle = gca;
    if ~isempty( opt.caxis )
        caxis( currentHangle, opt.caxis )
    end
    colorbarHandle = cbar;
    title(colorbarHandle,'dB');
    axes(currentHangle);
    return;
end;

% plotting paramters
% ------------------
if ng > 1 && ~isempty(opt.groupstats), addc = 1; else addc = 0; end;
if nc > 1 && ~isempty(opt.condstats  ), addr = 1; else addr = 0; end;

% compute significance mask
% --------------------------
if ~isempty(opt.interstats), pinter = opt.interstats{3}; end;

if ~isnan(opt.threshold) && ( ~isempty(opt.groupstats) || ~isempty(opt.condstats) )
    pcondplot  = opt.condstats;
    pgroupplot = opt.groupstats;
    pinterplot = pinter;
    maxplot = 1;
else
    warning off;
    for ind = 1:length(opt.condstats),  pcondplot{ind}  = -log10(opt.condstats{ind}); end;
    for ind = 1:length(opt.groupstats), pgroupplot{ind} = -log10(opt.groupstats{ind}); end;
    if ~isempty(pinter), pinterplot = -log10(pinter); end;
    maxplot = 3;
    warning on;
end;

% -------------------------------
% masking for significance of not
% -------------------------------
statmask = 0;
if strcmpi(opt.maskdata, 'on') && ~isnan(opt.threshold) && ...
        (~isempty(opt.condstats) || ~isempty(opt.condstats))
    addc = 0; addr = 0; statmask = 1;
end;

% -------------------------
% plot time/frequency image
% -------------------------
options = { 'chanlocs', opt.chanlocs, 'electrodes', 'off', 'cbar', 'off', ...
            'cmode', 'separate', opt.tftopoopt{:} };
if strcmpi(opt.freqscale, 'log'), options = { options{:} 'logfreq', 'native' }; end;

% adjust figure size
% ------------------
fig = figure('color', 'w');
pos = get(fig, 'position');
set(fig, 'position', [ pos(1)+15 pos(2)+15 pos(3)/2.5*(nc+addr), pos(4)/2*(ng+addc) ]);
pos = get(fig, 'position');
if strcmpi(opt.transpose, 'off'), set(gcf, 'position', [ pos(1) pos(2) pos(4) pos(3)]);
else                              set(gcf, 'position', pos);
end;

tmpc = [inf -inf];
for c = 1:nc
    for g = 1:ng
        hdl(c,g) = mysubplot(nc+addr, ng+addc, g + (c-1)*(ng+addc), opt.transpose);
        if ~isempty(data{c,g})
            tmpplot = mean(data{c,g},3);
            if ~isreal(tmpplot(1)), tmpplot = abs(tmpplot); end;
            if statmask,
                if ~isempty(opt.condstats), tmpplot(find(pcondplot{g}(:) == 0)) = 0;
                else                        tmpplot(find(pgroupplot{c}(:) == 0)) = 0;
                end;
            end;

            tftopo( tmpplot', timevals, freqs, 'title', opt.titles{c,g}, options{:});

            if isempty(opt.caxis) && ~isempty(tmpc)
                warning off;
                tmpc = [ min(min(tmpplot(:)), tmpc(1)) max(max(tmpplot(:)), tmpc(2)) ];
                warning on;
            else
                if ~isempty(opt.caxis)
                    caxis(opt.caxis);
                end;
            end;

            if c > 1
                ylabel('');
            end;
        end;

        % statistics accross groups
        % -------------------------
        if g == ng && ng > 1 && ~isempty(opt.groupstats) && ~isinf(pgroupplot{c}(1)) && ~statmask
            hdl(c,g+1) = mysubplot(nc+addr, ng+addc, g + 1 + (c-1)*(ng+addc), opt.transpose);
            tftopo( pgroupplot{c}', timevals, freqs, 'title', opt.titles{c,g+1}, options{:});
            caxis([-maxplot maxplot]);
        end;
    end;
end;

for g = 1:ng
    % statistics accross conditions
    % -----------------------------
    if ~isempty(opt.condstats) && ~isinf(pcondplot{g}(1)) && ~statmask && nc > 1
        hdl(nc+1,g) = mysubplot(nc+addr, ng+addc, g + c*(ng+addc), opt.transpose);
        tftopo( pcondplot{g}', timevals, freqs, 'title', opt.titles{nc+1,g}, options{:});
        caxis([-maxplot maxplot]);
    end;
end;

% color scale
% -----------
if isempty(opt.caxis)
    tmpc = [-max(abs(tmpc)) max(abs(tmpc))];
    for c = 1:nc
        for g = 1:ng
            axes(hdl(c,g));
            if ~isempty(tmpc)
                caxis(tmpc);
            end;
        end;
    end;
end;

% statistics accross group and conditions
% ---------------------------------------
if ~isempty(opt.groupstats) && ~isempty(opt.condstats) && ng > 1 && nc > 1
    hdl(nc+1,ng+1) = mysubplot(nc+addr, ng+addc, g + 1 + c*(ng+addr), opt.transpose);
    tftopo( pinterplot',  timevals, freqs, 'title', opt.titles{nc+1,ng+1}, options{:});
    caxis([-maxplot maxplot]);
    ylabel('');
end;

% color bars
% ----------
axes(hdl(nc,ng));
cbar_standard(opt.datatype, ng);
if isnan(opt.threshold) && (nc ~= size(hdl,1) || ng ~= size(hdl,2))
    ind = find(hdl(end:-1:1));
    axes(hdl(end-ind(1)+1));
    cbar_signif(ng, maxplot);
end;

% mysubplot (allow to transpose if necessary)
% -------------------------------------------
function hdl = mysubplot(nr,nc,ind,transp);

    r = ceil(ind/nc);
    c = ind -(r-1)*nc;
    if strcmpi(transp, 'on'), hdl = subplot(nc,nr,(c-1)*nr+r);
    else                      hdl = subplot(nr,nc,(r-1)*nc+c);
    end;

% colorbar for ERSP and scalp plot
% --------------------------------
function cbar_standard(datatype, ng);
    pos = get(gca, 'position');
    tmpc = caxis;
    fact = fastif(ng == 1, 40, 20);
    tmp = axes('position', [ pos(1)+pos(3)+max(pos(3)/fact,0.006) pos(2) max(pos(3)/fact,0.01) pos(4) ]);
    set(gca, 'unit', 'normalized');
    if strcmpi(datatype, 'itc')
         cbar(tmp, 0, tmpc, 10); ylim([0.5 1]);
         title('ITC');
    else cbar(tmp, 0, tmpc, 5);title('dB');
    end;

% colorbar for significance
% -------------------------
function cbar_signif(ng, maxplot);
    pos = get(gca, 'position');
    tmpc = caxis;
    fact = fastif(ng == 1, 40, 20);
    tmp = axes('position', [ pos(1)+pos(3)+max(pos(3)/fact,0.006) pos(2) max(pos(3)/fact,0.01) pos(4) ]);
    map = colormap;
    n = size(map,1);
    cols = [ceil(n/2):n]';
    image([0 1],linspace(0,maxplot,length(cols)),[cols cols]);
    %cbar(tmp, 0, tmpc, 5);
    tick = linspace(0, maxplot, maxplot+1);
    set(gca, 'ytickmode', 'manual', 'YAxisLocation', 'right', 'xtick', [], ...
        'ytick', tick, 'yticklabel', round(10.^-tick*1000)/1000);
    xlabel('');

% rapid filtering for ERP
% -----------------------
function tmpdata2 = myfilt(tmpdata, lowpass, highpass, factor, filtertype)

    tmpdata2 = reshape(tmpdata, size(tmpdata,1), size(tmpdata,2)*size(tmpdata,3)*size(tmpdata,4));
    tmpdata2 = eegfiltfft(tmpdata2',lowpass, highpass, factor, filtertype)';
    tmpdata2 = reshape(tmpdata2, size(tmpdata,1), size(tmpdata,2), size(tmpdata,3), size(tmpdata,4));

author: PatternRecognition

repo: OpenBMI


text:
 function y = mtimes(ob, x)
%function y = mtimes(ob, x)	for Gnufft object
% y = G * x	or x = G' * y

if ob.apower ~= 1, error 'power not done', end

%
% partial projection or backprojection (for ordered subsets)
%
if ob.is_subref
	error 'subref not done'

%
% full "projection"
%
elseif ~ob.is_transpose
	if ob.is_masked
		x = embed(x, ob.mask);
	end

	%
	% full dsft, which doesn't have shift implemented anyway
	%
%	if ob.is.dsft
%		x = reshape(x, ob.Nd);
%		y = dtft_mex('forward', ob.omega', x, int32(ob.nthread));

	%
	% nufft
	%
%	else
		x = reshape(x, [ob.Nd numel(x)/prod(ob.Nd)]);
		y = nufft(x, ob.st);
%	end

%
% full "back-projection"
%
else

%	if ob.is.dsft
%		if isreal(y)
%			y = complexify(y);
%			warning 'faking complex'
%		end
%		y = dtft_mex('adjoint', ob.omega', y, ...
%			int32([ob.nx ob.ny]'), int32(ob.nthread));
%		y = ob.pixel_size * y;

%	else
		y = nufft_adj(x, ob.st);
%	end

	if ob.is_masked
		y = y(ob.mask);
	end

	y = reshape(y, [ob.dims(1) numel(y) / ob.dims(1)]); % [*Nd,L]
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function plotx2east
% set the default plot direction of the x-axis

setMTEXpref('xAxisDirection','east');

author: mtex-toolbox

repo: mtex

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


text:
function [M,S]=rigidbody_transform(X)

% function [M]=rigidbody_transform(X)
% ------------------------------------------------------------------------
%
% Kevin Mattheus Moerman
% kevinmoerman@hotmail.com
% 22/04/2011
% ------------------------------------------------------------------------

%% Determine translation matrix
OR=mean(X,1);

%Defining translation matrix
T  = [1 0 0 OR(1);...
      0 1 0 OR(2);...
      0 0 1 OR(3);...
      0 0 0 1];

%% Determine rotation matrix

X=[X(:,1)-OR(1) X(:,2)-OR(2) X(:,3)-OR(3)]; %Centre points around mean
[~,S,V]=svd(X,0); %Singular value decomposition

%Defining direction cosine matrix

if 1-V(3,3)<eps('double')
    DCM=eye(3,3);
else
    rz=V(:,3); rz=rz./sqrt(sum(rz.^2)); %surface normal
    r=V(:,2); r=r./sqrt(sum(r.^2));
    rx=cross(rz,r);rx=rx./sqrt(sum(rx.^2));
    ry=cross(rx,rz);ry=ry./sqrt(sum(ry.^2));
    DCM=[rx(:) ry(:) rz(:)];
end

% N=-V(:,3)./V(3,3); %Surface normal
% rx=[1 0 N(1)]; rx=rx./sqrt(sum(rx.^2));
% ry=[0 1 N(2)]; ry=ry./sqrt(sum(ry.^2));
% rz=cross(rx,ry); rz=rz./sqrt(sum(rz.^2));
% DCM=[rx(:) ry(:) rz(:)];

R  = eye(4,4);
R(1:3,1:3)=DCM;

%% Create translation rotation matrix
M  = T * R * eye(4,4);

end

%%
% _*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:
classdef TitleColorable < handle
    % Mixin for component with Title font color

    % Copyright 2020-2021 The MathWorks Inc.

    %% Properties
    properties (AbortSet)

        % Title color
        TitleColor (1,3) double {wt.validators.mustBeBetweenZeroAndOne} = [0 0 0]

    end %properties

    %% Internal properties
    properties (AbortSet, Transient, NonCopyable, ...
            Access = {?wt.abstract.BaseWidget, ?wt.test.BaseWidgetTest} )

        % List of graphics controls to apply to
        TitleColorableComponents (:,1) matlab.graphics.Graphics

    end %properties

    %% Accessors
    methods

        function set.TitleColor(obj,value)
            obj.TitleColor = value;
            obj.updateTitleColorableComponents()
        end

        function set.TitleColorableComponents(obj,value)
            obj.TitleColorableComponents = value;
            obj.updateTitleColorableComponents()
        end

    end %methods

    %% Methods
    methods (Access = protected)

        function updateTitleColorableComponents(obj)

            wt.utility.fastSet(obj.TitleColorableComponents,"FontColor",obj.TitleColor)

        end %function

    end %methods

end %classdef

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%DEMO_HURDLE Demonstration of Logit Negative-binomial hurdle model
%            using Gaussian process prior
%
%  Description
%    Hurdle models can be used to model excess number of zeros
%    compared to usual Poisson and negative binomial count models.
%    Hurdle models assume a two-stage process, where the first
%    process determines whether the count is larger than zero, and
%    the second process determines the non-zero count. Both
%    processes are given a zero mean Gaussian process prior. The
%    two stage model formulation makes it possible to make
%    inference for the two latent processes separately.
%
%    In this demo we construct logit negative binomial hurdle model
%    by using logit and zero truncated negative binomial models.
%    The posterior inference is made with parallel-EP
%    approximation.
%
%    See also  DEMO_SPATIAL2, DEMO_CLASSIFIC1
%
% Copyright (c) 2008-2010 Jarno Vanhatalo
% Copyright (c) 2010-2011 Aki Vehtari
% Copyright (c) 2011 Jaakko Riihimäki

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

% load the data
S = which('demo_spatial1');
data = load(strrep(S,'demo_spatial1.m','demodata/spatial1.txt'));

x = data(:,1:2);
ye = data(:,3);
y = data(:,4);

x0=x;
x=bsxfun(@rdivide,bsxfun(@minus,x0,mean(x0)),std(x0));

% for zero-process
yz=double(y>0)*2-1;
% for count-process
ci=find(y>0);
yc=y(ci);
xc=x(ci,:);
yec=ye(ci,:);

% Create the covariance functions
pl = prior_t('s2',10);
pm = prior_sqrtunif();
cf = gpcf_matern32('lengthScale', 5, 'magnSigma2', 0.05, ...
                   'lengthScale_prior', pl, 'magnSigma2_prior', pm);

% Create the zero part
likz=lik_probit();
gpz=gp_set('lik',likz,'cf',cf,'jitterSigma2',1e-6,'latent_method','EP','latent_opt',struct('parallel','on'));
% Create the count part
likc=lik_negbinztr();
gpc=gp_set('lik',likc,'cf',cf,'jitterSigma2',1e-6,'latent_method','EP','latent_opt',struct('parallel','on'));
% Note that although both parts use 'cf', they have separate hyperparameters

% Set the options for the quasi-Newton optimization
opt=optimset('TolFun',1e-2,'TolX',1e-2,'Display','iter');
gpz=gp_optim(gpz,x,yz,'opt',opt,'optimf',@fminlbfgs);
gpc=gp_optim(gpc,xc,yc,'z',yec,'opt',opt,'optimf',@fminlbfgs);

% make prediction to the data points
[Efz, Varfz] = gp_pred(gpz, x, yz, x);
[Efc, Varfc] = gp_pred(gpc, xc, yc, xc, 'z', yec);

% Define help parameters for plotting
xii=sub2ind([60 35],x0(:,2),x0(:,1));
[X1,X2]=meshgrid(1:35,1:60);

% Plot the figures
figure
subplot(1,2,1)
G=NaN(size(X1));
G(xii)=Efz(:);
pcolor(X1,X2,G),shading flat
colorbar
axis equal
axis([0 35 0 60])
title('Posterior mean of latent zero process')

subplot(1,2,2)
G=NaN(size(X1));
G(xii(ci))=Efc(:);
pcolor(X1,X2,G),shading flat
colorbar
axis equal
axis([0 35 0 60])
title('Posterior mean of latent count process')

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function saveOBBs(obblist,filename)

fid = fopen(filename,'wt');

for i = 1:size(obblist,1)
    p = obblist(i,:);

    center = p(1:3);
    dir_1 = p(4:6);
    dir_2 = p(7:9);
    lengths = p(10:12);

    dir_1 = dir_1/norm(dir_1);
    dir_2 = dir_2/norm(dir_2);
    dir_3 = cross(dir_1,dir_2);
    dir_3 = dir_3/norm(dir_3);
    cornerpoints = zeros(8,3);

    d1 = 0.5*lengths(1)*dir_1;
    d2 = 0.5*lengths(2)*dir_2;
    d3 = 0.5*lengths(3)*dir_3;
    cornerpoints(1,:) = center-d1-d2-d3;
    cornerpoints(2,:) = center-d1+d2-d3;
    cornerpoints(3,:) = center+d1-d2-d3;
    cornerpoints(4,:) = center+d1+d2-d3;
    cornerpoints(5,:) = center-d1-d2+d3;
    cornerpoints(6,:) = center-d1+d2+d3;
    cornerpoints(7,:) = center+d1-d2+d3;
    cornerpoints(8,:) = center+d1+d2+d3;

    for j = 1:size(cornerpoints,1)
        fprintf(fid,'v  %f  %f  %f  \n',cornerpoints(j,1),cornerpoints(j,2),cornerpoints(j,3));
    end
end

N = 8;
for i = 1:size(obblist,1)
    fprintf(fid,'f  %d  %d  %d\n',(i-1)*N+6,(i-1)*N+4,(i-1)*N+2);
    fprintf(fid,'f  %d  %d  %d\n',(i-1)*N+6,(i-1)*N+8,(i-1)*N+4);
    fprintf(fid,'f  %d  %d  %d\n',(i-1)*N+8,(i-1)*N+7,(i-1)*N+4);
    fprintf(fid,'f  %d  %d  %d\n',(i-1)*N+7,(i-1)*N+3,(i-1)*N+4);
    fprintf(fid,'f  %d  %d  %d\n',(i-1)*N+7,(i-1)*N+5,(i-1)*N+3);
    fprintf(fid,'f  %d  %d  %d\n',(i-1)*N+5,(i-1)*N+2,(i-1)*N+1);
    fprintf(fid,'f  %d  %d  %d\n',(i-1)*N+5,(i-1)*N+6,(i-1)*N+2);
    fprintf(fid,'f  %d  %d  %d\n',(i-1)*N+4,(i-1)*N+3,(i-1)*N+1);
    fprintf(fid,'f  %d  %d  %d\n',(i-1)*N+2,(i-1)*N+4,(i-1)*N+1);
    fprintf(fid,'f  %d  %d  %d\n',(i-1)*N+5,(i-1)*N+8,(i-1)*N+6);
    fprintf(fid,'f  %d  %d  %d\n',(i-1)*N+5,(i-1)*N+7,(i-1)*N+8);
end

fclose(fid);
% plot3([cornerpoints(1,1),cornerpoints(2,1)],[cornerpoints(1,2),cornerpoints(2,2)],[cornerpoints(1,3),cornerpoints(2,3)],col);hold on;
% plot3([cornerpoints(1,1),cornerpoints(3,1)],[cornerpoints(1,2),cornerpoints(3,2)],[cornerpoints(1,3),cornerpoints(3,3)],col);hold on;
% plot3([cornerpoints(2,1),cornerpoints(4,1)],[cornerpoints(2,2),cornerpoints(4,2)],[cornerpoints(2,3),cornerpoints(4,3)],col);hold on;
% plot3([cornerpoints(3,1),cornerpoints(4,1)],[cornerpoints(3,2),cornerpoints(4,2)],[cornerpoints(3,3),cornerpoints(4,3)],col);hold on;
% plot3([cornerpoints(5,1),cornerpoints(6,1)],[cornerpoints(5,2),cornerpoints(6,2)],[cornerpoints(5,3),cornerpoints(6,3)],col);hold on;
% plot3([cornerpoints(5,1),cornerpoints(7,1)],[cornerpoints(5,2),cornerpoints(7,2)],[cornerpoints(5,3),cornerpoints(7,3)],col);hold on;
% plot3([cornerpoints(6,1),cornerpoints(8,1)],[cornerpoints(6,2),cornerpoints(8,2)],[cornerpoints(6,3),cornerpoints(8,3)],col);hold on;
% plot3([cornerpoints(7,1),cornerpoints(8,1)],[cornerpoints(7,2),cornerpoints(8,2)],[cornerpoints(7,3),cornerpoints(8,3)],col);hold on;
% plot3([cornerpoints(1,1),cornerpoints(5,1)],[cornerpoints(1,2),cornerpoints(5,2)],[cornerpoints(1,3),cornerpoints(5,3)],col);hold on;
% plot3([cornerpoints(2,1),cornerpoints(6,1)],[cornerpoints(2,2),cornerpoints(6,2)],[cornerpoints(2,3),cornerpoints(6,3)],col);hold on;
% plot3([cornerpoints(3,1),cornerpoints(7,1)],[cornerpoints(3,2),cornerpoints(7,2)],[cornerpoints(3,3),cornerpoints(7,3)],col);hold on;
% plot3([cornerpoints(4,1),cornerpoints(8,1)],[cornerpoints(4,2),cornerpoints(8,2)],[cornerpoints(4,3),cornerpoints(8,3)],col);hold on;

end

repo: GRAINS


text:
function th_img = mmse7(inp,inp_map,s,n)
[row col]=size(inp);
th_img=zeros(row,col);
th_img(1:(row/2^n),1:(col/2^n))=inp(1:(row/2^n),1:(col/2^n));
for i=n:-1:1
    a=inp((row/2^i)+1:(row/2^(i-1)),(col/2^i)+1:(col/2^(i-1))); %diagonal coefficients
    b=inp_map((row/2^i)+1:(row/2^(i-1)),(col/2^i)+1:(col/2^(i-1)));
    th_img((row/2^i)+1:(row/2^(i-1)),(col/2^i)+1:(col/2^(i-1)))=th(a,b,s);
    a=inp(1:(row/2^(i)),(col/2^i)+1:(col/2^(i-1)));   %horizontal coefficients
    b=inp(1:(row/2^(i)),(col/2^i)+1:(col/2^(i-1)));
    th_img(1:(row/2^(i)),(col/2^i)+1:(col/2^(i-1)))=th(a,b,s);
    a=inp((row/2^i)+1:(row/2^(i-1)),1:(col/2^(i)));     %vertical coefficients
    b=inp((row/2^i)+1:(row/2^(i-1)),1:(col/2^(i)));
    th_img((row/2^i)+1:(row/2^(i-1)),1:(col/2^(i)))=th(a,b,s);
end
end
function y=th(inp,inp1,s)
[row col]=size(inp);
lambda=1/var(inp1(:));
test1=zeros(row+6,col+6);
test2=zeros(size(test1));
test1(4:end-3,4:end-3)=inp;
[row col]=size(test1);
for i=4:row-3
    for j=4:col-3
        b=test1(i-3:i+3,j-3:j+3);
%         sigsq=max(0,((1/numel(b))*sum(sum(b.^2))-s));
%         lambda=1/sqrt(sigsq);
        A=numel(b)*(-1+sqrt(1+(8*lambda/numel(b)^2)*sum(sum(b.^2))));
        sig_cap_sq=max(0,(A/(4*lambda)-s));
        %test2(i,j)=(sqrt(2)*s)./sqrt(sigsq);
        test2(i,j)=(sig_cap_sq.*test1(i,j))./(sig_cap_sq+s);
    end
end
y=test2(4:row-3,4:col-3);
% % c=(abs(inp)-test2(4:row-3,4:col-3));
% % c(c<0)=0;
% % y=sign(inp).*c;
% [row col]=size(inp);
% y(1:row/2^n,1:col/2^n)=inp(1:row/2^n,1:col/2^n);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function r=rotpl2ro(u,v,t)
%ROTPL2RO find matrix to rotate in the plane containing u and v r=[u,v,t]
% Inputs:
%
%     U(n,1) and V(n,1) define a plane in n-dimensional space
%     T is the rotation angle in radians from U towards V. If T
%       is omitted it will default to the angle between U and V
%
% Outputs:
%
%     R(n,n)   Rotation matrix

%
%      Copyright (C) Mike Brookes 2007
%      Version: $Id: rotpl2ro.m,v 1.1 2007/11/21 16:50:32 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

u=u(:);
    n=length(u);
v=v(:);
l=sqrt(u'*u);
if l==0, error('input u is a zero vector'); end
u=u/l;      % normalize
q=v-v'*u*u;        % q is orthogonal to x
l=sqrt(q'*q);
if l==0          % u and v are colinear or v=zero
    [m,i]=max(abs(u));
    q=zeros(n,1);
    q(1+mod(i(1),n))=1;  % choose next available dimension
    q=q-q'*u*u;  % q is orthogonal to x
    l=sqrt(q'*q);
end
q=q/l;          % normalize
if nargin<3
    [s,c]=atan2sc(v'*q,v'*u);
    r=eye(n)+(c-1)*(u*u'+q*q')+s*(q*u'-u*q');
else
    r=eye(n)+(cos(t)-1)*(u*u'+q*q')+sin(t)*(q*u'-u*q');
end

author: decouples


text:
%  Figure 3.4     Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%script to generate Fig. 3.4
%% fig3_04.m
%Satellite Pulse response Example. 3.16
clf;
numG=[0 0 0.0002];
denG=[1 0 0];
sysG=tf(numG,denG);
t=0:0.01:10;
%pulse input
u1=[zeros(1,500) 25*ones(1,10) zeros(1,491)];
[y1]=lsim(sysG,u1,t);
plot(t,u1);
grid;
xlabel('Time (sec)');
ylabel('Thrust Fc');
title('Fig. 3.4(a): Thrust input')
axis([0 10 0 26]);
pause;
% conversion to degrees
ff=180/pi;
y1=ff*y1;
plot(t,y1);
grid;
xlabel('Time (sec)');
ylabel('\theta (deg)');
title('Fig. 3.4(b): Satellite attitude')

author: Sable

repo: mcbench-benchmarks


text:
function geometry_test03323 ( )

%*****************************************************************************80
%
%% TEST0323 tests I4COL_FIND_PAIR_WRAP.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  m = 5;
  n = 4;
  test_num = 5;

  item1_test = [ 22, 32, 22, 54, 54 ];
  item2_test = [ 32, 22, 23, 14, 11 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0323\n' );
  fprintf ( 1, '  I4COL_FIND_PAIR_WRAP finds the first occurrence of\n' );
  fprintf ( 1, '  a pair of item in an I4COL.\n' );
  fprintf ( 1, '  Items in the array are ordered by column, and\n' );
  fprintf ( 1, '  wraparound is allowed.\n' );

  for i = 1 : m
    for j = 1 : n
      a(i,j) = 10 * i + j;
    end
  end

  i4mat_print ( m, n, a, '  The matrix of columns:' );

  for test = 1 : test_num

    item1 = item1_test(test);
    item2 = item2_test(test);

    [ row, col ] = i4col_find_pair_wrap ( m, n, a, item1, item2 );

    fprintf ( 1, ...
      '  Item %d followed by item %d occurs in row %d and column %d\n', ...
      item1, item2, row, col );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [Vo, Ro, No, V] = generate_syntheticdata(F, N, K, noise_level, rho)
% Function to generate synthetic dataset
%
% This file is part of NMFLibrary.
%
% Created by H.Kasai
%
% Change log:
%
%

%
%     % Vo
%     W = randn(F,K);
%     W = max(W, 0);
%     W = min(W, 1);
%
%     H = randn(K,N);
%     H = max(H, 0);
%     H = min(H, 1);
%
%     Vo = W * H;
%
%     % noise
     No = noise_level*randn(F,N);
%
%     % add outlier
%     [V, Ro] = add_outlier(rho, F, N, Vo);
%
%     % V
%     V = V + No;
%     %V = max(V, 0);
%     %V = min(V,1);

    %
    sigma2 = 1 / sqrt(K);
    HN = makedist('Normal', 'mu', 0, 'sigma', sqrt(sigma2));
    Vo_n = random(HN, F, N);
    Vo = abs(Vo_n) ;
    Vo = min(Vo, 1);

    nu = rho;
    nu_tilda = 0.1;
    I = nu * N;
    card = nu_tilda * F;
    Ro = zeros(F,N);
    if rho > 0
        for i = 1 : N
            n_before = 0;
            if i < I
                for f = 1 : card
                    c = randi(F);
                     Ro(c,i) =  1 + (1+1)*rand(1, 1);
                     n = nnz(Ro(:,i));
                     if n_before == n
                         f = f - 1;
                     end
                     n_before = n;
                end

            end
        end
    end

    % V
    V = Vo + Ro + No;
    %V = max(V, 0);
    %V = min(V,1);

    index = find(V<0);
    V(index) = 0;

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function [curCosts] = scaleCosts(curCosts, counts)
  curCosts.total = curCosts.total/counts.total;
  curCosts.word = curCosts.word/counts.word;
end

author: lmthang


text:
function p = perms1(n,b)

% PERMS1(n,b), where n is a positive integer, creates a vector with n
% columns containing one permutation of the numbers from 1 to n. Input b is
% a number between 1 and n! which is the index of the required permutation,
% such that the output of PERMS1(n,b) for b from 1 to n! is equivalent to
% the output from perms(1:n). The order of the permutations is, however,
% different (e.g.,PERMS1(n,5) does not give the fifth row of the PERMS(1:n) output).

% The algorithm is based on imagining that the permutations are obtained by
% using n nested loops. The counter in each inner loop skips the values of all
% outer loops to generate the required permutations. A counter (index) in the
% innermost loop would count n! permutations. Working backwards, given the value
% of an index between 1 and n! this algorithm figures out the state of the
% counter in each these imaginary loops, which are actually the required
% permutation indices.

%K. H. Hamed,
%Sultan Qaboos University, Muscat, Oman
%khhamed@squ.edu.om
%21 November 2006

if nargin <2
    error('PERMS1:minrhs','Two Input Arguments are required.')
end;
if fix(n) ~= n || fix(b) ~= b ||n < 0 || b < 0 || ~isa(n,'double') || ~isa(b,'double') || ~isreal(n) || ~isreal(b)
  error('PERMS1:nnegint','Input Parameters n and b Must be Non-negative Integers');
end
a=factorial(n:-1:1);
if b>a;
    error('PERMS1:brange','Input Parameter b Should be between 1 and n!')
end;

p=zeros(1,n);                       %initialize output variable
idx=ones(1,n);                      %initialize loop counters
id=find(a<b,1);                     %locate the level at which the given index stands
for i=id-1:n-1                      %work down to figure out the counter states (idx) in each inner loop
    d=mod(b,a(i+1));                %how far down this loop?
    idx(i)=floor(b/a(i+1))+(d~=0);  %set counter, account for counter resetting
    b=d;
    if d==0;                        %adjust indices of inner loops if current loop is reset
        for j=i+1:n-1;
            idx(j)=n-j+1;
        end;
        break;
    end;
end;
m=1:n;                  %all possible indices
for i=1:n;
    p(i)=m(idx(i));     %set permutation index according to current counter index
    m=m(m~=p(i));       %current premutation index is no longer available
end;


text:
function exact = p47_exact ( )

%*****************************************************************************80
%
%% P47_EXACT returns the exact integral for problem 47.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the value of the integral.
%
  exact = - 4.0 / 9.0;

  return
end


text:
classdef matRad_ConstantRBEProjection < matRad_BackProjection
% matRad_BackProjection for optimization based on RBExDose
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2019 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    methods
        function obj = matRad_ConstantRBEProjection()
        end
    end

    methods
        function RBExD = computeSingleScenario(~,dij,scen,w)
            if ~isempty(dij.physicalDose{scen})
                RBExD = dij.physicalDose{scen} * (dij.RBE * w);
            else
                RBExD = [];
                matRad_cfg = MatRad_Config.instance();
                matRad_cfg.dispWarning('Empty scenario in optimization detected! This should not happen...\n');
            end
        end

        function wGrad = projectSingleScenarioGradient(~,dij,doseGrad,scen,~)
            if ~isempty(dij.physicalDose{scen})
                wGrad = ((dij.RBE * doseGrad{scen})' * dij.physicalDose{scen})';
            else
                wGrad = [];
                matRad_cfg = MatRad_Config.instance();
                matRad_cfg.dispWarning('Empty scenario in optimization detected! This should not happen...\n');
            end
        end
    end

end

repo: matRad


text:
%test script fgg_2D_experiment.m for the 2D NFFT based on Fast Gaussian
%Gridding.
%
%NOTE: In order for this FGG_2D to work, the C
%file "FGG_Convolution2D.c" must be compiled into a Matlab executable
%(cmex) with the following command: mex FGG_Convolution2D.c
%
%Code by (send correspondence to):
%Matthew Ferrara, Research Mathematician
%AFRL Sensors Directorate Innovative Algorithms Branch (AFRL/RYAT)
%Matthew.Ferrara@wpafb.af.mil

clear all;
close all;
%clc
%path(path, './nfftde');%NFFT mex files from Potts, et al.

% make an "image"
%Note for even lengths, the Nfft routine defines the image-space origin
% at pixel location [Nx/2 + 1, Ny/2 + 1].
% Convention:  x counts down rows, y counts L to R columns.
N=16;%even length assumed below...tedious to generalize...
z = zeros(N,N);%
%Make a smiley face:
z(N/2+3,N/2-1 : N/2+1 ) = 1;
z(N/2+2,N/2-2 ) = 1;
z(N/2+2,N/2+2 ) = 1;
z(N/2-1,N/2-1) = 1;
z(N/2+1,N/2) = 1;
z(N/2-1,N/2+1) = 1;
%imagesc(z)

N=[N,N];
img=double(z);
% Now, let's compute a matlab DFT in d dimensions using the "nfft" command.
% Note, use fftshifts to match indexing convention as used above in Pott's
% nfft.
data=fftshift(ifftn(ifftshift(img)));
DFTout = fftshift(fftn(ifftshift(data),N));
z=z(:);

% We need knots on [-1/2, 1-1/Nx]x[-1/2, 1-1/Ny] as fundamental period.
% make square grid of knots for exact comparison to fft
tmpx = linspace(-1/2,1/2 -1/N(1), N(1));% tmpx(end)=[];
tmpy = linspace(-1/2,1/2 -1/N(2), N(2));% tmpy(end)=[];

%this creates N+1 points, then discards point at +0.5.
%store my K knots as a d-by-K matrix (d=2 dimension here)
%...following four lines could be cleverly vectorized, avoiding loop.
[Y,X]=meshgrid(tmpy,tmpx);

knots=[X(:),Y(:)];

Nx=N(1);
Ny=N(2);
%set the desired number of digits of accuracy desired from the NUFFT
%routine:
Desired_accuracy = 6;%6=single precision, 12=double precision.
tic
MattOut_Gauss=FGG_2d_type1(data(:),knots,[Nx,Ny],Desired_accuracy);
imagesc(abs(MattOut_Gauss))
colorbar
title('(Type-I Fast Gaussian Gridding) NFFT output')
disp(['NUFFT evaluated in ',num2str(toc),' seconds'])
MattOut_t2=iFGG_2d_type2(MattOut_Gauss,knots,Desired_accuracy);
MattOut_Gauss=FGG_2d_type1(MattOut_t2(:),knots,[Nx,Ny],Desired_accuracy);
figure
imagesc(abs(MattOut_Gauss))
colorbar
title('(Type-I Fast Gaussian Gridding) NFFT output from Type-II-generated data')
norm(MattOut_t2-z(:))

figure
imagesc(abs(img-MattOut_Gauss))
colorbar
title('Error between DFT and NFFT')
Mean_Error=mean(abs(MattOut_Gauss(:)-DFTout(:)))


text:
function [i,pC,pE,Np] = spm_find_pC(varargin)
% Utility routine that finds the indices of non-zero covariance
% FORMAT [i,pC,pE,Np] = spm_find_pC(pC,pE,fields)
% FORMAT [i,pC,pE,Np] = spm_find_pC(DCM,fields)
% FORMAT [i,pC,pE,Np] = spm_find_pC(DCM)
%
% pC     - covariance matrix or variance stucture
% pE     - parameter structure
% fields - desired fields of pE
%
% or
%
% DCM    - DCM structure
%
% i      - find(diag(pC) > TOL)
% rC     - reduced covariances
% rE     - reduced expectation
%
%__________________________________________________________________________
% Copyright (C) 2015-2019 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_find_pC.m 7714 2019-11-26 11:25:50Z spm $

%-parse input arguments
%--------------------------------------------------------------------------
if nargin > 2
    pC     = varargin{1};
    pE     = varargin{2};
    fields = varargin{3};
elseif numel(varargin) > 1
    DCM    = varargin{1};
    fields = varargin{2};
else
    DCM    = varargin{1};
end

%-get prior density from DCM
%--------------------------------------------------------------------------
if nargin < 3
    if ischar(DCM)
        DCM = load(DCM,'DCM');
        DCM = DCM.DCM;
    end
    if any(isfield(DCM,{'options','M'}))
        try, [pC,pE] = spm_find_rC(DCM); end
    end
end

%-Deal with variance structures
%--------------------------------------------------------------------------
if isstruct(pC)
    q = spm_vec(pC);
else
    q = diag(pC);
end

%-Get indices
%--------------------------------------------------------------------------
i  = find(q > mean(q(q < 1024))/1024);
Np = numel(q);

%-subsample fields if necessary
%--------------------------------------------------------------------------
if nargin > 1
    if ischar(fields), fields = {fields}; end
    if isstruct(pE)
        j = spm_fieldindices(pE,fields{:});
        if isempty(j) && ~(~isempty(fields) && strcmp(fields{1},'none'))
            warning('%s not found. Returning all fields',...
                strjoin(cellstr(fields),','));
        else
            i = j(ismember(j,i));
        end
    end
end

return

function [pC,pE] = spm_find_rC(DCM)
% FORMAT [pC,pE] = spm_find_rC(DCM)
% model priors
%__________________________________________________________________________

% Get full priors and posteriors
%--------------------------------------------------------------------------
try
    pC  = DCM.M.pC;
    pE  = DCM.M.pE;
    return
end

% get priors from model specification
%--------------------------------------------------------------------------
if isfield(DCM.options,'spatial')

    % EEG or MEG
    %----------------------------------------------------------------------
    if strcmpi(DCM.options.analysis,'IND')
        [pE,dummy,pC] = spm_ind_priors(DCM.A,DCM.B,DCM.C,DCM.Nf);
    else
        [pE,  pC] = spm_dcm_neural_priors(DCM.A,DCM.B,DCM.C,DCM.options.model);

        try                                                     %#ok<TRYNC>
            try model   = DCM.options.model;   catch, model    = 'NMM'; end
            try spatial = DCM.options.spatial; catch, spatial  = 'LFP'; end
            DCM.M.dipfit.model = model;
            DCM.M.dipfit.type  = spatial;
            [pE,  pC] = spm_L_priors(DCM.M.dipfit,pE,pC);
            [pE,  pC] = spm_ssr_priors(pE,pC);
        end
    end

else

    % fMRI
    %----------------------------------------------------------------------
    if ~isfield(DCM,'d')
        DCM.d = zeros(size(DCM.a,1),size(DCM.a,1),0);
    end
    [pE,pC]   = spm_dcm_fmri_priors(DCM.a,DCM.b,DCM.c,DCM.d,DCM.options);
end

author: spm


text:
function view_bvt(lab1,valueMap)
    % plots the maxz LTA values calculated
    % with maxzlta.m or other similar values as a color map
    % needs valueMap, gx, gy, stri
    %
    % define size of the plot etc.
    %

    report_this_filefun();
    ZG=ZmapGlobal.Data;
    ZG.someColor = 'k';

    % Find out if figure already exists
    %
    bmapc=findobj('Type','Figure','-and','Name','b-value cross-section');

    % This is the info window text
    %
    ttlStr='The Z-Value Map Window                        ';
    hlpStr1zmap= ...
        ['                                                '
        ' This window displays seismicity rate changes   '
        ' as z-values using a color code. Negative       '
        ' z-values indicate an increase in the seismicity'
        ' rate, positive values a decrease.              '
        ' Some of the menu-bar options are               '
        ' described below:                               '
        '                                                '
        ' Threshold: You can set the maximum size that   '
        '   a volume is allowed to have in order to be   '
        '   displayed in the map. Therefore, areas with  '
        '   a low seismicity rate are not displayed.     '
        '   edit the size (in km) and click the mouse    '
        '   outside the edit window.                     '
        'FixAx: You can chose the minimum and maximum    '
        '        values of the color-legend used.        '
        'Polygon: You can select earthquakes in a        '
        ' polygon either by entering the coordinates or  '
        ' defining the corners with the mouse            '];
    hlpStr2zmap= ...
        ['                                                '
        'Circle: Select earthquakes in a circular volume:'
        '      Ni, the number of selected earthquakes can'
        '      be edited in the upper right corner of the'
        '      window.                                   '
        ' Refresh Window: Redraws the figure, erases     '
        '       selected events.                         '

        ' zoom: Selecting Axis -> zoom on allows you to  '
        '       zoom into a region. Click and drag with  '
        '       the left mouse button. type <help zoom>  '
        '       for details.                             '
        ' Aspect: select one of the aspect ratio options '
        ' Text: You can select text items by clicking.The'
        '       selected text can be rotated, moved, you '
        '       can change the font size etc.            '
        '       Double click on text allows editing it.  '
        '                                                '
        '                                                '];

    % Set up the Seismicity Map window Enviroment
    %
    if isempty(bmapc)
        bmapc = figure_w_normalized_uicontrolunits( ...
            'Name','b-value cross-section',...
            'NumberTitle','off', ...
            'backingstore','on',...
            'Visible','off', ...
            'Position',position_in_current_monitor(ZG.map_len(1), ZG.map_len(2)));
        % make menu bar

        lab1 = 'b-value';

        uicontrol('Units','normal',...
            'Position',[.0 .93 .08 .06],'String','Print ',...
            'callback',@callbackfun_001)

        uicontrol('Units','normal',...
            'Position',[.0 .75 .08 .06],'String','Close ',...
            'callback',@callbackfun_002)

        uicontrol('Units','normal',...
            'Position',[.0 .85 .08 .06],'String','Info ',...
            'callback',@callbackfun_003)

        uicontrol('Units','normal',...
            'Position',[.92 .80 .08 .05],'String','set ni',...
            'callback',@callbackfun_024)

        set_nia = uicontrol('style','edit','value',ni,'string',num2str(ni));
        set(set_nia, 'Callback', @callbackfun_025);
        set(set_nia,'units','norm','pos',[.94 .85 .06 .05],'min',10,'max',10000);
        nilabel = uicontrol('style','text','units','norm','pos',[.90 .85 .04 .05]);
        set(nilabel,'string','ni:','background',[.7 .7 .7]);

        ZG.tresh_km = max(r(:)); re4 = valueMap;
        nilabel2 = uicontrol('style','text','units','norm','pos',[.60 .92 .25 .06]);
        set(nilabel2,'string','MinRad (in km):','background',color_fbg);
        set_ni2 = uicontrol('style','edit','value',ZG.tresh_km,'string',num2str(ZG.tresh_km),...
            'background','y');
        set(set_ni2, 'Callback', @callbackfun_026)
        set(set_ni2,'units','norm','pos',[.85 .92 .08 .06],'min',0.01,'max',10000);

        uicontrol('Units','normal',...
            'Position',[.95 .93 .05 .05],'String','Go ',...
            'callback',@callbackfun_027)

        create_my_menu();

        colormap(jet)
    end   % This is the end of the figure setup

    % Now lets plot the color-map of the z-value
    %
    figure(bmapc);
    delete(findobj(bmapc,'Type','axes'));
    % delete(sizmap);
    reset(gca);
    cla
    set(gca,'NextPlot','replace')
    watchon;
    set(gca,'visible','off','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold',...
        'FontWeight','bold','LineWidth',1.5,...
        'Box','on','SortMethod','childorder')

    % set values greater ZG.tresh_km = nan
    %
    re4 = valueMap;
    re4(r > ZG.tresh_km) = nan;

    % plot image
    %
    orient portrait
    set(gcf,'PaperPosition', [2. 1 7.0 5.0])

    set(gca,'NextPlot','add')
    pco1 = pcolor(gx,gy,db12);
    axis([ min(gx) max(gx) min(gy) max(gy)])
    set(gca,'NextPlot','add');  shading flat; axis image

    hocm = hot;

    hocm(33,:) = [0.4 0.4 0.4];
    colormap(hocm)
    %brighten(0.7)
    caxis([-0.3 0.3])
    set(gca,'Color',[1.0 1.0 1.0]  )
    set(gca,'visible','on','FontSize',10,...
        'FontWeight','normal','LineWidth',1.5,...
        'Box','on','TickDir','out')
    set(gca,'NextPlot','add')

    % plot overlay
    %
    ploeqc = plot(newa(:,end),-newa(:,7),'.k');
    set(ploeqc,'Tag','eqc_plot','MarkerSize',ZG.ms6,'Marker',ty,'Color',ZG.someColor,'Visible','on')

    if exist('vox', 'var')
        plovo = plot(vox,voy,'*b');
        set(plovo,'MarkerSize',6,'LineWidth',1)
    end

    if exist('maix', 'var')
        pl = plot(maix,maiy,'*k');
        set(pl,'MarkerSize',12,'LineWidth',2)
    end

    if exist('maex', 'var')
        pl = plot(maex,-maey,'*k');
        set(pl,'MarkerSize',8,'LineWidth',2)
    end
    %overlay

    set(gca,'visible','on','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold',...
        'FontWeight','bold','LineWidth',1.5,...
        'Box','on','TickDir','out')
    h1 = gca;
    hzma = gca;

    % Create a colorbar
    %
    h5 = colorbar('horz');
    set(h5,'Pos',[0.35 0.05 0.4 0.04],...
        'FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.s,'TickDir','out')

    rect = [0.00,  0.0, 1 1];
    axes('position',rect)
    axis('off')
    %  Text Object Creation
    txt1 = text(...
        'Position',[ 0.33 0.07 0 ],...
        'HorizontalAlignment','right',...
        'FontSize',ZmapGlobal.Data.fontsz.s,....
        'FontWeight','bold',...
        'String',lab1);

    % Make the figure visible
    %
    axes(h1)
    set(gca,'visible','on','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold',...
        'FontWeight','bold','LineWidth',1.5,...
        'Box','on','TickDir','out')
    figure(bmapc);
    watchoff(bmapc)

    %% ui functions
    function create_my_menu()
        add_menu_divider();

    end

    %% callback functions

    function callbackfun_001(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        printdlg;
    end

    function callbackfun_002(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));

        f1=gcf;
        f2=gpf;
        set(f1,'Visible','off');
        close(bmapc);
        if f1~=f2
            figure(map);
        end
    end

    function callbackfun_003(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        helpdlg(help('view_bvt'),'Help for view_bvt')
    end

    function callbackfun_004(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        view_bv2(lab1,valueMap);
    end

    function callbackfun_005(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        h1 = gca;
        ZG=ZmapGlobal.Data;
        ZG.hold_state=false;
        cicros(1);
    end

    function callbackfun_006(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        h1 = gca;
        ZG=ZmapGlobal.Data;
        ZG.hold_state=false;
        cicros(2);
    end

    function callbackfun_007(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        h1 = gca;
        ZG=ZmapGlobal.Data;
        ZG.hold_state=true;
        cicros(0);
    end

    function callbackfun_008(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        h1 = gca;
        ZG=ZmapGlobal.Data;
        ZG.hold_state=false;
        cicros(3);
    end

    function callbackfun_011(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lab1 ='b-value change';
        valueMap = dbperc;
        view_bvt(lab1,valueMap);
    end

    function callbackfun_012(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lab1='b-value';
        valueMap = old;
        view_bvt(lab1,valueMap);
    end

    function callbackfun_013(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lab1='b-value';
        valueMap = meg;
        view_bvt(lab1,valueMap);
    end

    function callbackfun_014(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lab1='a-value';
        valueMap = avm;
        view_bvt(lab1,valueMap);
    end

    function callbackfun_015(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lab1='error in b';
        valueMap = stanm;
        view_bvt(lab1,valueMap);
    end

    function callbackfun_016(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lab1='difference in b';
        valueMap = old-meg;
        view_bvt(lab1,valueMap);
    end

    function callbackfun_017(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lab1='Mmax';
        valueMap = maxm;
        view_bvt(lab1,valueMap);
    end

    function callbackfun_018(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lab1='dM ';
        valueMap = maxm-magco;
        view_bvt(lab1,valueMap);
    end

    function callbackfun_019(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        def = {'6'};
        m = inputdlg('Magnitude of projected mainshock?','Input',1,def);
        m1 = m{:};
        m = str2num(m1);
        lab1 = 'Tr in yrs. (only smallest values shown)';
        valueMap =(teb - t0b)./(10.^(avm-m*old));
        view_bvt(lab1,valueMap);
    end

    function callbackfun_020(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lab1='Probability';
        valueMap = pro;
        view_bvt(lab1,valueMap);
    end

    function callbackfun_021(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lab1='Mcomp';
        valueMap = old1;
        view_bv2(lab1,valueMap);
    end

    function callbackfun_022(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lab1='Radius in [km]';
        valueMap = r;
        view_bvt(lab1,valueMap);
    end

    function callbackfun_024(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        ni=str2num(set_nia.String);
        num2str(ni);
    end

    function callbackfun_025(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
    end

    function callbackfun_026(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        ZG.tresh_km=str2double(set_ni2.String);
        set_ni2.String=num2str(ZG.tresh_km);
    end

    function callbackfun_027(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));

        pause(1);
        re4 =valueMap;
        view_bv2(lab1,valueMap);
    end
end

author: CelsoReyes


text:
function pyra = project_pyramid(model, pyra)
% pyra = project_pyramid(model, pyra)
%
% Project feature pyramid pyra onto PCA eigenvectors stored
% in model.coeff.

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2009-2012 Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

for i = 1:pyra.num_levels
  pyra.feat{i} = project(pyra.feat{i}, model.pca_coeff);
end

author: rbgirshick

repo: voc-dpm

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


text:
% SIN    Sine of argument in radians.
%    SIN(X) is the sine of the elements of X.
%
%    See also ASIN, SIND.
%
%    Reference page in Doc Center
%       doc sin
%
%    Other functions named sin
%
%       codistributed/sin    gpuArray/sin    sym/sin    ts/sin
%

author: jmaih


text:
classdef PTKGui < MimGuiBase
    % PTKGui. The user interface for the TD Pulmonary Toolkit.
    %
    %     To start the user interface, run ptk.m.
    %
    %     You do not need to modify this file. To add new functionality, create
    %     new plguins in the Plugins and GuiPlugins folders.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    methods
        function obj = PTKGui(splash_screen)
            obj = obj@MimGuiBase(PTKAppDef, splash_screen);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


function out=airdata(uu)
%
% Fake air data - this will be replaced with real air data in Chapter 4
%
% modified 12/11/2009 - RB

    % process inputs to function
    pn          = uu(1);             % North position (meters)
    pe          = uu(2);             % East position (meters)
    h           = -uu(3);            % altitude (meters)
    u           = uu(4);             % body velocity along x-axis (meters/s)
    v           = uu(5);             % body velocity along y-axis (meters/s)
    w           = uu(6);             % body velocity along z-axis (meters/s)
    phi         = 180/pi*uu(7);      % roll angle (degrees)
    theta       = 180/pi*uu(8);      % pitch angle (degrees)
    psi         = 180/pi*uu(9);      % yaw angle (degrees)
    p           = 180/pi*uu(10);     % body angular rate along x-axis (degrees/s)
    q           = 180/pi*uu(11);     % body angular rate along y-axis (degrees/s)
    r           = 180/pi*uu(12);     % body angular rate along z-axis (degrees/s)

    Va = sqrt(u^2+v^2+w^2);
    alpha = atan2(w,u);
    beta  = asin(v);
    wn    = 0;
    we    = 0;
    wd    = 0;

    out = [Va; alpha; beta; wn; we; wd];


text:
function test_bug1359

% MEM 3gb
% WALLTIME 00:10:00
% DEPENDENCY ft_read_header ft_read_data loadcnt

% the data is shared with test script 1490
filename = dccnpath('/home/common/matlab/fieldtrip/data/test/bug1490/sub1E3a.cnt');

hdr = ft_read_header(filename);
dat = ft_read_data(filename);

% The report is that the data looks weird and does not contain any
% negative values. I have not checked with old reading functions, but
% with FieldTrip revision 5913 (just after fixing the related bug 1490)
% the data looks meaningful. The channels have large offsets, some appear
% to have clipped, but there are channels that show a physiological signal
% with electrode drift. So it seems that the problem has disappeared.

assert(any(dat(:)>0) && any(dat(:)<0));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% Gabor Filter demo
%
% A GUI to interact with the 5 different Gabor filter parameters, while
% visualizing the resulting filter.
%

function varargout = gabor_filter_gui(ksize)
    % create the UI
    if nargin < 1, ksize = [121 121]; end
    h = buildGUI(ksize);
    if nargout > 0, varargout{1} = h; end
end

function onChange(~,~,h)
    %ONCHANGE  Event handler for UI controls

    % retrieve current values from UI controls
    sigma  = get(h.slid(5), 'Value') / 10;
    theta  = get(h.slid(4), 'Value') * pi/180;
    lambda = get(h.slid(3), 'Value');
    gamma  = get(h.slid(2), 'Value') / 100;
    psi    = get(h.slid(1), 'Value') * pi/180;

    % create Gabor filter
    kernel = cv.getGaborKernel('KSize',[h.ksize(2) h.ksize(1)], ...
        'Sigma',sigma, 'Theta',theta, 'Lambda',lambda, ...
        'Gamma',gamma, 'Psi',psi);

    % normalize filter to [0,1] range and resize it
    kernel = cv.normalize(kernel, 'NormType','MinMax');
    kernel = cv.resize(kernel, [h.sz(2) h.sz(1)]);

    % show result
    set(h.img, 'CData',kernel)
    set(h.txt(5), 'String',sprintf('Sigma  = %.2f',sigma))
    set(h.txt(4), 'String',sprintf('Theta  = %.2f',theta))
    set(h.txt(3), 'String',sprintf('Lambda = %.2f',lambda))
    set(h.txt(2), 'String',sprintf('Gamma  = %.2f',gamma))
    set(h.txt(1), 'String',sprintf('Psi    = %.2f',psi))
    drawnow
end

function onType(~,e,h)
    %ONTYPE  Event handler for key press on figure

    % handle keys
    switch e.Key
        case {'q', 'escape'}
            close(h.fig)
            return
        case 'h'
            onHelp([],[]);
        case 'r'
            onReset([],[],h);
    end
end

function onReset(~,~,h)
    set(h.slid(5), 'Value',400);    % sigma
    set(h.slid(4), 'Value',0);      % theta
    set(h.slid(3), 'Value',11);     % lambda
    set(h.slid(2), 'Value',100);    % gamma
    set(h.slid(1), 'Value',90);     % psi
    onChange([],[],h);
end

function onHelp(~,~)
    %ONHELP  Display usage help dialog

    helpdlg({
        'This GUI allows to interact with the 5 different Gabor filter'
        'parameters, while visualizing the resulting filter.'
        ''
        'Hot keys:'
        'ESC, q - quit the program'
        'r - reset parameters to original values'
        'h - this help dialog'
    });
end

function h = buildGUI(ksize)
    %BUILDGUI  Creates the UI

    % parameters
    sigma = 400; sigma_max = 1000;
    theta = 0;   theta_max = 180;
    lambda = 11; lambda_max = 100;
    gamma = 100; gamma_max = 200;
    psi = 90;    psi_max = 180;
    sz = [512 512];

    % build the user interface (no resizing to keep it simple)
    h = struct();
    h.ksize = ksize;  % size of the filter
    h.sz = sz;        % size of the image to show
    h.fig = figure('Name','Gabor Filter Demo', ...
        'NumberTitle','off', 'Menubar','none', 'Resize','off', ...
        'Position',[200 200 sz(2) sz(1)+129]);
    if ~mexopencv.isOctave()
        %HACK: not implemented in Octave
        movegui(h.fig, 'center');
    end
    h.ax = axes('Parent',h.fig, 'Units','pixels', 'Position',[1 130 sz(2) sz(1)]);
    if ~mexopencv.isOctave()
        h.img = imshow(zeros(sz), 'Parent',h.ax);
    else
        %HACK: https://savannah.gnu.org/bugs/index.php?45473
        axes(h.ax);
        h.img = imshow(zeros(sz));
    end
    text(5, 5, sprintf('KSize = %dx%d', ksize(2), ksize(1)), ...
        'Color','y', 'VerticalAlignment','top');

    props = {'Parent',h.fig, 'Style','text', 'String','', ...
        'FontSize',11, 'HorizontalAlignment','left'};
    h.txt(5) = uicontrol(props{:}, 'Position',[5   5 120 20]);
    h.txt(4) = uicontrol(props{:}, 'Position',[5  30 120 20]);
    h.txt(3) = uicontrol(props{:}, 'Position',[5  55 120 20]);
    h.txt(2) = uicontrol(props{:}, 'Position',[5  80 120 20]);
    h.txt(1) = uicontrol(props{:}, 'Position',[5 105 120 20]);

    props = {'Parent',h.fig, 'Style','slider', 'Min',0};
    h.slid(5) = uicontrol(props{:}, 'Position',[125   5 sz(2)-125-5 20], ...
        'Value',sigma, 'Max',sigma_max, 'SliderStep',[10 100]./(sigma_max-0));
    h.slid(4) = uicontrol(props{:}, 'Position',[125  30 sz(2)-125-5 20], ...
        'Value',theta, 'Max',theta_max, 'SliderStep',[2 20]./(theta_max-0));
    h.slid(3) = uicontrol(props{:}, 'Position',[125  55 sz(2)-125-5 20], ...
        'Value',lambda, 'Max',lambda_max, 'SliderStep',[1 10]./(lambda_max-0));
    h.slid(2) = uicontrol(props{:}, 'Position',[125  80 sz(2)-125-5 20], ...
        'Value',gamma, 'Max',gamma_max, 'SliderStep',[2 20]./(gamma_max-0));
    h.slid(1) = uicontrol(props{:}, 'Position',[125 105 sz(2)-125-5 20], ...
        'Value',psi, 'Max',psi_max, 'SliderStep',[2 20]./(psi_max-0));

    % hook event handlers, and trigger default start
    opts = {'Interruptible','off', 'BusyAction','cancel'};
    set(h.slid, 'Callback',{@onChange,h}, opts{:});
    set(h.fig, 'WindowKeyPressFcn',{@onType,h}, opts{:});
    onChange([],[],h);
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% out2 = readtimeseries('/afs/inf.ed.ac.uk/user/s08/s0880377/project/data/qdots/R1/100222/7W5UG85Y_F00000014.tif','',[],1,1);
% out14=squeeze(out2(:,:,:,0));
% o=dip_image(out14);
o=dip_image(data_out);
box1=4;
box2=2;
nclick = 10;
dipshow(7,mean(o,[],3));
coords = dipgetcoords(7,nclick);
box1x1 = coords(:,1) - box1;
box1x2 = coords(:,1) + box1;
box1y1 = coords(:,2) - box1;
box1y2 = coords(:,2) + box1;
for ii=1: 7%length(coords)
    oroi1 = o(box1x1(ii):box1x2(ii), box1y1(ii):box1y2(ii),:);
    [maxval, maxcoord] = max(max(oroi1,[],3));
    box2x1 = maxcoord(1)-box2;
    box2x2 = maxcoord(1)+box2;
    box2y1 = maxcoord(2)-box2;
    box2y2 = maxcoord(2)+box2;
    oroi2= oroi1(box2x1:box2x2,box2y1:box2y2,:);
    f(ii,:) = squeeze(double(sum(oroi2,[],[1 2])));

end

index=1;
figure(1)
for ii=1:N
    subplot(N,2,index)
    plot(f(ii,:))
      if ii<N
        set(gca,'xtick',[],'ytick',[])
    else
        set(gca,'ytick',[])
    end
    index = index +1;
    subplot(N,2,index)
    hist(f(ii,:),50);
    xlim([1.5,4]*10^5)
    if ii<N
        set(gca,'xtick',[],'ytick',[])
    else
        set(gca,'ytick',[])
    end
    index = index +1;
end
SaveImageFULL(['block' num2str(q) '_timeprofilesingledot'])
figure(2)
bar(max(f,[],2),'r')
hold on
bar(mean(f,2),'g')
bar(min(f,[],2),'b')
errorbar([1:8], mean(f,2), sqrt(var(f,[],2)),'+k')
l{1}='max';
l{2}='mean';
l{3}='min';
l{4}='sqrt(var)';
legend(l);
grid on
hold off
SaveImageFULL(['block' num2str(q) '_statistics'])

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function rmReconCrossCorrelate(numOfStimuli,roi,method,datatype)

    if ~exist('HOMEDIR','var') || isempty(HOMEDIR)
        homedir = pwd;
    else
        homedir = HOMEDIR;
    end

    if ~exist('numOfStimuli','var') || isempty(numOfStimuli)
        error('Need the number of different stimuli.');
    end

    if ~exist('method','var') || isempty(method)
        method = 'betasum';
        fprintf('[%s] No method supplied, using default (betasum) \n',mfilename);
    end

    if ~exist('roi','var') || isempty(roi)
        roi ='V1';
        fprintf('[%s] No ROI supplied, using default (V1) \n',mfilename);
    end

    numstim = numOfStimuli;

    correl = zeros(numstim,numstim);
    pval = zeros(numstim,numstim);
    cond = single(zeros(128*128,numstim));

    for i = 1:numstim
        stimuli(i) = load([homedir '/Recon/RFcovs/' datatype '/' roi '/' method '/RFcov-stim' num2str(i) '.mat']);

        covstim = load([homedir '/Gray/' datatype '/methods/' method '/CovPlot_' roi '_stim' num2str(i) '-' method '.mat']);

        cond(:,i) = single(reshape(covstim.CovPlot,1,128*128));
    end

    for si = 1 : numstim
        for si2 = 1 : numstim

            [cor, p] = corr(cond(:,si), stimuli(si2).RFcov(:));

            correl(si,si2) = cor;
            pval(si,si2) = p;
        end
    end

    save([homedir '/Recon/Correlations/CrossCorr-' roi '-' method '.mat'],'correl');
    save([homedir '/Recon/Correlations/CrossPval-' roi '-' method '.mat'],'pval');
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function format_colorbar(axes_handle, scalar)
chkarg(istypesizeof(scalar, 'Scalar2d') || istypesizeof(scalar, 'Scalar3d'), ...
	'"scalar" shuold be instancef of either Scalar2d or Scalar3d.');

% axes_handle = colorbar;

% Make sure that 0 is ticked when 0 is the minimum of the colorbar.
crange = caxis(gca);  % caxis is the property of the axes handle for a figure
ticks = get(axes_handle, 'YTick');
if crange(1) == 0 && ticks(1) > 0
	crange(1) = -eps * crange(2);
	caxis(gca, crange);
	ticks = [0 ticks];
end
set(axes_handle, 'YTick', ticks);

% Prevent 'x10^n' from showing on top of the colorbar.
ticklabel = textscan(num2str(ticks),'%s');
ticklabel = ticklabel{1};
set(axes_handle, 'YTickLabel', ticklabel);

% if scalar.unitvalue ~= 1e0
% 	str = [' (', char(hex2dec('00D7')), num2str(scalar.unitvalue, '%.2e'), ')'];
% 	title(axes_handle, str, 'HorizontalAlignment', 'left', 'VerticalAlignment', 'baseline');
% end
str = [' (', char(hex2dec('00D7'))];
if scalar.unitvalue < 1e4 && scalar.unitvalue > 1e-3
	str = [str, num2str(scalar.unitvalue), ')'];
else
	str = [str, num2str(scalar.unitvalue, '%.2e'), ')'];
end
title(axes_handle, str, 'HorizontalAlignment', 'left', 'VerticalAlignment', 'baseline');

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function varargout = line(rot,varargin)
% draw rotations connected by lines
%
% Syntax
%   line(rot,'linecolor','r','linewith',2)
%
% Input
%  rot - list of @rotation
%
% Example
%   cs1 = crystalSymmetry('321')
%   cs2 = crystalSymmetry('432')
%   oR = fundamentalRegion(cs1,cs2)
%   plot(oR)
%
%   % connect to vertices of the fundamental region
%   f = fibre(oR.V(1),oR.V(2))
%
%   % determine some orientations along the fibre
%   o = f.orientation
%
%   % plot the line
%   hold on
%   line(o,'color','r','linewidth',2)
%   hold off
%
% See also
%

[varargout{1:nargout}] = plot(reshape(rot,[],1),'all','edgecolor',...
  get_option(varargin,{'color','linecolor'},'k'),'Marker','none',varargin{:});

author: mtex-toolbox

repo: mtex


text:
% *************************************************************************
% Video Super-Resolution with Convolutional Neural Networks
%
% helper function to generate mat file from video frames
% --> feel free to modify and use for your own purpose
%
% Version 1.0
%
% Created by:   Armin Kappeler
% Date:         02/19/2016
%
% *************************************************************************
clearvars

OUTPUTPATH = 'data/lr_video/u3/foreman_u4_LR.mat';

DParam.prescaling = 1;
DParam.upscaleFactor = 4;
DParam.patchSize = 36;
DParam.stride = 36;

inputBasename = {'foreman_short'};
inputNrFrames = 20;

inputVidBasePath = '/home/armin/caffe/data/Superresolution/4vid/';
inputVidPath{1} = [inputVidBasePath inputBasename{1} '/original/Frame '];
inputVidPath{2} = '.png';

for i=1:inputNrFrames
    fileIdxStr = sprintf('%03d',i);
    filename = [inputVidPath{1} fileIdxStr inputVidPath{2}];

    imgRGB = im2double(imread(filename));

    %% only get Y-channel (gray) from YCbCr
    if (ndims(imgRGB)==3)
        img_ycbcr = rgb2ycbcr(imgRGB);
        imgHi = img_ycbcr(:,:,1);
    else
        imgHi = imgRGB;
    end

    %% preprocess image
    imgHi = crop2divisibleSize(imgHi,DParam.upscaleFactor);
    imgLo = imresize(imgHi,1/DParam.upscaleFactor,'bicubic');

    if i==1
        frames = zeros([size(imgLo),inputNrFrames] );
    end

    frames(:,:,i) = imgLo;

end

save(OUTPUTPATH,'frames')

author: thomas-koehler

repo: SupER


text:
function plotData(X, y)
%PLOTDATA Plots the data points X and y into a new figure
%   PLOTDATA(x,y) plots the data points with + for the positive examples
%   and o for the negative examples. X is assumed to be a Mx2 matrix.

% Create New Figure
figure; hold on;

% ====================== YOUR CODE HERE ======================
% Instructions: Plot the positive and negative examples on a
%               2D plot, using the option 'k+' for the positive
%               examples and 'ko' for the negative examples.
%

% Find Indices of Positive and Negative Examples
pos = find(y==1); neg = find(y==0);

% Plot Examples
plot(X(pos, 1), X(pos, 2), 'k+', 'LineWidth', 2, 'MarkerSize', 7);
plot(X(neg, 1), X(neg, 2), 'ko', 'MarkerFaceColor', 'y', 'MarkerSize', 7);

% =========================================================================

hold off;

end

author: rmarquis

repo: coursera-machinelearning


text:
function VO = spm_smoothto8bit(V,fwhm)
% 3 dimensional convolution of an image to 8bit data in memory
% FORMAT VO = spm_smoothto8bit(V,fwhm)
% V     - mapped image to be smoothed
% fwhm  - FWHM of Guassian filter width in mm
% VO    - smoothed volume in a form that can be used by the
%         spm_*_vol.mex* functions.
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_smoothto8bit.m 4310 2011-04-18 16:07:35Z guillaume $

if nargin>1 && fwhm>0,
    VO = smoothto8bit(V,fwhm);
else
    VO = V;
end
return;
%_______________________________________________________________________

%_______________________________________________________________________
function VO = smoothto8bit(V,fwhm)
% 3 dimensional convolution of an image to 8bit data in memory
% FORMAT VO = smoothto8bit(V,fwhm)
% V     - mapped image to be smoothed
% fwhm  - FWHM of Guassian filter width in mm
% VO    - smoothed volume in a form that can be used by the
%         spm_*_vol.mex* functions.
%_______________________________________________________________________

vx   = sqrt(sum(V.mat(1:3,1:3).^2));
s    = (fwhm./vx./sqrt(8*log(2)) + eps).^2;
r    = cell(1,3);
for i=1:3,
    r{i}.s = ceil(3.5*sqrt(s(i)));
    x      = -r{i}.s:r{i}.s;
    r{i}.k = exp(-0.5 * (x.*x)/s(i))/sqrt(2*pi*s(i));
    r{i}.k = r{i}.k/sum(r{i}.k);
end

buff = zeros([V.dim(1:2) r{3}.s*2+1]);

VO        = V;
VO.dt     = [spm_type('uint8') spm_platform('bigend')];
V0.dat    = uint8(0);
V0.dat(VO.dim(1:3)) = uint8(0);
VO.pinfo  = [];

for i=1:V.dim(3)+r{3}.s,
    if i<=V.dim(3),
        img      = spm_slice_vol(V,spm_matrix([0 0 i]),V.dim(1:2),0);
        msk      = find(~isfinite(img));
        img(msk) = 0;
        buff(:,:,rem(i-1,r{3}.s*2+1)+1) = ...
            conv2(conv2(img,r{1}.k,'same'),r{2}.k','same');
    else
        buff(:,:,rem(i-1,r{3}.s*2+1)+1) = 0;
    end

    if i>r{3}.s,
        kern    = zeros(size(r{3}.k'));
        kern(rem((i:(i+r{3}.s*2))',r{3}.s*2+1)+1) = r{3}.k';
        img     = reshape(buff,[prod(V.dim(1:2)) r{3}.s*2+1])*kern;
        img     = reshape(img,V.dim(1:2));
        ii      = i-r{3}.s;
        mx      = max(img(:));
        mn      = min(img(:));
        if mx==mn, mx=mn+eps; end
        VO.pinfo(1:2,ii) = [(mx-mn)/255 mn]';
        VO.dat(:,:,ii)   = uint8(round((img-mn)*(255/(mx-mn))));
    end
end


text:
function disp_(p)
%DISP_        Display of polynom in "disp_" mode
%

% written  07/22/02     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  olddisp = intvalinit('display',0);
  intvalinit('display_',0);
  display(p,inputname(1));
  intvalinit(olddisp,0);

repo: OpenMAS


text:
function rows = sb_sparse_to_mat(diinsy);

% SB_SPARSE_TO_MAT
%
% $Id$

rows = zeros(max(diinsy),1);
rows(diinsy) = 1;
rows = [1;rows];
rows = cumsum(rows);
rows = rows(1:end-1);
end


text:
function [ output_meta ] = meta2sicd_sio( sio_hdr, user_data )
%META2SICD_SIO Converts SIO header into SICD-style structure
%
% Really not much to do here since SIO header is so minimal.  Just fill in
% image size and datatype.
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

if nargin<2
    user_data=struct();
end

% Check to see if we have a SICD-in-SIO type file.  If there's a SICDMETA
% field in the user_data structure we'll assume that has all the valid
% metadata for this file and we'll use that.  Otherwise we'll just return a
% small 'stub' SICD metadata structure (not much information in the SIO
% header).
if (isfield(user_data,'SICDMETA'))
  output_meta = sicdxml2struct( xmlread( java.io.StringBufferInputStream(user_data.SICDMETA) ) );
else
  output_meta.ImageData.NumRows=sio_hdr(2);
  output_meta.ImageData.NumCols=sio_hdr(3);
  output_meta.ImageData.FullImage=output_meta.ImageData; % Assume full image, but no way to know for sure
  output_meta.ImageData.FirstRow=uint32(0); output_meta.ImageData.FirstCol=uint32(0);
  if (sio_hdr(4)==13)&&(sio_hdr(5)==8), output_meta.ImageData.PixelType='RE32F_IM32F'; end;
  output_meta.ImageData.SCPPixel.Row=fix(output_meta.ImageData.NumRows/2); % Not usually given explicitly in accompanying metadata, so just guess center
  output_meta.ImageData.SCPPixel.Col=fix(output_meta.ImageData.NumCols/2);
end
end

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

author: ngageoint


text:
function [mm]=misssum(X,def)
%MISSSUM sum of a matrix X with NaN's
%
%[mm]=misssum(X,def)
%
%This function calculates the sum of a matrix X.
%X may hold missing elements denoted by NaN's which
%are ignored.
%
%The result is standardized, that is, corrected for the lower
%number of contributing terms.
%
%Check that for no column of X, all values are missing

% Copyright (C) 1995-2006  Rasmus Bro & Claus Andersson
% Copenhagen University, DK-1958 Frederiksberg, Denmark, rb@life.ku.dk
%
% 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
% Street, Fifth Floor, Boston, MA  02110-1301, USA.

%Insert zeros for missing, correct afterwards
missidx = isnan(X);
i = find(missidx);
if ~isempty(i),
    X(i) = zeros(size(i));
end;

%Find the number of real(non-missing objects)
if min(size(X))==1,
   n_real=length(X)-sum(missidx);
   weight=length(X);
else
   n_real=size(X,1)-sum(missidx);
   weight=size(X,1);
end

i=find(n_real==0);
if isempty(i) %All values are real and can be corrected
   mm=weight*sum(X)./n_real;
else %There are columns with all missing, insert missing
   n_real(i)=1;
   mm=weight*sum(X)./n_real;
   mm(i)=i + NaN;
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function [all_theta] = oneVsAll(X, y, num_labels, lambda)
%ONEVSALL trains multiple logistic regression classifiers and returns all
%the classifiers in a matrix all_theta, where the i-th row of all_theta
%corresponds to the classifier for label i
%   [all_theta] = ONEVSALL(X, y, num_labels, lambda) trains num_labels
%   logisitc regression classifiers and returns each of these classifiers
%   in a matrix all_theta, where the i-th row of all_theta corresponds
%   to the classifier for label i

% Some useful variables
m = size(X, 1);
n = size(X, 2);

% You need to return the following variables correctly
all_theta = zeros(num_labels, n + 1);

% Add ones to the X data matrix
X = [ones(m, 1) X];

% ====================== YOUR CODE HERE ======================
% Instructions: You should complete the following code to train num_labels
%               logistic regression classifiers with regularization
%               parameter lambda.
%
% Hint: theta(:) will return a column vector.
%
% Hint: You can use y == c to obtain a vector of 1's and 0's that tell use
%       whether the ground truth is true/false for this class.
%
% Note: For this assignment, we recommend using fmincg to optimize the cost
%       function. It is okay to use a for-loop (for c = 1:num_labels) to
%       loop over the different classes.
%
%       fmincg works similarly to fminunc, but is more efficient when we
%       are dealing with large number of parameters.
%
% Example Code for fmincg:
%
%     % Set Initial theta
%     initial_theta = zeros(n + 1, 1);
%
%     % Set options for fminunc
%     options = optimset('GradObj', 'on', 'MaxIter', 50);
%
%     % Run fmincg to obtain the optimal theta
%     % This function will return theta and the cost
%     [theta] = ...
%         fmincg (@(t)(lrCostFunction(t, X, (y == c), lambda)), ...
%                 initial_theta, options);
%

initial_theta = zeros(n + 1, 1);

options = optimset('GradObj', 'on', 'MaxIter', 50);

for c = 1:num_labels
    all_theta(c,:) = fmincg (@(t)(lrCostFunction(t, X, (y == c), lambda)), initial_theta, options);

% =========================================================================

end

author: Borye

path: github-repos/MATLAB/Borye-machine-learning-coursera-1/machine-learning-coursera-1-033fdc2e6da393eeb1179a09aafe92362021effb/Week 4 Assignments/Multi-class Classification and Neural Networks/mlclass-ex3/oneVsAll.m


text:
function a = cauchy ( n, x, y )

%*****************************************************************************80
%
%% CAUCHY returns the CAUCHY matrix.
%
%  Formula:
%
%    A(I,J) = 1.0 / ( X(I) + Y(J) )
%
%  Example:
%
%    N = 5, X = ( 1, 3, 5, 8, 7 ), Y = ( 2, 4, 6, 10, 9 )
%
%    1/3  1/5  1/7  1/11 1/10
%    1/5  1/7  1/9  1/13 1/12
%    1/7  1/9  1/11 1/15 1/14
%    1/10 1/12 1/14 1/18 1/17
%    1/9  1/11 1/13 1/17 1/16
%
%    or, in decimal form,
%
%    0.333333      0.200000      0.142857      0.0909091     0.100000
%    0.200000      0.142857      0.111111      0.0769231     0.0833333
%    0.142857      0.111111      0.0909091     0.0666667     0.0714286
%    0.100000      0.0833333     0.0714286     0.0555556     0.0588235
%    0.111111      0.0909091     0.0769231     0.0588235     0.0625000
%
%  Properties:
%
%    A is generally not symmetric: A' /= A.
%
%    A is totally positive if 0 < X(1) < ... < X(N) and 0 < Y1 < ... < Y(N).
%
%    A will be singular if any X(I) equals X(J), or
%    any Y(I) equals Y(J), or if any X(I)+Y(J) equals zero.
%
%    A is generally not normal: A' * A /= A * A'.
%
%    The Hilbert matrix is a special case of the Cauchy matrix.
%
%    The Parter matrix is a special case of the Cauchy matrix.
%
%    The Ris or "ding-dong" matrix is a special case of the Cauchy matrix.
%
%    det ( A ) = product ( 1 <= I < J <= N ) ( X(J) - X(I) )* ( Y(J) - Y(I) )
%           / product ( 1 <= I <= N, 1 <= J <= N ) ( X(I) + Y(J) )
%
%    The inverse of A is
%
%      INVERSE(A)(I,J) = product ( 1 <= K <= N ) [ (X(J)+Y(K)) * (X(K)+Y(I)) ] /
%            [ (X(J)+Y(I)) * product ( 1 <= K <= N, K /= J ) (X(J)-X(K))
%                          * product ( 1 <= K <= N, K /= I ) (Y(I)-Y(K)) ]
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2015
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Gregory, David Karney,
%    Example 3.26,
%    A Collection of Matrices for Testing Computational Algorithms,
%    Wiley, New York, 1969, page 54,
%    LC: QA263.G68.
%
%    Nicholas Higham,
%    Accuracy and Stability of Numerical Algorithms,
%    SIAM, 1996.
%
%    Donald Knuth,
%    The Art of Computer Programming,
%    Volume 1, Fundamental Algorithms, Second Edition
%    Addison-Wesley, Reading, Massachusetts, 1973, page 36.
%
%    Olga Taussky, Marvin Marcus,
%    Eigenvalues of finite matrices,
%    in Survey of Numerical Analysis,
%    Edited by John Todd,
%    McGraw-Hill, New York, pages 279-313, 1962.
%
%    Evgeny Tyrtyshnikov,
%    Cauchy-Toeplitz matrices and some applications,
%    Linear Algebra and Applications,
%    Volume 149, 1991, pages 1-18.
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real X(N), Y(N), vectors that determine A.
%
%    Output, real A(N,N), the matrix.
%
  for i = 1 : n

    for j = 1 : n

      if ( x(i) + y(j) == 0.0 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'CAUCHY - Fatal error!\n' );
        fprintf ( 1, '  The denominator X(I)+Y(J) was zero\n' );
        fprintf ( 1, '  for I = %d\n', i );
        fprintf ( 1, '  X(I)  = %g\n', x(i) );
        fprintf ( 1, '  and J = %d\n', j );
        fprintf ( 1, '  Y(J)  = %g\n', y(j) );
        error ( 'CAUCHY - Fatal error!' );
      end

      a(i,j) = 1.0 / ( x(i) + y(j) );

    end
  end

  return
end

author: johannesgerer


text:
function [At,b,c,K] = opti2sedumi(prob)
%OPTI2SEDUMI  Converts an OPTI Problem (f,A,b,sdcone) to SeDuMi Format
% [At,b,c,K] = opti2sedumi(prob)

%   Copyright (C) 2013 Jonathan Currie (IPL)

%Check not already in sedumi format
if(isfield(prob,'sdcone') && ~isempty(prob.sdcone) && isstruct(prob.sdcone))
    At = getSDMIField(prob.sdcone,'At','A');
    b  = getSDMIField(prob.sdcone,'b','B');
    c  = getSDMIField(prob.sdcone,'c','C');
    K  = getSDMIField(prob.sdcone,'K','k');
    return;
end

%Check we have a linear f
if(~isfield(prob,'f') || isempty(prob.f) || ~isnumeric(prob.f))
    error('The SDP to convert must have a linear objective f in prob.f');
end
b = -prob.f; ndec = length(b);

%Convert linear constraints to general format
if(~isempty(prob.rl) || ~isempty(prob.ru))
    [AA,bb,Aeq,beq] = row2gen(prob.A,prob.rl,prob.ru);
else
    AA = prob.A; bb = prob.b; Aeq = prob.Aeq; beq = prob.beq;
end
%Convert all constraints and bounds to linear inequalities
[At,c] = inequal(AA,bb,Aeq,beq,prob.lb,prob.ub);
K.l = size(At,1);

%Check for Semidefinite constraints, concat as required
if(~isempty(prob.sdcone))
    if(iscell(prob.sdcone))
        for i = 1:length(prob.sdcone)
            chkSDDim(prob.sdcone{i},ndec,i);
            K.s(i) = sqrt(size(prob.sdcone{i},1));
            c = [c;-prob.sdcone{i}(:,1)]; %#ok<AGROW>
            At = [At;-prob.sdcone{i}(:,2:end)]; %#ok<AGROW>
        end
    else
        chkSDDim(prob.sdcone,ndec,1);
        K.s = sqrt(size(prob.sdcone,1));
        c = [c;-prob.sdcone(:,1)];
        At = [At;-prob.sdcone(:,2:end)];
    end
end

%Check for quadratic constraints
if(~isempty(prob.Q)), optiwarn('opti:sedumiQC','Currently quadratic constraints are not converted to SeDuMi format'); end

function chkSDDim(cone,ndec,i)
if(size(cone,2) ~= ndec+1)
    error('Semidefinite Cone %d does not have the correct number of columns',i);
end
dim = sqrt(size(cone,1));
if(floor(dim) ~= dim)
    error('Semidefinite Cone %d does not have the correct number of rows to form a square matrix');
end

function val = getSDMIField(conestr,f1,f2)
if(isfield(conestr,f1))
    val = conestr.(f1);
elseif(isfield(prob.sdcone,f2))
    val = conestr.(f2);
else
    error('It appears the problem is already in SeDuMi format, but OPTI could not find prob.sdcone.%s',f1);
end

author: Arrowstar

repo: ksptot


text:
function [isValid, errMsg] = validate_patient_module(dcmobj)
%"validate_patient_module"
%   Takes a Java dcmobj that contains the fields representing a patient,
%   and checks the validity of the fields and their contents.  This
%   includes a general check of fields based on type (Class 1,2,3,1C,2C)
%   and individual field check for any fields with special conditions on
%   them.
%
%JRA 06/12/06
%
%Usage:
%   [isValid, errMsg] = validate_patient_module(dcmobj)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

%Set cache buffer to true.
buf = 1;

%Set character set to [] for now.
cs = [];

%Default to valid, no error message.
isValid = 1; errMsg = '';

%Get the parent of the passed dcmobj, in case the validation process
%requires checking tags in other modules.
parent = dcmobj.getParent;
if ~isempty(parent)
    dcmobj = parent;
else
    %Do nothing.  This is the root object.
end

%Get the tags for the patient module.
[Class_1, Class_2, Class_3, Class_1C, Class_2C] = patient_module_tags;

%Send tags and dcmobj for general checking of field format and presence.
%See documentation of validate_general_fields.
[isValid, errMsg] = validate_general_fields(dcmobj, Class_1, Class_2, Class_3, Class_1C, Class_2C);
if ~isValid
    return;
end

%Create an iterator over the dcmobj to inspect each element
it = dcmobj.iterator;

%-----------------------Individual Field Checking-------------------------%
%Switch over top level tags in patient module, with individual handling for
%each tag that requires verification beyond the general field validation
%performed above.
while it.hasNext

    el  = it.next;
    tag = el.tag;
    switch tag
        case hex2dec('00100040')
            %"Patient's Sex" must contain M, F, O or null.
            validValues = {'M', 'F', 'O', ''};
            val = char(el.getString(cs, buf));
            if ~ismember(val, validValues);
               isValid      = 0;
               [name, hex]  = tag2strings(tag, dcmobj);
               errMsg       = [hex ' "' name '" must contain M, F, O or null.'];
               return;
            end

        case hex2dec('00081120')
            %Handle Referenced Patient Sequence here.
            [isValid, errMsg] = validate_referenced_patient_sequence(el, dcmobj);
            if ~isValid
                return;
            end

        case hex2dec('00120062')
            %"Patient Identity Removed" must contain YES, NO or null.
            validValues = {'YES', 'NO', ''};
            val = char(el.getString(cs, buf));
            if ~ismember(val, validValues);
               isValid      = 0;
               [name, hex]  = tag2strings(tag, dcmobj);
               errMsg       = [hex ' "' name '" must contain YES, NO or null.'];
               return;
            end

            %In addition, if YES, at least one of the two "de-identification
            %method" elements must be present.
            if strcmpi(val, 'YES') & ~dcmobj.contains(hex2dec('00120063')) & ~dcmobj.contains(hex2dec('00120064'))
               isValid      = 0;
               [name, hex]  = tag2strings(tag, dcmobj);
               errMsg       = [hex ' "' name '" contains YES, but no de-identification element is present.'];
               return;
            end

        case hex2dec('00120064')
            %Handle de-identification Method Code Sequence here.
            [isValid, errMsg] = validate_deidentification_method_code_sequence(el, dcmobj);
            if ~isValid
                return;
            end

        otherwise
            %Tag does not require special handling or is not in the patient
            %module.
    end

end

function [isValid, errMsg] = validate_deidentification_method_code_sequence(el, dcmobj)
%"validate_deidentification_method_code_sequence"
%   Validates the passed sequence element.
n   = el.countItems;
tag = el.tag;
isValid = 1; errMsg = '';

%Deidentification method code sequence must contain one or more child elements.
if n < 1
    isValid     = 0;
    [name, hex] = tag2strings(tag, dcmobj);
    errMsg      = [hex ' "' name '" must contain one or more elements.'];
end

for i=0:n-1
   childobj = el.getDicomObject(i);
   %Send each child object for independent validation by the macro's
   %validator.
   [isValid, errMsg] = validate_code_sequence_macro(dcmobj);
   if ~isValid
       return;
   end
end

function [isValid, errMsg] = validate_referenced_patient_sequence(el, dcmobj)
%"validate_referenced_patient_sequence"
%   Validates the passed sequence element.
n   = el.countItems;
tag = el.tag;
isValid = 1; errMsg = '';

%Referenced patient sequence must contain exactly one child element.
if n ~= 1
    isValid     = 0;
    [name, hex] = tag2strings(tag, dcmobj);
    errMsg      = [hex ' "' name '" must contain exactly one element.'];
end

%Pass the child elements for general validation, as Class 1 since
%conditions for their existance have been met.
Class_1 = hex2dec({'00081150','00081155'});
for i=0:n-1
    childobj = el.getDicomObject(i);
    [isValid, errMsg] = validate_general_fields(childobj, Class_1, [], [], [], []);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [fid, sens, label, sens_label] = read_zebris(Fname_zeb,skip)

% Reads Zebris files:
%   fiducials locations, and
%   either sensor file or headshape file or both
%
% FORMAT [fid, sens, label] = read_zebris(Fname_zeb,skip)
% Input:
% Fname_zeb  - Zebris ASCII file containing sensor locations (mm)
%             (headshape can also be considered here instead of sensors)
% skip       - first channels to skip
%
% Output:
% fid        - fiducial         locations (mm) in rows
% sens       - sensor/headshape locations (mm) in rows
% label      - labels of the fiducials
% sens_label - labels of the surface points, electrodes + headshape
%
% IMPORTANT: Note that Zebris data files should be -ASCII files with
% extension .sfp
% It is assumed that the .sfp file contains the location (mm) of fiducials
% (possibly twice), possibly followed by some additional named points for
% the electrodes, and then so more named location starting with 'sfl' for
% headshape locations.
% In some instances the first few channel locations may pertain to
% reference channels; the skip variable allows these to be skipped if
% necessary.
% The fiducial locations are flaged with the strings 'fidt9','fidnz' and
% 'fidt10'; indicating the leaft ear, nasion, and right ear, respectively.
% _________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Christophe Phillips
% $Id$

% checks and assigments
%--------------------------------------------------------------------------
try, skip;    catch, skip    = 0; end

[pth,nam,ext] = fileparts(Fname_zeb);
if ~strcmp(ext,'.sfp')
  warndlg(sprintf('Wrong input file format\n'));
  return
end

% --- READ Zebris Sensor + fiducial locations ---
%==========================================================================
try
  file = textread(Fname_zeb,'%s');
catch
  file = textread(fullfile(pwd,[nam ext]),'%s');
end
% remove zeros at the end
bool = 0;
while bool == 0
  if strcmp(file{end},'0')
    file = file(1:end-1);
  else
    bool = 1;
  end
end

% read fiducials
%--------------------------------------------------------------------------
NZ   = [];
LE   = [];
RE   = [];
temp = 0;
nl   = 1;

while temp == 0
  if strcmp(file{nl},'fidnz')
    NZ = [NZ ; str2num(file{nl+1}) str2num(file{nl+2}) str2num(file{nl+3}) ];
    nl = nl + 4;
  elseif strcmp(file{nl},'fidt9')
    LE = [LE ; str2num(file{nl+1}) str2num(file{nl+2}) str2num(file{nl+3}) ];
    nl = nl + 4;
  elseif strcmp(file{nl},'fidt10')
    RE = [RE ; str2num(file{nl+1}) str2num(file{nl+2}) str2num(file{nl+3}) ];
    nl = nl + 4;
  else
    temp = 1;
  end
end

% convert from cm to mm
%--------------------------------------------------------------------------
NZ    = mean(NZ,1); LE = mean(LE,1); RE = mean(RE,1);
fid   = [NZ; LE; RE];

label = [{'nas', 'lpa', 'rpa'}];

% read sensor locations or headshape locations
%--------------------------------------------------------------------------
sens  = [];
sens_label = [];

start = nl + skip*3;
for ii = start:4:length(file)
  sens_label = [sens_label, file(ii)];
  sens = [sens; ...
    str2num(file{ii+1}) str2num(file{ii+2}) str2num(file{ii+3})];
end


text:
function l = l1pp ( n, h )

%*****************************************************************************80
%
%% L1PP stores the 1D PP Laplacian as a full matrix.
%
%  Discussion:
%
%    The N grid points are assumed to be evenly spaced by H.
%
%    For N = 5, the discrete Laplacian with periodic boundary conditions
%    has the matrix form L:
%
%       2 -1  0  0 -1
%      -1  2 -1  0  0
%       0 -1  2 -1  0
%       0  0 -1  2 -1
%      -1  0  0 -1  2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 October 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points.
%    N must be at least 3.
%
%    Input, real H, the spacing between points.
%
%    Output, real L(N,N), the Laplacian matrix.
%
  if ( n < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'L1PP - Fatal error!\n' );
    fprintf ( 1, '  N < 3.\n' );
    error ( 'L1PP - Fatal error!' );
  end

  l = zeros ( n, n );

  i = 1;
  l(1,1) =  2.0 / h / h;
  l(1,2) = -1.0 / h / h;
  l(1,n) = -1.0 / h / h;

  for i = 2 : n - 1
    l(i,i-1) = -1.0 / h / h;
    l(i,i) =    2.0 / h / h;
    l(i,i+1) = -1.0 / h / h;
  end

  i = n;
  l(n,1) =   -1.0 / h / h;
  l(n,n-1) = -1.0 / h / h;
  l(n,n) =    2.0 / h / h;

  return
end

repo: jburkardt-m