text:
%% Figure 1 of RSDS paper
% ``Rotation, Scaling and Deformation Invariant Scattering
% for Texture Discrimination"
% Laurent Sifre, Stephane Mallat.
% Proc. IEEE CVPR 2013 Portland, Oregon
%
%%

clear; close all;

%% build a filter bank of oriented and dilated morlet wavelets
size_in = [512, 512];
options.min_margin = [0, 0];
filters = morlet_filter_bank_2d(size_in, options);

%% build two grids of dirac
dirac_step = 64;
dirac_grid_1 = zeros(size_in);
dirac_grid_1(1:dirac_step:end, 1:dirac_step:end) = 1;
dirac_grid_2 = zeros(size_in);
dirac_grid_2(dirac_step/2+1:dirac_step:end, ...
    dirac_step/2+1:dirac_step:end) = 1;

%% convolves the grid with morlet filters
dirac_grid_1_f = fft2(dirac_grid_1);
dirac_grid_2_f = fft2(dirac_grid_2);

psi_1 = filters.psi.filter{17};
psi_2 = filters.psi.filter{17+4};

texture_1 = conv_sub_2d(dirac_grid_1_f, psi_1, 0);
texture_1 = texture_1 + conv_sub_2d(dirac_grid_1_f, psi_2, 0);

texture_2 = conv_sub_2d(dirac_grid_1_f, psi_1, 0);
texture_2 = texture_2 + conv_sub_2d(dirac_grid_2_f, psi_2, 0);

%% take real part
texture_1 = real(texture_1);
texture_2 = real(texture_2);

%% display the two textures
subplot(121);
imagesc(texture_1);
subplot(122);
imagesc(texture_2);
colormap gray;

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function a = gridsel(alg,hyper)

%==============================================================
% GRIDSEL model selection object via grid search
%==============================================================
% A=GRIDSEL(A,H) returns a gridsel object initialized with
% algorithms A and hyperparameters H.
%
% Finds the best of the algorithm from the set A, and trains
% and stores that model.
%
% Hyperparameters.
%
%  a.child=A                -- methods to evaluate
%  a.loss='class_loss'      -- loss measure to use
%  a.score=cv('folds=5')    -- method of evaluating algorithms
%
% Model
%  a.scores=[]      % score of all of methods tried
%  a.best_score=[]  % score of best methods tried
%  a.best_index=[]  % index of best methods tried
%  a.best=[]        % learnt model of best method tried
%
% Example:
% % train 3 svms with C=1,2,3 and validate with 3 fold cross validation
% [r,a]=train(gridsel(param(svm,'C',[1,2,3]),{'score=cv;score.folds=3'}),gen(toy)) ;
%
% Methods:
%  train, test
%=============================================================
% Reference :
% Author    :
% Link      :
%=============================================================

  % hypers
  a.child=alg;     % original algorithm to start with
  a.score=cv; a.score.folds=5;
  a.loss='class_loss';

  % model
  a.all_algs=[];   % all trained models
  a.scores=[];     % score of all of methods tried
  a.best_score=[]; % score of best method tried
  a.best_index=[]; % index of best method tried
  a.best=[];       % learnt model of best method tried (without cv, etc.)

  p=algorithm('gridsel');
  a= class(a,'gridsel',p);

  if nargin==2,
    eval_hyper;
  end;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [data,mask,test,masktest]= genadding(problem)
global   mzeros convert
data=mzeros(problem.numsamples, 4 , ceil(problem.T*1.1));
mask=mzeros(problem.numsamples, 1 , ceil(problem.T*1.1));
test=mzeros(problem.numtest, 4 , ceil(problem.Ttest*1.1));
masktest=mzeros(problem.numtest, 1 , ceil(problem.Ttest*1.1));

T=fix(problem.T/1.1);
for i = 1 : problem.numsamples

    length = T+ fix(rand  * T / 10) ;
    mask(i,1,length)=1;

    data(i,3,1:length)=ones(1, length); %bias
    a=2*rand(1,length    ) - 1;
    data(i,1,1:length)=a;
    a=a';

    b=zeros(length,1);  %flag
    c=20*ones(length,1);

    data(i,4,1:length)=ones(1,length);%bias

    if T>10
        tmp1= ceil( rand * 10);
    else
        tmp1= ceil( rand * T/2);
    end

    data(i,2,tmp1)=1;
    b(tmp1,1)=1;

    tmp2=tmp1;
    while tmp2 == tmp1
        tmp2=  ceil(rand * T/2);
    end

    data(i,2,tmp2) = 1;
    b(tmp2,1)=1;

    if tmp2 ==1  || tmp1 ==1
        data(i,2,1)=0;
        b(1,1)=0;
    else
        data(i,2,1)=-1;
        b(1,1)=-1;
    end

    data(i,2,length-1)=-1;
    b(end,1)=-1;
    data(i,1,length) = 0.5+(data(i,1,tmp1)+data(i,1,tmp2))/4;
    c(end,1)= 0.5+(a(tmp1,1)+a(tmp2,1))/4;

end

for i =   size(data,3) :-1:1
    if max( max( data(:,:,i))) ~=0
        data=data(:,:,1:i);
        break
    end
end

for i = 1 : problem.numtest

    length = T+ fix(rand  * T / 10) ;
    masktest(i,1,length)=1;

    test(i,3,1:length)=ones(1, length); %bias
    test(i,1,1:length)=2*rand(1,length    ) - 1;

    test(i,4,1:length) = ones(1,length);%bias

    if T>10
        tmp1= ceil( rand * 10);
    else
        tmp1= ceil( rand * T/2);
    end

    test(i,2,tmp1)=1;

    tmp2=tmp1;
    while tmp2 == tmp1
        tmp2=  ceil(rand * T/2);
    end

    test(i,2,tmp2) = 1;

    if tmp2 ==1  || tmp1 ==1
        test(i,2,1)=0;
    else
        test(i,2,1)=-1;
    end

    test(i,2,length-1)=-1;
    test(i,1,length) = 0.5+(test(i,1,tmp1)+test(i,1,tmp2))/4;
end
for i =   size(test,3) :-1:1
    if max( max( test(:,:,i))) ~=0
        test=test(:,:,1:i);
        break
    end
end

tmp=test;
test=mzeros(size(test,1) ,  size(test,2) ,  size(test,3) +1  );
test(:,:,1)=tmp(:,:,1);
test(:,:,2:end)=tmp;

masktest=masktest(:,:,1:size(test,3));
tmp=masktest;
masktest=mzeros(size(masktest,1) ,  size(masktest,2) ,  size(masktest,3) +1  );
masktest(:,:,1)=tmp(:,:,1);
masktest(:,:,2:end)=tmp;

data=convert(data);
mask=convert(mask);
test=convert(test);
masktest=convert(masktest);

author: huashiyiqike

repo: LSTM-MATLAB


text:
function [GaitDBDir,probes,numprobe]=testData

GaitDBDir=['USFGait17_32x22x10/'];

%If you want to test on higher resolution, please download them from
%https://sites.google.com/site/tensormsl/#data
%or http://www.dsp.toronto.edu/~haiping/MSL#data
%and change the directory in this script.
%If you are using 32-bit Matlab, you may have out of memory problem. Please
%try:
%1. clear unused variables
%2. write matrix operations into for-loops
%3. use integer data type instead of single/double float point
%4. when passing variables to functions: save variables to mat files, clear
%the varialbes, load the variables in the calling function so Matlab won't
%keep two copies of the same data.

%GaitDBDir=['USFGait17_64x44x20/'];
%GaitDBDir=['USFGait17_128x88x20/'];

probes=['A','B','C','D','E','F','G'];
numprobe=length(probes);

author: Sable

repo: mcbench-benchmarks


text:
function [Xopt,yopt,Sopt,obj] = recover_mosek_sol(res,n)
idxtril         = tril(true(n,n));
Xopt            = zeros(n,n);
Xopt(idxtril)   = res.sol.itr.barx;
Xopt            = Xopt + Xopt';
Xopt            = Xopt - 0.5*diag(diag(Xopt));
yopt            = res.sol.itr.y;
Sopt            = zeros(n,n);
Sopt(idxtril)   = res.sol.itr.bars;
Sopt            = Sopt + Sopt';
Sopt            = Sopt - 0.5*diag(diag(Sopt));
obj             = [res.sol.itr.pobjval;res.sol.itr.dobjval];
end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
function y = log_normcdf( x, approx )

%LOG_NORMCDF   Logarithm of the cumulative normal distribution.
%   Y = LOG_NORMCDF(X) is the logarithm of the CDF of the normal
%   distribution at the point X.
%
%                                1    / x
%       LOG_NORMCDF(X) = LOG( ------- |   exp(-t^2/2) dt )
%                             sqrt(2) / -Inf
%
%   For numeric X, LOG_NORMCDF(X) is computed using the equivalent
%   expression LOG(0.5*ERFC(-X*SQRT(0.5))). When X is a CVX variable, a
%   a piecewise quadratic *approximation* is employed instead. This
%   approximation gives good results when -4 <= x <= 4, and will be
%   improved in future releases of CVX.
%
%   For array values of X, the LOG_NORMCDF returns an array of identical
%   size with the calculation applied independently to each element.
%
%   X must be real.
%
%   Disciplined convex programming information:
%       LOG_NORMCDF is concave and nondecreasing in X. Therefore, when used
%       in CVX specifications, X must be concave.

error(nargchk(1,2,nargin));
if ~isreal( x ),
    error( 'Argument must be real.' );
end
if nargin > 1,
    % For debugging purposes only
    y = cvx_constant(log_normcdf(cvx(x)));
else
    y = log(0.5*erfc(-x*sqrt(0.5)));
end

repo: REST


text:
function [varargout] = minDose(planC, varargin)
%Calculate a minDose metric.
%Stand alone metric :       minDose(planC, structNum, doseNum, 'Percent' or 'Absolute');
%Request m object   :       minDose(planC, 'getnewmetric');
%Evaluate m object  :       minDose(planC, m, 'evaluate');
%LM:  4 Jan 06, JOD: changed interface.
%
% 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/>.

indexS = planC{end};
optS = planC{indexS.CERROptions};

if(nargin > 1)
    call = varargin{end};
else
    warning('Incorrect Usage, try: minDose(planC, structNum, doseNum, ''Percent'' or ''Absolute'');');
    return;
end

switch upper(call)
    case 'GETNEWMETRIC'
        m.name = 'minDose';
        m.valueV = [];
        m.description = 'Returns the min dose in specified structure';
        m.functionName = @minDose;
        m.note = '';
        m.params(1).name = 'Structure';
        m.params(1).type = 'DropDown';
        m.params(1).list = {planC{indexS.structures}.structureName};
        m.params(1).value = 1;
        m.params(2).name = 'IsTarget?';
        m.params(2).type = 'DropDown';
        m.params(2).list ={'Target', 'Not Target'};
        m.params(2).value = 2;
        m.units = [];
        m.range = [0 inf];
        m.doseSets = [1];
        varargout = {planC,m};
    return;

    case 'EVALUATE'
        m = varargin{1};
        structName = planC{indexS.structures}(m.params(1).value).structureName;
        m.note = structName;
        m.valueV = [];

        maxD = -inf;
		for i=1:length(m.doseSets)
            [planC, doseBinsV, volsHistV] = getDVHMatrix(planC, m.params(1).value, m.doseSets(i));
            ans = calc_minDose(doseBinsV, volsHistV, m.params(2).value);
            m.valueV = [m.valueV ans];
            doseArray = getDoseArray(planC{indexS.dose}(m.doseSets(i)));
            maxD = max(maxD, max(doseArray(:)));
        end
        if m.params(2).value == 1
            %Is target, high dose is best.
            m.range = [0 max(m.valueV)];
        elseif m.params(2).value == 2;
            m.range = [maxD min(m.valueV)];
        end
        m.units = getDoseUnitsStr(i,planC);
        varargout = {planC,m};
        return;

    otherwise  %No specific call has been made, assume user just wants raw metric.
        if(nargin > 3)
            structNum = varargin{1};
            doseNum = varargin{2};
            volumeTypeString = varargin{3};
            if(strcmp(upper(volumeTypeString), 'ABSOLUTE'))
                volumeType = 2;
            else
                volumeType = 1;
            end
            [planC, doseBinsV, volsHistV] = getDVHMatrix(planC, structNum, doseNum);
            varargout = {calc_minDose(doseBinsV, volsHistV, volumeType)};
        else
            error('Not enough parameters, try: minDose(planC, structNum, doseNum, ''Percent'' or ''Absolute'');')
            return;
        end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
