text:
function [err,time,solver,eqn] = femPoisson3(mesh,pde,option,varargin)
%% FEMPOISSON3 solve Poisson equation by various finite element methods
%
%   FEMPOISSON3 computes approximations to the Poisson equation on a
%   sequence of meshes obtained by uniform refinement of a input mesh.
%
% See also Poisson, crack, Lshape
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

%% Check input arguments
if isfield(mesh,'node') && isfield(mesh,'elem')
    node = mesh.node;
    elem = double(mesh.elem);
else
    [node,elem] = cubemesh([0,1,0,1,0,1],0.25); % default mesh is a cube
end
if isfield(mesh,'bdFlag')
    bdFlag = mesh.bdFlag;
else
    bdFlag = setboundary3(node,elem,'Dirichlet');
end
if ~exist('option','var'), option = []; end
if ~exist('pde','var')
    pde = sincosdata3;                          % default data
end

%% Parameters
option = femoption(option);
maxIt = option.maxIt;
maxN = option.maxN;
L0 = option.L0;
elemType = option.elemType;
refType = option.refType;

%% Generate an initial mesh
for k = 1:L0
    if strcmp(refType,'red')
        [node,elem,bdFlag] = uniformrefine3(node,elem,bdFlag);
    elseif strcmp(refType,'bisect')
        [node,elem,bdFlag] = uniformbisect3(node,elem,bdFlag);
    end
end

%% Initialize err
errL2 = zeros(maxIt,1);   errH1 = zeros(maxIt,1);
erruIuh = zeros(maxIt,1); errMax = zeros(maxIt,1);
errTime = zeros(maxIt,1); solverTime = zeros(maxIt,1);
assembleTime = zeros(maxIt,1); meshTime = zeros(maxIt,1);
itStep = zeros(maxIt,1);  stopErr = zeros(maxIt,1); flag = zeros(maxIt,1);
N = zeros(maxIt,1); h = zeros(maxIt,1);

