text:
function [VO,M] = pm_segment(VF,PG,flags)
% Segment an MR image into Gray, White & CSF.
%
% FORMAT VO = pm_segment(PF,PG,flags)
% PF    - name(s) of image(s) to segment (must have same dimensions).
% PG    - name(s) of template image(s) for realignment.
%       - or a 4x4 transformation matrix which maps from the image to
%         the set of templates.
% flags - a structure normally based on defaults.segment
% VO    - optional output volume
% M     - affine transformation between template and image to segment
%
%                      The algorithm is four step:
%
% 1) Determine the affine transform which best matches the image with a
%    template image. If the name of more than one image is passed, then
%    the first image is used in this step. This step is not performed if
%    no template images are specified.
%
% 2) Perform Cluster Analysis with a modified Mixture Model and a-priori
%    information about the likelihoods of each voxel being one of a
%    number of different tissue types. If more than one image is passed,
%    then they they are all assumed to be in register, and the voxel
%    values are fitted to multi-normal distributions.
%
% 3) Perform morphometric operations on the grey and white partitions
%    in order to more accurately identify brain tissue. This is then used
%    to clean up the grey and white matter segments.
%
% 4) If no or 2 output arguments is/are specified, then the segmented
%    images are written to disk. The names of these images have "c1",
%    "c2" & "c3" appended to the name of the first image passed. The
%    'brainmask' is also created with "BrMsk_" as an appendix.
%
%_______________________________________________________________________
% Refs:
%
% Ashburner J & Friston KJ (1997) Multimodal Image Coregistration and
% Partitioning - a Unified Framework. NeuroImage 6:209-217
%
%_______________________________________________________________________
%
% The template image, and a-priori likelihood images are modified
% versions of those kindly supplied by Alan Evans, MNI, Canada
% (ICBM, NIH P-20 project, Principal Investigator John Mazziotta).
%_______________________________________________________________________
%
% This is a renamed version of the original spm_segment which has been
% removed from the main spm distribution, but copied into the FieldMap
% toolbox where it is still used.
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: pm_segment.m 4873 2012-08-30 19:06:26Z john $

% Create some suitable default values
%-----------------------------------------------------------------------

dr = fileparts(mfilename('fullpath'));
def_flags.estimate.priors = char(...
        fullfile(dr,'grey.nii'),...
        fullfile(dr,'white.nii'),...
        fullfile(dr,'csf.nii'));