classdef TukeyLoss < dagnn.Loss

  properties
    sigma = 1.
  end

  methods
    function outputs = forward(obj, inputs, params)
      outputs{1} = vl_nntukeyloss(inputs{1}, inputs{2}, ...
                       'instanceWeights', inputs{3}, ...
                       obj.opts{:}) ;
      n = obj.numAveraged ;
      m = n + size(inputs{1}, 4) ;
      obj.average = (n * obj.average + gather(outputs{1})) / m ;
      obj.numAveraged = m ;
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      derInputs{1} = vl_nntukeyloss(inputs{1}, inputs{2}, derOutputs{1}, ...
                                'instanceWeights', inputs{3}, ...
                                obj.opts{:}) ;
      derInputs{2} = [] ;
      derInputs{3} = [] ;
      derParams = {} ;
    end

    function obj = TukeyLoss(varargin)
      obj.load(varargin) ;
    end
  end
end

repo: MFT


text:
%function particles= resampleParticles(particles, Nmin, fp_debug, resample_randstream)
function particles= resampleParticles(particles, Nmin) %, resample_randstream
%
% Resample particles if their weight variance is such that N-effective
% is less than Nmin.
%

N = length(particles);
w = [particles.w];
ws = sum(w);
w  = w/ws;

Neff = 1 / sum(w .^ 2);

if Neff < Nmin
    %fprintf(fp_debug, 'Particles resampled, Neff = %d\n', Neff);
    %disp(['Particles resampled, Neff = ', Neff]);

    %[keep]    = stratified_resample(w, resample_randstream);%, resample_randstream
    [keep]    = stratified_resample(w);
    %fprintf(fp_debug, 'Keep particles: ');
    particles = particles(keep);
%     for i=1:N
%         particles(i).w= 1/N;
%         fprintf(fp_debug, '%d ', keep(i));
%
%     end
    %fprintf(fp_debug, '\n');
    %disp('Keep particles: ');disp(keep);
else
    for i=1:N
        particles(i).w= particles(i).w / ws;
    end
end

repo: Wheel-SLAM


text:
% Source code directory
addpath('source');

% Add path for using interpolation code vgg_interp2
addpath('external\imrender\vgg');

% Add VLFeat
addpath('external\vlfeat-0.9.20\toolbox\mex\mexw64');

% Add Mean-shift
addpath('external\MeanShift');

% Clear the screen
clc

repo: StructCompletion


text:
function [headinfo,nav,obs,tobs,fid]=decode_obsh(headinfo,nav,obs,fid)

global glc gls
tobs(1:glc.MAXOBSTYPE,1:3,glc.NSYS)=' '; del=zeros(3,1);
sta=gls.sta;

while ~feof(fid)
    line=fgets(fid);label=line(61:end);
    if size(line,2)<=60
        continue;
    elseif ~isempty(strfind(label,'MARKER NAME'))
        str=line(1:60);idx=find(str~=' ');
        if any(idx)
            str(idx(end)+1:end)=[];
            sta.name=str;
        end
    elseif ~isempty(strfind(label,'MARKER NUMBER'))
        str=line(1:20);idx=find(str~=' ');
        if any(idx)
            str(idx(end)+1:end)=[];
            sta.maker=str;
        end
    elseif ~isempty(strfind(label,'MAKER TYPE'))
    elseif ~isempty(strfind(label,'OBSERVER / AGENCY'))
    elseif ~isempty(strfind(label,'REC # / TYPE / VERS'))
        str=line(1:20);idx=find(str~=' ');
        if any(idx)
            str(idx(end)+1:end)=[];
            sta.recsno=str;
        end
        str=line(21:40);idx=find(str~=' ');
        if any(idx)
            str(idx(end)+1:end)=[];
            sta.rectype=str;
        end
        str=line(41:60);idx=find(str~=' ');
        if any(idx)
            str(idx(end)+1:end)=[];
            sta.recver=str;
        end
    elseif ~isempty(strfind(label,'ANT # / TYPE'))
        str=line(1:20);idx=find(str~=' ');
        if any(idx)
            str(idx(end)+1:end)=[];
            sta.antno=str;
        end
        str=line(21:40);idx=find(str~=' ');
        if any(idx)
            str(idx(end)+1:end)=[];
            sta.antdes=str;
        end
    elseif ~isempty(strfind(label,'APPROX POSITION XYZ'))
        j=1;
        for i=1:3
            sta.pos(i)=str2double(line(j:j+14));
            j=j+14;
        end
    elseif ~isempty(strfind(label,'ANTENNA: DELTA H/E/N'))
        j=1;
        for i=1:3
            del(i)=str2double(line(j:j+14));
            j=j+14;
        end
        sta.del(1)=del(2);sta.del(2)=del(3);sta.del(3)=del(1);
    elseif ~isempty(strfind(label,'ANTENNA: DELTA X/Y/Z'))
    elseif ~isempty(strfind(label,'ANTENNA: PHASECENTER'))
    elseif ~isempty(strfind(label,'ANTENNA: B.SIGHT XYZ'))
    elseif ~isempty(strfind(label,'ANTENNA: ZERODIR AZI'))
    elseif ~isempty(strfind(label,'ANTENNA: ZERODIR XYZ'))
    elseif ~isempty(strfind(label,'CENTER OF MASS: XYZ'))
    elseif ~isempty(strfind(label,'SYS / # / OBS TYPES' ))%ver 3.0
        syscode=line(1);
        switch syscode
            case 'G',i=glc.SYS_GPS;
            case 'R',i=glc.SYS_GLO;
            case 'E',i=glc.SYS_GAL;
            case 'C',i=glc.SYS_BDS;
            case 'J',i=glc.SYS_QZS;
            otherwise,continue;
        end
        n=str2double(line(4:7));
        k=7;nt=1;
        for j=1:n
            if k>58
                line=fgetl(fid);
                k=7;
            end
            str=line(k+1:k+3);%idx=find(tmp~=' ');tmp(idx(end)+1:end)=[];
            tobs(nt,:,i)=str;
            nt=nt+1;
            k=k+4;
        end
        if i==4&&(headinfo.ver-3.02)<1e-3
            for j=1:nt
                if tobs(j,2,i)=='1'
                    tobs(j,2,i)='2';
                end
            end
        end
    elseif ~isempty(strfind(label,'WAVELENGTH FACT L1/2'))
    elseif ~isempty(strfind(label,'# / TYPES OF OBSERV'))%ver 2.0
        n=str2double(line(1:7));% number of obs type
        j=10;nt=1;
        for i=1:n
            if j>58
                line=fgetl(fid);
                j=10;
            end
            if headinfo.ver<=2.99
                str=line(j+1:j+2);idx=find(str~=' ');str(idx(end)+1:end)=[];
                tobs(nt,:,1)=convcode(headinfo.ver,str,glc.SYS_GPS);
                tobs(nt,:,2)=convcode(headinfo.ver,str,glc.SYS_GLO);
                tobs(nt,:,3)=convcode(headinfo.ver,str,glc.SYS_GAL);
                tobs(nt,:,4)=convcode(headinfo.ver,str,glc.SYS_BDS);
                tobs(nt,:,5)=convcode(headinfo.ver,str,glc.SYS_QZS);
            end
            nt=nt+1;
            j=j+6;
        end
    elseif ~isempty(strfind(label,'SIGNAL STRENGTH UNIT'))
    elseif ~isempty(strfind(label,'INTERVAL'))
    elseif ~isempty(strfind(label,'TIME OF FIRST OBS'))
        if     strcmp(line(49:51),'GPS'), headinfo.tsys=glc.TSYS_GPS;
        elseif strcmp(line(49:51),'GLO'), headinfo.tsys=glc.TSYS_GLO;
        elseif strcmp(line(49:51),'GAL'), headinfo.tsys=glc.TSYS_GAL;
        elseif strcmp(line(49:51),'BDT'), headinfo.tsys=glc.TSYS_BDS;
        elseif strcmp(line(49:51),'QZS'), headinfo.tsys=glc.TSYS_QZS;
        end
    elseif ~isempty(strfind(label,'TIME OF LAST OBS'))
    elseif ~isempty(strfind(label,'RCV CLOCK OFFS APPL'))
    elseif ~isempty(strfind(label,'SYS / DCBS APPLIED'))
    elseif ~isempty(strfind(label,'SYS / PCVS APPLIED'))
    elseif ~isempty(strfind(label,'SYS / SCALE FACTOR'))
    elseif ~isempty(strfind(label,'SYS / PHASE SHIFTS'))
    elseif ~isempty(strfind(label,'GLONASS SLOT / FRQ #'))
        p=5;
        for i=1:8
            prn=str2double(line(p+1:p+2));
            fcn=str2double(line(p+4:p+5));
            if prn>=1&&prn<=glc.MAXPRNGLO
                nav.glo_fcn(prn)=fcn+8;
                p=p+7;
            end
        end
    elseif ~isempty(strfind(label,'GLONASS COD/PHS/BIS'))
        p=1;
        for i=1:4
            if     strcmp(line(p+1:p+3),'C1C'),nav.glo_cpbias(1)=str2double(line(p+4:p+4+8));
            elseif strcmp(line(p+1:p+3),'C1P'),nav.glo_cpbias(2)=str2double(line(p+4:p+4+8));
            elseif strcmp(line(p+1:p+3),'C2C'),nav.glo_cpbias(3)=str2double(line(p+4:p+4+8));
            elseif strcmp(line(p+1:p+3),'C2P'),nav.glo_cpbias(4)=str2double(line(p+4:p+4+8));
            end
            p=p+13;
        end
    elseif ~isempty(strfind(label,'LEAP SECONDS'))
        nav.leaps=str2double(line(1:6));
    elseif ~isempty(strfind(label,'# OF SALTELLITES'))
    elseif ~isempty(strfind(label,'PRN / # OF OBS'))
    elseif ~isempty(strfind(label,'END OF HEADER'))
        break;
    end
end
obs.sta=sta;

return

repo: GINav


text:
function [stack,dim]=read_LSstack_fast_float(fpath,dim)
%% reads whole stack

info=dir(fpath);
if (isempty(info))
    error('File Does not exist');
end

stack=read_LSstack_fast_float_mex64(fpath,int32(dim(1:2)));

size(stack)

stack=reshape(stack,dim);

author: xiuyechen

repo: FishExplorer


text:
function ebsd = flipud(ebsd)
% flip spatial ebsd-data from upside down
%
% Input
%  ebsd - @EBSD
%
% Output
%  flipped ebsd - @EBSD

m(1) = max(ebsd.prop.x);
m(2) = max(ebsd.prop.y);
ebsd = affinetrans(ebsd,[],[0 -m(2)]);

ebsd = rotate(ebsd,rotation.byAxisAngle(xvector,pi));

author: mtex-toolbox

repo: mtex

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


text:
function c8_print_test ( )