%% Finite Element Method
for k = 1:maxIt
%     bdFlag = sparse(double(bdFlag));
    % solve the equation
    switch elemType
        case 'P1'     % piecewise linear function P1 element
            [soln,eqn,info] = Poisson3(node,elem,bdFlag,pde,option);
        case 'CR'     % piecewise linear function CR element
            [soln,eqn,info] = Poisson3CR(node,elem,bdFlag,pde,option);
        case 'P2'     % piecewise quadratic function
            [soln,eqn,info] = Poisson3P2(node,elem,bdFlag,pde,option);
        case 'WG'     % weak Galerkin element
            [soln,eqn,info] = Poisson3WG(node,elem,bdFlag,pde,option);
    end
    uh = soln.u;
    % compute error
    t = cputime;
    if isfield(pde,'Du')
        if isfield(soln,'Du') && ~isempty(soln.Du) % Du is in the output
            errH1(k) = getH1error3(node,elem,pde.Du,soln.Du);
        else
            errH1(k) = getH1error3(node,elem,pde.Du,uh);
        end
    end
    if isfield(pde,'exactu')
        errL2(k) = getL2error3(node,elem,pde.exactu,uh);
        % interpolation
        switch elemType
            case 'P1'
                uI = Lagrangeinterpolate(pde.exactu,node,elem);
            case 'CR'
                uI = Lagrangeinterpolate(pde.exactu,node,elem,'CR',eqn.face);
            case 'P2'
                uI = Lagrangeinterpolate(pde.exactu,node,elem,'P2',eqn.edge);
            case 'WG'
                uI = Lagrangeinterpolate(pde.exactu,node,elem,'WG',eqn.face);
        end
        erruIuh(k) = sqrt((uh-uI)'*eqn.A*(uh-uI));
        errMax(k) = max(abs(uh-uI));
    end
    errTime(k) = cputime - t;
    % record time
    solverTime(k) = info.solverTime;
    assembleTime(k) = info.assembleTime;
    if option.printlevel>1
        fprintf('Time to compute the error %4.2g s \n H1 err %4.2g    L2err %4.2g \n',...
                 errTime(k), errH1(k), errL2(k));
    end
    % record solver information
    itStep(k) = info.itStep;
    stopErr(k) = info.stopErr;
    flag(k) = info.flag;
    % plot
    N(k) = length(soln.u);
    h(k) = 1./(size(node,1)^(1/3)-1);
    if  strcmp(elemType,'WG') % modify size for WG
        if ~isfield(option,'reducesystem') || (option.reducesystem == 1)
            N(k) = N(k) - size(elem,1); % reduced system
        end
    end
    if option.plotflag && N(k) < 2e4 % show mesh and solution for small size
        switch elemType
        case 'P1'     % piecewise linear function P1 element
            figure(1);  showresult3(node,elem,uh);
        case 'CR'     % piecewise linear function CR element
            continue;
        case 'P2'     % piecewise quadratic function
            figure(1);  showresult3(node,elem,uh(1:size(node,1)));
        case 'WG'     % weak Galerkin element
            continue;
        end
    end
    if N(k) > maxN
        break;
    end
    % refine mesh
    t = cputime;
    if strcmp(refType,'red')
        [node,elem,bdFlag] = uniformrefine3(node,elem,bdFlag);
    elseif strcmp(refType,'bisect')
        [node,elem,bdFlag] = uniformbisect3(node,elem,bdFlag);
    end
    meshTime(k) = cputime - t;
end

%% Plot convergence rates
if option.rateflag
    figure;
    set(gcf,'Units','normal');
    set(gcf,'Position',[0.25,0.25,0.55,0.4]);
    subplot(1,2,1)
    showrateh2(h(1:k),errH1(1:k),1,'-*','||Du-Du_h||',...
               h(1:k),errL2(1:k),1,'k-+','||u-u_h||');
    title(['Convergence Rate of 3D ', elemType, '-element']);
    subplot(1,2,2)
    showrateh2(h(1:k),erruIuh(1:k),1,'m-+','||Du_I-Du_h||',...
               h(1:k),errMax(1:k),1,'r-*','||u_I-u_h||_{\infty}');
    title(['Convergence Rate of 3D ', elemType, '-element']);
end

%% Output
err = struct('h',h(1:k),'N',N,'H1',errH1(1:k),'L2',errL2(1:k),...
             'uIuhH1',erruIuh(1:k),'uIuhMax',errMax(1:k));
time = struct('N',N,'err',errTime(1:k),'solver',solverTime(1:k), ...
              'assemble',assembleTime(1:k),'mesh',meshTime(1:k));
solver = struct('N',N(1:k),'itStep',itStep(1:k),'time',solverTime(1:k),...
                'stopErr',stopErr(1:k),'flag',flag(1:k));

%% Display error and time
if option.dispflag
    disp('Table: Error')
    colname = {'#Dof','h','||u-u_h||','||Du-Du_h||','||DuI-Du_h||','||uI-u_h||_{max}'};
    disptable(colname,err.N,[],err.h,'%0.3e',err.L2,'%0.5e',err.H1,'%0.5e',...
              err.uIuhH1,'%0.5e',err.uIuhMax,'%0.5e');

    disp('Table: CPU time')
    colname = {'#Dof','Assemble','Solve','Error','Mesh'};
    disptable(colname,time.N,[],time.assemble,'%0.2e',time.solver,'%0.2e',...
                      time.err,'%0.2e',time.mesh,'%0.2e');
end

repo: ifem


text:
function status = testYamlConversion
% test the functions for yaml import/export see if the conversion process
% changes the model content
%
% Usage: status = testYamlConversion
%

% Get model path
[ST, I]=dbstack('-completenames');
modelPath=fileparts(fileparts(fileparts(ST(I).file)));

% Import yaml model
ymlFile=fullfile(modelPath,'model','Human-GEM.yml');
model = importYaml(ymlFile, true);

% make sure there is no intermediate Yaml file under the current folder
warning('off', 'MATLAB:DELETE:FileNotFound')
if exist('testYamlConversion.yml','file')
    delete testYamlConversion.yml;
end

% export to yml and then import back
exportYaml(model,'testYamlConversion.yml');
importedHumanGEM = importYaml('testYamlConversion.yml', true);

% remove intermediate Yaml file
delete testYamlConversion.yml;

% compare the imported model from yaml with the original one
if isequal(model, importedHumanGEM)
    % model conversion between Matlab and Yaml files is successful
    disp('The conversion was successful.')
    status = 1;
else
    error('There are problems during the conversion between Matlab and Yaml files');
end

author: SysBioChalmers

repo: Human-GEM


text:
%PARSC Parse classifier
%
% 	PARSC(W)
%
% Displays the type and, for combining classifiers, the structure of the
% mapping W.
%
% See also MAPPINGS

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Physics, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: parsc.m,v 1.2 2006/03/08 22:06:58 duin Exp $

function parsc(w,space)

		% If W is not a mapping, do not generate an error but simply return.

	if (~isa(w,'prmapping')) return; end
	if (nargin == 1)
		space = '';
	end

	% Display W's type.

	display(w,space);

	% If W is a combining classifier, recursively call PARSC to plot the
	% combined classier.

	pars = getdata(w);
	if (iscell(pars))
		space = [space '  '];
		for i=1:length(pars)
			parsc(pars{i},space)
		end
	end

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [Ax,Ay,b,K] = convexhullConvex(varargin)
% Two lower bounds from tangents
% y > f(xL) + (x-xL)*df(xL)
% y > f(xU) + (x-xL)*df(xU)
% Upper bound from connecting extreme points
% y < f(xU)(x-xL)/(xU-xL) +  f(xL)(xU-x)/(xU-xL)
% can be written as
% Ax*x + Ay*y <= b

if rem(nargin,3)
    error('The convex hull generator assumes n triplets (x,f,df).')
end
m = nargin/3;
x = [varargin{(1:m)}]';
f = [varargin{(1:m)+m}]';
df = [varargin{(1:m)+2*m}]';

if df(1) == df(end)
    % This is a linear function
    % return y = f(1) + df*(x-x(1))
    Ax = df(1);
    Ay = -1;
    b = df(1)*x(1)-f(1);
    K.f = 1;
    K.l = 0;
    return
end

if all(diff(x)>=0) | all(diff(x(~isinf(x))))>=0 % Support [0 inf inf]
    Ay = [-ones(m,1);1];
    b  = [-f + x.*df; -f(end)*x(1)/(x(end)-x(1)) +  f(1)*x(end)/(x(end)-x(1))];
    Ax  = [df;-f(end)/(x(end)-x(1)) + f(1)/(x(end)-x(1))];

     % Don't use ill-conditioned cuts
    if df(1)<-1000
        Ax(1)=[];
        Ay(1) = [];
        b(1) = [];
    end
    if df(end-1)>1000
        Ax(end-1)=[];
        Ay(end-1) = [];
        b(end-1) = [];
    end
else
    Ax = [];
    Ay = [];
    b = [];
end
j = find(any(isnan([Ax Ay b]),2));
if ~isempty(j)
    b(j)=[];
    Ax(j,:)=[];
    Ay(j,:)=[];
end
K.f = 0;
K.l = length(b);
if 0
  sdpvar x y;plot(Ax*x+Ay*y <= b);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function varargout=lpcar2cc(varargin)
%
% For calling details please see v_lpcar2cc.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_lpcar2cc.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: lpcar2cc.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_lpcar2cc(varargin{:});
else
    v_lpcar2cc(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
%% DEMO_febio_0005_cube_tension_compression_shear
% Below is a demonstration for:
%
% * Building geometry for a cube with hexahedral elements
% * Defining the boundary conditions
% * Coding the febio structure
% * Running the model
% * Importing and visualizing the displacement and stress results

%% Keywords
%
% * febio_spec version 3.0
% * febio, FEBio
% * compression, tension, compressive, tensile, shear
% * displacement control, displacement boundary condition
% * hexahedral elements, hex8
% * cube, box, rectangular
% * static, solid
% * hyperelastic, Ogden
% * displacement logfile
% * stress logfile

%%

clear; close all; clc;

%% Plot settings
fontSize=20;
faceAlpha1=0.8;
markerSize=40;
markerSize2=20;
lineWidth=3;

%% Control parameters

% Path names
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'data','temp');

% Defining file names
febioFebFileNamePart='tempModel';
febioFebFileName=fullfile(savePath,[febioFebFileNamePart,'.feb']); %FEB file name
febioLogFileName=[febioFebFileNamePart,'.txt']; %FEBio log file name
febioLogFileName_disp=[febioFebFileNamePart,'_disp_out.txt']; %Log file name for exporting displacement
febioLogFileName_stress=[febioFebFileNamePart,'_stress_out.txt']; %Log file name for exporting stress

%Specifying dimensions and number of elements
cubeSize=10;
sampleWidth=cubeSize; %Width
sampleThickness=cubeSize; %Thickness
sampleHeight=cubeSize; %Height
pointSpacings=2*ones(1,3); %Desired point spacing between nodes
numElementsWidth=round(sampleWidth/pointSpacings(1)); %Number of elemens in dir 1
numElementsThickness=round(sampleThickness/pointSpacings(2)); %Number of elemens in dir 2
numElementsHeight=round(sampleHeight/pointSpacings(3)); %Number of elemens in dir 3

%Define applied displacement
stretchLoad=1.3;
displacementMagnitude=(stretchLoad*sampleHeight)-sampleHeight; %The displacement magnitude

%Material parameter set
c1=1e-3; %Shear-modulus-like parameter
m1=8; %Material parameter setting degree of non-linearity
k_factor=1e2; %Bulk modulus factor
k=c1*k_factor; %Bulk modulus

% FEA control settings
numTimeSteps=10; %Number of time steps desired
max_refs=25; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=6; %Optimum number of iterations
max_retries=5; %Maximum number of retires
dtmin=(1/numTimeSteps)/100; %Minimum time step size
dtmax=1/numTimeSteps; %Maximum time step size

%% Creating model geometry and mesh
% A box is created with tri-linear hexahedral (hex8) elements using the
% |hexMeshBox| function. The function offers the boundary faces with
% seperate labels for the top, bottom, left, right, front, and back sides.
% As such these can be used to define boundary conditions on the exterior.

% Create a box with hexahedral elements
cubeDimensions=[sampleWidth sampleThickness sampleHeight]; %Dimensions
cubeElementNumbers=[numElementsWidth numElementsThickness numElementsHeight]; %Number of elements
outputStructType=2; %A structure compatible with mesh view
[meshStruct]=hexMeshBox(cubeDimensions,cubeElementNumbers,outputStructType);

%Access elements, nodes, and faces from the structure
E=meshStruct.elements; %The elements
V=meshStruct.nodes; %The nodes (vertices)
Fb=meshStruct.facesBoundary; %The boundary faces
Cb=meshStruct.boundaryMarker; %The "colors" or labels for the boundary faces
elementMaterialIndices=ones(size(E,1),1); %Element material indices

%%
% Plotting model boundary surfaces and a cut view

hFig=cFigure;

subplot(1,2,1); hold on;
title('Model boundary surfaces and labels','FontSize',fontSize);
gpatch(Fb,V,Cb,'k',faceAlpha1);
colormap(gjet(6)); icolorbar;
axisGeom(gca,fontSize);

hs=subplot(1,2,2); hold on;
title('Cut view of solid mesh','FontSize',fontSize);
optionStruct.hFig=[hFig hs];
meshView(meshStruct,optionStruct);
axisGeom(gca,fontSize);

drawnow;

%% Defining the boundary conditions
% The visualization of the model boundary shows colors for each side of the
% cube. These labels can be used to define boundary conditions.

%Define supported node sets
logicFace=Cb==5; %Logic for current face set
Fr=Fb(logicFace,:); %The current face set
bcSupportList=unique(Fr(:)); %Node set part of selected face

%Prescribed displacement nodes
logicPrescribe=Cb==6; %Logic for current face set
Fr=Fb(logicPrescribe,:); %The current face set
bcPrescribeList=unique(Fr(:)); %Node set part of selected face

%%
% Visualizing boundary conditions. Markers plotted on the semi-transparent
% model denote the nodes in the various boundary condition lists.

hf=cFigure;
title('Boundary conditions','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
hold on;

gpatch(Fb,V,'kw','k',0.5);

hl(1)=plotV(V(bcSupportList,:),'k.','MarkerSize',markerSize);
hl(2)=plotV(V(bcPrescribeList,:),'r.','MarkerSize',markerSize);

legend(hl,{'BC support','BC prescribe'});

axisGeom(gca,fontSize);
camlight headlight;
drawnow;

%% Defining the FEBio input structure
% See also |febioStructTemplate| and |febioStruct2xml| and the FEBio user
% manual.

%Get a template with default settings
[febio_spec]=febioStructTemplate;

%febio_spec version
febio_spec.ATTR.version='3.0';

%Module section
febio_spec.Module.ATTR.type='solid';

%Create control structure for use by all steps
stepStruct.Control.time_steps=numTimeSteps;
stepStruct.Control.step_size=1/numTimeSteps;
stepStruct.Control.solver.max_refs=max_refs;
stepStruct.Control.solver.max_ups=max_ups;
stepStruct.Control.time_stepper.dtmin=dtmin;
stepStruct.Control.time_stepper.dtmax=dtmax;
stepStruct.Control.time_stepper.max_retries=max_retries;
stepStruct.Control.time_stepper.opt_iter=opt_iter;

%Add template based default settings to proposed control section
[stepStruct.Control]=structComplete(stepStruct.Control,febio_spec.Control,1); %Complement provided with default if missing

%Remove control field (part of template) since step specific control sections are used
febio_spec=rmfield(febio_spec,'Control');

febio_spec.Step.step{1}.Control=stepStruct.Control;
febio_spec.Step.step{1}.ATTR.id=1;
febio_spec.Step.step{2}.Control=stepStruct.Control;
febio_spec.Step.step{2}.ATTR.id=2;
febio_spec.Step.step{3}.Control=stepStruct.Control;
febio_spec.Step.step{3}.ATTR.id=3;
febio_spec.Step.step{4}.Control=stepStruct.Control;
febio_spec.Step.step{4}.ATTR.id=4;
febio_spec.Step.step{5}.Control=stepStruct.Control;
febio_spec.Step.step{5}.ATTR.id=5;

%Material section
materialName1='Material1';
febio_spec.Material.material{1}.ATTR.name=materialName1;
febio_spec.Material.material{1}.ATTR.type='Ogden';
febio_spec.Material.material{1}.ATTR.id=1;
febio_spec.Material.material{1}.c1=c1;
febio_spec.Material.material{1}.m1=m1;
febio_spec.Material.material{1}.c2=c1;
febio_spec.Material.material{1}.m2=-m1;
febio_spec.Material.material{1}.k=k;

% Mesh section
% -> Nodes
febio_spec.Mesh.Nodes{1}.ATTR.name='Object1'; %The node set name
febio_spec.Mesh.Nodes{1}.node.ATTR.id=(1:size(V,1))'; %The node id's
febio_spec.Mesh.Nodes{1}.node.VAL=V; %The nodel coordinates

% -> Elements
partName1='Part1';
febio_spec.Mesh.Elements{1}.ATTR.name=partName1; %Name of this part
febio_spec.Mesh.Elements{1}.ATTR.type='hex8'; %Element type
febio_spec.Mesh.Elements{1}.elem.ATTR.id=(1:1:size(E,1))'; %Element id's
febio_spec.Mesh.Elements{1}.elem.VAL=E; %The element matrix

% -> NodeSets
nodeSetName1='bcSupportList';
nodeSetName2='bcPrescribeList';

febio_spec.Mesh.NodeSet{1}.ATTR.name=nodeSetName1;
febio_spec.Mesh.NodeSet{1}.node.ATTR.id=bcSupportList(:);

febio_spec.Mesh.NodeSet{2}.ATTR.name=nodeSetName2;
febio_spec.Mesh.NodeSet{2}.node.ATTR.id=bcPrescribeList(:);

%MeshDomains section
febio_spec.MeshDomains.SolidDomain.ATTR.name=partName1;
febio_spec.MeshDomains.SolidDomain.ATTR.mat=materialName1;

%Boundary condition section
%-> Fix boundary conditions
febio_spec.Boundary.bc{1}.ATTR.type='fix';
febio_spec.Boundary.bc{1}.ATTR.node_set=nodeSetName1;
febio_spec.Boundary.bc{1}.dofs='x,y,z';

%STEP 1 Tension
febio_spec.Step.step{1}.Boundary.bc{1}.ATTR.type='prescribe';
febio_spec.Step.step{1}.Boundary.bc{1}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{1}.Boundary.bc{1}.dof='z';
febio_spec.Step.step{1}.Boundary.bc{1}.scale.ATTR.lc=1;
febio_spec.Step.step{1}.Boundary.bc{1}.scale.VAL=displacementMagnitude;
febio_spec.Step.step{1}.Boundary.bc{1}.relative=1;

febio_spec.Step.step{1}.Boundary.bc{2}.ATTR.type='fix';
febio_spec.Step.step{1}.Boundary.bc{2}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{1}.Boundary.bc{2}.dofs='x,y';

%STEP 2 Return form tension
febio_spec.Step.step{2}.Boundary.bc{1}.ATTR.type='prescribe';
febio_spec.Step.step{2}.Boundary.bc{1}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{2}.Boundary.bc{1}.dof='z';
febio_spec.Step.step{2}.Boundary.bc{1}.scale.ATTR.lc=2;
febio_spec.Step.step{2}.Boundary.bc{1}.scale.VAL=-displacementMagnitude;
febio_spec.Step.step{2}.Boundary.bc{1}.relative=1;

febio_spec.Step.step{2}.Boundary.bc{2}.ATTR.type='fix';
febio_spec.Step.step{2}.Boundary.bc{2}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{2}.Boundary.bc{2}.dofs='x,y';

%STEP 3 Compression
febio_spec.Step.step{3}.Boundary.bc{1}.ATTR.type='prescribe';
febio_spec.Step.step{3}.Boundary.bc{1}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{3}.Boundary.bc{1}.dof='z';
febio_spec.Step.step{3}.Boundary.bc{1}.scale.ATTR.lc=3;
febio_spec.Step.step{3}.Boundary.bc{1}.scale.VAL=-displacementMagnitude;
febio_spec.Step.step{3}.Boundary.bc{1}.relative=1;

febio_spec.Step.step{3}.Boundary.bc{2}.ATTR.type='fix';
febio_spec.Step.step{3}.Boundary.bc{2}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{3}.Boundary.bc{2}.dofs='x,y';

%STEP 4 Return from compression
febio_spec.Step.step{4}.Boundary.bc{1}.ATTR.type='prescribe';
febio_spec.Step.step{4}.Boundary.bc{1}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{4}.Boundary.bc{1}.dof='z';
febio_spec.Step.step{4}.Boundary.bc{1}.scale.ATTR.lc=4;
febio_spec.Step.step{4}.Boundary.bc{1}.scale.VAL=displacementMagnitude;
febio_spec.Step.step{4}.Boundary.bc{1}.relative=1;

febio_spec.Step.step{4}.Boundary.bc{2}.ATTR.type='fix';
febio_spec.Step.step{4}.Boundary.bc{2}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{4}.Boundary.bc{2}.dofs='x,y';

%STEP 5 Shear
febio_spec.Step.step{5}.Boundary.bc{1}.ATTR.type='prescribe';
febio_spec.Step.step{5}.Boundary.bc{1}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{5}.Boundary.bc{1}.dof='x';
febio_spec.Step.step{5}.Boundary.bc{1}.scale.ATTR.lc=5;
febio_spec.Step.step{5}.Boundary.bc{1}.scale.VAL=displacementMagnitude;
febio_spec.Step.step{5}.Boundary.bc{1}.relative=1;

febio_spec.Step.step{5}.Boundary.bc{2}.ATTR.type='fix';
febio_spec.Step.step{5}.Boundary.bc{2}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{5}.Boundary.bc{2}.dofs='y,z';

%LoadData section
% -> load_controller
febio_spec.LoadData.load_controller{1}.ATTR.id=1;
febio_spec.LoadData.load_controller{1}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{1}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{1}.points.point.VAL=[0 0; 1 1];

febio_spec.LoadData.load_controller{2}.ATTR.id=2;
febio_spec.LoadData.load_controller{2}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{2}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{2}.points.point.VAL=[1 0; 2 1];

febio_spec.LoadData.load_controller{3}.ATTR.id=3;
febio_spec.LoadData.load_controller{3}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{3}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{3}.points.point.VAL=[2 0; 3 1];

febio_spec.LoadData.load_controller{4}.ATTR.id=4;
febio_spec.LoadData.load_controller{4}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{4}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{4}.points.point.VAL=[3 0; 4 1];

febio_spec.LoadData.load_controller{5}.ATTR.id=5;
febio_spec.LoadData.load_controller{5}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{5}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{5}.points.point.VAL=[4 0; 5 1];

%Output section
% -> log file
febio_spec.Output.logfile.ATTR.file=febioLogFileName;
febio_spec.Output.logfile.node_data{1}.ATTR.file=febioLogFileName_disp;
febio_spec.Output.logfile.node_data{1}.ATTR.data='ux;uy;uz';
febio_spec.Output.logfile.node_data{1}.ATTR.delim=',';
febio_spec.Output.logfile.node_data{1}.VAL=1:size(V,1);

febio_spec.Output.logfile.element_data{1}.ATTR.file=febioLogFileName_stress;
febio_spec.Output.logfile.element_data{1}.ATTR.data='s1';
febio_spec.Output.logfile.element_data{1}.ATTR.delim=',';
febio_spec.Output.logfile.element_data{1}.VAL=1:size(E,1);

%% Quick viewing of the FEBio input file structure
% The |febView| function can be used to view the xml structure in a MATLAB
% figure window.

%%
% |febView(febio_spec); %Viewing the febio file|

%% Exporting the FEBio input file
% Exporting the febio_spec structure to an FEBio input file is done using
% the |febioStruct2xml| function.

febioStruct2xml(febio_spec,febioFebFileName); %Exporting to file and domNode

%% Running the FEBio analysis
% To run the analysis defined by the created FEBio input file the
% |runMonitorFEBio| function is used. The input for this function is a
% structure defining job settings e.g. the FEBio input file name. The
% optional output runFlag informs the user if the analysis was run
% succesfully.

febioAnalysis.run_filename=febioFebFileName; %The input file name
febioAnalysis.run_logname=febioLogFileName; %The name for the log file
febioAnalysis.disp_on=1; %Display information on the command window
febioAnalysis.runMode='external';%'internal';

[runFlag]=runMonitorFEBio(febioAnalysis);%START FEBio NOW!!!!!!!!

%% Import FEBio results

if runFlag==1 %i.e. a succesful run

    %%
    % Importing nodal displacements from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_disp),1,1);

    %Access data
    N_disp_mat=dataStruct.data; %Displacement
    timeVec=dataStruct.time; %Time

    %Create deformed coordinate set
    V_DEF=N_disp_mat+repmat(V,[1 1 size(N_disp_mat,3)]);

    %%
    % Importing element stress from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_stress),1,1);

    %Access data
    E_stress_mat=dataStruct.data;

    %%
    % Plotting the simulated results using |anim8| to visualize and animate
    % deformations

    DN_magnitude=sqrt(sum(N_disp_mat(:,:,end).^2,2)); %Current displacement magnitude

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);
    title('Displacement magnitude [mm]','Interpreter','Latex')
    hp=gpatch(Fb,V_DEF(:,:,end),DN_magnitude,'k',1); %Add graphics object to animate
    hp.Marker='.';
    hp.MarkerSize=markerSize2;
    hp.FaceColor='interp';
    gpatch(Fb,V,0.5*ones(1,3),'k',0.25); %A static graphics object

    axisGeom(gca,fontSize);
    colormap(gjet(250)); colorbar;
    caxis([0 max(DN_magnitude)]);
    axis(axisLim(V_DEF)); %Set axis limits statically
        camlight headlight;

    % Set up animation features
    animStruct.Time=timeVec; %The time vector
    for qt=1:1:size(N_disp_mat,3) %Loop over time increments
        DN_magnitude=sqrt(sum(N_disp_mat(:,:,qt).^2,2)); %Current displacement magnitude

        %Set entries in animation structure
        animStruct.Handles{qt}=[hp hp]; %Handles of objects to animate
        animStruct.Props{qt}={'Vertices','CData'}; %Properties of objects to animate
        animStruct.Set{qt}={V_DEF(:,:,qt),DN_magnitude}; %Property values for to set in order to animate
    end
    anim8(hf,animStruct); %Initiate animation feature
    drawnow;

        %%
    % Plotting the simulated results using |anim8| to visualize and animate
    % deformations

    [CV]=faceToVertexMeasure(E,V,E_stress_mat(:,:,end));

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);
    title('$\sigma_{1}$ [MPa]','Interpreter','Latex')
    hp=gpatch(Fb,V_DEF(:,:,end),CV,'k',1); %Add graphics object to animate
    hp.Marker='.';
    hp.MarkerSize=markerSize2;
    hp.FaceColor='interp';
    gpatch(Fb,V,0.5*ones(1,3),'k',0.25); %A static graphics object

    axisGeom(gca,fontSize);
    colormap(gjet(250)); colorbar;
