text:
%% DEMO_febio_00014_cube_varying_material
% 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
% * displacement control, displacement boundary condition
% * hexahedral elements, hex8
% * cube, box, rectangular
% * static, solid
% * spatially varying material properties
% * hyperelastic, Ogden
% * displacement logfile
% * stress logfile

%%

clear; close all; clc;

%% Plot settings
fontSize=15;
faceAlpha1=0.8;
markerSize=40;
lineWidth=3;
cMap=spectral(250);

%% 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 force
febioLogFileName_sed=[febioFebFileNamePart,'_sed_out.txt']; %Log file name for exporting strain energy density

%Specifying dimensions and number of elements
pointSpacings=1*ones(1,3); %Desired point spacing between nodes
cubeSize=10;
sampleWidth=cubeSize; %Width
sampleThickness=cubeSize; %Thickness
sampleHeight=cubeSize; %Height
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
appliedStrain=0.3; %Linear strain (Only used to compute applied stretch)
loadingOption='tension'; % or 'compression'
switch loadingOption
    case 'compression'
        stretchLoad=1-appliedStrain; %The applied stretch for uniaxial loading
    case 'tension'
        stretchLoad=1+appliedStrain; %The applied stretch for uniaxial loading
end
displacementMagnitude=(stretchLoad*sampleHeight)-sampleHeight; %The displacement magnitude

%Material parameter sets
testOpt=2; %1=Linear gradient of material, or 2=gyroid based material distribution
E_youngs_min=1e-3; %Lowest Youngs modulus
E_youngs_max=1; %Highest Youngs modulus
nu_min=0.3; %Lowest Poissons ratio
nu_max=0.45; %Lowest Poissons ratio

% 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
runMode='external';

%% 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

%% Define spatially varying material distribution data
% Here the element centre coordinates are used assign the material
% stiffness based on a particular function on these coordinates.

VE=patchCentre(E,V); %Element centres
switch testOpt
    case 1 %linear gradient in X direction
        S=VE(:,1);
    case 2 %gyroid
        %Scale coordinates for gyroid
        VE=VE-min(VE(:));
        VE=VE./max(VE(:));
        VE=(VE.*2*pi)-pi;

        %Evaluate gyroid
        S=triplyPeriodicMinimal(VE(:,1),VE(:,2),VE(:,3),'g');
end

%Normalize data
S=S-min(S(:)); %Subtract minimum -> range [0-...]
S=S./max(S(:)); %Devide by max -> range [0-1]

%Use scaling data S to generate element Youngs moduli
E_youngs_elem=S.*(E_youngs_max-E_youngs_min)+E_youngs_min;
nu_elem=S.*(nu_max-nu_min)+nu_min;

%Fix mesh struct for plotting
meshStruct.elements=E;
meshStruct.elementData=E_youngs_elem;

%%
% 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(gca,gjet(6)); icolorbar;
axisGeom(gca,fontSize);

hs=subplot(1,2,2); hold on;
title('Cut view of solid mesh and materials','FontSize',fontSize);
optionStruct.hFig=[hFig hs];
meshView(meshStruct,optionStruct);
colormap(gca,cMap); colorbar; caxis([E_youngs_min E_youngs_max]);
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
bcSupportList=unique(Fb(Cb==5,:)); %Node set part of selected face

%Prescribed displacement nodes
bcPrescribeList=unique(Fb(Cb==6,:)); %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 full support','BC z 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';

%Control section
febio_spec.Control.analysis='STATIC';
febio_spec.Control.time_steps=numTimeSteps;
febio_spec.Control.step_size=1/numTimeSteps;
febio_spec.Control.solver.max_refs=max_refs;
febio_spec.Control.solver.max_ups=max_ups;
febio_spec.Control.time_stepper.dtmin=dtmin;
febio_spec.Control.time_stepper.dtmax=dtmax;
febio_spec.Control.time_stepper.max_retries=max_retries;
febio_spec.Control.time_stepper.opt_iter=opt_iter;

%Material section
materialName1='Material1';
dataMapName1='MaterialParameterMap1';
dataMapName2='MaterialParameterMap2';
febio_spec.Material.material{1}.ATTR.name=materialName1;
febio_spec.Material.material{1}.ATTR.type='neo-Hookean';
febio_spec.Material.material{1}.ATTR.id=1;
febio_spec.Material.material{1}.E.ATTR.type='map'; %Calls for mapping of parameter
febio_spec.Material.material{1}.E.VAL=dataMapName1; %Calls for mapping of parameter
febio_spec.Material.material{1}.v.ATTR.type='map'; %Calls for mapping of parameter
febio_spec.Material.material{1}.v.VAL=dataMapName2; %Calls for mapping of parameter

% 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(:);

%MeshData secion
%-> Element data
febio_spec.MeshData.ElementData{1}.ATTR.name=dataMapName1;
febio_spec.MeshData.ElementData{1}.ATTR.elem_set=partName1;
febio_spec.MeshData.ElementData{1}.elem.ATTR.lid=(1:1:size(E,1))';
febio_spec.MeshData.ElementData{1}.elem.VAL=E_youngs_elem;

febio_spec.MeshData.ElementData{2}.ATTR.name=dataMapName2;
febio_spec.MeshData.ElementData{2}.ATTR.elem_set=partName1;
febio_spec.MeshData.ElementData{2}.elem.ATTR.lid=(1:1:size(E,1))';
febio_spec.MeshData.ElementData{2}.elem.VAL=nu_elem;

%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';

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

febio_spec.Boundary.bc{3}.ATTR.type='prescribe';
febio_spec.Boundary.bc{3}.ATTR.node_set=nodeSetName2;
febio_spec.Boundary.bc{3}.dof='z';
febio_spec.Boundary.bc{3}.scale.ATTR.lc=1;
febio_spec.Boundary.bc{3}.scale.VAL=displacementMagnitude;
febio_spec.Boundary.bc{3}.relative=0;

%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];

%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_sed;
febio_spec.Output.logfile.element_data{1}.ATTR.data='sed';
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
% febView(febioFebFileName);

%% 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=runMode;

[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_sed),1,1);

    %Access data
    E_sed_mat=dataStruct.data;

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

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

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);
    title('$\Psi$ $[J/m^3]$','Interpreter','Latex')

    hp=gpatch(Fb,V_DEF(:,:,end),CV,'k',1); %Add graphics object to animate
    hp.FaceColor='interp';

    axisGeom(gca,fontSize);
    colormap(cMap); colorbar;
    caxis([min(E_sed_mat(:)) max(E_sed_mat(:))]/3);
    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_sed_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 [image, imageCha, objOut] = CS_reconstruction(path, para)
% CS_RECONSTRUCTION main function to select the appropriate Compressed
% Sensing reconstruction
%
% input:
% path      - path to file without specifying file extension
%           - path to directory containing files
%           - empty (for file dialog)
%           - kSpace to be reconstructed with dimensionalities:
%             cell array: NSlice x NChannels x NRepetitions x NAverages
%             each cell: 2D: yPhase x xFreq x (nTime/nPhases)
%                        3D: yPhase x xFreq x zPhase
%                        4D: yPhase x xFreq x zPhase x nTime/nPhases
% para      - struct containing to be modified reconstruction parameters
%           and/or struct with handover parameter
%           - path to parameter file
%
%
% output:
% image     reconstructed and sum-of-squared channel combined image
% imageCha  reconstructed individual coil images
% objOut    reconstruction parameters
%
%
% (c) copyright 2012 - 2016 under BSD license
% Thomas Kuestner, University of Tuebingen and University of Stuttgart, Germany
% (thomas.kuestner@{med.uni-tuebingen.de,iss.uni-stuttgart.de})
%
% Please see license files (./licenses) for accompanied codes (which also
% include own modifications) of:
% - ESPIRiT v0.1.8  in @ESPIRiT, utils/utils_ESPIRiT     under LICENSE_ESPIRiT
% - bart v0.2.09    in utils/utils_BART                  under LICENSE_BART
% - SPIRiT v0.3     in @SPIRiT, utils/utils_GRAPPA,
%                   utils/utils_SPIRiT                   under LICENSE_SPIRiT
% - sparseMRI v0.2  in @sparseMRI, @p2DFT, @A_operator,
%                   @TVOP                                under LICENSE_SPARSEMRI
% - L1Magic v1.11   in @L1_Magic                         under LICENSE_L1MAGIC
% - GIST            in utils/utils_Proximal/GIST         under LICENSE_GIST
% - Wavelab850      in @Wavelet, utils/utils_Wavelet     under LICENSE_WAVELAB
% - Rice Wavelet v3.0 in utils/utils_WaveletRice         under LICENSE_WAVELETRICE
% - Curvelab v2.1.3 in utils/utils_TRAFO/CurveLab-2.1.3  under LICENSE_CURVELAB
% - NUFFT           in @NUFFT, utils/utils_NUFFT         under LICENSE_NUFFT
% - DTCWT           in utils/utils_TRAFO/DTCWT           under LICENSE_DTCWT
% - Shearlet        in utils/utils_TRAFO/Shearlet        under LICENSE_SHEARLET
% - Denoising       in utils/utils_Proximal/Denoising    under LICENSE_DENOISING

warning('off','MATLAB:dispatcher:nameConflict');
warning('off','MATLAB:mat2cell:TrailingUnityVectorArgRemoved');

% global prop

% prepare m-files
currpath = fileparts(mfilename('fullpath'));
addpath(genpath([currpath,filesep,'utils']));
addpath(genpath([currpath,filesep,'io']));
addpath(genpath([currpath,filesep,'preproc']));
addpath(genpath([currpath,filesep,'postproc']));

% load ADCs
if(nargin == 0 || ~exist('path','var') || isempty(path))
    if(strcmp(getenv('computername'),'M123PC'))
        defpath = 'K:\Compressed Sensing\Rohdaten\3D';
    else
        defpath = 'S:\Rohdaten\Thomas\ktFOCUSS\3D';
    end
    [file, path] = uigetfile({'*.mat;*.dat;', 'ADC measdata (*.mat, *.dat)'; '*.h5', 'ISMRMD (*.h5)'} ,'Select measdata file', 'MultiSelect', 'off', defpath);
    if(isequal(file,0))
        error('CS_reconstruction(): User Termination');
    end
    [~,filename,ext] = fileparts(file);
    lReadFromFile = true;
    % in general: [dData, iLC] = fMeasRead(datapath, 'Set', 0, 'Smp', [21 inf]);