%*****************************************************************************80
%
%% C8_PRINT_TEST tests C8_PRINT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 February 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'C8_PRINT_TEST\n' );
  fprintf ( 1, '  C8_PRINT prints a C8.\n' );
  fprintf ( 1, '\n' );

  c1 = 0.0 + 0.0 * i;
  c2 = 1.0 + 0.0 * i;
  c3 = 3.141592653589793 + 0.0 * i;
  c4 =  0.0 + 1.0 * i;
  c5 = 1.0 + 2.0 * i;
  c6 = -12.34 + 56.78 * i;
  c7 = 0.001 + 0.000002 * i;
  c8 = 3.0E+08 - 4.5E+09 * i;

  c8_print ( c1, '  Zero:' );
  c8_print ( c2, '  One:' );
  c8_print ( c3, '  Pi:' );
  c8_print ( c4, '  i:' );
  c8_print ( c5, '  1+2i:' );
  c8_print ( c6, ' -12.34 + 56.78i:' );
  c8_print ( c7, '  1E-3 + 2E-6i' );
  c8_print ( c8, '  3E8 - 4.5E9i:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [theta, J_history] = gradientDescent(X, y, theta, alpha, num_iters)
%GRADIENTDESCENT Performs gradient descent to learn theta
%   theta = GRADIENTDESCENT(X, y, theta, alpha, num_iters) updates theta by
%   taking num_iters gradient steps with learning rate alpha

% Initialize some useful values
m = length(y); % number of training examples
J_history = zeros(num_iters, 1);

for iter = 1:num_iters

    % ====================== YOUR CODE HERE ======================
    % Instructions: Perform a single gradient step on the parameter vector
    %               theta.
    %
    % Hint: While debugging, it can be useful to print out the values
    %       of the cost function (computeCost) and gradient here.
    %

    predictions =  X * theta;
    updates = X' * (predictions - y);
    theta = theta - alpha * (1/m) * updates;

    %theta = theta - alpha * (1/m) * sum(sqerrors) * X;

    %theta - (alpha/m) * (X' * (X * theta - y));

    %theta = theta - (alpha/m) * (X' * (X * theta - y));

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

    % Save the cost J in every iteration
    J_history(iter) = computeCost(X, y, theta);

end

end

author: vkosuri

repo: CourseraMachineLearning


text:
function scale = p40_scale ( neqn )

%*****************************************************************************80
%
%% P40_SCALE returns scale factors for problem p40.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real SCALE(NEQN), the scaling factors.
%
  scale(1:neqn,1) = 1.0;

  return
end


text:
function label = litekmeans(X, k)
% Perform k-means clustering.
%   X: d x n data matrix
%   k: number of seeds
% Written by Michael Chen (sth4nth@gmail.com).
n = size(X,2);
last = 0;
label = ceil(k*rand(1,n));  % random initialization
while any(label ~= last)
    [u,~,label] = unique(label);   % remove empty clusters
    k = length(u);
    E = sparse(1:n,label,1,n,k,n);  % transform label into indicator matrix
    m = X*(E*spdiags(1./sum(E,1)',0,k,k));    % compute m of each cluster
    last = label;
    [~,label] = max(bsxfun(@minus,m'*X,dot(m,m,1)'/2),[],1); % assign samples to the nearest centers
end
[~,~,label] = unique(label);


text:
function [M, movie] = tSeriesMovie(vw, scans, anatFlag, varargin)
% [M movie] = tSeriesMovie(vw, scans, anatFlag, [options]):
%
% Make a movie of the selected tSeries, using the current view's
% settings, and play using mplay viewer. Intended as a replacement
% for makeTSeriesMovie.
%
% NOTE: The code auto-resizes the movie to fit the screen nicely. Although
% this is a benefit, it usually causes the lower line of the mplay tool to
% turn black; this line contains some useful information like the current
% frame and image size. To get this info back, press the second button
% from the left ("truesize") in the movie window's toolbar.
%
% vw: the view struct. Currently supports inplane and flat level views.
%
% scans: set of scans to view from current data type.
%
% anatFlag: if 1, will threshold the tSeries and superimpose
% the parts above threshold on top of the view's anat image.
% Defaults to 0.
%
% Returns M, a matrix of the tSeries, and movie, an object-
% oriented struct generated by mplay used for programmatic
% control of the movie. If anatFlag is 0, M is size
% x by y by nFrames, where x and y are the view's slice dims.
% If anatFlag is 1, M is a 4-D array of size x by y by 3 by nFrames,
% where each 3-D subvolume is a truecolor image.
%
% Options:  [case-insensitive]
%   slices, [val]:          set slice numbers (defaults to current slice).
%                           WARNING: >1 slice can be very memory-hungry
%   zoom, [val]:            set zoom area (as with axis command)
%   detrend, [val]:         set detrend flag (see detrendTSeries for vals)
%   convertToPct,[flag]:    if 1, will convert to % signal change.
%   applyHisto,[flag]:      flag to apply standard histogram criterion (on
%                           by default) in which the tails of the intensity
%                           histogram are clipped.
%   compareFrames,[frame]:  if nonzero, will make a truecolor move in which
%                           the red channel is the tSeries, and the
%                           blue-green channel is a reference frame --
%                           e.g., if 1, the first frame of the tSeries
%   meanThresh,[val]:       set threshold for raw functionals -- any parts
%                           of the tSeries below this value (normalized
%                           from 0 to 1 --> min to max of tSeries) will
%                           be zeroed.
%   funcClip,[val]:         set clip values for functionals, if
%                           superimposing tSeries on anatomicals
%                           (i.e., anatFlag = 1). Values outside
%                           this range won't be overlaid. Default [0.2 1]
%   fps,[val]:              Set default frames per second of playback
%                           in the movie player.
%   saveAvi,[path]:         export the tSeries movie to an .AVI file
%                           in the specified path. If path is set to
%                           'dialog', will pop up a dialog to interactively
%                           set the target file. Note that because of
%                           codecs, exporting on Windows produces smaller
%                           files than on linux (at least in my
%                           experience). Uses movie2avi.
%
%
% More info on the mplay movie player for matlab:
%
% http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=3250&objectType=file
%
% Thanks Don Orofino!
%
% ras 12/04: wrote it.
% ras 03/05: added avi export.
% ras 11/06: added multiple slices option.
if notDefined('scans'),    scans = viewGet(vw, 'curScan'); end
if notDefined('anatFlag'), anatFlag = 0;                   end

%%%%% default params
slices = viewGet(vw, 'curSlice');
sDims  = viewGet(vw, 'sliceDims',scans(1));
if checkfields(vw,'ui','zoom')
    zoom = vw.ui.zoom;
    rsFactor = upSampleFactor(vw,scans(1));
    zoom(3:4) = zoom(3:4) ./ rsFactor(1);
    zoom(1:2) = zoom(1:2) ./ rsFactor(2);
    zoom = round(zoom);
    zoom(zoom==0) = 1;
else
    zoom = [1 sDims(2) 1 sDims(1)];
end
detrend         = 0;
convertToPct    = 0;
applyHisto      = 0; % use histogram-based thresholding like makeTSeriesMovie
compareFrame    = 0; % compare to a reference frame (if 0, don't compare)
meanThresh      = 0; % normalized raw intensity threshold
funcClip        = [0.2 1]; % clipping values for overlaying functionals, if anatFlag==1
autoCmap        = 1; % flag to use the full cmap range (if anatFlag==1)
movieFlag       = 1; % call the mplay GUI
fps             = 15;
saveAvi         = '';
viewType        = viewGet(vw,'viewType');
frames       = true; % true = use all frames from each scan; false = use 1st / last frame only

%%%%% parse the option flags
opts = unNestCell(varargin);  % to enable recursive passing of options
for i = 1:length(opts)
    if ischar(opts{i})
        switch lower(opts{i})
            case 'slices',          slices      = opts{i+1};
            case 'zoom',            zoom        = opts{i+1};
            case 'detrend',         detrend     = opts{i+1};
            case 'converttopct',    convertToPct= opts{i+1};
            case 'applyhisto',      applyHisto  = opts{i+1};
            case 'compareframes',   compareFrame= opts{i+1};
            case 'meanthresh',      meanThresh  = opts{i+1};
            case 'funcclip',        funcClip    = opts{i+1};
            case 'fps',             fps         = opts{i+1};
            case 'saveavi',         saveAvi     = opts{i+1};
            case 'frames',          frames      = opts{i+1};
            case 'autocmap',        autoCmap    = opts{i+1};
            case 'nomovie',         movieFlag   = 0;
        end
    end
end

%%%%% check view type
if isequal(viewType,'Flat')
    [M, movie] = flatLevelMovie(vw,scans,saveAvi);
    return
elseif isequal(viewType,'Volume') || isequal(viewType,'Gray')
    error('Volume/Gray support not provided yet.')
end

%%%%% init M intensity array
M = [];

%%%%% if >1 slices, run recursively
if length(slices) > 1
    [M, movie] = tSeriesMovieSlices(vw, scans, slices, anatFlag, opts);
    return
end

%%%%% get tSeries from selected scans
nFrames = 0;
hwait = mrvWaitbar(0,'Loading tSeries for movie...');
for scan = scans
    tSeries = loadtSeries(vw,scan,slices(1));
    if frames,  % get just the requested frames, throw out the rest
        % make sure we don't try to get more frames than exist (could
        % happen if the number of frames differs across scans)
        validframes = frames <= viewGet(vw, 'nframes', scan);
        % now restrict the tseries
        tSeries = tSeries(frames(validframes), :);
    end

    % update frame count
    nFrames = nFrames + size(tSeries,1);

    % detrend if selected
    if detrend ~= 0
        tSeries = detrendTSeries(tSeries, detrend, 20);
    end

    % append to M
    M = [M; tSeries];

    mrvWaitbar(find(scans==scan)/length(scans),hwait);
end
close(hwait);

%%%%% apply histogram criterion
if applyHisto==1
    M = histoThresh(M);
end

% if a mean intensity threshold is set, mask out stuff
% that doesn't have strong mean signal
if meanThresh > 0
    meanImg = mean(M,1);
    meanImg = rescale2(meanImg,[],[0 1000]);
    mask = (meanImg > 1000*meanThresh);
end

% convert to pct if selected
if convertToPct==1
    M = M ./ repmat(mean(M,1),[nFrames 1]);
end

% now apply mask if selected
if meanThresh > 0
    mask = repmat(mask,[nFrames 1]);
    M(mask==0) = min(M(:))-1;
end

%%%%% reshape and apply zoom
M = reshape(M', [sDims(1) sDims(2) nFrames]);
M = M(zoom(3):zoom(4),zoom(1):zoom(2),:);

%%%%% if anat flag set, convert M to 4D truecolor x time array
if anatFlag==1
	hwait = mrvWaitbar(0,'Superimposing over anat image...');

    % get cmap from map mode
    modeInfo    = viewGet(vw,'mapMode');
    numColors   = modeInfo.numColors;
    numGrays    = modeInfo.numGrays;
    cmap        = 255 .* modeInfo.cmap(numGrays+1:end,:);

    % make an overlay, indexing into the rows of cmap
    clipVals = min(M(:)) + funcClip .* (max(M(:))-min(M(:)));
    overlay  = uint8(rescale2(M, clipVals, [0 numColors]));
    if ~isequal(dataSize(vw,1),viewGet(vw,'Size'))

        % adjust zoom to be of view size again, not data size
        rsFactor  = upSampleFactor(vw,scans(1));
        zoom(3:4) = zoom(3:4) .* rsFactor(1);
        zoom(1:2) = zoom(1:2) .* rsFactor(2);

        % now up-sample the overlay to match the view resolution
        newSize = [zoom(4)-zoom(3) zoom(2)-zoom(1) nFrames];
        overlay = upSampleRep(overlay,newSize);

    end
    mrvWaitbar(1/4,hwait);

    % re-initialize M as a 4D array w/ the anat background
    anatImg = recomputeAnatImage(vw);
    anatImg = anatImg(zoom(3):zoom(4),zoom(1):zoom(2));
    anatImg = normalize(anatImg,0,255);
    M       = uint8(repmat(anatImg,[1 1 nFrames 3]));

    mrvWaitbar(1/2,hwait);

    % get the values for which the overlay > 0
    locs3D = find(overlay>0);
    [yy, xx, zz] = ind2sub(size(overlay),locs3D);

    % auto-scale the colormap
    if autoCmap==1
        overlay(locs3D) = rescale2(overlay(locs3D),[],[1 numColors]);
    end

    % for each color channel, plug in the appropriate
    % values from the cmap into the overlay locations
    for col = 1:3
        locs4D      = sub2ind(size(M),yy,xx,zz,repmat(col,size(yy)));
        M(locs4D)   = cmap(overlay(locs3D),col);
        mrvWaitbar(1/2+col/6,hwait);
    end
    close(hwait);

    % mplay likes the 4D truecolor array rows x cols x 3 x frames,
    % so permute to this order:
    M = permute(M,[1 2 4 3]);
elseif compareFrame > 0
    % make an R G B movie in which the red channel is the tSeries,
    % and the G,B channels are the reference frame
    ref = repmat(M(:,:,compareFrame), [1 1 nFrames 2]);
    M = cat(4, M, ref);
    M = permute(M, [1 2 4 3]);
    M = uint8(rescale2(M,[],[0 2^8-1]));
else
    % just reduce it to a uint8, for the purpose of rendering
    % (seems tricky to find a way around this bottleneck):
    M = uint8(rescale2(M,[],[0 2^8-1]));
end

%%%%% plug into mplayer
if movieFlag==1
    global mrSESSION dataTYPES;
    movie = implay(M,fps);

    ttltxt = sprintf('Movie: %s %s %s, Slice %i', ...
                        mrSESSION.sessionCode, ...
                        dataTYPES(vw.curDataType).name, ...
                        num2str(scans), slices(1));
    %set(get(movie,'hfig'),'Name',ttltxt);

    % play the movie
    %movie.play
end

%%%%% export to AVI if selected
if ~isempty(saveAvi)
    if isequal(saveAvi,'dialog')
        % put up a dialog to get path
        [f p] = uiputfile('*.avi','Save .avi movie as...');
        saveAvi = fullfile(p,f);
    end

    % check for .avi extension
    [~, p, ext] = fileparts(saveAvi);
    if ~isequal(lower(ext),'.avi')
        saveAvi = [saveAvi '.avi'];
    end

    % construct a matlab movie struct for the export
    mov = repmat(struct('cdata',[],'colormap',[]),1,nFrames);
    for frame = 1:nFrames
        if ndims(M)==4
            mov(frame).cdata = M(:,:,:,frame);
        else
            mov(frame).cdata = repmat(M(:,:,frame),[1 1 3]);
        end
    end

    % now export the movie
    if ispc
        codec = 'Indeo3';
    else
        codec = 'None';
    end
    movie2avi(mov,saveAvi,'Compression',codec,'Quality',100,'FPS',fps);
    fprintf('Exported movie to %s.\n',saveAvi);
end

return
% /-------------------------------------------------------------------/ %

% /-------------------------------------------------------------------/ %
function [M, movie] = tSeriesMovieSlices(vw, scans, slices, anatFlag, varargin)
% code to recursively call the main function (tSeriesMovie)
% for each slice, and assemble a montage across slices.
% Can be very memory-intensive.
M = []; movie = [];

%%%%% params
movieFlag = 1; % call the mplay GUI
fps = 15;
saveAvi = '';

% remove the 'slices' specification for recursion
opts = unNestCell(varargin);
% for i=1:length(opts), opts{i} = lower(opts{i}); end
ii = cellfind(opts, 'slices');
ok = setdiff(1:length(opts), [ii ii+1]);
opts = opts(ok);

% parse movie options
for i = 1:length(opts)
    if ischar(opts{i})
        switch lower(opts{i})
            case 'fps', fps = opts{i+1};
            case 'saveavi',
				% if there's a flag to save an .AVI file, we only want to
				% do that saving here, not in the recursive call to
				% tSeriesMovie below. So, mark it here, but removeit from
				% the options list
				saveAvi = opts{i+1};
				opts{i} = []; opts{i+1} = [];
            case 'nomovie', movieFlag = 0;
        end
    end
end

% get movie data for each slice
hwait = mrvWaitbar(0, 'Creating Montage Across Slices...');
set(hwait, 'Position', get(hwait, 'Position')+[0 100 0 0]);

nSlices = length(slices);
for i = 1:nSlices
    slice{i} = tSeriesMovie(vw, scans, anatFlag, 'slices', slices(i), ...
                            'nomovie', opts);
    mrvWaitbar(.8 * (i/nSlices), hwait);
end

%%%%%% make montage
nRows = ceil(sqrt(nSlices));
nCols = ceil(nSlices/nRows);
if nSlices < nRows*nCols    % pad out extra images
    for i = nSlices+1:nRows*nCols
        slice{i} = zeros(size(slice{1}));
    end
end
slice = reshape(slice, [nCols nRows])';

mrvWaitbar(.95, hwait);

for r = 1:nRows
    row = [];
    for c = 1:nCols
        row = cat(2, row, slice{1,c});
    end

    % be memory-efficient: clear the row we just grabbed
    slice = slice(2:end,:);

    M = cat(1, M, row);
end

close(hwait);

%%%%% show movie if selected
if movieFlag==1
    global mrSESSION dataTYPES;
    movie = mplay(M,fps);
    ttltxt = sprintf( 'Movie: %s %s %s, Slices %s', ...
                       mrSESSION.sessionCode, ...
                       dataTYPES(vw.curDataType).name, ...
                       num2str(scans), num2str(slices) );
    set(get(movie,'hfig'), 'Name', ttltxt);
end

%%%%% export to AVI if selected
if ~isempty(saveAvi)
    if isequal(saveAvi,'dialog')
        % put up a dialog to get path
        [f p] = uiputfile('*.avi','Save .avi movie as...');
        saveAvi = fullfile(p,f);
    end

    % check for .avi extension
    [f p ext] = fileparts(saveAvi);
    if ~isequal(lower(ext),'.avi')
        saveAvi = [saveAvi '.avi'];
    end

    % construct a matlab movie struct for the export
    nFrames = size(M, ndims(M));
    mov = repmat( struct('cdata',[], 'colormap',[]), 1, nFrames );
    for frame = 1:nFrames
        if ndims(M)==4
            mov(frame).cdata = M(:,:,:,frame);
        else
            mov(frame).cdata = repmat(M(:,:,frame),[1 1 3]);
        end
	end

	% allow the movie path to refer to directories that have not yet been
	% created (like 'Movies/')
	ensureDirExists( fileparts(fullpath(saveAvi)) );

    % now export the movie
    if ispc
        codec = 'Indeo5';
    else
        codec = 'None';
    end
    movie2avi(mov, saveAvi, 'Compression', codec, 'Quality', 100, 'FPS', fps);
    fprintf('Exported movie to %s.\n', saveAvi);
end

%% lastly, play the movie if it exists
if movieFlag==1
    movie.play
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%--- help for set_simulation_initial_conditions ---
%
%  Set the initial conditions for forecasting, simulation and irfs
%
%  ::
%
%    Initcond=set_simulation_initial_conditions(obj)
%
%  Args:
%
%     obj (rise | dsge | svar | rfvar): model object
%
%  Returns:
%     :
%
%     - **Initcond** [struct]: Initial conditions for simulation
%
%  Note:
%
%     - if future values of endogenous are found, they are discarded the
%       variables are not explicitly declared as conditioning endogenous
%       variables.
%
%     - if either endogenous or exogenous variables are declared as
%       conditioning variables, we automatically have a conditional forecasting
%       exercise.
%
%     - if no variable is declared as conditioning, we have a simulation
%       exercise, in which case, there is a burn-in period.
%
%     - Initialization is always done at the steady state and the declared
%       endogenous variables merely override the steady state values. This means
%       that if a variable is not found in the database, it is initialized at its
%       steady state.
%
%  Example:
%
%  See also:
%

author: jmaih


text:
function [lat, long ] = gcj2wgs(gcjLat, gcjLng)
% Not exact
% Fully vectorized

    lat = gcjLat;
    long = gcjLng;

    inChina = ~outOfChina(gcjLat, gcjLng);

    [dlat, dlng] = delta(gcjLat(inChina), gcjLng(inChina));

    lat(inChina) = gcjLat(inChina) - dlat;
    long(inChina) = gcjLng(inChina) - dlng;

end

author: googollee

repo: eviltransform


text:
function issquare=issquare(X)
%ISSQUARE Check if variable is square

n = X.dim(1);
m = X.dim(2);
issquare = (n==m);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%function dtiGenerateNormalizedData
% Generate a spatially-normalized set of summary data for two groups (boys
% and girls).

bd = '/biac3/wandell4/data/reading_longitude/dti_y1234/';
[sd,sc] = findSubjects([bd '*'],'dti06trilinrt');
behaveData = dtiGetBehavioralDataStruct(sc, '/biac3/wandell4/data/reading_longitude/read_behav_measures_longitude.csv', sd);
clear sd sc;
outDir = '/biac3/wandell4/data/reading_longitude/logNorm_analysis';
sumFile = fullfile(outDir,'sum_090715');

n = numel(behaveData);
templateName = 'SIRL54';
templateDir = fullfile(fileparts(which('mrDiffusion.m')),'templates');
template = fullfile(templateDir, [templateName '_EPI.nii.gz']);
t1Template = fullfile(templateDir, [templateName '_T1.nii.gz']);
desc = ['Normalized to ' templateName ' using B0.'];
clear templateName templateDir;
mT1 = 0;
mB0 = 0;
mBm = 0;
mDt6 = 0;
clear sn t1Sn;
for(ii=116:n)
    tic;
    fprintf('Processing %d of %d (%s)...\n',ii, n, behaveData(ii).dataDir);
    dtiRawFixDt6File(behaveData(ii).dataDir);
    [dt,t1] = dtiLoadDt6(behaveData(ii).dataDir);
    im = mrAnatHistogramClip(double(dt.b0),0.4,0.98);
    dt.b0 = im;
    im(bwareaopen(im==0,10000,6)) = NaN;
    evalc('sn{ii} = mrAnatComputeSpmSpatialNorm(im, dt.xformToAcpc, template);');
    evalc('dt_sn = dtiSpmDeformer(dt,sn{ii});');
    sn{ii}.VG.dat = [];
    dtXformToAcpc = dt_sn.xformToAcpc;
    adcUnits = dt_sn.adcUnits;
    mB0 = mB0 + dt_sn.b0;
    % We'll exclude non-PD tensors by removing them from the brain mask
    [vec,val] = dtiEig(dt_sn.dt6);
    mBm = mBm + double(dt_sn.brainMask & all(val>0,4));
    mDt6 = mDt6 + dt_sn.dt6;
    dt6(:,:,:,:,ii) = single(dt_sn.dt6);
    b0(:,:,:,ii) = single(dt_sn.b0);
    clear dt im dt_sn vec val;
    im = mrAnatHistogramClip(double(t1.img),0.4,0.98);
    evalc('t1Sn{ii} = mrAnatComputeSpmSpatialNorm(im, t1.xformToAcpc, t1Template);');
    t1Sn{ii}.VG.dat = [];
    [im,t1XformToAcpc] = mrAnatResliceSpm(im, t1Sn{ii}, mrAnatXformCoords(t1Sn{ii}.VG.mat,[1 1 1; t1Sn{ii}.VG.dim]), [1 1 1], 7, false);
    % mrAnatOverlayMontage(dt_sn.b0,dt_sn.xformToAcpc, im, t1Sn.VG.mat, autumn(256), [.8 1], [-20:2:60])
    mT1 = mT1 + im;
    clear im t1;
    toc
end
tmpFile = [tempname '.mat'];
save(tmpFile,'-v7.3');

clear vec val im ii tmpFile;
brainMask = mBm>=ceil(n*.95) & all(dt6(:,:,:,1,:)~=0,5);
b0 = dtiImgToInd(b0,brainMask);
dt6 = dtiImgToInd(dt6, brainMask);
mBm = single(mBm);
mB0 = single(mB0);
mDt6 = single(mDt6);
mT1 = single(mT1);
clear ans;

save(sumFile);

%return

dtiWriteTensorsToNifti(mDt6, dtXformToAcpc, desc, adcUnits, fullfile(outDir,'dti06','bin','tensors.nii.gz'));
dtiWriteNiftiWrapper(mB0, dtXformToAcpc, fullfile(outDir,'dti06','bin','b0.nii.gz'), 1, desc);
dtiWriteNiftiWrapper(uint8(brainMask), dtXformToAcpc, fullfile(outDir,'dti06','bin','brainMask.nii.gz'), 1, desc);
dtiWriteNiftiWrapper(mT1, t1XformToAcpc, fullfile(outDir,'dti06','t1.nii.gz'), 1, desc);


text:
function I = randpermNK(N,K)
% N is the number of trials/subjects
% K is the number of groups
if K==1
    I = cell(1);
    I{1} = randperm(N);
    return;
end
np=(N-mod(N,K))/K; % number of elements per group
[~,idx]=sort(rand(N,1));
i=1;
j=1;
I={};
while 1
    I{j}=idx(i:i+np-1,1);
    if N-(i+np)+1 < np
        I{j} = [I{j}; idx(i+np:end,1)];
        break
    end
    i=i+np;
    j=j+1;
end
end

author: OHBA-analysis

repo: HMM-MAR


text:
function [xy,scalFact]=ellips2Mercator(latLon,a,f)
%%ELLIPS2MERCATOR Convert from ellipsoidal (geodetic) latitude and
%               longitude into Mercator Easting and Northing. The Mercator
%               projection is useful, because lines of constant heading
%               (i.e. degrees East of North) are straight lines on the map
%               and the projection is conformal. However, the poles are
%               placed at +/-Inf. The Mercator projection is a type of
%               cylindrical projection.
%
%INPUTS: latLon A 2XN matrix of ellipsoidal [latitude;longitude] points
%               given in radians.
%             a The semi-major axis of the reference ellipsoid. If this
%               argument is omitted, the value in
%               Constants.WGS84SemiMajorAxis is used.
%             f The flattening factor of the reference ellipsoid. If this
%               argument is omitted, the value in Constants.WGS84Flattening
%               is used.
%
%OUTPUTS: xy A 2XN matrix of [Easting;Northing] points in the Mercator
%            projection. Easting is related to longitude and Northing is
%            related to latitude. The units are those of the semi-major
%            axis and are typically meters.
%   scalFact The scale factor (both parallel and meridian) at the given
%            points.
%
%The Mercator projection is discussed in detail in Chapter 7 of [1].
%Equations for the conversion along with a commonly used but very
%inaccurate spherical Mercator conversion (the "Web" Mercator map) are
%given in [2]. The "web" Mercator conversion can be obtained by setting
%f to zero.
%
%REFERENCES:
%[1] J. P. Snyder, "Map projections- a working manual," U.S. Geological
%    Survey, Tech. Rep. 1395, 1987.
%[2] Office of Geomatics, "National geospatial-intelligence agency
%   standardization document: Implementation practice: Web mercator map
%   projection," National Geospatial-Intelligence Agency, Tech. Rep.
%   NGA.SIG.0011 1.0.0 WEBMERC, 18 Feb. 2014. [Online]. Available:
%   http://earth-info.nga.mil/GandG/wgs84/web_mercator/(U)%20NGA_SIG_0011_1.0.0_WEBMERC.pdf
%
%August 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(a))
    a=Constants.WGS84SemiMajorAxis;
end

if(nargin<3||isempty(f))
    f=Constants.WGS84Flattening;
end

e2=2*f-f^2;
%The first numerical eccentricity of the ellipsoid.
e=sqrt(e2);

sinPhi=sin(latLon(1,:));
cosPhi=cos(latLon(1,:));
lambda=latLon(2,:);

x=a*lambda;
y=a*(atanh(sinPhi)-e*atanh(e*sinPhi));

xy=[x;y];

scalFact=sqrt(1-e2*sinPhi.^2)./cosPhi;

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function varargout = dim(varargin)
% file_array's dimension property
% For getting the value
% dat = dim(obj)
%
% For setting the value
% obj = dim(obj,dat)
%__________________________________________________________________________
% Copyright (C) 2005-2017 Wellcome Trust Centre for Neuroimaging

%
% $Id: dim.m 7147 2017-08-03 14:07:01Z spm $

if nargin==2
    varargout{1} = asgn(varargin{:});
elseif nargin==1
    varargout{1} = ref(varargin{:});
else
    error('Wrong number of arguments.');
end

%==========================================================================
% function dat = ref(obj)
%==========================================================================
function dat = ref(obj)
dat = obj.dim;

%==========================================================================
% function obj = asgn(obj,dat)
%==========================================================================
function obj = asgn(obj,dat)
if isnumeric(dat) && all(dat>=0) && all(rem(dat,1)==0)
    dat = [double(dat(:)') 1 1];
    lim = max([2 find(dat~=1)]);
    dat = dat(1:lim);
    obj.dim = dat;
else
    error('"dim" must be a vector of positive integers.');
end

author: spm


text:
function boundary_display ( v, bi )

%% BOUNDARY_DISPLAY displays a boundary.
%
%  Discussion:
%
%    The boundary is assumed to consist of one or more closed curves.
%
%    Each curve is represented by a sequence of line segments, traced in
%    counterclockwise order. (Interior holes go in clockwise order.)
%
%  Parameters:
%
%    Input, real V(V_NUM,2), the coordinates of vertices used to define the curves.
%
%    Input, integer BI(BI_NUM), a sequence of indices into V.  Each closed curve
%    is defined by giving a sequence of indices, which is terminated by
%    repeating the starting index.  Thus, BI = { 3, 1, 5, 3, 4, 2, 9, 7, 4 }
%    describes two curves: ( 3, 1, 5, 3 ) and (4, 2, 9, 7, 4 ).
%
  bi_num = length ( bi );

  clf
  hold on;
  next = 1;
  s = bi(1);
  t2 = s;
  draw = 1;

  while ( next < bi_num )
    t1 = t2;
    next = next + 1;
    t2 = bi(next);
    if ( draw )
      line ( [ v(t1,1), v(t2,1) ], [ v(t1,2), v(t2,2) ], 'LineWidth', 2, 'Color', 'k' );
      if ( t2 == s )
        draw = 0;
      end
    else
      s = t2;
      draw = 1;
    end

  end

  plot ( v(:,1), v(:,2), 'r.', 'MarkerSize', 15 );

  axis equal
  grid on
  hold off

  return
end


text:
function chEOG=identify_eog_channels(fn,x);
% IDENTIFY_EOG_CHANNELS returns bipolar EOG channels for
%  correcting of EOG artifacts using regression analysis
%
%  EOGchan = IDENTIFY_EOG_CHANNELS(...)
%
% EOGchan is a sparse matrix of size number_of_channels x 2.
% The sparsity ensures that missing samples of unrelated channels
% do not affect the data.
%
%  [...] = IDENTIFY_EOG_CHANNELS(filename)
%  [...] = IDENTIFY_EOG_CHANNELS(HDR)
%	filename or HDR struct can be used
%  [...] = IDENTIFY_EOG_CHANNELS(...,'x')
%     looks for EOG channels whos Label start with x
%
% see also: GET_REGRESS_EOG, SLOAD
%
% Reference(s):
% [1] Schlogl A, Keinrath C, Zimmermann D, Scherer R, Leeb R, Pfurtscheller G.
%	A fully automated correction method of EOG artifacts in EEG recordings.
%	Clin Neurophysiol. 2007 Jan;118(1):98-104. Epub 2006 Nov 7.
% 	http://dx.doi.org/10.1016/j.clinph.2006.09.003
%       http://pub.ist.ac.at/~schloegl/publications/schloegl2007eog.pdf

%	$Id: identify_eog_channels.m 2649 2011-03-09 09:52:44Z schloegl $
%	Copyright (C) 2006,2007,2009,2010 by Alois Schloegl
%    	This is part of the BIOSIG-toolbox http://biosig.sf.net/

% Biosig is free software; you can redistribute it and/or
% modify it under the terms of the GNU Library General Public
% License as published by the Free Software Foundation; either
% Version 3 of the License, or (at your option) any later version.
%
% This library 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
% Library General Public License for more details.
%
% You should have received a copy of the GNU Library General Public
% License along with this library; if not, write to the
% Free Software Foundation, Inc., 59 Temple Place - Suite 330,
% Boston, MA  02111-1307, USA.

if ischar(fn),
	HDR=sopen(fn);
	HDR=sclose(HDR);
elseif isstruct(fn)
	HDR=fn;
end;

% graz
g1 = strmatch('EOG-left',HDR.Label);
g2 = strmatch('EOG-central',HDR.Label);
g3 = strmatch('EOG-right',HDR.Label);
if isempty([g1,g2,g3])
	g1 = strmatch('EOG:ch01',HDR.Label);
	g2 = strmatch('EOG:ch02',HDR.Label);
	g3 = strmatch('EOG:ch03',HDR.Label);
end;

% berlin
if nargin<2,
	v1 = strmatch('eogv1',lower(HDR.Label));
	v2 = strmatch('eogv2',lower(HDR.Label));
	v0 = strmatch('eogv',lower(HDR.Label));
	v3 = strmatch('eogvp',lower(HDR.Label));
	v4 = strmatch('eogvn',lower(HDR.Label));

	h1 = strmatch('eogh1',lower(HDR.Label));
	h2 = strmatch('eogh2',lower(HDR.Label));
	h0 = strmatch('eogh' ,lower(HDR.Label));
	h3 = strmatch('eoghp',lower(HDR.Label));
	h4 = strmatch('eoghn',lower(HDR.Label));
else
	v1 = [];
	v2 = [];
	v0 = [];
	v3 = strmatch('xeogvp',lower(HDR.Label));
	v4 = strmatch('xeogvn',lower(HDR.Label));

	h1 = [];
	h2 = [];
	h0 = [];
	h3 = strmatch('xeoghp',lower(HDR.Label));
	h4 = strmatch('xeoghn',lower(HDR.Label));
end;

g = [g1;g2;g3];
v = [v1,v2,v3,v4];
if isempty(v), v=v0; end;
h = [h1,h2,h3,h4];
if isempty(h), h=h0; end;
if length(g)==3,
	chEOG = sparse([g1,g2,g2,g3],[1,1,2,2],[1,-1,1,-1],HDR.NS,2);
elseif length(g)==2,
	chEOG = sparse([g1,g2,g3],[1,1],[1,-1],HDR.NS,1);
else
	c = (length(v)>0);
	sz2 = (length(v)>0) + (length(h)>0);
	if length(v)==1,
		chEOG = sparse(v,c,1,HDR.NS,sz2);
	elseif length(v)==2,
		chEOG = sparse(v,[c,c],[1,-1],HDR.NS,sz2);
	else
		chEOG = 0;
	end;
	if length(h)==1,
		chEOG = chEOG+sparse(h,1+c,1,HDR.NS,1+c);
	elseif length(h)==2,
		chEOG = chEOG+sparse(h,[1,1]+c,[1,-1],HDR.NS,1+c);
	end;
end;

if size(chEOG,2)<2,
	warning('EOG channels are missing, or were not recognized');
end;

author: PatternRecognition

repo: OpenBMI


text:
function varargout = showufmf(varargin)
% showufmf
%
% Opens a GUI to show UFMF-compressed video
% Optional inputs:
%
% UFMFName: Name of input video.
%
%

% Edit the above text to modify the response to help showufmf

% Last Modified by GUIDE v2.5 18-Jul-2014 11:42:42

% Begin initialization code - DO NOT EDIT
gui_Singleton = 0;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @showufmf_OpeningFcn, ...
                   'gui_OutputFcn',  @showufmf_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1}) && exist(varargin{1},'file'),
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before showufmf is made visible.
function showufmf_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to showufmf (see VARARGIN)

% Choose default command line output for showufmf
handles.output = hObject;

handles.MAXSLIDERRES = 10^6;

handles.figpos = get(handles.figure1,'Position');

% set up path
if isempty(which('myparse')),
  if exist('../misc','file'),
    addpath('../misc');
  end
  while isempty(which('myparse')),
    miscdir = uigetdir('.','Choose "misc" folder to add to path');
    try
      addpath(miscdir);
    catch %#ok<CTCH>
    end
  end
end

% load previous values
handles.rcfilename = '.showufmfrc.mat';
handles.previous_values = struct;
if exist(handles.rcfilename,'file')
  handles.previous_values = load(handles.rcfilename);
end

if isfield(handles.previous_values,'filename'),
  handles.filename = handles.previous_values.filename;
  if ~exist(handles.filename,'file'),
    handles.filename = '';
  end
else
  handles.filename = '';
end

handles.menu_Colormaps = get(handles.menu_Colormap,'Children');
if isfield(handles.previous_values,'Colormap'),
  handles.Colormap = handles.previous_values.Colormap;
else
  handles.Colormap = 'gray';
end

for h = handles.menu_Colormaps',
  set(h,'Checked','off');
end
switch handles.Colormap,
  case 'gray',
    set(handles.menu_ColormapGray,'Checked','on');
  case 'jet',
    set(handles.menu_ColormapJet,'Checked','on');
end

if isfield(handles.previous_values,'MaxFPS'),
  handles.MaxFPS = handles.previous_values.MaxFPS;
else
  handles.MaxFPS = 0;
  handles.MinSPF = 0;
end
if isfield(handles.previous_values,'BackSubThresh'),
  handles.BackSubThresh = handles.previous_values.BackSubThresh;
else
  handles.BackSubThresh = 10;
end

[handles.filename,handles.MaxFPS,handles.FMFName,handles.BackSubThresh,figpos,handles.FirstFrame] = myparse(varargin,...
  'UFMFName',handles.filename,'MaxFPS',handles.MaxFPS,'FMFName','','BackSubThresh',handles.BackSubThresh,...
  'FigPos',[],'FirstFrame',1);
if ~isempty(figpos),
  % nan means auto-set
  figpos(isnan(figpos)) = handles.figpos(isnan(figpos));
  handles.figpos = figpos;
  set(handles.figure1,'Units','Pixels','Position',figpos);
end
if handles.MaxFPS == 0,
  handles.MinSPF = 0;
else
  handles.MinSPF = 1/handles.MaxFPS;
end
handles.UFMFNameIsInput = ismember('UFMFName',varargin(1:2:end));
handles.FMFNameIsInput = ~isempty(handles.FMFName) && exist(handles.FMFName,'file');

s1 = {'Compressed Frame'
  'Background Model'
  'Background Difference'
  'Is Foreground'};
s2 = {'Raw Frame'
  'Compression Error'
  'Thresholded Compression Error'};
if handles.FMFNameIsInput,
  set(handles.popupmenu_Show,'String',[s1;s2]);
else
  set(handles.popupmenu_Show,'String',s1);
end

% set slider listener
fcn = get(handles.slider_Frame,'Callback');
if verLessThan('matlab','8.4.0')
  handles.hslider_listener = handle.listener(handles.slider_Frame,...
    'ActionEvent',fcn);
  set(handles.slider_Frame,'Callback','');
else
  handles.hslider_listener = addlistener(handles.slider_Frame,...
    'ContinuousValueChange',fcn);
  set(handles.slider_Frame,'Callback','');
end

%handles.slider_listener = handle.listener(handles.slider_Frame,'ActionEvent',fcn);

% open video
handles = open_fmf(handles);

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes showufmf wait for user response (see UIRESUME)
% uiwait(handles.figure1);

function play(hObject)

handles = guidata(hObject);
set(hObject,'String','Stop','BackgroundColor',[.5,0,0]);
handles.isplaying = true;
guidata(hObject,handles);
tic;
t0 = now;
fnext = handles.f;

while true,
  handles = guidata(hObject);
  if ~handles.isplaying,
    break;
  end
  f = fnext;
  handles.f = f;
  handles = update_frame(handles);
  guidata(hObject,handles);
  if handles.MaxFPS > 0,
    tmp = toc;
    t1 = now;
    if tmp < handles.MinSPF
      pause(handles.MinSPF - t1);
      fnext = f + 1;
    elseif tmp > handles.MinSPF,
      fnext = max(1,round((t1-t0)*24*3600/handles.MinSPF));
      %fprintf('setting fnext to %d\n',fnext);
      if fnext > handles.nframes,
        break;
    end
      drawnow;
  else
    drawnow;
  end
  else
    fnext = f+1;
    drawnow;
  end

%   if handles.MaxFPS > 0,
%     tmp = toc;
%     if tmp < handles.MinSPF
%       pause(handles.MinSPF - tmp);
%     end
%   else
%     drawnow;
%   end
  tic;
end

stop(hObject);

function stop(hObject)

handles = guidata(hObject);
handles.isplaying = false;
set(hObject,'String','Play','BackgroundColor',[0,.5,0]);
guidata(hObject,handles);

% --- Outputs from this function are returned to the command line.
function varargout = showufmf_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

function handles = update_frame(handles)

try
  [handles.im,handles.timestamp,handles.cachedidx,handles.bb,handles.mu] = handles.readframe(handles.f);
catch ME
  warning([sprintf('Error reading frame %d\n',handles.f),getReport(ME)]);
  return;
end
s = get(handles.popupmenu_Show,'String');
v = get(handles.popupmenu_Show,'Value');
showType = s{v};
switch showType,
  case 'Compressed Frame',
    set(handles.himage,'CData',handles.im);
    set(handles.axes_Video,'CLim',[0,255]);
  case 'Background Model',
    set(handles.himage,'CData',handles.mu);
    set(handles.axes_Video,'CLim',[0,255]);
  case 'Background Difference',
    if handles.FMFNameIsInput,
      im = handles.readframe_raw(handles.f);
      if size(im,3) > 1,
        im = rgb2gray(im);
      end
      im = cast(im,handles.headerinfo.dataclass);
    else
      im = handles.im;
    end
    diff = max(im - handles.mu,handles.mu - im);
    set(handles.himage,'CData',diff);
    set(handles.axes_Video,'CLim',[0,max(1,max(diff(:)))]);
  case 'Is Foreground',
    if handles.FMFNameIsInput,
      im = handles.readframe_raw(handles.f);
      if size(im,3) > 1,
        im = rgb2gray(im);
      end
      im = cast(im,handles.headerinfo.dataclass);
    else
      im = handles.im;
    end
    if size(im,3) > 1,
      im = rgb2gray(im);
    end
    diff = max(im - handles.mu,handles.mu - im);
    isfore = diff >= handles.BackSubThresh;
    set(handles.himage,'CData',isfore);
    set(handles.axes_Video,'CLim',[0,1]);
  case 'Raw Frame',
    im = handles.readframe_raw(handles.f);
    if size(im,3) > 1,
      im = rgb2gray(im);
    end
    im = cast(im,handles.headerinfo.dataclass);
    handles.rawim = im;
    set(handles.himage,'CData',handles.rawim);
    set(handles.axes_Video,'CLim',[0,255]);
  case 'Compression Error',
    im = handles.readframe_raw(handles.f);
    if size(im,3) > 1,
      im = rgb2gray(im);
    end
    im = cast(im,handles.headerinfo.dataclass);
    handles.rawim = im;
    diff = max(handles.im - handles.rawim,handles.rawim - handles.im);
    set(handles.himage,'CData',diff);
    set(handles.axes_Video,'CLim',[0,max(1,max(diff(:)))]);
  case 'Thresholded Compression Error',
    im = handles.readframe_raw(handles.f);
    if size(im,3) > 1,
      im = rgb2gray(im);
    end
    im = cast(im,handles.headerinfo.dataclass);
    handles.rawim = im;
    diff = max(handles.im - handles.rawim,handles.rawim - handles.im);
    isfore = diff >= handles.BackSubThresh;
    set(handles.himage,'CData',isfore);
    set(handles.axes_Video,'CLim',[0,1]);
end
if(handles.headerinfo.is_fixed_size)
  y = cat(2,handles.bb(:,1),...
    handles.bb(:,1),...
    handles.bb(:,1)+handles.headerinfo.max_width,...
    handles.bb(:,1)+handles.headerinfo.max_width,...
    handles.bb(:,1),...
    nan(size(handles.bb,1),1))'-.5;
  x = cat(2,handles.bb(:,2),...
    handles.bb(:,2)+handles.headerinfo.max_height,...
    handles.bb(:,2)+handles.headerinfo.max_height,...
    handles.bb(:,2),...
    handles.bb(:,2),...
    nan(size(handles.bb,1),1))'-.5;
else
  y = cat(2,handles.bb(:,1),...
    handles.bb(:,1),...
    handles.bb(:,1)+handles.bb(:,3),...
    handles.bb(:,1)+handles.bb(:,3),...
    handles.bb(:,1),...
    nan(size(handles.bb,1),1))'-.5;
  x = cat(2,handles.bb(:,2),...
    handles.bb(:,2)+handles.bb(:,4),...
    handles.bb(:,2)+handles.bb(:,4),...
    handles.bb(:,2),...
    handles.bb(:,2),...
    nan(size(handles.bb,1),1))'-.5;
end
set(handles.hboxes,'xdata',x(:),'ydata',y(:));
set(handles.edit_Frame,'String',num2str(handles.f));
set(handles.slider_Frame,'Value',handles.f);

% --- Executes on slider movement.
function slider_Frame_Callback(hObject, eventdata, handles)
% hObject    handle to slider_Frame (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
v = get(hObject,'Value');
roundv = round(v);
if roundv == handles.f,
  if v > handles.f,
    handles.f = min(handles.nframes,handles.f+1);
  elseif v < handles.f,
    handles.f = max(1,handles.f-1);
  end
else
  handles.f = roundv;
end

handles = update_frame(handles);
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function slider_Frame_CreateFcn(hObject, eventdata, handles)
% hObject    handle to slider_Frame (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on button press in pushbutton_PlayStop.
function pushbutton_PlayStop_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton_PlayStop (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if ~handles.isplaying,
  play(hObject);
else
  stop(hObject);
end

function edit_Frame_Callback(hObject, eventdata, handles)
% hObject    handle to edit_Frame (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit_Frame as text
%        str2double(get(hObject,'String')) returns contents of edit_Frame as a double
f = str2double(get(hObject,'String'));
if isnan(f),
  set(hObject,'String',num2str(handles.f));
  return;
end
handles.f = max(1,min(f,handles.nframes));
handles = update_frame(handles);
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function edit_Frame_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit_Frame (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --------------------------------------------------------------------
function menu_File_Callback(hObject, eventdata, handles)
% hObject    handle to menu_File (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function menu_File_Open_Callback(hObject, eventdata, handles)
% hObject    handle to menu_File_Open (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

handles = open_fmf(handles);
guidata(hObject,handles);

function handles = open_fmf(handles)

if isfield(handles,'isplaying') && handles.isplaying,
  stop(handles.figure1);
  handles = guidata(handles.figure1);
end

if ~handles.UFMFNameIsInput,
  handles.filterspec = {'*.ufmf','MicroFlyMovieFormat (*.ufmf)'};

  [filename, pathname] = uigetfile(handles.filterspec, 'Choose UFMF video to play',handles.filename);
  if ~ischar(filename),
    return;
  end
  handles.filename = fullfile(pathname,filename);
else
  handles.UFMFNameIsInput = false;
end

if isfield(handles,'fid') && ~isempty(fopen(handles.fid)),
  fclose(handles.fid);
end
if isfield(handles,'fid_raw') && ~isempty(fopen(handles.fid_raw)),
  fclose(handles.fid_raw);
end
if isfield(handles,'himage') && ishandle(handles.himage),
  delete(handles.himage);
end

try
  [handles.readframe,handles.nframes,handles.fid,handles.headerinfo] = ...
    get_readframe_fcn(handles.filename);
catch ME
  s = sprintf('Could not read video %s.',handles.filename);
  uiwait(errordlg(s,'Error opening video'));
  rethrow(ME);
end

if handles.FMFNameIsInput,
  [handles.readframe_raw,handles.nframes_raw,handles.fid_raw,handles.headerinfo_raw] = ...
    get_readframe_fcn(handles.FMFName);
end

% set slider steps
% this seems to be the limit to slider step resolution
handles.stepsize = ceil(handles.nframes/handles.MAXSLIDERRES);
step1 = handles.stepsize/(handles.nframes-1);
sliderstep = [step1,min(1,100*step1)];
set(handles.slider_Frame,'Value',0,'SliderStep',sliderstep,'Min',1,'Max',handles.nframes);

% show first image
if ischar(handles.FirstFrame) && strcmpi(handles.FirstFrame,'maxnboxes'),
  nboxes = ufmf_read_nboxes(handles.headerinfo,1:handles.nframes);
  [~,handles.FirstFrame] = max(nboxes);
end
handles.f = handles.FirstFrame;
[handles.im,handles.timestamp,handles.cachedidx,handles.bb,handles.mu] = handles.readframe(handles.f);
hold(handles.axes_Video,'off');
if size(handles.im,3) == 1,
  handles.himage = imagesc(handles.im,'Parent',handles.axes_Video,[0,255]);
else
  handles.himage = image(uint8(handles.im),'Parent',handles.axes_Video);
end
colormap(handles.axes_Video,handles.Colormap);
axis(handles.axes_Video,'image','off');
axis(handles.axes_Video,[0,size(handles.im,2)+1,0,size(handles.im,1)+1]);
hold(handles.axes_Video,'on');
handles.hboxes = plot(handles.axes_Video,nan(1,2),nan(1,2),'r-');

handles = update_frame(handles);

handles.isplaying = false;

% --------------------------------------------------------------------
function menu_File_Quit_Callback(hObject, eventdata, handles)
% hObject    handle to menu_File_Quit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
figure1_CloseRequestFcn(handles.figure1, eventdata, handles);

% --- Executes when user attempts to close figure1.
function figure1_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: delete(hObject) closes the figure
if isfield(handles,'fid') && ~isempty(fopen(handles.fid)),
  fclose(handles.fid);
end
if isfield(handles,'fid_raw') && ~isempty(fopen(handles.fid_raw)),
  fclose(handles.fid_raw);
end
savefns = {'filename','FMFName','MaxFPS','BackSubThresh','Colormap'};
save(handles.rcfilename,'-struct','handles',savefns{:});
delete(hObject);

% --- Executes when figure1 is resized.
function figure1_ResizeFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if ~isfield(handles,'figpos'),
  return;
end

sliderpos = get(handles.slider_Frame,'Position');
editpos = get(handles.edit_Frame,'Position');
pushbuttonpos = get(handles.pushbutton_PlayStop,'Position');
axespos = get(handles.axes_Video,'Position');
newfigpos = get(handles.figure1,'Position');
popupmenupos = get(handles.popupmenu_Show,'Position');

woff = newfigpos(3) - handles.figpos(3);
hoff = newfigpos(4) - handles.figpos(4);

sliderpos(3) = sliderpos(3) + woff;
editpos(1) = editpos(1) + woff;
pushbuttonpos(1) = pushbuttonpos(1) + woff;
popupmenupos(1) = popupmenupos(1) + woff;

axespos(3) = axespos(3) + woff;
axespos(4) = axespos(4) + hoff;

set(handles.slider_Frame,'Position',sliderpos);
set(handles.edit_Frame,'Position',editpos);
set(handles.pushbutton_PlayStop,'Position',pushbuttonpos);
set(handles.axes_Video,'Position',axespos);
set(handles.popupmenu_Show,'Position',popupmenupos);

handles.figpos = newfigpos;
guidata(hObject,handles);

% --------------------------------------------------------------------
function menu_Edit_Preferences_Callback(hObject, eventdata, handles)
% hObject    handle to menu_Edit_Preferences (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

MAXFPS_IDX = 1;
BACKSUBTHRESH_IDX = 2;
prompts = {};
defAns = {};
prompts{MAXFPS_IDX} = 'Max FPS: ';
defAns{MAXFPS_IDX} = num2str(handles.MaxFPS);
prompts{BACKSUBTHRESH_IDX} = 'BackSub Thresh: ';
defAns{BACKSUBTHRESH_IDX} = num2str(handles.BackSubThresh);
while true,
  answer = inputdlg(prompts,'showufmf Preferences',1,defAns,'on');
  if isempty(answer),
    return;
  end
  v = str2double(answer{MAXFPS_IDX});
  iserror = false;
  if isnan(v),
    iserror = true;
  else
    defAns{MAXFPS_IDX} = num2str(v);
    handles.MaxFPS = v;
    handles.MinSPF = 1 / handles.MaxFPS;
  end
  v = str2double(answer{BACKSUBTHRESH_IDX});
  if isnan(v),
    iserror = true;
  else
    defAns{BACKSUBTHRESH_IDX} = num2str(v);
    handles.BackSubThresh = v;
  end

  if iserror,
    uiwait(warndlg('Illegal values entered. Max FPS and BackSub Thresh must be numbers','Bad Preferences'));
  else
    break;
  end
end
guidata(hObject,handles);

% --------------------------------------------------------------------
function menu_Help_About_Callback(hObject, eventdata, handles)
% hObject    handle to menu_Help_About (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

s = {};
s{end+1} = 'ShowUFMF';
s{end+1} = '';
s{end+1} = 'Kristin Branson';
s{end+1} = 'bransonk@janelia.hhmi.org';
s{end+1} = '';
s{end+1} = 'This is a GUI for showing UFMF compression results. The maximum frame rate can be set through the "Preferences..." menu. Set to <= 0 for no maximum frame rate. Control the frame shown with the slider or editable text box. The Play/Stop button does what you think it does.';
msgbox(s,'About ShowUFMF','help','Replace');

% --------------------------------------------------------------------
function menu_Edit_Callback(hObject, eventdata, handles)
% hObject    handle to menu_Edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function menu_Help_Callback(hObject, eventdata, handles)
% hObject    handle to menu_Help (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes on selection change in popupmenu_Show.
function popupmenu_Show_Callback(hObject, eventdata, handles)
% hObject    handle to popupmenu_Show (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns popupmenu_Show contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popupmenu_Show
handles = update_frame(handles);
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function popupmenu_Show_CreateFcn(hObject, eventdata, handles)
% hObject    handle to popupmenu_Show (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --------------------------------------------------------------------
function menu_ShowForegroundBoxes_Callback(hObject, eventdata, handles)
% hObject    handle to menu_ShowForegroundBoxes (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

v = get(hObject,'Checked');
if strcmpi(v,'on'),
  v = 'off';
else
  v = 'on';
end
set(hObject,'Checked',v);

if isfield(handles,'hboxes'),
  set(handles.hboxes,'Visible',v);
end

% --------------------------------------------------------------------
function menu_Colormap_Callback(hObject, eventdata, handles)
% hObject    handle to menu_Colormap (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function menu_ColormapGray_Callback(hObject, eventdata, handles)
% hObject    handle to menu_ColormapGray (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if strcmpi(get(hObject,'Checked'),'off'),
  for h = handles.menu_Colormaps',
    if h ~= hObject,
      set(h,'Checked','off');
    end
  end
  set(hObject,'Checked','on');
  set(handles.hboxes,'color','r');
  handles.Colormap = 'gray';
  colormap(handles.axes_Video,handles.Colormap);
  guidata(hObject,handles);
end

% --------------------------------------------------------------------
function menu_ColormapJet_Callback(hObject, eventdata, handles)
% hObject    handle to menu_ColormapJet (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if strcmpi(get(hObject,'Checked'),'off'),
  for h = handles.menu_Colormaps',
    if h ~= hObject,
      set(h,'Checked','off');
    end
  end
  set(hObject,'Checked','on');
  set(handles.hboxes,'color','k');
  handles.Colormap = 'jet';
  colormap(handles.axes_Video,handles.Colormap);
  guidata(hObject,handles);
end

% --- Executes on key press with focus on figure1 or any of its controls.
function figure1_WindowKeyPressFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  structure with the following fields (see FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)

switch eventdata.Key,
  case {'leftarrow','downarrow'},
    handles.f = max(1,handles.f-1);
    handles = update_frame(handles);
    guidata(hObject,handles);
  case {'rightarrow','uparrow'},
    handles.f = min(handles.nframes,handles.f+1);
    handles = update_frame(handles);
    guidata(hObject,handles);
end

author: kristinbranson

repo: JAABA


text:
function gammaM = createGammaDose(doseNum1,doseNum2,strNum,doseAgreement,...
    distAgreement,thresholdPercentMax,maxDistanceFactor,doseDiffMethod,planC)
% function gammaM = createGammaDose(doseNum1,doseNum2,strNum,...
% doseAgreement,distAgreement,thresholdPercentMax,maxDistanceFactor,doseDiffMethod,planC);
%
% This function creates and adds a gamma dose to planC.
% 3D Gamma is calculated between doseNum1 and doseNum2.
% dosePercent is the dose criteria based on max(doseNum1).
% doseAgreement = max(doseNum1)*dosePercent/100;
% distAgreement is in cm.
%
% APA, 04/26/2012

if ~exist('planC','var')
    global planC
end

indexS = planC{end};

% Prepare inputs for gamma calculation
assocScan = getDoseAssociatedScan(doseNum1,planC);
if isempty(assocScan)
    assocScan = 1;
end
[newXgrid, newYgrid, newZgrid, doseArray1, doseArray2, interpMask3M] = ...
    prepareDosesForGamma(doseNum1,doseNum2,strNum,assocScan,planC);

deltaX = abs(newXgrid(2) - newXgrid(1));
deltaY = abs(newYgrid(2) - newYgrid(1));
deltaZ = abs(newZgrid(2) - newZgrid(1));

if doseDiffMethod == 1
    doseAgreement = doseAgreement*max(planC{indexS.dose}(doseNum1).doseArray(:))/100;
end

thresholdAbsolute = thresholdPercentMax*max(planC{indexS.dose}(doseNum1).doseArray(:))/100;

maxCalcDistance = distAgreement * maxDistanceFactor;

gammaM = gammaDose3d(doseArray1, doseArray2, interpMask3M, ...
    [deltaX deltaY deltaZ], doseAgreement, distAgreement, ...
    maxCalcDistance, thresholdAbsolute, doseDiffMethod);

% Assume doses within the filter threshold pass gamma
%gammaM(isnan(gammaM)) = 0;

newDoseNum = length(planC{indexS.dose}) + 1;

%Remove old caching info.
planC{indexS.dose}(newDoseNum).cachedMask = [];
planC{indexS.dose}(newDoseNum).cachedColor = [];
planC{indexS.dose}(newDoseNum).cachedTime = [];

%Set coordinates.
planC{indexS.dose}(newDoseNum).sizeOfDimension1 = length(newXgrid);
planC{indexS.dose}(newDoseNum).sizeOfDimension2 = length(newYgrid);
planC{indexS.dose}(newDoseNum).sizeOfDimension3 = length(newZgrid);
planC{indexS.dose}(newDoseNum).horizontalGridInterval = newXgrid(2)-newXgrid(1);
planC{indexS.dose}(newDoseNum).verticalGridInterval = newYgrid(2)-newYgrid(1);
planC{indexS.dose}(newDoseNum).depthGridInterval = newZgrid(2)-newZgrid(1);
planC{indexS.dose}(newDoseNum).coord1OFFirstPoint = newXgrid(1);
planC{indexS.dose}(newDoseNum).coord2OFFirstPoint = newYgrid(1);
planC{indexS.dose}(newDoseNum).coord3OfFirstPoint = newZgrid(1);
planC{indexS.dose}(newDoseNum).zValues = newZgrid;
planC{indexS.dose}(newDoseNum).doseUnits = 'Gy';
planC{indexS.dose}(newDoseNum).doseArray = gammaM;
planC{indexS.dose}(newDoseNum).doseUID = createUID('dose');
strName = planC{indexS.structures}(strNum).structureName;
planC{indexS.dose}(newDoseNum).fractionGroupID = ['Gamma_',...
    num2str(doseAgreement),'%_',num2str(distAgreement*10),'mm','_',strName];

%Switch to new dose
sliceCallBack('selectDose', num2str(newDoseNum));


text:
function yp = fpcube ( x )

%*****************************************************************************80
%
%% FPCUBE sets the derivative of the cubic function.
%
%  Discussion:
%
%    Y(X) = ( ( X + 2 ) * X + 3 ) * X + 4
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real YP, the value of the derivative of the cubic function.
%
  yp = ( 3.0E+00 * x + 4.0E+00 ) * x + 3.0E+00;

  return
end


text:
function [ node_xyz, tet_node ] = tet_mesh_order4_example_set ( node_num, ...
  tet_num )

%*****************************************************************************80
%
%% TET_MESH_ORDER4_EXAMPLE_SET sets an example linear tet mesh.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer TET_NUM, the number of tetrahedrons.
%
%    Output, real NODE_XYZ(3,NODE_NUM), the node coordinates.
%
%    Output, integer TET_NODE(4,TET_NUM), the nodes
%    forming each tet.
%
  tet_node(1:4,1:tet_num) = [ ...
         1,   2,   4,  10; ...
         2,   4,   5,  10; ...
         2,   5,  10,  11; ...
         2,   3,   5,  11; ...
         4,   5,  10,  13; ...
         3,   5,   6,  11; ...
         5,  10,  11,  13; ...
         4,   5,   7,  13; ...
         5,   6,   8,  14; ...
         5,   7,   8,  13; ...
         6,   8,   9,  14; ...
        11,  13,  14,  19; ...
        12,  14,  15,  20; ...
         3,   6,  11,  12; ...
         5,   6,  11,  14; ...
         6,   9,  14,  15; ...
         6,  11,  12,  14; ...
         6,  12,  14,  15; ...
         7,   8,  13,  16; ...
         5,   8,  13,  14; ...
        10,  11,  13,  19; ...
         8,   9,  14,  17; ...
        11,  12,  14,  20; ...
         5,  11,  13,  14; ...
         8,  13,  14,  16; ...
         9,  14,  15,  17; ...
        13,  14,  16,  22; ...
         8,  14,  16,  17; ...
        14,  15,  17,  23; ...
        14,  16,  17,  22; ...
         9,  15,  17,  18; ...
        15,  17,  18,  23; ...
        14,  17,  22,  23; ...
        13,  14,  19,  22; ...
        11,  14,  19,  20; ...
        14,  15,  20,  23; ...
        15,  20,  21,  23; ...
        21,  23,  24,  29; ...
        20,  22,  23,  28; ...
        14,  19,  20,  22; ...
        15,  18,  23,  24; ...
        12,  15,  20,  21; ...
        15,  21,  23,  24; ...
        16,  17,  22,  25; ...
        19,  20,  22,  28; ...
        17,  18,  23,  26; ...
        20,  21,  23,  29; ...
        14,  20,  22,  23; ...
        17,  22,  23,  25; ...
        18,  23,  24,  26; ...
        22,  23,  25,  31; ...
        17,  23,  25,  26; ...
        23,  24,  26,  32; ...
        23,  25,  26,  31; ...
        18,  24,  26,  27; ...
        24,  26,  27,  32; ...
        23,  26,  31,  32; ...
        22,  23,  28,  31; ...
        20,  23,  28,  29; ...
        23,  24,  29,  32; ...
        24,  29,  30,  32; ...
        30,  32,  33,  38; ...
        29,  31,  32,  37; ...
        23,  28,  29,  31; ...
        24,  27,  32,  33; ...
        21,  24,  29,  30; ...
        24,  30,  32,  33; ...
        25,  26,  31,  34; ...
        28,  29,  31,  37; ...
        26,  27,  32,  35; ...
        29,  30,  32,  38; ...
        23,  29,  31,  32; ...
        26,  31,  32,  34; ...
        27,  32,  33,  35; ...
        31,  32,  34,  40; ...
        26,  32,  34,  35; ...
        32,  33,  35,  41; ...
        32,  34,  35,  40; ...
        27,  33,  35,  36; ...
        33,  35,  36,  41; ...
        32,  35,  40,  41; ...
        31,  32,  37,  40; ...
        29,  32,  37,  38; ...
        32,  33,  38,  41; ...
        33,  38,  39,  41; ...
        39,  41,  42,  47; ...
        38,  40,  41,  46; ...
        32,  37,  38,  40; ...
        33,  36,  41,  42; ...
        30,  33,  38,  39; ...
        33,  39,  41,  42; ...
        34,  35,  40,  43; ...
        37,  38,  40,  46; ...
        35,  36,  41,  44; ...
        38,  39,  41,  47; ...
        32,  38,  40,  41; ...
        35,  40,  41,  43; ...
        36,  41,  42,  44; ...
        40,  41,  43,  49; ...
        35,  41,  43,  44; ...
        41,  42,  44,  50; ...
        41,  43,  44,  49; ...
        36,  42,  44,  45; ...
        42,  44,  45,  50; ...
        41,  44,  49,  50; ...
        40,  41,  46,  49; ...
        38,  41,  46,  47; ...
        41,  42,  47,  50; ...
        42,  47,  48,  50; ...
        48,  50,  51,  56; ...
        47,  49,  50,  55; ...
        41,  46,  47,  49; ...
        42,  45,  50,  51; ...
        39,  42,  47,  48; ...
        42,  48,  50,  51; ...
        43,  44,  49,  52; ...
        46,  47,  49,  55; ...
        44,  45,  50,  53; ...
        47,  48,  50,  56; ...
        41,  47,  49,  50; ...
        44,  49,  50,  52; ...
        45,  50,  51,  53; ...
        49,  50,  52,  58; ...
        44,  50,  52,  53; ...
        50,  51,  53,  59; ...
        50,  52,  53,  58; ...
        45,  51,  53,  54; ...
        51,  53,  54,  59; ...
        50,  53,  58,  59; ...
        49,  50,  55,  58; ...
        47,  50,  55,  56; ...
        50,  51,  56,  59; ...
        51,  56,  57,  59; ...
        50,  55,  56,  58; ...
        51,  54,  59,  60; ...
        48,  51,  56,  57; ...
        51,  57,  59,  60; ...
        52,  53,  58,  61; ...
        53,  54,  59,  62; ...
        50,  56,  58,  59; ...
        53,  58,  59,  61; ...
        54,  59,  60,  62; ...
        53,  59,  61,  62; ...
        54,  60,  62,  63 ]';

  node_xyz(1:3,1:node_num) = [ ...
  0.0,  0.0,  0.0; ...
  0.0,  0.0,  0.5; ...
  0.0,  0.0,  1.0; ...
  0.0,  0.5,  0.0; ...
  0.0,  0.5,  0.5; ...
  0.0,  0.5,  1.0; ...
  0.0,  1.0,  0.0; ...
  0.0,  1.0,  0.5; ...
  0.0,  1.0,  1.0; ...
  0.5,  0.0,  0.0; ...
  0.5,  0.0,  0.5; ...
  0.5,  0.0,  1.0; ...
  0.5,  0.5,  0.0; ...
  0.5,  0.5,  0.5; ...
  0.5,  0.5,  1.0; ...
  0.5,  1.0,  0.0; ...
  0.5,  1.0,  0.5; ...
  0.5,  1.0,  1.0; ...
  1.0,  0.0,  0.0; ...
  1.0,  0.0,  0.5; ...
  1.0,  0.0,  1.0; ...
  1.0,  0.5,  0.0; ...
  1.0,  0.5,  0.5; ...
  1.0,  0.5,  1.0; ...
  1.0,  1.0,  0.0; ...
  1.0,  1.0,  0.5; ...
  1.0,  1.0,  1.0; ...
  1.5,  0.0,  0.0; ...
  1.5,  0.0,  0.5; ...
  1.5,  0.0,  1.0; ...
  1.5,  0.5,  0.0; ...
  1.5,  0.5,  0.5; ...
  1.5,  0.5,  1.0; ...
  1.5,  1.0,  0.0; ...
  1.5,  1.0,  0.5; ...
  1.5,  1.0,  1.0; ...
  2.0,  0.0,  0.0; ...
  2.0,  0.0,  0.5; ...
  2.0,  0.0,  1.0; ...
  2.0,  0.5,  0.0; ...
  2.0,  0.5,  0.5; ...
  2.0,  0.5,  1.0; ...
  2.0,  1.0,  0.0; ...
  2.0,  1.0,  0.5; ...
  2.0,  1.0,  1.0; ...
  2.5,  0.0,  0.0; ...
  2.5,  0.0,  0.5; ...
  2.5,  0.0,  1.0; ...
  2.5,  0.5,  0.0; ...
  2.5,  0.5,  0.5; ...
  2.5,  0.5,  1.0; ...
  2.5,  1.0,  0.0; ...
  2.5,  1.0,  0.5; ...
  2.5,  1.0,  1.0; ...
  3.0,  0.0,  0.0; ...
  3.0,  0.0,  0.5; ...
  3.0,  0.0,  1.0; ...
  3.0,  0.5,  0.0; ...
  3.0,  0.5,  0.5; ...
  3.0,  0.5,  1.0; ...
  3.0,  1.0,  0.0; ...
  3.0,  1.0,  0.5; ...
  3.0,  1.0,  1.0 ]';

  return
end


text:
function handle = xyzhumanevaVisualise2(pos,fid)

% XYZHUMANEVAVISUALISE2
%
% COPYRIGHT : Carl Henrik Ek and Neil Lawrence, 2008

% MOCAP

if(nargin<2)
  fid = 2;
end

figure(fid);

% determine parametrization
if(max(max(abs(pos)))<100)
  type = 'jon';
else
  type = 'raquel';
end

joint = xyzhumaneva2joint(pos);
handle{1} = xyzhumanevaDraw(joint);axis equal;view([1 0 0]);
switch type
 case 'raquel'
  set(gca,'XLim',[-300 300],'YLim',[-300 300],'ZLim',[-800 800]);
 case 'jon'
  set(gca,'XLim',[-0.35 0.35],'YLim',[-0.35 0.35],'ZLim',[-1 0.3]);
end

figure(fid+2);
handle{2} = xyzhumanevaDraw(joint);axis equal;view([1 1 1]);
switch type
 case 'raquel'
  set(gca,'XLim',[-300 300],'YLim',[-300 300],'ZLim',[-800 800]);
 case 'jon'
  set(gca,'XLim',[-0.35 0.35],'YLim',[-0.35 0.35],'ZLim',[-1 0.3]);
end

return

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
classdef PTKAirwaysLungsAndVessels < PTKPlugin
    % PTKAirwaysLungsAndVessels. Plugin for displaying segented airways, lungs and vessels together
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %
    %
    %     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.
    %

    properties
        ButtonText = 'Airways, Lungs<BR> and vesssels'
        ToolTip = 'Shows a segmentation of the airways and lungs'
        Category = 'Lungs'

        AllowResultsToBeCached = true
        AlwaysRunPlugin = true
        PluginType = 'ReplaceOverlay'
        HidePluginInDisplay = false
        FlattenPreviewImage = true
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = true
        Visibility = 'Developer'
    end

    methods (Static)
        function results = RunPlugin(dataset, reporting)

            % Get the airway results
            [airways_results, airways_image] = dataset.GetResult('PTKAirways');

            % Get the results for the left and right lung segmentations
            [left_and_right_lungs_results, left_and_right_lungs_image] = dataset.GetResult('PTKLeftAndRightLungs');

            % Get the vessel results
            vesselness_image = dataset.GetResult('PTKVesselness');

            % Make our output image the lugn segmentation
            results = left_and_right_lungs_image.BlankCopy;

            % Make an image matrix with 1 for right lung, 2 for left lung and 3
            % for airways
            combined_image_raw = uint8(left_and_right_lungs_image.RawImage > 0);
            combined_image_raw(vesselness_image.RawImage > 0.3) = 3;
            combined_image_raw(airways_image.RawImage == 1) = 6;

            % Change the output image
            results.ChangeRawImage(combined_image_raw);

        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [y] = spm_int_D(P,M,U)
% integrates a MIMO bilinear system dx/dt = f(x,u) = A*x + B*x*u + Cu + D;
% FORMAT [y] = spm_int_D(P,M,U)
% P   - model parameters
% M   - model structure
%   M.delays - sampling delays (s); a vector with a delay for each output
%   M.states - a vector of indices if M.x(:) to be used in updating df/dx
%   M.nsteps - increase number of time steps by this number (default = 1)
%
% U   - input structure or matrix
%
% y   - response y = g(x,u,P)
%__________________________________________________________________________
% Integrates the bilinear approximation to the MIMO system described by
%
%    dx/dt = f(x,u,P) = A*x + u*B*x + C*u + D
%    y     = g(x,u,P) = L*x;
%
% at v = M.ns is the number of samples [default v = size(U.u,1)]
%
% spm_int_D will also handle static observation models by evaluating
% g(x,u,P).  It will also handle timing delays if specified in M.delays
%
%--------------------------------------------------------------------------
%
% SPM solvers or integrators
%
% spm_int_ode:  uses ode45 (or ode113) which are one and multi-step solvers
% respectively.  They can be used for any ODEs, where the Jacobian is
% unknown or difficult to compute; however, they may be slow.
%
% spm_int_J: uses an explicit Jacobian-based update scheme that preserves
% nonlinearities in the ODE: dx = (expm(dt*J) - I)*inv(J)*f.  If the
% equations of motion return J = df/dx, it will be used; otherwise it is
% evaluated numerically, using spm_diff at each time point.  This scheme is
% infallible but potentially slow, if the Jacobian is not available (calls
% spm_dx).
%
% spm_int_E: As for spm_int_J but uses the eigensystem of J(x(0)) to eschew
% matrix exponentials and inversion during the integration. It is probably
% the best compromise, if the Jacobian is not available explicitly.
%
% spm_int_B: As for spm_int_J but uses a first-order approximation to J
% based on J(x(t)) = J(x(0)) + dJdx*x(t).
%
% spm_int_L: As for spm_int_B but uses J(x(0)).
%
% spm_int_U: like spm_int_J but only evaluates J when the input changes.
% This can be useful if input changes are sparse (e.g., boxcar functions).
% It is used primarily for integrating EEG models
%
% spm_int_D: Fast integrator that uses a bilinear approximation to the
% Jacobian evaluated using spm_soreduce. This routine will also allow for
% sparse sampling of the solution and delays in observing outputs. It is
% used primarily for integrating fMRI models
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_int_D.m 5667 2013-10-02 18:26:06Z karl $

% convert U to U.u if necessary
%--------------------------------------------------------------------------
if ~isstruct(U), u.u = U; U = u; end
try, dt = U.dt; catch, U.dt = 1; end

 % number of times to sample (v) and number of microtime bins (u)
%--------------------------------------------------------------------------
u      = size(U.u,1);
try, v = M.ns;  catch, v = u;    end

% get expansion point
%--------------------------------------------------------------------------
x = [1; spm_vec(M.x)];

% add [0] states if not specified
%--------------------------------------------------------------------------
if ~isfield(M,'f')
    M.f = inline('sparse(0,1)','x','u','P','M');
    M.n = 0;
    M.x = sparse(0,0);
end

% output nonlinearity, if specified
%--------------------------------------------------------------------------
try
    g   = fcnchk(M.g,'x','u','P','M');
catch
    g   = inline('x','x','u','P','M');
    M.g = g;
end

% Bilinear approximation (1st order)
%--------------------------------------------------------------------------
[M0,M1,M2] = spm_soreduce(M,P);
n          = length(M2);                     % n states
m          = length(M1);                     % m inputs

% delays
%--------------------------------------------------------------------------
try
    D  = max(round(M.delays/U.dt),1);
catch
    D  = ones(M.l,1)*round(u/v);
end

% state-dependent effects to include during integration
%--------------------------------------------------------------------------
try
    M2 = M2(M.states);
    n  = length(M2);
end

% decrease integration time steps
%--------------------------------------------------------------------------
try
    N  = max(M.nsteps,1);
catch
    N  = 1;
end

% Evaluation times (t) and indicator array for inputs (su) and output (sy)
%==========================================================================

% get times that the input changes
%--------------------------------------------------------------------------
i     = [1 (1 + find(any(diff(U.u),2))')];
su    = sparse(1,i,1,1,u);

% get times that the response is sampled
%--------------------------------------------------------------------------
s     = ceil([0:v - 1]*u/v);
for j = 1:M.l
    i       = s + D(j);
    sy(j,:) = sparse(1,i,1:v,1,u);
end

% get (N) intervening times to evaluate
%--------------------------------------------------------------------------
i     = ceil((0:(v - 1)*N)*u/v/N) + D(1);
sx    = sparse(1,i,1,1,u);

% time in seconds
%--------------------------------------------------------------------------
t     = find(su | any(sy,1) | sx);
sy    = full(sy(:,t));
dt    = [diff(t) 0]*U.dt;

% Integrate
%--------------------------------------------------------------------------
y     = zeros(M.l,v);
J     = M0;
U.u   = full(U.u);
for i = 1:length(t)

    % input dependent changes in Jacobian
    %----------------------------------------------------------------------
    u     = U.u(t(i),:);
    J     = M0;
    for j = 1:m
        J = J + u(j)*M1{j};
    end

    % state dependent changes in Jacobian
    %----------------------------------------------------------------------
    for j = 1:n
        J = J + (x(j + 1) - M.x(j))*M2{j};
    end

    % output sampled
    %----------------------------------------------------------------------
    if any(sy(:,i))
        q      = spm_unvec(x(2:end),M.x);
        q      = spm_vec(feval(g,q,u,P,M));
        j      = find(sy(:,i));
        s      = sy(j(1),i);
        y(j,s) = q(j);
    end

    % compute updated states x = expm(J*dt)*x;
    %----------------------------------------------------------------------
    x  = spm_expm(J*dt(i),x);

    % check for convergence
    %----------------------------------------------------------------------
    if norm(x,1) > 1e6, break, end

end
y      = real(y');