%     caxis([min(E_stress_mat(:)) max(E_stress_mat(:))]);
    axis(axisLim(V_DEF)); %Set axis limits statically
    camlight headlight;

    % Set up animation features
    animStruct.Time=timeVec; %The time vector
    for qt=1:1:size(N_disp_mat,3) %Loop over time increments

        [CV]=faceToVertexMeasure(E,V,E_stress_mat(:,:,qt));

        %Set entries in animation structure
        animStruct.Handles{qt}=[hp hp]; %Handles of objects to animate
        animStruct.Props{qt}={'Vertices','CData'}; %Properties of objects to animate
        animStruct.Set{qt}={V_DEF(:,:,qt),CV}; %Property values for to set in order to animate
    end
    anim8(hf,animStruct); %Initiate animation feature
    drawnow;

end

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function dsq=define_distances_btw_control_points()

% Copyright (C) <2007>  <Francesc Moreno-Noguer, Vincent Lepetit, Pascal Fua>
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the version 3 of the GNU General Public License
% as published by the Free Software Foundation.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
% Francesc Moreno-Noguer, CVLab-EPFL, September 2007.
% fmorenoguer@gmail.com, http://cvlab.epfl.ch/~fmoreno/

%relative coordinates of the control points
c1=[1,0,0];
c2=[0,1,0];
c3=[0,0,1];
c4=[0,0,0];

d12=(c1(1)-c2(1))^2 + (c1(2)-c2(2))^2 + (c1(3)-c2(3))^2;
d13=(c1(1)-c3(1))^2 + (c1(2)-c3(2))^2 + (c1(3)-c3(3))^2;
d14=(c1(1)-c4(1))^2 + (c1(2)-c4(2))^2 + (c1(3)-c4(3))^2;
d23=(c2(1)-c3(1))^2 + (c2(2)-c3(2))^2 + (c2(3)-c3(3))^2;
d24=(c2(1)-c4(1))^2 + (c2(2)-c4(2))^2 + (c2(3)-c4(3))^2;
d34=(c3(1)-c4(1))^2 + (c3(2)-c4(2))^2 + (c3(3)-c4(3))^2;

dsq=[d12,d13,d14,d23,d24,d34]';

author: cvlab-epfl

repo: EPnP

save_path: github-repos/MATLAB/cvlab-epfl-EPnP


text:
function [E,J]=SynthMeasWatsonSHCylNeuman_PGSE(x, grad_dirs, G, delta, smalldel, fibredir, roots)
% Substrate: Impermeable cylinders with one radius in an empty background.
% Orientation distribution: Watson's distribution with SH approximation
% Pulse sequence: Pulsed gradient spin echo
% Signal approximation: Gaussian phase distribution.
%
% [E,J]=SynthMeasWatsonSHCylNeuman_PGSE(x, grad_dirs, G, delta, smalldel, fibredir, roots)
% returns the measurements E according to the model and the Jacobian J of the
% measurements with respect to the parameters.  The Jacobian does not
% include derivates with respect to the fibre direction.
%
% x is the list of model parameters in SI units:
% x(1) is the diffusivity of the material inside the cylinders.
% x(2) is the radius of the cylinders.
% x(3) is the concentration parameter of the Watson's distribution
%
% grad_dirs is the gradient direction for each measurement.  It has size [N
% 3] where N is the number of measurements.
%
% G, delta and smalldel are the gradient strength, pulse separation and
% pulse length of each measurement in the protocol.  Each has
% size [N 1].
%
% fibredir is a unit vector along the symmetry axis of the Watson's
% distribution.  It must be in Cartesian coordinates [x y z]' with size [3 1].
%
% roots contains solutions to the Bessel function equation from function
% BesselJ_RootsCyl.
%
% author: Gary Hui Zhang (gary.zhang@ucl.ac.uk)
%

if length(x) ~= 3
    error('the first argument should have exactly three parameters');
end

d=x(1);
R=x(2);
kappa=x(3);

l_q = size(grad_dirs,1);

% Parallel component
if nargout > 1
    [LePar, J_LePar] = CylNeumanLePar_PGSE(d, G, delta, smalldel);
else
    LePar = CylNeumanLePar_PGSE(d, G, delta, smalldel);
end

% Perpendicular component
if nargout > 1
    [LePerp, J_LePerp] = CylNeumanLePerp_PGSE(d, R, G, delta, smalldel, roots);
else
    LePerp = CylNeumanLePerp_PGSE(d, R, G, delta, smalldel, roots);
end
ePerp = exp(LePerp);

% Compute the Legendre weighted signal
Lpmp = LePerp - LePar;
if nargout > 1
    [lgi, J_lgi] = LegendreGaussianIntegral(Lpmp, 6);
else
    lgi = LegendreGaussianIntegral(Lpmp, 6);
end

% Compute the spherical harmonic coefficients of the Watson's distribution
if nargout > 1
    [coeff, J_coeff] = WatsonSHCoeff(kappa);
else
    coeff = WatsonSHCoeff(kappa);
end
coeffMatrix = repmat(coeff, [l_q, 1]);

% Compute the dot product between the symmetry axis of the Watson's distribution
% and the gradient direction
%
% For numerical reasons, cosTheta might not always be between -1 and 1
% Due to round off errors, individual gradient vectors in grad_dirs and the
% fibredir are never exactly normal.  When a gradient vector and fibredir are
% essentially parallel, their dot product can fall outside of -1 and 1.
%
% BUT we need make sure it does, otherwise the legendre function call below
% will FAIL and abort the calculation!!!
%
cosTheta = grad_dirs*fibredir;
badCosTheta = find(abs(cosTheta)>1);
cosTheta(badCosTheta) = cosTheta(badCosTheta)./abs(cosTheta(badCosTheta));

