text:
%kPalette 'Palette Object Menuform'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros Palette.pane file
%
% Parameters:
% MultiChoice: paletteType 'Palette Type', default: 1: 'Type of palette to display'
%    Choices are:
%   1: 'Color Bar'
%   2: 'Color Cell'
%   3: 'Color Wheel'
% InputFile: colorColorfile ' InputFile', required: 'input file'
%
% Example: kPalette(colorColorfile, {'paletteType',1;'colorColorfile',''})
%
% Khoros helpfile follows below:
% .begin tagged
%
% .item "The Palette Visual Object"
% The palette visual object provides a visual display of the colors defined by a
% colormap.
%
% .item "Palette Type"
% This attribute specifies how the palette should be displayed; the palette
% may be displayed as a Color Bar, Color Cell, or Color Wheel.
%
% .item "Input File"
% This input file selection corresponds to the XVW_COLOR_COLORFILE attribute
% used by the palette object.  It is used to specify which input file should
% be used as the data object for which the palette will display colors.
%
% .end tagged

function varargout = kPalette(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kPalette(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'paletteType', 1;'colorColorfile', '__input'};
maxval={0,0};
minval={0,0};
istoggle=[0,0];
was_set=istoggle * 0;
paramtype={'MultiChoice','InputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(0);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
callKhoros([w 'Palette" '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% ----------------------------------------------------------------------------
% function hfssAssignMaster(fid, Name, ObjName, iUStart, iUEnd, Units,
%                               [ReverseV])
%
% Description :
% -------------
% Creates the necessary VB Script to assign a Master Boundary to a given Object.
%
% Parameters :
% ------------
% fid     - file identifier of the HFSS script file.
% Name    - name of the master boundary (appears under 'Boundaries' in HFSS).
% ObjName - name of the (sheet-like) object to which the master boundary is to
%           be assigned.
% iUStart - (vector) starting point of the U vector. Specify as
%           [x, y, z].
% iUEnd   - (vector) ending point of the U vector. Specify as
%           [x, y, z].
% Units   - specify as 'meter', 'in', 'cm' (defined in HFSS).
% [ReverseV] - (boolean, optional) reverses vector V (defaults to
%                true if the U vector points +y, and to false elsewise)
% Example :
% ---------
% fid = fopen('myantenna.vbs', 'wt');
% ...
% hfssAssignMaster(fid, 'Master', 'Sheet', [-width/2, 0, 0], ...
%	                 [width/2, 0, 0], 'meter', false);
%
% ----------------------------------------------------------------------------

% ----------------------------------------------------------------------------
% CHANGELOG
%
% 20-May-2013: *Initial release.
% ----------------------------------------------------------------------------

% ----------------------------------------------------------------------------
% Written by Pablo Alcon Garcia
% pabloalcongarcia@gmail.com / palcon@tsc.uniovi.es
% 20 May 2013
% ----------------------------------------------------------------------------
function hfssAssignMaster(fid, Name, ObjName, iUStart, iUEnd, Units, ReverseV)

% arguments processor.
if (nargin < 6)
	error('Insufficient # of arguments !');
elseif (nargin < 7)
    if iUEnd(2)~=iUStart(2)
        ReverseV = true;
    else
        ReverseV = false;
    end
end

if ReverseV
    ReverseV = 'true';
else
    ReverseV = 'false';
end

% Preamble.
fprintf(fid, '\n');
fprintf(fid, 'Set oModule = oDesign.GetModule("BoundarySetup")\n');

% Parameters
fprintf(fid, 'oModule.AssignMaster _\n');
fprintf(fid, 'Array("NAME:%s", _\n', Name);
fprintf(fid, '\t"Objects:=", Array("%s"), _\n', ObjName);
fprintf(fid, '\tArray("NAME:CoordSysVector", "Origin:=", _\n');
fprintf(fid, '\t\tArray("%f%s", "%f%s", "%f%s"), _\n', ...
        iUStart(1), Units, iUStart(2), Units, iUStart(3), Units);
fprintf(fid, '\t\t"UPos:=", Array("%f%s", "%f%s", "%f%s") _\n', ...
        iUEnd(1), Units, iUEnd(2), Units, iUEnd(3), Units);
fprintf(fid, '\t\t), _\n');
fprintf(fid, '\t"ReverseV:=", %s _\n',ReverseV);
fprintf(fid, '\t)\n');

author: yuip

repo: hfss-api

save_path: github-repos/MATLAB/yuip-hfss-api


text:
function json=moxunit_fieldtrip_util_struct2json(s)
% convert struct to json string representation
%
% json=moxunit_fieldtrip_util_struct2json(s)
%
% Inputs:
%   s                   struct for which the values must all be strings,
%                       or scalar logical or double arrays
%
% Output:
%   json                JSON string representation of the input
%

    if ~isstruct(s)
        error('input must be a struct');
    end

    keys=fieldnames(s);
    n=numel(keys);
    args=cell(2,n);
    for k=1:n
        key=keys{k};
        raw_value=s.(key);

        if ischar(raw_value)
            value=sprintf('"%s"',escape(raw_value));

        elseif isnumeric(raw_value) && numel(raw_value)==1
            value=sprintf('%d',raw_value);

        elseif islogical(raw_value) && numel(raw_value)==1
            if raw_value
                value='true';
            else
                value='false';
            end

        else
            error('Unsupported input for key %s: <class %s>',...
                            key, class(raw_value));
        end

        args{1,k}=sprintf('"%s":%s',key,value);

        if k<n
            sep=',';
        else
            sep='';
        end
        args{2,k}=sep;
    end

    json=sprintf('{%s}',sprintf('%s',args{:}));

function s=escape(s)
    s=regexprep(s,'"','\\"');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [R_ned_2_frame, ned_x, ned_y, ned_z] = computeNedFrameInFrame(rVect)
    %Source: https://en.wikipedia.org/wiki/North_east_down

    rNorm = norm(rVect);
    lambda = AngleZero2Pi(atan2(rVect(2),rVect(1)));
    phi = pi/2 - acos(rVect(3)/rNorm);

    RFrame2Ned = [-sin(phi)*cos(lambda), -sin(lambda), -cos(phi)*cos(lambda);
                 -sin(phi)*sin(lambda),  cos(lambda), -cos(phi)*sin(lambda);
                  cos(phi),              0,           -sin(phi)]';

    R_ned_2_frame = RFrame2Ned';
    ned_x = R_ned_2_frame(:,1);
    ned_y = R_ned_2_frame(:,2);
    ned_z = R_ned_2_frame(:,3);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function Z = conj(Y)
%CONJ (overloaded)

Z = Y;
Z.basis = conj(Z.basis);
% Reset info about conic terms
Z.conicinfo = [0 0];

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function demoEquivalentEllipsoid
%DEMOEQUIVALENTELLIPSOID Demo program for the use of ellipsoids
%
%   Example
%     demoEquivalentEllipsoid
%
%   See also
%     demoRevolutionSurface, demoDrawTubularMesh

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2011-06-21,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

%% Data generation

% Generate gaussian 3D data
nPoints = 1000;
points = randn(nPoints, 3);

% point clouds parameters
center = [20 30 40];
sizes  = [70 40 10];
orient = [50 30 30];

% transform points to make a gaussian cloud
transfo = composeTransforms3d(...
    createScaling3d(sizes), ...
    eulerAnglesToRotation3d(orient), ...
    createTranslation3d(center));
points = transformPoint3d(points, transfo);

% display data
figure;
drawPoint3d(points, '.');
hold on;
axis equal;
view([80 -10]);

%% Equivalent ellipsoid computation and display

% Fit a 3D equivalent ellipsoid to data
elli = equivalentEllipsoid(points);

% draw the ellipsoid with transparency
drawEllipsoid(elli, 'FaceColor', 'g', 'FaceAlpha', .5);

%% Add ellipses and main axes

drawEllipsoid(elli, 'FaceColor', 'g', 'FaceAlpha', .5, ...
          'drawEllipses', true, 'EllipseColor', 'b', 'EllipseWidth', 2, ...
          'drawAxes', true);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function[varargout]=cellsum(varargin)
%CELLSUM  Sum of each element a cell array, possibly weighted.
%
%   S=CELLSUM(X) where X is a cell array of N arrays, is equivalent to
%
%      S(1,1)=SUM(X{1}(:)),  S(2,1)=SUM(X{2}(:)), ...,  S(N,1)=SUM(X{N}(:))
%
%   thus returning an N x 1 array containing the sum over all values of
%   each element in the cell array.
%
%   In taking the sum, non-finite values are ignored, as in VSUM.
%   S is a column vector of the same length as X.
%
%   [S1,S2,...,SP]=CELLSUM(X1,X2,...,XP) also works for P different
%   input arguments.
%
%   CELLSUM(X1,X2,...XP);  with no output arguments overwrites the
%   original input variables.
%   __________________________________________________________________
%
%   Weighted sums
%
%   CELLSUM(X,'weight',W) or CELLSUM(X1,X2,...,XP,'weight',W) where W is
%   a cell array of the same size as the other input variables, computes
%   the weighted sum, using the weighting factor ABS(W).^2.
%   __________________________________________________________________
%
%   Parallelization
%
%   CELLSUM(...,'parallel') parallelizes the computation using a PARFOR
%   loop.  This requires that Matlab's Parallel Computing Toolbox be
%   installed, and is useful for very large datasets.
%   __________________________________________________________________
%
%   See also CELLSTD, CELLMED, JCELL.
%
%   Usage: s=cellsum(x);
%          [s1,s2,s3]=cellsum(x1,x2,x3);
%          cellsum(x1,x2,x3,'weight',w);
%          cellsum(x1,x2,x3);
%          cellsum(x1,x2,x3,'parallel');
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2014--2019 J.M. Lilly --- type 'help jlab_license' for details

if ~iscell(varargin{1})
    error('X must be a cell array.')
end

weight=[];
cores='serial';

if ischar(varargin{end})
    if strcmpi(varargin{end}(1:3),'par')||strcmpi(varargin{end}(1:3),'ser')
        cores=varargin{end};
    end
    varargin=varargin(1:end-1);
end

if length(varargin)>1
    if ischar(varargin{end-1})
        if strcmpi(varargin{end-1}(1:3),'wei')
            weight=varargin{end};
        end
        varargin=varargin(1:end-2);
    end
end

if strcmpi(cores(1:3),'par')
    if exist('parpool')~=2
        disp('Sorry, parallel algorithm requires the Parallel Computing Toolbox.')
        disp('Defaulting to the standard algorithm.')
        cores='serial';
    end
end

for i=1:length(varargin)
    if isempty(weight)
        varargout{i}=cellsum_one_unweighted(varargin{i},cores);
    else
        varargout{i}=cellsum_one_weighted(varargin{i},weight,cores);
    end
end

eval(to_overwrite(length(varargin)))

function[y]=cellsum_one_unweighted(x,cores)

y=nan*zeros(length(x),1);
if strcmpi(cores(1:3),'par')
    parfor i=1:length(x)
        y(i)=sum(x{i}(:),1);
        if ~isfinite(y(i))
            y(i)=vsum(x{i}(:),1);
        end
    end
else
    for i=1:length(x)
        y(i)=sum(x{i}(:),1);
        if ~isfinite(y(i))
            y(i)=vsum(x{i}(:),1);
        end
    end
end

function[y]=cellsum_one_weighted(x,weight,cores)

y=nan*zeros(length(x),1);
if strcmpi(cores(1:3),'par')
    parfor i=1:length(x)
        w=squared(weight{i}(:));
        y(i)=sum(w.*x{i}(:),1);
        if ~isfinite(y(i))
            y(i)=vsum(w.*x{i}(:),1);
        end
    end
else
    for i=1:length(x)
        w=squared(weight{i}(:));
        y(i)=sum(w.*x{i}(:),1);
        if ~isfinite(y(i))
            y(i)=vsum(w.*x{i}(:),1);
        end
    end
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [ sepScene, wallPanoNormal] = compSurfaceLabel( rotImg )%compSurfaceLabel( data, folderName )
%COMPSURFACELABEL Compute geometric context surface label for panorama
%   sepScene: surface label for each seprate view, 6 channels
%   wallPanoNormal: surface label for panorama, 7 channels

% [rotImg, data, R] = loadBuffer( data, folderName );
%% emprical good separate of panorama
cutSize = 640;
% fov = pi*2/3;
fov = [pi/2*ones(1,8) pi/2*ones(1,8)];
xh = [-pi:(pi/4):(3/4*pi) -pi:(pi/4):(3/4*pi)];
yh = [zeros(1, 8) -1/6*pi*ones(1,8)];

% fov = [pi/2*ones(1,36) pi/2*ones(1,36)];
% xh = [-pi:(pi/18):(17/18*pi) -pi:(pi/18):(17/18*pi)];
% yh = [zeros(1,36) -1/6*pi*ones(1,36)];

% fov = [pi/2*ones(1,8)];
% xh = [-pi:(pi/4):(3/4*pi)];
% yh = [zeros(1,8)];

% fov = [pi/2*ones(1,36)];
% xh = [-pi:(pi/18):(17/18*pi)];
% yh = [zeros(1,36)];
[sepScene] = separatePano( rotImg, fov, xh, yh, cutSize);

%%
% imdir = './roomModel/SpatialLayout_shrink/Images_resized/';
% workspcdir = './roomModel/SpatialLayout_shrink/tempworkspace/';
% imsegdir = './roomModel/SpatialLayout_shrink/Imsegs/';
% moveFolder = './roomModel/SpatialLayout_shrink/';
sepScene(1).cimages = [];
sepScene(1).maxCimage = [];
numScene = length(sepScene);
for sid = 1:numScene
    img = sepScene(sid).img;
    img = imresize(img, [640 640]);
    try
        [ cimages ] = getSurfaceLabelInitial( img );
    catch
        cimages{1} = 1/7*ones(cutSize,cutSize,7);
    end
%     [cimages, Polyg] = getSurfaceLabelLayout( sepScene(sid), data  )
    sepScene(sid).cimages = cimages;
    [~, indd] = max(cimages{1}, [], 3);
    sepScene(sid).maxCimage = indd;
end

%% combine
% 1 middle wall=yellow; 2 left wall=red; 3 right wall=cyan; 4 floor=green; 5 ceiling=blue; 6 objects=pink
width = 2048;   height = 1024;
wallPano = zeros(height, width, 7);
wallPanoWeight = zeros(height, width);
for vid = 1:length(sepScene)
    vx = sepScene(vid).vx;

    [curPano, curPanoValid] = ...
    im2Sphere( sepScene(vid).cimages{1}, sepScene(vid).fov, width, height, ...
                sepScene(vid).vx, sepScene(vid).vy);
    curPano(~curPanoValid) = 0;
    wallPanoWeight(curPanoValid(:,:,1)) = wallPanoWeight(curPanoValid(:,:,1)) + 1;

    normMiddle = vx+pi;
    normMLR = normValue([normMiddle normMiddle-pi/2 normMiddle+pi/2], -pi, pi);
    normGlobal = [-pi -pi/2 0 pi/2];
    for i = 1:3
        angle = abs(normValue(normGlobal - normMLR(i), -pi, pi));
        angle(angle>pi/2) = pi/2; % cos^2 weight
%         [~,I] = sort(angle, 'ascend'); % nearest neighbor
%         angle(I(1)) = 0;
%         angle(I(2:4)) = pi/2;
        for j = find(angle<pi/2)
            wallPano(:,:,j) = wallPano(:,:,j) + curPano(:,:,i)*(cos(angle(j)).^2);
        end
    end

    wallPano(:,:,5) = wallPano(:,:,5) + curPano(:,:,4); % floor
    wallPano(:,:,6) = wallPano(:,:,6) + curPano(:,:,5); % ceiling
    wallPano(:,:,7) = wallPano(:,:,7) + curPano(:,:,6); % object
end
wallPanoNormal = wallPano./repmat(wallPanoWeight+0.000001,[1 1 7]);

%% more informative combine
% % Left wall: view+(90-FOV/2) ~ view+(90-FOV/2)+180
% % Right wall: view-(90-FOV/2)-180 ~ view-(90-FOV/2)
% % Middle wall: Left wall + Right wall
% width = 1024;   height = 512;
% wallPano = zeros(512, 1024, 4);
% for vid = 1:length(sepScene)
%     vx = sepScene(vid).vx;
%     fov = sepScene(vid).fov;
%     [curPano, curPanoValid] = ...
%         im2Sphere( sepScene(vid).cimages{1}, sepScene(vid).fov, width, height, ...
%                     sepScene(vid).vx, sepScene(vid).vy);
%     curPano(~curPanoValid) = 0;
%     %Left
%     lowBound = vx+(pi/2-fov/2);
%     uppBound = vx+(pi/2-fov/2)+pi;
%     range = normRange(lowBound, uppBound, -pi, pi);
%     Linside = insideRange([-pi -pi/2 0 pi/2], range);
% %     Lweight = [1 1 1 1];
%     %Right
%     lowBound = vx-(pi/2-fov/2)-pi;
%     uppBound = vx-(pi/2-fov/2);
%     range = normRange(lowBound, uppBound, -pi, pi);
%     Rinside = insideRange([-pi -pi/2 0 pi/2], range);
% %     Rweight = [1 1 1 1];
%     %Middle
%     lowBound = vx+(pi/2-fov/2);
%     uppBound = vx-(pi/2-fov/2);
%     range = normRange(lowBound, uppBound, -pi, pi);
%     Minside = insideRange([-pi -pi/2 0 pi/2], range);
% %     Mweight = [1 1 1 1];
%     %Vote Left
%     score = curPano(:,:,2);
%     voteMatrix = repmat(reshape(Linside, [1 1 4]), [512 1024 1]);
%     voteMatrix = voteMatrix .* repmat(score, [1 1 4]);
%     wallPano = wallPano + voteMatrix;
%
%     %Vote Right
%     score = curPano(:,:,3);
%     voteMatrix = repmat(reshape(Linside, [1 1 4]), [512 1024 1]);
%     voteMatrix = voteMatrix .* repmat(score, [1 1 4]);
%     wallPano = wallPano + voteMatrix;
%
%     %Vote Middle
%     score = curPano(:,:,1);
%     voteMatrix = repmat(reshape(Linside, [1 1 4]), [512 1024 1]);
%     voteMatrix = voteMatrix .* repmat(score, [1 1 4]);
%     wallPano = wallPano + voteMatrix;
%
% end

%% heuristic combine
% MASK = zeros(640,640);
% MASK(213:427,1:320) = 1;
% MASK(213:427,321:640) = 2;
% All_MLR2ABCD = zeros(8,3);
% for vid = [2 4 6 8]
%     histLabelL = hist(sepScene(vid).maxCimage(MASK<1.5 & MASK>0.5), 1:3);
%     histLabelR = hist(sepScene(vid).maxCimage(MASK>1.5 & MASK<2.5), 1:3);
%     histLabelL = histLabelL./sum(histLabelL) + 0.0001;
%     histLabelR = histLabelR./sum(histLabelR) + 0.0001;
%     entropyL = -sum(histLabelL.*log2(histLabelL));
%     entropyR = -sum(histLabelR.*log2(histLabelR));
%     MLR2ABCD = [0 0 0]; %[M L R], 1234 for ABCD, 0 for Null
%     histLabel = hist(sepScene(vid).maxCimage(MASK<2.5 & MASK>0.5), 1:3);
%     [~, IW] = sort(histLabel, 'descend');
%     if entropyL<entropyR
%         [~, I] = max(histLabelL);
%         MLR2ABCD(I) = vid/2;
%         if ~ismember(I, IW(1:2))
%             fprintf('Probably some problems\n');
%         end
%         if I==1
%             MLR2ABCD(3) = rem(vid/2,4)+1;
%         else
%             MLR2ABCD(setdiff(IW(1:2),I)) = rem(vid/2,4)+1;
%         end
%     else
%         [~, I] = max(histLabelR);
%         MLR2ABCD(I) = rem(vid/2,4)+1;
%         if ~ismember(I, IW(1:2))
%             fprintf('Probably some problems\n');
%         end
%         if I==1
%             MLR2ABCD(2) = vid/2;
%         else
%             MLR2ABCD(setdiff(IW(1:2),I)) = vid/2;
%         end
%     end
%     All_MLR2ABCD(vid,:) = MLR2ABCD;
% end
% width = 1024;   height = 512;
% for vid = [1 3 5 7]
%     preID = vid - 1; preID(preID==0) = 8;
%     latID = vid + 1;
%
%     curMaxCimage = sepScene(vid).maxCimage;
%     curMaxCimage(MASK<0.5) = 0;
%     [curPano, curPanoValid] = ...
%         im2Sphere( curMaxCimage, sepScene(vid).fov, width, height, ...
%                     sepScene(vid).vx, sepScene(vid).vy);
%     curPano(~curPanoValid) = 0;
%
%     preMaxCimage = sepScene(preID).maxCimage;
%     preMaxCimage(MASK<0.5) = 0;
%     [prePano, prePanoValid] = ...
%         im2Sphere( preMaxCimage, sepScene(preID).fov, width, height, ...
%                     sepScene(preID).vx, sepScene(preID).vy);
%     prePano(~prePanoValid) = 0;
%
%     latMaxCimage = sepScene(latID).maxCimage;
%     latMaxCimage(MASK<0.5) = 0;
%     [latPano, latPanoValid] = ...
%         im2Sphere( latMaxCimage, sepScene(latID).fov, width, height, ...
%                     sepScene(latID).vx, sepScene(latID).vy);
%     latPano(~latPanoValid) = 0;
%
%     preConsist = [0 0 0];
%     preNumber = [0 0 0];
%     for i = 1:3
%         mask = curPano>i-0.5 & curPano<i+0.5 & prePano>0.5 & prePano<3.5;
%         h = hist(prePano(mask),1:3);
%         [B,I] = max(h);
%         J = All_MLR2ABCD(preID,I);
%         preConsist(i) = J;
%         preNumber(i) = B;
%
%         mask = curPano>i-0.5 & curPano<i+0.5;
%         if sum(mask(:))<250 || B/sum(mask(:))<0.05
%             preConsist(i) = 0;
%         end
%     end
%     latConsist = [0 0 0];
%     latNumber = [0 0 0];
%     for i = 1:3
%         mask = curPano>i-0.5 & curPano<i+0.5 & latPano>0.5 & latPano<3.5;
%         h = hist(latPano(mask),1:3);
%         [B,I] = max(h);
%         J = All_MLR2ABCD(latID,I);
%         latConsist(i) = J;
%         latNumber(i) = B;
%
%         mask = curPano>i-0.5 & curPano<i+0.5;
%         if sum(mask(:))<250 || B/sum(mask(:))<0.05
%             latConsist(i) = 0;
%         end
%     end
%     for i = 1:3
%         if preNumber(i)>latNumber(i)
%             All_MLR2ABCD(vid, i) = preConsist(i);
%         else
%             All_MLR2ABCD(vid, i) = latConsist(i);
%         end
%     end
% end
% wallPano = zeros(512, 1024, 4);
% for vid = [2 4 6 8]
%     [curPano, curPanoValid] = ...
%     im2Sphere( sepScene(vid).cimages{1}, sepScene(vid).fov, width, height, ...
%                 sepScene(vid).vx, sepScene(vid).vy);
%     curPano(~curPanoValid) = 0;
%     voteMat = zeros(512, 1024, 4);
%     for i = find(All_MLR2ABCD(vid,:)~=0)
%         voteMat(:,:,All_MLR2ABCD(vid,i)) = curPano(:,:,i);
%     end
%     wallPano = wallPano + voteMat;
% end

end

function [nv] = normValue( v, minBound, maxBound)
nv = rem((v-minBound)+2*(maxBound-minBound), maxBound-minBound) + minBound;
end

function [range] = normRange( low, upp, minBound, maxBound)
normLow = rem((low-minBound)+2*(maxBound-minBound), maxBound-minBound) + minBound;
normUpp = rem((upp-minBound)+2*(maxBound-minBound), maxBound-minBound) + minBound;
if normLow<=normUpp
    range = [normLow normUpp];
else
    range = [normLow maxBound; minBound normUpp];
end
end

function [inside] = insideRange(values, range)
inside = false(size(values));
for i = 1:size(range,1)
    inside = inside | (values>=range(i,1) & values<=range(i,2));
end
end

author: zouchuhang

repo: LayoutNet


text:
function out = spm_deformations(job)
% Various deformation field utilities
% FORMAT out = spm_deformations(job)
% job - a job created via spm_cfg_deformations.m
% out - a struct with fields
%       .def    - file name of created deformation field
%       .warped - file names of warped images
%
% See spm_cfg_deformations.m for more information.
%__________________________________________________________________________
% Copyright (C) 2005-2015 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_deformations.m 6577 2015-10-15 15:22:11Z volkmar $

[Def,mat] = get_comp(job.comp);
out = struct('def',{{}},'warped',{{}},'surf',{{}},'jac',{{}});
for i=1:numel(job.out)
    fn = fieldnames(job.out{i});
    fn = fn{1};
    switch fn
    case 'savedef'
        out.def    = [out.def;    save_def(Def,mat,job.out{i}.(fn))];
    case 'pull'
        out.warped = [out.warped; pull_def(Def,mat,job.out{i}.(fn))];
    case 'push'
        out.warped = [out.warped; push_def(Def,mat,job.out{i}.(fn))];
    case 'surf'
        out.surf   = [out.surf;   surf_def(Def,mat,job.out{i}.(fn))];
    case 'savejac'
        out.jac    = [out.jac;    jac_def(Def,mat,job.out{i}.(fn))];
    otherwise
        error('Unknown option');
    end
end

%==========================================================================
% function [Def,mat] = get_comp(job)
%==========================================================================
function [Def,mat] = get_comp(job)
% Return the composition of a number of deformation fields.

if isempty(job)
    error('Empty list of jobs in composition');
end
[Def,mat] = get_job(job{1});
for i=2:numel(job)
    Def1         = Def;
    mat1         = mat;
    [Def,mat]    = get_job(job{i});
    M            = inv(mat1);
    tmp          = zeros(size(Def),'single');
    tmp(:,:,:,1) = M(1,1)*Def(:,:,:,1)+M(1,2)*Def(:,:,:,2)+M(1,3)*Def(:,:,:,3)+M(1,4);
    tmp(:,:,:,2) = M(2,1)*Def(:,:,:,1)+M(2,2)*Def(:,:,:,2)+M(2,3)*Def(:,:,:,3)+M(2,4);
    tmp(:,:,:,3) = M(3,1)*Def(:,:,:,1)+M(3,2)*Def(:,:,:,2)+M(3,3)*Def(:,:,:,3)+M(3,4);
    Def(:,:,:,1) = single(spm_diffeo('bsplins',Def1(:,:,:,1),tmp,[1,1,1,0,0,0]));
    Def(:,:,:,2) = single(spm_diffeo('bsplins',Def1(:,:,:,2),tmp,[1,1,1,0,0,0]));
    Def(:,:,:,3) = single(spm_diffeo('bsplins',Def1(:,:,:,3),tmp,[1,1,1,0,0,0]));
    clear tmp
end

%==========================================================================
% function [Def,mat] = get_job(job)
%==========================================================================
function [Def,mat] = get_job(job)
% Determine what is required, and pass the relevant bit of the
% job out to the appropriate function.

fn = fieldnames(job);
fn = fn{1};
switch fn
    case {'comp'}
        [Def,mat] = get_comp(job.(fn));
    case {'def'}
        [Def,mat] = get_def(job.(fn));
    case {'dartel'}
        [Def,mat] = get_dartel(job.(fn));
    case {'sn2def'}
        [Def,mat] = get_sn2def(job.(fn));
    case {'inv'}
        [Def,mat] = get_inv(job.(fn));
    case {'id'}
        [Def,mat] = get_id(job.(fn));
    case {'idbbvox'}
        [Def,mat] = get_idbbvox(job.(fn));
    otherwise
        error('Unrecognised job type');
end

%==========================================================================
% function [Def,mat] = get_sn2def(job)
%==========================================================================
function [Def,mat] = get_sn2def(job)
% Convert a SPM _sn.mat file into a deformation field, and return it.

vox = job.vox;
bb  = job.bb;
sn  = load(job.matname{1});

if any(isfinite(bb(:))) || any(isfinite(vox))
    [bb0, vox0] = spm_get_bbox(sn.VG(1));
    if any(~isfinite(vox)), vox = vox0; end
    if any(~isfinite(bb)),  bb  = bb0;  end
    bb  = sort(bb);
    vox = abs(vox);

    % Adjust bounding box slightly - so it rounds to closest voxel.
    bb(:,1) = round(bb(:,1)/vox(1))*vox(1);
    bb(:,2) = round(bb(:,2)/vox(2))*vox(2);
    bb(:,3) = round(bb(:,3)/vox(3))*vox(3);

    M   = sn.VG(1).mat;
    vxg = sqrt(sum(M(1:3,1:3).^2));
    ogn = M\[0 0 0 1]';
    ogn = ogn(1:3)';

    % Convert range into range of voxels within template image
    x   = (bb(1,1):vox(1):bb(2,1))/vxg(1) + ogn(1);
    y   = (bb(1,2):vox(2):bb(2,2))/vxg(2) + ogn(2);
    z   = (bb(1,3):vox(3):bb(2,3))/vxg(3) + ogn(3);

    og  = -vxg.*ogn;
    of  = -vox.*(round(-bb(1,:)./vox)+1);
    M1  = [vxg(1) 0 0 og(1) ; 0 vxg(2) 0 og(2) ; 0 0 vxg(3) og(3) ; 0 0 0 1];
    M2  = [vox(1) 0 0 of(1) ; 0 vox(2) 0 of(2) ; 0 0 vox(3) of(3) ; 0 0 0 1];
    mat = sn.VG(1).mat*(M1\M2);
    % dim = [length(x) length(y) length(z)];
else
    dim = sn.VG(1).dim;
    x   = 1:dim(1);
    y   = 1:dim(2);
    z   = 1:dim(3);
    mat = sn.VG(1).mat;
end

[X,Y] = ndgrid(x,y);

st = size(sn.Tr);

if (prod(st) == 0)
    affine_only = true;
    basX = 0;
    basY = 0;
    basZ = 0;
else
    affine_only = false;
    basX = spm_dctmtx(sn.VG(1).dim(1),st(1),x-1);
    basY = spm_dctmtx(sn.VG(1).dim(2),st(2),y-1);
    basZ = spm_dctmtx(sn.VG(1).dim(3),st(3),z-1);
end

Def = zeros([numel(x),numel(y),numel(z),3],'single');

for j=1:length(z)
    if (~affine_only)
        tx = reshape( reshape(sn.Tr(:,:,:,1),st(1)*st(2),st(3)) *basZ(j,:)', st(1), st(2) );
        ty = reshape( reshape(sn.Tr(:,:,:,2),st(1)*st(2),st(3)) *basZ(j,:)', st(1), st(2) );
        tz = reshape( reshape(sn.Tr(:,:,:,3),st(1)*st(2),st(3)) *basZ(j,:)', st(1), st(2) );

        X1 = X    + basX*tx*basY';
        Y1 = Y    + basX*ty*basY';
        Z1 = z(j) + basX*tz*basY';
    end

    Mult = sn.VF.mat*sn.Affine;
    if (~affine_only)
        X2= Mult(1,1)*X1 + Mult(1,2)*Y1 + Mult(1,3)*Z1 + Mult(1,4);
        Y2= Mult(2,1)*X1 + Mult(2,2)*Y1 + Mult(2,3)*Z1 + Mult(2,4);
        Z2= Mult(3,1)*X1 + Mult(3,2)*Y1 + Mult(3,3)*Z1 + Mult(3,4);
    else
        X2= Mult(1,1)*X + Mult(1,2)*Y + (Mult(1,3)*z(j) + Mult(1,4));
        Y2= Mult(2,1)*X + Mult(2,2)*Y + (Mult(2,3)*z(j) + Mult(2,4));
        Z2= Mult(3,1)*X + Mult(3,2)*Y + (Mult(3,3)*z(j) + Mult(3,4));
    end

    Def(:,:,j,1) = single(X2);
    Def(:,:,j,2) = single(Y2);
    Def(:,:,j,3) = single(Z2);
end

%==========================================================================
% function [Def,mat] = get_def(job)
%==========================================================================
function [Def,mat] = get_def(job)
% Load a deformation field saved as an image
Nii = nifti(job{1});
Def = single(Nii.dat(:,:,:,1,:));
d   = size(Def);
if d(4)~=1 || d(5)~=3, error('Deformation field is wrong!'); end
Def = reshape(Def,[d(1:3) d(5)]);
mat = Nii.mat;

%==========================================================================
% function [Def,mat] = get_dartel(job)
%==========================================================================
function [Def,mat] = get_dartel(job)

Nii    = nifti(job.flowfield{1});

if ~isempty(job.template{1})
   %Nt        = nifti(job.template{1});
    [pth,nam] = fileparts(job.template{1});
    if exist(fullfile(pth,[nam '_2mni.mat']),'file')
        load(fullfile(pth,[nam '_2mni.mat']),'mni');
    else
        % Affine registration of Dartel Template with MNI space.
        %------------------------------------------------------------------
        fprintf('** Affine registering "%s" with MNI space **\n', nam);
        tpm        = fullfile(spm('Dir'),'tpm','TPM.nii');
        Mmni       = spm_get_space(tpm);
        Nt         = nifti(job.template{1});
        mni.affine = Mmni/spm_klaff(Nt,tpm);
        mni.code   = 'MNI152';
        save(fullfile(pth,[nam '_2mni.mat']),'mni', spm_get_defaults('mat.format'));
    end
    Mat    = mni.affine;
   %do_aff = true;
else
    Mat    = Nii.mat;
   %do_aff = false;
end

% Integrate a Dartel flow field
y0      = spm_dartel_integrate(Nii.dat,job.times,job.K);
if all(job.times == [0 1]),
    mat = Nii.mat0;
    Def = affine(y0,single(Mat));
else
    mat = Mat;
    Def = affine(y0,single(Nii.mat0));
end

%==========================================================================
% function [Def,mat] = get_id(job)
%==========================================================================
function [Def,mat] = get_id(job)
% Get an identity transform based on an image volume
[mat, dim] = spm_get_matdim(job.space{1});
Def = identity(dim, mat);

%==========================================================================
% function [Def,mat] = get_idbbvox(job)
%==========================================================================
function [Def,mat] = get_idbbvox(job)
% Get an identity transform based on bounding box and voxel size.
% This will produce a transversal image.
[mat, dim] = spm_get_matdim('', job.vox, job.bb);
Def = identity(dim, mat);

%==========================================================================
% function [Def,mat] = get_inv(job)
%==========================================================================
function [Def,mat] = get_inv(job)
% Invert a deformation field (derived from a composition of deformations)

NT          = nifti(job.space{:});
[Def0,mat0] = get_comp(job.comp);
M0          = mat0;
M1          = inv(NT.mat);
M0(4,:)     = [0 0 0 1];
M1(4,:)     = [0 0 0 1];
Def         = spm_diffeo('invdef',Def0,NT.dat.dim(1:3),M1,M0);
mat         = NT.mat;
Def         = spm_extrapolate_def(Def,mat);

%==========================================================================
% function fname = save_def(Def,mat,job)
%==========================================================================
function fname = save_def(Def,mat,job)
% Save a deformation field as an image

ofname = job.ofname;
if isempty(ofname), fname = {}; return; end;

[pth,nam] = fileparts(ofname);
if isfield(job.savedir,'savepwd')
    wd = pwd;
elseif isfield(job.savedir,'saveusr')
    wd = job.savedir.saveusr{1};
else
    wd = pwd;
end

fname = {fullfile(wd,['y_' nam '.nii'])};
dim   = [size(Def,1) size(Def,2) size(Def,3) 1 3];
dtype = 'FLOAT32-LE';
off   = 0;
scale = 1;
inter = 0;
dat   = file_array(fname{1},dim,dtype,off,scale,inter);

N      = nifti;
N.dat  = dat;
N.mat  = mat;
N.mat0 = mat;
N.mat_intent  = 'Aligned';
N.mat0_intent = 'Aligned';
N.intent.code = 'VECTOR';
N.intent.name = 'Mapping';
N.descrip     = 'Deformation field';
create(N);
N.dat(:,:,:,1,1) = Def(:,:,:,1);
N.dat(:,:,:,1,2) = Def(:,:,:,2);
N.dat(:,:,:,1,3) = Def(:,:,:,3);

%==========================================================================
% function fname = jac_def(Def,mat,job)
%==========================================================================
function fname = jac_def(Def,mat,job)
% Save Jacobian determinants of deformation field

ofname = job.ofname;
if isempty(ofname), fname = {}; return; end;

[pth,nam] = fileparts(ofname);
if isfield(job.savedir,'savepwd')
    wd = pwd;
elseif isfield(job.savedir,'saveusr')
    wd = job.savedir.saveusr{1};
else
    wd = pwd;
end

Dets  = spm_diffeo('def2det',Def)/det(mat(1:3,1:3));
Dets(:,:,[1 end]) = NaN;
Dets(:,[1 end],:) = NaN;
Dets([1 end],:,:) = NaN;

fname = {fullfile(wd,['j_' nam '.nii'])};
dim   = [size(Def,1) size(Def,2) size(Def,3) 1 1];
dtype = 'FLOAT32-LE';
off   = 0;
scale = 1;
inter = 0;
dat   = file_array(fname{1},dim,dtype,off,scale,inter);

N      = nifti;
N.dat  = dat;
N.mat  = mat;
N.mat0 = mat;
N.mat_intent  = 'Aligned';
N.mat0_intent = 'Aligned';
N.intent.code = 0;
N.intent.name = 'Mapping';
N.descrip     = 'Jacobian Determinants';
create(N);
N.dat(:,:,:,1,1) = Dets;

%==========================================================================
% function out = pull_def(Def,mat,job)
%==========================================================================
function out = pull_def(Def,mat,job)

PI      = job.fnames;
intrp   = job.interp;
intrp   = [intrp*[1 1 1], 0 0 0];
out     = cell(numel(PI),1);

if numel(PI)==0, return; end

if job.mask
    oM  = zeros(4,4);
    odm = zeros(1,3);
    dim = size(Def);
    msk = true(dim);
    for m=1:numel(PI)
        [pth,nam,ext,num] = spm_fileparts(PI{m});
        NI = nifti(fullfile(pth,[nam ext]));
        dm = NI.dat.dim(1:3);
        if isempty(num)
            j_range = 1:size(NI.dat,4);
        else
            num     = sscanf(num,',%d');
            j_range = num(1);
        end
        for j=j_range

            M0 = NI.mat;
            if ~isempty(NI.extras) && isstruct(NI.extras) && isfield(NI.extras,'mat')
                M1 = NI.extras.mat;
                if size(M1,3) >= j && sum(sum(M1(:,:,j).^2)) ~=0
                    M0 = M1(:,:,j);
                end
            end
            M   = inv(M0);
            if ~all(M(:)==oM(:)) || ~all(dm==odm)
                tmp = affine(Def,M);
                msk = tmp(:,:,:,1)>=1 & tmp(:,:,:,1)<=size(NI.dat,1) ...
                    & tmp(:,:,:,2)>=1 & tmp(:,:,:,2)<=size(NI.dat,2) ...
                    & tmp(:,:,:,3)>=1 & tmp(:,:,:,3)<=size(NI.dat,3);
            end
            oM  = M;
            odm = dm;
        end
    end
end

oM = zeros(4,4);
spm_progress_bar('Init',numel(PI),'Resampling','volumes completed');
for m=1:numel(PI)

    % Generate headers etc for output images
    %----------------------------------------------------------------------
    [pth,nam,ext,num] = spm_fileparts(PI{m});
    NI = nifti(fullfile(pth,[nam ext]));

    j_range = 1:size(NI.dat,4);
    k_range = 1:size(NI.dat,5);
    l_range = 1:size(NI.dat,6);
    if ~isempty(num)
        num = sscanf(num,',%d');
        if numel(num)>=1, j_range = num(1); end
        if numel(num)>=2, k_range = num(2); end
        if numel(num)>=3, l_range = num(3); end
    end

    NO = NI;
    if isfield(job.savedir,'savepwd')
        wd = pwd;
    elseif isfield(job.savedir,'saveusr')
        wd = job.savedir.saveusr{1};
    elseif isfield(job.savedir,'savesrc')
        wd = pth;
    else
        wd = pwd;
    end

    if sum(job.fwhm.^2)==0
        newprefix  = spm_get_defaults('normalise.write.prefix');
        NO.descrip = sprintf('Warped');
    else
        newprefix  = [spm_get_defaults('smooth.prefix') spm_get_defaults('normalise.write.prefix')];
        NO.descrip = sprintf('Smoothed (%gx%gx%g subopt) warped',job.fwhm);
    end
    if isfield(job,'prefix') && ~isempty(job.prefix)
        NO.dat.fname = fullfile(wd,[job.prefix nam ext]);
    else
        NO.dat.fname = fullfile(wd,[newprefix nam ext]);
    end
    dim            = size(Def);
    dim            = dim(1:3);
    NO.dat.dim     = [dim NI.dat.dim(4:end)];
    NO.dat.offset  = 0; % For situations where input .nii images have an extension.
    NO.mat         = mat;
    NO.mat0        = mat;
    NO.mat_intent  = 'Aligned';
    NO.mat0_intent = 'Aligned';
    if isempty(num)
        out{m}     = NO.dat.fname;
    else
        out{m}     = [NO.dat.fname, ',', num2str(num(1))];
    end
    NO.extras      = [];
    create(NO);

    % Smoothing settings
    vx  = sqrt(sum(mat(1:3,1:3).^2));
    krn = max(job.fwhm./vx,0.25);

    % Loop over volumes within the file
    %----------------------------------------------------------------------
    %fprintf('%s',nam);
    for j=j_range

        M0 = NI.mat;
        if ~isempty(NI.extras) && isstruct(NI.extras) && isfield(NI.extras,'mat')
            M1 = NI.extras.mat;
            if size(M1,3) >= j && sum(sum(M1(:,:,j).^2)) ~=0
                M0 = M1(:,:,j);
            end
        end
        M  = inv(M0);
        if ~all(M(:)==oM(:))
            % Generate new deformation (if needed)
            Y     = affine(Def,M);
        end
        oM = M;
        % Write the warped data for this time point
        %------------------------------------------------------------------
        for k=k_range
            for l=l_range
                C   = spm_diffeo('bsplinc',single(NI.dat(:,:,:,j,k,l)),intrp);
                dat = spm_diffeo('bsplins',C,Y,intrp);
                if job.mask
                    dat(~msk) = NaN;
                end
                if sum(job.fwhm.^2)~=0
                    spm_smooth(dat,dat,krn); % Side effects
                end
                NO.dat(:,:,:,j,k,l) = dat;
                %fprintf('\t%d,%d,%d', j,k,l);
            end
        end
    end
    %fprintf('\n');
    spm_progress_bar('Set',m);
end
spm_progress_bar('Clear');

%==========================================================================
% function out = push_def(Def,mat,job)
%==========================================================================
function out = push_def(Def,mat,job)
% Generate deformation, which is the inverse of the usual one (it is for "pushing"
% rather than the usual "pulling"). This deformation is affine transformed to
% allow for different voxel sizes and bounding boxes, and also to incorporate
% the affine mapping between MNI space and the population average shape.
%--------------------------------------------------------------------------

% Deal with desired bounding box and voxel sizes.
%--------------------------------------------------------------------------
if isfield(job.fov,'file')
    N1   = nifti(job.fov.file);
    mat0 = N1.mat;
    dim  = N1.dat.dim(1:3);
else
    bb   = job.fov.bbvox.bb;
    vox  = job.fov.bbvox.vox;
    [mat0, dim] = spm_get_matdim('', vox, bb);
end

M   = inv(mat0);
y0  = affine(Def,M);

if isfield(job,'weight') && ~isempty(job.weight) && ~isempty(job.weight{1})
    wfile = job.weight{1};
    Nw    = nifti(wfile);
    Mw    = Nw.mat;
    wt    = Nw.dat(:,:,:,1,1,1);
else
    wt    = [];
end

odm = zeros(1,3);
oM  = zeros(4,4);
PI  = job.fnames;
out = cell(numel(PI),1);
for m=1:numel(PI)

    % Generate headers etc for output images
    %----------------------------------------------------------------------
    [pth,nam,ext,num] = spm_fileparts(PI{m});
    NI = nifti(fullfile(pth,[nam ext]));
    j_range = 1:size(NI.dat,4);
    k_range = 1:size(NI.dat,5);
    l_range = 1:size(NI.dat,6);
    if ~isempty(num)
        num = sscanf(num,',%d');
        if numel(num)>=1, j_range = num(1); end
        if numel(num)>=2, k_range = num(2); end
        if numel(num)>=3, l_range = num(3); end
    end

    if isfield(job.savedir,'savepwd')
        wd = pwd;
    elseif isfield(job.savedir,'saveusr')
        wd = job.savedir.saveusr{1};
    elseif isfield(job.savedir,'savesrc')
        wd = pth;
    else
        wd = pwd;
    end

    NO = NI;
    if job.preserve
        NO.dat.scl_slope = 1.0;
        NO.dat.scl_inter = 0.0;
        NO.dat.dtype     = 'float32-le';
        if sum(job.fwhm.^2)==0
            newprefix  = [spm_get_defaults('deformations.modulate.prefix') spm_get_defaults('normalise.write.prefix')];
            NO.descrip = sprintf('Warped & Jac scaled');
        else
            newprefix  = [spm_get_defaults('smooth.prefix') spm_get_defaults('deformations.modulate.prefix') spm_get_defaults('normalise.write.prefix')];
            NO.descrip = sprintf('Smoothed (%gx%gx%g) warped Jac scaled',job.fwhm);
        end
    else
        if sum(job.fwhm.^2)==0
            newprefix  = spm_get_defaults('normalise.write.prefix');
            NO.descrip = sprintf('Warped');
        else
            newprefix  = [spm_get_defaults('smooth.prefix') spm_get_defaults('normalise.write.prefix')];
            NO.descrip = sprintf('Smoothed (%gx%gx%g opt) warped',job.fwhm);
        end
    end
    if isfield(job,'prefix') && ~isempty(job.prefix)
        NO.dat.fname = fullfile(wd,[job.prefix nam ext]);
    else
        NO.dat.fname = fullfile(wd,[newprefix nam ext]);
    end
    NO.dat.dim     = [dim NI.dat.dim(4:end)];
    NO.dat.offset  = 0; % For situations where input .nii images have an extension.
    NO.mat         = mat0;
    NO.mat0        = mat0;
    NO.mat_intent  = 'Aligned';
    NO.mat0_intent = 'Aligned';

    if isempty(num)
        out{m}     = NO.dat.fname;
    else
        out{m}     = [NO.dat.fname, ',', num2str(num(1))];
    end

    NO.extras      = [];
    create(NO);

    % Smoothing settings
    vx  = sqrt(sum(mat0(1:3,1:3).^2));
    krn = max(job.fwhm./vx,0.25);

    % Loop over volumes within the file
    %----------------------------------------------------------------------
    fprintf('%s',nam); drawnow;
    for j=j_range

        % Need to resample the mapping by an affine transform
        % so that it maps from voxels in the native space image
        % to voxels in the spatially normalised image.
        %------------------------------------------------------------------
        M0 = NI.mat;
        if ~isempty(NI.extras) && isstruct(NI.extras) && isfield(NI.extras,'mat')
            M1 = NI.extras.mat;
            if size(M1,3) >= j && sum(sum(M1(:,:,j).^2)) ~=0
                M0 = M1(:,:,j);
            end
        end

        M   = mat\M0;
        dm  = [size(NI.dat),1,1,1,1];
        if ~all(dm(1:3)==odm) || ~all(M(:)==oM(:))
            % Generate new deformation (if needed)
            y   = zeros([dm(1:3),3],'single');
            for d=1:3
                yd = y0(:,:,:,d);
                for x3=1:size(y,3)
                    y(:,:,x3,d) = single(spm_slice_vol(yd,M*spm_matrix([0 0 x3]),dm(1:2),[1 NaN]));
                end
            end
        end

        odm = dm(1:3);
        oM  = M;
        % Write the warped data for this time point.
        %------------------------------------------------------------------
        for k=k_range
            for l=l_range
                f  = single(NI.dat(:,:,:,j,k,l));
                if isempty(wt)
                    if ~job.preserve
                        % Unmodulated - note the slightly novel procedure
                        [f,c] = spm_diffeo('push',f,y,dim);
                        spm_smooth(f,f,krn); % Side effects
                        spm_smooth(c,c,krn); % Side effects
                        f = f./(c+0.001);
                    else
                        % Modulated, by pushing
                        scal = abs(det(NI.mat(1:3,1:3))/det(NO.mat(1:3,1:3))); % Account for vox sizes
                        f    = spm_diffeo('push',f,y,dim)*scal;
                        spm_smooth(f,f,krn); % Side effects
                    end
                else
                    if isequal(size(wt),size(f)) && sum((Mw(:)-M0(:)).^2)<1e-6
                        wtw = single(wt);
                        f   = single(f.*wt);
                    else
                        wtw = zeros(size(f),'single');
                        for z=1:size(wt,3)
                            Mz = Mw\M0*[1 0 0 0; 0 1 0 0; 0 0 1 z; 0 0 0 1];
                            wtw(:,:,z) = single(spm_slice_vol(wt,Mz,[size(f,1),size(f,2)],1));
                        end
                    end
                    if ~job.preserve
                        % Unmodulated - note the slightly novel procedure
                        f = spm_diffeo('push',f.*wtw,y,dim);
                        c = spm_diffeo('push',wtw,y,dim);
                        spm_smooth(f,f,krn); % Side effects
                        spm_smooth(c,c,krn); % Side effects
                        f = f./(c+0.001);
                    else
                        % Modulated, by pushing
                        scal = abs(det(NI.mat(1:3,1:3))/det(NO.mat(1:3,1:3))); % Account for vox sizes
                        f    = spm_diffeo('push',f.*wtw,y,dim)*scal;
                        spm_smooth(f,f,krn); % Side effects
                    end
                    clear wtw
                end
                NO.dat(:,:,:,j,k,l) = f;
                fprintf('\t%d,%d,%d', j,k,l); drawnow;
            end
        end
    end
    fprintf('\n'); drawnow;
end

%==========================================================================
% function out = surf_def(Def,mat,job)
%==========================================================================
function out = surf_def(Def,mat,job)
filenames = job.surface;
out       = cell(numel(filenames),1);
for i=1:numel(filenames)
    fname     = deblank(job.surface{i});
    [pth,nam] = fileparts(fname);
    fprintf('%s\n', nam);
    d      = size(Def);
    tmp    = double(reshape(Def,[d(1:3) 1 d(4)]));
    Tmesh  = spm_swarp(fname, tmp,mat);
    if isfield(job.savedir,'savepwd')
        wd = pwd;
    elseif isfield(job.savedir,'saveusr')
        wd = job.savedir.saveusr{1};
    elseif isfield(job.savedir,'savesrc')
        wd = pth;
    else
        wd = pwd;
    end
    filename = fullfile(wd,[nam,'_warped', '.gii']);
    save(gifti(Tmesh), filename);
    out{i}   = filename;
end

%==========================================================================
% function Def = affine(y,M)
%==========================================================================
function Def = affine(y,M)
Def          = zeros(size(y),'single');
Def(:,:,:,1) = y(:,:,:,1)*M(1,1) + y(:,:,:,2)*M(1,2) + y(:,:,:,3)*M(1,3) + M(1,4);
Def(:,:,:,2) = y(:,:,:,1)*M(2,1) + y(:,:,:,2)*M(2,2) + y(:,:,:,3)*M(2,3) + M(2,4);
Def(:,:,:,3) = y(:,:,:,1)*M(3,1) + y(:,:,:,2)*M(3,2) + y(:,:,:,3)*M(3,3) + M(3,4);

%==========================================================================
% function Def = identity(d,M)
%==========================================================================
function Def = identity(d,M)
[y1,y2]   = ndgrid(single(1:d(1)),single(1:d(2)));
Def       = zeros([d 3],'single');
for y3=1:d(3)
    Def(:,:,y3,1) = y1*M(1,1) + y2*M(1,2) + (y3*M(1,3) + M(1,4));
    Def(:,:,y3,2) = y1*M(2,1) + y2*M(2,2) + (y3*M(2,3) + M(2,4));
    Def(:,:,y3,3) = y1*M(3,1) + y2*M(3,2) + (y3*M(3,3) + M(3,4));
end


text:
% @Author: aaronmishkin
% @Date:   2018-07-26T13:13:16-07:00
% @Email:  amishkin@cs.ubc.ca
% @Last modified by:   aaronmishkin
% @Last modified time: 2018-07-26T13:22:15-07:00

% Load the toy example experiment data.
file_name = strcat('./toy_example_experiment_data.mat');
load(file_name)

method = {'Vadam', 'VOGN-1'}
colors = [0 0.8 0.5; 0 0 1; 1 0 0];
f = figure('Position', [50,50,1000,715]); clf;
% exact post
contourf(w1,w2,reshape(post,[n,n]),5);
cbh=colorbar;
colormap(gray);
hold on
% map estimate
h(1) = plot(wmap(1),wmap(2),'+','color', 'k', 'MarkerSize',7, 'linewidth', 8, 'markersize', 7);
h(2) = plot_gaussian_ellipsoid(w_exact_vi, C_exact_vi, 1);
set(h(2), 'color', 0*[1 1 1], 'linestyle', '-.', 'linewidth', 12);

for m = 1:length(method)
    plot(w_all(1,t,m), w_all(2,t,m), 'or', 'color', colors(m,:), 'linewidth', 8, 'markerfacecolor', colors(m,:), 'markersize', 5);
    h(m+2) = plot_gaussian_ellipsoid(w_all(:,t,m), Sigma_all(:,:,t,m), 1);
    set(h(m+2), 'color', colors(m,:), 'linewidth', 12);
end
% mf-exact
axis([0 20 0 12]);
plot(w_exact_vi(1), w_exact_vi(2), 'o', 'color', 0*[1 1 1], 'linewidth', 8, 'markerfacecolor', 0*[1 1 1], 'markersize', 5);
      t = maxIters;

hl = legend(h, {'MAP', 'VI-Exact', 'Vadam', 'VOGN-1'}, 'location', 'northwest');
hx = xlabel('\theta_1');
hy = ylabel('\theta_2');
set(gca, 'fontsize', 24);
set([hx,hy], 'fontsize', 24, 'fontname', 'helvetica');
set(hl, 'fontsize', 24, 'fontname', 'helvetica');
set(gca, 'xtick', [0:5:20], 'ytick', [0:5:10], 'tickdir', 'out');
set(cbh,'YTick',[0:1e-3:5e-3])

h = gcf;

f.Position(3) = 1000;
f.Position(4) = 800;

axesObjs = get(h, 'Children');  %axes handles
dataObjs = get(axesObjs, 'Children'); %handles to low-level graphics objects in axes

cell_obj = dataObjs(3, 1);
cell_obj = cell_obj{1};

vogn_dot = cell_obj(3);
vogn_line = cell_obj(2);

vadam_dot = cell_obj(5);
vadam_line = cell_obj(4);

vi_dot = cell_obj(1);
vi_line = cell_obj(6);
map_dot = cell_obj(7);

linewidth = 12;

set(vogn_dot,'color', 'b')
set(vogn_line,'color', 'b')
set(vogn_line,'linewidth', linewidth)

set(vadam_dot, 'color', 'r')
set(vadam_line,'color', 'r')
set(vadam_line,'linewidth', linewidth)

set(vi_dot,'color', [0, 1.0, 0])
set(vi_line,'color', [0, 1.0, 0])
set(vi_line,'linewidth', linewidth)

set(map_dot, 'markerSize', 22)
set(map_dot, 'color', 'black')

%%%%
fontsize = 30;
ax = gca;
set(ax, 'fontsize', fontsize);

ax.XLim = [0, 20];
ax.YLim = [0, 11];
ax.XLabel.String = 'Weight 1';
ax.YLabel.String = 'Weight 2';
%%%
hLegend = findobj(gcf, 'Type', 'Legend');
hLegend.FontSize = 30;
%set(ax, 'LooseInset', get(ax, 'TightInset'));

hLegend.String{2} = 'MF-Exact';

grid off

set(f,'Units','Inches');
pos = get(f,'Position');
set(f,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])

savefig('figures/figure_two_a.fig')
saveas(f, 'figures/figure_two_a.pdf')
close all

author: emtiyaz

repo: vadam

save_path: github-repos/MATLAB/emtiyaz-vadam


text:
function [sFile, ChannelMat] = in_fopen_manscan(DataFile)
% IN_FOPEN_MANSCAN: Open a MANSCAN file (continuous recordings)
%
% USAGE:  [sFile, ChannelMat] = in_fopen_manscan(DataFile)

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

%% ===== READ HEADER =====
% Get text file (.mbi)
MbiFile = strrep(DataFile, '.mb2', '.mbi');
% If doesn't exist: error
if ~file_exist(MbiFile)
    error('Cannot open file: missing text file .mbi');
end
% Initialize header
iEpoch = 1;
hdr.epoch(iEpoch).Comment = {};
hdr.epoch(iEpoch).Channel = [];
hdr.Events = [];
curBlock = '';
% Read file line by line
fid = fopen(MbiFile,'r');
while(1)
    % Reached the end of the file: exit the loop
    if feof(fid)
        break;
    end
    % Read one line
    buf = fgetl(fid);
    if isempty(buf)
        curBlock = '';
        continue;
    end
    % Split line based on space characters
    splitBuf = str_split(buf, [0 9 32]);
    if isempty(splitBuf)
        continue;
    end
    % Current block
    if ~isempty(curBlock)
        switch curBlock
            case 'event'
                % Ignore MANSCAN auto events
                if any(strcmpi(evtProp, 'Channel'))
                    continue;
                end
                % New event
                iEvt = length(hdr.Events) + 1;
                % Read event name and epoch
                hdr.Events(iEvt).Name   = splitBuf{1};
                hdr.Events(iEvt).iEpoch = iEpoch;
                % Read each property
                for iProp = 2:length(evtProp)
                    if strcmpi(evtProp{iProp}, 'Channel')
                        hdr.Events(iEvt).Name = [hdr.Events(iEvt).Name '_' splitBuf{iProp}];
                    elseif strcmpi(evtType{iProp}, 'String')
                        hdr.Events(iEvt).(evtProp{iProp}) = splitBuf{iProp};
                    else
                        hdr.Events(iEvt).(evtProp{iProp}) = str2num(splitBuf{iProp});
                    end
                end

            case 'channel'
                % New channel
                iChan = length(hdr.epoch(iEpoch).Channel) + 1;
                % Read channel name
                hdr.epoch(iEpoch).Channel(iChan).Name = splitBuf{1};
                % Read each property
                for iProp = 2:length(chanProp)
                    if strcmpi(chanProp{iProp}, 'UnitDisplayToFile')
                        propValue = str2num(splitBuf{iProp});
                    else
                        propValue = splitBuf{iProp};
                    end
                    hdr.epoch(iEpoch).Channel(iChan).(chanProp{iProp}) = propValue;
                end
        end
        % Next line
        continue;
    end
    % Detect keyword
    switch lower(splitBuf{1})
        case 'datatypeid'
            iEpoch = iEpoch + 1;
            hdr.epoch(iEpoch).Comment = {};
            hdr.epoch(iEpoch).Channel = [];
        case 'comment'
            hdr.epoch(iEpoch).Comment{end+1} = strtrim(strrep(buf, splitBuf{1}, ''));
        case 'event'
            curBlock = 'event';
            evtProp = splitBuf;
            evtType = str_split(fgetl(fid), [0 9 32]);
        case 'channel'
            curBlock = 'channel';
            chanProp = splitBuf;
        case 'offsetdisplayexperiment'
            hdr.epoch(iEpoch).OffsetDisplay = str2num(splitBuf{2});
        case 'experimenttime'
            hdr.epoch(iEpoch).ExperimentTime = strtrim(strrep(buf, splitBuf{1}, ''));
        case 'worddatafile'
            % Read channel order
            hdr.epoch(iEpoch).ChannelOrder = str_split(fgetl(fid), [0 9 32]);
            % Read record info
            datainfo = str_split(fgetl(fid), [0 9 32]);
            hdr.epoch(iEpoch).StartData1 = str2num(datainfo{1});
            hdr.epoch(iEpoch).StartData2 = str2num(datainfo{2});
            hdr.epoch(iEpoch).Sweeps     = str2num(datainfo{3});
            hdr.epoch(iEpoch).BinFile    = datainfo{4};
    end
end

%% ===== CREATE BRAINSTORM SFILE STRUCTURE =====
% Initialize returned file structure
sFile = db_template('sfile');
% Add information read from header
sFile.byteorder  = 'l';
sFile.filename   = DataFile;
sFile.format = 'EEG-MANSCAN';
sFile.device = 'MANSCAN';
% Comment: short filename
[tmp__, sFile.comment, tmp__] = bst_fileparts(DataFile);
% Consider that the sampling rate of the file is the sampling rate of the first signal
sFile.prop.sfreq = 256;
sFile.prop.nAvg  = 1;
% No info on bad channels
sFile.channelflag = ones(length(hdr.epoch(iEpoch).ChannelOrder), 1);
% Acquisition date
try
    sFile.acq_date = str_date(hdr.epoch.ExperimentTime, 'mm/dd/yy');
catch
end

%% ===== EPOCHS =====
if (length(hdr.epoch) <= 1)
    sFile.prop.times   = [0, hdr.epoch(1).Sweeps-1] ./ sFile.prop.sfreq;
else
    % Build epochs structure
    for iEpoch = 1:length(hdr.epoch)
        sFile.epochs(iEpoch).label   = sprintf('Epoch #%02d', iEpoch);
        sFile.epochs(iEpoch).times   = [0, hdr.epoch(iEpoch).Sweeps-1] ./ sFile.prop.sfreq;
        sFile.epochs(iEpoch).nAvg    = 1;
        sFile.epochs(iEpoch).select  = 1;
        sFile.epochs(iEpoch).bad     = 0;
        sFile.epochs(iEpoch).channelflag = [];

        % Check if all the epochs have the same channel list
        if (iEpoch > 1) && ~isequal(hdr.epoch(iEpoch).ChannelOrder, hdr.epoch(1).ChannelOrder)
            error('Channel list must remain constant across epochs.');
        end
    end
    % Extract global min/max for time and samples indices
    sFile.prop.times   = [min([sFile.epochs.times]),   max([sFile.epochs.times])];
end

%% ===== CREATE EMPTY CHANNEL FILE =====
nChannels = length(hdr.epoch(1).ChannelOrder);
ChannelMat = db_template('channelmat');
ChannelMat.Comment = [sFile.device ' channels'];
ChannelMat.Channel = repmat(db_template('channeldesc'), [1, nChannels]);
hdr.Gains = [];
% For each channel
for iChan = 1:nChannels
    % Find channel in description list
    chName = hdr.epoch(1).ChannelOrder{iChan};
    iDesc = find(strcmpi({hdr.epoch(1).Channel.Name}, chName));
    sDesc = hdr.epoch(1).Channel(iDesc);
    % Type
    if isfield(sDesc, 'IsEEG') && ~isempty(sDesc.IsEEG)
        if strcmpi(sDesc.IsEEG, 'Yes')
            ChannelMat.Channel(iChan).Type = 'EEG';
        else
            ChannelMat.Channel(iChan).Type = 'Misc';
        end
    elseif isfield(sDesc, 'IsEEG')
        ChannelMat.Channel(iChan).Type = 'EEG REF';
    else
        ChannelMat.Channel(iChan).Type = 'EEG';
    end
    % Create structure
    ChannelMat.Channel(iChan).Name    = chName;
    ChannelMat.Channel(iChan).Loc     = [0; 0; 0];
    ChannelMat.Channel(iChan).Orient  = [];
    ChannelMat.Channel(iChan).Weight  = 1;
    ChannelMat.Channel(iChan).Comment = '';
    % Check that this channel has a gain
    if isempty(hdr.epoch(1).Channel(iDesc).UnitDisplayToFile)
        hdr.Gains(iChan,1) = 1;
    else
        hdr.Gains(iChan,1) = hdr.epoch(1).Channel(iDesc).UnitDisplayToFile;
    end
end

%% ===== FORMAT EVENTS =====
if ~isempty(hdr.Events)
    % Get all the epochs names
    uniqueEvt = unique({hdr.Events.Name});
    % Create one category for each event
    for iEvt = 1:length(uniqueEvt)
        % Get all the occurrences
        iOcc = find(strcmpi({hdr.Events.Name}, uniqueEvt{iEvt}));
        % Get the samples for all the occurrences
        sample   = [hdr.Events(iOcc).BeginSample];
        duration = [hdr.Events(iOcc).DurationSamples];
        % Extended event = duration is not 1 for all the markers
        if ~all(duration == 1)
            sample = [sample; sample + duration];
        end
        % Create event structure
        sFile.events(iEvt).label    = uniqueEvt{iEvt};
        sFile.events(iEvt).times    = sample ./ sFile.prop.sfreq;
        sFile.events(iEvt).epochs   = [hdr.Events(iOcc).iEpoch];
        sFile.events(iEvt).select   = 1;
        sFile.events(iEvt).channels = [];
        sFile.events(iEvt).notes    = [];
    end
end
% Save file header
sFile.header = hdr;

author: brainstorm-tools


text:
function tblstring = makeSAMMIdataString(tbl)
    tblstring = '[[''';

    %Headers
    tblstring = strcat(tblstring,strjoin(tbl.Properties.VariableNames,''','''),'''],');

    %Data
    dat = table2array(tbl);

    %Make key string
    fun = @(a) strjoin(sprintfc('%g',a{1}),'","');
    fd = arrayfun(fun,num2cell(dat,2),'UniformOutput',false);
    fd = strcat('["',tbl.Properties.RowNames,'","',fd,'"]');
    tblstring = strcat(tblstring,strjoin(fd,','),']');
end

author: opencobra

repo: cobratoolbox


text:
background.marker_pos(:,1+mod(background.smax:background.smax+size(background_marker_pos,2)-1,background.buffer_len)) = background_marker_pos+logical(background_marker_pos)*background.mmax;
background.mmax = background.mmax + nnz(background_marker_pos);

repo: REST


text:
function pck = computePCK(dist,range)

pck = zeros(numel(range),size(dist,2)+1);

for jidx = 1:size(dist,2)
    % compute PCK for each threshold
    for k = 1:numel(range)
        pck(k,jidx) = 100*mean(squeeze(dist(1,jidx,:))<=range(k));
    end
end

% compute average PCK
for k = 1:numel(range)
    pck(k,end) = 100*mean(reshape(squeeze(dist(1,:,:)),size(dist,2)*size(dist,3),1)<=range(k));
end

end

author: Guanghan

repo: GNet-pose

save_path: github-repos/MATLAB/Guanghan-GNet-pose


text:
%CALCOPTICALFLOWSPARSETODENSE  Fast dense optical flow based on PyrLK sparse matches interpolation
%
%     flow = cv.calcOpticalFlowSparseToDense(from, to)
%     flow = cv.calcOpticalFlowSparseToDense(from, to, 'OptionName',optionValue, ...)
%
% ## Input
% * __from__ first 8-bit 3-channel or 1-channel image.
% * __to__ second 8-bit 3-channel or 1-channel image of the same size as
%   `from`.
%
% ## Output
% * __flow__ computed flow image that has the same size as `from` and type
%   `single` (2-channels). Flow for `(x,y)` is stored in the third dimension.
%
% ## Options
% * __GridStep__ stride used in sparse match computation. Lower values usually
%   result in higher quality but slow down the algorithm. default 8
% * __K__ number of nearest-neighbor matches considered, when fitting a
%   locally affine model. Lower values can make the algorithm noticeably
%   faster at the cost of some quality degradation. default 128
% * __Sigma__ parameter defining how fast the weights decrease in the
%   locally-weighted affine fitting. Higher values can help preserve fine
%   details, lower values can help to get rid of the noise in the output flow.
%   default 0.05
% * __UsePostProcessing__ defines whether the
%   cv.FastGlobalSmootherFilter.fastGlobalSmootherFilter is used for
%   post-processing after interpolation. default true
% * __FGSLambda__ see the respective parameter of
%   cv.FastGlobalSmootherFilter.fastGlobalSmootherFilter, default 500.0
% * __FGSSigma__ see the respective parameter of
%   cv.FastGlobalSmootherFilter.fastGlobalSmootherFilter, default 1.5
%
% See also: cv.calcOpticalFlowPyrLK, cv.SparsePyrLKOpticalFlow
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function match = at_dense_tc(desc1,desc2)

[idx12, dis12] = yael_nn(desc2, desc1, 1);
[idx21, dis21] = yael_nn(desc1, desc2, 1);

Ndesc = length(idx12);
match = NaN(3,Ndesc);
for ii=1:Ndesc
  if ~isnan(idx12(1,ii))
    if idx21(1,idx12(1,ii)) == ii
      match(:,ii) = [ii; single(idx12(1,ii)); dis12(1,ii)];
    end
  end
end
match = match(:,~isnan(match(1,:)));

author: HajimeTaira


text:
function B = transform_image_LO_TO(I, corners, interporation)

% B = transform_image3(I, corners, interporation)
% calculate the transformation matrix for each area

B = nan(128,128);

if notDefined('interporation'), interporation = [];  end;

% LO1
input_points = corners{1};
% base_points =  [40 20;40 35; 80 35; 80 20];
base_points =  [30 80; 45 80; 45 40; 30 40];
udata = [1 121];  vdata = [1 121];  % input coordinate system
tform = maketform('projective',input_points, base_points);
[tmp,xdata,ydata] = imtransform(I,tform,interporation,'udata',udata,...
                                                'vdata',vdata,...
                                                'xdata',udata,...
                                                'ydata',vdata,...
                                                'size',size(I),...
                                                'fill',128);

B(round(20/120*128):round(100/120*128), round(20/120*128):round(45/120*128))=tmp(round(20/120*128):round(100/120*128),round(20/120*128):round(45/120*128));

% LO2
input_points = corners{2};
% base_points =  [40 50;40 35; 80 35; 80 50];
base_points =  [60 80; 45 80; 45 40; 60 40];
udata = [1 121];  vdata = [1 121];  % input coordinate system
tform = maketform('projective',input_points, base_points);
[tmp,xdata,ydata] = imtransform(I,tform,interporation,'udata',udata,...
                                                'vdata',vdata,...
                                                'xdata',udata,...
                                                'ydata',vdata,...
                                                'size',size(I),...
                                                'fill',128);
B(round(20/120*128):round(100/120*128), round(45/120*128):round(60/120*128))=tmp(round(20/120*128):round(100/120*128),round(45/120*128):round(60/120*128));

% TO1
input_points = corners{3};
% base_points =  [40 50;40 65; 80 65; 80 50];
base_points =  [60 80; 75 80; 75 40; 60 40];
udata = [1 121];  vdata = [1 121];  % input coordinate system
tform = maketform('projective',input_points, base_points);
[tmp,xdata,ydata] = imtransform(I,tform,interporation,'udata',udata,...
                                                'vdata',vdata,...
                                                'xdata',udata,...
                                                'ydata',vdata,...
                                                'size',size(I),...
                                                'fill',128);

B(round(20/120*128):round(100/120*128), round(60/120*128):round(75/120*128))=tmp(round(20/120*128):round(100/120*128),round(60/120*128):round(75/120*128));

% TO2
input_points = corners{4};
base_points =  [90 80; 75 80; 75 40; 90 40];
udata = [1 121];  vdata = [1 121];  % input coordinate system
tform = maketform('projective',input_points, base_points);
[tmp,xdata,ydata] = imtransform(I,tform,interporation,'udata',udata,...
                                                'vdata',vdata,...
                                                'xdata',udata,...
                                                'ydata',vdata,...
                                                'size',size(I),...
                                                'fill',128);

B(round(20/120*128):round(100/120*128), round(75/120*128):round(90/120*128))=tmp(round(20/120*128):round(100/120*128),round(75/120*128):round(90/120*128));

% TO2 side (for the control of noise average)
input_points = corners{6};
base_points =  [90 80; 105 80; 105 40; 90 40];
udata = [1 121];  vdata = [1 121];  % input coordinate system
tform = maketform('projective',input_points, base_points);
[tmp,xdata,ydata] = imtransform(I,tform,interporation,'udata',udata,...
                                                'vdata',vdata,...
                                                'xdata',udata,...
                                                'ydata',vdata,...
                                                'size',size(I),...
                                                'fill',128);

B(round(20/120*128):round(100/120*128), round(90/120*128):round(105/120*128))=tmp(round(20/120*128):round(100/120*128),round(90/120*128):round(105/120*128));

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function fgplvmCmu35Animate(experimentNo, missingData)

% FGPLVMCMU35ANIMATE Animate the test data jointly with predictions.
% FORMAT
% DESC animates the stick figure using 'fill in' from the GP-LVM models.
% ARG experimentNo : the experiment number (1 is four dimensions, 2 is three
% dimensions, 3 is five dimensions). Default value is 1.
% ARG missingData : either 'leg' or 'body' for missing data being the leg
% or the upper body. Default is 'body'.
%
% COPYRIGHT : Neil D. Lawrence, 2008
%
% SEEALSO : demCmu35gplvm1, demCmu35SequenceOptimisef

% FGPLVM

if nargin < 2
  missingData = 'body';
  if nargin < 1
    experimentNo = 1;
  end
end

dataSetName = 'cmu35gplvm';

% load skeleton
skel = acclaimReadSkel('35.asf');
[tmpchan, skel] = acclaimLoadChannels('35_01.amc', skel);

[Y, lbls, Ytest, lblstest] = lvmLoadData(dataSetName);
channelTest = Ytochannels(Ytest);
capName = dataSetName;
capName(1) = upper(dataSetName(1));
missingData(1) = upper(missingData(1));
load(['dem' capName 'Yvals' missingData num2str(experimentNo) '.mat'])
channelPred = Ytochannels(Ypred);

function channels = Ytochannels(Y)

% YTOCHANNELS Convert Y to channel values.

xyzInd = [2];
xyzDiffInd = [1 3];
rotInd = [4 6];
rotDiffInd = [5];
generalInd = [7:38 41:47 49:50 53:59 61:62];
startInd = 1;
endInd = length(generalInd);
channels(:, generalInd) = 180*Y(:, startInd:endInd)/pi;
startInd = endInd + 1;
endInd = endInd + length(xyzDiffInd);
channels(:, xyzDiffInd) = cumsum(Y(:, startInd:endInd), 1);
startInd = endInd + 1;
endInd = endInd + length(xyzInd);
channels(:, xyzInd) = Y(:, startInd:endInd);
startInd = endInd + 1;
endInd = endInd + length(xyzDiffInd);
channels(:, xyzDiffInd) = 0; %cumsum(Y(:, startInd:endInd), 1);
startInd = endInd + 1;
endInd = endInd + length(rotInd);
channels(:, rotInd) = real(asin(Y(:, startInd:endInd))*180/pi);
channels(:, rotInd(end)) = channels(:, rotInd(end))+270;
startInd = endInd + 1;
endInd = endInd + length(rotDiffInd);
channels(:, rotDiffInd) = 0;%cumsum(asin(Y(:, startInd:endInd)), 1))*180/pi;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function output = callquadprog(interfacedata)

options = interfacedata.options;
model = yalmip2quadprog(interfacedata);

if options.savedebug
    save quadprogdebug model
end

if options.showprogress;showprogress(['Calling ' interfacedata.solver.tag],options.showprogress);end
solvertime = tic;
solveroutput = callsolver(model,options);
solvertime = toc(solvertime);
solution = quadprogsol2yalmipsol(solveroutput,model);

% Save all data sent to solver?
if ~options.savesolverinput
    model = [];
end

% Save all data from the solver?
if ~options.savesolveroutput
    solveroutput = [];
end

% Standard interface
Primal      = solution.x(:);
Dual        = solution.D_struc;
problem     = solution.problem;
infostr     = yalmiperror(solution.problem,interfacedata.solver.tag);
if ~options.savesolverinput
    solverinput = [];
else
    solverinput = model;
end
if ~options.savesolveroutput
    solveroutput = [];
else
    solveroutput = solveroutput;
end
% Standard interface
output = createOutputStructure(Primal,Dual,[],problem,infostr,solverinput,solveroutput,solvertime);

function solveroutput = callsolver(model,options)
x = [];
fmin = [];
flag = [];
output = [];
lambda = [];
if nnz(model.Q) == 0
    % BUG in LIN/QUADPROG, computation of lambda crashes in some rare
    % cases. To avoid seeing this when we don't want the lambdas anyway, we
    % don't ask for it
    if options.saveduals
        [x,fmin,flag,output,lambda] = linprog(model.c, model.A, model.b, model.Aeq, model.beq, model.lb, model.ub, model.x0,model.ops);
    else
        lambda = [];
        [x,fmin,flag,output] = linprog(model.c, model.A, model.b, model.Aeq, model.beq, model.lb, model.ub, model.x0,model.ops);
    end
else
    if options.saveduals
        [x,fmin,flag,output,lambda] = quadprog(model.Q, model.c, model.A, model.b, model.Aeq, model.beq, model.lb, model.ub, model.x0,model.ops);
    else
        lambda = [];
        if options.savesolveroutput
            [x,fmin,flag,output] = quadprog(model.Q, model.c, model.A, model.b, model.Aeq, model.beq, model.lb, model.ub, model.x0,model.ops);
        else
            [x,fmin,flag] = quadprog(model.Q, model.c, model.A, model.b, model.Aeq, model.beq, model.lb, model.ub, model.x0,model.ops);
        end
    end
    if flag==5
        [x,fmin,flag,output,lambda] = quadprog(model.Q, model.c, model.A, model.b, model.Aeq, model.beq, model.lb, model.ub, [],model.ops);
    end
end
if isempty(x)
    x = nan(length(model.c),1);
end
solveroutput.x = x;
solveroutput.fmin = fmin;
solveroutput.flag = flag;
solveroutput.output = output;
solveroutput.lambda = lambda;

function solution = quadprogsol2yalmipsol(solveroutput,model)

solution.x = solveroutput.x(:);
% Internal format for duals
if ~isempty(solveroutput.lambda)
    solution.D_struc = [solveroutput.lambda.eqlin;solveroutput.lambda.ineqlin];
else
    solution.D_struc = [];
end

% Check, currently not exhaustive...
problem = 0;
if solveroutput.flag==0
    solution.problem = 3;
else
    if solveroutput.flag==-3
        solution.problem = 2;
    elseif solveroutput.flag==-2
        solution.problem = 1;
    else
        if solveroutput.flag>0
            solution.problem = 0;
        else
            if isempty(solveroutput.x)
                solution.x = repmat(nan,length(model.f),1);
            end
            if any((model.A*solveroutput.x-model.b)>sqrt(eps)) | any( abs(model.Aeq*solveroutput.x-model.beq)>sqrt(eps))
                solution.problem = 1; % Likely to be infeasible
            else
                if model.f'*solveroutput.x<-1e10 % Likely unbounded
                    solution.problem = 2;
                else          % Probably convergence issues
                    solution.problem = 5;
                end
            end
        end
    end
end


text:
% op_concatSubspecs.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% out=op_concatSubspecs(in1,in2);
%
% DESCRIPTION:
% Concatenate two scans along the subspecs dimension.  Two scans with 50
% averages each will now look like a single scan with 100 averages.
%
% INPUTS:
% in1    = first input in matlab structure format.
% in2    = second input in matlab structure format.
%
% OUTPUTS:
% out    = Output following concatenation along the subspecs dimension.

function out=op_concatSubspecs(in1,in2);

if in1.dims.subSpecs ~= in2.dims.subSpecs || in1.dims.t ~= in2.dims.t || in1.dims.coils ~= in2.dims.coils || in1.dims.averages ~=in2.dims.averages
    error('subSpecs dimensions must be the same for both inputs');
end

%if subspecs dimension is zero, make a new dimension for them
if in1.dims.subSpecs==0
    newSubspecsDim=max([in1.dims.t in1.dims.coils in1.dims.averages in1.dims.subSpecs in1.dims.extras])+1;
else
    newSubspecsDim=in1.dims.subSpecs;
end

fids=cat(newSubspecsDim,in1.fids,in2.fids);
specs=cat(newSubspecsDim,in1.specs,in2.specs);
sz=size(fids);

%FILLING IN DATA STRUCTURE
out=in1;
out.fids=fids;
out.specs=specs;
out.sz=sz;
out.dims=in1.dims;
out.dims.subSpecs=newSubspecsDim;

out.rawSubspecs=in1.rawSubspecs+in2.rawSubspecs;
out.subspecs=in1.subspecs+in2.subspecs;

%FILLING IN THE FLAGS
out.flags=in1.flags;
out.flags.writtentostruct=1;
out.flags.subtracted=0;

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
% PLOTPDFTEX prints a PDF figure for inclusion in LaTeX documents, with
%    LaTeX fonts. It creates a complete pdf-file with annotation of the
%    figure (titles, labels and texts) and graphics (lines,arrows,markers,
%    ...). This function uses the LAPRINT function and makes simply a
%    concatenation of the text part and the graphical part in a single PDF
%    file. It requires:
%
%                  MATLAB        |      OTHER
%               -----------------|----------------------------
%                 laprint.m      |      latex
%                                |      dvips
%                                |      ps2epsi or eps2eps
%                                |      epstopdf
%                                |      GhostScript
%
%    The diference between ps2epsi and eps2eps is that the first is a
%    standard EPS converter, but creates very large file size. The commande
%    eps2eps can be replace by ps2epsi at line 191
%
%    USAGE:
%    ------------------------
%
%    PLOTPDFTEX creates the PDF file of the current graphical figure (CGF)
%    named LaTeXfile.pdf
%
%    PLOTPDFTEX(H) creates the PDF file from the graphical figure with
%    handle H, named LaTeXfile.pdf
%
%    PLOTPDFTEX(H,FILENAME) creates the PDF file from the graphical figure
%    with handle H, named FILENAME.PDF (FILENAME is a character array of
%    the filename, with or without the extension .PDF)
%
%    PLOTPDFTEX(H,FILENAME,[FONTSCALING LINESCALING]) scales the font and
%    line sizes relative to the figure size (values greater than 1 scale
%    up, less than 1 scale down).  If nan, the default of 1 is used.
%
%    PLOTPDFTEX(H,FILENAME,[FONTSCALING LINESCALING],PACKAGES) to adding
%    packages such as fonts, mathfont,...
%    ex: PACKAGES = 'amssymb, times'
%
% See also LAPRINT:
%          http://www.uni-kassel.de/fb16/rat/matlab/laprint/laprintdoc.ps
%
%_______________________________
% Adjusted from PLOTEPSTEX by
% J.C. Olivier                                    09/2007
% http://www.jc-olivier.2007.fr
%
%_______________________________
%
% Emmett Lalish                                   06/2008

function plotpdftex(h,FileName,FontLineSize,Packages)

if nargin >= 3
    fontsize=FontLineSize(1);
    linesize=FontLineSize(2);

    id=find(isnan(FontLineSize));
    if length(find(id==1))==1
        fontsize=1;
    end

    if length(find(id==2))==1
        linesize=1;
    end
end

if nargin < 3
   fontsize=1; linesize=1;
end

if nargin < 2
    FileName='LaTeXfile';
end

if nargin < 1
    h=gcf;
end

TempName = strcat('TEMP',num2str(round(rand*10000))); %Generate random file name

DTI = get(0,'defaulttextinterpreter');   %Backup of the text interpreter
set(0,'defaulttextinterpreter','none');  %Delete the TeX interpreter

id=findstr(FileName,'.');
if ~isempty(id)
    FileName(id:end)=[];
end

if ~exist('laprint.m')
    disp('Laprint M-file does not exist or is not in the MATLAB''s search path.');
    disp('This file can be downloaded at: http://www.mathworks.com/matlabcentral/fileexchange');
    disp('                                                         Try again...');
    return;
end

laprint(h,TempName,'width',12/fontsize,'factor',linesize/fontsize,'scalefonts','off','asonscreen','on'); % Generate the EPS/LaTeX file

disp(sprintf('\n-------------------------------------------------------'));
disp(sprintf('PlotEpsTeX                                      02/2008'));
disp(sprintf('-------------------------------------------------------'));
disp(sprintf('Exporting figure in EPS format, including LaTeX strings'));
disp(sprintf('This function is based on laprint.m:'));
disp(sprintf('http://www.uni-kassel.de/fb16/rat/matlab/laprint/laprintdoc.ps\n'));
%-------------------------------------------------------
% Temporary LaTeX file
%-------------------------------------------------------
fid = fopen(strcat(TempName,'2.tex'),'w');

fprintf(fid,'\\documentclass[11pt, oneside]{article}\n');
fprintf(fid,'\\usepackage{graphicx}\n');
fprintf(fid,'\\usepackage{amsmath}\n');
fprintf(fid,'\\usepackage[T1]{fontenc}\n');
fprintf(fid,'\\usepackage[latin1]{inputenc}\n');
fprintf(fid,'\\usepackage{ae}\n');
fprintf(fid,'\\usepackage{psfrag}\n');
fprintf(fid,'\\usepackage{color}\n');

if nargin==4
    fprintf(fid,'\\usepackage{%s}\n',Packages); % Suplementary packages
end
fprintf(fid,'\\pagestyle{empty}\n');
fprintf(fid,' \n');
fprintf(fid,'\\begin{document}\n');
fprintf(fid,'    \\begin{figure}\n');
fprintf(fid,'        \\centering\n');
fprintf(fid,'        \\input{%s}\n',TempName);
fprintf(fid,'    \\end{figure}\n');
fprintf(fid,' \n');
fprintf(fid,'\\end{document}\n');
fclose(fid);

%-------------------------------------------------------
% LaTeX Command
%-------------------------------------------------------
Str=sprintf('latex --src -interaction=nonstopmode %s2.tex',TempName);
disp(sprintf('\n[LaTeX Command] %s',Str));
[hdos,wdos]=system(Str);

if hdos ~=0
    if isunix
        dos(sprintf('rm %s*',TempName));
    else
        dos(sprintf('del %s*',TempName));
    end
    error('Error %d -- LATEX:\n%s',hdos ,wdos);
    return;
end

%-------------------------------------------------------
% DviPs Command
%-------------------------------------------------------
Str=sprintf('dvips %s2.dvi -o %s2.ps',TempName,TempName);
disp(sprintf('[DVIPS Command] %s',Str));
[hdos,wdos]=system(Str);
if hdos ~=0
    if isunix
        dos(sprintf('rm %s*',TempName));
    else
        dos(sprintf('del %s*',TempName));
    end
    error('Error %d -- DVIPS:\n%s', hdos , wdos);
    return;
end

%-------------------------------------------------------
% eps2eps Command (can be changed with PS2EPSI)
%-------------------------------------------------------
Str = sprintf('eps2eps -dNOCACHE %s2.ps %s.eps',TempName,FileName);
disp(sprintf('[EPS2EPS Command] %s',Str));
[hdos,wdos]=system(Str);
if hdos ~=0
    if isunix
        dos(sprintf('rm %s*',TempName));
    else
        dos(sprintf('del %s*',TempName));
    end
    error('Error %d -- PS2EPSI:\n%s',hdos ,wdos);
    return;
end

%-------------------------------------------------------
% epstopdf Command
%-------------------------------------------------------
Str = sprintf('epstopdf %s.eps',FileName);
disp(sprintf('[EPSTOPDF Command] %s',Str));
[hdos,wdos]=system(Str);
if hdos ~=0
    if isunix
        dos(sprintf('rm %s*',TempName));
        dos(sprintf('rm %s.eps',FileName));
    else
        dos(sprintf('del %s*',TempName));
        dos(sprintf('del %s.eps',FileName));
    end
    error('Error %d -- EPSTOPDF:\n%s',hdos ,wdos);
    return;
end

disp(sprintf('... OK!\nPDF file [%s.pdf] has been created in the current directory\n',FileName));
%-------------------------------------------------------
% Delete all the temporary files
%-------------------------------------------------------
if isunix
    dos(sprintf('rm %s*',TempName));
    dos(sprintf('rm %s.eps',FileName));
else
    dos(sprintf('del %s*',TempName));
    dos(sprintf('del %s.eps',FileName));
end

set(0,'defaulttextinterpreter',DTI);

return;

author: Sable

repo: mcbench-benchmarks


text:
function matVer = getMLVersion
%"MLVersion"
%   Returns the numeric version of Matlab currently running.  Uses the first
%   3 digits of the Version field returned from ver('matlab').
%
%JRA 11/15/04
%
%Usage:
%   version = MLVersion
%
% 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/>.

verOutput = ver('MATLAB');
if isempty(verOutput)
    verOutput(1).Version = '';
end
verString = verOutput.Version;
matVer = str2num(verString);

% dotInd = strfind(verString, '.');
% startInd = [1 dotInd+1];
% endInd = [dotInd - 1 length(verString)];
% for i = 1:length(startInd)
%     versionCell{i} = str2num(verString(startInd(i):endInd(i)));
% end
%
% if length(versionCell) == 2
%     matVer = str2num([num2str(versionCell{1}),num2str(versionCell{2})]);
% elseif length(versionCell) == 3
%     matVer = str2num([num2str(versionCell{1}),num2str(versionCell{2}),'.',num2str(versionCell{3})]);
% end

% verOutput = ver('MATLAB');
% verString = verOutput.Version;
%
% matVer = str2num(verString);

% dotInd = strfind(verString, '.');
% startInd = [1 dotInd+1];
% endInd = [dotInd - 1 length(verString)];
%
% for i = 1:length(startInd)
%     versionCell{i} = str2num(verString(startInd(i):endInd(i)));
% end
%
% if length(versionCell) > 1
%     version = versionCell{1} + versionCell{2}*1/(10^length(versionCell{2}));
% else
%     version = versionCell{1};
% end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
classdef DesignVariableAndMeshSaverFromRes < handle

    properties (Access = private)
        levelSet
        backgroundMesh
        iterationsNumber
    end

    properties (Access = private)
        testPath
        testName
    end

    methods (Access = public)

        function obj = DesignVariableAndMeshSaverFromRes(cParams)
            obj.init(cParams)
        end

        function save(obj)
            obj.readIterNumbers()
            for iter = 1:length(obj.iterationsNumber)
                obj.obtainLevelSetAndMesh(iter);
                obj.saveLevelSetAndMesh(iter);
            end
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.testPath  = cParams.testPath;
            obj.testName  = cParams.testName;
        end

        function readIterNumbers(obj)
            files = dir(fullfile(obj.testPath,'*.res'));
            nFiles = numel(files);
            iterN  = zeros(nFiles,1);
            for ifiles = 1:nFiles
                file = files(ifiles);
                iter = obj.getIter(file);
                iterN(ifiles) = iter;
            end
            obj.iterationsNumber = sort(iterN);
        end

        function obtainLevelSetAndMesh(obj,iter)
            it = obj.iterationsNumber(iter);
            s.fileName = [obj.testName,num2str(it)];
            s.folderPath = fullfile(obj.testPath);
            wM = WrapperMshResFiles(s);
            wM.compute();
            obj.levelSet       = wM.dataRes.DesignVar1;
            obj.backgroundMesh = wM.mesh;
        end

        function saveLevelSetAndMesh(obj,iter)
            x = obj.levelSet;
            mesh = obj.backgroundMesh;
            outputName = ['DesignVariable',num2str(iter)];
            outputFile = fullfile(obj.testPath,outputName);
            save(outputFile,'x','mesh');
        end

    end

    methods (Access = private, Static)

        function iter = getIter(file)
            fName = file.name;
            iterCell = regexp(fName,'[0-9]','match');
            iterM    = cell2mat(iterCell);
            iter     = str2double(iterM);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function iter = newIterator(max_its,varargin)
%function iter = newIterator(max_its,varargin)
%
% create new iterator.  assumes the value is always increasing
    iter.max_its = max_its;
    [iter.thresh iter.min_its iter.debug] = process_options(varargin,'thresh',0,'min-its',0,'debug',false);
    iter.prev = -inf;
    iter.done = false;
    iter.its = 0;
end

author: jacobeisenstein

repo: SAGE


text:
function sClusters = in_clusters(ClusterFile, FileFormat)
% IN_CLUSTERS: Read clusters of channels from a file

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

% Read file
switch (FileFormat)
    case 'BST'
        ClusterMat = load(ClusterFile);
        if ~isfield(ClusterMat, 'Clusters')
            error('Invalid Brainstorm clusters file: missing field "Clusters".');
        end
    case 'MNE'
        sMontage = in_montage_mne(FileName);
        for i = 1:length(sMontage)
            ClusterMat.Clusters(i).Label   = sMontage(i).Name;
            ClusterMat.Clusters(i).Sensors = sMontage(i).ChanNames;
        end
    otherwise
        error('Invalid file format.');
end
if isempty(ClusterMat.Clusters)
    error(['No clusters available in file: ' FileName]);
end

% Create standardized Brainstorm structure
sClusters = repmat(db_template('cluster'), 1, length(ClusterMat.Clusters));
% Loop on all the new clusters
for i = 1:length(ClusterMat.Clusters)
    sClusters(i).Sensors = ClusterMat.Clusters(i).Sensors;
    if isfield(ClusterMat.Clusters(i), 'Label') && ~isempty(ClusterMat.Clusters(i).Label)
        sClusters(i).Label = ClusterMat.Clusters(i).Label;
    end
    if isfield(ClusterMat.Clusters(i), 'Color') && ~isempty(ClusterMat.Clusters(i).Color)
        sClusters(i).Color = ClusterMat.Clusters(i).Color;
    end
    if isfield(ClusterMat.Clusters(i), 'Function') && ~isempty(ClusterMat.Clusters(i).Function)
        sClusters(i).Function = ClusterMat.Clusters(i).Function;
    end
end


text:
classdef EnergyComputerForVoigt < EnergyComputer

    properties
    end

    methods (Access = public)

        function obj = EnergyComputerForVoigt(C,s)
            obj.init(C,s)
            obj.computeEnergy()
        end
    end

    methods (Access = protected)

        function computeEnergy(obj)
            C = obj.fourthOrder.getValue();
            s = obj.secondOrder.getValue();
            d = obj.secondOrder.getVoigtDimension();
            en = 0;
            for i = 1:d
                for j = 1:d
                    si  = s(i);
                    Cij = C(i,j);
                    sj = s(j);
                    e = si*Cij*sj;
                    en = en + e;
                end
            end
            obj.energy = en;
        end

    end
end


text:
function [bcx,bcy] = specific_flow(xbd,ybd)
%forwardstep_flow   Forward facing step inflow condition
%   [bcx,bcy] = specific_flow(xbd,ybd);
%   input
%          xbd          x coordinate vector
%          ybd          y coordinate vector
%
%   specifies forward step flow boundary condition
%   IFISS function: DJS; 6 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
bcx=0*xbd; bcy=0*xbd;
k=find(xbd==5); bcx(k)=-0.5*(1+ybd(k)).*(1-ybd(k));
return

author: johannesgerer

repo: jburkardt-m


text:
function [ score ] = findNearbyObject( all_xyz, obj_xyz )
%FINDNEARBYOBJECT Summary of this function goes here
%   Detailed explanation goes here
N = size(all_xyz,1);

lowxyz = max(all_xyz(:,1:3), repmat(obj_xyz(1:3),N,1));
topxyz = min(all_xyz(:,4:6), repmat(obj_xyz(4:6),N,1));
inter = prod(max(topxyz-lowxyz, 0.01), 2);
area1 = prod(max(all_xyz(:,4:6)-all_xyz(:,1:3), 0.01), 2);
area2 = repmat(prod(max(obj_xyz(4:6)-obj_xyz(1:3), 0.01), 2), N, 1);

score = inter./(area1+area2-inter);
valid = any(topxyz-lowxyz<=0,2);
score(valid) = 0;

end


text:
function tt() % autogenerated function wrapper
 % turned into function by Celso G Reyes 2017

ZG=ZmapGlobal.Data; % used by get_zmap_globals
report_this_filefun(mfilename('fullpath'));

catSave =...
    [ 'zmap_message_center.set_info(''Save Data'',''  '');think;',...
    '[file1,path1] = uigetfile(fullfile(ZmapGlobal.Data.data_dir, ''*.mat''), ''Earthquake Datafile'');end;',...
    'if length(file1) > 1 , sapa2 = [''save '' path1 file1 '' a faults main mainfault coastline  ''],',...
    'eval(sapa2) ,end, done'];

end

author: CelsoReyes


function Y = vl_nnloss(X,c,dzdy)
% VL_NNLOSS  CNN log-loss
%    Y = VL_NNLOSS(X, C) applies the the logistic loss to the data
%    X. X has dimension H x W x D x N, packing N arrays of W x H
%    D-dimensional vectors.
%
%    C contains the class labels, which should be integers in the range
%    1 to D. C can be an array with either N elements or with H x W x
%    1 x N dimensions. In the fist case, a given class label is
%    applied at all spatial locations; in the second case, different
%    class labels can be specified for different locations.
%
%    D can be thought of as the number of possible classes and the
%    function computes the softmax along the D dimension. Often W=H=1,
%    but this is not a requirement, as the operator is applied
%    convolutionally at all spatial locations.
%
%    DZDX = VL_NNLOSS(X, C, DZDY) computes the derivative DZDX of the
%    CNN with respect to the input X given the derivative DZDY with
%    respect to the block output Y. DZDX has the same dimension as X.

% Copyright (C) 2014 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

% no division by zero
X = X + 1e-4 ;
sz = [size(X,1) size(X,2) size(X,3) size(X,4)] ;

% index from 0
c = c - 1 ;

if numel(c) == sz(4)
  % one label per image
  c = reshape(c, [1 1 1 sz(4)]) ;
  c = repmat(c, [sz(1) sz(2)]) ;
else
  % one label per spatial location
  sz_ = [size(c,1) size(c,2) size(c,3) size(c,4)] ;
  assert(isequal(sz_, [sz(1) sz(2) 1 sz(4)])) ;
end

% convert to indeces
c_ = 0:numel(c)-1 ;
c_ = 1 + ...
  mod(c_, sz(1)*sz(2)) + ...
  (sz(1)*sz(2)) * c(:)' + ...
  (sz(1)*sz(2)*sz(3)) * floor(c_/(sz(1)*sz(2))) ;

n = sz(1)*sz(2) ;
if nargin <= 2
  Y = - sum(log(X(c_))) / n ;
else
  Y_ = - (1./X) * (dzdy/n) ;
  Y = Y_*0 ;
  Y(c_) = Y_(c_) ;
end


text:
function qwv_2d_test ( )

%*****************************************************************************80
%
%% QWV_2D_TEST tests the QWV_2D library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 May 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'QWV_2D_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the QWV_2D library.\n' );

  qwv_2d_test01 ( );
  qwv_2d_test02 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'QWV_2D_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
% 近似计算 算法稳定性的演示
clear;

% 解法一   n从小到大
S0 = 0.182;
S(1) = 1 - 5*S0;
for n = 2 : 8
    S(n) = 1/n - 5*S(n-1);
    S(n) = vpa(S(n),3);
end
disp('解法一：')
disp(S)

% 解法二   n从大到小
S(8) = 0.0204;
for n = 8: -1 : 2
    S(n-1) = 1/(5*n) - S(n)/5;
    S(n-1) = vpa(S(n-1),3);
end
disp('解法二：')
disp(S)

% 精确值
syms x
for n = 1 : 5 % 8
    S(n) = int(x^n/(x+5),0,1);
    S(n) = vpa(S(n),3);
end
disp('精确值：')
disp(S)

repo: NumericalAnalysis