text:
function view = addROIline(view,sgn)
%
% view = addROIline(view,[sgn])
%
% Click on two points in the image and find an ROI along a line
% between them.  What does line mean?  Geodesic or screen line?
%
% If sgn~=0, adds user-specified line to selected ROI in
% current slice. If sgn==0, removes the line from the ROI.
%
% If you change this function make parallel changes in:
%   all addROI*.m functions
%
% bw, 4/30/99

% error if no current ROI
if view.selectedROI == 0
  myErrorDlg('No current ROI');
  return
end

if ~exist('sgn','var')
  disp('Default:  adding coords')
  sgn = 1;
end
% This is only really meaningful in the flat view... DO a check to make
% sure it's a FLAT
% Get curSlice
curSlice = viewGet(vw, 'Current Slice');

if view.rotateImageDegrees(curSlice)>0  %detects if the Flat your are currently selecting the points ist rotated
            errordlg('Flat must not be rotated','Warning!'); %Results would be wrong!
end

% Get current ROI coords
curCoords = getCurROIcoords(view);

% Save prevCoords for undo
view.prevCoords = curCoords;

% -------- Copied and edited from measureCOrticalDistance...

    % Select flat figure and get a single point from the user
    figure(view.ui.figNum)
    disp('Click left to add points, right to quit');
    button = 1;
    count = 0;
    w = 0.5;
    coords = [];

    while(button~=3)
        [x,y,button] = ginput(1);
        if(button==3)
            break;
        end
        count = count+1;
        coords = [coords,round([y;x;curSlice])];
        h(count) = line([x-w,x-w,x+w,x+w,x-w],[y-w,y+w,y+w,y-w,y-w],'Color','w');

    end
    % Delete the temporarily drawn squares
    for ii=1:length(h)
        delete(h(ii));
    end
    clear h;

nPoints=size(coords,2);

% Check if outside image
%
% dims=size(view.ui.image);
% if (min(rgn(:,1))< 1 | max(rgn(:,1))>dims(1) | ...
%       min(rgn(:,2))< 1 | max(rgn(:,2))>dims(2))
%   myWarnDlg('Must choose line endpoints within image boundaries');
%   return;
% end

% In findLinePoints, if y1 == y2, we draw the horizontal line.
% if x1 == x2 we draw a vertical line.
% otherwise, we sample along the longer direction and find the
% appropriate value along the shorter direction.
finalCoords=[];

for thisLineSegment=1:nPoints-1

y1 = coords(2,thisLineSegment); y2 = coords(2,thisLineSegment+1);
x1 = coords(1,thisLineSegment); x2 = coords(1,thisLineSegment+1);

[thisXList, thisYList] = findLinePoints([x1 y1], [x2 y2]);

if ieNotDefined('newCoords')
    newCoords=[thisXList(:),thisYList(:)];
else
    newCoords=[newCoords;[thisXList(:),thisYList(:)]];
end

end

% Do an (inverse) rotation if necessary
if (strcmp(view.viewType,'Flat'))
    newCoords=(rotateCoords(view,newCoords,1));
end

% Add in the slice index
newCoords=[newCoords,ones(size(newCoords,1),1)*curSlice]';

% Convert coords to canonical frame of reference
newCoords = curOri2CanOri(view,newCoords);

% Merge/remove coordinates
if sgn
  disp('Merging Coords')
  coords = mergeCoords(curCoords,newCoords);
else
  disp('Removing Coords')
  coords = removeCoords(newCoords,curCoords);
end

view.ROIs(view.selectedROI).coords = coords;

return;


text:
function R = slsumstat_blks(data, statfunc, varargin)
%SLSUMSTAT_BLKS Sums up statistics on all blocks for partitioned data
%
% $ Syntax $
%   - R = slsumstat_blks(data, statfunc, ...)
%
% $ Arguments $
%   - data:         the data array of cell array of data array filenames
%   - statfunc:     the function to compute statistics on data array
%   - R:            the sum statistics
%
% $ Description $
%   - R = slsumstat_blks(data, statfunc, ...) computes statistics on
%     all data blocks and sums them by plus. When data is an array,
%     it just invokes statfunc on it and return, if data is a cell array
%     of filenames, it loads data of each block, computes statistics
%     blockwise and sums them up to give the final result.
%
% $ Remarks $
%   - The function is widely applicable for diverse types of computation.
%     The only conditions is that the function values on the whole matrix
%     is equivalent to the sum of function values on all blocks.
%
%   - The function values can be either a scalar or an array of any
%     dimensions, provided that the values produced on all blocks have
%     equal size.
%
% $ History $
%   - Created by Dahua Lin, on Aug 8th, 2006
%

%% parse and verify input arguments

if nargin < 2
    raise_lackinput('slsumstat_blks', 2);
end
if ~isnumeric(data) && ~iscell(data)
    error('sltoolbox:invalidargs', ...
        'data should be either a numeric array or a cell array of strings');
end

%% compute

if isnumeric(data)
    R = feval(statfunc, data, varargin{:});
else
    n = numel(data);

    % first section
    curdata = slreadarray(data{1});
    R = feval(statfunc, curdata, varargin{:});

    % other section
    if n > 1
        for i = 2 : n
            curdata = slreadarray(data{i});
            curR = feval(statfunc, curdata, varargin{:});
            R = R + curR;
        end
    end
end

author: lmthang


text:
function [dLdp,iCpY,L] = mci_discount_deriv (P,M,U,Y)
% Gradient of likelihood for discounting model
% FORMAT [dLdp,iCpY,L] = mci_discount_deriv (P,M,U,Y)
%
% P         parameters
% M         model structure
% U         contains rewards and times
% Y         data
%
% dLdp      gradient of log joint
% iCpY      curvature (Fisher Information)
% L         log joint
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: mci_discount_deriv.m 6548 2015-09-11 12:39:47Z will $

dLdp = spm_diff(M.L,P,M,U,Y,1);
dLdp = full(dLdp(:));

X = spm_diff('mci_discount_act',P,M,U,1);
g = mci_discount_gen (P,M,U);
Lambda=diag(g.*(1-g));
iCpY=X'*Lambda*X;

if nargout > 2
    L=mci_discount_like (P,M,U,Y);
end

author: spm


text:
%%%%%script to generate figure 1 of Invariant Scattering Convolution Networks
%%%%%feb 2012%%%%%%%
%%%%%%Joan Bruna and Stephane Mallat$$$$$$$$$

%%%%%This figure shows the Morlet Wavelet used throughout the numerical
%%%%%experiments of the paper%%%%%%%%

N=256;
%options.wavelet_family='morlet';
%options=configure_wavelet_options(options);
%filters=options.filter_bank_name([N N], options);
foptions.J=2;
foptions.L=6;
soptions.M=2;
[Wop,filters]=wavelet_factory_2d([N N], foptions, soptions);

R=8;

Nc=55;

figure
imagesc(flipud(crop(real(fftshift(ifft2(filters.psi.filter{R}.coefft{1}))),Nc,N/2)))
colormap gray
axis square
set(gca,'YTick',[])
set(gca,'XTick',[])

figure
imagesc(flipud(crop(imag(fftshift(ifft2(filters.psi.filter{R}.coefft{1}))),Nc,N/2)))
colormap gray
axis square
set(gca,'YTick',[])
set(gca,'XTick',[])

figure
imagesc(add_circular_mask(flipud(abs(fftshift(filters.psi.filter{R}.coefft{1}))),0.98,1))
%imagesc((flipud(abs(fftshift(filters.psi{1}{2}{2})))))
cc=colormap(gray);
colormap(1-cc)
axis square
axis off
set(gca,'YTick',[])
set(gca,'XTick',[])

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [Psamp,noise,M] = spm_mci_fixed (mcmc,w,fixed,noise,M,U,Y)
% Group fixed effects estimation
% FORMAT [Psamp,noise,M] = spm_mci_fixed (mcmc,w,fixed,noise,M,U,Y)
%
% mcmc      Sampling parameters
% w(:,n)    Random effects for nth subject
% fixed     [fixed.pE, fixed.pC] prior over fixed effects
% noise     noise structure
% M,U,Y     Model, input, data structures
%
% Psamp     Samples, [maxits x M{1}.Np]
% noise     updated noise model
% M         updated model structures
%
% Uses Langevin Monte Carlo
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

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

try, verbose=mcmc.verbose; catch, verbose=0; end
try, maxits=mcmc.maxits; catch, maxits=64; end
try, plot_int=mcmc.plot_int; catch, plot_int=1; end
try, update_obs_step=mcmc.update_obs_step; catch, update_obs_step=64; end
try, h=mcmc.h; catch, h=0.5; end

assign=mcmc.assign;

% Prior over fixed effects
pE=fixed.pE;
Np=length(pE);
V=eye(Np);
ipC=inv(fixed.pC);
log_prior_t2=spm_logdet(ipC)/2-0.5*Np*log(2*pi);

N=length(M);
for n=1:N
    M{n}.logdet_Ce=spm_logdet(M{n}.Ce);
    M{n}.iCe = inv(M{n}.Ce);
end

try, init=mcmc.init; catch, init=fixed.vpE; end
xinit = init;
x = zeros(maxits,Np);
x(1,:) = xinit';

if verbose figure; end

% Tune h by monitoring acceptance rate
tune_h=1;
acc_block=32;
acc_low=0.3;
acc_high=0.7;
total_acc_target=64; % Number of accepted samples to get
acc=zeros(maxits,1);

i=1;
while (i < maxits) && (sum(acc) < total_acc_target),

    if verbose
        if mod(i,plot_int) == 0 && i > 2
            spm_mci_progress (x,E,i);
        end
    end

    if mod(i,acc_block)==0 && tune_h
        % Change step size h ?
        Nacc=sum(acc(i-acc_block+1:i-1));
        prop_acc=Nacc/acc_block;
        if prop_acc < acc_low
            if verbose, disp('Decreasing step size ...'); end
            h=h/2;
        elseif prop_acc > acc_high
            if verbose, disp('Increasing step size ...'); end
            h=h*2;
        end
    end

    % Proposal (first proposal always accepted)
    if i==1
        pos=x(1,:)';
        curr=[];
    else
        pos=spm_normrnd(curr.mu,curr.Cp,1);
    end

    % Quantities re proposal
    prop.pos=pos;

    % Get gradient and curvature of log like
    dLdp=zeros(1,Np);
    iCpY=zeros(Np,Np);

    for n=1:N,
        if ~isempty(w), wsub=w(:,n); else, wsub=[]; end
        [dLdp_n,iCpY_n,st] = spm_mci_grad_curve (assign,wsub,pos,M{n},U{n},Y{n},'fixed');

        dLdp = dLdp + dLdp_n;
        iCpY = iCpY + iCpY_n;
    end
    dlogprior = - (pos-pE)'*ipC;
    j=dLdp+dlogprior;

    % Log like
    log_like=0;
    for n=1:N,
        if ~isempty(w), wsub=w(:,n); else, wsub=[]; end
        [Pinit,Pflow] = spm_mci_init_flow (assign,wsub,pos,M{n});
        like_n = spm_mci_like_ind (Pflow,Pinit,M{n},U{n},Y{n});
        log_like = log_like + like_n;
    end
    e=pos-pE;
    log_prior = - e'*ipC*e/2 + log_prior_t2;
    prop.L=log_like+log_prior;

    if st==-1
        disp('Integration problem in spm_mci_random.m');
        keyboard
    end

    % Posterior covariance under local linear approximation
    prop.Cp = h*inv(iCpY+ipC);
    prop.mu = pos+0.5*h*prop.Cp*j(:);

    prop.iCp = inv(prop.Cp);
    prop.logdetCp = spm_logdet(prop.Cp);

    % Accept proposal ?
    [curr,accepted,bayes_fb(i),dL(i)] = spm_mci_mh_update(curr,prop,verbose);
    acc(i)=accepted;
    E(i) = -curr.L;
    if i > 1
        dEdit(i-1)=100*(E(i)-E(i-1))/E(i-1);
    end
    x(i,:)=curr.pos;

    i=i+1;