else
    if(ischar(path))
        if(isdir(path))
            [file, path] = uigetfile({'*.mat;*.dat;', 'ADC measdata (*.mat, *.dat)'; '*.h5', 'ISMRMD (*.h5)'} ,'Select measdata file', 'MultiSelect', 'off', path);
            if(isequal(file,0))
                error('CS_reconstruction(): User Termination');
            end
            [~,filename,ext] = fileparts(file);
        else
            [path, filename, ext] = fileparts(path);
        end
        lReadFromFile = true;
    else
        lReadFromFile = false;
    end
end

% read inputs
if(lReadFromFile)
    if(strcmp(ext,'h5')) % ISMRMD
        data = h5read([path,filesep,filename,ext],'dataset/data');
        % convert to processable kSpace
        [ kSpace, measPara ] = convertKSpace_main( data );
    else % ADC measdata (Siemens)
        if(exist(fullfile(path,[filename,'.mat']),'file'))
            matfile = whos('-file', fullfile(path,[filename,'.mat']));
        end
        if(exist(fullfile(path,[filename,'.mat']),'file') && ~ismember('iLC', {matfile.name}))
            % load mat file to check if it contains iLC/iSP or something else
            data = load(fullfile(path,[filename,'.mat']));
            kSpace = data.kSpace;
            measPara.dim = zeros(1,5);
            measPara.LCall = zeros(1,4);
            dimNames = {'y (phase)', 'x (frequency)', 'z (phase)', 't (time)', 'channels'};
            LCallNames = {'slices', 'acquisitions', 'echos', 'repetitions'};
            dimDefault = ones(1,5);
            LCallDefault = [size(kSpace,1), size(kSpace,4), size(kSpace,5), size(kSpace,3)];
            tmp = size(kSpace{1});
            if(LCallDefault(1) > 1 && length(tmp) > 2)
                dimDefault(1:2) = tmp(1:2);
                dimDefault(4) = tmp(3);
            else
                dimDefault(1:length(tmp)) = tmp;
            end
            dimDefault(5) = size(kSpace,2);
            if(isfield(data,'dim') && isfield(data,'LCall'))
                measPara.dim = data.dim;
                measPara.LCall = data.LCall;
            elseif(isfield(data,'dim') && ~isfield(data,'LCall'))
                measPara.dim = data.dim;
                fprintf('Please specify amount of:\n');
                for i=1:length(measPara.LCall)
                    helper = input(sprintf('%s [%d]: ',LCallNames{i}, LCallDefault(i)),'s');
                    if(isempty(helper)), helper = LCallDefault(i); end;
                    measPara.LCall(i) = helper;
                end

            elseif(~isfield(data,'dim') && isfield(data,'LCall'))
                measPara.LCall = data.LCall;
                fprintf('Please specify dimensions:\n');
                for i=1:length(measPara.dim)
                    helper = input(sprintf('%s [%d]: ',dimNames{i}, dimDefault(i)),'s');
                    if(isempty(helper)), helper = dimDefault(i); end;
                    measPara.dim(i) = helper;
                end

            else
                fprintf('Please specify dimensions:\n');
                for i=1:length(measPara.dim)
                    helper = input(sprintf('%s [%d]: ',dimNames{i}, dimDefault(i)),'s');
                    if(isempty(helper)), helper = dimDefault(i); end;
                    measPara.dim(i) = helper;
                end
                fprintf('Please specify amount of:\n');
                for i=1:length(measPara.LCall)
                    helper = input(sprintf('%s [%d]: ',LCallNames{i}, LCallDefault(i)),'s');
                    if(isempty(helper)), helper = LCallDefault(i); end;
                    measPara.LCall(i) = helper;
                end
            end
            clear 'LCallDefault' 'dimDefault';
            loadedVars = {'dimension', 'oversampling', 'measPara', 'iLC', 'iLCPositions', 'drecksMDH'};
            for i=1:length(loadedVars)
                if(isfield(data,loadedVars{i}))
                    if(any(strcmp(loadedVars{i},{'dimension','oversampling'})))
                        eval(sprintf('measPara.%s = data.%s;', loadedVars{i}, loadedVars{i}));
                    else
                        eval(sprintf('%s = data.%s;', loadedVars{i}, loadedVars{i}));
                    end
                else
                    if(strcmp(loadedVars{i},'measPara'))
                        continue;
                    end
                    eval(sprintf('%s = [];', loadedVars{i}));
                end
            end
            clear 'data';
        else
            [dData, iLC, iEvalInfoMask, drecksMDH, iLCPositions] = fMeas_main(fullfile(path, [filename, ext]));
            measPara = [];
        end
    end
else
    % kSpace was directly handed over
    kSpace = path;
    path = currpath; % for compatibility later
    filename = '';
    if(isstruct(para))
        loadedVars = {'measPara', 'dimension', 'oversampling', 'iLC', 'iLCPositions', 'drecksMDH', 'savePath'};
        for i=1:length(loadedVars)
            if(isfield(para,loadedVars{i}))
                if(any(strcmp(loadedVars{i},{'dimension','oversampling'})))
                    eval(sprintf('measPara.%s = para.%s;', loadedVars{i}, loadedVars{i}));
                else
                    eval(sprintf('%s = para.%s;', loadedVars{i}, loadedVars{i}));
                end
                eval(sprintf('para = rmfield(para,''%s'');', loadedVars{i}));
            else
                eval(sprintf('%s = [];', loadedVars{i}));
            end
        end
    end
end

%% save output
prop.flagSave = false;

%% show output in imagine
prop.flagPlot = true;

%% display output
prop.disp.names = {'Repetitions', 'Averages', 'Slices', 'Channels', 'Time', 'Frequency', 'Extracting K-Space', 'Phase Correction', 'GRAPPA calibration', 'FOCUSS', 'CG', 'ESPReSSo', 'Trafo', 'POCS', 'MC', 'Line Search', 'Log barrier', 'Newton', 'Proximal Average', 'ADMM', 'SENSE - Slice', 'SENSE - Time', 'SENSE - Calibration', 'SENSE - Kernel', 'SENSE - Concatenate'};
if(usejava('jvm') && ~feature('ShowFigureWindows'))
    % use text-based alternative
    prop.flagDisp = false;
    prop.flagPlot = prop.flagPlot & false;

    % initialize some console variables
    prop.maxvals = zeros(1,length(prop.disp.names));
    prop.lastOut = '';
else
    % use GUI dialogs
    prop.flagDisp = true;
    prop.flagPlot = prop.flagPlot & true;

    % initialize some GUI variables
    prop.disp.openBar = false(1,size(prop.disp.names,2));
    prop.disp.colors = colorGradient(rgb('DarkRed'),rgb('ForestGreen'),100);
    prop.disp.ppm = cell(1,length(prop.disp.names));
    % initialize some console variables
    prop.maxvals = zeros(1,length(prop.disp.names));
    prop.lastOut = '';
end

%% parallel computing
prop.flagParallel = false;
prop.openPool = false;
helper = ver;
lFound = false;
for i=1:length(helper)
    if(strcmp(helper(i).Name,'Parallel Computing Toolbox'))
        lFound = true;
        prop.flagParallel = prop.flagParallel & true;
        if(verLessThan('matlab','8.4'))
            if(matlabpool('size') > 0)
    %             matlabpool close
                prop.openPool = true;
            end
        else
            if(~isempty(gcp('nocreate')))
                prop.openPool = true;
            end
        end
        break;
    end
end
if(~lFound), prop.flagParallel = false; end

%% load reconstruction parameter
if(nargin > 1 && exist('para', 'var') && ischar(para))
    [pathPara, filenamePara, ~] = fileparts(para);
    cd(pathPara);
    eval([filenamePara,';']);
    cd(currpath);
else
    % load default parameter set
    eval('parameters_default;');
end

%% replace parameters (if necessary)
if(nargin > 1 && exist('para', 'var') && isstruct(para))
    names = fieldnames(para);
    for i=1:length(names)
        if(isstruct(para.(names{i})))
            innerNames = fieldnames(para.(names{i}));
            for j=1:length(innerNames)
                eval(sprintf('%s.(innerNames{j}) = para.(names{i}).(innerNames{j});',names{i}));
                if(strcmp(innerNames{j},'trafodim') && ~any(strcmp(innerNames,'scrambledim')))
                    eval(sprintf('%s.scrambledim = para.(names{i}).(innerNames{j})(1,:);',names{i}));
                end
            end
            clear 'innerNames'
        else
            eval(sprintf('%s = para.%s;',names{i},names{i}));
        end
    end
    clear 'names'
end
clear 'para';
% save MC parameters
if(exist('savePath','var'))
    mc.savePath = savePath;
else
    mc.savePath = '';
end
% check for missing necessary parameters
loadedVars = {'measPara', 'iLC', 'iLCPositions', 'drecksMDH'};
for iI=1:length(loadedVars)
    if(eval(sprintf('~exist(''%s'',''var'')',loadedVars{iI})))
        eval(sprintf('%s = [];', loadedVars{iI}));
    end
end

%% initialize dispProgress
dispProgress('InitDispProgress',prop);

%% get data dimension and loop counters
[measPara, espresso, postproc] = extractMeasPara( iLC, measPara, espresso, postproc, flagOversampling, drecksMDH, iLCPositions );

%% check variable consistency and existence
allPara = { cstype,   measPara.dimension,   measPara.dim,   iNOUTER,   iNINNER,   p,   lambda,   epsilon,   measPara.oversampling,   lambdaCalib,   calibTyk,   reconTyk,   kernelSize,   calibSize,   opt_problem,   postproc,   trafo,   FFTwindow,   espresso,   measPara.LCall, lbtol,   tvmu,   cgtol,   cgmaxiter,   flagToolbox,   flagZeropadding,  measPara.aniso,  measPara.precision, reconDIM,   flags,   lambdaGroup,   lambdaNLTV,   lambdaTV;...
           'cstype', 'dimension',           'dim',         'iNOUTER', 'iNINNER', 'p', 'lambda', 'epsilon', 'oversampling',          'lambdaCalib', 'calibTyk', 'reconTyk', 'kernelSize', 'calibSize', 'opt_problem', 'postproc', 'trafo', 'FFTwindow', 'espresso', 'LCall',        'lbtol', 'tvmu', 'cgtol', 'cgmaxiter', 'flagToolbox', 'flagZeropadding', 'aniso',        'precision',        'reconDIM', 'flags', 'lambdaGroup', 'lambdaNLTV', 'lambdaTV'}; %#ok<*NODEF>

[consExist, msg, msgWarn, kernelSize, calibSize, FFTwindow, trafo, flagToolbox, measPara.precision, flags, reconDIM] = checkConsExist(allPara);
if(~consExist)
    error('CS_reconstruction:checkConsExist','%s', msg);
end
if(~isempty(msgWarn))
    warning('CS_reconstruction:checkConsExist\n%s', msgWarn);