% Compute the SH values at cosTheta
sh = zeros(size(coeff));
shMatrix = repmat(sh, [l_q, 1]);
for i = 1:7
    shMatrix(:,i) = sqrt((i - .75)/pi);
    % legendre function returns coefficients of all m from 0 to l
    % we only need the coefficient corresponding to m = 0
    % WARNING: make sure to input ROW vector as variables!!!
    % cosTheta is expected to be a COLUMN vector.
    tmp = legendre(2*i - 2, cosTheta');
    tmp = tmp';
    shMatrix(:,i) = shMatrix(:,i) .* tmp(:,1);
end

E = sum(lgi.*coeffMatrix.*shMatrix, 2);
% with the SH approximation, there will be no guarantee that E will be positive
% but we need to make sure it does!!! replace the negative values with 10% of
% the smallest positive values
E(find(E<=0)) = min(E(find(E>0)))*0.1;
E = 0.5*E.*ePerp;

% Compute the Jacobian matrix
if(nargout>1)
    % dePerp/dd
    dePerpdd = E.*J_LePerp(1);
    % dePar/dd
    dElgi = sum(J_lgi.*coeffMatrix.*shMatrix, 2);
    dePardd = 0.5*dElgi.*(J_LePerp(:,:,1) - J_LePar).*ePerp;
    % dE/dd
    dEdd = dePardd + dePerpdd;

	% dePerp/dR
    dePerpdR = E.*J_LePerp(2);
    % dePar/dR
    dePardR = 0.5*dElgi.*J_LePerp(:,:,2).*ePerp;
    % dE/dR
    dEdR = dePardR + dePerpdR;

    % dE/dK
    J_coeffMatrix = repmat(J_coeff, [l_q, 1]);
    dEdk = sum(lgi.*J_coeffMatrix.*shMatrix,2);
    dEdk = 0.5*dEdk.*ePerp;

    % Construct the jacobian matrix.
    J = zeros(length(E), 3);
    J(:,1) = dEdd;
    J(:,2) = dEdR;
    J(:,3) = dEdk;
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [Vm]=patchCentre(F,V)

% function [Vm]=patchCentre(F,V)
% ------------------------------------------------------------------------
%
%
%
% Change log:
% 2014
% 2019/04/22 Fixed such that output matches dimentionality of input (output
% used to always be 3D even if input is 2D)
% 2019/04/22 Updated to handle cell input
% ------------------------------------------------------------------------

%%

if isa(F,'cell') %Cell of faces (e.g. potentially of different types)
    Vm=repmat({zeros(size(F,1),3)},size(F)); %Allocate Vm to be a cell like F
    for q=1:1:numel(F) %Loop over face sets
        Vm{q}=patchCentre(F{q},V); %Override cell entries by face centres
    end
else
    Vm=zeros(size(F,1),size(V,2)); %Allocate memory for
    for q=1:1:size(V,2)
        X=V(:,q);
        if size(F,1)==1 %Treat single element case since MATLAB is not consistent
            Vm(:,q)=mean(X(F));
        else
            Vm(:,q)=mean(X(F),2);
        end
    end
end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.


text:
function [SE_MR,SE_MR_asymptotic] = functionComputeSE_DL_MR_impairments(channelGainOverNoise,tau_c,M,K,L,p,rho,f,kappatUE,kapparBS,kappatBS,kapparUE)
%Compute DL SE for different transmit precoding schemes using Corollary
%6.6.
%
%INPUT:
%channelGainOverNoise = K x L x L matrix containing the average channel
%                       gains normalized by the noise variance of al
%                       channels (in dB). channelGainOverNoise(k,j,l) is
%                       the channel gain between between UE k in cell j
%                       and BS l
%tau_c                = Length of coherence block
%M                    = Number of antennas per BS
%K                    = Number of UEs per cell
%L                    = Number of BSs and cells
%p                    = Uplink transmit power per UE (same for everyone)
%rho                  = Downlink transmit power per UE (same for everyone)
%f                    = Pilot reuse factor
%kappatUE             = Hardware quality of the UEs' transmitters
%kapparBS             = Hardware quality of the BSs' receivers
%kappatBS             = Hardware quality of the BSs' transmitters
%kapparUE             = Hardware quality of the UEs' receivers
%
%OUTPUT:
%SE_MR            = K x L matrix where element (k,l) is the downlink SE of
%                   UE k in cell l achieved with MR precoding
%SE_MR_asymptotic = K x L matrix where element (k,l) is the asymptotic
%                   downlink SE of UE k in cell l achieved with MR precoding
%
%
%This Matlab function was developed to generate simulation results to:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Initiate the hardware qualities as perfect if these are not given as input
if nargin<9
    kappatUE = 1;
end

if nargin<10
    kapparBS = 1;
end

if nargin<11
    kappatBS = 1;
end

if nargin<12
    kapparUE = 1;
end

%Transform the (normalized) average channel gains from dB to linear scale
betaValues = 10.^(channelGainOverNoise/10);

%Compute length of pilot sequences
tau_p = f*K;

%Generate pilot pattern
if f == 1

    pilotPattern = ones(L,1);

elseif f == 2 %Only works in the running example with its 16 BSs

    pilotPattern = kron(ones(2,1),[1; 2; 1; 2; 2; 1; 2; 1]);

elseif f == 4 %Only works in the running example with its 16 BSs

    pilotPattern = kron(ones(2,1),[1; 2; 1; 2; 3; 4; 3; 4]);

elseif f == 16 %Only works in the running example with its 16 BSs

    pilotPattern = (1:L)';

end

%Prepare to store signal and interference terms for all UEs
signalTerm = zeros(K,L);
interferenceTerm = zeros(K,L);

%If the asymptotic SE should be computed
if nargout>1

    signalTerm_asymptotic = zeros(K,L);
    interferenceTerm_asymptotic = zeros(K,L);

end

%% Go through all cells
for j = 1:L

    %Compute the parameter in (6.50)
    G = (1 + kappatBS*(M-1))/(M*kappatBS*kapparUE);

    %Go through all UEs
    for k = 1:K

        %Extract cells that use same pilots as cell j
        groupMembers = find(pilotPattern==pilotPattern(j))';

        %Compute the parameter in (6.38)
        psi_jk = 1/( p*kappatUE*kapparBS*tau_p*sum(betaValues(k,groupMembers,j)) + p*(1-kappatUE*kapparBS)*sum(sum(betaValues(:,:,j))) + 1 );

        %Compute the numerator of (6.47)
        signalTerm(k,j) = rho*p*(betaValues(k,j,j))^2*tau_p*psi_jk*M;

        %Compute parameter in (6.49)
        Flijk = (1+p*betaValues(:,:,j)*psi_jk*( 1 - kappatUE*kapparBS + (1-kappatUE)*kappatBS*kapparBS*(M-1) ))/(kappatBS*kapparUE*kappatUE*kapparBS);

        %Compute the non-coherent interference that UE k in cell j causes
        %to all the UEs, which is a part of the first term in the
        %denominator of (6.47)
        interferenceTerm = interferenceTerm + rho*betaValues(:,:,j).*Flijk;

        %Compute the coherent interference that UE k in cell j causes to
        %all the pilot-sharing UEs, which is a part of the second term in
        %the denominator of (6.47)
        interferenceTerm(k,groupMembers) = interferenceTerm(k,groupMembers) + rho*p*(betaValues(k,groupMembers,j)).^2*tau_p*psi_jk*M*G;

        %Compute the last two terms in the denominator of (6.47)
        interferenceTerm(k,j) = interferenceTerm(k,j) - signalTerm(k,j) + 1/(kappatUE*kapparBS*kappatBS*kapparUE);

        %If the asymptotic SE in (6.51) should be computed
        if nargout>1

            %Compute the numerator of (6.51), multiplied with psi_jk
            signalTerm_asymptotic(k,j) = rho*(betaValues(k,j,j))^2*psi_jk;

            %Compute the part of the first term in the denominator of
            %(6.51) that UE k in cell j contributes to
            interferenceTerm_asymptotic = interferenceTerm_asymptotic + rho*(betaValues(:,:,j)).^2*psi_jk*(1-kappatUE)/(kappatUE*kapparUE*tau_p);

            %Compute the part of the second term in the denominator of
            %(6.51) that UE k in cell j contributes to
            interferenceTerm_asymptotic(k,groupMembers) = interferenceTerm_asymptotic(k,groupMembers) +  rho*(betaValues(k,groupMembers,j)).^2*psi_jk/kapparUE;

            %Subtract the numerator from the denominator, which is an
            %alternative way of obtaining the third term in the denominator
            interferenceTerm_asymptotic(k,j) = interferenceTerm_asymptotic(k,j) - rho*(betaValues(k,j,j))^2*psi_jk;

        end

    end

end

%Compute the prelog factor assuming only downlink transmission
prelogFactor = (tau_c-tau_p)/tau_c;

%Compute the final SE expression based on Corollary 6.6
SE_MR = prelogFactor * log2(1+signalTerm ./ real(interferenceTerm));

if nargout>1

    %Compute the final asymptotic SE expression based on Corollary 6.7
    SE_MR_asymptotic = prelogFactor * log2(1+signalTerm_asymptotic ./ real(interferenceTerm_asymptotic));

end

author: emilbjornson

repo: massivemimobook


text:
function [CSF, seg]=mtv_mrQ_Seg_kmeans_simple(T1,BM,M0,mmPerVox,CSFVal)
% function [CSF, seg]=mtv_mrQ_Seg_kmeans_simple(T1,BM,M0)
%
% Clips brain to ventricles area, then segments using k-means with three
% clusters.
%
% This function uses FSL to segment into three tissues. It takes the CSF
% tissue and restricts it by the T1 values. The CSF is also restricted to
% be in the center of the brain in a box of approximately 60mm x 80mm x
% 40mm. Assuming that the brain is in AC-PC space, this is where the
% ventricles should be.
% Assumption:
%       CSFVal: R1 for water at body temperature (in 1/sec).
%                  Default is 0.35.
%
% ~INPUTS~
%       T1: VFA T1.
%       BM: Brain Mask.
%       M0: VFA M0.
%       CSFVal: Threshold on the R1 value of CSF (default=0.35s-1)
%
% ~OUTPUTS~
%       CSF: Refined mask of ventricules (box in the center + smooth)
%       seg: brain clustering (GM=1;DeepGray=2;WM=3;CSF=4)
%
% Adapted from:
% (C) Mezer lab, the Hebrew University of Jerusalem, Israel
%   2015
%
%

%% I. Loading and definitions
if ~exist('CSFVal','var')
    CSFVal=1/0.35; % R1 (in 1/sec) of water at body temp  (minimum value)
end
R1=1./T1;
BM = BM & ~isinf(R1);
%% II. Perform k-means in a "while" loop

    fprintf('\n Performing segmentation for CSF file ...              \n');

seg=zeros(size(R1));

mask= R1>CSFVal & BM;
notdone=0;

while notdone==0

    [IDX,C] =kmeans(R1(mask),3);
    seg(mask)=IDX;

    notdone=1;

    % Check we don't get a strange cluster that is very small in size.
    %(if so, this might be just noise)
    if  length(find(IDX==1))/length(IDX)<0.05
        mask=mask & seg~=1;
        notdone=0;
    end

    if  length(find(IDX==2))/length(IDX)<0.05
        mask=mask & seg~=2;
        notdone=0;
    end

    if  length(find(IDX==3))/length(IDX)<0.05
        mask=mask & seg~=3;
        notdone=0;
    end
end

% check if the clusters' means are too similar
if abs(1-C(1)/C(2))<0.1  &&  abs(1-C(1)/C(3))<0.1
    [IDX,C] =kmeans(R1(mask),1);
elseif abs(1-C(1)/C(2))<0.1
    [IDX,C] =kmeans(R1(mask),2);
elseif abs(1-C(1)/C(3))<0.1
    [IDX,C] =kmeans(R1(mask),2);
elseif abs(1-C(2)/C(3))<0.1
    [IDX,C] =kmeans(R1(mask),2);
end

%% III. Create segmentation file of the clipped brain
seg=zeros(size(R1));
seg(mask)=IDX;

% The tissue with the highest value is white matter (WM), the tissue with
% the lowest value is gray matter (GM), and the tissue with the
% intermediate value is the deep nuclei and the tissue between the WM and
% the GM.
%
% In some segmentations the lowest value is is air, intermediate is
% GM, and highest is WM.
%
% Either way, the order is maintained and we get a segmentation of GM, WM
% and CSF.

[val,idx]=sort(C);
GMclass=1;DEEPclass=2;WMclass=3; CSFclass=4;

seg(seg==idx(1))=4; seg(seg==idx(2))=5;seg(seg==idx(3))=6;
seg(seg==4)=GMclass; seg(seg==5)=DEEPclass; seg(seg==6)=WMclass;

CSF= R1<CSFVal & BM;
seg(CSF)=CSFclass; % any region that is mostly water.

% Clip mask size
boxsize = [30 40 20];
sz=size(CSF); szH=round(sz./2);
XX=boxsize(1)./round(mmPerVox(1));
YY=boxsize(2)./round(mmPerVox(2));
ZZ=boxsize(3)./round(mmPerVox(3));

CSF(szH(1)+XX:end,:,:)=0;
CSF(1:szH(1)-XX,:,:)=0;

CSF(:,1:szH(2)-YY,:)=0;
CSF(:,szH(2)+YY:end,:,:)=0;
if numel(szH)>2
    CSF(:,:,1:max(1,szH(3)-ZZ))=0;
    CSF(:,:,min(end,szH(3)+ZZ):end)=0;
end
%% IV. Smoothe in space
[CSF1] = ordfilt3D(CSF,6);
CSF1=CSF &  CSF1;

CSF2= CSF1 & R1<0.25 & R1>0.2 & M0<prctile(M0(BM),99);

%% V. Some issues

if length(find(CSF2))<200
           fprintf(['\n Warning: We could find only ' num2str(length(find(CSF2))) ' csf voxels. This makes the CSF WF estimation noisy. Consider reducing the CSF T1 threshold.              \n']);
end

% Larger ventricles region.
% This mask is good for cases when CSF ventricles are hard to identify.
% It may reduce the accuracy.
CSF= CSF & R1<0.25 & R1>0.2 & M0<prctile(M0(BM),99);


text:
function varargout=spm_XYZreg_Ex2(varargin)
% Example of Registry enabled XYZ GUI control / function
% FORMAT...
%_______________________________________________________________________
%
% Help goes here...
%
% Object must be indentifiable via a unique HandleGraphics object.
% In this code, this handle is called hMe.
%
% This HandleGraphics objects 'UserData' *must* be a structure.
% The structure must have a field called 'hReg', which stores the handle
% of the registry when linked, and is empty when not. Some utility features
% of spm_XYZreg will set/delete the handle directly...
%
% There must be a 'SetCoords' function for this object, with call:
%   spm_res_ui('SetCoords',xyz,hMe,hC)
% ...this can handle interna, co-ordinate setting (as in this example), but
% must also call the registry.
%
% The registry update function is:
%   spm_XYZreg('SetCoords',xyz,hReg,hMe);
% ...which must be called at all points where the local co-ordinates can be
% changed. It is robust to invalid or empty hReg.
%
% It's *vital* to specify caller handles (hC), so that the registry doesn't
% end up in an infinite loop of updating!
%
% Hey, if your function has multiple places where you can change the XYZ,
% you could use an ``internal'' registry locally, with the external registry
% as one of it's entries! (I think?)
%
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Andrew Holmes
% $Id: spm_XYZreg_Ex2.m 5219 2013-01-29 17:07:07Z spm $

%=======================================================================
switch lower(varargin{1}), case 'create'
%=======================================================================
% hMe = spm_XYZreg_Ex2('Create',M,D,xyz)
if nargin<4, xyz=[0;0;0]; else, xyz=varargin{4}; end
if nargin<3, error('Insufficient arguments'), end
D = varargin{3};
M = varargin{2};

xyz = spm_XYZreg('RoundCoords',xyz,M,D);

F = figure;

%-Create control:
% - note the UserData structure: hReg is necessary for cross-registration,
%   M, D, & xyz are used internally...
% - note how the DeleteFcn closes the objects figure...
hMe = uicontrol(F,'Style','Pushbutton',...
    'String',sprintf('[%.3f, %.3f, %.3f]',xyz),...
    'Position',[100 100 300 30],...
    'FontSize',14,'FontWeight','Bold',...
    'Callback',...
    'spm_XYZreg_Ex2(''SetCoords'',input(''Enter [x,y,z]'''' vector: ''),gcbo);',...
    'UserData',struct(...
        'hReg', [],...
        'M',    M,...
        'D',    D,...
        'xyz',  xyz ),...
    'DeleteFcn','delete(gcbf)');

varargout = {hMe};

%=======================================================================
case 'setcoords'    % Set co-ordinates
%=======================================================================
% [xyz,d] = spm_XYZreg_Ex2('SetCoords',xyz,hMe,hC)
if nargin<4, hC=0; else hC=varargin{4}; end
if nargin<3, error('Insufficient arguments'), else hMe=varargin{3}; end
%-Or possibly some clever code to guess the handle
if nargin<2, error('Set co-ords to what!'), else xyz=varargin{2}; end

UD = get(hMe,'UserData');

%-Check validity of coords only when called without a caller handle
%-----------------------------------------------------------------------
if hC<=0
    [xyz,d] = spm_XYZreg('RoundCoords',xyz,UD.M,UD.D);
    if d>0 && nargout<2, warning(sprintf(...
        '%s: Co-ords rounded to nearest voxel center: Discrepancy %.2f',...
        mfilename,d)), end
else
    d = [];
end

%-Update local XYZ information
%-----------------------------------------------------------------------
UD.xyz = xyz;
set(hMe,'UserData',UD)
set(hMe,'String',sprintf('[%.3f, %.3f, %.3f]',xyz))

%-Tell the registry, if we've not been called by the registry...
%-----------------------------------------------------------------------
if (~isempty(UD.hReg) && UD.hReg~=hC)
    spm_XYZreg('SetCoords',xyz,UD.hReg,hMe);
end

%-Return arguments
%-----------------------------------------------------------------------
varargout = {xyz,d};

%=======================================================================
otherwise
%=======================================================================
error('Unknown action string')

%=======================================================================
end

author: spm


text:
function r83row_print_part_test ( )

%*****************************************************************************80
%
%% R83ROW_PRINT_PART_TEST tests R83ROW_PRINT_PART.
%
%  Discussion:
%
%    An R83ROW is a (3,N) array of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2015
%
%  Author:
%
%    John Burkardt
%
  n = 10;
%
%  Set up the array by rows.
%
  v = [ ...
    11.0,  12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0; ...
    21.0,  22.0, 23.0, 24.0, 25.0, 26.0, 27.0, 28.0, 29.0, 30.0; ...
    31.0,  32.0, 33.0, 34.0, 35.0, 36.0, 37.0, 38.0, 39.0, 40.0  ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R83ROW_PRINT_PART_TEST\n' );
  fprintf ( 1, '  R83ROW_PRINT_PART prints part of an R83ROW.\n' );

  max_print = 2;
  r83row_print_part ( n, v, max_print, '  Output with MAX_PRINT = 2' );

  max_print = 5;
  r83row_print_part ( n, v, max_print, '  Output with MAX_PRINT = 5' );

  max_print = 25;
  r83row_print_part ( n, v, max_print, '  Output with MAX_PRINT = 25' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function nodesInFaces=findFaceNodes(mesh,group)
% Here's a little problem - outsideNodes can contain groups linked by long, thin chains
% These mess up the perimeter-finder since each group contains a closed perim.
% Get rid of them here by finding only those nodes that are in faces

facesInGroup=findFacesInGroup(mesh,group);
nodesInFaces=unique(mesh.uniqueFaceIndexList(facesInGroup,:));

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function T = sptoeplitz(col, row)
%SPTOEPLITZ   Sparse Toeplitz matrix.
%   SPTOEPLITZ(C, R) produces a sparse nonsymmetric Toeplitz matrix having
%   C as its first column and R as its first row. Neither C nor R needs to
%   be sparse. No full-size dense matrices are formed.
%
%   SPTOEPLITZ(R) is a sparse symmetric/Hermitian Toeplitz matrix.
%
%   Examples:
%     sptoeplitz( real( (1i).^(0:8) ) )   % 9x9, 41 nonzeros
%     sptoeplitz( [-2 1 zeros(1,9998)] ); % classic 2nd difference
%
% See also TOEPLITZ, SPDIAGS.

% Based on SPTOEPLITZ.M on the Mathworks File Exchange,
% Copyright (c) 2006 by Tobin Driscoll (tobin.driscoll@gmail.com).

% Developer note: This is needed by multiplication operators.

% This part is borrowed from built-in Toeplitz.
if ( nargin < 2 ) % Symmetric case
  col(1) = conj(col(1));
  row = col;
  col = conj(col);
else
  if ( col(1) ~= row(1) )
    warning('MATLAB:sptoeplitz:DiagonalConflict',['First element of ' ...
      'input column does not match first element of input row. ' ...
      '\n         Column wins diagonal conflict.'])
  end
end

% Size of result.
m = length(col(:));
n = length(row(:));

% Only use toeplitz if you have too... fairly slow.
if ( (m < 2e3) && (n < 2e3) )
    % Note: nnz(col) -- Number of nonzero elements in col
    if ( nnz(col) == 1 && nnz(row) == 1 )
        Ic = find(col);
        Ir = find(row);
        if ( Ic == 1 )
            T = spdiags(col(Ic)*ones(m, 1), 0, m, n);
        else
            T = spdiags([col(Ic)*ones(m, 1) row(Ir)*ones(n , 1)], ...
                [-Ic + 1, Ir - 1], m, n);
        end
    else
        T = toeplitz(col, row);
        T = sparse(T);
    end
else
    % Locate the nonzero diagonals.
    [ic, jc, sc] = find(col(:));
    row(1) = 0;  % not used
    [ir, jr, sr] = find(row(:));

    % Use spdiags for construction.
    d = [ ir - 1; 1 - ic ];
    B = repmat( [ sr; sc ].', min(m, n), 1 );
    T = spdiags(B, d, m, n);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% dftfilt2() - discrete complex wavelet filters
%
% Usage:
%   >> wavelet = dftfilt2( freqs, cycles, srate, cyclefact)
%
% Inputs:
%   freqs    - frequency array
%   cycles   - cycles array. If one value is given, all wavelets have
%              the same number of cycles. If two values are given, the
%              two values are used for the number of cycles at the lowest
%              frequency and at the highest frequency, with linear
%              interpolation between these values for intermediate
%              frequencies
%   srate    - sampling rate (in Hz)
%
%   cycleinc - ['linear'|'log'] increase mode if [min max] cycles is
%              provided in 'cycle' parameter. {default: 'linear'}
%   type     - ['sinus'|'morlet'] wavelet type is a sinusoid with
%              cosine (real) and sine (imaginary) parts tapered by
%              a Hanning or Morlet function. 'morlet' is a typical Morlet
%              wavelet (with p=2*pi and sigma=0.7) best matching the
%              'sinus' Hanning taper) {default: 'morlet'}
% Output:
%   wavelet - cell array of wavelet filters
%
% Note: The length of the window is automatically computed from the
%       number of cycles and is always made odd.
%
% Authors: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 3/28/2003

% Copyright (C) 3/28/2003 Arnaud Delorme 8, SCCN/INC/UCSD, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function wavelet = dftfilt2( freqs, cycles, srate, cycleinc, type);

    if nargin < 3
        error('3 arguments required');
    end;
    if nargin < 5
        type = 'morlet';
    end;

    % compute number of cycles at each frequency
    % ------------------------------------------
    if length(cycles) == 1
        cycles = cycles*ones(size(freqs));
    elseif length(cycles) == 2
        if nargin == 4 & strcmpi(cycleinc, 'log') % cycleinc
            cycles = linspace(log(cycles(1)), log(cycles(2)), length(freqs));
            cycles = exp(cycles);
        else
            cycles = linspace(cycles(1), cycles(2), length(freqs));
        end;
    end;

    % compute wavelet
    for index = 1:length(freqs)

        % number of cycles depend on window size
        % number of cycles automatically reduced if smaller window
        % note: as the number of cycle changes, the frequency shifts a little
        %       this has to be fixed

        winlen = cycles(index)*srate/freqs(index);
        winlenint = floor(winlen);
        if mod(winlenint,2) == 1, winlenint = winlenint+1; end;
        winval = linspace(winlenint/2, -winlenint/2, winlenint+1);

        if strcmpi(type, 'sinus') % Hanning
            win = exp(2i*pi*freqs(index)*winval/srate);
            wavelet{index} = win .* hanning(length(winval))';

        else % Morlet
            t = freqs(index)*winval/srate;
            p = 2*pi;
            s = cycles(index)/5;
            wavelet{index} = exp(j*t*p)/sqrt(2*pi) .* ...
                (exp(-t.^2/(2*s^2))-sqrt(2)*exp(-t.^2/(s^2)-p^2*s^2/4));
        end;
    end;

    return;

    % testing
    % -------
    wav1 = dftfilt2(5, 5, 256); wav1 = wav1{1};
    abs1 = linspace(-floor(length(wav1)),floor(length(wav1)), length(wav1));
    figure; plot(abs1, real(wav1), 'b');

    wav2 = dftfilt2(5, 3, 256); wav2 = wav2{1};
    abs2 = linspace(-floor(length(wav2)),floor(length(wav2)), length(wav2));
    hold on; plot(abs2, real(wav2), 'r');

    wav3 = dftfilt2(5, 1.4895990, 256); wav3 = wav3{1};
    abs3 = linspace(-floor(length(wav3)),floor(length(wav3)), length(wav3));
    hold on; plot(abs3, real(wav3), 'g');

    wav4 = dftfilt2(5, 8.73, 256); wav4 = wav4{1};
    abs4 = linspace(-floor(length(wav4)),floor(length(wav4)), length(wav4));
    hold on; plot(abs4, real(wav4), 'm');

    % more testing
    % ------------
    freqs = exp(linspace(0,log(10),10));
    win = dftfilt2(freqs, [3 3], 256, 'linear', 'sinus');
    win = dftfilt2(freqs, [3 3], 256, 'linear', 'morlet');

    freqs = [12.0008   13.2675   14.5341   15.8007   17.0674   18.3340   19.6007   20.8673   22.1339   23.4006   24.6672   25.9339   27.2005 28.4671   29.7338   31.0004   32.2670   33.5337   34.8003   36.0670   37.3336   38.6002   39.8669   41.1335   42.4002   43.6668 44.9334   46.2001   47.4667 ...
             48.7334   50.0000];

    win = dftfilt2(freqs, [3 12], 256, 'linear'); size(win)

    winsize = 0;
    for index = 1:length(win)
        winsize = max(winsize,length(win{index}));
    end;
    allwav = zeros(winsize, length(win));
    for index = 1:length(win)
        wav1 = win{index};
        abs1 = linspace(-(length(wav1)-1)/2,(length(wav1)-1)/2, length(wav1));
        allwav(abs1+(winsize+1)/2,index) = wav1(:);
    end;
    figure; imagesc(imag(allwav));

% syemtric hanning function
function w = hanning(n)
if ~rem(n,2)
   w = .5*(1 - cos(2*pi*(1:n/2)'/(n+1)));
   w = [w; w(end:-1:1)];
else
   w = .5*(1 - cos(2*pi*(1:(n+1)/2)'/(n+1)));
   w = [w; w(end-1:-1:1)];
end

author: PatternRecognition

repo: OpenBMI


text:
function [phi,theta] = sphCoord(Coord3D)
% [phi,theta] = sphCoord(Coord3D)
%
% This function transforms spherical 3D XYZ spherical coordinates to 2D
% (theta,phi) spherical coordinates.
%
% Input
% Coord3D: Data Matrix with the spherical XYZ coordinates at the first
% three columns respectively.
%
% Output
% theta: Latitude
% phi: Longitud.
%
% $Revision: 1.2 $  $Date: 2015/01/06 17:14:55 $
% Original Author: Jorge Luis Bernal Rusiel
% CVS Revision Info:
%    $Author: mreuter $
%    $Date: 2015/01/06 17:14:55 $
%    $Revision: 1.2 $

% [theta,phi,r] = cart2sph(overlayData(:,1),overlayData(:,2),overlayData(:,3));
% //It was confusion

X=Coord3D(:,1); Y=Coord3D(:,2); Z=Coord3D(:,3);
phi=atan2(Y,X);
theta=atan2(sqrt(X.^2+Y.^2),Z);


text:
function [servicesUrl, options] = getHostOptions(host)
%% Set the options associated with the services for host.
    csrfUrl = [host '/services'];
    servicesUrl = [host '/services_submit'];
    [cookie, csrftoken] = getSessionInfo(csrfUrl);
    header = ["Content-Type" "application/json"; ...
              "Accept" "application/json"; ...
              "X-CSRFToken" csrftoken; "Cookie" cookie];

    options = weboptions('MediaType', 'application/json', ...
                         'Timeout', 120, 'HeaderFields', header);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%TIME_8
%
%   The purpose of this test is to determine whether it pays of to do
%   special optimizations for p=1 q=2 and p=2 q=3%

Lr=[480000*sf^2,262144*sf^2,900*sf^2];
ar=[     600*sf,        512,       2];
Mr=[     800*sf,       1024,  600*sf];

for ii=1:length(Lr)

  L=Lr(ii);

  M=Mr(ii);
  a=ar(ii);

  g=rand(L,1);
  f=rand(L,1);

  [L, a, M]

  gf=comp_wfac(g,a,M);

  c1=mex_dgt_fac_1(f,gf,a,M,1);

  % This is an attempt of cacheflushing
  f=rand(L,1); gf=rand(size(gf));

  c2=mex_dgt_fac_7(f,gf,a,M,1);

end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [isrcs, idets, measures, measure_type] = nst_unformat_channels(channel_labels, warn_bad_channels)
% NST_UNFORMAT_CHANNELS extract sources, dectectors and measures information
% from channel labels with *homogeneous* type (eg wavelength or Hb).
%
%   [ISRCS, IDETS, MEAS, CHAN_TYPE] = NST_UNFORMAT_CHANNELS(CHANNEL_LABELS)
%       CHANNEL_LABELS (cell array of str):
%           each str is formatted as 'SxDyWLz' or 'SxDyHbt', where:
%               x: source index
%               y: detector index
%               z: wavelength
%               t: Hb type (O, R, T).
%           Consistency is asserted so that:
%               - channels are unique
%               - channel type is homogeneous
%           Warning is issued if:
%               - the number of measures per pair is not homogeneous
%                 eg one pair has only one wavelength
%
%           Examples: S1D2WL685, S01D7WL830, S3D01HbR
%
%
%        ISRCS (array of int): extracted source indexes
%        IDETS (array of int): extracted detector indexes
%        MEAS (array of int | cell array of str): extracted measure values
%        CHAN_TYPE (array of int): channel type (see NST_CHANNEL_TYPES for enum).
%
%   See also NST_UNFORMAT_CHANNEL, NST_CHANNEL_TYPES, NST_FORMAT_CHANNEL

assert(iscellstr(channel_labels));

if nargin < 2
    warn_bad_channels = 0;
end

nb_channels = length(channel_labels);
isrcs = zeros(1, nb_channels);
idets = zeros(1, nb_channels);
measures = cell(1, nb_channels);
mtypes = zeros(1, nb_channels);
reformated_channels = cell(1, nb_channels);
for ichan=1:nb_channels
    [isrc, idet, meas, mtype] = nst_unformat_channel(channel_labels{ichan}, warn_bad_channels);
    isrcs(ichan) = isrc;
    idets(ichan) = idet;
    measures{ichan} = meas;
    mtypes(ichan) = mtype;
    if ~isnan(isrc) && ~isnan(idet)
        % reformat channels to make sure formatting is consistent
        % -> allow proper duplicate detection after
        reformated_channels{ichan} = nst_format_channel(isrc, idet, meas);
    else
        reformated_channels{ichan} = '';
    end
end

% check uniqueness
[~, i_unique] = unique(reformated_channels);
duplicates = reformated_channels;
duplicates(i_unique) = [];
duplicates(strcmp(duplicates, '')) = []; %remove unrecognized channels
i_duplicates = ismember(reformated_channels, unique(duplicates));
if ~isempty(duplicates)
    msg = sprintf('Duplicated channels: "%s". Indexes: %s', ...
                  strjoin(channel_labels(i_duplicates), ', '), num2str(i_duplicates));
    throw(MException('NIRSTORM:NonUniqueChannels', msg));
end

% check homogeneity of channel type
mtypes = mtypes(~isnan(mtypes));
if length(unique(mtypes)) > 1
    throw(MException('NIRSTORM:NonHomogeneousMeasure', 'Measure type is not homogeneous.'));
end

measure_types = nst_measure_types();
measure_type = mtypes(1);
if measure_type == measure_types.WAVELENGTH
   measures = cell2mat(measures);
end

%% Check number of measures per pair
max_idet = (max(idets)+1);
pairs_hash = isrcs(~isnan(isrcs)) * max_idet + idets(~isnan(isrcs));
unique_pairs_hash = unique(pairs_hash);
mcounts = sparse(ones(1, length(unique_pairs_hash)), unique_pairs_hash, ...
                      ones(1, length(unique_pairs_hash)), 1, max(unique_pairs_hash), ...
                      nb_channels);
if ~iscell(measures)
    all_measures = unique(measures(~isnan(measures)));
else
    all_measures = unique(measures(cellfun(@(e) ischar(e) || ~isnan(e), measures)));
end
for ichan=1:nb_channels
    if ~isnan(isrcs(ichan))
        if measure_type == measure_types.WAVELENGTH
            measure_hash = find(measures(ichan)==all_measures);
        elseif measure_type == measure_types.HB
            measure_hash = find(strcmp(measures(ichan),all_measures));
        end
        h_idx = isrcs(ichan) * max_idet + idets(ichan);
        mcounts(h_idx) = mcounts(h_idx) + measure_hash;
    end
end

% pairs with too many measures:
i_inconsistant_measures = (mcounts ~= sum(1:length(all_measures)) + 1) & (mcounts ~= 0);

if nnz(i_inconsistant_measures) > 0
    [isrc_incon,idet_incon,counts] = find(i_inconsistant_measures);
    inconsistent_pairs = cell(1, length(isrc_incon));
    for ipair=1:length(isrc_incon)
        inconsistent_pairs{ipair} = sprintf('S%dD%d', isrc_incon(ipair), idet_incon(ipair));
    end
    warning('Inconsistent measure(s) for pair(s): %s', strjoin(inconsistent_pairs, ', '));
end

end

author: brainstorm-tools


text:
classdef DASCMOP8 < PROBLEM
% <multi> <real> <large/none> <constrained>
% Difficulty-adjustable and scalable constrained benchmark MOP

%------------------------------- Reference --------------------------------
% Z. Fan, W. Li, X. Cai, H. Li, C. Wei, Q. Zhang, K. Deb, and E. Goodman,
% Difficulty adjustable and scalable constrained multi-objective test
% problem toolkit, Evolutionary Computation, 2020, 28(3): 339-378.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Wenji Li

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 3;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            x_all = X(:,3:1:end);
            g_1   = sum((x_all - 0.5) .* (x_all - 0.5) - cos(20.0 * pi .* ( x_all - 0.5)),2);
            sum1  = obj.D - 2 + g_1;
            PopObj(:,1) = cos(0.5 * pi * X(:,1)) .* cos(0.5 * pi * X(:,2)) + sum1;
            PopObj(:,2) = cos(0.5 * pi * X(:,1)) .* sin(0.5 * pi * X(:,2)) + sum1;
            PopObj(:,3) = sin(0.5 * pi * X(:,1)) + sum1;
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,X)
            x_all  = X(:,3:1:end);
            g_1    = sum((x_all - 0.5) .* (x_all - 0.5) - cos(20.0 * pi .* ( x_all - 0.5)),2);
            sum1   = obj.D - 2 + g_1;
            PopCon = Constraint(X(:,1),X(:,2),sum1);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R = UniformPoint(N,3);
            R = R./repmat(sqrt(sum(R.^2,2)),1,3);
            X(:,2) = atan(R(:,2)./R(:,1))/0.5/pi;
            X(:,1) = acos(R(:,1)./cos(0.5*pi*X(:,2)))/0.5/pi;
            C(:,1) = -sin(20*pi*X(:,1));
            C(:,2) = -cos(20*pi*X(:,2));
            R(any(C>1e-2,2),:) = [];
            R = R + 0.5;
        end
        %% Generate the feasible region
        function R = GetPF(obj)
            [X1,X2] = meshgrid(linspace(0,1,100));
            x = cos(0.5*pi*X1).*cos(0.5*pi*X2);
            y = cos(0.5*pi*X1).*sin(0.5*pi*X2);
            z = sin(0.5*pi*X1);
            fes = all(Constraint(X1(:),X2(:),0.5)<=0,2);
            z(reshape(~fes,size(z))) = nan;
            R = {x+0.5,y+0.5,z+0.5};
        end
    end
end

function PopCon = Constraint(X1,X2,sum1)
    % set the parameters of constraints
    DifficultyFactors = [0.5,0.5,0.5];
    % Type-I parameters
    a = 20;
    b = 2 * DifficultyFactors(1) - 1;
    % Type-II parameters
    d = 0.5;
    if DifficultyFactors(2) == 0.0
        d = 0.0;
    end
    e = d - log(DifficultyFactors(2));
    if isfinite(e) == 0
        e = 1e+30;
    end
    % Type-III parameters
    r = 0.5 * DifficultyFactors(3);
    % Calculate objective values
    PopObj(:,1) = cos(0.5 * pi * X1) .* cos(0.5 * pi * X2) + sum1;
    PopObj(:,2) = cos(0.5 * pi * X1) .* sin(0.5 * pi * X2) + sum1;
    PopObj(:,3) = sin(0.5 * pi * X1) + sum1;
    % Type-I constraints
    PopCon(:,1) = b - sin(a * pi * X1);
    PopCon(:,2) = b - cos(a * pi * X2);
    % Type-II constraints
    PopCon(:,3) = -(e - sum1) .* (sum1 - d);
    if DifficultyFactors(2) == 1.0
        PopCon(:,3) = 1e-4 - abs(sum1 - e);
    end
    % Type-III constraints
    x_k = [1.0, 0.0, 0.0, 1.0 / sqrt(3.0)];
    y_k = [0.0, 1.0, 0.0, 1.0 / sqrt(3.0)];
    z_k = [0.0, 0.0, 1.0, 1.0 / sqrt(3.0)];
    for k=1:length(x_k)
        PopCon(:,3+k) = r * r - ((PopObj(:,1) - x_k(k))).^2  -...
            ((PopObj(:,2) - y_k(k))).^2 - ((PopObj(:,3) - z_k(k))).^2;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%% Example 2.5: Viscous Conservation Law
% In this example, we consider the viscous conservation law
%
% $$ u_t + f(u)_x = \mu u_{xx}, \qquad u(x,0)=u_0(x).$$
%
% A common method to derive efficient methods for this equation is to split
% the evolution into two different operators
%
% $$ S(t): \quad v_t + f(v)_x = 0, \qquad v(x,0)=v_0(x)$$
%
% $$ H(t): \quad w_t = \mu w_{xx}, \qquad w(y,0)=w_0(y)$$
%
% The approximate solution is constructed from the formula
%
% $$u(x,t)\approx [H(\Delta t)\circ S(\Delta t) ]^n u_0(x)$$
%
% This way, one can utilize highly efficient solvers for each of the two
% subequations. Here, we use a spectral difference method for the heat
% equation. Moreover, because the purpose here is simply to demonstrate the
% operator splitting, we use the simple first-order Lax-Friedrichs scheme
% for the hyperbolic conservation law ($r=\Delta t/\Delta x$)
%
% $$u_i^{n+1} = \frac12\bigl(u_{i-1}^n + u_{i+1}^n\bigr)
%     - \frac12 r \bigl[ f(u_{i+1}^n) - f(u_{i-1}^n)\bigr]$$

%% Initial setup
N  = 512; h = 2*pi/N;
T  = 6;
x  = -pi+(0:N)*h; x = 0.5*(x(1:end-1)+x(2:end));
u0 = exp(-4*sin((x+2)/2).^2);

%% Burgers' equation
% A classical example of a viscous conservation law is the so-called
% Burgers' equation
%
% $$ u_t + ( 0.5 u^2)_x = \mu u_{xx}
%
% which can be considered a simplified model of the momentum equation from
% the Navier-Stokes equations. Here, we will consider the evolution of a
% smooth profile for \mu=0.01 using periodic boundary conditions.
nsplit = 50;
u = consheat('flux', u0, x, T, nsplit, 0.01);
subplot(2,1,1), plot(x,u(:,1:10:nsplit+1))
subplot(2,1,2), surf(x,linspace(0,T,nsplit+1),u')
shading interp; view(-10,50), axis tight

%%
% As we see from the figure, the nonlinear convective forces sharpen the
% smooth profile into a shock layer whose width is determined by the
% balance of the diffusion term and the self-sharpening mechanisms in the
% nonlinear flux function.

%% Compare different time steps
% Next, we increase \mu to 0.1 and consider the effect of different number
% of splitting steps.
for n=1:4,
	nsplit = 4*3.^(n-1);
	u = consheat('flux', u0, x, T, nsplit);
   subplot(2,2,n); contourf(x, linspace(0,T,nsplit+1), u');
   xlabel('x'), ylabel('t'), colorbar, title([num2str(nsplit) ' steps']);
end;
%%
% The approximations computed with four and twelve splitting steps clearly
% capture the dominant evolutionary behavior of the solution, but also
% contain nonphysical wiggles caused by splitting errors. As we increase
% the number of splitting steps, the splitting errors decrease and the
% approximate solutions are able to capture the whole evolution without
% significant artifacts created by our splitting strategy.

%% Comparing different shock widths
% As a last example, we will look at the evolution of the smooth profile
% for different balances between the convective and diffusive forces.
mu = 4; nsplit = 50;
t  = linspace(0,T,nsplit+1);
for n=1:4,
   mu = 0.25*mu;
	u  = consheat('flux',u0,x,T,nsplit,mu);
	subplot(2,2,n); contourf(x,t,u');
   xlabel('x'), ylabel('t'), colorbar, title(['mu = ', num2str(mu)]);
end;
%%
% For large values of \mu, the diffusive forces dominate and the solution
% profile evolves toward a flat surface. As \mu decreases, the convective
% forces take over and the solution develops into a viscous shock that
% travels to the right.

repo: Aero-matlab


text:
function gram_vector = find_gram_vector(samplesf, new_sample, num_training_samples, params)
% Find the inner product of the new sample with the existing samples. TO be
% used for distance calculation

% Note: Since getting the 'exact' distance between the samples is not important,
% the inner product computation is done by using only the half spectrum for
% efficiency. In practice the error incurred by this is negligible. Also
% since we wannt to merge samples that are similar, and finding the best
% match is not important, small error in the distance computation doesn't
% matter

gram_vector = inf(params.nSamples,1,'like', params.data_type);

num_feature_blocks = numel(new_sample);

if num_training_samples == params.nSamples
    % This if statement is only for speed
    ip = zeros(1,'like', params.data_type);
    for k = 1:num_feature_blocks
        ip_block = 2*reshape(samplesf{k}, num_training_samples, []) * conj(new_sample{k}(:));
        ip = ip + real(ip_block);
    end

    gram_vector = ip;
elseif num_training_samples > 0
    ip = zeros(1,'like', params.data_type);
    for k = 1:num_feature_blocks
        ip_block = 2*reshape(samplesf{k}(1:num_training_samples,:,:,:),num_training_samples, []) * conj(new_sample{k}(:));
        ip = ip + real(ip_block);
    end

    gram_vector(1:num_training_samples) = ip;
end
end

author: martin-danelljan

repo: ECO


text:
%  Figure 7.18      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
% fig7_18.m
% script to plot the step position responses of the tape
% drive servo with dominant second order, and LQR designs
clf;

f =[0    2.0000         0         0         0;
   -0.1000   -0.3500    0.1000    0.1000    0.7500;
         0         0         0    2.0000         0 ;
    0.4000    0.4000   -0.4000   -1.4000         0  ;
         0   -0.0300         0         0   -1.0000   ];

g =[0;
     0;
     0 ;
     0  ;
     1   ];
h3 =[0.5000         0    0.5000         0         0];
ht =[-0.2000   -0.2000    0.2000    0.2000         0];

p2 =[-0.7070+0.7070*i;
  -0.7070-0.7070*i    ;
  -4.0000               ;
  -4.0000                ;
  -4.0000                 ];

j=0;
p22 =p2/1.5;
k2=acker(f,g,p22);

k2 =[8.5123   20.3457   -1.4911   -7.8821    6.1927];

f2c=f-g*k2;

s=[f g;h3 j];
r=[0 0 0 0 0 1]';
n=s\r;
nu=n(6)
nx=n(1:5)
nbar2=k2*nx+nu
t=0:.2:12;
sys2=ss(f2c,g*nbar2,h3,j);
y2=step(sys2,t);
% lqr
kqr=lqr(f,g,h3'*h3,1)

klqr =[0.6526    2.1667    0.3474    0.5976    1.0616];

flqrc=f-g*klqr;
nbarlqr = klqr*nx+nu;

syslqr=ss(flqrc,g*nbarlqr,h3,j);
ylqr=step(syslqr,t);

plot(t,y2,t,ylqr,'-.','LineWidth',2);
xlabel('Time (msec)');
ylabel('Tape position, x_3');
text(5,1.1, 'LQR');
text(4.5,.9,'Dominant second-order');
title('Fig. 7.18: Step responses of tape servo designs');
nicegrid;

author: Sable

repo: mcbench-benchmarks


text:
function out=tanh(x)

out=sinh(x)./cosh(x);

author: opencobra

repo: cobratoolbox


text:
function [zd,stat]=zdres(flag,obs,sv,nav,rr,opt,rtk)

global glc
nobs=size(obs,1); nf=rtk.NF;
zazel=[0 pi/2]; stat=1;
zd0=struct('y',zeros(2*nf,1),'LOS',zeros(1,3),'azel',zeros(1,2));
zd=repmat(zd0,nobs,1);

if norm(rr)<=0,stat=0;return;end

%compute earth tide correction
if opt.tidecorr
    dr=tidedisp(gpst2utc(obs(1).time),rr,7,nav.erp,nav.otlp(:,:,flag));
    rr=rr+dr;
end

pos=ecef2pos(rr);

for i=1:nobs

    rs = sv(i).pos; dts = sv(i).dts;

    [r,LOS]=geodist(rs,rr); azel=satazel(pos,LOS);
    if r<=0||azel(2)<opt.elmin,continue;end

    if satexclude(obs(i).sat,sv(i).vars,sv(i).svh)==0,continue;end

    r=r-dts*glc.CLIGHT;

    zhd=trop_saa(pos,zazel,0);
    r=r+tropmapf(obs(i).time,pos,azel)*zhd;

    dantr=antmodel(obs(i).sat,opt.pcvr(flag),opt.antdel(flag,:),azel,opt.posopt(2));

    zd(i)=zdres_sat(flag,r,obs(i),nav,azel,dantr,opt,rtk,zd(i));
    zd(i).LOS=LOS; zd(i).azel=azel;
end

return

repo: GINav


text:
function f = fint(x)
%fint   subdivision function
%   f = fint(x);
%   called by fitint
%   IFISS function: DJS; 28 February 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
      global global_N global_INTL global_LASTDL
      f=(global_INTL*(x-1.0)-global_LASTDL)*(x^global_N) + global_LASTDL;
%     fprintf('\n %10.4e   %10.4e',x,f)
      return


text:
%% string get_algorithm_name_by_number(int)
% algorithm_id - int
% algorithm_name - string
%
function [algorithm_name] = get_algorithm_name_by_id(algorithm_id)
  %clear all; clc; algorithm_id = 'FPCP';
  L = load_algorithm_list();
  [m,n] = size(L{1,1});
  algorithm_name = '';
  for i = 1 : m
    a_id = strtrim(char(L{1,2}(i)));
    a_name = strtrim(char(L{1,3}(i)));
    %disp([algorithm_id '<->' a_id]);
    if(strcmp(algorithm_id,a_id))
      algorithm_name = a_name;
      break;
    end
  end
end

author: andrewssobral

repo: lrslibrary


text:
function poseRefPoint = getPoseRefPoint(keypointsAll,pidxsHead)

poseRefPoint = zeros(size(keypointsAll,2),2);
for clusidx = 1:size(keypointsAll,2)
    pp = [];
    % reference point is head center
%     pidxHC = 8; % head center
%     pidxN = 15; % neck
%     pidxHT = 16; % head top
    for pidx = pidxsHead
        pp = [pp; keypointsAll{pidx,clusidx}];
    end
    if (isempty(pp))
        pp = [-inf -inf];
    end
    poseRefPoint(clusidx,:) = mean(pp,1);
end

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
classdef WFG1 < PROBLEM
% <multi/many> <real> <large/none> <expensive/none>
% Benchmark MOP proposed by Walking Fish Group
% K --- --- The position parameter, which should be a multiple of M-1

%------------------------------- Reference --------------------------------
% S. Huband, P. Hingston, L. Barone, and L. While, A review of
% multiobjective test problems and a scalable test problem toolkit, IEEE
% Transactions on Evolutionary Computation, 2006, 10(5): 477-506.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties(Access = private)
        K;  % Position parameter
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            if isempty(obj.M); obj.M = 3; end
            obj.K = obj.ParameterSet(obj.M-1);
            if isempty(obj.D); obj.D = obj.K + 10; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = 2 : 2 : 2*obj.D;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            [N,D] = size(PopDec);
            M = obj.M;
            K = obj.K;
            L = D - K;
            D = 1;
            S = 2 : 2 : 2*M;
            A = ones(1,M-1);

            z01 = PopDec./repmat(2:2:size(PopDec,2)*2,N,1);

            t1 = zeros(N,K+L);
            t1(:,1:K)     = z01(:,1:K);
            t1(:,K+1:end) = s_linear(z01(:,K+1:end),0.35);

            t2 = zeros(N,K+L);
            t2(:,1:K)     = t1(:,1:K);
            t2(:,K+1:end) = b_flat(t1(:,K+1:end),0.8,0.75,0.85);

            t3 = zeros(N,K+L);
            t3 = b_poly(t2,0.02);

            t4 = zeros(N,M);
            for i = 1 : M-1
                t4(:,i) = r_sum(t3(:,(i-1)*K/(M-1)+1:i*K/(M-1)),2*((i-1)*K/(M-1)+1):2:2*i*K/(M-1));
            end
            t4(:,M) = r_sum(t3(:,K+1:K+L),2*(K+1):2:2*(K+L));

            x = zeros(N,M);
            for i = 1 : M-1
                x(:,i) = max(t4(:,M),A(i)).*(t4(:,i)-0.5)+0.5;
            end
            x(:,M) = t4(:,M);

            h      = convex(x);
            h(:,M) = mixed(x);
            PopObj = repmat(D*x(:,M),1,M) + repmat(S,N,1).*h;
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            M = obj.M;
            R = UniformPoint(N,M);
            c = ones(size(R,1),M);
            for i = 1 : size(R,1)
                for j = 2 : M
                    temp = R(i,j)/R(i,1)*prod(1-c(i,M-j+2:M-1));
                    c(i,M-j+1) = (temp^2-temp+sqrt(2*temp))/(temp^2+1);
                end
            end
            x = acos(c)*2/pi;
            temp = (1-sin(pi/2*x(:,2))).*R(:,M)./R(:,M-1);
            a = 0 : 0.0001 : 1;
            E = abs(temp*(1-cos(pi/2*a))-1+repmat(a+cos(10*pi*a+pi/2)/10/pi,size(x,1),1));
            [~,rank] = sort(E,2);
            for i = 1 : size(x,1)
                x(i,1) = a(min(rank(i,1:10)));
            end
            R      = convex(x);
            R(:,M) = mixed(x);
            R      = repmat(2:2:2*M,size(R,1),1).*R;
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            if obj.M == 2
                R = obj.GetOptimum(100);
            elseif obj.M == 3
                a = linspace(0,pi/2,10)';
                x = (1-cos(a))*(1-cos(a'));
                y = (1-cos(a))*(1-sin(a'));
                z = 1 - a*ones(size(a'))*2/pi - cos(20*a*ones(size(a'))+pi/2)/10/pi;
                R = {x*2,y*4,z*6};
            else
                R = [];
            end
        end
    end
end

function Output = s_linear(y,A)
    Output = abs(y-A)./abs(floor(A-y)+A);
end

function Output = b_flat(y,A,B,C)
    Output = A+min(0,floor(y-B))*A.*(B-y)/B-min(0,floor(C-y))*(1-A).*(y-C)/(1-C);
    Output = round(Output*1e4)/1e4;
end

function Output = b_poly(y,a)
    Output = y.^a;
end

function Output = r_sum(y,w)
    Output = sum(y.*repmat(w,size(y,1),1),2)./sum(w);
end

function Output = convex(x)
    Output = fliplr(cumprod([ones(size(x,1),1),1-cos(x(:,1:end-1)*pi/2)],2)).*[ones(size(x,1),1),1-sin(x(:,end-1:-1:1)*pi/2)];
end

function Output = mixed(x)
    Output = 1-x(:,1)-cos(10*pi*x(:,1)+pi/2)/10/pi;
end


text:
function [ feature ] = block_dct( im)
%
%   I: INPUT IMAGE
%   nSc: number of scale
%   feature: Returned feature on Block_DCT statistics

% shape parameter estimation on blocks and pool them together,
% select the first 10th,last 10th and mean as featurs

gama_L1 = blkproc(im,[3,3],[2,2],@gama_dct);
gama_sorted_temp = sort(gama_L1(:),'ascend');
gama_count = length(gama_sorted_temp);
p10_gama_L1=mean(gama_sorted_temp(1:ceil(gama_count*0.1)));
% p10_last_gama_L1=mean(gama_sorted_temp(fix(gama_count*0.9):end));
p100_gama_L1=mean(gama_sorted_temp(:));
clear gama_sorted_temp gama_count

feature=[p10_gama_L1 p100_gama_L1];

% coefficient variation estimation on blocks and pool them together,
% select the first 10th,last 10th and mean as featurs
coeff_var_L1 = blkproc(im,[3,3],[2,2],@coeff_var_dct);
cv_sorted_temp = sort(coeff_var_L1(:),'ascend');
cv_count = length(cv_sorted_temp);
% p10_cv_L1=mean(cv_sorted_temp(1:ceil(cv_count*0.1)));
p10_last_cv_L1=mean(cv_sorted_temp(fix(cv_count*0.9):end));
p100_cv_L1=mean(cv_sorted_temp(:));
clear cv_sorted_temp cv_count

feature=[feature p10_last_cv_L1 p100_cv_L1];

ori1_rho_L1 = blkproc(im,[3,3],[2,2],@oriented1_dct_rho_config3);
ori2_rho_L1 = blkproc(im,[3,3],[2,2],@oriented2_dct_rho_config3);
ori3_rho_L1 = blkproc(im,[3,3],[2,2],@oriented3_dct_rho_config3);
temp_size=size(ori1_rho_L1);
var_temp=zeros(temp_size);

for i=1:temp_size(1)
    for j=1:temp_size(2)
        var_temp(i,j)=var([ori1_rho_L1(i,j) ori2_rho_L1(i,j) ori3_rho_L1(i,j)]);
   end
end
ori_rho_L1=var_temp;

ori_sorted_temp = sort(ori_rho_L1(:),'ascend');
ori_count = length(ori_sorted_temp);
% p10_orientation_L1=mean(ori_sorted_temp(1:ceil(ori_count*0.1)));
p10_last_orientation_L1=mean(ori_sorted_temp(fix(ori_count*0.9):end));
p100_orientation_L1=mean(ori_sorted_temp(:));
clear var_ori_sorted_temp rho_count

feature=[feature p10_last_orientation_L1 p100_orientation_L1];

end

repo: sr-metric


text:
function [X,Y,Z]=cone0(z,t,p)
% [X,Y,Z]=cone0 defines a non-orthogonal
% coordinate system using a conical surface
% with a planar base
X=z.*tan(t).*cos(p); Y=z.*tan(t).*sin(p); Z=z;


text:
% Earth Observing System Data Visualization
%
%   Part 2 (Using the Mapping Toolbox)
%
%   This M-file contains the sequence of MATLAB commands listed in
%   the November 2002 MATLAB Digest article on "Reading and Visualizing
%   Data from the Earth Observing System (EOS)," from the section that
%   requires the Mapping Toolbox.
%
%   Rob Comer and Chris Lawton
%   Copyright 2002 The MathWorks, Inc.

%== MAP PROJECTION OF COMBINED SST AND NDVI ==============================

% Set up figure, title, colormap, and map axes. Use a Mollweide projection
% with Greenwich as its central meridian.
figure('Color','k');
title({'SST from MODIS and NDVI from AVHRR, July 2001'},'Color','w')
colormap(cmap);
ax = axesm('MapProjection','mollweid','Origin',[0 0 0]);
set(ax,'Color',[0 0 0])

% Display the global NDVI on the map axes. The CDataMapping property
% must be set to 'direct' so that the data values directly index into the
% colormap. Add 1 because indexing for uint8s is zero-based but indexing
% for doubles is one-based. FLIPUD is required because columns in a
% "regular matrix map" must run from south to north.
ndviNorthwestCorner = [89.5 -179.5];  % Lat/lon of upper-left-most cell center
ndviCellSize = 1;                     % Size of grid cells in degrees
hmesh = meshm(flipud(double(ndviMapped) + 1),...  % Converting uint8 image
              [ndviCellSize ndviNorthwestCorner],... % The 'legend'
              'CDataMapping','direct');  % Go directly into the colormap

% Display the MODIS SST swath over the global NDVI grid.
hsurf = surfm(sstLat,sstLon,double(sstMapped) + 1,...
              'CDataMapping','direct');

% Show 2 color scale bars, one for SST and one for NDVI.
% Move the map to make room, then add the color scales.
set(ax,'Position',get(ax,'Position') + [0 0.12 0 0])
colorscale(sstCmapLim,sstDataLim,5,'horiz',...
    'Position',[0.1 0.10, 0.8, 0.03],'XColor','w','YColor','w')
title('Sea Surface Temperature','Color','w')
xlabel('degrees Celsius')
colorscale(ndviCmapLim,ndviDataLim,0.2,'horiz',...
    'Position',[0.1 0.28, 0.8, 0.03],'XColor','w','YColor','w')
title('Normalized Difference Vegetation Index','Color','w')


text:
classdef GridOptions < handle
    properties
        dx              double {isfinite}                   % east-west distance between grid points
        dy              double                              % north-south distance between grid points
        dz              double                              % vertical distance between grid points
        % Defines whether horizontal distances are constant, or whether they scale as the grid
        % deviates from the equator.
        followMeridians     matlab.lang.OnOffSwitchState  	= 'off'
        gridEntireArea      matlab.lang.OnOffSwitchState    = 'off'
        FixedAnchorPoint    double                          = []
        % grid cannot be used past these limits [xmin xmax ymin ymax]
        AbsoluteGridLimits  (1,4) double
        GridType            GridTypes                       = GridTypes.XY;
        RefEllipsoid        referenceEllipsoid
    end

    properties(Dependent, Hidden)
        horizUnits
        dzUnits
    end

    methods
        function obj = GridOptions(grid_type, dx_dy_dz, RefEllipsoid, varargin)
            % old usage: obj = GridOptions(dx, dy, dz_km, horiz_units, follow_meridians, gridEntireArea)
            % GRIDOPTIONS defines parameters that are used to create ZmapGrids
            %
            % obj = GRIDOPTIONS( GRID_TYPE, [dx,dy,dz], ELLPSOID ) where GRID_TYPE is a GridType,
            % and defines the spatial orientation of the grid.  [dx,dy,dz] are a doublet or triplet
            % describing the grid spacing, and ELLIPSOID defines how everything will be interpreted.
            % the ellipsoid's LengthUnit determines the units of this grid.  (See the option for
            % FollowMeridians, below). An specifying the Ellipsoid value of nonEllipsoid() will
            % declare that this is a cartesian grid instead of a geodetic one.
            %
            % obj = GRIDOPTIONS( ..., 'FollowMeridians',true) Horizontal Units are degrees and grid
            % therefore narrows toward the poles.
            %
            % obj = GRIDOPTIONS(...,'GridEntireArea',true) grids entire area.
            %
            % For a horizontal grid, the first parameter should be [dx, dy] or [dx,dy,nan]
            % For a cross-sectional grid, the first parameter should be [dS nan dz], where dS is
            % a linear distance.
            % For a 3-d grid, the first parameters should be [dx,dy,dz]
            %
            %  see also ZmapGrid, referenceEllipsoid, GridType, nonEllipsoid
            p = inputParser;
            p.addRequired('GridType');
            p.addRequired('dx_dy_dz');
            p.addRequired('RefEllipsoid');
            p.addParameter('FollowMeridians', false);
            p.addParameter('GridEntireArea', false);
            p.parse(grid_type, dx_dy_dz, RefEllipsoid,varargin{:})

            obj.GridType        = p.Results.GridType;
            obj.RefEllipsoid    = p.Results.RefEllipsoid;
            obj.followMeridians = p.Results.FollowMeridians;
            obj.gridEntireArea  = p.Results.GridEntireArea;

            if iscartesian(obj.RefEllipsoid)
                assert(p.Results.FollowMeridians == false, 'Cannot follow Meridians for a cartesian grid');
            end

            switch obj.GridType
                case 'XY'
                    obj.dx = dx_dy_dz(1);
                    obj.dy = dx_dy_dz(2);
                    obj.dz = 1;

                case 'XZ'
                    obj.dx = dx_dy_dz(1);
                    obj.dy = nan;
                    obj.dz = dx_dy_dz(end);

                case 'XYZ'
                    obj.dx = dx_dy_dz(1);
                    obj.dy = dx_dy_dz(2);
                    obj.dz = dx_dy_dz(3);
            end
        end

        function u = get.horizUnits(obj) % for backwards compatibility
            if obj.followMeridians
                u = 'degree';
            else
                u = obj.RefEllipsoid.LengthUnit;
            end
        end

        function u = get.dzUnits(obj) % for backwards compatibility
            u = obj.RefEllipsoid.LengthUnit;
        end

    end % methods section

    methods(Static)
        function [mygrid, mygridopts] = fromDialog(existing_gridopt, ellipsoid, shape)
            % FROMDIALOG shows an interactive dialog box allowing user to choose grid
            mygrid = [];
            mygridopts = [];
            if ~exist('ellipsoid', 'var')
                ellipsoid = getappdata(groot, 'ZmapDefaultReferenceEllipsoid');
            end

            if ~exist('shape', 'var')
                shape = ShapeGeneral();
            end
            if exist('existing_gridopt', 'var') && ~isempty(existing_gridopt)
                gc = grid_chooser(ellipsoid, existing_gridopt, shape);
            else
                gc = grid_chooser(ellipsoid, [], shape);
            end
            gc.ResultDump = @set_values;
            waitfor(gc)
            function set_values(g, gop)
                mygrid = g;
                mygridopts = gop;
            end
        end
    end %static methods
end

author: CelsoReyes


text:
function d = dsphere ( p, xc, yc, zc, r )

%*****************************************************************************80
%
%% DSPHERE returns the signed distance of one or more points to a sphere.
%
%  Discussion:
%
%    The corresponding routine in 2D is called DCIRCLE.
%
%  Licensing:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, real P(NP,3), the point coordinates.
%
%    Input, real XC, YC, ZC, the coordinates of the center of the sphere.
%
%    Input, real R, the radius of the sphere.
%
%    Output, real D(NP), the signed distance of each point to the
%    sphere.  The point is inside, on, or outside the sphere depending
%    on whether D is negative, zero, or positive.
%
  d = sqrt ( ( p(:,1) - xc ).^2 ...
           + ( p(:,2) - yc ).^2 ...
           + ( p(:,3) - zc ).^2 ) - r;

  return
end