end

Psamp=x(1:i-1,:);


text:
% Display ASCII waiting bar
%
% Inputs
%  max_index : maximum index
%  -- or --
%  pb        : progress bar structure
%  i         : current
%  title     : text shown after the remaining time
%
% Outputs
%  pb        : progress bar structure
%
% Example
%  pb = progressbar(max_index);
%  for index=1:max_index;
%    [...]
%    pb = progressbar(pb, index);
%  end
%
% Copyright (c) 2008 Ircam-CNRS UMR9912-STMS
%
% License
%  This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file is
%  distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%  PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <gilles.degottex@ircam.fr>
%

function pb = progressbar(pbin, i, text)

	if nargin==1
		pb.max_n = pbin;
		pb.history = [];
		pb.i_old = 0;
        pb.finished = false;
        pb.nbp = 0;
		tic;
		return;
	else
		pb = pbin;
	end

	if toc<1/16 && i<pb.max_n; return; end

%  	if nargin~=2
%  		return;
%  	end

    if ~pb.finished

	    pb.history = [pb.history, toc];
	    pb.iteration_time = mean(pb.history);

	    pos = i / pb.max_n;
	    str='|=';
	    for c = 1:floor(pos*50)
		    str = [str, '='];
	    end
	    for c = floor(pos*50):50-1
		    str = [str, ' '];
	    end

	    di = i - pb.i_old;
	    nrest = (pb.max_n-i)/di;
	    minutes = floor(sum(pb.history)/60);
	    seconds = rem(sum(pb.history), 60);
	    str_times = ['  elapsed: ', num2str(minutes), ':', num2str(floor(seconds))];
	    if nrest>0
		    rem_minutes = floor(pb.iteration_time*nrest/60);
		    rem_seconds = rem(pb.iteration_time*nrest, 60);
		    str_times = [str_times, ' remaining: ', num2str(rem_minutes), ':', num2str(round(rem_seconds))];
	    else
            str_times = [str_times, '               '];
        end

		str_full = [num2str(floor(pos * 100)),'%% ', str_times];
		if 1
			str_full = [repmat('\b',1,pb.nbp),str_full];
		else
			str_full = ['\r',str_full];
		end

        if nargin>2
        	str_full = [str_full, ': ', text];
        end
       	str_full = [str_full, '       '];
        nbp = fprintf(str_full);
        pb.nbp = nbp - pb.nbp;

	    pb.i_old = i;
	    tic;

        if i>=pb.max_n;
            pb.finished = true;
            disp(' ');
        end
    end

return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
% Copyright 2019, The MathWorks, Inc.
%% animateRotorPosition
%
% This script will run a Simulink model of a BLDC that is energized in
% one coil and animate the movement of the rotor. Once the animation figure
% is rendered, there is a 5 second pause before the animation begins.

close all
clear

% Below parameters are defined in the Simulink model
% Sample time
% Ts = 2e-5;
% Number of pole pairs
% p = 1;
% Initial rotor angle in degrees
% th0 = 0;
% Sector
% sector = 6;

mdl = 'Modeling_commutation_logic.slx';
open_system(mdl);

sim(mdl)
try

    %
    r = 1.2;
    theta = linspace(0,2*pi);
    x = cos(theta);
    y = sin(theta);

    x1 = 0.8*cos(theta);
    y1 = 0.8*sin(theta);

    xa = cos(0);
    ya = sin(0);

    xb = cos(2*pi/3);
    yb = sin(2*pi/3);

    xc = cos(-2*pi/3);
    yc = sin(-2*pi/3);

    xat = r*cos(0);
    yat = r*sin(0);

    xbt = r*cos(2*pi/3);
    ybt = r*sin(2*pi/3);

    xct = r*cos(-2*pi/3);
    yct = r*sin(-2*pi/3);

    hf = figure(1);h = plot(x,y,'k-',[0 xa],[0 ya],'k-',[0 xb],[0 yb],'k-',[0 xc],[0  yc],'k-',xa,ya,'ko',xb,yb,'ko',xc,yc,'ko',x1,y1,'k-');grid

    set(hf,'Color',[1 1 1])

    ha = gca;

    set(ha,'Visible','off')

    ht1 = text(xat,yat,'A');
    ht2 = text(xbt,ybt,'B');
    ht3 = text(xct,yct,'C');

    set(ht1,'FontSize',14);
    set(ht2,'FontSize',14);
    set(ht3,'FontSize',14);

    axis([-1.2 1.2 -1.2 1.2])
    axis equal

    set(h(5),'MarkerSize',20)
    set(h(6),'MarkerSize',20)
    set(h(7),'MarkerSize',20)

    switch_pattern = switchPattern.signals.values(1:40:end,:);

    init_switch_pattern = num2str(switch_pattern(1,:));

    switch init_switch_pattern

        case num2str([1 0 0 0 0 1])

            set(h(5),'MarkerFaceColor',[1 0 0])
            set(h(7),'MarkerFaceColor',[0 0 1])
            set(h(6),'MarkerFaceColor',[1 1 1])

        case num2str([0 0 1 0 0 1])

            set(h(6),'MarkerFaceColor',[1 0 0])
            set(h(7),'MarkerFaceColor',[0 0 1])
            set(h(5),'MarkerFaceColor',[1 1 1])

        case num2str([0 1 1 0 0 0])

            set(h(5),'MarkerFaceColor',[0 0 1])
            set(h(6),'MarkerFaceColor',[1 0 0])
            set(h(7),'MarkerFaceColor',[1 1 1])

        case num2str([0 1 0 0 1 0])

            set(h(5),'MarkerFaceColor',[0 0 1])
            set(h(7),'MarkerFaceColor',[1 0 0])
            set(h(6),'MarkerFaceColor',[1 1 1])

        case num2str([0 0 0 1 1 0])

            set(h(6),'MarkerFaceColor',[0 0 1])
            set(h(7),'MarkerFaceColor',[1 0 0])
            set(h(5),'MarkerFaceColor',[1 1 1])

        case num2str([1 0 0 1 0 0])

            set(h(5),'MarkerFaceColor',[1 0 0])
            set(h(6),'MarkerFaceColor',[0 0 1])
            set(h(7),'MarkerFaceColor',[1 1 1])

    end

    l1 = 0.7;
    ro = thetaSim.signals.values(1)*pi/180-pi/2;

    xr1 = l1*cos(ro);
    yr1 = l1*sin(ro);
    xr2 = -l1*cos(ro);
    yr2 = -l1*sin(ro);

    hold on,hl1 = plot([xr1 xr2],[yr1 yr2],'m-');

    set(hl1,'LineWidth',3)

    hpr = plot(xr2,yr2,'ko');
    hpb = plot(xr1,yr1,'ko');

    set(hpr,'MarkerSize',20)
    set(hpb,'MarkerSize',20)
    set(hpr,'MarkerFaceColor',[1 0 0])
    set(hpb,'MarkerFaceColor',[0 0 1])

    pause

    ro1 = thetaSim.signals.values(1:40:end)-90;

    for l = 1:numel(ro1)

        switch_pattern1 = num2str(switch_pattern(l,:));
        ro = ro1(l)*pi/180;

        xr1 = l1*cos(ro);
        yr1 = l1*sin(ro);
        xr2 = -l1*cos(ro);
        yr2 = -l1*sin(ro);

        set(hpr,'XData',xr2);
        set(hpr,'YData',yr2);
        set(hpb,'XData',xr1);
        set(hpb,'YData',yr1);

        set(hl1,'XData',[xr1 xr2]);
        set(hl1,'YData',[yr1 yr2]);

        switch switch_pattern1

            case num2str([1 0 0 0 0 1])

                set(h(5),'MarkerFaceColor',[1 0 0])
                set(h(6),'MarkerFaceColor',[1 1 1])
                set(h(7),'MarkerFaceColor',[0 0 1])

            case num2str([0 0 1 0 0 1])

                set(h(6),'MarkerFaceColor',[1 0 0])
                set(h(7),'MarkerFaceColor',[0 0 1])
                set(h(5),'MarkerFaceColor',[1 1 1])

            case num2str([0 1 1 0 0 0])

                set(h(5),'MarkerFaceColor',[0 0 1])
                set(h(6),'MarkerFaceColor',[1 0 0])
                set(h(7),'MarkerFaceColor',[1 1 1])

            case num2str([0 1 0 0 1 0])

                set(h(5),'MarkerFaceColor',[0 0 1])
                set(h(7),'MarkerFaceColor',[1 0 0])
                set(h(6),'MarkerFaceColor',[1 1 1])

            case num2str([0 0 0 1 1 0])

                set(h(6),'MarkerFaceColor',[0 0 1])
                set(h(7),'MarkerFaceColor',[1 0 0])
                set(h(5),'MarkerFaceColor',[1 1 1])

            case num2str([1 0 0 1 0 0])

                set(h(5),'MarkerFaceColor',[1 0 0])
                set(h(6),'MarkerFaceColor',[0 0 1])
                set(h(7),'MarkerFaceColor',[1 1 1])

        end

        drawnow

    end
catch

    disp("The animation was closed before the end of the simulated data")

end

author: mathworks


text:
%% RESTICTED geometry DEFINITION (OMAS_restrictedGeometry.m) %%%%%%%%%%%%%%
% This function computes the geometry that can be observed within a defined
% radius.

function [restrictedGeometry] = OMAS_restrictedGeometry(point,radius,geometry)
% INPUTS:
% centroid - The 3D cartesian reference point.
% radius   - The detection radius of the observational sphere.
% geometry - The geometry of the second object in the reference frame of
%            the centroid.

% TO ENSURE FUNCTIONALITY
sendCompleteGeometryOnDetection = 1;

newSurfaceMembers = [];
% EVALUATE THE geometry'S FACE MEMBERS
for face = 1:size(geometry.faces,1)
    % NOTE:
    % - If all vertices are outside of the sphere, the sphere may still
    %   intersect the plane at an edge or mid-plane.

    % MEMBERS IF THE PLANE
    memberID_A = geometry.faces(face,1);
    memberID_B = geometry.faces(face,2);
    memberID_C = geometry.faces(face,3);
    % TEST THE FACE AGAINST THE SPHERE
    faceIsInsideRadius = OMAS_sphereTriangleIntersection(point,radius,...
                                                         geometry.vertices(memberID_A,:)',...
                                                         geometry.vertices(memberID_B,:)',...
                                                         geometry.vertices(memberID_C,:)');
    if sendCompleteGeometryOnDetection && faceIsInsideRadius
    % IF ANY FACE IS VISIBLE, SEND THE COMPLETE GEOMETRY
        restrictedGeometry = geometry;
        return
    elseif faceIsInsideRadius
    % IF THE FACE IS FOUND TO VIOLATE THE CONSTRAINT
%         fprintf('\t face %d interacts with the sphere\n',face);
        % add all vertices belonging to that face to the matrix
        newSurfaceMembers = vertcat(newSurfaceMembers,memberID_A,memberID_B,memberID_C);
    end
end
% PARSE COMMON POINTS
newSurfaceMembers = unique(newSurfaceMembers,'rows');   % The unique vertices
newVertices = geometry.vertices(newSurfaceMembers,:);   % Extract the common vertices

% THE CONSTRAINT SURFACE
% sphereGeometry = OMAS_graphics.defineSphere(zeros(3,1),radius,6);
% [A B] = SurfaceIntersection(sphereGeometry,geometry);

% TRIANGULATE THE NEW VOLUME
triangulatedVolume = triangulation(boundary(newVertices(:,1),newVertices(:,2),newVertices(:,3),1),newVertices);
restrictedGeometry = geometry;
restrictedGeometry.vertices = triangulatedVolume.Points;
restrictedGeometry.faces = triangulatedVolume.ConnectivityList;
restrictedGeometry.normals = OMAS_graphics.normals(restrictedGeometry);    % Compute the surface normals

% DEBUG
% fig = figure(2);
% ax = axes(fig);
% axis square vis3d;
% % THE OBJECT GEOMETRY
% patch(ax,...
%     'Faces',geometry.faces,...
%     'Vertices',geometry.vertices,...
%     'faceColor','g');
% % THE CONSTRAINT VOLUME
% patch(ax,...
%     'Faces',sphereGeometry.faces,...
%     'Vertices',sphereGeometry.vertices,...
%     'faceColor','b',...
%     'faceAlpha',0.2);
% % THE RESULTING GEOMETRY
% patch(ax,...
%     'Faces',restrictedGeometry.faces,...
%     'Vertices',restrictedGeometry.vertices,...
%     'faceColor','b',...
%     'faceAlpha',0.2);
% end

end

% THE SQUARED DISTANCE TO ALL VERTICES
% verticesSq = (geometry.vertices(:,1) - point(1)).^2 ...
%            + (geometry.vertices(:,2) - point(2)).^2 ...
%            + (geometry.vertices(:,3) - point(3)).^2;
% % VERTEX INSIDE RADIUS LOGICALS
% vertexDetectionLogicals = verticesSq <= radius^2;
% vertexIDset = 1:1:numel(vertexDetectionLogicals);
%
% % VISIBLE POINT ID's
% visibleIDs = vertexIDset(vertexDetectionLogicals);
% if numel(visibleIDs) < 1
%     disp('no vertices are visible');
% end

% [A,~] = ismember(relativeFaces,vertIndices);
% completeFaceIndices = sum(A,2) > 2;                        % Minimum number of connections

repo: OpenMAS


text:
function [MPa] = inH2O2MPa(inH2O)
% Convert pressure from inches of water column at 4 degrees to megapascals
% Chad Greene 2012
MPa = inH2O*0.000249089;

author: Sable

repo: mcbench-benchmarks


text:
function [metabolite_structure] = getMissingHMDBMolForm(metabolite_structure,molFileDirectory,retrievePotHMDB,startSearch,endSearch)
% This function uses getMolFileFromHMDB.m to obtain mol files for those
% metabolite structure entries that have hmdb id's specified. It also has
% the option find novel hmdb entries by querying hmdb for the metabolite
% names. Only perfect matches are considered. See retrievePotHitsHMDB.m for
% more details.
%
% INPUT
% metabolite_structure  metabolite structure
% molFileDirectory      directory where mol files should be stored
% retrievePotHMDB       default: true (attention: this could be time
%                       consuming)
% startSearch           specify where the search should start in the
%                       metabolite structure. Must be numeric (optional, default: all metabolites
%                       in the structure will be search for)
% endSearch             specify where the search should end in the
%                       metabolite structure. Must be numeric (optional, default: all metabolites
%                       in the structure will be search for)
%
% OUTPUT
% metabolite_structure  updated metabolite structure
%
% Ines Thiele 09/21

if ~exist('retrievePotHMDB','var')
    retrievePotHMDB = 1;
end
annotationSource = 'HMDB';
annotationType = 'automatic';

if ~exist('startSearch','var')
    startSearch = 1;
end
if ~exist('endSearch','var')
    endSearch = length(F);
end
%% obtain inchiString, charged formula, mol files for each metabolite
F = fieldnames(metabolite_structure);
for i = startSearch : endSearch
    if retrievePotHMDB == 1
        if isempty(metabolite_structure.(F{i}).hmdb) || ~isempty(find(isnan(metabolite_structure.(F{i}).hmdb)))
            % try to find the hmdb entry based on name
            met = metabolite_structure.(F{i}).metNames;
             % remove abbr with non-breaking space
             if any(ismember(met,char(160)))
                 met=strsplit(met,char(160));
                 met= met{1};
             end
            % remove parenthesis from metabolite names
            met = regexprep(met,'\(','');
            met = regexprep(met,'\)','');
            met = regexprep(met,'\[','');
            met = regexprep(met,'\]','');
            % retrieve potential hits
            try
                hmdb =retrievePotHitsHMDB(met); % the catching should avoid that script fails when connection to HMDB fails
            catch
                hmdb = [];
            end
            if ~isempty(hmdb)
                metabolite_structure.(F{i}).hmdb = hmdb;
                metabolite_structure.(F{i}).hmdb_source = ['Metabolite searched in HMDB by name',':',annotationType,':',datestr(now)];
            end
        end
    end
    % check that hmdb id is defined
    if ~isempty(metabolite_structure.(F{i}).hmdb) && length(find(isnan(metabolite_structure.(F{i}).hmdb)))==0
        % get inchistring only if there is no inchistring for the metabolite
        % defined
        %         if (isempty(metabolite_structure.(F{i}).inchiString)) || length(find(isnan(metabolite_structure.(F{i}).inchiString)))> 0
        %             % input into this function is the HMDB ID
        %             [inchiString] = getInchiStringFromHMDB( metabolite_structure.(F{i}).hmdb);
        %             if ~isempty(inchiString)
        %                 metabolite_structure.(F{i}).inchiString = inchiString;
        %                 metabolite_structure.(F{i}).inchiString_source = [annotationSource,':',annotationType,':',datestr(now)];
        %             end
        %         end
        % get mol files for each entry
        if isempty(metabolite_structure.(F{i}).hasmolfile) || length(find(isnan(metabolite_structure.(F{i}).hasmolfile)))>0 ...
                || strcmp(metabolite_structure.(F{i}).hasmolfile,'0')
            %             [outFile] = getMolFileFromHMDB(metabolite_structure.(F{i}).VMHId, metabolite_structure.(F{i}).hmdb,molFileDirectory);
            %             if ~isempty(outFile)
            %                 metabolite_structure.(F{i}).hasmolfile = num2str(1);
            %                 metabolite_structure.(F{i}).hasmolfile_source = [annotationSource,':',annotationType,':',datestr(now)];
            %
            [metabolite_structure,molCollectionReport] = getMolFilesMultipleSources(metabolite_structure, molFileDirectory,i,i,'hmdb');
            [metabolite_structure,IDsAdded,InchiKeyList,InchiStringList] = generateInchiFromMol(metabolite_structure,molFileDirectory, 1, 1,1,i,i);
            % end
        end

    end
end
[metabolite_structure] = addMetFormulaCharge(metabolite_structure,startSearch,endSearch);

author: opencobra

repo: cobratoolbox


text:
function subaperturedemo(varargin)
%SUBAPERTUREDEMO Read image, generate, and view subaperture-processed image using MATLAB GUIs.
%    subaperturedemo('PropertyName',PropertyValue,...)
%
% Calculates the subaperture-processed image of a complex image (selected
% through a MATLAB dialog box) using the properties specified. The AOI over
% which to compute the subaperture-processed image is selected through a MATLAB
% GUI. This version of subaperture processing does NOT require that the complete
% data fit into memory. It processes from any format handled by OPEN_READER and
% outputs to files in SIO format.
%
%       Property name     Description
%       frames            number of frames (default = 7)
%       apfraction        fraction of aperture for each subaperture
%                            (default = .25)
%       method            'normal' (default), 'fullpixel', or 'minimal'
%       platformdir       platform direction, 'right' (default) or 'left'
%       dim               dimension over which to split subaperture
%                            (default = 1)
%       fill              fill factor
%
% Output frames are stored in one frame per file, where the first frame has
% the filename OUTFILE, and each consecutive frame has the frame number
% appended onto the filename.
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

demo_core(@subaperturefile, varargin{:});

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

author: ngageoint


text:
function x = mtimes( a, x )
    %MINUS Multiplication of TT/MPS block tensor by scalar or vector
    %   X = MTIMES(A, X) multiplies the TT/MPS tensor X by A.
    %       If A is a scalar, all blocks are multiplied by the x.order-th root of A.
    %       if A is a vector of size X.p, all X.p slices of the supercore X.U{x.mu} are multiplied
    %       by the corresponding entry in A.
    %
    %   See also PLUS, MINUS.

    %   TTeMPS Toolbox.
    %   Michael Steinlechner, 2013-2016
    %   Questions and contact: michael.steinlechner@epfl.ch
    %   BSD 2-clause license, see LICENSE.txt

    %x.U{1} = a*x.U{1};

    if length(a) == 1
        c = a^(1/x.order);
        for i = 1:x.order
            x.U{i} = c*x.U{i};
        end
    elseif length(a) == x.p
        for i = 1:x.p
            x.U{x.mu}(:,:,:,i) = a(i) * x.U{x.mu}(:,:,:,i);
        end
    else
        error('Dimension mismatch! Can only multiply block tensor X by scalar (whole tensor) or by a vector of length X.p')
    end
end

author: NicolasBoumal

repo: manopt


text:
function [day] = us2day(us)
% Convert time from microseconds to days.
% Chad Greene 2012
day = us*1.157407407407e-11 ;


text:
function [U,G] = unzip_corners(A)
  % UNZIP_CORNERS Given a triangle mesh where corners of each triangle index
  % different matrices of attributes (e.g. read from an OBJ file), unzip the
  % corners into unique efficiently: attributes become properly vertex valued
  % (usually creating greater than #V but less than #F*3 vertices).
  %
  % Inputs:
  %   A  #F by 3 by #A, typically A = cat(3,F,FTC,FN)
  % Outputs:
  %   U  #U by #A list of indices into each attribute for each unique mesh
  %     vertex: U(v,a) is the attribute index of vertex v in attribute a.
  %   G  #F by 3 list of triangle indices into U
  % Example:
  %   [V,F,TC,FTC] = readOBJ('~/Downloads/kiwis/kiwi.obj');
  %   [U,G] = unzip_corners(cat(3,F,FTC));
  %   % display mesh
  %   tsurf(G,V(U(:,1),:));
  %   % display texture coordinates
  %   tsurf(G,TC(U(:,2),:));
  %

  % C  #F*3 by 2 list of vertex indices and uv indices
  C = reshape(A,[],size(A,3));
  [U,~,J] = unique(C,'rows');
  ss = size(A,2);
  G = J(reshape(1:size(A,1)*ss,[],ss));
end

author: alecjacobson

repo: gptoolbox


text:
display('Testing data')
total_test = cell(5+1,1);
test_folder = 'C:\Users\brend_000\Documents\MATLAB\Machine Learning\Ghost Target\Images\cambridge';
for x = 1:5
    display(['Reading img', num2str(x)])
    img = imread([test_folder,'\img', num2str(x),'.jpg']);
    data = cell(8,8);
    mynxel = nxel(8,img);
    [a,b] = size(mynxel);
    for i = 1:a
        for j = 1:b
            data{i,j} = readImage(mynxel{i,j});
        end
    end
    total_test{1,1} = 'Image';
    total_test{x+1,1} = ['img',num2str(x)];
    total_test{1,2} = 'Feature';
    total_test{x+1,2} = data;

    load([test_folder,'_labels\img',num2str(x),'.mat']);
    total_test{1,3} = 'Label Set';
    total_test{x+1,3} = imgLabel;
end
load('model.mat')
load('pca_model')

[x_test, y_test, accuracy, test_labels, confusion_matrix] = validate(total_test,model,'naive');
save('Results\confusion_matrix_naive.mat','confusion_matrix');

[x_test, y_test, accuracy, test_labels, pca_confusion_matrix] = validate(total_test,pca_model,'pca');
save('Results\confusion_matrix_pca.mat','pca_confusion_matrix');

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [ a, b ] = p13_lim ( dim_num )

%*****************************************************************************80
%
%% P13_LIM returns the integration limits for problem 13.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Output, real A(DIM_NUM), B(DIM_NUM), the lower and upper
%    limits of integration.
%
  a(1:dim_num) = 0.0;
  b(1:dim_num) = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [] = show_map(ax, bounds, map_img_filename)
% plot raster map in figure and fix plot bounds
%
% dependency
%   lat_lon_proportions, File Exchange ID = 32462,
%   (c) 2011 by Jonathan Sullivan
%   http://www.mathworks.com/matlabcentral/fileexchange/32462-correctly-proportion-a-latlon-plot
%
% 2010.11.21 (c) Ioannis Filippidis, jfilippidis@gmail.com
%
% See also PLOT_WAY.

hold(ax, 'on')

% image provided ?
if ~isempty(map_img_filename)
    map_img = imread(map_img_filename);
    image('Parent', ax, 'CData', flipdim(map_img,1),...
          'XData', bounds(1,1:2), 'YData', bounds(2,1:2))
end

plot(ax, [bounds(1,1), bounds(1,1), bounds(1,2), bounds(1,2), bounds(1,1)],...
         [bounds(2,1), bounds(2,2), bounds(2,2), bounds(2,1), bounds(2,1)],...
         'ro-')

xlabel(ax, 'Longitude (^o)')
ylabel(ax, 'Latitude (^o)')
title(ax, 'OpenStreetMap osm file')

axis(ax, 'image')
axis(ax, [bounds(1, :), bounds(2, :) ] )
lat_lon_proportions(ax)

author: johnyf

repo: openstreetmap


text:
function [b a]=get_peak_filter(g,Q,f,Fs)

A=10^(g/40);
w=2*pi*f/Fs;
sn=sin(w);
cs=cos(w);
al=sn/(2*Q);

b=[1+al*A   -2*cs   1-al*A];
a=[1+al/A   -2*cs   1-al/A];

b=b/a(1);
a=a/a(1);


text:
%Voice Based Biometric System
%By Ambavi K. Patel.

function mfcc3=melcep(mel1,powspc1,noc1,fsize1,nwin1)
mfcc1(1:noc1,1:nwin1)=0;
for j=1:nwin1
 n=floor(fsize1/2);
 mfcc1(:,j)=dct(log(mel1*abs(powspc1(1:n,j).^2)));
 %converting the power-spectrum to a melfrequency spectrum &taking the logarithm of that spectrum and by computing its inverse using
 %discrete cosine transform
end;
mfcc2(1:13,1:nwin1)=0;
mfcc2(1:13,1:nwin1)=mfcc1(3:15,1:nwin1);
mfcc3=nanclr(mfcc2);


text:
%%******************************************************************
%% ops:
%%
%%   Z = ops(X,operand,Y,alpha);
%%
%%  INPUT:        X = a matrix or a scalar
%%                    or a CELL ARRAY consisting only of matrices
%%          operand = sym, transpose, triu, tril,
%%                    real, imag, sqrt, abs, max, min, nnz,
%%                    spdiags, ones, zeros, norm, sum, row-norm, blk-norm
%%                    rank1, rank1inv, inv
%%                    +,  -, *, .*,  ./, .^
%%     Y (optional) = a matrix or a scalar
%%                    or a CELL ARRAY consisting only of matrices
%% alpha (optional) = a scalar
%%                    or the variable blk.
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

   function Z = ops(X,operand,Y,alpha);

   spdensity = 0.4;

   if (nargin == 2)
      if strcmp(operand,'sym');
         if ~iscell(X);
            [m,n] = size(X);
            if (m == n);
               Z = (X+X')/2;
            elseif (n == 1);
               Z = X;
            else;
               error('X must be square matrix or a column vector');
            end;
         else
            Z = cell(size(X));
            for p = 1:length(X);
                [m,n] = size(X{p});
                if (m == n);
                   Z{p} = (X{p}+X{p}')/2;
                elseif (n == 1);
                   Z{p} = X{p};
                else;
                   error('X{p} must be square matrix or a column vector');
                end;
            end;
         end;
      elseif strcmp(operand,'sqrt') | strcmp(operand,'abs') | ...
         strcmp(operand,'real') | strcmp(operand,'imag');
         if ~iscell(X);
            eval(['Z = ',operand,'(X);']);
         else;
            Z = cell(size(X));
            for p = 1:length(X);
               eval(['Z{p} = ',operand,'(X{p});']);
            end;
         end;
      elseif strcmp(operand,'max') | strcmp(operand,'min') | ...
             strcmp(operand,'sum');
         if ~iscell(X);
            eval(['Z = ',operand,'(X);']);
         else;
            Z = [];
            for p = 1:length(X);
                eval(['Z = [Z  ',operand,'(X{p})','  ];']);
            end;
         end;
         eval(['Z = ',operand,'(Z);']);
      elseif strcmp(operand,'transpose') | strcmp(operand,'triu') | ...
             strcmp(operand,'tril');
         if ~iscell(X);
            if (size(X,1) == size(X,2));
               eval(['Z = ',operand,'(X);']);
            elseif (size(X,2) == 1);
               eval(['Z = X;']);
            else;
               error('X must be square matrix or a column vector');
            end;
         else
            Z = cell(size(X));
            for p = 1:length(X);
                if (size(X{p},1) == size(X{p},2));
                   eval(['Z{p} = ',operand,'(X{p});']);
                elseif (size(X{p},2) == 1);
                   eval(['Z{p} = X{p};']);
                else;
                   error('X{p} must be square matrix or a column vector');
                end;
            end;
         end;
      elseif strcmp(operand,'norm');
         if ~iscell(X);
            Z = full(sqrt(sum(sum(X.*X))));
         else
            Z = 0;
            for p = 1:length(X); Z = Z + sum(sum(X{p}.*X{p})); end;
            Z = sqrt(Z);
         end;
      elseif strcmp(operand,'blk-norm');
         if ~iscell(X);
            Z = full(sqrt(sum(sum(X.*X))));
         else
            Z = zeros(length(X),1);
            for p = 1:length(X); Z(p) = sum(sum(X{p}.*X{p})); end;
            Z = sqrt(Z);
         end;
      elseif strcmp(operand,'inv');
         if ~iscell(X);
            [m,n] = size(X); n2 = n*n;
            if (m==n)
               Z = inv(X);
               if (nnz(Z) > spdensity*n2)
                  Z = full(Z);
               else
                  Z = sparse(Z);
               end
            elseif (m > 1 & n == 1);
               Z = 1./X;
               if (nnz(Z) > spdensity*n)
                  Z = full(Z);
	       else
	          Z = sparse(Z);
               end
            end
         else
            Z = cell(size(X));
            for p = 1:length(X);
               [m,n] = size(X{p}); n2 = n*n;
               if (m==n)
                  Z{p} = inv(X{p});
                  if (nnz(Z{p}) > spdensity*n2)
                     Z{p} = full(Z{p});
		  else
                     Z{p} = sparse(Z{p});
                  end
               elseif (m > 1 & n == 1);
                  Z{p} = 1./X{p};
                  if (nnz(Z{p}) > spdensity*n)
                     Z{p} = full(Z{p});
		  else
                     Z{p} = sparse(Z{p});
                  end
               end
            end
         end
      elseif strcmp(operand,'getM');
         if ~iscell(X);
            Z = size(X,1);
         else
            for p = 1:length(X); Z(p) = size(X{p},1); end;
            Z = sum(Z);
         end;
      elseif strcmp(operand,'nnz');
         if ~iscell(X);
            Z = nnz(X);
         else;
            for p = 1:length(X);
               Z(p) = nnz(X{p});
            end;
            Z = sum(Z);
         end;
      elseif strcmp(operand,'ones');
         if ~iscell(X);
            Z = ones(size(X));
         else
            Z = cell(size(X));
            for p = 1:length(X);
                Z{p} = ones(size(X{p}));
            end
         end
      elseif strcmp(operand,'zeros');
         if ~iscell(X);
            [m,n] = size(X);
            Z = sparse(m,n);
         else
            Z = cell(size(X));
            for p = 1:length(X);
               [m,n] = size(X{p});
               Z{p} = sparse(m,n);
            end
         end
      elseif strcmp(operand,'identity');
         blk = X;
         Z = cell(size(blk,1),1);
         for p = 1:size(blk,1)
            pblk = blk(p,:); n = sum(pblk{2});
            if strcmp(pblk{1},'s')
               Z{p} = speye(n,n);
            elseif strcmp(pblk{1},'q')
               s = 1+[0, cumsum(pblk{2})];
               len = length(pblk{2});
               Z{p} = zeros(n,1);
               Z{p}(s(1:len)) = ones(len,1);
            elseif strcmp(pblk{1},'l')
               Z{p} = ones(n,1);
            elseif strcmp(pblk{1},'u')
               Z{p} = zeros(n,1);
            end
         end
      elseif strcmp(operand,'row-norm');
         if ~iscell(X);
            if (size(X,2) == size(X,1));
               Z = sqrt(sum((X.*conj(X))'))';
            elseif (size(X,2) == 1);
               Z = abs(X);
            end
         else
            Z = cell(size(X));
            for p = 1:length(X);
               if (size(X{p},2) == size(X{p},1));
                  Z{p} = sqrt(sum((X{p}.*conj(X{p}))'))';
               elseif (size(X{p},2) == 1);
                  Z{p} = abs(X{p});
               end
            end
         end
      end
   end
%%
   if (nargin == 3)
      if strcmp(operand,'spdiags');
         if ~iscell(Y);
            [m,n] = size(Y);
            if (m == n);
               Z = spdiags(X,0,m,n);
            else
               Z = X;
            end
         else
            Z = cell(size(Y));
            for p = 1:length(Y);
               [m,n] = size(Y{p});
               if (m == n);
                  Z{p} = spdiags(X{p},0,m,n);
               else;
                  Z{p} = X{p};
               end
            end
         end
      elseif strcmp(operand,'inprod')
         if ~iscell(X) & ~iscell(Y)
    	    Z = (Y'*X)';
	 elseif iscell(X) & iscell(Y)
    	    Z = zeros(size(X{1},2),1);
   	    for p=1:length(X)
	        Z = Z + (Y{p}'*X{p})';
            end
         end
      elseif strcmp(operand,'+') | strcmp(operand,'-') | ...
             strcmp(operand,'/') | strcmp(operand,'./') | ...
             strcmp(operand,'*') | strcmp(operand,'.*') | ...
             strcmp(operand,'.^');
         if (~iscell(X) & ~iscell(Y));
            eval(['Z = X',operand,'Y;']);
         elseif (iscell(X) & iscell(Y))
            Z = cell(size(X));
            for p = 1:length(X);
 	       if (size(X{p},2) == 1) & (size(Y{p},2) == 1) & ...
                  (strcmp(operand,'*') | strcmp(operand,'/'));
                  eval(['Z{p} = X{p}.',operand,'Y{p};']);
               else
                  eval(['Z{p} = X{p} ',operand,'Y{p};']);
               end
            end
         elseif (iscell(X) & ~iscell(Y));
	    if (length(Y) == 1); Y = Y*ones(length(X),1); end
            Z = cell(size(X));
            for p = 1:length(X);
                eval(['Z{p} = X{p}',operand,'Y(p);']);
            end
         elseif (~iscell(X) & iscell(Y));
            Z = cell(size(Y));
	    if (length(X) == 1); X = X*ones(length(Y),1); end
            for p = 1:length(Y);
                eval(['Z{p} = X(p)',operand,'Y{p};']);
            end
         end
      else
         error([operand,' is not available, check input arguments']);
      end
   end
%%
   if (nargin == 4)
      if strcmp(operand,'rank1') | strcmp(operand,'rank1inv');
         Z = cell(size(alpha,1),1);
         for p = 1:size(alpha,1);
            if ~strcmp(alpha{p,1},'diag');
               blktmp = alpha{p,2};
               if (length(blktmp) == 1);
                  if strcmp(operand,'rank1');
                     Z{p} = (X{p}*Y{p}' + Y{p}*X{p}')/2;
                  else;
                     Z{p} = 2./(X{p}*Y{p}' + Y{p}*X{p}');
                  end
               else
                  Xp = X{p};
                  Yp = Y{p};
                  n = sum(blktmp);
                  Zp = sparse(n,n);
                  s = [0 cumsum(blktmp)];
                  if strcmp(operand,'rank1');
                     for i = 1:length(blktmp)
                        pos = [s(i)+1 : s(i+1)];
                        x = Xp(pos);
                        y = Yp(pos);
                        Zp(pos,pos) = sparse((x*y' + y*x')/2);
                     end;
                     Z{p} = Zp;
                  else
                     for i = 1:length(blktmp)
                        pos = [s(i)+1 : s(i+1)];
                        x = Xp(pos);
                        y = Yp(pos);
                        Zp(pos,pos) = sparse(2./(x*y' + y*x'));
                     end
                     Z{p} = Zp;
                  end
               end
            elseif strcmp(alpha{p,1},'diag');
               if strcmp(operand,'rank1');
                  Z{p} = X{p}.*Y{p};
               else
                  Z{p} = 1./(X{p}.*Y{p});
               end
            end
         end
      elseif strcmp(operand,'+') | strcmp(operand,'-');
         if ~iscell(X) & ~iscell(Y);
            eval(['Z = X',operand,'alpha*Y;']);
         elseif (iscell(X) & iscell(Y));
            Z = cell(size(X));
	    if (length(alpha) == 1);
               alpha = alpha*ones(length(X),1);
            end
            for p = 1:length(X);
               eval(['Z{p} = X{p}',operand,'alpha(p)*Y{p};']);
            end
         else
            error('X, Y are different objects');
         end
      else
         error([operand,' is not available']);
      end
   end
%%============================================================

author: yu-jiang

repo: radpbook

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


text:
%-----------------------------------------------------------------------
% Job saved on 18-Nov-2014 09:28:14 by cfg_util (rev $Rev: 6134 $)
% spm SPM - SPM12 (6225)
% cfg_basicio BasicIO - Unknown
%-----------------------------------------------------------------------
matlabbatch{2}.spm.stats.fmri_est.spmmat(1) = cfg_dep('fMRI model specification: SPM.mat File', substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','spmmat'));
matlabbatch{2}.spm.stats.fmri_est.write_residuals = 0;
matlabbatch{2}.spm.stats.fmri_est.method.Classical = 1;

author: translationalneuromodeling

repo: tapas


text:
function [Gain, errMsg] = bst_duneuro(cfg)
% BST_DUNEURO: Call Duneuro to compute a FEM solution for Brainstorm.
%
% USAGE:  [Gain, errMsg] = bst_duneuro(cfg)

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

% Initialize returned values
Gain = [];
% Install/load duneuro plugin
[isInstalled, errMsg, PlugDesc] = bst_plugin('Install', 'duneuro', cfg.Interactive);
if ~isInstalled
    return;
end
bst_plugin('SetProgressLogo', 'duneuro');
% Get DUNEuro executable
DuneuroExe = bst_fullfile(PlugDesc.Path, PlugDesc.SubFolder, 'bin', ['bst_duneuro_meeg_', bst_get('OsType')]);
if ispc
    DuneuroExe = [DuneuroExe, '.exe'];
else
    DuneuroExe = [DuneuroExe, '.app'];
end

%% ===== SENSORS =====
% Select modality for DUNEuro
isEeg  = strcmpi(cfg.EEGMethod, 'duneuro')  && ~isempty(cfg.iEeg);
isMeg  = strcmpi(cfg.MEGMethod, 'duneuro')  && ~isempty(cfg.iMeg);
isEcog = strcmpi(cfg.ECOGMethod, 'duneuro') && ~isempty(cfg.iEcog);
isSeeg = strcmpi(cfg.SEEGMethod, 'duneuro') && ~isempty(cfg.iSeeg);

% Get the modality
if ((isEeg || isEcog || isSeeg) && isMeg)
    dnModality = 'meeg';
elseif (isEeg || isEcog || isSeeg)
    dnModality = 'eeg';
elseif  isMeg
    dnModality = 'meg'; % from DUNEuro side, EEG, sEEG, ECOG  uses the same process
else
    errMsg = 'No valid modality available.';
    return;
end

% Get EEG positions
% Combined modalities for EEG/sEEG/EcoG as EEG
if (isEeg || isEcog || isSeeg)
    cfg.iEeg = [cfg.iEeg, cfg.iSeeg, cfg.iEcog];
    EegLoc = cat(2, cfg.Channel(cfg.iEeg).Loc);
end

% Get MEG positions/orientations
if isMeg
    MegChannels = [];
    for iChan = 1:length(cfg.iMeg)
        sChan = cfg.Channel(cfg.iMeg(iChan));
        for iInteg = 1:size(sChan.Loc, 2)
            MegChannels = [MegChannels; iChan, sChan.Loc(:,iInteg)', sChan.Orient(:,iInteg)', sChan.Weight(iInteg)];
        end
    end
    % In the case where the MEG integration points are used
    if cfg.UseIntegrationPoint == 0
        % loop over the integration Points
        % chan_loc = figure_3d('GetChannelPositions', cfg, cfg.iMeg); % <= this function is not sufficient, we need also the weights.
        MegChannelsTemp = [];
        for iChan = 1 : length(cfg.iMeg)
            group = MegChannels(MegChannels(:,1) == iChan,:);
            groupPositive = group(group(:,end)>0,:);
            groupNegative = group(group(:,end)<0,:);
            if ~isempty(groupPositive)
                %equivalentPositionPostive = sum(repmat(abs(groupPositive(:,end)),[1 3])  .* groupPositive(:,2:4));
                equivalentPositionPostive = mean(groupPositive(:,2:4));
                MegChannelsTemp = [MegChannelsTemp; iChan  equivalentPositionPostive groupPositive(1,5:7)  sum(groupPositive(:,end))];
            end
            if ~isempty(groupNegative)
                %equivalentPositionNegative = sum(repmat(abs(groupNegative(:,end)),[1 3])  .* groupNegative(:,2:4));
                equivalentPositionNegative = mean(groupNegative(:,2:4));
                MegChannelsTemp = [MegChannelsTemp; iChan  equivalentPositionNegative groupNegative(1,5:7)  sum(groupNegative(:,end))];
            end
        end
        MegChannels = MegChannelsTemp;
    end
end

%% ===== HEAD MODEL =====
% Load FEM mesh
FemMat = load(cfg.FemFile);
% Get mesh type
switch size(FemMat.Elements,2)
    case 4,  ElementType = 'tetrahedron';
    case 8,  ElementType = 'hexahedron';
end
% Remove unselected tissue from the head model (MEG only), this could be also done for sEEG later
if strcmp(dnModality, 'meg')
    % Remove the elements corresponding to the unselected tissues
    iRemove = find(~ismember(FemMat.Tissue, find(cfg.FemSelect)));
    if ~isempty(iRemove)
        FemMat = fem_remove_elem(FemMat, iRemove);
    end
elseif strcmp(dnModality,'meeg') && (sum(cfg.FemSelect) ~= length(unique(FemMat.Tissue)))
    errMsg = 'Reduced head model cannot be used when computing MEG+EEG simultaneously.';
    return;
end
% Hexa mesh: detect whether the geometry was adapted
if strcmpi(ElementType, 'hexahedron')
    GeometryAdapted = [];
    % Detect using the options of the Brainstorm process that created the file
    if isfield(FemMat, 'History') && ~isempty(FemMat.History) && ~isempty(strfind([FemMat.History{:,3}], 'NodeShift'))
        strOptions = [FemMat.History{:,3}];
        iTag = strfind(strOptions, 'NodeShift');
        val = sscanf(strOptions(iTag:end), 'NodeShift=%f');
        if ~isempty(val)
            GeometryAdapted = (val > 0);
        end
    end
    % Otherwise, try to guess based on the geometry
    if isempty(GeometryAdapted)
        % Compute the distance between the first two nodes of each element
        dist = sqrt(sum([FemMat.Vertices(FemMat.Elements(:,1),1) - FemMat.Vertices(FemMat.Elements(:,2),1), ...
         FemMat.Vertices(FemMat.Elements(:,2),2) - FemMat.Vertices(FemMat.Elements(:,2),2), ...
         FemMat.Vertices(FemMat.Elements(:,2),3) - FemMat.Vertices(FemMat.Elements(:,2),3)] .^ 2, 2));
        % If the distance is not constant: then the geomtry is adapted
        GeometryAdapted = (max(abs(dist - dist(1))) > 1e-9);
    end
    % Copy value in DUNEuro options
    if ~isempty(GeometryAdapted)
        cfg.GeometryAdapted = GeometryAdapted;
        disp(['DUNEURO> Detected parameter: GeometryAdapted=', bool2str(GeometryAdapted)]);
    end
end
% Isotropic
if (cfg.Isotropic)
    % Isotropic without tensor
    if ~cfg.UseTensor
        if strcmp(ElementType,'hexahedron')
            MeshFile = 'head_model.dgf';
        else
            MeshFile = 'head_model.msh';
        end
    % Isotropic with tensor
    else
        if strcmp(ElementType,'hexahedron')
            errMsg = 'Using the tensor model with hexahedral mesh is not supported for now.';
            return;
        end
        MeshFile = 'head_model.geo';
    end
% Anisotropic (with tensor)
else
    if strcmp(ElementType,'hexahedron')
        errMsg = 'Using the anisotropy model with hexahedral mesh is not supported for now.';
        return;
    end
    MeshFile = 'head_model.geo';
    cfg.UseTensor = true;
end

% Get temp folder
TmpDir = bst_get('BrainstormTmpDir', 0, 'duneuro');
% Display message
bst_progress('text', 'DUNEuro: Writing temporary files...');
disp(['DUNEURO> Writing temporary files to: ' TmpDir]);

% Write mesh model
MeshFile = fullfile(TmpDir, MeshFile);
out_fem(FemMat, MeshFile);

%% ====== SOURCE SPACE =====
% Source space type
switch (cfg.HeadModelType)
    case 'volume'
        % TODO or keep it as it's now....
    case 'surface'
        bst_progress('text', 'DUNEuro: Fixing source space...');
        % Read cortex file
        sCortex = bst_memory('LoadSurface', cfg.CortexFile);
        cfg.GridLoc = sCortex.Vertices;
        % Shrink the cortex surface by XX mm
        if (cfg.SrcShrink > 0)
            % Get spherical coordinates of the surface normals
            [azimuth, elevation] = cart2sph(sCortex.VertNormals(:,1), sCortex.VertNormals(:,2), sCortex.VertNormals(:,3));
            % Find components to shrink the surface in the three dimensions
            depth = cfg.SrcShrink ./ 1000 .* [cos(elevation) .* cos(azimuth), cos(elevation) .* sin(azimuth), sin(elevation)];
            % Apply to the cortex surface
            cfg.GridLoc = cfg.GridLoc - depth;
        end
        % Force all the dipoles within the GM layer
        iGM = find(panel_duneuro('CheckType', FemMat.TissueLabels, 'gray'), 1);
        iWM = find(panel_duneuro('CheckType', FemMat.TissueLabels, 'white'), 1);
        if cfg.SrcForceInGM && ~isempty(iGM)
            % Install/load iso2mesh plugin
            [isInstalled, errMsg] = bst_plugin('Install', 'iso2mesh', cfg.Interactive);
            if ~isInstalled
                % Delete the temporary files
                file_delete(TmpDir, 1, 1);
                return;
            end
            % Extract GM vertices and elements
            [gmVert, gmElem] = removeisolatednode(FemMat.Vertices, FemMat.Elements(FemMat.Tissue == iGM,:));
            % Compute the centroid of the GM elements
            nElem = size(gmElem, 1);
            nMesh = size(gmElem, 2);
            ElemCenter = zeros(nElem, 3);
            for i = 1:3
                ElemCenter(:,i) = sum(reshape(gmVert(gmElem,i), nElem, nMesh)')' / nMesh;
            end

            % Extract GM envelope
            envFaces = volface(FemMat.Elements(FemMat.Tissue <= iGM,:));
            [envVert, envFaces] = removeisolatednode(FemMat.Vertices, envFaces);
            % Find the dipoles outside of the GM envelope
            iVertOut = find(~inpolyhedron(envFaces, envVert, cfg.GridLoc));
            if ~isempty(iVertOut)
                disp(['DUNEURO> Warning: ' num2str(length(iVertOut)) ' dipole(s) outside of the GM.']);
            end

            % If there is a white matter layer: find the dipoles inside the WM and move them outside to the GM
            if ~isempty(iWM)
                % Extract GM envelope
                wmFaces = volface(FemMat.Elements(FemMat.Tissue == iWM,:));
                [wmVert, wmFaces] = removeisolatednode(FemMat.Vertices, wmFaces);
                % Find the dipoles outside of the GM envelope
                iVertWM = find(inpolyhedron(wmFaces, wmVert, cfg.GridLoc));
                if ~isempty(iVertWM)
                    disp(['DUNEURO> Warning: ' num2str(length(iVertWM)) ' dipole(s) inside the WM.']);
                    iVertOut = union(iVertOut, iVertWM);
                end
            end

            % Move each vertex towards the centroid of the closest GM element
            % view_surface_matrix(cfg.GridLoc, sCortex.Faces)
            for i = 1:length(iVertOut)
                bst_progress('text', sprintf('DUNEuro: Fixing dipole %d/%d...', i, length(iVertOut)));
                % Find the closest GM centroid
                iTarget = dsearchn(ElemCenter, cfg.GridLoc(iVertOut(i),:));
                targetFaces = volface(gmElem(iTarget,:));

                % OPTION #1: Replace the vertex position directly with the centroid.
                % => Problem: might project multiple vertices on the same centroid...
                % cfg.GridLoc(iVertOut ,:) = ElemCenter(iTarget,:);

                % OPTION #2: Gradually move the vertex towards the center of the centroid, until it is located inside the element
                % Move the vertex towards the center until it is inside the element
                %nFix = 10;
                %for iFix = 1:nFix
                %    tmpVert = (nFix - iFix)/nFix * cfg.GridLoc(iVertOut(i),:) + iFix/nFix * ElemCenter(iTarget,:);
                %   if inpolyhedron(targetFaces, gmVert, tmpVert)
                %        distMove = sqrt(sum((cfg.GridLoc(iVertOut(i),:) - tmpVert) .^ 2)) * 1000;
                %        disp(sprintf('DUNEURO> Dipole #%d moved inside the GM (%1.2fmm)', iVertOut(i), distMove));
                %        cfg.GridLoc(iVertOut(i),:) = tmpVert;
                %        break;
                %    end
                %end

                % OPTION #3: move the vertex towards the centroid of the element, and then place the final dipole
                % in the symetric point to the center, as the image of the computed vertex
                % x-----o-----x'
                % ^      ^      ^____ : x' the image of x, or the final dipole position
                % |       |_________ : o is the center of the elem, and middle of [x,x']
                % |_____________ : the point inside the element determined by the tmpVert in the following equation
                nFix = 10; % divid into 10 segments
                iFix = 7; % ratio of the distance tmpVert from the centroid
                tmpVert = (nFix - iFix)/nFix * cfg.GridLoc(iVertOut(i),:) + iFix/nFix * ElemCenter(iTarget,:);
                distPoint = ElemCenter(iTarget,:) - tmpVert;
                newPoint = ElemCenter(iTarget,:) + (distPoint);
                distMove = sqrt(sum((cfg.GridLoc(iVertOut(i),:) - newPoint) .^ 2)) * 1000;
                % use the image/symeric point if it's inside GM
                if inpolyhedron(targetFaces, gmVert, newPoint)
                    disp(sprintf('DUNEURO> iDipole %d/%d : Dipole #%d moved inside the GM (%1.2fmm) (option 3 :as image)',i,length(iVertOut), iVertOut(i), distMove));
                    cfg.GridLoc(iVertOut(i),:) = newPoint;
                elseif inpolyhedron(targetFaces, gmVert, tmpVert) % use the original point if it's inside GM
                    disp(sprintf('DUNEURO> iDipole %d/%d : Warning Dipole #%d moved outside the GM (%1.2fmm) (option 3 :as image)', i,length(iVertOut),iVertOut(i), distMove));
                    % use the original distance unstead of the image
                    distMove = sqrt(sum((cfg.GridLoc(iVertOut(i),:) - tmpVert) .^ 2)) * 1000;
                    disp(sprintf('DUNEURO> iDipole %d/%d : Correction 1: Dipole #%d moved inside the GM (%1.2fmm) (option 3: not image)', i,length(iVertOut),iVertOut(i), distMove));
                    cfg.GridLoc(iVertOut(i),:) = tmpVert;
                else % Use the option 2 defined by Francois
                    disp(sprintf('DUNEURO> iDipole %d/%d : Warning Dipole #%d moved outside the GM (%1.2fmm) (option 3 :as image)', i,length(iVertOut),iVertOut(i), distMove));
                    nFix = 20; % with 10 it's not working for some extrem case, then I upgrade it to 20
                    for iFix = 1: nFix
                       tmpVert = (nFix - iFix)/nFix * cfg.GridLoc(iVertOut(i),:) + iFix/nFix * ElemCenter(iTarget,:);
                      if inpolyhedron(targetFaces, gmVert, tmpVert)
                           distMove = sqrt(sum((cfg.GridLoc(iVertOut(i),:) - tmpVert) .^ 2)) * 1000;
                           disp(sprintf('DUNEURO> iDipole %d/%d : Correction 2: Dipole #%d moved inside the GM (%1.2fmm) (option2)', i,length(iVertOut),iVertOut(i), distMove));
                           cfg.GridLoc(iVertOut(i),:) = tmpVert;
                           break;
                       end
                    end
                end
            end
            % view_surface_matrix(cfg.GridLoc, sCortex.Faces)

%             %%%% =============================================
%             % Now similar process for the WM with an extra and unexpedted step ...
%             if ~isempty(iWM)
%                 disp('Checking the the WM ...');
%                 %Extract WM surface
%                 wm_tetra = FemMat.Elements(FemMat.Tissue <= iWM,:);
%                 wm_face = volface(wm_tetra);
%                 [nwmf, ewmf] = removeisolatednode(FemMat.Vertices,wm_face);
%                 % check if any dipoles is inside the WM
%                 wMfv.vertices = nwmf;
%                 wMfv.faces = ewmf;
%                 tic
%                 wMin = inpolyhedron(wMfv, sCortex.Vertices);
%                 wMindex_in = find(wMin);
%                 disp(['There are ' num2str(sum(wMin)) ' dipoles inside the WM']);
%                 disp('Moving these dipoles to the GM tissues ...');
%                 twm = toc;
%                 if ~isempty(wMindex_in)
%                     % 1- move the dipole from inside the WM to the GM surface
%                     % ==> this is for testing, when we use the centroide
%                     % directely, some dipole remains within the WM ...
%                     GMcentroide = 0; % just for testing, to use directely the GM centroides
%                     if GMcentroide == 1
%                         k = dsearchn(elem_centroide,sCortex.Vertices(wMindex_in,:));
%                         NewVertices(wMindex_in ,:) = ElemCenter(k,:);
%                         wMoutFinal = inpolyhedron(wMfv, NewVertices);
%                         disp(['Now, there are ' num2str(sum(wMoutFinal)) ' dipoles inside  the WM']);
%                     else
%                         k = dsearchn(gMfv.vertices,sCortex.Vertices(wMindex_in,:));
%                         NewVertices(wMindex_in ,:) = gMfv.vertices(k,:);
%                         wMoutFinal = inpolyhedron(wMfv, NewVertices);
%                         disp(['There are ' num2str(sum(wMoutFinal)) ' dipoles inside  the WM']);
%                         disp(['These dipoles are moved to the GM outer surface nodes']);
%                         disp(['Moving these dipoles to the nearest centroide of the GM element...']);
%                         % Now, we move these same dipole from the surface to the
%                         % nearest centroide ==> to be sure it's inside the GM
%                         % ==> fill partially the FEM condition
%                         % move again  to the centoide
%                         k = dsearchn(ElemCenter,NewVertices(wMindex_in,:));
%                         NewVertices(wMindex_in ,:) = ElemCenter(k,:);
%                         % just to check
%                         wMoutFinal = inpolyhedron(wMfv, NewVertices);
%                         disp(['Now, there are ' num2str(sum(wMoutFinal)) ' dipoles inside  the WM']);
%                         disp(['All the dipoles are moved to nearest centroid of the GM']);
%                     end
%                 end
%             end

        end
    case 'mixed'
        % TODO : not used ?
end

%% ===== SOURCE MODEL =====
bst_progress('text', 'DUNEuro: Writing temporary files...');
% Write the source/dipole file
DipoleFile = fullfile(TmpDir, 'dipole_model.txt');
% Unconstrained orientation for each dipole. ie the output file have 3*N dipoles (3 directions for each)
%     [x1,y1,z1, 1, 0, 0,;
%      x1,y1,z1, 0, 1, 0;
%      x1,y1,z1, 0, 0, 1;
%      x2,y2,z2, 1, 0, 0;
%      .... ];
dipoles = [kron(cfg.GridLoc, ones(3,1)), kron(ones(size(cfg.GridLoc,1), 1), eye(3))];
fid = fopen(DipoleFile, 'wt+');
fprintf(fid, '%d %d %d %d %d %d \n', dipoles');
fclose(fid);

%% ===== SENSOR MODEL =====
% Write the EEG electrode file
ElecFile = 'electrode_model.txt';
if isEeg || isEcog || isSeeg
    fid = fopen(fullfile(TmpDir, ElecFile), 'wt+');
    fprintf(fid, '%d %d %d  \n', EegLoc);
    fclose(fid);
end
% Write the MEG sensors file
CoilFile = fullfile(TmpDir, 'coil_model.txt');
ProjFile = fullfile(TmpDir, 'projection_model.txt');
if isMeg
    % Write coil file
    CoilsLoc = MegChannels(:,2:4);
    fid = fopen(CoilFile, 'wt+');
    fprintf(fid, '%d %d %d  \n', CoilsLoc');
    fclose(fid);
    % Write projection file
    CoilsOrient = MegChannels(:,5:7);
    fid = fopen(ProjFile, 'wt+');
    fprintf(fid, '%d %d %d  \n', CoilsOrient');
    fclose(fid);
end

%% ===== CONDUCTIVITY MODEL =====
% Isotropic without tensor
if ~cfg.UseTensor
    CondFile = fullfile(TmpDir, 'conductivity_model.con');
    fid = fopen(CondFile, 'w');
    fprintf(fid, '%d\t', cfg.FemCond);
    fclose(fid);
% With tensor (isotropic or anisotropic)
else
    CondFile = fullfile(TmpDir, 'conductivity_model.knw');
    % Transformation matrix  and tensor mapping on each direction
    CondTensor = zeros(length(FemMat.Elements),6) ;
    for ind =1 : length(FemMat.Elements)
        temp0 = reshape(FemMat.Tensors(ind,:),3,[]);
        T1 = temp0(:,1:3); % get the 3 eigen vectors
        l =  diag(temp0(:,4)); % get the eigen value as 3x3
        temp = T1 * l * T1'; % reconstruct the tensors
        CondTensor(ind,:) = [temp(1) temp(5) temp(9) temp(4) temp(8) temp(7)]; % this is the right order
    end
    % write the tensors
    out_fem_knw(FemMat, CondTensor, CondFile);
end

%% ===== WRITE MINI FILE =====
% Open the mini file
IniFile = fullfile(TmpDir, 'duneuro_minifile.mini');
fid = fopen(IniFile, 'wt+');
% General setting
fprintf(fid, '__name = %s\n\n', IniFile);
if strcmp(cfg.SolverType, 'cg')
    fprintf(fid, 'type = %s\n', cfg.FemType);
end
fprintf(fid, 'element_type = %s\n', ElementType);
fprintf(fid, 'solver_type = %s\n', cfg.SolverType);
fprintf(fid, 'geometry_adapted = %s\n', bool2str(cfg.GeometryAdapted));
fprintf(fid, 'tolerance = %d\n', cfg.Tolerance);
% [electrodes]
if isEcog || isSeeg
    % Instead of selecting the electrode on the outer surface,
    % uses the nearest FEM node as the electrode location
    cfg.ElecType = 'closest_subentity_center';
end
if strcmp(dnModality, 'eeg') || strcmp(dnModality, 'meeg')
    fprintf(fid, '[electrodes]\n');
    fprintf(fid, 'filename = %s\n', fullfile(TmpDir, ElecFile));
    fprintf(fid, 'type = %s\n', cfg.ElecType);
    fprintf(fid, 'codims = %s\n', '3');
end
% [meg]
if strcmp(dnModality, 'meg') || strcmp(dnModality, 'meeg')
    fprintf(fid, '[meg]\n');
    fprintf(fid, 'intorderadd = %d\n', cfg.MegIntorderadd);
    fprintf(fid, 'type = %s\n', cfg.MegType);
    fprintf(fid, 'cache.enable = %s\n',bool2str(cfg.EnableCacheMemory) );
    % [coils]
    fprintf(fid, '[coils]\n');
    fprintf(fid, 'filename = %s\n', CoilFile);
    % [projections]
    fprintf(fid, '[projections]\n');
    fprintf(fid, 'filename = %s\n', ProjFile);
end
% [dipoles]
fprintf(fid, '[dipoles]\n');
fprintf(fid, 'filename = %s\n', DipoleFile);
% [volume_conductor.grid]
fprintf(fid, '[volume_conductor.grid]\n');
fprintf(fid, 'filename = %s\n', MeshFile);
% [volume_conductor.tensors]
fprintf(fid, '[volume_conductor.tensors]\n');
fprintf(fid, 'filename = %s\n', CondFile);
% [solver]
fprintf(fid, '[solver]\n');
fprintf(fid, 'solver_type = %s\n', cfg.SolvSolverType);
fprintf(fid, 'preconditioner_type = %s\n', cfg.SolvPrecond);
if strcmp(cfg.SolverType, 'cg')
    fprintf(fid, 'cg_smoother_type = %s\n', cfg.SolvSmootherType);
end
fprintf(fid, 'intorderadd = %d\n', cfg.SolvIntorderadd);
% Discontinuous Galerkin
if strcmp(cfg.SolverType, 'dg')
    fprintf(fid, 'dg_smoother_type = %s\n', cfg.DgSmootherType);
    fprintf(fid, 'scheme = %s\n', cfg.DgScheme);
    fprintf(fid, 'penalty = %d\n', cfg.DgPenalty);
    fprintf(fid, 'edge_norm_type = %s\n', cfg.DgEdgeNormType);
    fprintf(fid, 'weights = %s\n', bool2str(cfg.DgWeights));
    fprintf(fid, 'reduction = %s\n', bool2str(cfg.DgReduction));
end
% [solution]
fprintf(fid, '[solution]\n');
fprintf(fid, 'post_process = %s\n', bool2str(cfg.SolPostProcess)); % true/false
fprintf(fid, 'subtract_mean = %s\n', bool2str(cfg.SolSubstractMean)); % boolean
% [solution.solver]
fprintf(fid, '[solution.solver]\n');
fprintf(fid, 'reduction = %d\n', cfg.SolSolverReduction);
% [solution.source_model]
fprintf(fid, '[solution.source_model]\n');
fprintf(fid, 'type = %s\n', cfg.SrcModel);
fprintf(fid, 'intorderadd = %d\n', cfg.SrcIntorderadd);
fprintf(fid, 'intorderadd_lb = %d\n', cfg.SrcIntorderadd_lb);
fprintf(fid, 'numberOfMoments = %d\n', cfg.SrcNbMoments);
fprintf(fid, 'referenceLength = %d\n', cfg.SrcRefLen);
fprintf(fid, 'weightingExponent = %d\n', cfg.SrcWeightExp);
fprintf(fid, 'relaxationFactor = %e\n', 10^(-cfg.SrcRelaxFactor));
fprintf(fid, 'mixedMoments = %s\n', bool2str(cfg.SrcMixedMoments));
fprintf(fid, 'restrict = %s\n', bool2str(cfg.SrcRestrict));
fprintf(fid, 'initialization = %s\n', cfg.SrcInit);
% [brainstorm]
fprintf(fid, '[brainstorm]\n');
fprintf(fid, 'modality = %s\n', dnModality);
fprintf(fid, 'output_folder = %s\n', [TmpDir, filesep]);
fprintf(fid, 'save_eeg_transfer_file = %s\n', bool2str(cfg.BstSaveTransfer));
fprintf(fid, 'save_meg_transfer_file = %s\n', bool2str(cfg.BstSaveTransfer));
fprintf(fid, 'save_meeg_transfer_file = %s\n', bool2str(cfg.BstSaveTransfer));
fprintf(fid, 'eeg_transfer_filename = %s\n', cfg.BstEegTransferFile);
fprintf(fid, 'meg_transfer_filename = %s\n', cfg.BstMegTransferFile);
fprintf(fid, 'eeg_leadfield_filename = %s\n', cfg.BstEegLfFile);
fprintf(fid, 'meg_leadfield_filename = %s\n', cfg.BstMegLfFile);
% Close file
fclose(fid);

%% ===== RUN DUNEURO ======
% Assemble command line
callStr = ['"' DuneuroExe '"' ' ' '"' IniFile '"'];
bst_progress('text', 'DUNEuro: Computing leadfield...');
disp(['DUNEURO> System call: ' callStr]);
tic;
% Call DUNEuro
status = system(callStr)
if (status ~= 0)
    errMsg = 'Error during the DUNEuro computation, see logs in the command window.';
    return;
end
disp(['DUNEURO> FEM computation completed in: ' num2str(toc) 's']);

%% ===== READ LEADFIELD ======
bst_progress('text', 'DUNEuro: Reading leadfield...');
% EEG
if (isEeg || isEcog || isSeeg)
    GainEeg = in_duneuro_bin(fullfile(TmpDir, cfg.BstEegLfFile))';
end

%MEG
if isMeg
    GainMeg = in_duneuro_bin(fullfile(TmpDir, cfg.BstMegLfFile))';

    % === POST-PROCESS MEG LEADFIELD ===
    % Compute the total magnetic field
    dipoles_pos_orie = [kron(cfg.GridLoc,ones(3,1)), kron(ones(length(cfg.GridLoc),1), eye(3))];

    % a- Compute the MEG Primary Magnetic B-field analytically (formula of Sarvas)
    dip_pos = dipoles_pos_orie(:,1:3);
    dip_mom = dipoles_pos_orie(:,4:6);
    Bp = zeros(size(MegChannels,1), size(dip_pos,1));
    for i = 1:size(CoilsLoc,1)
        for j = 1 : size(dip_pos,1)
            R = CoilsLoc(i,:);
            R_0 = dip_pos(j,:);
            A = R - R_0;
            a = norm(A);
            aa = A./(a^3);
            BpHelp = cross(dip_mom(j,:),aa);
            Bp(i,j) = BpHelp * CoilsOrient(i, :)'; % projection of the primary B-field along the coil orientations
        end
    end

    % b- The total magnetic field B = Bp + Bs;
    %  full B-field
    Bs = GainMeg;
    mu = 4*pi*1e-4; % check the value of the units maybe it needs to be mu = 4*pi*1e-7
    Bfull = (mu/(4*pi)) * (Bp - Bs);

    % c- Apply the weight :
    [channelIndex] = unique(MegChannels(:,1));
    nbChannel = length(channelIndex);
    weighted_B = zeros(nbChannel,size(Bfull,2));
    for iCh = 1 : nbChannel
        communChannel = find(iCh==MegChannels(:,1));
        BcommunChannel = Bfull(communChannel(:),:);
        WcommunChannel =  MegChannels(communChannel(:), 8: end);
        weighted_B(iCh,:) = sum (BcommunChannel.*WcommunChannel,1);
    end
    GainMeg = weighted_B;
end

% Fill the unused channels with NaN
Gain = NaN * zeros(length(cfg.Channel), 3 * length(cfg.GridLoc));
if isMeg
    Gain(cfg.iMeg,:) = GainMeg;
    % scaling the MEG Gain matrix
    Gain(cfg.iMeg,:) = GainMeg/1000;
end
if (isEeg || isEcog || isSeeg)
    Gain(cfg.iEeg,:) = GainEeg;
end

% Delete the temporary files
file_delete(TmpDir, 1, 1);

%% ===== SAVE TRANSFER MATRIX ======
disp('DUNEURO> TODO: Save transferOut.dat to database.')

% Remove logo
bst_plugin('SetProgressLogo', []);

end

%% =================================================================================
%  === SUPPORT FUNCTIONS  =========================================================
%  =================================================================================

%% ===== BOOL => STR =====
function str = bool2str(bool)
    if bool
        str = 'true';
    else
        str = 'false';
    end
end

author: brainstorm-tools


text:
function value = r4_bide ( x )

%*****************************************************************************80
%
%% R4_BIDE: exponentially scaled derivative, Airy function Bi of an R4 argument.
%
%  Discussion:
%
%    if X < 0,
%      R4_BIDE ( X ) = R4_BID ( X )
%    else
%      R4_BIDE ( X ) = R4_BID ( X ) * exp ( - 2/3 * X^(3/2) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 October 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the exponentially scaled derivative of
%    the Airy function Bi of X.
%
  persistent atr
  persistent bif2cs
  persistent bifcs
  persistent big2cs
  persistent bigcs
  persistent bip1cs
  persistent bip2cs
  persistent btr
  persistent nbif
  persistent nbif2
  persistent nbig
  persistent nbig2
  persistent nbip1
  persistent nbip2
  persistent x2sml
  persistent x3sml
  persistent x32sml
  persistent xbig

  atr =  8.7506905708484345;
  btr = -2.0938363213560543;

  if ( isempty ( nbif ) )

    bif2cs = [ ...
       0.323493987603522033521E+00, ...
       0.086297871535563559139E+00, ...
       0.002994025552655397426E+00, ...
       0.000051430528364661637E+00, ...
       0.000000525840250036811E+00, ...
       0.000000003561751373958E+00, ...
       0.000000000017146864007E+00, ...
       0.000000000000061663520E+00, ...
       0.000000000000000171911E+00, ...
       0.000000000000000000382E+00 ]';
    bifcs = [ ...
       0.1153536790828570243E+00, ...
       0.0205007894049192875E+00, ...
       0.0002135290278902876E+00, ...
       0.0000010783960614677E+00, ...
       0.0000000032094708833E+00, ...
       0.0000000000062930407E+00, ...
       0.0000000000000087403E+00, ...
       0.0000000000000000090E+00 ]';
    big2cs = [ ...
       1.6062999463621294578E+00, ...
       0.7449088819876088652E+00, ...
       0.0470138738610277380E+00, ...
       0.0012284422062548239E+00, ...
       0.0000173222412256624E+00, ...
       0.0000001521901652368E+00, ...
       0.0000000009113560249E+00, ...
       0.0000000000039547918E+00, ...
       0.0000000000000130017E+00, ...
       0.0000000000000000335E+00 ]';
    bigcs = [ ...
      -0.097196440416443537390E+00, ...
       0.149503576843167066571E+00, ...
       0.003113525387121326042E+00, ...
       0.000024708570579821297E+00, ...
       0.000000102949627731379E+00, ...
       0.000000000263970373987E+00, ...
       0.000000000000458279271E+00, ...
       0.000000000000000574283E+00, ...
       0.000000000000000000544E+00 ]';
    bip1cs = [ ...
      -0.1729187351079553719E+00, ...
      -0.0149358492984694364E+00, ...
      -0.0005471104951678566E+00, ...
       0.0001537966292958408E+00, ...
       0.0000154353476192179E+00, ...
      -0.0000065434113851906E+00, ...
       0.0000003728082407879E+00, ...
       0.0000002072078388189E+00, ...
      -0.0000000658173336470E+00, ...
       0.0000000074926746354E+00, ...
       0.0000000011101336884E+00, ...
      -0.0000000007265140553E+00, ...
       0.0000000001782723560E+00, ...
      -0.0000000000217346352E+00, ...
      -0.0000000000020302035E+00, ...
       0.0000000000019311827E+00, ...
      -0.0000000000006044953E+00, ...
       0.0000000000001209450E+00, ...
      -0.0000000000000125109E+00, ...
      -0.0000000000000019917E+00, ...
       0.0000000000000015154E+00, ...
      -0.0000000000000004977E+00, ...
       0.0000000000000001155E+00, ...
      -0.0000000000000000186E+00 ]';
    bip2cs = [ ...
      -0.13269705443526630495E+00, ...
      -0.00568443626045977481E+00, ...
      -0.00015643601119611610E+00, ...
      -0.00001136737203679562E+00, ...
      -0.00000143464350991284E+00, ...
      -0.00000018098531185164E+00, ...
       0.00000000926177343611E+00, ...
       0.00000001710005490721E+00, ...
       0.00000000476698163504E+00, ...
      -0.00000000035195022023E+00, ...
      -0.00000000058890614316E+00, ...
      -0.00000000006678499608E+00, ...
       0.00000000006395565102E+00, ...
       0.00000000001554529427E+00, ...
      -0.00000000000792397000E+00, ...
      -0.00000000000258326243E+00, ...
       0.00000000000121655048E+00, ...
       0.00000000000038707207E+00, ...
      -0.00000000000022487045E+00, ...
      -0.00000000000004953477E+00, ...
       0.00000000000004563782E+00, ...
       0.00000000000000332998E+00, ...
      -0.00000000000000921750E+00, ...
       0.00000000000000094157E+00, ...
       0.00000000000000167154E+00, ...
      -0.00000000000000055134E+00, ...
      -0.00000000000000022369E+00, ...
       0.00000000000000017487E+00, ...
       0.00000000000000000207E+00 ]';

    eta = 0.1 * r4_mach ( 3 );
    nbif = r4_inits ( bifcs, 8, eta );
    nbig = r4_inits ( bigcs, 9, eta );
    nbif2 = r4_inits ( bif2cs, 10, eta );
    nbig2 = r4_inits ( big2cs, 10, eta );
    nbip1 = r4_inits ( bip1cs, 24, eta );
    nbip2 = r4_inits ( bip2cs, 29, eta );
    x2sml = sqrt ( eta );
    x3sml = eta^0.3333;
    x32sml = 1.3104 * x3sml * x3sml;
    xbig = r4_mach ( 2 )^0.6666;

  end

  if ( x <= - 1.0 )
    [ xn, phi ] = r4_admp ( x );
    value = xn * sin ( phi );
  elseif ( 0.0 <= x && x <= x32sml )
    x2 = 0.0;
    x3 = 0.0;
    value = x2 * ( r4_csevl ( x3, bifcs, nbif ) + 0.25 ) ...
      + r4_csevl ( x3, bigcs, nbig ) + 0.5;
  elseif ( abs ( x ) <= x2sml )
    x2 = 0.0;
    x3 = 0.0;
    value = x2 * ( r4_csevl ( x3, bifcs, nbif ) + 0.25 ) ...
      + r4_csevl ( x3, bigcs, nbig ) + 0.5;
    value = value * exp ( - 2.0 * x * sqrt ( x ) / 3.0 );
  elseif ( x <= x3sml )
    x2 = x * x;
    x3 = 0.0;
    value = x2 * ( r4_csevl ( x3, bifcs, nbif ) + 0.25 ) ...
      + r4_csevl ( x3, bigcs, nbig ) + 0.5;
    value = value * exp ( - 2.0 * x * sqrt ( x ) / 3.0 );
  elseif ( x <= 1.0 )
    x2 = x * x;
    x3 = x * x * x;
    value = x2 * ( r4_csevl ( x3, bifcs, nbif ) + 0.25 ) ...
      + r4_csevl ( x3, bigcs, nbig ) + 0.5;
    value = value * exp ( - 2.0 * x * sqrt ( x ) / 3.0 );
  elseif ( x <= 2.0 )
    z = ( 2.0 * x * x * x - 9.0 ) / 7.0;
    value = exp ( - 2.0 * x * sqrt ( x ) / 3.0 ) ...
      * ( x * x * ( 0.25 + r4_csevl ( z, bif2cs, nbif2 ) ) ...
      + 0.5 + r4_csevl ( z, big2cs, nbig2 ) );
  elseif ( x <= 4.0 )
    sqrtx = sqrt ( x );
    z = atr / ( x * sqrtx ) + btr;
    value = ( 0.625 ...
      + r4_csevl ( z, bip1cs, nbip1 ) ) * sqrt ( sqrtx );
  elseif ( x < xbig )
    sqrtx = sqrt ( x );
    z = 16.0 / ( x * sqrtx ) - 1.0;
    value = ( 0.625 + r4_csevl ( z, bip2cs, nbip2 ) ) ...
      * sqrt ( sqrtx );
  else
    sqrtx = sqrt ( x );
    z = - 1.0;
    value = ( 0.625 + r4_csevl ( z, bip2cs, nbip2 ) ) * sqrt ( sqrtx );
  end

  return
end


text:
function [nV, eV, aV, alt, lat, long, type, id] = DecodeBits_ADI(bits, currentLat, currentLong)
% Copyright 2015, The MathWorks, Inc.

% Read message bits and decode valid messages for position, velocity and
% altitude data

% Initialize data
nV = 0;
eV = 0;
aV = 0;
alt = 0;
lat = 0;
long = 0;
type = 'X';
id = [bits(9:12)'*[8;4;2;1] bits(13:16)'*[8;4;2;1] bits(17:20)'*[8;4;2;1] bits(21:24)'*[8;4;2;1] bits(25:28)'*[8;4;2;1] bits(29:32)'*[8;4;2;1]]

% Check 9th and 10th hex characters for mesasge type
tf1 = bits(33:36)'*[8;4;2;1];
tf2 = bits(37:40)'*[8;4;2;1];
if tf1 == 9 && tf2 == 9
    [nV, eV, aV] = AltVelCalc_ADI(bits);
    type = 'A';
elseif tf1 == 5 || tf1 == 6
    [lat, long, alt] = LatLongCalcSingle_ADI(bits, currentLat, currentLong);
    type = 'L';
elseif tf1 == 9 && tf2 == 0
    [lat, long, alt] = LatLongCalcSingle_ADI(bits, currentLat, currentLong);
    type = 'L';
end

author: analogdevicesinc