end
clear 'allPara' 'msg' 'msgWarn' 'consExist';

%% initialize recon objects
if(strcmp(cstype,'FOCUSS'))
    obj = FOCUSS(iNOUTER, iNINNER, p, lambda, epsilon, measPara, lambdaCalib, lambdaTV, lambdaESPReSSo, lambdaMC, kernelSize, calibTyk, calibSize, flagZeropadding, FFTwindow, trafo, espresso, sScaling, mc, lSense);

elseif(strcmp(cstype,'SPIRiT_CG') || strcmp(cstype,'SPIRiT_POCS'))
    obj = SPIRiT_Wrapper(cstype(8:length(cstype)), iNINNER, measPara, lambda, kernelSize, calibTyk, reconTyk, trafo, espresso, FFTwindow, calibSize); %#ok<*COLND>

elseif(strcmp(cstype,'ESPIRiT_CG') || strcmp(cstype,'ESPIRiT_L1'))
    obj = ESPIRiT_Wrapper(cstype(9:length(cstype)), iNINNER, iNIterSplit, measPara, kernelSize, calibTyk, reconTyk, eigThresh_k, eigThresh_im, n_maps, splitWeight, trafo, flagToolbox, path, espresso, FFTwindow, calibSize);

elseif(strcmp(cstype,'BART'))
    obj = BART_Wrapper(lambda, lambdaCalib, lambdaTV, lambdaMC, measPara, kernelSize, n_maps, trafo, espresso, FFTwindow, currpath, mc, calibSize);

elseif(strcmp(cstype,'sparseMRI'))
    obj = sparseMRI(iNOUTER, iNINNER, measPara, lambda, lambdaTV, lambdaESPReSSo, p, trafo, espresso, FFTwindow, l1Smooth, lineSearchItnlim, lineSearchAlpha, lineSearchBeta, lineSearchT0, gradToll);

elseif(strcmp(cstype,'L1_Magic_TV') || strcmp(cstype,'L1_Magic_L1') || strcmp(cstype,'L1_Magic_TVDantzig') || strcmp(cstype,'L1_Magic_L1Dantzig'))
    obj = L1_Magic(cstype(10:length(cstype)), measPara, epsilon, lbtol, tvmu, cgtol, cgmaxiter, pdmaxiter, espresso);

elseif(strcmp(cstype,'Zero'))
    obj = Zero(measPara, espresso);

elseif(strcmp(cstype,'FCSA_WaTMRI') || strcmp(cstype,'FCSA_SLEP') || strcmp(cstype,'FCSA_proxA') || strcmp(cstype,'BFCSA_proxA') || strcmp(cstype,'ADMM_proxA') || strcmp(cstype,'SB_proxA') || ...
        strcmp(cstype,'A_PFISTA') || strcmp(cstype,'A_TDIHT') || strcmp(cstype,'A_ADMM_L1') || strcmp(cstype,'A_ADMM_L0') || strcmp(cstype,'A_ADMM_SCAD') || strcmp(cstype,'A_ADMM_MCP') || strcmp(cstype,'A_ADMM_ATAN') || strcmp(cstype,'A_ADMM_PSHRINK'))
    obj = Proximal(cstype, iNINNER, itrNLTV, iNOUTER, measPara, trafo, lambda, lambdaTV, lambdaGroup, lambdaNLTV, lambdaNLTV_h, mue, regularizerWeights, flags, reconDIM, espresso);

else
    error('CS_reconstruction(): Undefined reconstruction cstype');
end

%% optimize fft command
if(prop.flagOptim)
    optimizeFFT(measPara.dimension, measPara.dim, kernelSize, fft_planner_method, trafo.fftdim, lambdaCalib, flagZeropadding, measPara.precision);
end

%% extract kSpace data from measurement file
if(exist('iLC','var') && ~isempty(iLC) && ~exist('kSpace', 'var'))
    [obj.kSpace, evalMask] = extractKSpace_main(dData, iLC, iEvalInfoMask, measPara);
else
    obj.kSpace = kSpace;
    clear 'kSpace';
end

%% correct phase perturbations (echo alignment)
if(isfield(measPara,'sequenceName') && strcmp(measPara.sequenceName, 'CS_EPI'))
    obj.kSpace = alignEchos(obj.kSpace, evalMask, measPara, dData, iLC, iEvalInfoMask);
end

% clear all unnecessary variables
clear 'dData' 'iLC' 'dDataRep' 'iLCRep' 'dDataSli' 'iLCSli' 'drecksMDH' 'iLCPositions' 'evalMask' 'helper' 'espresso' 'trafo' 'flagOversampling' 'mc';
sCallingStack = dbstack;
if(strcmp(sCallingStack(end).name,'fRetroGateAndRecon') || strcmp(sCallingStack(end).name,'MotionCorrection')) % for 4D CS_Retro
    evalin('caller', 'clear ''kSpace'';'); % ''para''
end
clear 'sCallingStack';

%% start reconstruction
cs_time = tic;
imageCha = cell(size(obj.kSpace));

if(lAutoCheck && ( ((strcmp(measPara.dimension,'3D') || strcmp(measPara.dimension,'4D')) && nnz(obj.kSpace{1,1}) >= 0.95*prod(measPara.dim(1:4))) || ...
        (strcmp(measPara.dimension,'2D') && nnz(obj.kSpace{1,1}) >= 0.95*prod(measPara.dim([1:2,4]))) || ...
        (obj.espresso.state && ~isempty(regexp(filename,'\w*x1\D\w*','once'))) || (obj.espresso.state && isfield(measPara, 'CSAcceleration') && measPara.CSAcceleration == 1)))
    % check for acceleration factor == 1 ==> full acquisition
    if(strcmp(measPara.dimension,'2D'))
        fftdim = 1:2;
        pfdim = [3 1 2];
        subs = {':', ':', ':', ':'}; % y-x-cha
        subsid = 4;
        if(measPara.dim(4) > 1)
            pfdim = [4 1 2 3];
            subs = {':', ':', ':', ':'}; % y-x-t-cha
            subsid = 3;
        end
    elseif(strcmp(measPara.dimension,'3D'))
        fftdim = 1:3;
        pfdim = [4 1 2 3];
        subs = {':', ':', ':', ':', ':'}; % y-x-z-cha
        subsid = 5;
    elseif(strcmp(measPara.dimension,'4D'))
        fftdim = 1:3;
        pfdim = [5 1 2 3 4];
        subs = {':', ':', ':', ':', ':'}; % y-x-z-t-cha
        subsid = 4;
    end
    totalCount = measPara.LCall(2)*measPara.LCall(4)*measPara.LCall(1)*measPara.dim(4);
    dispProgress('Repetitions', 0, measPara.LCall(4));
    dispProgress('Averages', 0, totalCount);
    dispProgress('Time', 0, totalCount);
    for iRep=1:measPara.LCall(4)
        for iAvg=1:measPara.LCall(2)
            if(obj.espresso.state)
                for iSli=1:measPara.LCall(1)
                    kSpaceTmp = squeeze(cell2mat(shiftdim(obj.kSpace(iSli,:,iRep,iAvg),-2)));
                    imgTmp = zeros(size(kSpaceTmp));
                    for iTime=1:measPara.dim(4)
                        subs{subsid} = iTime;
                        imgTmp(subs{:}) = ipermute(espresso_recon(permute(kSpaceTmp(subs{:}),pfdim), obj.espresso.iter, false),pfdim);
                        dispProgress('Time', ((iRep-1)*measPara.LCall(2)*measPara.LCall(1)*measPara.dim(4) + (iAvg-1)*measPara.LCall(1)*measPara.dim(4) + (iSli-1)*measPara.dim(4) +iTime)/totalCount);
                    end
                    if(strcmp(measPara.dimension,'2D') && measPara.dim(4) == 1) % 2D
                        imageCha(iSli,:,iRep,iAvg) = shiftdim(squeeze(mat2cell(imgTmp, size(obj.kSpace{1,1},1),size(obj.kSpace{1,1},2),ones(1,measPara.dim(5)))),-1);
                    else % 2Dt, 3D, 4D
                        imageCha(iSli,:,iRep,iAvg) = shiftdim(squeeze(mat2cell(imgTmp, size(obj.kSpace{1,1},1),size(obj.kSpace{1,1},2),size(obj.kSpace{1,1},3),ones(1,measPara.dim(5)))),-1);
                    end
                    clear 'imgTmp' 'kSpaceTmp';
                end
            else
                imageCha(:,:,iRep,iAvg) = cellfun(@(x) ifftnshift(x,fftdim), obj.kSpace(:,:,iRep,iAvg), 'UniformOutput', false);
            end
            dispProgress('Averages', ((iRep-1)*measPara.LCall(2)*measPara.LCall(1)*measPara.dim(4) + iAvg*measPara.LCall(1)*measPara.dim(4))/totalCount);
        end
        dispProgress('Repetitions', iRep/measPara.LCall(4));
    end
    dispProgress('Time','Close');
else
    % Compressed Sensing reconstruction
    dispProgress('Repetitions', 0, measPara.LCall(4));
    dispProgress('Averages', 0, measPara.LCall(4)*measPara.LCall(2));
    for iRep=1:measPara.LCall(4)
        for iAvg=1:measPara.LCall(2)
            imageCha(:,:,iRep,iAvg) = obj.recon_main(iRep,iAvg); % obj.kSpace(:,:,iRep,iAvg)
            dispProgress('Averages', ((iRep-1)*measPara.LCall(2) + iAvg)/(measPara.LCall(4)*measPara.LCall(2)));
        end
        dispProgress('Repetitions', iRep/measPara.LCall(4));
    end
end
if(exist('cutOutMask','var')), imageCha = imageCha(cutOutMask); end;
dispProgress('Averages', 'Close');
dispProgress('Repetitions', 'Close');
obj.kSpace = []; % clear kSpace variable

cs_time = toc(cs_time);
dRecontime(1) = mod(cs_time,60);
dRecontime(2) = mod((cs_time - dRecontime(1))/60,60);
dRecontime(3) = mod(((cs_time - dRecontime(1))/60 - dRecontime(2))/60, 24);
if(cs_time >= 3600)
    fprintf('Execution time: %.2dh %.2dmin %.2ds\n', dRecontime(3), dRecontime(2), round(dRecontime(1)));
elseif(cs_time >= 60)
    fprintf('Execution time: %.2dmin %.2ds\n', dRecontime(2), round(dRecontime(1)));
else
    fprintf('Execution time: %.2ds\n', round(dRecontime(1)));
end
obj.dRecontime = dRecontime;

%% postprocessing
if(isempty(obj.kernelSize))
    obj.kernelSize = kernelSize;
end