def_flags.estimate.reg    = 0.01;
def_flags.estimate.cutoff = 30;
def_flags.estimate.samp   = 3;
def_flags.estimate.bb     =  [[-88 88]' [-122 86]' [-60 95]'];
def_flags.estimate.affreg.smosrc = 8;
def_flags.estimate.affreg.regtype = 'mni';
def_flags.estimate.affreg.weight = '';
def_flags.write.cleanup   = 1;
def_flags.write.wrt_cor   = 1;
def_flags.write.wrt_brV   = 0;
def_flags.graphics        = 1;

if nargin<3, flags = def_flags; end;
if ~isfield(flags,'estimate'),        flags.estimate        = def_flags.estimate;        end;
if ~isfield(flags.estimate,'priors'), flags.estimate.priors = def_flags.estimate.priors; end;
if ~isfield(flags.estimate,'reg'),    flags.estimate.reg    = def_flags.estimate.reg;    end;
if ~isfield(flags.estimate,'cutoff'), flags.estimate.cutoff = def_flags.estimate.cutoff; end;
if ~isfield(flags.estimate,'samp'),   flags.estimate.samp   = def_flags.estimate.samp;   end;
if ~isfield(flags.estimate,'bb'),     flags.estimate.bb     = def_flags.estimate.bb;     end;
if ~isfield(flags.estimate,'affreg'), flags.estimate.affreg = def_flags.estimate.affreg; end;
if ~isfield(flags.estimate.affreg,'smosrc'),
    flags.estimate.affreg.smosrc = def_flags.estimate.affreg.smosrc;
end;
if ~isfield(flags.estimate.affreg,'regtype'),
    flags.estimate.affreg.regtype = def_flags.estimate.affreg.regtype;
end;
if ~isfield(flags.estimate.affreg,'weight'),
    flags.estimate.affreg.weight = def_flags.estimate.affreg.weight;
end;
if ~isfield(flags,'write'),         flags.write         = def_flags.write;         end;
if ~isfield(flags.write,'cleanup'), flags.write.cleanup = def_flags.write.cleanup; end;
if ~isfield(flags.write,'wrt_cor'), flags.write.wrt_cor = def_flags.write.wrt_cor; end;
if ~isfield(flags.write,'wrt_brV'), flags.write.wrt_brV = def_flags.write.wrt_brV; end;
if ~isfield(flags,'graphics'),      flags.graphics      = def_flags.graphics;      end;

%-----------------------------------------------------------------------

if ischar(VF), VF= spm_vol(VF); end;

SP         = init_sp(flags.estimate,VF,PG);
[x1,x2,x3] = get_sampling(SP.MM,VF,flags.estimate.samp,flags.estimate.bb);
BP         = init_bp(VF, flags.estimate.cutoff, flags.estimate.reg);
CP         = init_cp(VF,x3);
sums       = zeros(8,1);

for pp=1:length(x3),
    [raw,msk] = get_raw(VF,x1,x2,x3(pp));
    s         = get_sp(SP,x1,x2,x3(pp));
    CP        = update_cp_est(CP,s,raw,msk,pp);
    sums      = sums + reshape(sum(sum(s,1),2),8,1);
end;
sums = sums/sum(sums);
CP   = shake_cp(CP);

[CP,BP,SP] = run_segment(CP,BP,SP,VF,sums,x1,x2,x3);

%save segmentation_results.mat CP BP SP VF sums

[g,w,c] = get_gwc(VF,BP,SP,CP,sums,flags.write.wrt_cor);
if flags.write.cleanup, [g,w,c,b] = clean_gwc(g,w,c); end;

% Create the segmented images + the brain mask.
%-----------------------------------------------------------------------
%offs  = cumsum(repmat(prod(VF(1).dim(1:2)),1,VF(1).dim(3)))-prod(VF(1).dim(1:2));
%pinfo = [repmat([1/255 0]',1,VF(1).dim(3)) ; offs];
[pth,nm,xt] = fileparts(deblank(VF(1).fname));

if flags.write.wrt_brV
    Nwrt = 4;
else
    Nwrt = 3;
end
for j=1:Nwrt,
    tmp   = fullfile(pth,['c', num2str(j), nm, xt]);
    if j==4, tmp = fullfile(pth,['BrMsk_', nm  xt]); end

    VO(j) = struct(...
        'fname',tmp,...
        'dim',    VF(1).dim(1:3),...
        'dt',     [spm_type('uint8'), spm_platform('bigend')],...
        'mat',    VF(1).mat,...
        'pinfo',  [1/255 0 0]',...
        'descrip','Segmented image');
end;

if nargout==0 || nargout==2,
    VO = spm_create_vol(VO);

    spm_progress_bar('Init',VF(1).dim(3),'Writing Segmented','planes completed');
    for pp=1:VF(1).dim(3),
        VO(1) = spm_write_plane(VO(1),double(g(:,:,pp))/255,pp);
        VO(2) = spm_write_plane(VO(2),double(w(:,:,pp))/255,pp);
        VO(3) = spm_write_plane(VO(3),double(c(:,:,pp))/255,pp);
        if flags.write.wrt_brV
            VO(4) = spm_write_plane(VO(4),double(b(:,:,pp))/255,pp);
        end
        spm_progress_bar('Set',pp);
    end;
    spm_progress_bar('Clear');
end;

VO(1).dat = g; VO(1).pinfo = VO(1).pinfo(1:2,:);
VO(2).dat = w; VO(2).pinfo = VO(2).pinfo(1:2,:);
VO(3).dat = c; VO(3).pinfo = VO(3).pinfo(1:2,:);
if flags.write.wrt_brV
    VO(4).dat = b; VO(4).pinfo = VO(4).pinfo(1:2,:);
end
if nargout==2
    M = SP.MM/VF(1).mat;
end

if flags.graphics, display_graphics(VF,VO,CP.mn,CP.cv,CP.mg); end;

return;
%=======================================================================

%=======================================================================
function [y1,y2,y3] = affine_transform(x1,x2,x3,M)
    y1 = M(1,1)*x1 + M(1,2)*x2 + M(1,3)*x3 + M(1,4);
    y2 = M(2,1)*x1 + M(2,2)*x2 + M(2,3)*x3 + M(2,4);
    y3 = M(3,1)*x1 + M(3,2)*x2 + M(3,3)*x3 + M(3,4);
return;
%=======================================================================

%=======================================================================
function display_graphics(VF,VS,mn,cv,mg)
% Do the graphics
nb = 3;
spm_figure('Clear','Graphics');
fg = spm_figure('FindWin','Graphics');
if ~isempty(fg),
    % Show some text
    %-----------------------------------------------------------------------
    ax = axes('Position',[0.05 0.8 0.9 0.2],'Visible','off','Parent',fg);
    text(0.5,0.80, 'Segmentation','FontSize',16,'FontWeight','Bold',...
        'HorizontalAlignment','center','Parent',ax);

    text(0,0.65, ['Image:  ' spm_file(VF(1).fname,'short50')],...
        'FontSize',14,'FontWeight','Bold','Parent',ax);

    text(0,0.40, 'Means:','FontSize',12,'FontWeight','Bold','Parent',ax);
    text(0,0.30, 'Std devs:' ,'FontSize',12,'FontWeight','Bold','Parent',ax);
    text(0,0.20, 'N vox:','FontSize',12,'FontWeight','Bold','Parent',ax);
    for j=1:nb,
        text((j+0.5)/(nb+1),0.40, num2str(mn(1,j)),...
            'FontSize',12,'FontWeight','Bold',...
            'HorizontalAlignment','center','Parent',ax);
        text((j+0.5)/(nb+1),0.30, num2str(sqrt(cv(1,1,j))),...
            'FontSize',12,'FontWeight','Bold',...
            'HorizontalAlignment','center','Parent',ax);
        text((j+0.5)/(nb+1),0.20, num2str(mg(1,j)/sum(mg(1,:))),...
            'FontSize',12,'FontWeight','Bold',...
            'HorizontalAlignment','center','Parent',ax);
    end;
    if length(VF) > 1,
        text(0,0.10,...
        'Note: only means and variances for the first image are shown',...
        'Parent',ax,'FontSize',12);
    end;

    M1 = VS(1).mat;
    M2 = VF(1).mat;
    for i=1:5,
        M   = spm_matrix([0 0 i*VF(1).dim(3)/6]);
        img = spm_slice_vol(VF(1),M,VF(1).dim(1:2),1);
        img(1,1) = eps;
        ax = axes('Position',...
            [0.05 0.75*(1-i/5)+0.05 0.9/(nb+1) 0.75/5],...
            'Visible','off','Parent',fg);
        imagesc(rot90(img), 'Parent', ax);
        set(ax,'Visible','off','DataAspectRatio',[1 1 1]);

        for j=1:3,
            img = spm_slice_vol(VS(j),M2\M1*M,VF(1).dim(1:2),1);
            ax  = axes('Position',...
                [0.05+j*0.9/(nb+1) 0.75*(1-i/5)+0.05 0.9/(nb+1) 0.75/5],...
                'Visible','off','Parent',fg);
            image(rot90(img*64), 'Parent', ax);
            set(ax,'Visible','off','DataAspectRatio',[1 1 1]);
        end;
    end;

    spm_print;
    drawnow;
end;
return;
%=======================================================================

%=======================================================================
function M = get_affine_mapping(VF,VG,aflags)

if ~isempty(VG) && ischar(VG), VG = spm_vol(VG); end;

if ~isempty(VG) && isstruct(VG),
    % Affine registration so that a priori images match the image to
    % be segmented.
    %-----------------------------------------------------------------------

    VFS = spm_smoothto8bit(VF(1),aflags.smosrc);

    % Scale all images approximately equally
    % ---------------------------------------------------------------
    for i=1:length(VG),
        VG(i).pinfo(1:2,:) = VG(i).pinfo(1:2,:)/spm_global(VG(i));
    end;
    VFS(1).pinfo(1:2,:) = VFS(1).pinfo(1:2,:)/spm_global(VFS(1));

    spm_plot_convergence('Init','Affine Registration','Mean squared difference','Iteration');
    flags     = struct('sep',aflags.smosrc, 'regtype',aflags.regtype,'WG',[],'globnorm',0,'debug',0);
    M         = eye(4);
    [M,scal]  = spm_affreg(VG, VFS, flags, M);

    if ~isempty(aflags.weight), flags.WG = spm_vol(aflags.weight); end;

    flags.sep = aflags.smosrc/2;
    M         = spm_affreg(VG, VFS, flags, M,scal);
    spm_plot_convergence('Clear');

elseif all(size(VG) == [4 4])
    % Assume that second argument is a matrix that will do the job
    %-----------------------------------------------------------------------
    M = VG;
else
    % Assume that image is normalized
    %-----------------------------------------------------------------------
    M = eye(4);
end
return;
%=======================================================================

%=======================================================================
function [x1,x2,x3] = get_sampling(MM,VF,samp,bb1)
% Voxels to sample during the cluster analysis
%-----------------------------------------------------------------------

% A bounding box for the brain in Talairach space.
%bb = [ [-88 88]' [-122 86]' [-60 95]'];
%c = [bb(1,1) bb(1,2) bb(1,3) 1
%     bb(1,1) bb(1,2) bb(2,3) 1
%     bb(1,1) bb(2,2) bb(1,3) 1
%     bb(1,1) bb(2,2) bb(2,3) 1
%     bb(2,1) bb(1,2) bb(1,3) 1
%     bb(2,1) bb(1,2) bb(2,3) 1
%     bb(2,1) bb(2,2) bb(1,3) 1
%     bb(2,1) bb(2,2) bb(2,3) 1]';
%tc = MM\c;
%tc = tc(1:3,:)';
%mx = max(tc);
%mn = min(tc);
%bb = [mn ; mx];
%vx = sqrt(sum(VF(1).mat(1:3,1:3).^2));
%samp = round(max(abs([4 4 4]./vx), [1 1 1]));
%x1 = bb(1,1):samp(1):bb(2,1);
%x2 = bb(1,2):samp(2):bb(2,2);
%x3 = bb(1,3):samp(3):bb(2,3);
%return;

% A bounding box for the brain in Talairach space.
if nargin<4, bb1 = [ [-88 88]' [-122 86]' [-60 95]']; end;

% A mapping from a unit radius sphere to a hyper-ellipse
% that is just enclosed by the bounding box in Talairach
% space.
M0 = [diag(diff(bb1)/2) mean(bb1)';[0 0 0 1]];

% The mapping from voxels to Talairach space is MM,
% so the ellipse in the space of the image becomes:
M0 = MM\M0;

% So to work out the bounding box in the space of the
% image that just encloses the hyper-ellipse.
tmp = M0(1:3,1:3);
tmp = diag(tmp*tmp'/diag(sqrt(diag(tmp*tmp'))));
bb  = round([M0(1:3,4)-tmp M0(1:3,4)+tmp])';
bb  = min(max(bb,[1 1 1 ; 1 1 1]),[VF(1).dim(1:3) ; VF(1).dim(1:3)]);

% Want to sample about every 3mm
tmp  = sqrt(sum(VF(1).mat(1:3,1:3).^2))';
samp = round(max(abs(tmp.^(-1)*samp), [1 1 1]'));

x1 = bb(1,1):samp(1):bb(2,1);
x2 = bb(1,2):samp(2):bb(2,2);
x3 = bb(1,3):samp(3):bb(2,3);

return;
%=======================================================================

%=======================================================================
function [CP,BP,SP] = run_segment(CP,BP,SP,VF,sums,x1,x2,x3)
oll = -Inf;
spm_plot_convergence('Init','Segmenting','Log-likelihood','Iteration #');

for iter = 1:64,
    ll= 0;
    for pp = 1:length(x3), % Loop over planes
        bf        = get_bp(BP,x1,x2,x3(pp));
        [raw,msk] = get_raw(VF,x1,x2,x3(pp));
        s         = get_sp(SP,x1,x2,x3(pp));
        cor       = bf.*raw;
        [P,ll0]   = get_p(cor,msk,s,sums,CP,bf);
        ll        = ll + ll0;
        CP        = update_cp_est(CP,P,cor,msk,pp);
        BP        = update_bp_est(BP,P,cor,CP,msk,x1,x2,x3(pp));
    end;

    BP = update_bp(BP);
    if iter>1, spm_plot_convergence('Set',ll); end;
    %fprintf('\t%g\n', ll);

    % Stopping criterion
    %-----------------------------------------------------------------------
    if iter == 2,
        ll2 = ll;
    elseif iter > 2 && abs((ll-oll)/(ll-ll2)) < 0.0001
        break;
    end;
    oll = ll;
end;
spm_plot_convergence('Clear');
return;
%=======================================================================

%=======================================================================
function BP = init_bp(VF,co,reg)
m        = length(VF);
tmp      = sqrt(sum(VF(1).mat(1:3,1:3).^2));
BP.nbas  = max(round((VF(1).dim(1:3).*tmp)/co),[1 1 1]);
BP.B1    = spm_dctmtx(VF(1).dim(1),BP.nbas(1));
BP.B2    = spm_dctmtx(VF(1).dim(2),BP.nbas(2));
BP.B3    = spm_dctmtx(VF(1).dim(3),BP.nbas(3));

nbas     = BP.nbas;
if prod(BP.nbas)>1,
    % Set up a priori covariance matrix
    vx = sqrt(sum(VF(1).mat(1:3,1:3).^2));
    kx=(pi*((1:nbas(1))'-1)*pi/vx(1)/VF(1).dim(1)*10).^2;
    ky=(pi*((1:nbas(2))'-1)*pi/vx(2)/VF(1).dim(2)*10).^2;
    kz=(pi*((1:nbas(3))'-1)*pi/vx(3)/VF(1).dim(3)*10).^2;

    % Cost function based on sum of squares of 4th derivatives
    IC0 =  (1*kron(kz.^4,kron(ky.^0,kx.^0)) +...
            1*kron(kz.^0,kron(ky.^4,kx.^0)) +...
            1*kron(kz.^0,kron(ky.^0,kx.^4)) +...
            4*kron(kz.^3,kron(ky.^1,kx.^0)) +...
            4*kron(kz.^3,kron(ky.^0,kx.^1)) +...
            4*kron(kz.^1,kron(ky.^3,kx.^0)) +...
            4*kron(kz.^0,kron(ky.^3,kx.^1)) +...
            4*kron(kz.^1,kron(ky.^0,kx.^3)) +...
            4*kron(kz.^0,kron(ky.^1,kx.^3)) +...
            6*kron(kz.^2,kron(ky.^2,kx.^0)) +...
            6*kron(kz.^2,kron(ky.^0,kx.^2)) +...
            6*kron(kz.^0,kron(ky.^2,kx.^2)) +...
           12*kron(kz.^2,kron(ky.^1,kx.^1)) +...
           12*kron(kz.^1,kron(ky.^2,kx.^1)) +...
           12*kron(kz.^1,kron(ky.^1,kx.^2)) )*reg;

    %IC0(1) = max(IC0);
    BP.IC0 = diag(IC0(2:end));

    % Initial estimate for intensity modulation field
    BP.T   = zeros(nbas(1),nbas(2),nbas(3),length(VF));
    %-----------------------------------------------------------------------
else
    BP.T   = zeros([1 1 1 length(VF)]);
    BP.IC0 = [];
end;
BP.Alpha = zeros(prod(BP.nbas(1:3)),prod(BP.nbas(1:3)),m);
BP.Beta  = zeros(prod(BP.nbas(1:3)),m);
return;
%=======================================================================

%=======================================================================
function BP = update_bp_est(BP,p,cor,CP,msk,x1,x2,x3)
if prod(BP.nbas)<=1, return; end;
B1 = BP.B1(x1,:);
B2 = BP.B2(x2,:);
B3 = BP.B3(x3,:);
for j=1:size(BP.Alpha,3),
    cr  = cor(:,:,j);
    w1 = zeros(size(cr));
    w2 = zeros(size(cr));
    for i=[1 2 3 4 5 6 7 8],
        tmp = p(:,:,i)*CP.cv(j,j,i)^(-1);
        w1  = w1 + tmp.*(CP.mn(j,i) - cr);
        w2  = w2 + tmp;
    end;
    wt1       = 1 + cr.*w1;
    wt2       = cr.*(cr.*w2 - w1);
    wt1(~msk) = 0;
    wt2(~msk) = 0;

    BP.Beta(:,j)    = BP.Beta(:,j)    + kron(B3',spm_krutil(wt1,B1,B2,0));
    BP.Alpha(:,:,j) = BP.Alpha(:,:,j) + kron(B3'*B3,spm_krutil(wt2,B1,B2,1));
end;
return;
%=======================================================================

%=======================================================================
function BP = update_bp(BP)
if prod(BP.nbas)<=1, return; end;
for j=1:size(BP.Alpha,3),
    x     = BP.T(:,:,:,j);
    x     = x(:);
    x     = x(2:end);
    Alpha = BP.Alpha(2:end,2:end,j);
    Beta  = BP.Beta(2:end,j);
    x     = (Alpha + BP.IC0)\(Alpha*x + Beta);

    BP.T(:,:,:,j) = reshape([0 ; x],BP.nbas(1:3));
    BP.Alpha      = zeros(size(BP.Alpha));
    BP.Beta       = zeros(size(BP.Beta));
end;
return;
%=======================================================================

%=======================================================================
function bf = get_bp(BP,x1,x2,x3)
bf = ones(length(x1),length(x2),size(BP.Alpha,3));
if prod(BP.nbas)<=1, return; end;
B1 = BP.B1(x1,:);
B2 = BP.B2(x2,:);
B3 = BP.B3(x3,:);
for i=1:size(BP.Alpha,3),
    t = reshape(reshape(BP.T(:,:,:,i),...
        BP.nbas(1)*BP.nbas(2),BP.nbas(3))*B3', BP.nbas(1), BP.nbas(2));
    bf(:,:,i) = exp(B1*t*B2');
end;
return;
%=======================================================================

%=======================================================================
function [dat,msk] = get_raw(VF,x1,x2,x3)
[X1,X2,X3] = ndgrid(x1,x2,x3);
for i=1:length(VF),
    [Y1,Y2,Y3] = affine_transform(X1,X2,X3,VF(i).mat\VF(1).mat);
    dat(:,:,i) = spm_sample_vol(VF(i),Y1,Y2,Y3,1);
end;
msk = all(dat,3) & all(isfinite(double(dat)),3);
return;
%=======================================================================

%=======================================================================
function CP = init_cp(VF,x3)
n = 8;
m = length(VF);
p = length(x3);
CP.mom0 = zeros(1,n,p)+eps;
CP.mom1 = zeros(m,n,p);
CP.mom2 = zeros(m,m,n,p)+eps;

% Occasionally the dynamic range of the images is such that many voxels
% all have the same intensity.  Adding cv0 is an attempt to improve the
% stability of the algorithm if this occurs. The value 0.083 was obtained
% from var(rand(1000000,1)).  It prbably isn't the best way of doing
% things, but it appears to work.
CP.cv0 = zeros(m,m);
for i=1:m,
    if spm_type(VF(i).dt(1),'intt'),
        CP.cv0(i,i)=0.083*mean(VF(i).pinfo(1,:));
    end;
end;
return;
%=======================================================================

%=======================================================================
function CP = shake_cp(CP)
CP.mom0(:,5,:)   = CP.mom0(:,1,:);
CP.mom0(:,6,:)   = CP.mom0(:,2,:);
CP.mom0(:,7,:)   = CP.mom0(:,3,:);
CP.mom1(:,5,:)   = CP.mom1(:,1,:);
CP.mom1(:,6,:)   = CP.mom1(:,2,:);
CP.mom1(:,7,:)   = CP.mom1(:,3,:);
CP.mom1(:,8,:)   = 0;
CP.mom2(:,:,5,:) = CP.mom2(:,:,1,:);
CP.mom2(:,:,6,:) = CP.mom2(:,:,2,:);
CP.mom2(:,:,7,:) = CP.mom2(:,:,3,:);
return;
%=======================================================================

%=======================================================================
function CP = update_cp_est(CP,P,dat,msk,p)
m   = size(dat,3);
d   = size(P);
P   = reshape(P,[d(1)*d(2),d(3)]);
dat = reshape(dat,[d(1)*d(2),m]);
P(~msk(:),:)   = [];
dat(~msk(:),:) = [];
for i=1:size(CP.mom0,2),
    CP.mom0(1,i,p)   = sum(P(:,i));
    CP.mom1(:,i,p)   = sum((P(:,i)*ones(1,m)).*dat)';
    CP.mom2(:,:,i,p) = ((P(:,i)*ones(1,m)).*dat)'*dat;
end;

for i=1:size(CP.mom0,2),
    CP.mg(1,i)   = sum(CP.mom0(1,i,:),3);
    CP.mn(:,i)   = sum(CP.mom1(:,i,:),3)/CP.mg(1,i);

    tmp          = (CP.mg(1,i).*CP.mn(:,i))*CP.mn(:,i)';
    tmp          = tmp-eye(size(tmp))*eps*1e6;
    CP.cv(:,:,i) = (sum(CP.mom2(:,:,i,:),4) - tmp)/CP.mg(1,i) + CP.cv0;
end;
CP.mg   = CP.mg/sum(CP.mg);
return;
%=======================================================================

%=======================================================================
function [p,ll] = get_p(cor,msk,s,sums,CP,bf)
d   = [size(cor) 1 1];
n   = size(CP.mg,2);
cor = reshape(cor,d(1)*d(2),d(3));
cor = cor(msk,:);
p   = zeros(d(1)*d(2),n);
if ~any(msk), p  = reshape(p,d(1),d(2),n); ll=0; return; end;

for i=1:n,
    amp       = 1/sqrt((2*pi)^d(3) * det(CP.cv(:,:,i)));
    dst       = (cor-ones(size(cor,1),1)*CP.mn(:,i)')/sqrtm(CP.cv(:,:,i));
    dst       = sum(dst.*dst,2);
    tmp       = s(:,:,i);
    p(msk,i)  = (amp*CP.mg(1,i)/sums(i))*exp(-0.5*dst).*tmp(msk) +eps;
end;
sp = sum(p,2);
ll = sum(log(sp(msk).*bf(msk)+eps));
sp(~msk) = Inf;
for i=1:n, p(:,i) = p(:,i)./sp; end;
p  = reshape(p,d(1),d(2),n);
return;
%=======================================================================

%=======================================================================
function SP = init_sp(flags,VF,PG)
SP.VB       = spm_vol(flags.priors);
MM          = get_affine_mapping(VF,PG,flags.affreg);
%VF          = spm_vol(PF);
SP.MM       = MM*VF(1).mat;
SP.w        = 0.98;
return;
%=======================================================================

%=======================================================================
function s = get_sp(SP,x1,x2,x3)
[X1,X2,X3] = ndgrid(x1,x2,x3);
[Y1,Y2,Y3] = affine_transform(X1,X2,X3,SP.VB(1).mat\SP.MM);
w1  = SP.w;
w2  = (1-w1)/2;
s   = zeros([size(Y1),4]);
for i=1:3,
    s(:,:,i) = spm_sample_vol(SP.VB(i),Y1,Y2,Y3,1)*w1+w2;
end;
s(:,:,4:8)   = repmat(abs(1-sum(s(:,:,1:3),3))/5,[1 1 5]);
return;
%=======================================================================

%=======================================================================
function [g,w,c] = get_gwc(VF,BP,SP,CP,sums,wc)

if wc,
    VC = VF;
    for j=1:length(VF),
        [pth,nm,xt,vr] = spm_fileparts(deblank(VF(j).fname));
        VC(j).fname    = fullfile(pth,['m' nm xt vr]);
        VC(j).descrip  = 'Bias corrected image';
    end;
    VC = spm_create_vol(VC);
end;

spm_progress_bar('Init',VF(1).dim(3),'Creating Segmented','planes completed');
x1 = 1:VF(1).dim(1);
x2 = 1:VF(1).dim(2);
x3 = 1:VF(1).dim(3);

g = uint8(0); g(VF(1).dim(1),VF(1).dim(2),VF(1).dim(3)) = 0;
w = uint8(0); w(VF(1).dim(1),VF(1).dim(2),VF(1).dim(3)) = 0;
c = uint8(0); c(VF(1).dim(1),VF(1).dim(2),VF(1).dim(3)) = 0;

for pp=1:length(x3),
    bf        = get_bp(BP,x1,x2,x3(pp));
    [raw,msk] = get_raw(VF,x1,x2,x3(pp));
    cor       = raw.*bf;
    if wc,
        for j=1:length(VC),
            VC(j) = spm_write_plane(VC(j),cor(:,:,j),pp);
        end;
    end;
    s         = get_sp(SP,x1,x2,x3(pp));
    p         = get_p(cor,msk,s,sums,CP,bf);
    g(:,:,pp) = uint8(round(p(:,:,1)*255));
    w(:,:,pp) = uint8(round(p(:,:,2)*255));
    c(:,:,pp) = uint8(round(p(:,:,3)*255));

    spm_progress_bar('Set',pp);
end;
spm_progress_bar('Clear');

return;
%=======================================================================

%=======================================================================
function [g,w,c,b] = clean_gwc(g,w,c)
b    = w;
b(1) = w(1);

% Build a 3x3x3 seperable smoothing kernel
%-----------------------------------------------------------------------
kx=[0.75 1 0.75];
ky=[0.75 1 0.75];
kz=[0.75 1 0.75];
sm=sum(kron(kron(kz,ky),kx))^(1/3);
kx=kx/sm; ky=ky/sm; kz=kz/sm;

% Erosions and conditional dilations
%-----------------------------------------------------------------------
niter = 32;
spm_progress_bar('Init',niter,'Extracting Brain','Iterations completed');
for j=1:niter,
    if j>2, th=0.15; else th=0.6; end; % Dilate after two its of erosion.
    for i=1:size(b,3),
        gp = double(g(:,:,i));
        wp = double(w(:,:,i));
        bp = double(b(:,:,i))/255;
        bp = (bp>th).*(wp+gp);
        b(:,:,i) = uint8(round(bp));
    end;
    spm_conv_vol(b,b,kx,ky,kz,-[1 1 1]);
    spm_progress_bar('Set',j);
end;
th = 0.05;
for i=1:size(b,3),
    gp       = double(g(:,:,i))/255;
    wp       = double(w(:,:,i))/255;
    cp       = double(c(:,:,i))/255;
    bp       = double(b(:,:,i))/255;
    bp       = ((bp>th).*(wp+gp))>th;
    g(:,:,i) = uint8(round(255*gp.*bp./(gp+wp+cp+eps)));
    w(:,:,i) = uint8(round(255*wp.*bp./(gp+wp+cp+eps)));
    c(:,:,i) = uint8(round(255*(cp.*bp./(gp+wp+cp+eps)+cp.*(1-bp))));
    b(:,:,i) = uint8(round(255*bp));
end;
spm_progress_bar('Clear');
return;


text:
function r=v_rotqr2ro(q)
%ROTQR2RO converts a real quaternion to a 3x3 rotation matrix
% Inputs:
%
%     Q(4,...)      Real-valued quaternion array (possibly unnormalized)
%
% Outputs:
%
%     R(3,3,...)    Rotation matrix array
%                   Plots a diagram if no output specified
%
% In the quaternion representation of a rotation, and q(1) = cos(t/2)
% where t is the angle of rotation in the range 0 to 2pi
% and q(2:4)/sin(t/2) is a unit vector lying along the axis of rotation
% a positive rotation about [0 0 1] takes the X axis towards the Y axis.
%
%      Copyright (C) Mike Brookes 2007-2018
%      Version: $Id: v_rotqr2ro.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

persistent a b c d e f g h m
if isempty(a)
    a=[1 5 9];
    b=[11 16 6];
    c=[16 6 11];
    d=[4 8 3];
    e=[10 15 14];
    f=[4 2 3];
    g=[2 6 7];
    h=[1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4]';
    m=[1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4]';
end
sz=size(q);
q=reshape(q,4,[]);      % convert to 2D matrix
nq =size(q,2);          % number of quaternions to convert
p=2*q(h,:).*q(m,:)./repmat(sum(q.^2,1),16,1); % force normalized and calculate quadratic terms
r=zeros(9,nq);          % space for nq rotation matrices
r(a,:)=1-p(b,:)-p(c,:);
r(d,:)=p(e,:)-p(f,:);
r(g,:)=p(e,:)+p(f,:);
r=reshape(r,[3 3 sz(2:end)]);
if ~nargout
    % display rotated cube
    q1=q(:,1);
    q1=q1/sqrt(q1'*q1); % normalize the first quaternion
    clf('reset'); % clear current axis
    v0=[-1 1 1 -1 -1 1 1 -1; -1 -1 1 1 -1 -1 1 1; -1 -1 -1 -1 1 1 1 1]*0.5; % unrotated coordinates
    v=r(:,:,1)*v0; % use first elemnt of r to rotate
    fv=[4 1 5 8; 2 3 7 6; 1 2 6 5; 3 4 8 7; 2 1 4 3; 5 6 7 8]; % verices for each face
    fc=[0 1 1; 1 0 0; 1 0 1; 0 1 0; 1 1 0; 0 0 1]; % colours for faces
    xc={[25 46 46 25; 55 55 49 49]/100,
        [25 33 33 39 39 46 46 39 39 33 33 25; 55 55 63 63 55 55 49 49 42 42 49 49]/100,
        [48 56 62 69 76 66 76 69 62 56 48 59; 68 68 58 68 68 53 37 37 47 37 37 53]/100,
        [48 55 62 69 77 65 65 59 59; 68 68 55 68 68 50 37 37 50]/100,
        [50 74 74 57 74 74 50 50 67 50; 68 68 62 43 43 37 37 43 62 62]/100};
    xf=[1 3; 2 3; 1 4; 2 4; 1 5; 2 5]; % characters to plot on each face
    nf=size(fv,1); % number of faces
    for i=1:6
        p(i)=patch(v(1,fv(i,:)),v(2,fv(i,:)),v(3,fv(i,:)),fc(i,:));
        set(p(i),'FaceAlpha',0.65);
        k=1.001; % factor to move out labels slightly to get correct depth ordering
        for j=1:2
            xij=xc{xf(i,j)}; % relative coordinates of character vertices
            patch(k*(v(1,fv(i,1))+(v(1,fv(i,2))-v(1,fv(i,1)))*xij(1,:)+(v(1,fv(i,4))-v(1,fv(i,1)))*xij(2,:)), ...
                k*(v(2,fv(i,1))+(v(2,fv(i,2))-v(2,fv(i,1)))*xij(1,:)+(v(2,fv(i,4))-v(2,fv(i,1)))*xij(2,:)), ...
                k*(v(3,fv(i,1))+(v(3,fv(i,2))-v(3,fv(i,1)))*xij(1,:)+(v(3,fv(i,4))-v(3,fv(i,1)))*xij(2,:)),1-fc(i,:));
        end
    end
    qa=q1(2:4);
    qm=max(abs(qa));
    if qm>1e-6*abs(q1(1))
        qa=qa/(qm/0.7); % scale so axis extends outside the cube
        hold on;
        plot3([1 -1]*qa(1),[1 -1]*qa(2),[1 -1]*qa(3),'-');
        hold off
    end
    th=360/pi*acos(abs(q1(1)));
    xlabel('x axis');
    ylabel('y axis');
    zlabel('z axis');
    q=q(:,1)*((2*(q(find(q1(:)~=0,1))>0)-1)); % force leading coefficient to be positive
    title(sprintf('%d^\\circ, qr'' = [%.2f,%.2f,%.2f,%.2f], eu_{xyzo}'' = [%d, %d, %d]^\\circ',round(th),q1(:),round(v_rotqr2eu('xyz',q1)*180/pi)));
    axis([-1 1 -1 1 -1 1 0 1]*sqrt(3)/2);
    axis equal
    grid on
    view(3);
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function varargout = online_ssvep_gui_analysis(varargin)
% ONLINE_SSVEP_GUI_ANALYSIS MATLAB code for online_ssvep_gui_analysis.fig
%      ONLINE_SSVEP_GUI_ANALYSIS, by itself, creates a new ONLINE_SSVEP_GUI_ANALYSIS or raises the existing
%      singleton*.
%
%      H = ONLINE_SSVEP_GUI_ANALYSIS returns the handle to a new ONLINE_SSVEP_GUI_ANALYSIS or the handle to
%      the existing singleton*.
%
%      ONLINE_SSVEP_GUI_ANALYSIS('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in ONLINE_SSVEP_GUI_ANALYSIS.M with the given input arguments.
%
%      ONLINE_SSVEP_GUI_ANALYSIS('Property','Value',...) creates a new ONLINE_SSVEP_GUI_ANALYSIS or raises
%      the existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before online_ssvep_gui_analysis_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to online_ssvep_gui_analysis_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

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

% Last Modified by GUIDE v2.5 14-Aug-2017 19:12:30

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

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

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

% Choose default command line output for online_ssvep_gui_analysis
% save('InitData.mat','handles');
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);
initialize_gui(hObject, handles, false);

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

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

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

function initialize_gui(fig_handle, handles, isreset)
% If the metricdata field is present and the reset flag is false, it means
% we are we are just re-initializing a GUI by calling it from the cmd line
% while it is up. So, bail out as we dont want to reset the data.
RESET(handles, true);
% Update handles structure
guidata(handles.figure1, handles);

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

% Hints: get(hObject,'String') returns contents of freq as text
%        str2double(get(hObject,'String')) returns contents of freq as a double
input = get(handles.freq,'string');
if regexp(input, '[^0-9., ;]')
    set(handles.noti_text, 'String', sprintf('[%s] is not acceptable', input))
    set(handles.freq,'string','6.67, 7.5, 8.57,10, 12');
end

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

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

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

% Hints: get(hObject,'String') returns contents of time_stimulus as text
%        str2double(get(hObject,'String')) returns contents of time_stimulus as a double
input = get(handles.time_stimulus,'string');
if regexp(input, '[^0-9]')
    set(handles.noti_text, 'String', sprintf('[%s] is not acceptable', input))
    set(handles.timeSti,'string','5');
end

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

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

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

% Hints: get(hObject,'String') returns contents of freq_range_min as text
%        str2double(get(hObject,'String')) returns contents of freq_range_min as a double
input = get(handles.freq_range_min,'string');
if regexp(input, '[^0-9.]')
    set(handles.noti_text, 'String', sprintf('[%s] is not acceptable', input))
    set(handles.freq_range_min,'string','0.1');
end

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

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

% --- Executes on button press in btnReset.
function btnReset_Callback(hObject, eventdata, handles)
% hObject    handle to btnReset (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% disp(get(handles.btwBox,'String'));
% load('InitData.mat');
% handles.output = hObject;
% disp(get(handles.btwBox,'String'));
% guidata(handles.figure1, handles);
RESET(handles,false);

% --- Executes on button press in btnStart.
function btnStart_Callback(hObject, eventdata, handles)
% hObject    handle to btnStart (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global sock state
try
    chan = get(handles.chan,'String');
    if isequal(chan,'all')
        chan = filt_EEG_CHANNEL(state.clab);
        chan = cell2mat(chan(1,:));
    else
        chan = str2double(strsplit(chan, ' ,'));
    end
    time_stimulus = str2double(get(handles.time_stimulus, 'String'));
    freq = str2double(strsplit(get(handles.freq, 'String'),','));
    fRange = [str2double(get(handles.freq_range_min, 'String')), ...
        str2double(get(handles.freq_range_max,'String'))];
    freq_class = strsplit(get(handles.freq_class,'String'),',');
    axes_fft = handles.axes_FFT;
    axes_cca = handles.axes_CCA;
    results_text = handles.results_text;
catch
    set(handles.noti_text, 'String', 'Somethings wrong...');
    return;
end
if ~isequal(sock.status,'closed')
    s = online_ssvep_Analysis(axes_fft, axes_cca, results_text, {'chan', chan;...
        'time_stimulus', time_stimulus; 'freq', freq;...
        'fRange', fRange; 'class', freq_class});
else
    set(handles.noti_text, 'String', 'check your socket');
    return;
end
bbci_acquire_bv('close');
set(handles.noti_text, 'String', s);

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

% Hints: get(hObject,'String') returns contents of chan as text
%        str2double(get(hObject,'String')) returns contents of chan as a double
input = get(handles.chan,'string');
if ~isequal(input,'all')
    if regexp(input, '[^0-9., ;]')
        set(handles.noti_text, 'String', sprintf('[%s] is not acceptable', input))
        set(handles.freq,'string','all');
    end
end

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

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

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

try
    bbci_acquire_bv('close');
    params = struct;
    state = bbci_acquire_bv('init', params);
catch
    set(handles.noti_text,'String', 'Please open the BrainVisino Recorder');
    return;
end

if isempty(sock) || isequal(sock.status,'closed')
    port = str2double(get(handles.tcpip_port,'String'));
    sock = tcpip('localhost', port,'timeout',2);
    try
        fopen(sock);
    catch
        set(handles.noti_text, 'String', 'First open the server');
        return;
    end
    while true
        fwrite(sock,19);
        [~, ~, markerdescr, ~] = bbci_acquire_bv(state);
        if ~isempty(markerdescr)&&(markerdescr == 19)
            flushoutput(sock);
            break;
        end
    end
    set(handles.closebtn,'Visible', 'On');
    set(handles.trig_button, 'String', 'OK!');
    set(handles.noti_text, 'String', 'Connect!');
end

% --- Executes on button press in selChan.
function selChan_Callback(hObject, eventdata, handles)
% hObject    handle to selChan (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global state
filt_EEG_CHANNEL(state.clab);
select_channel_v2({'chan', state.clab});
% figure;
% uitable('Data',handles.initialization.bv.clab);

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

% Hints: get(hObject,'String') returns contents of freq_class as text
%        str2double(get(hObject,'String')) returns contents of freq_class as a double

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

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

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

% Hints: get(hObject,'String') returns contents of tcpip_port as text
%        str2double(get(hObject,'String')) returns contents of tcpip_port as a double
input = get(handles.tcpip_port,'string');
if regexp(input, '[^0-9]')
    set(handles.noti_text, 'String', sprintf('[%s] is not acceptable', input))
    set(handles.timeSti,'string','12300');
end

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

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

function RESET(handles, init)
% handles    empty - handles not created until after all CreateFcns called
% init        initialization
set(handles.chan,'String','all');
set(handles.time_stimulus, 'String','5');
set(handles.freq, 'String','6.67, 7.5, 8.57,10, 12');
set(handles.freq_range_min, 'String', '0.1');
set(handles.freq_range_max, 'String', '50');
set(handles.freq_class,'String','up, left, center, right, down');
set(handles.trig_button,'String', 'Check');
set(handles.closebtn,'Visible', 'Off');
set(handles.results_text, 'String', '');
set(handles.noti_text, 'String', 'Welcome');
set(handles.axes_FFT, 'XTickLabel','');
set(handles.axes_CCA, 'XTickLabel', '');
set(handles.axes_FFT, 'YTickLabel', '');
set(handles.axes_CCA, 'YTickLabel', '');
delete(get(handles.axes_FFT,'Children'));
delete(get(handles.axes_CCA,'Children'));
global sock
try
    bbci_acquire_bv('close');
catch
    set(handles.noti_text,'String', 'Please open the BV');
end
try
    fclose(sock);
catch
end

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

% Hints: get(hObject,'String') returns contents of freq_range_max as text
%        str2double(get(hObject,'String')) returns contents of freq_range_max as a double
input = get(handles.freq_range_max,'string');
if regexp(input, '[^0-9.]')
    set(handles.noti_text, 'String',  sprintf('[%s] is not acceptable', input))
    set(handles.timeSti,'string','50');
end

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

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

% --- Executes on button press in closebtn.
function closebtn_Callback(hObject, eventdata, handles)
% hObject    handle to closebtn (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global sock;
fclose(sock);
set(handles.noti_text, 'String', 'Close the socket');
set(handles.trig_button,'String', 'Check');
set(handles.closebtn,'Visible', 'Off');

author: PatternRecognition

repo: OpenBMI


text:
% Digital circuit sizing example (GP)
% Boyd, Kim, Patil, and Horowitz, "Digital circuit optimization
% via geometric programming"
% Written for CVX by Almir Mutapcic 02/08/06
%
% Solves the problem of choosing gate scale factors x_i to give
% minimum ckt delay, subject to limits on the total area and power.
% Uses max gate arrival time T formulation that avoids evaluation
% of the delay over all paths in the circuit.
%
%   minimize   T_bar
%       s.t.   T_j <= T_bar      for j an output gate
%              T_j + d_i <= T_i  for j in FI(i)
%              P <= Pmax, A <= Amax
%              x >= 1
%
% where variables are x and T.
%
% We use the circuit topology presented in figure 1 (page 902),
% where we take gates 1, 3 and 6 to be inverters (INV),
% gates 2 and 7 to be three input NANDs (NAND3),
% and gates 4 and 5 to be two input NORs (NOR2).

%********************************************************************
% user specified data (specify problem constant and ckt topology)
%********************************************************************
m = 7;        % number of gates
Vdd = 5;      % supply voltage
Amax = 250;   % maximum area spec

% gate specs
INV   = struct('Cin',3, 'Cint',3, 'Rdrv',0.48, 'A',3,  'Ileak',0.006);
NAND3 = struct('Cin',4, 'Cint',6, 'Rdrv',0.48, 'A',8,  'Ileak',0.007);
NOR2  = struct('Cin',5, 'Cint',6, 'Rdrv',0.48, 'A',10, 'Ileak',0.009);

clear gates;
gates([1 3 6]) = INV;
gates([2 7])   = NAND3;
gates([4 5])   = NOR2;

% primary inputs and primary outputs labels (start with m+1)
primary_inputs = [8 9 10];
primary_outputs = [11 12];
M = m + length( primary_inputs ) + length( primary_outputs );

% fan-in cell array
FI = cell(M,1);
FI{1} = 8;
FI{2} = [8 9 10];
FI{3} = 10;
FI{4} = [1 2];
FI{5} = [2 3];
FI{6} = 4;
FI{7} = [3 4 5];
FI{8} = [];
FI{9} = [];
FI{10} = [];
FI{11} = 6;
FI{12} = 7;

% primary output has Cin capacitance (but has no Cload)
Cin_po = sparse(M,1);
Cin_po(primary_outputs) = [10 10];

% primary input has Cload capacitance (but has no Cin)
Cload_pi = sparse(M,1);
Cload_pi(primary_inputs) = [10 10 10];

% activity frequency of gates and primary inputs
f_gates = 0.001*ones(m,1);
f_pi = sparse(M,1);
f_pi(primary_inputs) = 0.001*[10 10 10];

%********************************************************************
% derived problem data (computed from user inputs)
%********************************************************************
% fan-out cell array (compute it from the fan-in cell array)
FO = cell(M,1);
for gate = [1:m primary_outputs]
  preds = FI{gate};
  for k = 1:length(preds)
    FO{preds(k)}(end+1) = gate;
  end
end

% input and internal capacitance of gates, and driving resistance
Cin_norm  = [gates.Cin]';
Cint_norm = [gates.Cint]';
Rdrv_norm = [gates.Rdrv]';

% area specification for each gate with unit scaling
A_norm = [gates.A]';

% leakage current of gate i with unit scaling
Ileak_norm = [gates.Ileak]';

%********************************************************************
% optimization (with tradeoff curve generation)
%********************************************************************
% objective is the upper bound on the overall delay
% and that is the max of arrival times for output gates
output_gates = [FI{primary_outputs}];

% varying parameters for the tradeoff curve
N = 25;
Pmax = linspace(10,20,N);
min_delay = zeros(N,1);
disp('Generating the optimal tradeoff curve...')
for n = 1:N
  fprintf('Pmax = %6.2f: ',Pmax(n));
  cvx_begin gp quiet
    % optimization variables
    variable x(m)                 % scale factor
    variable T(m)                 % arrival times

    % input capacitance is an affine function of sizes
    Cin  = Cin_norm.*x;
    Cint = Cint_norm.*x;

    % driving resistance is inversily proportional to sizes
    R = Rdrv_norm./x;

    % gate delay is the product of its driving resistance and load cap.
    Cload = cvx( zeros(m,1) );
    for gate = 1:m
      if ~ismember( FO{gate}, primary_outputs )
        Cload(gate) = sum( Cin(FO{gate}) );
      else
        Cload(gate) = Cin_po( FO{gate} );
      end
    end

    % delay
    D = 0.69*ones(m,1).*R.*( Cint + Cload );

    % total area
    area = A_norm'*x;

    % total power calculation
    Pdyn = Vdd^2*sum( f_pi(primary_inputs).*Cload_pi(primary_inputs) ) + ...
           Vdd^2*(f_gates'*(Cint + Cload));
    Pstat = Vdd*Ileak_norm'*x;
    power = Pdyn + Pstat;

    minimize( max( T(output_gates) ) )
    subject to
      % constraints
      x >= 1; %#ok
      area <= Amax; %#ok
      power <= Pmax(n); %#ok

      % create timing constraints
      for gate = 1:m
        if ~ismember( FI{gate}, primary_inputs )
          for j = FI{gate}
            % enforce T_j + D_j <= T_i over all gates j that drive i
            D(gate) + T(j) <= T(gate); %#ok
          end
        else
          % enforce D_i <= T_i for gates i connected to primary inputs
          D(gate) <= T(gate); %#ok
        end
      end
  cvx_end
  fprintf( 'delay = %3.2f\n', cvx_optval );
  min_delay(n) = cvx_optval;
end

% plot the tradeoff curve
figure, clf
plot(Pmax,min_delay);
xlabel('Pmax'); ylabel('Dmin');
title(['Tradeoff curve for Amax = ' num2str(Amax)])
disp('Optimal tradeoff curve plotted.')

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
% this script will plot the cumulative moment
% release as a function of time

%  Stefan Wiemer  2/95

report_this_filefun(mfilename('fullpath'));

% open a new figure
figure
set(gcf,'PaperPosition',[2 1 5.5 7.5])

% create the buttons

matdraw

%  Do the calculation
%  newt2 is the currently selected catalog, newt2.Magnitude is the
% vextor containing the magnitudes
c = cumsum( 10.^(1.5*newt2.Magnitude + 16.1));

% plot the results in an xy plot
pl = plot(newt2.Date,c);
set(pl,'LineWidth',2.0)
xlabel('Time in years ','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
ylabel('Cumulative Moment ','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)

% add text -  maybe
%te = text(0.1,0.9,'log10(Mo) = 1.5Ms + 16.1;','Units','normalized','FontWeight','bold')

% change the layout of the axes slightly
set(gca,'visible','on','FontSize',ZmapGlobal.Data.fontsz.m,'FontWeight','bold',...
    'FontWeight','bold','LineWidth',1.5,...
    'Box','on')

hold on
grid

author: CelsoReyes


text:
function h=m_text(long,lat,varargin)
% M_TEXT Text Annotation
%    M_TEXT(LONG,LAT,'string') adds the text in the quotes to location
%    (LONG,LAT) on the currently defined map projection. If LONG and LAT
%    are vectors, M_TEXT writes the text at all locations given.
%    If 'string' is an array the same number of rows as the
%    length of LONG and LAT, M_TEXT marks each point with the
%    corresponding row of the 'string' array.
%
%    M_TEXT returns a column vector of handles to TEXT objects, one
%    handle per text object. TEXT objects are children of AXES objects.
%
%    M_TEXT(LONG,LAT,'string',property/value pairs) can be used to
%    change fontsize, weight, color, etc using the standard TEXT
%    properties.
%
%    See also TEXT.

% Rich Pawlowicz (rich@ocgy.ubc.ca) 17/Jan/1998
%
% This software is provided "as is" without warranty of any kind. But
% it's mine, so you can't sell it.

% 31/Jul/99 - changed to allow for X/Y vectors.
% 6/Nov/00 - eliminate returned stuff if ';' neglected (thx to D Byrne)

global MAP_PROJECTION

% Have to have initialized a map first

if isempty(MAP_PROJECTION)
  disp('No Map Projection initialized - call M_PROJ first!');
  return;
end

[X,Y]=m_ll2xy(long,lat,'clip','off');
%h=text('position',[X(:) Y(:)],'tag','m_text','string',varargin{:});
% Fix to allow vectors of X/Y to work
h=text(X(:),Y(:),varargin{:});
if ~isempty(h) && isempty(get(h(1),'tag'))
 set(h,'tag','m_text');
end

if nargout==0
 clear h
end

author: goGPS-Project


text:
function E = energyVector(C,x,V,P,varargin)
% Calculates Energy velocity vector (km/s)
%
% Description
% Energy velocity for lossless elastic medium (i.e. no attenuation)
% Good proxy for group velocity, which typically has some energy loss
% The formula is given by
% F.I. Fedorov(1968)Theory of Elastic Waves in Crystals, 375 pp. New York: Penum Press.
%
% Ve_i = C_ijkl P_j P_l X_k / rho*V
%
% N.B. E_magnitude should be equal or more than plane wave velocity vp, vs1 or vs2
%
% David Mainprice 6/02/2018
%
% Syntax
%   E = energyVector(C,x,v,p)
%   E = energyVector(C,x,v,p,rho)
%   E = energyVector(C,[],vFun,pFun)
%   E = energyVector(C,[],vFun,pFun,rho)
%
% Input
%  C - @stiffnessTensor (units GPa)
%  x - @vector3d propagation direction
%  v - plane wave velocity (unit km/s) e.g. vp,vs1 or vs2
%  p - @vector3d plane wave polarization vector e.g. pp,ps1 or ps2
%  vFun - @S2Fun plane wave velocity (unit km/s) e.g. vp,vs1 or vs2
%  pFun - @S2AxisField plane wave polarization vector e.g. pp,ps1 or ps2
%  rho - density in g/cm3
%
% Output
%  E - Energy velocity vector (units km/s)
%

% return a function if required
if isempty(x)
  if isa(V,'S2FunTri')
    E = S2VectorFieldTri(V.tri,energyVector(C,V.vertices,V.values,P.values,varargin{:}));
  else
    E = S2VectorFieldHarmonic.quadrature(@(x) energyVector(C,x,V,P,varargin{:}),'bandwidth',128,C.CS);
  end
  return
end

% make X and P to be unit vectors
x = x.normalize;

if ~isa(P,'vector3d'), P = P.eval(x); end
if ~isnumeric(V), V = V.eval(x); end
P = P.normalize;

% get density
if nargin == 5
  rho = varargin{1};
elseif isfield(C.opt,'density')
  rho = C.opt.density;
else
  rho = 1;
  warning('No density given! I''m going to use the density rho=1.');
end

% E_vector
E = vector3d(EinsteinSum(C,[1 -2 -3 -4],P(:),-2,P(:),-4,x(:),-3))./rho ./V(:);

end

author: mtex-toolbox

repo: mtex

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


text:
function inpC = regCorrIntGradWiener(inp, Int, Noise);
% regCorrIntGradWiener - corrects for the intensity gradient applying a
%                     wiener-like filtering.
%
%    inpC = regCorrIntGradWiener(inp, Int, Noise);
%
% INPUT:
%   inp   - original inplanes
%   Int   - estimated intensity
%   Noise - estimated power-spatial distribution of the noise
%
% Oscar Nestares - 5/99
%

% minimum noise
sigma2 = mean(Noise(:))/2;

% correction
inpC = inp.*Int ./ (Int.^2 + Noise + sigma2);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%cells(.5*n-1,.5*n-1)=1;
%cells(.5*n-2,.5*n-2)=1;
%cells(.5*n-3,.5*n-3)=1;
cells = (rand(n,n))<.5 ;
%how long for each case to stability or simple oscillators
%build an image and display it
imh = image(cat(3,cells,z,z));
set(imh, 'erasemode', 'none')
axis equal
axis tight
%index definition for cell update
x = 2:n-1;
y = 2:n-1;
%Main event loop
stop= 0; %wait for a quit button push
run = 0; %wait for a draw
freeze = 0; %wait for a freeze
while (stop==0)

 if (run==1)
 %nearest neighbor sum
 sum(x,y) = cells(x,y-1) + cells(x,y+1) + ...
 cells(x-1, y) + cells(x+1,y) + ...
 cells(x-1,y-1) + cells(x-1,y+1) + ...
 cells(3:n,y-1) + cells(x+1,y+1);
 % The CA rule
 cells = (sum==3) | (sum==2 & cells);
 %draw the new image
 set(imh, 'cdata', cat(3,cells,z,z) )
 %update the step number diaplay
 stepnumber = 1 + str2num(get(number,'string'));
 set(number,'string',num2str(stepnumber))
 end

 if (freeze==1)
 run = 0;
 freeze = 0;
 end
 drawnow %need this in the loop for controls to work
end
cells = (rand(n,n))<.5 ;
%how long for each case to stability or simple oscillators
%build an image and display it
imh = image(cat(3,cells,z,z));
set(imh, 'erasemode', 'none')
axis equal
axis tight
%index definition for cell update
x = 2:n-1;
y = 2:n-1;
%Main event loop
stop= 0; %wait for a quit button push
run = 0; %wait for a draw
freeze = 0; %wait for a freeze
while (stop==0)

 if (run==1)
 %nearest neighbor sum
 sum(x,y) = cells(x,y-1) + cells(x,y+1) + ...
 cells(x-1, y) + cells(x+1,y) + ...
 cells(x-1,y-1) + cells(x-1,y+1) + ...
 cells(3:n,y-1) + cells(x+1,y+1)+...
 cells(x,y);
 % The CA rule
 cells = ~((sum<4) | (sum==5));
 %draw the new image
 set(imh, 'cdata', cat(3,cells,z,z) )
 %update the step number diaplay
 stepnumber = 1 + str2num(get(number,'string'));
 set(number,'string',num2str(stepnumber))
 end

 if (freeze==1)
 run = 0;
 freeze = 0;
 end

 drawnow %need this in the loop for controls to work

end

author: HuangCongQing


text:
function res = le(A,B)
%LE           Implements  A<=B  elementwise for long (refers only to midpoints)
%

% written  11/06/99     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  C = A - B;
  res = ( C.sign==-1 ) | all( C.mantissa==0 , 2);

repo: OpenMAS


text:
function K = sdlfmaXsdlfmaKernComputeBlock(lfmKern1, lfmKern2, t1, t2, ...
    kyy, kyv, kvy, kvv, i, j, generalConst)

% SDLFMAXSDLFMAKERNCOMPUTEBLOCK Computes SDLFM kernel matrix for block i,j
% FORMAT
% DESC computes the kernel matrix for the SDLFM kernel function in the
% block specified at indeces i,j. It assumes the computation for functions
% that describe accelerations (acceleration 1 and acceleration 2).
% ARG lfmKern1 : structure containing parameters for the system 1
% ARG lfmKern2 : structure containing parameters for the system 2
% ARG t1 : times at which the system 1 is evaluated
% ARG t2 : times at which the system 2 is evaluated
% ARG kyy : covariance for the initial conditions between position 1 and
% position 2 at block i,j
% ARG kyv : covariance for the initial conditions between position 1 and
% velocity 2 at block i,j
% ARG kvy : covariance for the initial conditions between velocity 1 and
% position 2 at block i,j
% ARG kvv : covariance for the initial conditions between velocity 1 and
% velocity 2 at block i,j
% ARG i : interval to be evaluated for system 1
% ARG j : interval to be evaluated for system 2
% ARG generalConstant : constants evaluated with sdlfmKernComputeConstant.m
% RETURN K : the kernel matrix portion of block i,j
%
% COPYRIGHT : Mauricio A. Alvarez, 2010.

% KERN

if nargin<11
    j = i;
    generalConst = [];
end

a1 = sdlfmaMeanCompute(lfmKern1(1), t1, 'Pos');
b1 = sdlfmaMeanCompute(lfmKern1(1), t1, 'Vel');
a2 = sdlfmaMeanCompute(lfmKern2(1), t2, 'Pos');
b2 = sdlfmaMeanCompute(lfmKern2(1), t2, 'Vel');

K = kyy*a1*a2.' + kyv*a1*b2.' + kvy*b1*a2.' + kvv*b1*b2.';

if i==j
    for k=1:length(lfmKern1)
        K  = K + lfmaXlfmaKernCompute(lfmKern1(k), lfmKern2(k), t1, t2);
    end
else
    if i>j
        AccelPos = zeros(1, length(t2));
        AccelVel = zeros(1, length(t2));
        for k=1:length(lfmKern1)
            AccelPos = AccelPos + lfmaXlfmKernCompute(lfmKern2(k), lfmKern1(k), t2, lfmKern2(k).limit).';
            AccelVel = AccelVel + lfmaXlfmvKernCompute(lfmKern2(k), lfmKern1(k), t2, lfmKern2(k).limit).';
        end
        if isempty(generalConst{i,j})
            K = K + a1*AccelPos + b1*AccelVel;
        else
            K = K + (generalConst{i,j}(1,1)*a1 + generalConst{i,j}(2,1)*b1)*AccelPos + ...
                (generalConst{i,j}(1,2)*a1 + generalConst{i,j}(2,2)*b1)*AccelVel;
        end
    else
        AccelPos = zeros(length(t1),1);
        AccelVel = zeros(length(t1),1);
        for k =1:length(lfmKern1)
            AccelPos = AccelPos + lfmaXlfmKernCompute(lfmKern1(k), lfmKern2(k), t1, lfmKern1(k).limit);
            AccelVel = AccelVel + lfmaXlfmvKernCompute(lfmKern1(k), lfmKern2(k), t1, lfmKern1(k).limit);
        end
        if isempty(generalConst{i,j})
            K = K + AccelPos*a2.' + AccelVel*b2.';
        else
            K = K + AccelPos*(generalConst{i,j}(1,1)*a2.' + generalConst{i,j}(2,1)*b2.') + ...
                AccelVel*(generalConst{i,j}(1,2)*a2.' + generalConst{i,j}(2,2)*b2.');
        end
    end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% ConvertToGray.m
% -------------------------------------------------------------------
%
% Authors: Sun Li
% Date:    05/05/2013
% Last modified: 05/05/2013
% -------------------------------------------------------------------

function imgGray = ConvertToGray(img)
    if size(img, 3) == 3,
        imgGray = RGBTOGRAY(img);
    else
        imgGray = img;
    end
end

author: thfylsty


text:
function [yi, xi] = downsample_scnlab(y, orig_samprate, new_samprate, varargin)
% Uses linear interpolation to resample a vector from one sampling
% rate to another
%
% :Usage:
% ::
%
%     [yi, xi] = downsample_scnlab(y, orig_samprate, new_samprate, [doplot])
%
% :Inputs:
%
%   **orig_samprate:**
%        sampling rate in Hz
%
%   **new_samprate:**
%        desired sampling rate in Hz
%
% ..
%    I prefer this to matlab's downsample.m because linear interpolation is
%    robust and does fewer weird things to the data.
%
%     Tor Wager, June 2009
% ..
%
% :Example:
% ::
%
%    % Downsample a 100 Hz signal to a scanning TR of 2 sec
%    % signal at 100 Hz, sample to low-freq TR of 0.5 hz (2 sec TR)
%    % every 100 / TR = 100/.5 = 200 samples
%
%    [yi, xi] = downsample_scnlab(y, 100, .5)
%

doplot = 0;
if length(varargin) > 0, doplot = varargin{1}; end

downsamplerate = orig_samprate ./ new_samprate;

nobs = length(y);

% observations in resampled output
nobs_out = round(nobs ./ downsamplerate);

x = 1:nobs;

xi = linspace(0, nobs, nobs_out);

yi = interp1(x, y, xi, 'linear', 'extrap');

if doplot

    create_figure('Downsample plot');
    plot(y)
    hold on; plot(xi, yi, 'r')

end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%% Demo: DISTATIS
%
% The data used here is available from http://cosmomvpa.org/datadb.zip
%
% It is based on the following work:
% * Connolly et al (2012), Representation of biological classes in the human
%   brain. Journal of Neuroscience, doi 10.1523/JNEUROSCI.5547-11.2012
%
% Six categories (monkey, lemur, mallard, warbler, ladybug, lunamoth)
% during ten runs in an fMRI study.
%
% This example shows the application of DISTATIS, which tries to find an
% optimal 'compromise' dissimilarity matrix across a set of observations
% (participants)
%
% Reference:
%   - Abdi, H., Valentin, D., O?Toole, A. J., & Edelman, B. (2005).
%     DISTATIS: The analysis of multiple distance matrices. In
%     Proceedings of the IEEE Computer Society: International conference
%     on computer vision and pattern recognition, San Diego, CA, USA
%     (pp. 42?47).
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

%% Set data paths
% The function cosmo_config() returns a struct containing paths to tutorial
% data. (Alternatively the paths can be set manually without using
% cosmo_config.)
config=cosmo_config();
study_path=fullfile(config.tutorial_data_path,'ak6');
output_path=config.output_data_path;

readme_fn=fullfile(study_path,'README');
cosmo_type(readme_fn);

% reset citation list
cosmo_check_external('-tic');

%% Preprocessing for DISTATIS: RSM analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subject_ids={'s01','s02','s03','s04','s05','s06','s07','s08'};
nsubjects=numel(subject_ids);

mask_label='vt_mask';

ds_rsms=cell(nsubjects,1); % allocate space for output
for subject_num=1:nsubjects
    subject_id=subject_ids{subject_num};

    % Code from here is pretty much identical to that above >>>

    % set path for this subject
    data_path=fullfile(study_path,subject_id);

    % Define data locations and load data from even and odd runs
    mask_fn=fullfile(data_path, [mask_label '.nii']); % vt mask

    % Use odd runs only
    data_fn=fullfile(data_path,'glm_T_stats_odd.nii');
    ds=cosmo_fmri_dataset(data_fn,'mask',mask_fn,...
                            'targets',1:6,'chunks',1);

    ds_rsm=cosmo_dissimilarity_matrix_measure(ds);

    % set chunks (one chunk per subject)
    ds_rsm.sa.chunks=subject_num*ones(size(ds_rsm.samples,1),1);
    ds_rsms{subject_num}=ds_rsm;
end

% combine data from all subjects
all_ds=cosmo_stack(ds_rsms);

%% Run DISTATIS
distatis=cosmo_distatis(all_ds);

%% show comprimise distance matrix
[compromise_matrix,dim_labels,values]=cosmo_unflatten(distatis,1);

labels={'monkey', 'lemur', 'mallard', 'warbler', 'ladybug', 'lunamoth'};
n_labels=numel(labels);
figure();
imagesc(compromise_matrix)
title('DSM');
set(gca,'YTick',1:n_labels,'YTickLabel',labels);
set(gca,'XTick',1:n_labels,'XTickLabel',labels);
ylabel(dim_labels{1});
xlabel(dim_labels{2});
colorbar

% skip if stats toolbox is not present
if cosmo_check_external('@stats',false)
    figure();
    hclus = linkage(compromise_matrix);
    dendrogram(hclus,'labels',labels,'orientation','left');
    title('dendrogram');

    figure();
    F = cmdscale(squareform(compromise_matrix));
    text(F(:,1), F(:,2), labels);
    title('2D MDS plot');
    mx = max(abs(F(:)));
    xlim([-mx mx]); ylim([-mx mx]);
end

%% show citation information
cosmo_check_external('-cite');

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function value = daub2_condition ( n )

%*****************************************************************************80
%
%% DAUB2_DETERMINANT returns the L1 condition of the DAUB2 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real VALUE, the L1 condition.
%
  c0 = sqrt ( 2.0 ) / 2.0;
  c1 = sqrt ( 2.0 ) / 2.0;

  a_norm = abs ( c0 ) + abs ( c1 );
  b_norm = a_norm;
  value = a_norm * b_norm;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function skypos = place (tjd, object, locatn, icoord, star, observ)

% this function computes the apparent direction of a star or solar
% system body at a specified time and in a specified coordinate
% system.  based on kaplan, et al. (1989), astronomical journal 97,
% 1197-1210, with some enhancements from klioner (2003),
% astronomical journal 125, 1580-1597.

%  tjd    = tt julian date for place (in)

%  object = character string identifying object of interest (in)

%           for solar system
%           body,             specify the name using all upper-
%                             case letters ('sun', 'moon',
%                             'jupiter', etc.),
%                             - or -
%                             specify the body id number
%                             in a 4-character string of the
%                             form '=nnn', where nnn is the
%                             body id number
%           for star,         provide a blank string, the word
%                             'star', or any string beginning
%                             with '*'

%  locatn = integer code specifying location of observer (in)

%           set locatn=0 for observer at geocenter
%           set locatn=1 for observer on surface of earth
%           set locatn=2 for observer on near-earth spacecraft

%  icoord = integer code specifying coordinate system of output
%           position (in)

%           set icoord=0 for gcrs (or 'local gcrs')
%           set icoord=1 for true equator and equinox of date
%           set icoord=2 for true equator and cio of date
%           set icoord=3 for astrometric coordinates, i.e.,
%                        without light deflection or aberration

%  star   = array of catalog data for star (in)

%           (not used if solar system body requested)

%           star(1) = icrs right ascension in hours
%           star(2) = icrs declination in degrees
%           star(3) = icrs proper motion in ra in
%                     milliarcseconds/year
%           star(4) = icrs proper motion in dec in
%                     milliarcseconds/year
%           star(5) = parallax in milliarcseconds
%           star(6) = radial velocity in kilometers/second
%           further star array elements are not used here
%           but are reserved for future use

%  observ = array of data specifying location of observer (in)

%           (not used if locatn=0)

%           for locatn = 1,

%           observ(1) = geodetic longitude (wgs-84) of observer
%                       (east +) in degrees
%           observ(2) = geodetic latitude (wgs-84) of observer
%                       (north +) in degrees
%           observ(3) = height of observer above ellipsoid
%                       in meters
%           observ(4) = value of delta-t in seconds
%                       (delta-t=tt-ut1)
%           observ(5) = (not used, reserved for future use)
%           observ(6) = (not used, reserved for future use)

%           for locatn = 2,

%           observ(1) = geocentric x in kilometers
%           observ(2) = geocentric y in kilometers
%           observ(3) = geocentric z in kilometers
%           observ(4) = geocentric x-dot in kilometers/second
%           observ(5) = geocentric y-dot in kilometers/second
%           observ(6) = geocentric z-dot in kilometers/second
%           with respect to true equator and equinox of date

%  skypos = array of output data specifying object's place
%           on the sky at time tjd, with respect to the
%           specified output coordinate system (out)

%           skypos(1) = x, dimensionless      unit vector
%           skypos(2) = y, dimensionless      toward object
%           skypos(3) = z, dimensionless
%           skypos(4) = apparent, topocentric, or astrometric
%                       right ascension in hours
%           skypos(5) = apparent, topocentric, or astrometric
%                       declination in degrees
%           skypos(6) = true (geometric, euclidian) distance
%                       to solar system body in au at time tjd,
%                       or 0.0d0 for star
%           skypos(7) = radial velocity in kilometers/second

%           further skypos array elements are not used here
%           but are reserved for future use

% note 1: values of locatn and icoord for various standard kinds of place:

%   locatn = 0 and icoord = 1 apparent place
%   locatn = 1 and icoord = 1 topocentric place
%   locatn = 0 and icoord = 0 virtual place
%   locatn = 1 and icoord = 0 local place
%   locatn = 0 and icoord = 3 astrometric place
%   locatn = 1 and icoord = 3 topocentric astrometric place

% note 2: arrays star and skypos may be expanded in the future, and
% this can be allowed for in the calling code by dimensioning
% these arrays with 20 and 10 elements, respectively, even though
% elements beyond star(6) and skypos(7) are not now referred to in
% this subroutine.

% note 3: if locatn = 1 and observ(4) = 0.0d0, the value of delta-t will
% be obtained from getdt, which provides the last value of delta-t
% defined by the user via call to setdt.

% note 4: skypos(7), the radial velocity, is the predicted
% radial velocity measure (z) times the speed of light, an
% inherently spectroscopic measure. for a star, it
% includes all effects, such as gravitational red shift,
% contained in the catalog barycentric radial velocity measure,
% which is assumed given in star(6). for a solar system
% body, it applies to a fictitious emitter at the center of the
% observed object, assumed massless (no gravitational red shift),
% and does not in general apply to reflected light.

% ported from NOVAS 3.1

%%%%%%%%%%%%%%%%%%%%%%%

% t0 = tdb julian date of epoch j2000.0 (tt)

t0 = 2451545.0d0;

% get body ids

iearth = idss_novas ('earth');

isun = idss_novas ('sun');

% get c, the speed of light in au/day

c = astcon ('c(au/day)', 1.0d0);

% check on earth as an observed object

if (strcmp(object, '=3') == 1 && locatn ~= 2)

    fprintf ('\n place: will not process earth as observed object except when locatn = 2');

    return

end

% compute tdbjd, the tdb julian date corresponding to ttjd

ttjd = tjd;

tdbjd = tjd;

[x, secdif] = novas_times (tdbjd);

tdbjd = ttjd + secdif / 86400.0d0;

% get position and velocity of the earth wrt barycenter of solar system, in icrs

% fprintf('\n\nearth state vector wrt barycenter\n');

[peb, veb, ierr] = solsys (tdbjd, iearth, 0);

if (ierr ~= 0)

    fprintf ('\nplace: cannot obtain coordinates of earth at jd %16.8f', tjd);

    return

end

% get position and velocity of the sun wrt barycenter of solar system, in icrs

% fprintf('\n\nsun state vector wrt barycenter\n');

[psb, vsb, ierr] = solsys (tdbjd, isun, 0);

if (ierr ~= 0)

    fprintf ('\nplace: cannot obtain coordinates of sun at jd %16.8f', tjd);

    return

end

% get position and velocity of observer

if (locatn == 1 || locatn == 2)

    % for topocentric place, get geocentric position and velocity
    % vectors of observer

    % fprintf('\n\nobserver state vector wrt geocenter\n');

    [pog, vog] = geopos (ttjd, locatn, observ);

    loc = 1;

else

    % for geocentric place, there is nothing to do

    for j = 1:3

        pog(j) = 0.0d0;

        vog(j) = 0.0d0;

    end

    loc = 0;

end

% compute position and velocity of observer wrt barycenter of
% solar system (galilean transformation from gcrs to bcrs)

for j = 1:3

    pob(j) = peb(j) + pog(j);

    vob(j) = veb(j) + vog(j);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find geometric position of observed object
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (strcmp(object, 'star') == 1 || strcmp(object, ' ') == 1 || strcmp(object(1:1), '*'))

    %%%%%%%%%%%%%%%%%%%%%%%%%
    % observed object is star
    %%%%%%%%%%%%%%%%%%%%%%%%%

    idbody = -9999;

    % get position of star updated for its space motion

    [pos1, vel1] = vectrs (star(1), star(2), star(3), star(4), star(5), star(6));

    dt = dlight (pos1, pob);

    pos2 = propmo (t0, pos1, vel1, tdbjd + dt);

    % get position of star wrt observer (corrected for parallax)

    [pos3, tlight] = geocen (pos2, pob);

    dis = 0.0d0;

else

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % observed object is solar system body
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % get id number of body

    if (strcmp(object(1:1), '=') == 1)

        idbody = object(2:length(object));

    else

        idbody = idss_novas (object);

        if (idbody == -9999)

            fprintf ('\nplace: cannot obtain coordinates of object at jd %16.8f', tjd);

            return

        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % get position of body wrt barycenter of solar system
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % fprintf('\n\nobject state vector wrt barycenter\n');

    [pos1, vel1, ierr] = solsys (tdbjd, str2num(idbody), 0);

    if (ierr ~= 0)

        fprintf ('\nplace: cannot obtain coordinates of object at jd %16.8f', tjd);

        return

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % get position of body wrt observer, and true (euclidian) distance
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % fprintf('\n\nobject state vector wrt observer\n');

    [pos2, tlight] = geocen (pos1, pob);

    dis = tlight * c;

    % get position of body wrt observer, antedated for light-time

    % fprintf('\n\nobject state vector wrt observer - littim\n');

    [pos3, tlight] = littim (tdbjd, idbody, pob, 0.0);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% apply gravitational deflection of light and aberration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (icoord == 3)

    % these calculations are skipped for astrometric place

    for j = 1:3

        pos5(j) = pos3(j);

    end

else

    % variable loc determines whether earth deflection is included

    if (loc == 1)

        [x, frlimb] = limang (pos3, pog);

        if (frlimb < 0.8d0)

            loc = 0;

        end

    end

    % compute gravitational deflection and aberration

    pos4 = grvdef (tdbjd, loc, pos3, pob);

    pos5 = aberat (pos4, vob, tlight);

    % position vector is now in gcrs

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% transform, if necessary, to output coordinate system
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (icoord == 1)

    % transform to equator and equinox of date

    pos6 = frame (pos5, 1);

    pos7 = preces (t0, pos6, tdbjd);

    pos8 = nutate (tdbjd, pos7);

elseif (icoord == 2)

    % transform to equator and cio of date

    % obtain the basis vectors, in the gcrs, of the celestial
    % intermediate system

    kcio = cioloc (tdbjd, rcio);

    [px, py, pz] = ciobas (tdbjd, rcio, kcio);

    % transform position vector to celestial intermediate system

    pos8(1) = px(1) * pos5(1) + px(2) * pos5(2) + px(3) * pos5(3);

    pos8(2) = py(1) * pos5(1) + py(2) * pos5(2) + py(3) * pos5(3);

    pos8(3) = pz(1) * pos5(1) + pz(2) * pos5(2) + pz(3) * pos5(3);

else

    % no transformation -- keep coordinates in gcrs (or icrs for astrometric place)

    for j = 1:3

        pos8(j) = pos5(j);

    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set up star data, if applicable
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (idbody == -9999)

    rvs(1) = star(1);

    rvs(2) = star(2);

    rvs(3) = star(6);

    if (star(5) <= 0.0d0)

        vel1(1) = 0.0d0;

        vel1(2) = 0.0d0;

        vel1(3) = 0.0d0;

    end

else

    rvs(1) = 0.0d0;

    rvs(2) = 0.0d0;

    rvs(3) = 0.0d0;

end

% compute distances: observer-geocenter, observer-sun, object-sun

rvd(1) = sqrt((pob(1) - peb(1))^2 + (pob(2) - peb(2))^2 + (pob(3) - peb(3))^2);

rvd(2) = sqrt((pob(1) - psb(1))^2 + (pob(2) - psb(2))^2 + (pob(3) - psb(3))^2);

rvd(3) = sqrt((pos1(1) - psb(1))^2 + (pos1(2) - psb(2))^2 + (pos1(3) - psb(3))^2);

rv = radvl (pos3, vel1, vob, rvs, rvd);

% finish up

[ra, dec] = angles (pos8);

x = sqrt (pos8(1)^2 + pos8(2)^2 + pos8(3)^2);

for j = 1:3

    skypos(j) = pos8(j) / x;

end

skypos(4) = ra;

skypos(5) = dec;

skypos(6) = dis;

skypos(7) = rv;


text:
function D=distMat(P1, P2)
%
% Euclidian distances between vectors
% each vector is one row

if nargin == 2
    P1 = double(P1);
    P2 = double(P2);

    X1=repmat(sum(P1.^2,2),[1 size(P2,1)]);
    X2=repmat(sum(P2.^2,2),[1 size(P1,1)]);
    R=P1*P2';
    D=real(sqrt(X1+X2'-2*R));
else
    P1 = double(P1);

    % each vector is one row
    X1=repmat(sum(P1.^2,2),[1 size(P1,1)]);
    R=P1*P1';
    D=X1+X1'-2*R;
    D = real(sqrt(D));
end

author: willard-yuan

repo: hashing-baseline-for-image-retrieval


text:
function [] = mci_plot_dist (dist,j,xlims)
% Plot probability density
% FORMAT [] = mci_plot_dist (dist,j,xlims)
%
% dist      struct with fields
%
% .Ep       posterior mean
% .P        [Np x Ns] sample matrix
% .ind      indices of samples dist burn-in
% .names
% .ks       set to 1 for kernel smoothing (default)
% j         jth variable
% xlims     xlims(1,2) for lower/upper limits
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: mci_plot_dist.m 6697 2016-01-27 14:57:28Z spm $

try, ks=dist.ks; catch, ks=1; end

if nargin > 2
    limits=1;
else
    limits=0;
end
lw=2;

switch lower(dist.type),
    case 'boxplot',
        boxplot(dist.P(j,dist.ind));
    case 'sample',
        if ks
            [g,xi]=ksdensity(dist.P(j,dist.ind));
            plot(xi,g,dist.color,'LineWidth',lw);set(gca,'YTick',[]);
        else
            hist(dist.P(j,dist.ind),20);
        end
    case 'gaussian',
        m=dist.Ep(j);s=sqrt(full(dist.Cp(j,j)));
        xi=linspace(m-4*s,m+4*s,100);
        g=spm_Npdf(xi,m,s^2);
        plot(xi,g,dist.color,'LineWidth',lw);set(gca,'YTick',[]);
end
set(gca,'FontSize',16);
xlabel(dist.names{j});
grid on

if limits
    xlim([xlims(1) xlims(2)]);
end

author: spm


text:
% trainTestSurfaceContextScript
%
% Assumes you have created:
%   imsegs  :  data structure that contains information on all training
%              and testing images
%   trainind:  the indices of the training images
%   testind:   the indices of the testing images

basedir = '???'; % set this appropriately
imdir = [basedir '/images/all_images'];
datadir = [basedir '/data'];
resultsdir = [basedir '/results'];

DO_TRAIN = 1;
DO_TEST = 1;

if DO_TRAIN

    nsegments = [10 20 30 40 50 60 80 100];

    % Computes the superpixel features
    if ~exist('spfeatures')
        spfeatures = mcmcGetAllSuperpixelData(imdir, imsegs);
        save([datadir '/superpixelData.mat'], 'spfeatures');
    end

    % Trains the superpixel classifiers
    if ~exist('vclassifierSP')
        [vclassifierSP, hclassifierSP] = ...
            mcmcTrainSuperpixelClassifier(spfeatures(trainind), imsegs(trainind));
        save([datadir '/superpixelClassifier.mat'], 'vclassifierSP', 'hclassifierSP');
    end

    % Computes the same-label features
    if ~exist('efeatures')
        [efeatures, adjlist] = mcmcGetAllEdgeData(spfeatures, imsegs);
        save([datadir '/edgeData.mat'], 'efeatures', 'adjlist');
    end

    % Trains the same-label classifier
    if ~exist('eclassifier')
        eclassifier = mcmcTrainEdgeClassifier(efeatures(trainind), ...
            adjlist(trainind), imsegs(trainind));
        ecal = calibrateEdgeClassifier(efeatures(trainind), adjlist(trainind), ...
            imsegs(trainind), eclassifier, 1);
        ecal = ecal{1};
        save([datadir '/edgeClassifier.mat'], 'eclassifier', 'ecal');
    end

    % Computes the multiple segmentations, the segment features, and the
    % ground truth labels for each segment
    if ~exist('labdata')
        % gather data
        for f = 1:numel(imsegs)

            disp([num2str(f) ': ' imsegs(f).imname])

            [pvSP{f}, phSP{f}, pE{f}] = mcmcInitialize(spfeatures{f}, efeatures{f}, ...
                adjlist{f}, imsegs(f), vclassifierSP, hclassifierSP, eclassifier, ecal, 'none');
            smaps{f} = generateMultipleSegmentations2(pE{f}, adjlist{f}, imsegs(f).nseg, nsegments);

            im = im2double(imread([imdir '/' imsegs(f).imname]));
            imdata = mcmcComputeImageData(im, imsegs(f));

            for k = 1:numel(nsegments)
                labdata{f, k} = mcmcGetSegmentFeatures(imsegs(f), spfeatures{f}, imdata, smaps{f}(:, k), (1:max(smaps{f}(:, k))));
                [mclab{f, k}, mcprc{f, k}, allprc{f, k}, trainw{f, k}] = segmentation2labels(imsegs(f), smaps{f}(:, k));
                unilabel{f, k} = mclab{f, k}.*(mcprc{f, k}>0.95);
                seglabel{f,k} =  1*(mcprc{f, k}>0.95) + (-1)*(mcprc{f, k}<0.95);
            end
        end
        save([datadir '/allData.mat'], 'smaps', 'labdata', 'segdata', 'mclab', 'mcprc', 'allprc', 'seglabel', 'unilabel', 'trainw', 'pvSP', 'phSP', 'pE');
    end

    % Trains the segment classifiers
    if ~exist('vclassifier')
        sclassifier = mcmcTrainSegmentationClassifier2(labdata(trainind, :), seglabel(trainind{k}, :), trainw(trainind, :));
        [vclassifier, hclassifier] = ...
            mcmcTrainSegmentClassifier2(labdata(trainind, :), unilabel(trainind, :), trainw(trainind, :), 50000);
    end

    save([datadir '/allClassifiers.mat'], 'vclassifier', 'hclassifier', 'sclassifier', 'eclassifier', 'vclassifierSP', 'hclassifierSP');

end

if DO_TEST

    % Computes the label confidences for each superpixel in the test images
    % and gives the final accuracy.
    % pg{image number}(superpixel number, [000 left center right porous solid sky])
    % gives the superpixel label confidences
    [vacc, hacc, vcm, hcm, pg] = ...
        testMultipleSegmentationsCV2(imsegs(testind), labdata(testind, :), ...
        labdata(testind, :), smaps(testind), ...
        vclassifier, hclassifier, sclassifier, pvSP(testind), phSP(testind), 1);
    save([datadir '/results.mat'], 'vacc', 'hacc', 'vcm', 'hcm', 'pg');

    % Computes and writes the labeled images
    for f = 1:testind
        im = im2double(imread([imdir '/' imsegs(f).imname]));
        [pv, ph] = splitpg(pg{f});
        lim = APPgetLabeledImage2(im, imsegs(f), pv, ph);
        imwrite(im, [resultsdir '/' imsegs(f).imname]);
        imwrite(lim, [resultsdir '/' strtok(imsegs(f).imname, '.') '.l.jpg']);
    end
end

author: zouchuhang

repo: LayoutNet


text:
% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

%%begin

% Generate a Reeds-Shepp path for 3-point turn
dl = 0.05;
q0 = [0 0 0]'; qf = [0 0 pi]';
maxcurv = 1/5;   % 5m turning circle
rs = ReedsShepp(q0, qf, maxcurv, dl)

% set up a vehicle model
[car.image,~,car.alpha] = imread('car2.png');
car.rotation = 180;
car.centre = [81,110];
car.centre = [648; 173];
car.length = 4.2;

% now animate
clf; plotvol([-4 8 -6 6])

a = gca;
a.XLimMode = 'manual';
a.YLimMode = 'manual';
set(gcf, 'Color', 'w')
grid on
a = gca;
xyzlabel

plot_vehicle(rs.path, 'model', car, 'trail', 'r:');

author: petercorke

repo: robotics-toolbox-matlab


text:
function [ x, more ] = sgmga_vcn_naive ( dim_num, level_weight, x_max, x, ...
  q_min, q_max, more )

%*****************************************************************************80
%
%% SGMGA_VCN_NAIVE returns the next constrained vector.
%
%  Discussion:
%
%    This function uses a naive algorithm, which quickly becomes unsuitable
%    for higher dimensions.  The function SGMGA_VCN is an attempt at a more
%    efficient calculation of the same quantities.
%
%    We consider vectors X of dimension DIM_NUM satisfying:
%
%      0 <= X(1:DIM_NUM) <= X_MAX(1:DIM_NUM).
%
%    and define
%
%      Q = sum ( 1 <= I <= DIM_NUM ) LEVEL_WEIGHT(I) * X(I)
%
%    and seek X satisfying the constraint:
%
%      Q_MIN < Q <= Q_MAX
%
%    For sparse grid applications, we compute
%
%      LEVEL_WEIGHT_MIN_POS = minimum positive entry in LEVEL_WEIGHT
%
%    and assume there is an underlying LEVEL used to index the sets of
%    constrained vectors, and that
%
%      Q_MAX = LEVEL * LEVEL_WEIGHT_MIN_POS
%      Q_MIN = LEVEL - LEVEL_WEIGHT_MIN_POS * sum ( LEVEL_WEIGHT(:) )
%      X_MAX(I) = LEVEL * LEVEL_WEIGHT_MIN_POS / LEVEL_WEIGHT(I)
%
%    This routine returns, one at a time exactly those X which satisfy
%    the constraint.  No attempt is made to return the X values in
%    any particular order as far as Q goes.
%
%  Example:
%
%    LEVEL_WEIGHT:          1.000000        1.000000
%
%    Q_MIN:        0.000000
%    Q_MAX:        2.000000
%    X_MAX:                         2         2
%
%         1        1.000000         1         0
%         2        2.000000         2         0
%         3        1.000000         0         1
%         4        2.000000         1         1
%         5        2.000000         0         2
%
%    LEVEL_WEIGHT:          1.000000        2.000000
%
%    Q_MIN:       -1.000000
%    Q_MAX:        2.000000
%    X_MAX:                         2         1
%
%         1        0.000000         0         0
%         2        1.000000         1         0
%         3        2.000000         2         0
%         4        2.000000         0         1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    An Anisotropic Sparse Grid Stochastic Collocation Method for Partial
%    Differential Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2411-2442.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the number of components in the vector.
%
%    Input, real LEVEL_WEIGHT(DIM_NUM), the anisotropic weights.
%
%    Input, integer X_MAX(DIM_NUM), the maximum
%    values allowed in each component of X.
%
%    Input, integer X(DIM_NUM).  On first call, with
%    MORE = FALSE, the input value of X is not important.  On subsequent calls,
%    the input value of X should be the output value from the previous call.
%
%    Input, real Q_MIN, Q_MAX, the lower and upper limits on the sum.
%
%    Input, logical MORE.  On input, if the user has set MORE
%    FALSE, the user is requesting the initiation of a new sequence
%    of values.  If MORE is TRUE, then the user is requesting "more"
%    values in the current sequence.
%
%    Output, integer X(DIM_NUM).
%    On output, (with MORE = TRUE), the value of X will be the "next"
%    vector in the reverse lexicographical list of vectors that satisfy
%    the condition.  However, on output with MORE = FALSE, the vector
%    X is meaningless, because there are no more vectors in the list.
%
%    Output, logical MORE.  If MORE is TRUE on output,
%    then another value was found and returned in X, but if MORE is
%    FALSE, then there are no more values in the sequence, and X is
%    NOT the next value.
%
  if ( ~more )

    more = 1;

    x = zeros ( dim_num, 1 );

    q = level_weight(1:dim_num)' * x(1:dim_num);

    if ( q_min < q && q <= q_max )
      return
    end

  end

  while ( 1 )

    i = 1;

    while ( 1 )

      if ( x(i) < x_max(i) )
        break
      end

      if ( dim_num <= i )
        more = 0;
        return
      end

      i = i + 1;

    end

    x(i) = x(i) + 1;
    x(1:i-1) = 0;

    q = level_weight(1:dim_num)' * x(1:dim_num);

    if ( q_min < q && q <= q_max )
      break
    end

  end

  return
end


text:
function v = r8vec_house_column ( n, a, k )

%*****************************************************************************80
%
%% R8VEC_HOUSE_COLUMN defines a Householder premultiplier that "packs" a column.
%
%  Discussion:
%
%    The routine returns a vector V that defines a Householder
%    premultiplier matrix H(V) that zeros out the subdiagonal entries of
%    column K of the matrix A.
%
%       H(V) = I - 2 * v * v'
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 January 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix A.
%
%    Input, real A(N,1), column K of the matrix A.
%
%    Input, integer K, the column of the matrix to be modified.
%
%    Output, real V(N,1), a vector of unit L2 norm which defines an
%    orthogonal Householder premultiplier matrix H with the property
%    that the K-th column of H*A is zero below the diagonal.
%

%
%  Destroy all row vectors!
%
  a = a(:);

  v(1:n,1) = 0.0;

  if ( k < 1 || n <= k )
    return
  end

  s = sqrt ( sum ( a(k:n,1).^2 ) );

  if ( s == 0.0 )
    return
  end

  if ( a(k,1) < 0.0 )
    v(k,1) = a(k,1) - abs ( s ) ;
  else
    v(k,1) = a(k,1) + abs ( s );
  end
  v(k+1:n,1) = a(k+1:n,1);

  s = sqrt ( sum ( v(k:n,1).^2 ) );
  v(k:n,1) = v(k:n,1) / s;

  return
end


text:
function v = cvx_value( x )
warning( 'CVX:IllegalUse', ...
    'CVX error: illegal use of a cvx problem object has been detected.\n   Please do not copy or manipulate the value of ''cvx_problem'' in any way.', 1 ); %#ok
v = [];

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function ind = gpBlockIndices(model, blockNo)

% GPBLOCKINDICES Return indices of given block.
% FORMAT
% DESC returns the indices of a given block for the PITC
% approximation.
% ARG model : the model for which the indices are being computed.
% ARG blockNum : the block number for which the indices are
% required.
% RETURN indices : the data indices associated with given block.
%
% SEEALSO : gpComputeAlpha, gpCovGrads, gpLogLikeGradients,
% gpLogLikelihood, gpUpdateAD
%
% COPYRIGHT : Neil D. Lawrence, 2006

% GP

if ~strcmp(model.approx, 'pitc')
  error('Block indices only relevant for pitc approximation.');
else
  if blockNo == 1
    startVal = 1;
  else
    startVal = model.blockEnd(blockNo-1)+1;
  end
  endVal = model.blockEnd(blockNo);
  ind = startVal:endVal;
end