image = cell(1,measPara.LCall(4));
for iRep=1:measPara.LCall(4)
    imageAvg = cell(1,measPara.LCall(2));
    for iAvg=1:measPara.LCall(2)
        % each element: 2D (y-x + one slice), 3D (y-x-t + one slice, y-x + several slices, y-x-z) or
        % 4D (y-x-t + several slices, y-x-z-t)
        imageAvg{1,iAvg} = postproc_main(obj, imageCha(:,:,iRep,iAvg), postproc);
    end
    image{iRep} = squeeze(mean(cell2mat(shiftdim(imageAvg,-4)),6));
end
image = squeeze(cell2mat(shiftdim(image,-3)));

% % average over repetitions
% image = squeeze(mean(cell2mat(shiftdim(image,-3)),5));

%% delete fftw_wisdom variable
if(prop.flagOptim)
    evalin('base', 'clear ''fftw_wisdom''');
end

%% save image
if(nargout > 2 || prop.flagSave)
    % convert obj->struct and reduce size
    objProp = properties(obj);
    for iObj = 1:numel(objProp)
        if(any(strcmp(objProp{iObj},{'kernel','kernelImg','kCalib','A','kSpace'})))
            continue;
        end
        eval(['objOut.',objProp{iObj},' = obj.',objProp{iObj},';']);
    end
end
if(prop.flagSave)
    if(exist('filename','var') && ~isempty(filename))
        savename = [filename,'_',cstype,'_',transformation];
    else
        savename = [cstype,'_',transformation];
    end
    dirContent = dir(path);
    occurence = 0;
    for i=1:length(dirContent)
        if(dirContent(i).isdir)
            continue;
        end
        tmp = regexp(dirContent(i).name,[savename,'[_]*\d*'],'once');
%         if(strcmpi(dirContent(i).name,savename))
        if(~isempty(tmp))
            occurence = occurence + 1;
        end
    end
    if(occurence > 0)
         savename = [savename,'_',num2str(occurence)];
    end
%     save([path,filesep,savename,'.mat'], 'image', 'imageCha', 'objOut');
    save([path,filesep,savename,'.mat'], 'image', 'objOut', '-v7.3');
end

%% plot resulting image
if(prop.flagPlot && exist('imagine','file'))
    if(iscell(image))
        imagePlot = cell2mat(shiftdim(image,-1));
    else
        imagePlot = image;
    end

    % set displayed name
    figname = [cstype,'_',transformation];
    tmp = functions(FFTwindow.type);
    if(~strcmp(tmp.function,'rectwin'))
        figname = [figname,'_',tmp.function];
    end
    if(obj.espresso.state)
        figname = [figname,'_espresso'];
    end
    figname = [figname,'_',postproc.type];

    % if imagine is already open, add new image in same window
    hfig = findobj('type','figure','name','IMAGINE 2.0 Belly Jeans');
    if(length(hfig) > 1)
        hfig = max(hfig);
    end

    if(~isempty(hfig) && ishandle(hfig))
        fh = get(hfig,'userdata');
        data = fh();

        oldImage = cell(1,size(data,2));
        oldTitle = oldImage;
        cmd = '';
        for i=1:size(data,2)
            oldImage{i} = data(i).dImg;
            oldTitle{i} = data(i).sName;
            cmd = [cmd, sprintf('oldImage{%d},''Name'',oldTitle{%d},',i,i)];
        end

        close(hfig);
        cmd = ['imagine(',cmd, 'imagePlot,''Name'',figname);'];
        eval(cmd);
        clear oldImage oldTitle cmd hfig;
    else
        imagine(imagePlot,'Name',figname);
    end
end

end

author: thomaskuestner


text:
classdef RiceD
%%RICED Functions to handle the Rice distribution.
%Implemented methods are: mean, var, PDF, CDF, rand
%
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

methods(Static)

function val=mean(s,sigma)
%%MEAN Obtain the mean of the Rice distribution for given noncentrality
%      and scale parameters.
%
%INPUTS: s The noncentrality parameter of the distribution.
%    sigma The scale parameter of the distribution.
%
%OUTPUTS: val The mean of the Rice distribution.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    param=-s^2/(2*sigma^2);
    val=sigma*sqrt(pi/2)*exp(param/2)*((1-param)*besseli(0,-param/2)-param*besseli(1,-param/2));
end

function val=var(s,sigma)
%%VAR Obtain the variance of the Rice distribution for given noncentrality
%     and scale parameters.
%
%INPUTS: s The noncentrality parameter of the distribution.
%    sigma The scale parameter of the distribution.
%
%OUTPUTS: val The variance of the Rice distribution.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    param=-s^2/(2*sigma^2);
    L=exp(param/2)*((1-param)*besseli(0,-param/2)-param*besseli(1,-param/2));
    val=2*sigma^2+s^2-pi*(sigma^2/2)*L^2;
end

function val=PDF(x,s,sigma)
%%PDF Evaluate the Rice probability distribution function at one or more
%     desired points.
%
%INPUTS: x The point(s) at which the Rice PDF is to be evaluated. Note that
%          x>=0.
%        s The noncentrality parameter of the distribution.
%    sigma The scale parameter of the distribution.
%
%OUTPUTS: val The value(s) of the Rice PDF with parameters s and sigma
%             evaluated at x.
%
%The PDF of the Rice distribution is given in Chapter 2.1.4 of [1].
%
%REFERENCES:
%[1] J. G. Proakis, Digital Communications. Ed. 4, Boston, MA:
%    McGraw Hill, 2001.
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    %We are evaluating
    %val=(x/sigma^2).*exp(-(x.^2+s^2)/(2*sigma^2)).*besseli(0,x*s/sigma^2);
    %However, if x*s/sigma^2 is large, then this function could return
    %NaNs. large doesn't have to be all that large. For example, s=30,
    %sigma=1, x=100. means that ratio is 3000. Thus, we use besseli with
    %the third argument to get a scaled value, then take the logarithm, do
    %everything else in the logarithmic domain and take the exponent to
    %undo the logarithm.

    besselArg=x*s/sigma^2;
    val=exp(log(besseli(0,besselArg,1))+besselArg-(x.^2+s^2)/(2*sigma^2)+log(x/sigma^2));
end

function val=CDF(x,s,sigma)
%%CDF Evaluate the cumulative distribution function of the Rice
%     distribution at desired points.
%
%INPUTS: x The point(s) at which the Rice CDF is to be evaluated. Note that
%          x>=0.
%        s The noncentrality parameter of the distribution.
%    sigma The scale parameter of the distribution.
%
%OUTPUTS: prob The value(s) of the CDF of the Rice distribution with
%              parameters k and theta evaluated at x.
%
%The CDF of the Rice distribution can be expressed in terms of Marcum's Q
%function as shown in Chapter 2.1.4 of [1].
%
%REFERENCES:
%[1] J. G. Proakis, Digital Communications. Ed. 4, Boston, MA:
%    McGraw Hill, 2001.
%
%February 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
    numPoints=length(x(:));
    val=zeros(size(x));

    for curPoint=1:numPoints
        val(curPoint)=1-MarcumQ(1,s/sigma,x(curPoint)/sigma);
    end
end

function vals=rand(N,s,sigma)
%%RAND Generate Rice-distributed random variables with the given
%      parameters.
%
%INPUTS: N If N is a scalar, then rand returns an NXN matrix of random
%          variables. If N=[M,N1] is a two-element row vector, then rand
%          returns an MXN1 matrix of random variables.
%        s The noncentrality parameter of the distribution.
%    sigma The scale parameter of the distribution.
%
%OUTPUTS: vals A matrix whose dimensions are determined by N of the
%              generated Rice random variables.
%
%This generates Rice distributed random variables by transforming normally
%distributed random variables using the identity given in Chapter 2.1.4 of
%[1].
%
%REFERENCES:
%[1] J. G. Proakis, Digital Communications. Ed. 4, Boston, MA:
%    McGraw Hill, 2001.
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(isscalar(N))
        dims=[N, N];
    else
        dims=N;
    end

    m=s/sqrt(2);

    X=sigma*randn(dims)+m;
    Y=sigma*randn(dims)+m;

    vals=sqrt(X.^2+Y.^2);
end

end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% Copyright 2018 Marc René Schädler
%
% This file is part of the mobile hearing aid prototype project
% The the mobile hearing aid prototype project 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.
%
% The mobile hearing aid prototype project 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 the mobile hearing aid prototype project. If not, see http://www.gnu.org/licenses/.

function writegaintable(filename, gt_data)
  fid = fopen(filename,'w');
  gt_min = -10;
  gt_step = 1;
  if fid > 0
    fprintf(fid,'gtdata = [');
    [num_rows, num_cols] = size(gt_data);
    for i=1:num_rows
      fprintf(fid,'[');
      for j=1:num_cols
        fprintf(fid,'%.2f',gt_data(i,j));
        if j < num_cols
          fprintf(fid,' ');
        else
          fprintf(fid,']');
        end
      end
      if i < num_rows
        fprintf(fid,';');
      else
        fprintf(fid,']');
      end
    end
    fprintf(fid,'\n');
    fprintf(fid,'gtmin = [');
    for i=1:num_rows
      fprintf(fid,'%.2f',gt_min);
      if i < num_rows
        fprintf(fid,' ');
      else
        fprintf(fid,']');
      end
    end
    fprintf(fid,'\n');
    fprintf(fid,'gtstep = [');
    for i=1:num_rows
      fprintf(fid,'%.2f',gt_step);
      if i < num_rows
        fprintf(fid,' ');
      else
        fprintf(fid,']');
      end
    end
    fprintf(fid,'\n');
    fclose(fid);
  end
end

author: m-r-s

repo: hearingaid-prototype


text:
%kkimportraw 'Low level import of raw (binary) data'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kimportraw.pane file
%
% Parameters:
% InputFile: i1 'Raw (binary) input', required: 'raw (binary) input'
% InputFile: i2 'Add-to object', optional: 'Object to which the new data should be added'
% Integer: wsize 'Width ', default: 1: 'Size of width dimension'
% Integer: hsize 'Height ', default: 1: 'Size of height dimension'
% Integer: dsize 'Depth ', default: 1: 'Size of depth dimension'
% Integer: tsize 'Time ', default: 1: 'Size of time dimension'
% Integer: esize 'Elements ', default: 1: 'Size of elements dimension'
% Integer: skip 'Offset or skip (bytes)', default: 0: 'skip this number of bytes before reading'
% OutputFile: o 'Output object', required: 'output object'
%
% Example: o = kkimportraw({i1, i2}, {'i1','';'i2','';'wsize',1;'hsize',1;'dsize',1;'tsize',1;'esize',1;'skip',0;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kimportraw - Low level import of raw (binary) data
%
%  DESCRIPTION
% .I kimportraw
% can be used to read raw binary data and import it to a data segment.  If
% the second input file is supplied, "Add-to Object", the new data
% will be appended onto that object.  Otherwise, an new object is created.
% The data is stored in the "Destination segment" specified.
%
% The source architecture, data type information, and data dimensions and
% structure  must be known to permit the data to be imported correctly.
%
% The -skip argument can be used to cause the specified number of bytes to be
% skipped before reading the data that will eventually end up in the output
% object. This is handy for skipping over a header or other information in
% a raw data file to get to the actual data.
%
% The index order of the data is assumed to be WHDTE.
%
% Access to the raw input is done using a special input mode (obtained
% by setting the KOBJ_RAW flag when the input file is opened) that forces
% the data to be interpreted in raw mode. This means that even if the
% data is actually in one of the supported formats, that fact will be
% ignored. This potentially lets you read the headers and other goodies
% and store them as data in the output object.
%
% If the destination of the raw data is specified to be the comment attribute
% of the output object, the raw data will be accessed as KUBYTE regardless of
% the setting of the -type argument. The number of characters given by the
% product of the sizes of the dimensions will be read in and appended to the
% comment attribute in the output object. Only those characters up to
% the first NULL in the data will actually end up in the comment
% since the comment is required to be a NULL terminated string.
%
%
%
%  EXAMPLES
% The following example shows how one might read a raw data file generated by a
% frame grabber. Here, we will assume that the data is a 512 col by 480 row
% unsigned byte color (RGB) image, with a 1024 byte header. We will assume that
% the color data is stored with the whole R plane, then the whole G plane, then
% the whole B plane.
% .DS
%   % kimportraw -o khoros.rgb -i1 framegrab.rgb -type 3 -wsize 512 -hsize 480 -esize 3 -segment 1 -skip 1024 -color RGB
% .DE
%
%  "SEE ALSO"
% kimportasc(1), kformats(1), kcomment(1), kcolorspace(1)
%
%  RESTRICTIONS
% The data is assumed to be stored in the raw file with a WHDTE index order.
% It would be nice to be able to directly handle mulitply interleaved data,
% but specifiying how that is to be done in 5D is a mess.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function SO3F = rotate_outer(SO3F,varargin)
% rotate ODF
%
% Syntax
%   SO3F = rotate_outer(SO3F,rot)
%
% Input
%  SO3F - @SO3Fun
%  rot  - @rotation
%
% Output
%  SO3F - @SO3FunHandle

SO3F = SO3FunHandle(@(rot) SO3F.eval(rot));

SO3F = rotate_outer(SO3F,varargin{:});

author: mtex-toolbox

repo: mtex

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


text:
function bessel_y0_int_values_test ( )

%*****************************************************************************80
%
%% BESSEL_YO_INT_VALUES_TEST demonstrates the use of BESSEL_YO_INT_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BESSEL_YO_INT_VALUES_TEST:\n' );
  fprintf ( 1, '  BESSEL_Y0_INT_VALUES stores values of \n' );
  fprintf ( 1, '  the integral of the Bessel function Y0.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X           FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = bessel_y0_int_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function status = test_exp_imp(modus)
%TEST_EXP_IMP tests the different expansions of sound fields in time-domain
%
%   Usage: status = test_exp_imp(modus)
%
%   Input parameters:
%       modus   - 0: numerical
%                 1: visual
%
%   Output parameters:
%       status  - true or false
%
%   TEST_EXP_IMP(modus) checks, if the circular basis expansions for plane waves
%   and point sources are working. The circular basis expansions are converted
%   to plane wave decompositions. Additionally, modal weighting functions are
%   tested. Optionally, sound field plots of the plane wave decompositions are
%   used for verification.

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

status = false;

%% ===== Checking of input  parameters ===================================
nargmin = 1;
nargmax = 1;
narginchk(nargmin,nargmax);

%% ===== Configuration ===================================================
% Parameters
conf = SFS_config;
conf.plot.loudspeakers = 0;  % do not plot loudspeakers
conf.plot.usedb = 1;  % use decibel for sound field plot
conf.modal_window_parameter = 2.0;  % parameter for kaiser window

c = conf.c;  % speed of sound

% for sound field plots
X = [-2, 2];
Y = [-2, 2];
Z = 0;

fhp = 1000;  % high-pass filter corner-frequency for pwd of point source

% plane waves with equi-angular distribution
Npw = 1024;
phi0 = (0:Npw-1).'*2*pi/Npw;
x0 = [cos(phi0) sin(phi0)];
x0(:,3) = 0;
x0(:,4:6) = x0(:,1:3);
x0(:,7) = 1;

% test scenarios
scenarios = { ...
  'pw', [ 0.0 -1.0 0.0],  5, 'rect'   , [0.0 0.0 0.0]
  'pw', [ 0.0 -1.0 0.0], 15, 'rect'   , [0.0 0.0 0.0]
  'pw', [ 0.0 -1.0 0.0], 15, 'kaiser' , [0.0 0.0 0.0]
  'pw', [ 0.0 -1.0 0.0],  5, 'rect'   , [0.5 1.0 0.0]
  'pw', [ 0.0 -1.0 0.0], 15, 'rect'   , [0.5 1.0 0.0]
  'pw', [ 0.0 -1.0 0.0], 15, 'kaiser' , [0.5 1.0 0.0]
  'ps', [ 0.0  2.5 0.0],  5, 'rect'   , [0.0 0.0 0.0]
  'ps', [ 0.0  2.5 0.0], 15, 'rect'   , [0.0 0.0 0.0]
  'ps', [ 0.0  2.5 0.0], 15, 'kaiser' , [0.0 0.0 0.0]
  'ps', [ 0.0  2.5 0.0],  5, 'rect'   , [0.5 1.0 0.0]
  'ps', [ 0.0  2.5 0.0], 15, 'rect'   , [0.5 1.0 0.0]
  'ps', [ 0.0  2.5 0.0], 15, 'kaiser' , [0.5 1.0 0.0]
  };

%% ===== Main ============================================================

for ii=1:size(scenarios)

    src = scenarios{ii,1};  % source type
    xs = scenarios{ii,2};  % source position / direction of plane wave
    Nce = scenarios{ii,3};  % modal order
    conf.modal_window = scenarios{ii,4};  % type of modal weighting function
    xq = scenarios{ii,5};  % expansion centre

    % circular expansion coefficients
    switch src
    case 'pw'
        [pm, delay_offset] = circexp_imp_pw(xs,Nce,xq,conf);
        t = xq*xs.'./c;
    case 'ps'
        [pm, delay_offset] = circexp_imp_ps(xs,Nce,xq,fhp,conf);
        t = norm(xq-xs)./c;
    end

    % modal weighting of coefficients
    wm = modal_weighting(Nce, conf);
    pm = bsxfun(@times, pm, wm);

    % conversion to plane wave decomposition
    ppwd = pwd_imp_circexp(pm, Npw);

    if modus
      % delay plane waves according to expansion centre
      delays = -xq*x0(:,1:3).'./c;
      delay_offset = delay_offset - min(delays);
      delays = delays - min(delays);
      ppwd = delayline(ppwd,delays,1,conf);

      % sound field plot
      sound_field_imp(X,Y,Z,x0,'pw',ppwd,t+delay_offset,conf);

      % title string
      str = 'Plane wave decompostion of modally bandlimited';
      switch src
      case 'pw'
          str = sprintf('%s plane wave', str);
      case 'ps'
          str = sprintf('%s point source', str);
      end
      str = sprintf(['%s ([%1.1f %1.1f %1.1f]):\n%s-window (M=%d), ' ...
        'center of modal expansion at [%1.1f %1.1f %1.1f]'], ...
        str, xs, conf.modal_window, Nce, xq);
      title(str);

    end
end

status = true;

author: sfstoolbox

repo: sfs-matlab


text:
% Euclidean distance between polyhedra
% Section 8.2.1, Boyd & Vandenberghe "Convex Optimization"
% Joelle Skaf - 10/09/05
%
% Given two polyhedra C = {x | A1*x <= b1} and D = {x | A2*x <= b2}, the
% distance between them is the optimal value of the problem:
%           minimize    || x - y ||_2
%               s.t.    A1*x <= b1
%                       A2*y <= b2

% Input data
randn('state',0);
rand('state',0);

n  = 5;
m1 = 2*n;
m2 = 3*n;
A1 = randn(m1,n);
A2 = randn(m2,n);
b1 = rand(m1,1);
b2 = rand(m2,1) + A2*randn(n,1);

% Solution via CVX
cvx_begin
    variables x(n) y(n)
    minimize (norm(x - y))
    A1*x <= b1;
    A2*y <= b2;
cvx_end

% Displaying results
disp('------------------------------------------------------------------');
disp('The distance between the 2 polyhedra C and D is: ' );
disp(['dist(C,D) = ' num2str(cvx_optval)]);

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [J grad] = nnCostFunction(nn_params, ...
                                   input_layer_size, ...
                                   hidden_layer_size, ...
                                   num_labels, ...
                                   X, y, lambda)
%NNCOSTFUNCTION Implements the neural network cost function for a two layer
%neural network which performs classification
%   [J grad] = NNCOSTFUNCTON(nn_params, hidden_layer_size, num_labels, ...
%   X, y, lambda) computes the cost and gradient of the neural network. The
%   parameters for the neural network are "unrolled" into the vector
%   nn_params and need to be converted back into the weight matrices.
%
%   The returned parameter grad should be a "unrolled" vector of the
%   partial derivatives of the neural network.
%

% Reshape nn_params back into the parameters Theta1 and Theta2, the weight matrices
% for our 2 layer neural network
Theta1 = reshape(nn_params(1:hidden_layer_size * (input_layer_size + 1)), ...
                 hidden_layer_size, (input_layer_size + 1));

Theta2 = reshape(nn_params((1 + (hidden_layer_size * (input_layer_size + 1))):end), ...
                 num_labels, (hidden_layer_size + 1));

% Setup some useful variables
m = size(X, 1);

% You need to return the following variables correctly
J = 0;
Theta1_grad = zeros(size(Theta1));
Theta2_grad = zeros(size(Theta2));

% ====================== YOUR CODE HERE ======================
% Instructions: You should complete the code by working through the
%               following parts.
%
% Part 1: Feedforward the neural network and return the cost in the
%         variable J. After implementing Part 1, you can verify that your
%         cost function computation is correct by verifying the cost
%         computed in ex4.m

% Add ones to the X data matrix
X = [ones(m, 1) X];
% Convert y from (1-10) class into num_labels vector
yd = eye(num_labels);
y = yd(y,:);

%%% Map from Layer 1 to Layer 2
a1=X;
% Coverts to matrix of 5000 examples x 26 thetas
z2=X*Theta1';
% Sigmoid function converts to p between 0 to 1
a2=sigmoid(z2);

%%% Map from Layer 2 to Layer 3
% Add ones to the h1 data matrix
a2=[ones(m, 1) a2];
% Converts to matrix of 5000 exampls x num_labels
z3=a2*Theta2';
% Sigmoid function converts to p between 0 to 1
a3=sigmoid(z3);

% Compute cost
%logisf=(-y)'*log(a3)-(1-y)'*log(1-a3);
logisf=(-y).*log(a3)-(1-y).*log(1-a3); % Becos y is now a matrix, so use dot product, unlike above
%J=((1/m).*sum(sum(logisf)));	% This line is correct if there is no regularization
% Try with ...
% J=((1/m).*sum((logisf)));
% That will give J in 10 columns (it has summed m samples), so need to sum again

%% Regularized cost
Theta1s=Theta1(:,2:end);
Theta2s=Theta2(:,2:end);
J=((1/m).*sum(sum(logisf)))+(lambda/(2*m)).*(sum(sum(Theta1s.^2))+sum(sum(Theta2s.^2)));

% Part 2: Implement the backpropagation algorithm to compute the gradients
%         Theta1_grad and Theta2_grad. You should return the partial derivatives of
%         the cost function with respect to Theta1 and Theta2 in Theta1_grad and
%         Theta2_grad, respectively. After implementing Part 2, you can check
%         that your implementation is correct by running checkNNGradients
%
%         Note: The vector y passed into the function is a vector of labels
%               containing values from 1..K. You need to map this vector into a
%               binary vector of 1's and 0's to be used with the neural network
%               cost function.
%
%         Hint: We recommend implementing backpropagation using a for-loop
%               over the training examples if you are implementing it for the
%               first time.
%
% Set all the D to zeros
tridelta_1=0;
tridelta_2=0;

% Compute delta, tridelta and big D
	delta_3=a3-y;
    z2=[ones(m,1) z2];
	delta_2=delta_3*Theta2.*sigmoidGradient(z2);
    delta_2=delta_2(:,2:end);
	tridelta_1=tridelta_1+delta_2'*a1; % Same size as Theta1_grad (25x401)
    tridelta_2=tridelta_2+delta_3'*a2; % Same size as Theta2_grad (10x26)
	Theta1_grad=(1/m).*tridelta_1;
    Theta2_grad=(1/m).*tridelta_2;
    %Theta1_grad=0;
	%Theta2_grad=0;
%end

% Part 3: Implement regularization with the cost function and gradients.
%
%         Hint: You can implement this around the code for
%               backpropagation. That is, you can compute the gradients for
%               the regularization separately and then add them to Theta1_grad
%               and Theta2_grad from Part 2.
%

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

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

% Unroll gradients
grad = [Theta1_grad(:) ; Theta2_grad(:)];

end

author: yhyap

repo: machine-learning-coursera


text:
function xferOrbits = getMultiFlybyXferOrbits(x, numRev, r1Sc, r2Sc, v1Sc, v2Sc, bodiesInfo, celBodyData)
%getMultiFlybyXferOrbits Summary of this function goes here
%   Detailed explanation goes here

%     [~, ~, secInDay, ~] = getSecondsInVariousTimeUnits();

    flybyBodies = bodiesInfo(2:end-1);
    numFB = length(flybyBodies);
    numBodies = length(bodiesInfo);
    numREVS = length(bodiesInfo) - 1;

    if(numFB == 0)
        tm = x(:,3);
        numRevInds = x(:,4);
    else
        tm = x(:,numBodies+1:numBodies+numREVS);
        numRevInds = x(:,numBodies+1+numREVS:end);
    end

    tm = round(tm);
    tm(tm==2) = -1;
    numTM = size(tm,2);

%     x = x(:,1:end-numTM);
%     daTimes = cumsum(x,2);
    x = x(:,1:end-numTM-numREVS);
    daTimes = cumsum(x,2);

    parentBodyInfo = bodiesInfo{1}.getParBodyInfo(celBodyData);
    gmu = parentBodyInfo.gm;

    xferOrbits = zeros(length(bodiesInfo)-1,10);
    for(i=1:length(bodiesInfo)-1) %#ok<*NO4LP>
        b1Info = bodiesInfo{i};
        b2Info = bodiesInfo{i+1};

        t1 = daTimes(i);
        t2 = daTimes(i+1);
        dt = t2 - t1;

        [r1, ~] = getStateAtTime(b1Info, t1, getParentGM(b1Info, celBodyData));
        [r2, ~] = getStateAtTime(b2Info, t2, getParentGM(b2Info, celBodyData));
        [v1, v2] = orbit.lambert(r1', r2', tm(i)*dt/86400, numRev(i), gmu);

        if(any(isnan(v1)) || any(isnan(v2))) %failsafe incase the solver does something dumb
            v1 = v1Sc(:,i)';
            v2 = v2Sc(:,i)';
        end

        [sma, ecc, inc, raan, arg, tru] = vect_getKeplerFromState(r1,v1',gmu);
        [~, ~, ~, ~, ~, tru2] = vect_getKeplerFromState(r2,v2',gmu);

        if(inc >= pi/2)
            [v1, v2] = orbit.lambert(r1', r2', -tm(i)*dt/86400, numRev(i), gmu);

            [sma, ecc, inc, raan, arg, tru] = vect_getKeplerFromState(r1,v1',gmu);
            [~, ~, ~, ~, ~, tru2] = vect_getKeplerFromState(r2,v2',gmu);
        end

        orbitData = [sma ecc inc raan arg tru tru2 t1 t2 gmu];
        xferOrbits(i,:) = orbitData;
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function g = p16_g ( n, x )

%*****************************************************************************80
%
%% P16_G evaluates the gradient for problem 16.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 May 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the values of the variables.
%
%    Output, real G(N), the gradient of the objective function.
%
  g = zeros ( n, 1 );

  f1 = 1.5   - x(1) * ( 1.0 - x(2)    );
  f2 = 2.25  - x(1) * ( 1.0 - x(2) * x(2) );
  f3 = 2.625 - x(1) * ( 1.0 - x(2) * x(2) * x(2) );

  df1dx1 = - ( 1.0 - x(2) );
  df1dx2 = x(1);
  df2dx1 = - ( 1.0 - x(2) * x(2) );
  df2dx2 = 2.0 * x(1) * x(2);
  df3dx1 = - ( 1.0 - x(2) * x(2) * x(2) );
  df3dx2 = 3.0 * x(1) * x(2) * x(2);

  g(1) = 2.0 * ( f1 * df1dx1 + f2 * df2dx1 + f3 * df3dx1 );
  g(2) = 2.0 * ( f1 * df1dx2 + f2 * df2dx2 + f3 * df3dx2 );

  return
end


text:
function type = r82poly2_type ( a, b, c, d, e, f )

%*****************************************************************************80
%
%% R82POLY2_TYPE analyzes a second order polynomial in two variables.
%
%  Discussion:
%
%    The polynomial has the form
%
%      A x^2 + B y^2 + C xy + Dx + Ey + F = 0
%
%    The possible types of the solution set are:
%
%     1: a hyperbola;
%        9x^2 -  4y^2       -36x - 24y -  36 = 0
%     2: a parabola;
%        4x^2 +  1y^2 - 4xy + 3x -  4y +   1 = 0;
%     3: an ellipse;
%        9x^2 + 16y^2       +36x - 32y -  92 = 0;
%     4: an imaginary ellipse (no real solutions);
%         x^2 +   y^2       - 6x - 10y + 115 = 0;
%     5: a pair of intersecting lines;
%                        xy + 3x -   y -   3 = 0
%     6: one point;
%         x^2 +  2y^2       - 2x + 16y +  33 = 0;
%     7: a pair of distinct parallel lines;
%                 y^2            -  6y +   8 = 0
%     8: a pair of imaginary parallel lines (no real solutions);
%                 y^2            -  6y +  10 = 0
%     9: a pair of coincident lines.
%                 y^2            -  2y +   1 = 0
%    10: a single line;
%                             2x -   y +   1 = 0;
%    11; all space;
%                                          0 = 0;
%    12; no solutions;
%                                          1 = 0;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Daniel Zwillinger, editor,
%    CRC Standard Mathematical Tables and Formulae,
%    CRC Press, 30th Edition, 1996, pages 282-284.
%
%  Parameters:
%
%    Input, real A, B, C, D, E, F, the coefficients.
%
%    Output, integer TYPE, indicates the type of the solution set.
%

%
%  Handle the degenerate case.
%
  if ( a == 0.0 && b == 0.0 && c == 0.0 )
    if ( d == 0.0 && e == 0.0 )
      if ( f == 0.0 )
        type = 11;
      else
        type = 12;
      end
    else
      type = 10;
    end
    return
  end

  delta = 8.0 * a * b * f ...
        + 2.0 * c * e * d ...
        - 2.0 * a * e * e ...
        - 2.0 * b * d * d ...
        - 2.0 * f * c * c;

  j = 4.0 * a * b - c * c;

  if ( delta ~= 0.0 )
    if ( j < 0.0 )
      type = 1;
    elseif ( j == 0.0 )
      type = 2;
    elseif ( 0.0 < j )
      if ( r8_sign ( delta ) ~= r8_sign ( a + b ) )
        type = 3;
      elseif ( r8_sign ( delta ) == r8_sign ( a + b ) )
        type = 4;
      end
    end
  elseif ( delta == 0.0 )
    if ( j < 0.0 )
      type = 5;
    elseif ( 0.0 < j )
      type = 6;
    elseif ( j == 0.0 )

      k = 4.0 * ( a + b ) * f - d * d - e * e;

      if ( k < 0.0 )
        type = 7;
      elseif ( 0.0 < k )
        type = 8;
      elseif ( k == 0.0 )
        type = 9;
      end

    end
  end

  return
end


text:
 function test = ir_read_op(arg1, arg2, arg3, arg4, arg5)
%function test = ir_read_op(dir, file, slice, frame, chat)
%function test = ir_read_op(file, slice, frame, chat)
%|
%| Use ASPIRE "op" command to read in a data file
%| from directory 'dir' (which is optional)
%| or optionally one slice,frame of a file
%| Jeff Fessler

slice = [];
frame = 0;
chat = 0;

if nargin < 1 || nargin > 5
	help read
	error args
elseif nargin == 1
	file = arg1;
else
	if ischar(arg2)
		file = [arg1 arg2];
		if nargin >= 3,	slice	= arg3;	end
		if nargin >= 4,	frame	= arg4;	end
		if nargin >= 5,	chat	= arg5;	end
	else
		file = arg1;
		if nargin >= 2,	slice	= arg2;	end
		if nargin >= 3,	frame	= arg3;	end
		if nargin >= 4,	chat	= arg4;	end
	end
end

if length(file) > 4
	ismat = strcmp(file((end-3):end), '.mat');
else
	ismat = 0;
end

op = sprintf('op -chat %d ', chat);

if exist(file) ~= 2
	fail(['File ' file ' not found!'])
end

tmp = '/tmp/tmp.mat';
if exist(tmp) == 2
	os_run(['/bin/rm -f ' tmp])
end

if ismat
	type = 'double';
else
%	type = '-';
	type = 'raw'; % ask op for raw file format
end

if ~isempty(slice)
	if length(slice) == 1
		choose = sprintf('%d %d %d %d 0 0', slice, slice, frame, frame);
	elseif length(slice) == 2
		choose = sprintf('%d %d %d %d 0 0', slice(1), slice(2), frame, frame);
	else
		error slice
	end
	comm = [op 'slice ' tmp ' ' file ' ' type ' ' choose];
else
	comm = [op 'convert ' tmp ' ' file ' ' type ' test'];
end

if chat
	disp(comm)
end

os_run(comm)

if ~exist(tmp, 'file')
	fail(['in read.m: Conversion error on: ' file])
end

load(tmp)
eval(['!/bin/rm -f ' tmp])

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
% MSH_REFINE: construct a refined mesh from a given one. The function only
%                refines the mesh, it does not refine any space.
%
%     msh_fine = msh_refine (msh, nsub);
%
% INPUTS:
%
%     msh:  an object of the msh_cartesian class (see msh_cartesian)
%     nsub: number of uniform subdivisions to apply on the elements, and for each direction
%
% OUTPUT:
%
%     msh_fine: an object of the class msh_cartesian (see msh_cartesian)
%
% Copyright (C) 2015 Rafael Vazquez
%
%    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/>.

function msh_fine = msh_refine (msh, nsub)

  rule = msh_gauss_nodes (msh.nqn_dir);
  [~, zeta] = kntrefine (msh.breaks, nsub-1, ones(1, msh.ndim), zeros(1, msh.ndim));
  [qn, qw] = msh_set_quad_nodes (zeta, rule);

  auxiliary_geometry.rdim = msh.rdim;
  auxiliary_geometry.map = msh.map;
  auxiliary_geometry.map_der = msh.map_der;
  if (isfield (struct (msh), 'map_der2'))
    auxiliary_geometry.map_der2 = msh.map_der2;
  end

  boundary = ~isempty (msh.boundary);
  if (msh.ndim > 1)
    bnd = [];
    for ii = 1:numel (msh.boundary)
      bnd(ii).rdim = msh.boundary(ii).rdim;
      bnd(ii).map = msh.boundary(ii).map;
      bnd(ii).map_der = msh.boundary(ii).map_der;
    end
    auxiliary_geometry.boundary = bnd;
  end
  msh_fine = msh_cartesian (zeta, qn, qw, auxiliary_geometry, 'boundary', boundary);

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function gB = flip(gB,ind)
% flip a boundary from phase1/phase2 to phase2/phase1
%
% Syntax
%   gB = flip(gB) % flip all boundaries
%   gB = flip(gB,ind) % flip the boundaries specified in ind
%
% Input
%  gB  - @grainBoundary
%  ind - indeces of the boundaries to flip
%
% Output
%  grains - @grainBoundary
%

if nargin == 2
  gB.ebsdId(ind,:) = fliplr(gB.ebsdId(ind,:));
  gB.grainId(ind,:) = fliplr(gB.grainId(ind,:));
  gB.phaseId(ind,:) = fliplr(gB.phaseId(ind,:));
  gB.misrotation(ind) = inv(gB.misrotation(ind));
else
  gB.eebsdId = fliplr(gB.ebsdId);
  gB.grainId = fliplr(gB.grainId);
  gB.phaseId = fliplr(gB.phaseId);
  gB.misrotation = inv(gB.misrotation);
end


text:
% Clear and close everything
clear all;
close all;

% Add paths needed
addpath('functions/matGeom/geom3d/')
addpath('functions/lips/')

% Load our data from file
path = '../input/floorplan_spencer_small.txt';
planes2d = load_2dplanedata(path);

% Convert to polygons with a height of 8ft
height = 8;
planes3d = planes2dtopolygons3d(planes2d, height);

% Plot the polyons (uses geom3d function)
figure;
for ii=1:size(planes3d,2)
	drawPolygon3d(planes3d{ii}(:,1),planes3d{ii}(:,2),planes3d{ii}(:,3),'b');
    hold on;
end
drawCoordinates3d([0,0,0],[0,0,0],10)
axis equal
xlabel('x-direction (ft)')
ylabel('y-direction (ft)')
zlabel('z-direction (ft)')

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
%% Copyright (C) 2015, 2016 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@symfun formula (@var{f})
%% Return a symbolic expression for this symfun.
%%
%% This returns the expression that defines the function given by
%% the symfun @var{f}.  Typically this is the RHS of the symfun:
%% @example
%% @group
%% syms x
%% f(x) = sin(x);
%% formula(f)
%%   @result{} ans = (sym) sin(x)
%% @end group
%% @end example
%%
%% The command @ref{@@symfun/argname} gives the independent variables
%% of the @@symfun.  Basically, @code{argname} for the independent
%% variables and @code{formula} for the dependent expression.
%%
%% If the symfun @var{f} is abstract, this returns @var{f} as a
%% sym:
%% @example
%% @group
%% syms f(x)
%% formula(f)   % but note it's a sym
%%   @result{} ans = (sym) f(x)
%% @end group
%% @end example
%%
%% @seealso{@@symfun/argnames}
%% @end defmethod

function g = formula(f)
  g = f.sym;
end

%!test
%! % simple
%! syms x
%! f(x) = sin(x);
%! g = formula(f);
%! assert (isequal (g, sin(x)));

%!test
%! % concrete: return is a sym, not a symfun
%! syms x
%! f(x) = sin(x);
%! g = formula(f);
%! assert (~isa(g, 'symfun'));

%!test
%! % abstract: return is a sym, not a symfun
%! syms f(x)
%! g = formula(f);
%! assert (~isa(g, 'symfun'));

repo: octsympy


text:
classdef MPDMP < PROBLEM
% <multi/many> <real>
% The multi-point distance minimization problem
% lower --- -100 --- Lower bound of decision variables
% upper ---  100 --- Upper bound of decision variables

%------------------------------- Reference --------------------------------
% M. Koppen and K. Yoshida, Substitute distance assignments in NSGA-II for
% handling many-objective optimization problems, Proceedings of the
% International Conference on Evolutionary Multi-Criterion Optimization,
% 2007, 727-741.
%------------------------------- 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)
        Points; % Vertexes
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            % Parameter setting
            [lower,upper] = obj.ParameterSet(-100,100);
            if isempty(obj.M); obj.M = 10; end
            obj.M        = max(obj.M,3);
            obj.D        = 2;
            obj.lower    = zeros(1,2) + lower;
            obj.upper    = zeros(1,2) + upper;
            obj.encoding = ones(1,obj.D);
            % Generate vertexes
            if mod(obj.M,2) == 0
                Angle = (2.*(1:obj.M)-3).*pi./obj.M;
            else
                Angle = (2.*(1:obj.M)-2).*pi./obj.M;
            end
            obj.Points = [sin(Angle)',cos(Angle)'];
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            PopObj = pdist2(PopDec,obj.Points);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            [X,Y] = ndgrid(linspace(-1,1,ceil(sqrt(N))));
            ND    = inpolygon(X(:),Y(:),obj.Points(:,1),obj.Points(:,2));
            R     = pdist2([X(ND),Y(ND)],obj.Points);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            if obj.M == 3
                [X,Y]    = ndgrid(linspace(-1,1,40));
                R        = pdist2([X(:),Y(:)],obj.Points);
                ND       = inpolygon(X(:),Y(:),obj.Points(:,1),obj.Points(:,2));
                R(~ND,:) = nan;
                R = {reshape(R(:,1),size(X)),reshape(R(:,2),size(X)),reshape(R(:,3),size(X))};
            else
                R = [];
            end
        end
        %% Display a population in the decision space
        function DrawDec(obj,Population)
            Draw(obj.Points([1:end,1],:),'-k','LineWidth',1.5,{'\it x\rm_1','\it x\rm_2',[]});
            Draw(obj.Points,'o','MarkerSize',6,'Marker','o','Markerfacecolor',[1 1 1],'Markeredgecolor',[.4 .4 .4]);
            Draw(Population.decs);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
addpath('functions/ht_pnp_function');
addpath('functions/at_netvlad_function');
addpath('functions/utils');
addpath('functions/wustl_function');
addpath('functions/relja_matlab');
addpath('functions/relja_matlab/matconvnet/');
addpath('functions/netvlad/');

addpath('functions/yael_matlab_linux64_v438');
addpath('functions/Inpaint_nans/Inpaint_nans');
run('functions/vlfeat/toolbox/vl_setup.m');
run('functions/matconvnet/matlab/vl_setupnn.m');

author: HajimeTaira


text:
function [u,dudx,dudp] = u_Fourier(x,P,t,in)
% input function for free-form (Fourier decomposition) deterministic DCM
% function [u,dudx,dudp] = u_Fourier(x,P,t,in)
% IN:
%   - x: [useless]
%   - P: nW*nuX1 vector of Fourier projection parameters, where nW is the
%   number of Fourier basis functions and nu is the dimensionality of the
%   output
%   - t: time at which the Fourier basis set is evaluated
%   - in: user-defined structure containing the set of Fourier frequencies
% OUT:
%   - u: the nuXnt matrix of outputs, where nt is the length of the
%   argument t.
%   - dudx: [useless]
%   - dudp: gradient wrt the parameters

% construct Fourier basis at time t
X = zeros(length(in.W),length(t));
for i = 1:length(t)
    X(:,i) = cos(in.W(:).*pi.*t(i)./in.T);
end

% project basis onto input space
P = reshape(P,length(in.W),[]);
u = (X'*P)';

% get gradients
dudx = [];
dudp = kron(eye(size(P,2)),X);

author: MBB-team

repo: VBA-toolbox


text:
% Copyright (C) 2013 Quan Wang <wangq10@rpi.edu>,
% Signal Analysis and Machine Perception Laboratory,
% Department of Electrical, Computer, and Systems Engineering,
% Rensselaer Polytechnic Institute, Troy, NY 12180, USA
%
% You are free to use this software for academic purposes if you cite our paper:
% Quan Wang and Kim L. Boyer.
% Feature Learning by Multidimensional Scaling and its Applications in Object Recognition.
% 2013 26th SIBGRAPI Conference on Graphics, Patterns and Images (Sibgrapi). IEEE, 2013.
%
% For commercial use, please contact the authors.

% Visualize the confusion matrix.

function draw_matrix(CM)
% CM: confusion matrix

figure;
imagesc(1-CM,[0 1]);
colormap gray;
N=size(CM,1);
for i=1:N
    for j=1:N
        if CM(i,j)>0.5
            the_color=[1 1 1];
        else
            the_color=[0 0 0];
        end
        text(j-0.2,i,[num2str(CM(i,j)*100,'%.2f') '%'],...
            'Color',the_color,'FontSize',15);
    end
end
title('confusion matrix');
xlabel('predicted label');
ylabel('true label');

author: Sable

repo: mcbench-benchmarks


text:
function pathG = build_initial_path()
global parameters
pathG=[];
j=1;
for i=1:size(parameters.trajectory,2)
    line = parameters.trajectory{i}.line;
    p0 = line(:,1);
    p1 = line(:,2);
    T = parameters.trajectory{i}.T0;

    partial_path = build_partial_path(p0, p1, T);
    %append
    for k=1:size(partial_path,2)
        pathG{j}.T=partial_path{k}.T;
        j=j+1;
    end
end

function partial_path = build_partial_path(p0, pf, T0)
global parameters
%this defines the line to follow in direction
deltaV = pf - p0;
deltaV = deltaV/norm(deltaV);
%and total length of the movement
total_length = norm(pf - p0); %m
n_steps = total_length/parameters.delta_mov;
n_steps = round(n_steps);
%linspace!!
%mov = linspace(0, total_length, parameters.N);
%line0 = [p0 pf];
%plot(line0(1,:),line0(2,:),'k')

partial_path=[];
Ti = T0;
for i=1:n_steps %length(mov)
    fprintf('Fast global plan %d out of %d\n', i, n_steps)
    %update to next point in trajectory
    %Ti(1:3,4) = p0 + mov(i)*deltaV;
    Ti(1:3,4) = p0 + i*parameters.delta_mov*deltaV;
    partial_path{i}.T = Ti;
end

repo: ARTE


text:
%
% Perform minmax norm
%
function [output, input_min, input_max, scale] = MinMaxNorm(input,min_val, max_val)
if nargin==1
    min_val = -1;
    max_val = 1;
end

input_min = min(input(:));
input_max = max(input(:));

input_range = input_max-input_min;
target_range = max_val - min_val;
scale = target_range ./ input_range;

output = bsxfun(@plus, input, -input_min);
output = bsxfun(@times, output, scale);
output = bsxfun(@plus, output, min_val);

end

author: singaxiong

repo: SignalGraph


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
% function [wOpt,his] = MLPIR(ML,varargin)
%
% Multi-Level Parametric Image Registration
%
% minimizes J^h(w) = D^h(T(Y(w)),R) + S^h(w) for h=coarse:fine
% see PIRobjFctn for a default objective and GaussNewton for a default optimizer
%
% Input:
%   ML        struct of coarse to fine representations of the data,
%           see getMultilevel.m
%  varargin optinonal parameters, see default parameter
%
% Output:
%  wOpt     optimal parameter for the parametric part
%  MLhis    iteration history
%
% for level=minLevel:maxLevel,
%   get data(level)
%   if level>minLevel, w0 = wOpt; end;
%   get wOpt running PIR  using w0 as starting guess
% end%For
%==============================================================================

function [wOpt,his] = MLPIR(ML,varargin)

if nargin == 0,
  help(mfilename);
  E6_HNSP_MLPIR_SSD_affine2D
  return;
end;

% setup default parameter for parametric pre-registration
PIRobj      = @PIRobjFctn;
PIRpara     = optPara('PIR-GN');
w0          = trafo('w0');  % starting guess
wStop       = w0;           % starting value (independent of level) used for stopping only

% configure regularizer
wRef        = w0;           % regularization: (w-wRef)'*M*(w-wRef)
M           = [];           %
beta        = 0;            % regularization: H -> H + beta*I
getGrid     = @getCellCenteredGrid;

% setup additional default parameter
pause       = 0;            % flag for pauses
plotIter    = 0;            % flag for output of iteration history each level
plotMLiter  = 0;            % flag for output of summarized iteration history
dimstr      = @(m) sprintf('m = [%s]',sprintf(' %d',m));

[ML,minLevel,maxLevel] = getMultilevel(ML);

for k=1:2:length(varargin),   % overwrites default parameter
  eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

% intialize
omega = ML{end}.omega;
wOpt  = w0;
his   = [];

FAIRmessage(sprintf('%s: MultiLevel Image Registration',mfilename),'#')
fprintf('>> %-20s : %s\n','omega',dimstr(omega));
fprintf('>> %-20s : %s\n','IMAGE MODEL',imgModel);
fprintf('>> %-20s : %s\n','DISTANCE',distance);
fprintf('>> %-20s : %s\n','TRAFO',trafo);
fprintf('>> %-20s :\n','PARARMETRIC-REGISTRATION');
fprintf('>> %-20s : objective=<%s>, method = <%s>\n',...
  'PIR',func2str(PIRobj),func2str(PIRpara.scheme));
fprintf('>> %-20s : minLevel=%d, maxLevel=%d\n','LEVEL',minLevel,maxLevel);
FAIRmessage('#')

FAIRmessage(sprintf('%s, minLevel=%d:maxLevel=%d',...
  'MultiLevel Parametric Image Registration',minLevel,maxLevel));

tic;
% -- for loop over all levels ---------------------------------------------
for level=minLevel:maxLevel;

  FAIRmessage(sprintf('%s: level %d from %d to %d, %s',...
    mfilename,level,minLevel,maxLevel,dimstr(ML{level}.m)));

  % get data for current level, compute interpolation coefficients
  m     = ML{level}.m;
  [T,R] = imgModel('coefficients',ML{level}.T,ML{level}.R,omega);

  % update transformation
  trafo('set','omega',omega,'m',m);

  % initialize plots
  PIRpara.Plots('reset','mode','PIR-multi level','fig',level);
  PIRpara.Plots('init',struct('Tc',T,'Rc',R,'omega',omega,'m',m));

  % ----- call PIR ------------------------------------
  xc   = getGrid(omega,m);
  Rc   = imgModel(R,omega,center(xc,m));
  fctn = @(wc) PIRobj(T,Rc,omega,m,beta,M,wRef,xc,wc);
  if level == minLevel,
    fctn([]);   % report status
  else
    w0 = wOpt;  % update starting guess
  end;
  PIRpara.yStop = wStop;
  f = fieldnames(PIRpara);
  v = struct2cell(PIRpara);
  temp = reshape({f{:};v{:}},1,[]);
  [wOpt,hisPIR] = PIRpara.scheme(fctn,w0,temp{:});
  % ----- end PIR --------------------------------------

  if plotIter,
    plotIterationHistory(hisPIR,'J',[1,2,5],'fig',20+level);
  end;

  % update iteration history
  if level == minLevel,
    his.str = hisPIR.str;
    his.his = hisPIR.his;
  else
    his.his = [his.his;hisPIR.his];
  end;
  doPause(pause)

end;%for level
% -- for loop over all levels ---------------------------------------------
his.time = toc;
if plotMLiter,
  plotMLIterationHistory(his,'fig',30);
end;
if isempty(wStop), wStop = w0; end
his.reduction = hisPIR.his(end,2) / (hisPIR.his(1,2)+(hisPIR.his(1,2)==0));
J = find(his.his(:,1)==-1);
his.iter(minLevel:maxLevel) = his.his([J(2:end)-1;size(his.his,1)],1)';

FAIRmessage([mfilename,' : done !']);

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

function doPause(p)
if strcmp(p,'on'),
  pause;
elseif p>0,
  pause(p);
end;
%==============================================================================


text:
function test_bug2464

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_selectdata

filename = dccnpath('/home/common/matlab/fieldtrip/data/test/bug2464.mat');
load(filename);

cfg = [];
cfg.foilim = [10 20];
sel = ft_selectdata(cfg, freq);

cfg = [];
cfg.frequency = [10 20];
sel2 = ft_selectdata(cfg, freq);

assert(isequal(rmfield(sel, 'cfg'), rmfield(sel2, 'cfg')));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [res, list] = modality(this, scalp, planar)
% Returns data modality
% FORMAT [res, list] = modality(this, scalp)
%
% scalp - 1 (default) only look at scalp modalities
%         0  look at all modalities
% planar - 1 distinguish between MEG planar and other MEG
%          0 (default) do not distinguish
% If more than one modality is found the function returns 'Multimodal'
% in res and a cell array of modalities in list.
% _______________________________________________________________________
% Copyright (C) 2008-2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: modality.m 6542 2015-09-09 11:48:34Z karl $

if nargin == 1
    scalp = 1;
end
if nargin < 3
    planar = 0;
end

list = {};

if ~isempty(indchantype(this, {'MEG', 'MEGPLANAR', 'MEGCOMB'}))
    if planar
        if ~isempty(indchantype(this, 'MEGPLANAR'))
            list = [list {'MEGPLANAR'}];
        end
        if ~isempty(indchantype(this, 'MEG'))
            list = [list {'MEG'}];
        end
        if ~isempty(indchantype(this, 'MEGCOMB'))
            list = [list {'MEGCOMB'}];
        end
    else
        list = [list {'MEG'}];
    end
end

if ~isempty(indchantype(this, 'EEG'))
    list = [list {'EEG'}];
end

if  ~isempty(indchantype(this, 'LFP')) && ~scalp
    list = [list {'LFP'}];
end

if  ~isempty(indchantype(this, 'ILAM')) && ~scalp
    list = [list {'ILAM'}];
end

switch numel(list)
    case 0
        res = 'Other';
    case 1
        res = list{1};
    otherwise
        res = 'Multimodal';
end

author: spm


text:
classdef AnimatorPointOfView < uint8
    % The list of animator point of view
    %
    % @author omar @date 2017-06-01
    %
    % Copyright (c) 2017, UMICH Biped Lab
    % All right reserved.
    %
    % Redistribution and use in source and binary forms, with or without
    % modification, are permitted only in compliance with the BSD 3-Clause
    % license, see
    % http://www.opensource.org/licenses/BSD-3-Clause

    enumeration
        North (0)
        South (1)
        East (2)
        West (3)
        Front (5)
        Back (6)
        Right (7)
        Left (8)
        TopSouthEast (9)
        TopFrontLeft (10)
        Free (256)
    end
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev