text:
function [fMc, fProbability, mResult] =  plot_SchusterMc(mCatalog)
% function [fMc, fProbability, mResult] = plot_SchusterMc(mCatalog)
% ----------------------------------------------------
% Determine the magnitude of completness using Schuster's Method
% objectively. Philosophy: Detect if walkout > 95%
% significance level two times, then set upper magnitude as Mc. If that
% happens more times, choose the smaller magnitude!
%
% Incoming variables:
% mCatalog : EQ catalog
%
% Outgoing variables:
% fMc          : Magnitude of completeness
% fProbability : Probability of exceeding the 95% level radius
% mResult      : Result matrix for all cases
%
% Author: J. Woessner, woessner@seismo.ifg.ethz.ch
% last update: 04.06.02

fMinMag = round(min(mCatalog(:,6))*10)/10;
fMaxMag = round(max(mCatalog(:,6))*10)/10;

prompt = {'Enter minimum magnitude:','Enter maximum magnitude:','Magnitude increment:'};
dlg_title = 'Search parameters';
num_lines= 1;
def     = {num2str(fMinMag),num2str(fMaxMag),'0.2'};
answer  = inputdlg(prompt,dlg_title,num_lines,def);
fMinMag = str2double(answer(1));
fMaxMag = str2double(answer(2));
fMagIncr = str2double(answer(3));

% Initialize
mResult = [];

fMag=fMinMag;
i=0;
figure_w_normalized_uicontrolunits('tag','schuster','Name','Schuster Walkout');
figure_w_normalized_uicontrolunits('tag','schuster2','Name','Schuster Walkout');
%n = floor((fMaxMag-fMag)/fMagIncr);
n = floor((fMaxMag-fMagIncr-fMag)/0.1)+2;
if mod(n,2) ~= 0
    fN=round(n/2);
else
    fN=floor(n/2);
end
while (fMag+fMagIncr) < (fMaxMag)
    % Calculate result matrix
    vSel = (mCatalog(:,6) >= fMag & mCatalog(:,6) < fMag+fMagIncr);
    mCat = mCatalog(vSel,:);
    [mWalkout, fR95, fProb, PHI, R] = calc_Schusterwalk(mCat);
    [vThetaWalkout,vRadWalkout] = cart2pol(mWalkout(:,1),mWalkout(:,2));
    fMaxRadius = max(abs(vRadWalkout(:,1)));
    mResult = [mResult; fMag+fMagIncr, fMaxRadius, max(abs(mWalkout(:,1))), max(abs(mWalkout(:,2))), fR95, fProb, R];

    % Subplot count
    i=i+1;
    figure_w_normalized_uicontrolunits(findobj('tag','schuster'));
    % Plot Schuster walkout
    subplot(fN,2,i)
    plot(mWalkout(:,1),mWalkout(:,2),'-','Color',[0.5 0.5 0.5]);
    hold on;
    polar([0:360]*pi/180,ones(1,361)*fR95,'--k');
    [x, y] =pol2cart(PHI,R);
    %plot([0,x],[0 y],'g');
    plot ([0 0],[-1/5*fR95,1/5*fR95],'k');
    plot ([-1/5*fR95,1/5*fR95],[0 0],'k');
%     hold off
%     text (0,2/5*fR95,'0:00','HorizontalAlignment','center');
%     text (2/5*fR95,0,'6:00','HorizontalAlignment','left');
%     text (0,-2/5*fR95,'12:00','HorizontalAlignment','center');
%     text (-2/5*fR95,0,'18:00','HorizontalAlignment','right');
%     text (0,fR95 ,'95KI','VerticalAlignment','bottom','HorizontalAlignment','center')
    axis square;
    sTitlestr = ['M = ' num2str(fMag) ' - ' num2str(fMag+fMagIncr)];
    title(sTitlestr);
    % Radius comparison for the plot
    sTextstr = [num2str(length(mCat(:,8))) ' / ' num2str(fMaxRadius/fR95)];
    text (-fR95,-fR95,sTextstr);

    % Plot hourly histogram
    figure_w_normalized_uicontrolunits(findobj('tag','schuster2'));
    subplot(n,1,i);
    vTime = (mCat(:,8)*60+mCat(:,9))/60; % Calculate decimal hour
    histogram(vTime,-0.5:1:24.5)
    fMag= fMag+0.1;
end % END of WHILE

% Select walkout > 95% level
try
    vSel = (mResult(:,2) >= mResult(:,5));
    mResOut = mResult(vSel,:);
    mDiffResOut=diff(mResOut);
    mDiffResOut(:,1) = round(mDiffResOut(:,1)*10)/10;
    [vIndice]=find(mDiffResOut(:,1) == 0.1);
    if isempty(vIndice)
        fMc = nan;
        fProbability = nan;
    else
        fMc = mResult(max(vIndice)+2,1); % Maximum value for Mc
        fProbability = mResult(max(vIndice)+2,6);
        for nCnt=length(vIndice):-1:2
            fdI = vIndice(nCnt)-vIndice(nCnt-1);
            if fdI > 1
                fMc = mResult(vIndice(nCnt-1)+2,1); % Find Mc
                fProbability = mResult(max(vIndice)+2,6);
            end
        end
    end
catch
    fMc = nan;
    fProbability = nan;
end

author: CelsoReyes


text:
filename='Gripping_quad_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'nonadjoint_compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
constraint_case = 'INEQUALITY';
optimizer = 'SLERP'; incrementFactor = 1;
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.5;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
classdef prtUiControl < hgsetget

    properties

        handle
    end

    properties (Dependent)
        backgroundColor
        beingDeleted
        busyAction
        buttonDownFcn
        callback
        cData
        children
        createFcn
        deleteFcn
        enable
        extent
        fontAngle
        fontName
        fontSize
        fontUnits
        fontWeight
        foregroundColor
        handleVisibility
        hitTest
        horizontalAlignment
        interruptible
        keyPressFcn
        listboxTop
        max
        min
        parent
        position
        selected
        selectionHighlight
        sliderStep
        string
        style
        tag
        tooltipString
        type
        uiContextMenu
        units
        userData
        value
        visible
    end
    methods
        function self = prtUiControl(varargin)
            self.handle = uicontrol(varargin{:});
        end

        function set.backgroundColor(self, val)
            set(self.handle,'backgroundColor',val);
        end
        function val = get.backgroundColor(self)
            val = get(self.handle,'backgroundColor');
        end
        function set.beingDeleted(self, val)
            set(self.handle,'beingDeleted',val);
        end
        function val = get.beingDeleted(self)
            val = get(self.handle,'beingDeleted');
        end
        function set.busyAction(self, val)
            set(self.handle,'busyAction',val);
        end
        function val = get.busyAction(self)
            val = get(self.handle,'busyAction');
        end
        function set.buttonDownFcn(self, val)
            set(self.handle,'buttonDownFcn',val);
        end
        function val = get.buttonDownFcn(self)
            val = get(self.handle,'buttonDownFcn');
        end
        function set.callback(self, val)
            set(self.handle,'callback',val);
        end
        function val = get.callback(self)
            val = get(self.handle,'callback');
        end
        function set.cData(self, val)
            set(self.handle,'cData',val);
        end
        function val = get.cData(self)
            val = get(self.handle,'cData');
        end
        function set.children(self, val)
            set(self.handle,'children',val);
        end
        function val = get.children(self)
            val = get(self.handle,'children');
        end
        function set.createFcn(self, val)
            set(self.handle,'createFcn',val);
        end
        function val = get.createFcn(self)
            val = get(self.handle,'createFcn');
        end
        function set.deleteFcn(self, val)
            set(self.handle,'deleteFcn',val);
        end
        function val = get.deleteFcn(self)
            val = get(self.handle,'deleteFcn');
        end
        function set.enable(self, val)
            set(self.handle,'enable',val);
        end
        function val = get.enable(self)
            val = get(self.handle,'enable');
        end
        function set.extent(self, val)
            set(self.handle,'extent',val);
        end
        function val = get.extent(self)
            val = get(self.handle,'extent');
        end
        function set.fontAngle(self, val)
            set(self.handle,'fontAngle',val);
        end
        function val = get.fontAngle(self)
            val = get(self.handle,'fontAngle');
        end
        function set.fontName(self, val)
            set(self.handle,'fontName',val);
        end
        function val = get.fontName(self)
            val = get(self.handle,'fontName');
        end
        function set.fontSize(self, val)
            set(self.handle,'fontSize',val);
        end
        function val = get.fontSize(self)
            val = get(self.handle,'fontSize');
        end
        function set.fontUnits(self, val)
            set(self.handle,'fontUnits',val);
        end
        function val = get.fontUnits(self)
            val = get(self.handle,'fontUnits');
        end
        function set.fontWeight(self, val)
            set(self.handle,'fontWeight',val);
        end
        function val = get.fontWeight(self)
            val = get(self.handle,'fontWeight');
        end
        function set.foregroundColor(self, val)
            set(self.handle,'foregroundColor',val);
        end
        function val = get.foregroundColor(self)
            val = get(self.handle,'foregroundColor');
        end
        function set.handleVisibility(self, val)
            set(self.handle,'handleVisibility',val);
        end
        function val = get.handleVisibility(self)
            val = get(self.handle,'handleVisibility');
        end
        function set.hitTest(self, val)
            set(self.handle,'hitTest',val);
        end
        function val = get.hitTest(self)
            val = get(self.handle,'hitTest');
        end
        function set.horizontalAlignment(self, val)
            set(self.handle,'horizontalAlignment',val);
        end
        function val = get.horizontalAlignment(self)
            val = get(self.handle,'horizontalAlignment');
        end
        function set.interruptible(self, val)
            set(self.handle,'interruptible',val);
        end
        function val = get.interruptible(self)
            val = get(self.handle,'interruptible');
        end
        function set.keyPressFcn(self, val)
            set(self.handle,'keyPressFcn',val);
        end
        function val = get.keyPressFcn(self)
            val = get(self.handle,'keyPressFcn');
        end
        function set.listboxTop(self, val)
            set(self.handle,'listboxTop',val);
        end
        function val = get.listboxTop(self)
            val = get(self.handle,'listboxTop');
        end
        function set.max(self, val)
            set(self.handle,'max',val);
        end
        function val = get.max(self)
            val = get(self.handle,'max');
        end
        function set.min(self, val)
            set(self.handle,'min',val);
        end
        function val = get.min(self)
            val = get(self.handle,'min');
        end
        function set.parent(self, val)
            set(self.handle,'parent',val);
        end
        function val = get.parent(self)
            val = get(self.handle,'parent');
        end
        function set.position(self, val)
            set(self.handle,'position',val);
        end
        function val = get.position(self)
            val = get(self.handle,'position');
        end
        function set.selected(self, val)
            set(self.handle,'selected',val);
        end
        function val = get.selected(self)
            val = get(self.handle,'selected');
        end
        function set.selectionHighlight(self, val)
            set(self.handle,'selectionHighlight',val);
        end
        function val = get.selectionHighlight(self)
            val = get(self.handle,'selectionHighlight');
        end
        function set.sliderStep(self, val)
            set(self.handle,'sliderStep',val);
        end
        function val = get.sliderStep(self)
            val = get(self.handle,'sliderStep');
        end
        function set.string(self, val)
            set(self.handle,'string',val);
        end
        function val = get.string(self)
            val = get(self.handle,'string');
        end
        function set.style(self, val)
            set(self.handle,'style',val);
        end
        function val = get.style(self)
            val = get(self.handle,'style');
        end
        function set.tag(self, val)
            set(self.handle,'tag',val);
        end
        function val = get.tag(self)
            val = get(self.handle,'tag');
        end
        function set.tooltipString(self, val)
            set(self.handle,'tooltipString',val);
        end
        function val = get.tooltipString(self)
            val = get(self.handle,'tooltipString');
        end
        function set.type(self, val)
            set(self.handle,'type',val);
        end
        function val = get.type(self)
            val = get(self.handle,'type');
        end
        function set.uiContextMenu(self, val)
            set(self.handle,'uiContextMenu',val);
        end
        function val = get.uiContextMenu(self)
            val = get(self.handle,'uiContextMenu');
        end
        function set.units(self, val)
            set(self.handle,'units',val);
        end
        function val = get.units(self)
            val = get(self.handle,'units');
        end
        function set.userData(self, val)
            set(self.handle,'userData',val);
        end
        function val = get.userData(self)
            val = get(self.handle,'userData');
        end
        function set.value(self, val)
            set(self.handle,'value',val);
        end
        function val = get.value(self)
            val = get(self.handle,'value');
        end
        function set.visible(self, val)
            set(self.handle,'visible',val);
        end
        function val = get.visible(self)
            val = get(self.handle,'visible');
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function dataOut=plotMultipleProjectedAmps_SingleCondition(view,scanToPlot,projectionPhase)
%
% plotMultipleProjectedAmps_SingleCondition(view,projectionPhase)
%
% Bar plot of the amplitudes for each scan, averaging across
% all pixels (in all slices) in a selection of ROIs. All y-axes are made the same. The bar heights
% and a coarse SEM can be obtained from get(gca,'UserData').
% Amplitudes are projected against a single phase (supplied).

% gmb  5/25/98
% bw   2/19/99  Added seY field to the UserData field.
%	    seY is an estimate of the variability in the
%      amplitudes.  It is the SEM of the in the complex
%      (amp*exp(-i*ph)) representation.  The values are
%      computed in vectorMean.m
% fwc   11/07/02 plots data relative to current view
%       added plotting of multiple ROIs
%       ROI selection copied from plotMultipleTSeries.m
% arw   042505 Correctly extracts projections phase from current scan if
%   none is sprcified.
%   Projection phase is in radians.
%

mrGlobals;

refScan = getCurScan(view);
if (~exist('projectionPhase','var'))
    computeProjPhase=1;
else
    if (isnan(projectionPhase))
        % User input requested
        projectionPhase=input('Enter projection phase in radians:');
    end

    computeProjPhase=0;
end

% Select ROIs
nROIs=size(view.ROIs,2);
roiList=cell(1,nROIs);
for r=1:nROIs
    roiList{r}=view.ROIs(r).name;
end
selectedROIs = find(buttondlg('ROIs to Plot',roiList));
nROIs=length(selectedROIs);
if (nROIs==0)
    error('No ROIs selected');
end

% Plot it
selectGraphWin
clf
fontSize = 8;
headerStr = ['Mean Amplitudes'];
set(gcf,'Name',headerStr);

minylim=0;
maxylim=0;
nrows=0;
ncols=0;

nscans = numScans(view);
ROIamps=zeros(nROIs,1);
ROIseZ=zeros(nROIs,1);
ROImeanPhs=zeros(nROIs,1);

for r=1:nROIs

    n=selectedROIs(r);
    view = selectROI(view,n); % is there another way? Well yes - we could be opening up a hidden window and doing all this invisibly.

    [meanAmps,meanPhs,seZ] = vectorMeans(view);
    meanAmps=meanAmps(scanToPlot);
    meanPhs=meanPhs(scanToPlot);
    seZ=seZ(scanToPlot);

    if (computeProjPhase)
        projectionPhase=meanPhs(refScan);
    end

    % Compute the amplitude projected onto the reference phase
    meanAmps = meanAmps.*cos(meanPhs-projectionPhase);

    ROIamps(r)=meanAmps;
    ROIseZ(r)=seZ;

    ROImeanPhs(r)=meanPhs;
    %xstr{r}=[view.ROIs(selectedROIs(r)).name];
    xstr{r}=int2str(r);
    roiName{r}=view.ROIs(selectedROIs(r)).name;
    fprintf(['\n#%d :',roiName{r}],r);
end

dataOut.ROIamps=ROIamps;
dataOut.ROIseZ=ROIseZ;
dataOut.ROIname=roiName;

% Now do the plotting

    nrows=1;
    ncols=1;

    subplot(nrows,ncols,1);

    %plot the bar graph
    size(ROIamps)
    size(ROIseZ)
    size(xstr)

    h=mybar(ROIamps,ROIseZ,xstr,'');

    xlabel('ROI','FontSize',fontSize);
    ylabel('Mean Amplitude','FontSize',fontSize);
    set(gca,'FontSize',ceil(fontSize*1.2));
    conditionName=dataTYPES(view.curDataType).scanParams(scanToPlot).annotation;
    fprintf('\nCondition #%d :',scanToPlot);

    title(['Condition #: ' int2str(scanToPlot)]);

end

%Save the data in gca('UserData')
data.y =ROIamps(:);
data.refScan = refScan;
data.seY = ROIseZ(:); % this should probably be adapted

set(gca,'UserData',data);

% give all plots same y-axis

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function dt_stability_limit = checkStability(kgrid, medium)
% CHECKSTABILITY   Compute maximum stable timestep for k-space models
%
% DESCRIPTION:
%       checkStability calculates the maximum time step for which the
%       k-space propagation models kspaceFirstOrder1D, kspaceFirstOrder2D
%       and kspaceFirstOrder3D are stable. These models are unconditionally
%       stable when the reference sound speed is equal to or greater than
%       the maximum sound speed in the medium and there is no absorption.
%       However, when the reference sound speed is less than the maximum
%       sound speed the model is only stable for sufficiently small time
%       steps. The criterion is more stringent (the time step is smaller)
%       in the absorbing case.
%
%       The time steps given are accurate when the medium properties are
%       homogeneous. For a heterogeneous media they give a useful, but not
%       exact, estimate.
%
%       The timesteps given are accurate when the medium properties are
%       homogeneous. For a heterogeneous media they give a useful, but not
%       exact, estimate.
%
% USAGE:
%       dt_stability_limit = checkStability(kgrid, medium)
%
% INPUTS:
%
%       kgrid              - structure returned by makeGrid holding grid parameters
%       medium             - structure holding the medium properties
%
% OUTPUTS:
%
%       dt_stability_limit - maximum timestep for stability. (Set to Inf
%                            when the model is unconditionally stable.)
%
% ABOUT:
%       author             - Ben Cox
%       date               - 12th August 2014
%       last update        - 25th August 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox
%
% See also kspaceFirstOrder1D, kspaceFirstOrder2D, kspaceFirstOrder3D,
% makeGrid, makeTime

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

% Find the maximum wavenumber
kmax = max(kgrid.k(:));

% calculate the reference sound speed for the fluid code, using the
% maximum by default which ensures the model is unconditionally stable
if isfield(medium, 'sound_speed_ref')
    if isnumeric(medium.sound_speed_ref)
        c_ref = medium.sound_speed_ref;
    elseif strcmp(medium.sound_speed_ref, 'min')
        c_ref = min(medium.sound_speed(:));
    elseif strcmp(medium.sound_speed_ref, 'mean')
        c_ref = mean(medium.sound_speed(:));
    else strcmp(medium.sound_speed_ref, 'max')
        c_ref = max(medium.sound_speed(:));
    end
else
    c_ref = max(medium.sound_speed(:));
end

% calculate the timesteps required for stability
if ~isfield(medium, 'alpha_coeff') || all(medium.alpha_coeff(:) == 0)

    % =====================================================================
    % NON-ABSORBING CASE
    % =====================================================================

    if c_ref >= max(medium.sound_speed(:))

        % set the timestep to Inf when the model is unconditionally stable.
        dt_stability_limit = Inf;

    else

        % set the timestep required for stability when c_ref~=max(medium.sound_speed(:))
        dt_stability_limit = 2/(c_ref * kmax) * asin(c_ref/max(medium.sound_speed(:)));

    end

else

    % =====================================================================
    % ABSORBING CASE
    % =====================================================================

    % convert the absorption coefficient to nepers.(rad/s)^-y.m^-1
    medium.alpha_coeff = db2neper(medium.alpha_coeff, medium.alpha_power);

    % calculate the absorption constant
    if ~(isfield(medium, 'alpha_mode') && strcmp(medium.alpha_mode, 'no_absorption'))
        absorb_tau = -2*medium.alpha_coeff.*medium.sound_speed.^(medium.alpha_power - 1);
    else
        absorb_tau = 0;
    end

    % calculate the dispersion constant
    if ~(isfield(medium, 'alpha_mode') && strcmp(medium.alpha_mode, 'no_dispersion'))
        absorb_eta = 2*medium.alpha_coeff.*medium.sound_speed.^(medium.alpha_power)*tan(pi*medium.alpha_power/2);
    else
        absorb_eta = 0;
    end

    % Estimate the timestep required for stability in the absorbing case by
    % assuming the k-space correction factor, kappa = 1;
    % (Note that absorb_tau and absorb_eta are negative quantities.)
    temp1 = max(medium.sound_speed(:)) * min(absorb_tau(:)) * kmax^(medium.alpha_power-1);
    temp2 = 1 - min(absorb_eta(:)) * kmax^(medium.alpha_power-1);
    dt_estimate = (temp1 + sqrt(temp1^2 + 4*temp2))/(temp2 * kmax * max(medium.sound_speed(:)));

    % Use a fixed point iteration to find the correct timestep, assuming
    % now that kappa = kappa(dt), using the previous estimate as a starting
    % point

    % First define the function to iterate
    kappa = @(dt) sinc(c_ref*kmax*dt/2);
    temp3 = @(dt) max(medium.sound_speed(:))*min(absorb_tau(:))*kappa(dt)*kmax^(medium.alpha_power-1);
    func_to_solve = @(dt) (temp3(dt) + sqrt((temp3(dt))^2 + 4*temp2 ))/(temp2*kmax*kappa(dt)*max(medium.sound_speed(:)));

    % run the fixed point iteration
    dt_stability_limit = dt_estimate;
    dt_old = 0;
    accuracy = 1e-12;
    while abs(dt_stability_limit-dt_old)>accuracy
        dt_old = dt_stability_limit;
        dt_stability_limit = func_to_solve(dt_stability_limit);
    end

end

repo: Aero-matlab


text:
function geneID = findGeneIDs(model, geneList)
% Finds gene numbers in a model
%
% USAGE:
%
%    geneID = findGeneIDs(model, geneList)
%
% INPUTS:
%    model:       COBRA model structure
%    geneList:    List of genes
%
% OUTPUT:
%    geneID:      List of gene IDs corresponding to `geneList`
%
% .. Author: - Jeff Orth 7/16/09

if (iscell(geneList))
    [tmp,geneID] = ismember(geneList,model.genes);
else
    geneID = find(strcmp(model.genes,geneList));
    if (isempty(geneID))
        geneID = 0;
    end
    if (length(geneID) > 1)
        geneID = geneID(1);
    end
end

author: opencobra

repo: cobratoolbox


text:
classdef TopOptComputer < handle

    properties (Access = public)
        computation
        settings
    end

    properties (Access = private)
        testName
    end

    methods (Access = public)
        function obj = TopOptComputer(cParams)
            obj.testName = cParams.testName;
        end

        function compute(obj)
            obj.createSettings();
            topOptSolver = TopOpt_Problem(obj.settings);
            topOptSolver.computeVariables();
            obj.computation = topOptSolver;
        end
    end

    methods (Access = private)

        function createSettings(obj)
            obj.createOldSettings();
            obj.translateToNewSettings();
        end

        function createOldSettings(obj)
            fileName = obj.testName;
            s = Settings(fileName);
            s.warningHoleBC = false;
            s.printIncrementalIter = false;
            s.printChangingFilter = false;
            s.printing = false;
            obj.settings = s;
        end

        function translateToNewSettings(obj)
            translator = SettingsTranslator();
            translator.translate(obj.settings);
            fileName = translator.fileName;
            obj.settings  = SettingsTopOptProblem(fileName);
        end

    end
end


text:
function test_suite = test_graphDiameter
%TEST_GRAPHDIAMETER  One-line description here, please.
%
%   output = test_graphDiameter(input)
%
%   Example
%   testGraphDiameter
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-05-19,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testGraph02(testCase) %#ok<*DEFNU>

[nodes, edges] = createTestGraph02;
diam = graphDiameter(nodes, edges);

exp = 6;
testCase.assertEqual(exp, diam);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function varargout = quickvis(varargin)
% VL_QUICKVIS Create an edge image from a Quickshift segmentation.
%   IEDGE = VL_QUICKVIS(I, RATIO, KERNELSIZE, MAXDIST, MAXCUTS) creates an edge
%   stability image from a Quickshift segmentation. RATIO controls the tradeoff
%   between color consistency and spatial consistency (See VL_QUICKSEG) and
%   KERNELSIZE controls the bandwidth of the density estimator (See VL_QUICKSEG,
%   VL_QUICKSHIFT). MAXDIST is the maximum distance between neighbors which
%   increase the density.
%
%   VL_QUICKVIS takes at most MAXCUTS thresholds less than MAXDIST, forming at
%   most MAXCUTS segmentations. The edges between regions in each of these
%   segmentations are labeled in IEDGE, where the label corresponds to the
%   largest DIST which preserves the edge.
%
%   [IEDGE,DISTS] = VL_QUICKVIS(I, RATIO, KERNELSIZE, MAXDIST, MAXCUTS) also
%   returns the DIST thresholds that were chosen.
%
%   IEDGE = VL_QUICKVIS(I, RATIO, KERNELSIZE, DISTS) will use the DISTS
%   specified
%
%   [IEDGE,DISTS,MAP,GAPS] = VL_QUICKVIS(I, RATIO, KERNELSIZE, MAXDIST, MAXCUTS)
%   also returns the MAP and GAPS from VL_QUICKSHIFT.
%
%   See Also: VL_QUICKSHIFT(), VL_QUICKSEG(), VL_HELP().
[varargout{1:nargout}] = vl_quickvis(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function insertPlan()
%
%Insert a plan from another study into the currently open study.
%A check is made that the CT matches and the plan is not already in the current study.
%
%Last Modified: 20 Mar 2006    KU
%                8 July 06,    KU,     Modified check and renaming of duplicate fractionGroupID names.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

global planC stateS

%Start diary and timer for import log.
startTime = now;
tmpFileName = tempname;
diary(tmpFileName);

%Load study containing the plan(s) to be inserted.
if isfield(stateS, 'CERRFile') && ~isempty(stateS.CERRFile)
    if stateS.workspacePlan
        %If workspace plan, ie no directory, use CERR root.
        stateS.CERRFile = fullfile(getCERRPath, 'workspacePlan');
    end
    dir = fileparts(stateS.CERRFile);
    wd = cd;
    cd(dir);
    [fname, pathname] = uigetfile('*.mat;*.bz2', ...
    'Select the study containing the plan(s) to insert.','Location',[100,100]);  %at position 100, 100 in pixels.
    cd(wd);
else
    [fname, pathname] = uigetfile('*.mat;*.bz2', ...
    'Select the study containing the plan(s) to insert.','Location',[100,100]);  %at position 100, 100 in pixels.
end

if fname == 0
    disp('No file selected.');
    return
end

file = [pathname fname];

[temp_planC] = loadPlanC_temp(fname,file);

%check if CT is the same as CT in current study.
try
    if ~isequal(temp_planC{temp_planC{end}.scan}(1).scanArray, planC{planC{end}.scan}(1).scanArray)
        sentence1 = 'CT of selected study does not match CT of current study. No plans can be inserted.';
        Zmsgbox=msgbox(sentence1, 'modal');
        waitfor(Zmsgbox);
        clear temp_planC;
        return
    end
catch
    sentence1 = 'CT of selected study does not match CT of current study. No plans can be inserted.';
    Zmsgbox=msgbox(sentence1, 'modal');
    waitfor(Zmsgbox);
    clear temp_planC;
    return
end

%Check that there are plans in temp_planC.
if isempty(temp_planC{temp_planC{end}.dose}) ||...
        (isempty(temp_planC{temp_planC{end}.dose}(1).doseArray) && isempty(temp_planC{temp_planC{end}.dose}(1,end).doseArray))
    sentence1 = 'There are no plans in the selected study.';
    Zmsgbox=msgbox(sentence1, 'modal');
    waitfor(Zmsgbox);
    clear temp_planC;
    return
end

addedplans = 0;
insert = 1;
while insert == 1

    %Select plans to be inserted.
    str = {temp_planC{temp_planC{end}.dose}.fractionGroupID};
    [planNum,ok] = listdlg('PromptString','Select a plan to insert:',...
                    'ListSize',[200 200],'SelectionMode','single',...
                    'ListString',str, 'OKString', 'Insert Plan');

    if ok == 0
        insert = 0;
        continue
    end

    %Check if selected plan is already in current study.
    match = 0;
    if length(planC{planC{end}.dose}) >=1
        for k=1:length(planC{planC{end}.dose})
            try
                if isequal(temp_planC{temp_planC{end}.dose}(1,planNum).doseArray, planC{planC{end}.dose}(1,k).doseArray)
                    match = 1;
                    sentence1 = 'Selected plan is already in current study.';
                    Zmsgbox=msgbox(sentence1, 'modal');
                    waitfor(Zmsgbox);
                end
            catch
            end
        end
    end

    %If selected plan is not yet in current study:
    if match ~= 1

        % Check dose and beam geometry for duplicate fractionGroupID names.
        if ischar(temp_planC{temp_planC{end}.dose}(1,planNum).fractionGroupID)
            ID = temp_planC{temp_planC{end}.dose}(1,planNum).fractionGroupID;
        else
            ID = num2str(temp_planC{temp_planC{end}.dose}(1,planNum).fractionGroupID);
        end
        newPlanName = ID;

        % Check dose for duplicate fractionGroupID names.
        if length(planC{planC{end}.dose}) >=1
            n = 0;
            matchname = 1;
            while matchname == 1
                matchname = 0;
                for k=1:length(planC{planC{end}.dose})
                    if ischar(planC{planC{end}.dose}(1,k).fractionGroupID)
                        idString = planC{planC{end}.dose}(1,k).fractionGroupID;
                    else
                        idString = num2str(planC{planC{end}.dose}(1,k).fractionGroupID);
                    end

                    if strcmpi(idString, newPlanName)
                        matchname = 1;
                        n=n+1;
                        newPlanName = [ID,' (',num2str(n),')'];
                    end
                end
            end
            temp_planC{temp_planC{end}.dose}(1,planNum).fractionGroupID = newPlanName;
            if n >= 1
                disp(['A plan with the name "', ID, '" already exists in the current study.']);
                disp(['Changing new plan name to "',  newPlanName, '".']);
                sentence1 = (['A plan with the name "', ID, '" already exists in the current study. ',...
                            'Changing new plan name to "',  newPlanName, '".']);
                Zmsgbox=msgbox(sentence1, 'modal');
                beep
                waitfor(Zmsgbox);
            end

           % Check beam geometry for duplicate fractionGroupID names.
            if ~strcmpi(ID, newPlanName)
                for i=1:length(temp_planC{temp_planC{end}.beamGeometry})
                    if ischar(temp_planC{temp_planC{end}.beamGeometry}(1,i).fractionGroupID)
                        tempID = temp_planC{temp_planC{end}.beamGeometry}(1,i).fractionGroupID;
                    else
                        tempID = num2str(temp_planC{temp_planC{end}.beamGeometry}(1,i).fractionGroupID);
                    end
                    if strcmpi(tempID, ID)
                        temp_planC{temp_planC{end}.beamGeometry}(1,i).fractionGroupID = newPlanName;
                    end
                end
            end
        end

        % Add associated scan and dose UID's to dose UID fields.
        scanUID = planC{planC{end}.scan}(1).scanUID;
        temp_planC{temp_planC{end}.dose}(1,planNum).assocScanUID = scanUID;
        temp_planC{temp_planC{end}.dose}(1,planNum).doseUID = createUID('dose');

        % Add new plan to planC.
        if isempty(planC{planC{end}.dose}) || isempty(planC{planC{end}.dose}(1,end).doseArray)     %no existing plan in planC

            planC{planC{end}.dose} = temp_planC{temp_planC{end}.dose}(1,planNum);    %add dose

            ok = 0;
            for i=1:length(temp_planC{temp_planC{end}.beamGeometry})     %add beams
                if ischar(temp_planC{temp_planC{end}.beamGeometry}(1,i).fractionGroupID)
                    tempID = temp_planC{temp_planC{end}.beamGeometry}(1,i).fractionGroupID;
                else
                    tempID = num2str(temp_planC{temp_planC{end}.beamGeometry}(1,i).fractionGroupID);
                end
                if strcmpi(tempID, newPlanName)
                    ok = ok+1;
                    if ok == 1
                        planC{planC{end}.beamGeometry} = temp_planC{temp_planC{end}.beamGeometry}(1,i);
                    else
                        planC{planC{end}.beamGeometry} = [planC{planC{end}.beamGeometry}, temp_planC{temp_planC{end}.beamGeometry}(1,i)];
                    end
                end
            end

        else        %at least one plan already in planC
            if isfield(planC{planC{end}.dose}(1),'cachedMask')
                temp_planC{temp_planC{end}.dose}(1,planNum).cachedMask = [];
            end
            planC{planC{end}.dose} = dissimilarInsert(planC{planC{end}.dose}, temp_planC{temp_planC{end}.dose}(1,planNum));     %add dose

            for i=1:length(temp_planC{temp_planC{end}.beamGeometry})     %add beams
                if ischar(temp_planC{temp_planC{end}.beamGeometry}(1,i).fractionGroupID)
                    tempID = temp_planC{temp_planC{end}.beamGeometry}(1,i).fractionGroupID;
                else
                    tempID = num2str(temp_planC{temp_planC{end}.beamGeometry}(1,i).fractionGroupID);
                end
                if strcmpi(tempID, newPlanName)
                    planC{planC{end}.beamGeometry} = dissimilarInsert(planC{planC{end}.beamGeometry}, temp_planC{temp_planC{end}.beamGeometry}(1,i));
                end
            end
        end

        try   %Need try statement in case original study was in RTOG format, not Dicom RT.
            try
                dose_info = temp_planC{temp_planC{end}.dose}(1,planNum).DICOMHeaders;
            catch
                dose_info = temp_planC{temp_planC{end}.dose}(1,planNum).DICOMHeaders{1,1};
            end
            for i = 1:length(temp_planC{temp_planC{end}.beams})
                plan_info = temp_planC{temp_planC{end}.beams}(1,i);
                if strcmpi(dose_info.ReferencedRTPlanSequence.Item_1.ReferencedSOPInstanceUID, plan_info.SOPInstanceUID)
                    if isempty(planC{planC{end}.beams})
                        planC{planC{end}.beams} = temp_planC{temp_planC{end}.beams}(1,i);
                    else
                        planC{planC{end}.beams} = dissimilarInsert(planC{planC{end}.beams}, temp_planC{temp_planC{end}.beams}(1,i));
                    end
                    break
                end
            end
        end

        planC{planC{end}.dose}(1,end).doseNumber = length(planC{planC{end}.dose});

        %switch to new dose after a short pause.
        pause(.1);
        sliceCallBack('selectDose', num2str(planC{planC{end}.dose}(1,end).doseNumber));

        %switch to scan #1 if not already selected.
        sliceCallBack('selectScan', num2str(1));

        %Turn on dose and switch to colorwash if not already on.
        if stateS.doseToggle == -1
            stateS.doseToggle = - stateS.doseToggle;
            set(findobj('tag', 'doseToggle'), 'checked', 'on');
        end
        stateS.optS.dosePlotType = 'colorwash';
        stateS.doseChanged = 1;
        stateS.doseSetChanged = 1;
        stateS.doseDisplayChanged = 1;
        stateS.structsChanged = 1;
        sliceCallBack('refresh')

        addedplans = addedplans + 1;
        display(['The plan ', ID, ' from study ', fname, ' has been added to the current study.'])

        sentence1 = horzcat('The selected plan has been added to the current study. ',...
            'The changes have not been saved.  Please remember to save changes ',...
            'before closing the study.');
        Zmsgbox=msgbox(sentence1, 'modal');
        waitfor(Zmsgbox);

    end

    answer = questdlg('Insert another plan?', 'Insert plan?', 'Yes', 'No', 'Yes');
    switch answer
        case 'Yes'
            insert = 1;
        case 'No'
            insert = 0;
        otherwise
    end

end

indexS = planC{end};
%Check dose-grid
for doseNum = 1:length(planC{indexS.dose})
    if planC{indexS.dose}(doseNum).zValues(2) - planC{indexS.dose}(doseNum).zValues(1) < 0
        planC{indexS.dose}(doseNum).zValues = flipud(planC{indexS.dose}(doseNum).zValues);
        planC{indexS.dose}(doseNum).doseArray = flipdim(planC{indexS.dose}(doseNum).doseArray,3);
    end
end

%Stop diary and write it to planC.
diary off;
endTime = now;
logC = file2cell(tmpFileName);
delete(tmpFileName);
indexS = planC{end};

if addedplans ~= 0
    temp_planC{temp_planC{end}.importLog} = [];
    temp_planC{temp_planC{end}.importLog}.importLog = logC;
    temp_planC{temp_planC{end}.importLog}.startTime = datestr(startTime);
    temp_planC{temp_planC{end}.importLog}.endTime = datestr(endTime);

    planC{indexS.importLog} = [planC{indexS.importLog}, temp_planC{temp_planC{end}.importLog}];
end
clear temp_planC;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [pnt, tri] = mesh_cube()

% MESH_CUBE creates a triangulated cube
%
% Use as
%   [pos, tri] = mesh_cube()
%
% See also MESH_TETRAHEDRON, MESH_OCTAHEDRON, MESH_ICOSAHEDRON, MESH_SPHERE, MESH_CONE

% Copyright (C) 2019, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.

pnt = [
  -1 -1 -1
  -1  1 -1
   1  1 -1
   1 -1 -1
  -1 -1  1
  -1  1  1
   1  1  1
   1 -1  1
  ];

tri = [
  1 2 4
  2 3 4
  1 5 6
  1 6 2
  2 6 7
  2 7 3
  3 7 8
  3 8 4
  4 8 5
  4 5 1
  5 6 8
  6 7 8
  ];

author: spm


text:
% TEST_MAXWELL_THICK_RING_G_NMNN: data function for Neumann boundary condition.

function g = test_maxwell_thick_ring_g_nmnn (x, y, z, ind)

  [theta, r] = cart2pol (x, y);
  g = zeros ([3, size(x)]);
  switch (ind)
    case 1
      g(1,:,:) = sin(theta) .* (cos(x) - cos(y));
      g(2,:,:) = cos(theta) .* (cos(y) - cos(x));
      g(3,:,:) = -x .* sin(theta) - y .* cos(theta);
    case 2
      g(1,:,:) = sin(theta) .* (cos(y) - cos(x));
      g(2,:,:) = cos(theta) .* (cos(x) - cos(y));
      g(3,:,:) = x .* sin(theta) + y .* cos(theta);
    case 3
      g(1,:,:) = cos(x) - cos(y);
    case 4
      g(2,:,:) = cos(y) - cos(x);
    case 5
      g(1,:,:) = y;
      g(2,:,:) = x;
    case 6
      g(1,:,:) = -y;
      g(2,:,:) = -x;
    otherwise
      error ('g_nmnn: unknown reference number')
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function out= bbci_calibrate_evalFeature(signal, bbci_feature, mrk)
%BBCI_CALIBRATE_EVALFEATURE - Perform feature extration
%
%Synopsis:
%  FV= bbci_calibrate_evalFeature(EPO, BBCI_FEATURE)
%  FV= bbci_calibrate_evalFeature(CNT, BBCI_FEATURE, MRK)
%
%Arguments:
%  EPO - Structure of epoched data.
%  BBCI_FEATURE - Struct (array) specifying the feature extraction,
%      subfield of 'bbci' structure of bbci_apply.
%  CNT - Structure of continuous data.
%  MRK - Structure of markers.
%
%Output:
%  FV - Structure holding the extracted feature vector

% 11-2011 Benjamin Blankertz

bbci_feature= bbciutil_transformProc2FcnParam(bbci_feature);
bbci_feature= opt_setDefaults(bbci_feature, {'signal', 1});

feature= repmat(struct, [1 length(bbci_feature)]);
for f= 1:length(bbci_feature),
  BF= bbci_feature(f);
  fv= signal(BF.signal);
  if nargin>=3,
    fv= proc_segmentation(fv, mrk, BF.ival);
  end

  if f==1 && isfield(fv, 'y'),
    out.y= fv.y;
  end

  for k= 1:length(BF.fcn),
    fv= BF.fcn{k}(fv, BF.param{k}{:});
  end
  % clash all feature dimensions (but preserve sample dimension)
  sz= size(fv.x);
  feature(f).x= reshape(fv.x, [prod(sz(1:end-1)) sz(end)]);
end
out.x= cat(1, feature.x);

author: bbci


text:
function formatsvnrevision(v1, v2)

disp('Getting revision information...');
system(sprintf('svn log --verbose -r%d:%d > tmprev.txt', v1, v2));

oldpwd = pwd;
newpwd = fileparts(which('eeglab.m'));
cd(newpwd);

fid = fopen('tmprev.txt', 'r');

state = 1;
while ~feof(fid)
    txt = fgetl(fid);
    if length(txt) > 4 && strcmpi(txt(1:4), '----')
        state = 1;
        if feof(fid), return; end;
    end;
    if state == 1,
        % get rev number
        txt  = fgetl(fid);
        rev = strtok(txt);
        revnum = rev(2:end);

        % get function name
        txt = fgetl(fid);
        txt = fgetl(fid);
        [tmp file] = strtok(txt);
        file = deblank(file);
        [tmp file ext] = fileparts(file);
        file = [ file ext ];
        state = 2;
    end;
    if isempty(txt)
        state = 3;
    end;
    if state == 3
        revinfo = fgetl(fid);
        revinfo = revinfo(1:end);
        state = 1;
        fprintf('** %s, %s (SVN %s - Arno)\n', file, revinfo, revnum);
    end;
end;
cd(oldpwd);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
clear all; close all;
I=imread('circles.png');
J=bwmorph(I, 'remove');
figure;
subplot(121);  imshow(I);
subplot(122);  imshow(J);

author: HuangCongQing


text:
function [shape] = ft_read_headshape(filename, varargin)

% FT_READ_HEADSHAPE reads the fiducials and/or the measured headshape
% from a variety of files (like CTF and Polhemus). The headshape and
% fiducials can for example be used for coregistration.
%
% Use as
%   [shape] = ft_read_headshape(filename)
%
% See also FT_READ_VOL, FT_READ_SENS, FT_WRITE_HEADSHAPE

% Copyright (C) 2008-2010 Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: ft_read_headshape.m 2914 2011-02-21 16:08:36Z jansch $

% check the input: if filename is a cell-array, call ft_read_headshape recursively and combine the outputs
if iscell(filename)
  for i=1:numel(filename)
    tmpshape = ft_read_headshape(filename{i}, varargin{:});
    if i==1,
      shape = tmpshape;
    else
      tmpshape  = ft_convert_units(tmpshape, shape.unit);
      npnt      = size(shape.pnt,1);
      shape.pnt = cat(1, shape.pnt, tmpshape.pnt);
      if isfield(shape, 'tri') && isfield(tmpshape, 'tri')
        shape.tri = cat(1, shape.tri, tmpshape.tri + npnt);
      elseif ~isfield(shape, 'tri') && ~isfield(tmpshape, 'tri')
        % this is ok
      else
        error('not all input files seem to contain a triangulation');
      end
    end
  end
  return
end

% test whether the file exists
if ~exist(filename)
  error(sprintf('file ''%s'' does not exist', filename));
end

% get the options
fileformat  = keyval('format',      varargin);
coordinates = keyval('coordinates', varargin); if isempty(coordinates), coordinates = 'head'; end

if isempty(fileformat)
  fileformat = ft_filetype(filename);
end

% start with an empty structure
shape           = [];
shape.pnt       = [];
shape.fid.pnt   = [];
shape.fid.label = {};

switch fileformat
  case {'ctf_ds', 'ctf_hc', 'ctf_meg4', 'ctf_res4', 'ctf_old'}
    [p, f, x] = fileparts(filename);

    if strcmp(fileformat, 'ctf_old')
      fileformat = ft_filetype(filename);
    end

    if strcmp(fileformat, 'ctf_ds')
      filename = fullfile(p, [f x], [f '.hc']);
    elseif strcmp(fileformat, 'ctf_meg4')
      filename = fullfile(p, [f '.hc']);
    elseif strcmp(fileformat, 'ctf_res4')
      filename = fullfile(p, [f '.hc']);
    end

    orig = read_ctf_hc(filename);
    switch coordinates
      case 'head'
        shape.fid.pnt = cell2mat(struct2cell(orig.head));
      case 'dewar'
        shape.fid.pnt = cell2mat(struct2cell(orig.dewar));
      otherwise
        error('incorrect coordinates specified');
    end
    shape.fid.label = fieldnames(orig.head);

  case 'ctf_shape'
    orig = read_ctf_shape(filename);
    shape.pnt = orig.pnt;
    shape.fid.label = {'NASION', 'LEFT_EAR', 'RIGHT_EAR'};
    for i = 1:numel(shape.fid.label)
      shape.fid.pnt = cat(1, shape.fid.pnt, ...
        getfield(orig.MRI_Info, shape.fid.label{i}));
    end

  case {'4d_xyz', '4d_m4d', '4d_hs', '4d', '4d_pdf'}
    [p, f, x] = fileparts(filename);
    if ~strcmp(fileformat, '4d_hs')
      filename = fullfile(p, 'hs_file');
    end
    [shape.pnt, fid] = read_bti_hs(filename);

    % I'm making some assumptions here
    % which I'm not sure will work on all 4D systems

    %fid = fid(1:3, :);

    [junk, NZ] = max(fid(1:3,1));
    [junk, L]  = max(fid(1:3,2));
    [junk, R]  = min(fid(1:3,2));
    rest       = setdiff(1:size(fid,1),[NZ L R]);

    shape.fid.pnt = fid([NZ L R rest], :);
    shape.fid.label = {'NZ', 'L', 'R'};
    if ~isempty(rest),
      for i = 4:size(fid,1)
        shape.fid.label{i} = ['fiducial' num2str(i)];
        % in a 5 coil configuration this corresponds with Cz and Inion
      end
    end

  case 'itab_asc'
    shape = read_itab_asc(filename);

  case 'neuromag_mex'
    [co,ki,nu] = hpipoints(filename);
    fid = co(:,find(ki==1))';

    [junk, NZ] = max(fid(:,2));
    [junk, L]  = min(fid(:,1));
    [junk, R]  = max(fid(:,1));

    shape.fid.pnt = fid([NZ L R], :);
    shape.fid.label = {'NZ', 'L', 'R'};

  case {'mne_source'}
    % read the source space from an MNE file
    ft_hastoolbox('mne', 1);

    src = mne_read_source_spaces(filename, 1);
    shape = [];
    % only keep the points that are in use
    inuse1 = src(1).inuse==1;
    inuse2 = src(2).inuse==1;
    shape.pnt=[src(1).rr(inuse1,:); src(2).rr(inuse2,:)];

    % only keep the triangles that are in use; these have to be renumbered
    newtri1 = src(1).use_tris;
    newtri2 = src(2).use_tris;
    for i=1:numel(src(1).vertno)
      newtri1(newtri1==src(1).vertno(i)) = i;
    end
    for i=1:numel(src(2).vertno)
      newtri2(newtri2==src(2).vertno(i)) = i;
    end
    shape.tri  = [newtri1; newtri2 + numel(src(1).vertno)];
    shape.orig.pnt = [src(1).rr; src(2).rr];
    shape.orig.tri = [src(1).tris; src(2).tris + src(1).np];
    shape.orig.inuse = [src(1).inuse src(2).inuse]';

  case {'neuromag_mne', 'neuromag_fif'}
    % read the headshape and fiducials from an MNE file
    hdr = ft_read_header(filename,'headerformat','neuromag_mne');
    nFid = size(hdr.orig.dig,2); %work out number of fiducials
    switch coordinates
      case 'head' % digitiser points should be stored in head coordinates by default

        fidN=1;
        pntN=1;
        for i=1:nFid %loop over fiducials
          %check this point is in head coordinates:
          if hdr.orig.dig(i).coord_frame~=4 % 4 is MNE constant for head coordinates
            error(['Digitiser point (' num2str(i) ') not stored in head coordinates!']);
          end

          switch hdr.orig.dig(i).kind % constants defined in MNE - see p.215 of MNE manual
            case 1 % Cardinal point (nasion, LPA or RPA)
              %get location of fiducial:
              shape.fid.pnt(fidN,1:3) = hdr.orig.dig(i).r*100; %multiply by 100 to convert to cm
              switch hdr.orig.dig(i).ident
                case 1 % LPA
                  shape.fid.label{fidN} = 'LPA';
                case 2 % nasion
                  shape.fid.label{fidN} = 'Nasion';
                case 3 % RPA
                  shape.fid.label{fidN} = 'RPA';
                otherwise
                  error('Unidentified cardinal point in file!');
              end
              fidN = fidN + 1;

            case 2 % HPI coil
              shape.pnt(pntN,1:3) = hdr.orig.dig(i).r*100;
              pntN = pntN + 1;
            case 3 % EEG electrode location (or ECG)
              shape.pnt(pntN,1:3) = hdr.orig.dig(i).r*100;
              pntN = pntN + 1;
            case 4 % Additional head point
              shape.pnt(pntN,1:3) = hdr.orig.dig(i).r*100;
              pntN = pntN + 1;
            otherwise
              warning('Unidentified digitiser point in file!');
          end

        end
        shape.fid.label=shape.fid.label';

      case 'dewar'
        error('Dewar coordinates not supported for headshape yet (MNE toolbox)');
      otherwise
        error('Incorrect coordinates specified');
    end

  case {'yokogawa_mrk', 'yokogawa_ave', 'yokogawa_con', 'yokogawa_raw' }

    hdr = read_yokogawa_header(filename);
    marker = hdr.orig.matching_info.marker;

    % markers 1-3 identical to zero: try *.mrk file
    if ~any([marker(:).meg_pos])
      [p, f, x] = fileparts(filename);
      filename = fullfile(p, [f '.mrk']);
      if exist(filename, 'file')
        hdr    = read_yokogawa_header(filename);
        marker = hdr.orig.matching_info.marker;
      end
    end

    % non zero markers 1-3
    if any([marker(:).meg_pos])
      shape.fid.pnt = cat(1, marker(1:5).meg_pos);
      sw_ind = [3 1 2];
      shape.fid.pnt(1:3,:)= shape.fid.pnt(sw_ind, :);
      shape.fid.label = {'nas'; 'lpa'; 'rpa'; 'Marker4'; 'Marker5'};
    else
      error('no coil information found in Yokogawa file');
    end

    % Convert to the units of the grad.
    shape = ft_convert_units(shape, 'cm');

  case 'yokogawa_coregis'

    in_str = textread(filename, '%s');
    nr_items = size(in_str,1);
    ind = 1;
    coil_ind = 1;
    while ind < nr_items
      if strcmp(in_str{ind},'MEG:x=')
        shape.fid.pnt = [shape.fid.pnt; str2num(strtok(in_str{ind+1},[',','['])) ...
          str2num(strtok(in_str{ind+3},[',','['])) str2num(strtok(in_str{ind+5},[',','[']))];
        shape.fid.label = [shape.fid.label ; ['Marker',num2str(coil_ind)]];
        coil_ind = coil_ind + 1;
        ind = ind + 6;
      else
        ind = ind +1;
      end
    end
    if size(shape.fid.label,1) ~= 5
      error('Wrong number of coils');
    end

    sw_ind = [3 1 2];

    shape.fid.pnt(1:3,:)= shape.fid.pnt(sw_ind, :);
    shape.fid.label(1:3)= {'nas', 'lpa', 'rpa'};

  case 'polhemus_fil'
    [shape.fid.pnt, shape.pnt, shape.fid.label] = read_polhemus_fil(filename, 0);

  case 'spmeeg_mat'
    tmp = load(filename);
    if isfield(tmp.D, 'fiducials') && ~isempty(tmp.D.fiducials)
      shape = tmp.D.fiducials;
    else
      error('no headshape found in SPM EEG file');
    end

  case 'matlab'
    tmp = load(filename);
    if isfield(tmp, 'shape')
      shape = tmp.shape;
    elseif isfield(tmp, 'bnd')
      % the variable in the file is most likely a precomputed triangulation of some
      % sort
      shape = tmp.bnd;
    elseif isfield(tmp, 'elec')
      shape.fid.pnt   = tmp.elec.pnt;
      shape.fid.label = tmp.elec.label;
    else
      error('no headshape found in Matlab file');
    end

  case {'freesurfer_triangle_binary', 'freesurfer_quadrangle'}
    % the freesurfer toolbox is required for this
    ft_hastoolbox('freesurfer', 1);
    [pnt, tri] = read_surf(filename);
    if min(tri(:)) == 0
      % start counting from 1
      tri = tri + 1;
    end
    shape.pnt = pnt;
    shape.tri = tri;
    shape = rmfield(shape, 'fid');

  case 'mne_tri'
    % FIXME this should be implemented, consistent with ft_write_headshape
    keyboard

  case 'mne_pos'
    % FIXME this should be implemented, consistent with ft_write_headshape
    keyboard

  otherwise
    % try reading it from an electrode of volume conduction model file
    success = false;

    if ~success
      % try reading it as electrode positions
      % and treat those as fiducials
      try
        elec = ft_read_sens(filename);
        if ~ft_senstype(elec, 'eeg')
          error('headshape information can not be read from MEG gradiometer file');
        else
          shape.fid.pnt   = elec.pnt;
          shape.fid.label = elec.label;
          success = 1;
        end
      catch
        success = false;
      end % try
    end

    if ~success
      % try reading it as volume conductor
      % and treat the skin surface as headshape
      try
        vol = ft_read_vol(filename);
        if ~ft_voltype(vol, 'bem')
          error('skin surface can only be extracted from boundary element model');
        else
          if ~isfield(vol, 'skin')
            vol.skin = find_outermost_boundary(vol.bnd);
          end
          shape.pnt = vol.bnd(vol.skin).pnt;
          shape.tri = vol.bnd(vol.skin).tri; % also return the triangulation
          success = 1;
        end
      catch
        success = false;
      end % try
    end

    if ~success
      error('unknown fileformat for head shape information');
    end
end

if issubfield(shape, 'fid.label')
  % ensure that it is a column
  shape.fid.label = shape.fid.label(:);
end

% this will add the units to the head shape
shape = ft_convert_units(shape);

author: PatternRecognition

repo: OpenBMI


text:
% Copyright (C) 2017-2018 Titus Cieslewski, RPG, University of Zurich,
%   Switzerland
%   You can contact the author at <titus at ifi dot uzh dot ch>
% Copyright (C) 2017-2018 Siddharth Choudhary, College of Computing,
%   Georgia Institute of Technology, Atlanta, GA, USA
% Copyright (C) 2017-2018 Davide Scaramuzza, RPG, University of Zurich,
%   Switzerland
%
% This file is part of dslam_open.
%
% dslam_open 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.
%
% dslam_open 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 dslam_open. If not, see <http://www.gnu.org/licenses/>.

function [informationExchangeGaussSeidel, nrGaussSeidelIterations, ...
    informationExchangeDDF, commLinks] = getDistGaussSeidelStats(dir)

filename = [dir '/trace'];
trace_file = sprintf('%s_overall_error.txt',filename);
trace = dlmread(trace_file, ' ');

end_index = find(trace(1,:) == -1);
rotation_iterations = end_index-1;

end_index = find(trace(2,:) == -1);
pose_iterations = end_index-1;

centralized_error = trace(3,1);
distributed_error = trace(4,1);

sizePose = 6*8;  % Assumed to be 6 doubles, 1 double = 8 byte
sizeRotation = 9*8;  % Assumed to be 9 doubles, 1 double = 8 byte

commLinks = 2 * trace(5,1);

% Gauss-Seidel Communication
nrGaussSeidelIterations = pose_iterations + rotation_iterations;
informationExchangeGaussSeidel = pose_iterations.*(commLinks*sizePose) + rotation_iterations.*(commLinks*sizeRotation);

% Compare with DDF-SAM Communication
nrGNIterations = 1; % typically 3
sizePoseDDFSAM = 6*8;  % Assumed to be 12 doubles, 1 double = 8 byte
informationExchangeDDF = nrGNIterations.*((commLinks*sizePoseDDFSAM) + (commLinks*sizePoseDDFSAM).^2);

author: uzh-rpg


text:
classdef DMUtilities
    % DMUtilities. Utility functions related to Dicom files
    %
    %
    %     Licence
    %     -------
    %     Part of DicoMat. https://github.com/tomdoel/dicomat
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the BSD 3-Clause license. Please see the file LICENSE for details.
    %

    methods (Static)

        function is_dicom = IsDicom(file_path, file_name)
            % Returns true if this is a Dicom file

            if strcmp(file_name, 'DICOMDIR')
                is_dicom = false;
                return
            end

            full_file_name = [file_path, filesep, file_name];

            is_dicom = DMFallbackDicomLibrary.getLibrary.isdicom(full_file_name);
        end

        function dicom_series_uid = GetDicomSeriesUid(fileName, dictionary)
            % Gets the series UID for a Dicom file

            if isempty(dictionary)
                dictionary = DMDicomDictionary.GroupingDictionary;
            end

            header = DMFallbackDicomLibrary.getLibrary.dicominfo(fileName, dictionary);

            if isempty(header)
                dicom_series_uid = [];
            else
                % If no SeriesInstanceUID tag then this is not a valid Dicom image (it
                % might be a DICOMDIR)
                if isfield(header, 'SeriesInstanceUID')
                    dicom_series_uid = header.SeriesInstanceUID;
                else
                    dicom_series_uid = [];
                end
            end
        end

        function metadata = ReadMetadata(fileName, dictionary, reporting)
            % Reads in Dicom metadata from the specified file
            try
                metadata = DMFallbackDicomLibrary.getLibrary.dicominfo(fileName, dictionary);
            catch exception
                reporting.Error('DMUtilities:MetaDataReadFail', ['Could not read metadata from the Dicom file ' file_name '. Error:' exception.message]);
            end
        end

        function metadata = ReadGroupingMetadata(fileName, reporting)
            % Reads in Dicom metadata from the specified file
            try
                metadata = DMFallbackDicomLibrary.getLibrary.dicominfo(fileName, DMDicomDictionary.GroupingDictionary);
            catch exception
                reporting.Error('DMUtilities:MetaDataReadFail', ['Could not read metadata from the Dicom file ' file_name '. Error:' exception.message]);
            end
        end

        function metadata = ReadEssentialMetadata(fileName, reporting)
            % Reads in Dicom metadata from the specified file
            try
                metadata = DMFallbackDicomLibrary.getLibrary.dicominfo(fileName);
            catch exception
                reporting.Error('DMUtilities:MetaDataReadFail', ['Could not read metadata from the Dicom file ' file_name '. Error:' exception.message]);
            end
        end

        function image_data = ReadDicomImageFromMetadata(metadata, reporting)
            % Reads in Dicom image data from the specified metadata

            try
                image_data = DMFallbackDicomLibrary.getLibrary.dicomread(metadata);
            catch exception
                reporting.Error('DMUtilities:DicomReadError', ['Error while reading the Dicom file. Error:' exception.message]);
            end
        end

        function ReadDicomImageIntoWrapperFromMetadata(metadata, image_wrapper, slice_index, reporting)
            % Reads in Dicom image data from the specified metadata. The image data
            % is stored directly into the RawImage matrix of a CoreWrapper object
            try
                image_wrapper.RawImage(:, :, slice_index, :) = DMFallbackDicomLibrary.getLibrary.dicomread(metadata);

            catch exception
                reporting.Error('DMUtilities:DicomReadError', ['Error while reading the Dicom file. Error:' exception.message]);
            end
        end

        function [name, short_name] = PatientNameToString(patient_name)
            if ischar(patient_name)
                name = patient_name;
            else
                name = '';
                short_name = '';
                if isstruct(patient_name)
                    name = DMUtilities.AddOptionalField(name, patient_name, 'FamilyName', false);
                    name = DMUtilities.AddOptionalField(name, patient_name, 'GivenName', false);
                    name = DMUtilities.AddOptionalField(name, patient_name, 'MiddleName', false);
                    name = DMUtilities.AddOptionalField(name, patient_name, 'NamePrefix', false);
                    name = DMUtilities.AddOptionalField(name, patient_name, 'NameSuffix', false);

                    short_name = DMUtilities.AddOptionalField(short_name, patient_name, 'FamilyName', true);
                    short_name = DMUtilities.AddOptionalField(short_name, patient_name, 'GivenName', true);
                    short_name = DMUtilities.AddOptionalField(short_name, patient_name, 'MiddleName', true);
                    short_name = DMUtilities.AddOptionalField(short_name, patient_name, 'NamePrefix', true);
                    short_name = DMUtilities.AddOptionalField(short_name, patient_name, 'NameSuffix', true);
                end
            end
        end

        function new_text = AddOptionalField(text, struct_name, field_name, only_if_nonempty)
            if isempty(text) || ~only_if_nonempty
                new_text = text;
                if isfield(struct_name, field_name) && ~isempty(struct_name.(field_name))
                    if isempty(text)
                        prefix = '';
                    else
                        prefix = ', ';
                    end
                    new_text = [text, prefix, struct_name.(field_name)];
                end
            else
                new_text = text;
            end
        end

        function uid = GetIdentifierFromFilename(file_name)
            [~, uid, ~] = fileparts(file_name);
        end

        function dicom_filenames = RemoveNonDicomFiles(image_path, filenames)
            dicom_filenames = [];
            for index = 1 : length(filenames)
                if (DMUtilities.IsDicom(image_path, filenames{index}))
                    dicom_filenames{end + 1} = filenames{index};
                end
            end
        end

        function tag_string = ConvertTag32ToString(tag_32)
            group_value = floor(tag_32/65536);
            tag_value = mod(tag_32, 65536);
            tag_string = ['(' dec2hex(group_value, 4) ',' dec2hex(tag_value, 4) ')'];
        end

        function tag_32 = ConvertStringToTag32(tag_string)
            group = tag_string(2:5);
            elem = tag_string(7:10);
            tag_32 = 65536*hex2dec(group) + hex2dec(elem);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
%{
    Copyright (C) Cambridge Electronic Design Limited 2014
    Author: James Thompson
    Web: www.ced.co.uk email: james@ced.co.uk, softhelp@ced.co.uk

    This file is part of CEDS64ML, a MATLAB interface to the SON64 library.

    CEDS64ML 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.

    CEDS64ML 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 CEDS64ML.  If not, see <http://www.gnu.org/licenses/>.
%}

function [ iOk, sUnitsOut ] = CEDS64ChanUnits( fhand, iChan, sUnitsIn )
%CEDS64CHANUNITS Gets and sets channels units
%   [ iOk {, sUnitsOut} ] = CEDS64ChanUnits( fhand, iChan{, sUnitsIn} )
%   Input
%   fhand - An integer handle to an open file
%   iChan - A channel number
%   sUnitsIn - (Optional) A string conatining the new units
%   Output
%   iOk - 0 if the units were set correctly otherwise a negative error code.
%   sUnitsOut - A string containing the old units

if (nargin == 2 || nargin == 3) % always get the old units
    %step 1 find out how big the comment is
    dummystring = blanks(1);
    [iSize] = calllib('ceds64int', 'S64GetChanUnits', fhand, iChan, dummystring, -1);
    %step 2 create a string buffer of the correct size
    stringptr = blanks(iSize+1);
    [iOk, sUnitsOut] = calllib('ceds64int', 'S64GetChanUnits', fhand, iChan, stringptr, 0);
else
    iOk = -22;
end

% has there been an error?
if iOk < 0
    return;
end

% if not set the new units if we're given some
if (nargin == 3)
    iOk = calllib('ceds64int', 'S64SetChanUnits', fhand, iChan, sUnitsIn);
end
end

author: brainstorm-tools


text:
function rmse_v = navego_rmse (nav, gnss, ref_n, ref_g)
% navego_rmse: calculates the Root Mean Squared Errors (RMSE) between
% a INS/GNSS system and a reference data structure, and between GNSS-only
% solution and a reference data structure.
%
% INPUT
%   nav_e, INS/GNSS integration data structure.
%   gnss,  GNSS data structure.
%   ref_n, Reference data structure ajusted for INS/GNSS estimations.
%   ref_g, Reference data structure ajusted for GNSS measurements.
%
% OUTPUT
%   rmse_v, vector with all RMSE.
%       RMSE_roll;  RMSE_pitch; RMSE_yaw; (degrees, degrees, degrees)
%       RMSE_vn;    RMSE_ve;    RMSE_vd;  (m/s, m/s, m/s)
%       RMSE_lat;   RMSE_lon;   RMSE_h;   (m, m, m)
%       RMSE_vn_g;  RMSE_ve_g;  RMSE_vd_g;(m/s, m/s, m/s)
%       RMSE_lat_g; RMSE_lon_g; RMSE_h_g; (m, m, m)
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% Version: 006
% Date:    2021/03/16
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

D2R = (pi/180);     % degrees to radians
R2D = (180/pi);     % radians to degrees

%% INS/GNSS ATTITUDE RMSE

RMSE_roll  = rmse (nav.roll , ref_n.roll)  .* R2D;
RMSE_pitch = rmse (nav.pitch, ref_n.pitch) .* R2D;

% Differences greater than 300 deg are avoided when comparing yaw angles.
% The idea is to avoid to compare values of yaw angles when, for example, the
% reference yaw is near pi and the nav yaw is near -pi. Both yaw angles, pi
% and -pi, represent the same heading angle (moving South).

nav.yaw = correct_yaw(nav.yaw);
ref_n.yaw = correct_yaw(ref_n.yaw);

idx = ( abs(nav.yaw - ref_n.yaw) < (300 * D2R) );
RMSE_yaw = rmse ( nav.yaw(idx), ref_n.yaw(idx) ) .* R2D;

%% INS/GNSS VELOCITY RMSE

if (isfield(nav, 'vel') && isfield(ref_n, 'vel'))
    RMSE_vn = rmse (nav.vel(:,1),  ref_n.vel(:,1));
    RMSE_ve = rmse (nav.vel(:,2),  ref_n.vel(:,2));
    RMSE_vd = rmse (nav.vel(:,3),  ref_n.vel(:,3));
else
    RMSE_vn = NaN;
    RMSE_ve = NaN;
    RMSE_vd = NaN;
    warning('navego_rmse: no NED velocity field was found in INS/GNSS data.');
end

%% INS/GNSS POSITION RMSE

[RM,RN] = radius(ref_n.lat);
LAT2M = (RM + ref_n.h);                     % Coefficient for lat, radians to meters
LON2M = (RN + ref_n.h) .* cos(ref_n.lat);   % Coefficient for lon, radians to meters

RMSE_lat = rmse (nav.lat.* LAT2M, ref_n.lat.* LAT2M) ;
RMSE_lon = rmse (nav.lon.* LON2M, ref_n.lon.* LON2M) ;
RMSE_h   = rmse (nav.h, ref_n.h);

%% GNSS VELOCITY RMSE

if (isfield(gnss, 'vel') && isfield( ref_g, 'vel'))
    RMSE_vn_g = rmse (gnss.vel(:,1), ref_g.vel(:,1));
    RMSE_ve_g = rmse (gnss.vel(:,2), ref_g.vel(:,2));
    RMSE_vd_g = rmse (gnss.vel(:,3), ref_g.vel(:,3));
else
    RMSE_vn_g = NaN;
    RMSE_ve_g = NaN;
    RMSE_vd_g = NaN;
    warning('navego_rmse: no NED velocity field was found in GNSS data.');
end

%% GNSS POSITION RMSE

[RMg,RNg] = radius(ref_g.lat);
LAT2Mg = (RMg + ref_g.h);                   % Coefficient for lat, radians to meters
LON2Mg = (RNg + ref_g.h) .* cos(ref_g.lat); % Coefficient for lon, radians to meters

RMSE_lat_g = rmse (gnss.lat.* LAT2Mg, ref_g.lat.* LAT2Mg) ;
RMSE_lon_g = rmse (gnss.lon.* LON2Mg, ref_g.lon.* LON2Mg) ;
RMSE_h_g   = rmse (gnss.h, ref_g.h);

%%

rmse_v = [  RMSE_roll;  RMSE_pitch; RMSE_yaw;
            RMSE_vn;    RMSE_ve;    RMSE_vd;
            RMSE_lat;   RMSE_lon;   RMSE_h;
            RMSE_vn_g;  RMSE_ve_g;  RMSE_vd_g;
            RMSE_lat_g; RMSE_lon_g; RMSE_h_g; ];
end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
%TRUNCQUAD_MODES  Compute colour modes for a truncated quadratic kernel
%
%   [modes depth energy inliers] = truncquad_modes(I, thresh[, use_variance[, search_width]])
%
% Computes the colour modes for each set of input vectors, using the
% truncated quadratic kernel.
%
% IN:
%   I - CxLxM array of L input colour vectors with C channels at M depths.
%   thresh - scalar cost truncation threshold.
%   use_variance - 0: sum cost over all input vectors; 1: sum cost over
%                  inlying input vectors / num_inliers; 2: sum cost over
%                  inliers / (num_inliers - 1). Default: 0.
%   search_width - scalar indicating depths above and below current depth
%                  to search for modes within. Default: M.
%
% OUT:
%   modes - CxN matrix of colour modes.
%   depth - 1xN list of indices of depth for each mode.
%   energy - 1xN list of costs for each mode.
%   inliers - LxN logical array indicating inliers used to calculate the
%             colour of each mode.

% $Id: truncquad_modes.m,v 1.2 2008/11/17 11:27:35 ojw Exp $

function varargout = truncquad_modes(varargin)
funcName = mfilename;
sourceList = {[funcName '.cxx']}; % Cell array of source files
vgg_mexcompile_script; % Compilation happens in this script
return

repo: SelfExSR


text:
function [w_out,ln_out] = straitln(term_mat,w,mode)
% STRAITLN Straight line approximation of bode response.
%          STRAITLN computes the straight line approximation of either the
%          magnitude response or phase response.  MODE = 1 computes the
%          magnitude response; MODE = 2 computes the phase response.

% Author: Craig Borghesani
% Date: 9/3/94
% Revised: 10/20/94
% Copyright (c) 1999, Prentice-Hall

% convert term_mat to extract lead/lag elements
term_mat = termcnvt(term_mat);
[num,den] = termextr(term_mat);
n=find(num~=0); d=find(den~=0);
dcgain = num(n(length(n)))/den(d(length(d)));
w_out = []; ln_out = [];

% pull out break frequency values and add to frequency vector
break_frq = [];
pole_zero = term_mat(find(term_mat(:,4)==4 | term_mat(:,4)==5),1);
wn = term_mat(find(term_mat(:,4)==6 | term_mat(:,4)==7),2);
break_frq = sort([abs([pole_zero(:)',wn(:)'])]);
if length(break_frq),
 break_frq = sort([break_frq,break_frq/10,break_frq*10]);
 break_frq(logical([0,diff(break_frq)==0]))=[];
end
wnew = sort([w,break_frq]);
wnew(logical([0,diff(wnew)==0]))=[];
[r,c]=size(term_mat);

lw = length(wnew);

% initialize magnitude straight line with gain
if mode == 1,
 ln_out = zeros(1,lw)+20*log10(abs(dcgain));
else
 ln_out = zeros(1,lw)-180*(term_mat(1,1)<0);
end

% handle integrators/differentiators first
first_w = wnew(1);
for k = 1:abs(term_mat(2,1)),
 if mode == 1,
  ln_out = ln_out - sign(term_mat(2,1))*20*log10(first_w);
  ln_out = ln_out - sign(term_mat(2,1))*20*(log10(wnew)-log10(first_w));
 else
  ln_out = ln_out - sign(term_mat(2,1))*90;
 end
end

% add in the rest of the terms
for k = 3:r,

% find break frequencies for both magnitude and phase plots
 if any(term_mat(k,4)==[4,5]),
  brk = find(wnew == abs(term_mat(k,1)));
 else
  brk = find(wnew == abs(term_mat(k,2)));
 end
 brk_w = wnew(brk);
 brkl = find(wnew == brk_w/10);
 brkl_w = wnew(brkl);
 brkr = find(wnew == brk_w*10);
 brkr_w = wnew(brkr);

 if any(term_mat(k,4)==[4,5]),

  if term_mat(k,4)==4,
   if mode == 1,
    ln_out(brk:lw) = ln_out(brk:lw) - 20*(log10(wnew(brk:lw))-log10(brk_w));
   else
    ln_out = ln_out - 180*(term_mat(k,1)<0);
    ln_out(brkl:brkr) = ln_out(brkl:brkr) - sign(term_mat(k,1))*45*(log10(wnew(brkl:brkr))-log10(brkl_w));
    ln_out((brkr+1):lw) = ln_out((brkr+1):lw) - 90*sign(term_mat(k,1));
   end
  else
   if mode == 1,
    ln_out(brk:lw) = ln_out(brk:lw) + 20*(log10(wnew(brk:lw))-log10(brk_w));
   else
    ln_out = ln_out - 180*(term_mat(k,1)<0);
    ln_out(brkl:brkr) = ln_out(brkl:brkr) + sign(term_mat(k,1))*45*(log10(wnew(brkl:brkr))-log10(brkl_w));
    ln_out((brkr+1):lw) = ln_out((brkr+1):lw) + 90*sign(term_mat(k,1));
   end
  end

 elseif any(term_mat(k,4)==[6,7]),
  if term_mat(k,4)==6,
   if mode == 1,
    ln_out(brk:lw) = ln_out(brk:lw) - 40*(log10(wnew(brk:lw))-log10(brk_w));
   else
    ln_out(brkl:brkr) = ln_out(brkl:brkr) - sign(term_mat(k,2))*90*(log10(wnew(brkl:brkr))-log10(brkl_w));
    ln_out((brkr+1):lw) = ln_out((brkr+1):lw) - 180*sign(term_mat(k,2));
   end
  else
   if mode == 1,
    ln_out(brk:lw) = ln_out(brk:lw) + 40*(log10(wnew(brk:lw))-log10(brk_w));
   else
    ln_out(brkl:brkr) = ln_out(brkl:brkr) + sign(term_mat(k,2))*90*(log10(wnew(brkl:brkr))-log10(brkl_w));
    ln_out((brkr+1):lw) = ln_out((brkr+1):lw) + 180*sign(term_mat(k,2));
   end
  end
 end

end

if mode == 1,

 w_out = wnew;

else

% re-adjust phase vector
 ph = phasecor(exp(i*ln_out*pi/180),[-360,0]);
 brkph=find(abs(diff(ph))>170);
 t=1; ln_out=[]; w_outt=[];
 for k=brkph,
  ln_out=[ln_out,ph(t:k),NaN];
  w_out=[w_out,wnew(t:k),NaN];
  t=k+1;
 end

 ln_out=[ln_out,ph(t:length(ph))];
 w_out=[w_out,wnew(t:length(wnew))];
end

author: Sable

repo: mcbench-benchmarks


text:
function varargout = cumsum2(varargin)
%CUMSUM2   Double indefinite integral of a CHEBFUN2.
%   F = CUMSUM2(F) returns the double indefinite integral of a CHEBFUN2. That is
%                   y  x
%                  /  /
%   CUMSUM2(F) =  |  |   f(x,y) dx dy   for  (x,y) in [a,b] x [c,d],
%                 /  /
%                c  a
%
%   where [a,b] x [c,d] is the domain of f.
%
% See also CUMSUM, SUM, SUM2.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

[varargout{1:nargout}] = cumsum2@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [node,face,elem]=meshunitsphere(tsize,maxvol)
%
% [node,face,elem]=meshunitsphere(tsize,maxvol)
%
% create the surface and/or volumetric mesh of a unit sphere
% centered at [0 0 0] and radius 1
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%   tsize: maximum size of the surface triangles (from 0 to 1)
%   maxvol: maximum volume of the tetrahedron; if one wants to return
%           elem without specifying maxvol, maxvol=tsize^3
%
% output:
%   node: node coordinates, 3 columns for x, y and z respectively
%   face: integer array with dimensions of NB x 3, each row represents
%         a surface mesh face element
%   elem: integer array with dimensions of NE x 4, each row represents
%         a tetrahedron. If ignored, this function only produces the surface
%
% example:
%   [node,face]=meshunitsphere(0.05);
%   [node,face,elem]=meshunitsphere(0.05,0.01);
%   plotmesh(node,elem,'x>0'); axis equal;
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

dim=60;
esize=tsize*dim;
thresh=dim/2-1;
[xi,yi,zi]=meshgrid(0:0.5:dim,0:0.5:dim,0:0.5:dim);
dist=thresh-sqrt((xi-30).^2+(yi-30).^2+(zi-30).^2);
dist(find(dist<0))=0;
clear xi yi zi;

% extract a level-set at v=thresh, being a sphere with R=thresh
% the maximum element size of the surface triangles is tsize*dim

[node,face]=vol2restrictedtri(dist,1,[dim dim dim],dim*dim*dim,30,esize,esize,40000);
node=(node-0.5)*0.5;

[node,face]=removeisolatednode(node,face);

node=(node-30)/28;
r0=sqrt(sum((node.*node)'));
node=node.*repmat(1./r0(:),1,3);

if(nargout==3)
   if(nargin==1) maxvol=tsize*tsize*tsize; end
   [node,elem,face]=surf2mesh(node,face,[-1 -1 -1]*1.1,[1 1 1]*1.1,1,maxvol,[],[]);
   elem=elem(:,1:4);
end
face=face(:,1:3);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ selvalid ] = softSelect( hypScores, thres, soft )
%SOFTSELECT Summary of this function goes here
%   Detailed explanation goes here
nm = (thres-hypScores)./soft;
rn = rand(length(nm),1);
selvalid = rn>nm;

end

author: zouchuhang

repo: LayoutNet


text:
function allPosVar2D = rmPositionVariance2D(view, voxel, samp)
%
% rmPositionVariance2D - calculate the visual area size covered by the pRF
% centers of surrounding nodes, normalized by cortical area size (per 1 mm^2)
%
%  allPosVar2D = rmPositionVariance2D(view, voxel, samp)
%
% INPUT
%  view: VOLUME view should be provided
%  voxel: voxel size in mm (default = 1.5)
%  samp: sampling of neurons within a voxel (default=1)
%        1: uniform sampling, 2: gaussian sampling
% OUTPUT
%  allPosVar2D: 2D position variance (sigma)
%
% KA wrote it 08/10

if ieNotDefined('view')
    view = getCurView;
end
if ieNotDefined('voxel')
    voxel = 1.5;
end
if ieNotDefined('samp')
    samp = 1;
end

nodes   = double(view.nodes);
edges   = double(view.edges);
numNeighbors = double(view.nodes(4,:));
edgeOffsets  = double(view.nodes(5,:));

allPosVar2D=zeros(1,size(nodes,2));

X = view.rm.retinotopyModels{1}.x0;
Y = view.rm.retinotopyModels{1}.y0;

for ii=1:size(nodes,2)
    tri=[];
    neighbors_pRF_area = 0;
    neighbors_cort_area = 0;
    X_tmp = X(ii);
    Y_tmp = Y(ii);
    tmp=1;
    if nodes(6,ii)==1   % use the nodes in layer 1
        neighbors = edges(:, edgeOffsets(ii):edgeOffsets(ii)+numNeighbors(ii)-1);
        % calculate neighbors only in layer 1
        remove_nodes=find(nodes(6,neighbors)~=1);
        neighbors(remove_nodes)=[];

        if size(neighbors,2)>=3
            % Calculate the visual field size covered by pRF centers of neighboring
            % nodes
            for jj=1:size(neighbors,2) % remove neighbors having the same pRF center
                if sum(X_tmp==X(neighbors(jj))|Y_tmp==Y(neighbors(jj)))==0
                    X_tmp = [X_tmp X(neighbors(jj))];
                    Y_tmp = [Y_tmp Y(neighbors(jj))];
                end
            end
            if size(X_tmp,2)>=3 % remove the case where all points are on a single line
                for kk=3:size(X_tmp,2)
                    tmp=tmp & (X_tmp(kk)-X_tmp(2))/(X_tmp(1)-X_tmp(2))*(Y_tmp(1)-Y_tmp(2))+Y_tmp(2)==Y_tmp(kk);
                end
                if tmp~=1
                    % get triangles in visual field to calculate visual
                    % field coverage
                    tri_visual = delaunay(X_tmp,Y_tmp);
                    % calculate the sum of the area size of all triangles
                    for num=1:size(tri_visual,1)
                        a = norm([X_tmp(tri_visual(num,1))-X_tmp(tri_visual(num,2)) Y_tmp(tri_visual(num,1))-Y_tmp(tri_visual(num,2))]);
                        b = norm([X_tmp(tri_visual(num,2))-X_tmp(tri_visual(num,3)) Y_tmp(tri_visual(num,2))-Y_tmp(tri_visual(num,3))]);
                        c = norm([X_tmp(tri_visual(num,1))-X_tmp(tri_visual(num,3)) Y_tmp(tri_visual(num,1))-Y_tmp(tri_visual(num,3))]);
                        s = (a+b+c)/2;

                        neighbors_pRF_area = neighbors_pRF_area + sqrt(s.*(s-a).*(s-b).*(s-c));
                    end
                end
            end

            % Calculate the cortical area size covered by neighboring nodes

            % get triangles in cortex to calculate cortical area size
            for jj=1:size(neighbors,2)
                neighbors_tmp = edges(:,edgeOffsets(neighbors(jj)):edgeOffsets(neighbors(jj))+numNeighbors(neighbors(jj))-1);
                remove_nodes_tmp = find(nodes(6,neighbors_tmp)~=1);
                neighbors_tmp(remove_nodes_tmp)=[];
                for kk=1:size(neighbors_tmp,2)
                    if size(find(neighbors_tmp(kk)==neighbors),2)~=0
                        tmp = sort([ii neighbors(jj) neighbors_tmp(kk)]);
    %                     if size(find(sum((tri-ones(size(tri,1),1)*sort([ii neighbors(jj) neighbors_tmp(kk)]))')'==0),1)
                        if size(tri,1)==0 | (size(tri,1)>0 & ~(sum(tri(:,1)==tmp(1))*sum(tri(:,2)==tmp(2))*sum(tri(:,3)==tmp(3))))
                            tri=[tri;tmp];
                        end
                    end
                end
            end
            % calculate the sum of the cortical area size
            for num=1:size(tri,1)
                a = norm([nodes(1,tri(num,1))-nodes(1,tri(num,2)) nodes(2,tri(num,1))-nodes(2,tri(num,2)) nodes(3,tri(num,1))-nodes(3,tri(num,2))]);
                b = norm([nodes(1,tri(num,2))-nodes(1,tri(num,3)) nodes(2,tri(num,2))-nodes(2,tri(num,3)) nodes(3,tri(num,2))-nodes(3,tri(num,3))]);
                c = norm([nodes(1,tri(num,3))-nodes(1,tri(num,1)) nodes(2,tri(num,3))-nodes(2,tri(num,1)) nodes(3,tri(num,3))-nodes(3,tri(num,1))]);
                s = (a+b+c)/2;

                neighbors_cort_area = neighbors_cort_area + sqrt(s.*(s-a).*(s-b).*(s-c));
            end
        end
        allPosVar2D(ii) = neighbors_pRF_area/neighbors_cort_area;
    end
end

% multiply voxel size and transform into standard deviation
% uniform sampling in voxel
if samp==1
    allPosVar2D(ii) = sqrt(allPosVar2D(ii))*voxel/sqrt(3);
% gaussian sampling in voxel (diameter of the scatter area is assumed to be fwhm)
elseif dist==2
    allPosVar2D(ii) = sqrt(allPosVar2D(ii))*voxel/sqrt(2*log(2));
end

return;

author: vistalab


text:
function model = make_empty_3D_model()

model.x = [];
model.y = [];
model.z = [];
model.nx = [];
model.ny = [];
model.nz = [];
model.r = [];
model.g = [];
model.b = [];

author: lmb-freiburg

repo: orion

save_path: github-repos/MATLAB/lmb-freiburg-orion


text:
function i4vec_order_type_test ( )

%*****************************************************************************80
%
%% I4VEC_ORDER_TYPE_TEST tests I4VEC_ORDER_TYPE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 4;
  test_num = 6;
%
%  Each ROW of the definition is a COLUMN of the matrix.
%
  x_test = [ ...
    1, 3, 2, 4; ...
    2, 2, 2, 2; ...
    1, 2, 2, 4; ...
    1, 2, 3, 4; ...
    4, 4, 3, 1; ...
    9, 7, 3, 0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_ORDER_TYPE_TEST\n' );
  fprintf ( 1, '  I4VEC_ORDER_TYPE classifies an integer vector as\n' );
  fprintf ( 1, '  -1: no order\n' );
  fprintf ( 1, '   0: all equal;\n' );
  fprintf ( 1, '   1: ascending;\n' );
  fprintf ( 1, '   2: strictly ascending;\n' );
  fprintf ( 1, '   3: descending;\n' );
  fprintf ( 1, '   4: strictly descending.\n' );

  for test = 1 : test_num

    x(1:n) = x_test(1:n,test);

    order = i4vec_order_type ( n, x );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  The following vector has order type %d\n', order );
    fprintf ( 1, '\n' );
    for j = 1 : n
      fprintf ( 1, '  %6d  %6d\n', j, x(j) );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function a = c8mat_nint ( m, n, a )

%*****************************************************************************80
%
%% C8MAT_NINT rounds the entries of a C8MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the order of the array.
%
%    Input, complex A(M,N), the array to be rounded.
%
%    Output, complex A(M,N), the rounded array.
%
  for j = 1 : n
    for i = 1 : m
      a(i,j) = c8_nint ( a(i,j) );
    end
  end

  return
end


text:
function [coef]=comp_dwiltii(coef2,a)
%COMP_DWILT  Compute Discrete Wilson transform.
%

M=size(coef2,1)/2;
N=size(coef2,2);
W=size(coef2,3);
L=N*a;

coef=zeros(2*M,N/2,W,assert_classname(coef2));

% ---- m is zero ---------
coef(1,:,:)=coef2(1,1:2:N,:);

% --- m is odd ----------
coef(2:2:M,:,:)    = i/sqrt(2)*(coef2(2:2:M,1:2:N,:)+coef2(2*M:-2:M+2,1:2:N,:));
coef(M+2:2:2*M,:,:)= 1/sqrt(2)*(coef2(2:2:M,2:2:N,:)-coef2(2*M:-2:M+2,2:2:N,:));

% --- m is even ---------
coef(3:2:M,:,:)=     1/sqrt(2)*(coef2(3:2:M,1:2:N,:)-coef2(2*M-1:-2:M+2,1:2:N,:));
coef(M+3:2:2*M,:,:)= i/sqrt(2)*(coef2(3:2:M,2:2:N,:)+coef2(2*M-1:-2:M+2,2:2:N,:));

% --- m is nyquest ------
if mod(M,2)==0
  coef(M+1,:,:) = i*coef2(M+1,2:2:N,:);
else
  coef(M+1,:,:) = i*coef2(M+1,1:2:N,:);
end;

coef=reshape(coef,M*N,W);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function fem2d_pack_test21 ( )

%*****************************************************************************80
%
%% TEST21 tests SPHERE_GRID_Q16.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  element_order = 16;
  nelemx = 2;
  nelemy = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST21\n' );
  fprintf ( 1, '  SPHERE_GRID_Q16_ELEMENT sets up a grid of\n' );
  fprintf ( 1, '    Q16 quadrilaterals on a sphere.\n' );
  fprintf ( 1, '  SPHERE_GRID_Q16_ELEMENT_NUM returns the number\n' );
  fprintf ( 1, '    of elements in the grid\n' );
  fprintf ( 1, '  SPHERE_GRID_Q16_NODE_NUM returns the number\n' );
  fprintf ( 1, '    of nodes in the grid.\n' );
  fprintf ( 1, '  SPHERE_GRID_Q16_NODE_XYZ returns the coordinates\n' );
  fprintf ( 1, '    of nodes in the grid.\n' );

  element_num = sphere_grid_q16_element_num ( nelemx, nelemy );
  node_num = sphere_grid_q16_node_num ( nelemx, nelemy );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Expected number of nodes =    %d\n', node_num );
  fprintf ( 1, '  Expected number of elements = %d\n', element_num );

  element_node = sphere_grid_q16_element ( nelemx, nelemy );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The elements and their nodes, listed in a way\n' );
  fprintf ( 1, '  that suggests their geometry:\n' );
  fprintf ( 1, '\n' );

  element = element_num;

  for j = 1 : nelemy
    for i = 1 : nelemx
      fprintf ( 1, '\n' );
      fprintf ( 1, '%4d  %4d%4d%4d%4d\n', element, element_node(13:16,element) );
      fprintf ( 1, '      %4d%4d%4d%4d\n',         element_node(9:12,element) );
      fprintf ( 1, '      %4d%4d%4d%4d\n',         element_node(5:8,element) );
      fprintf ( 1, '      %4d%4d%4d%4d\n',         element_node(1:4,element) );
      element = element - 1;
    end
  end

  node_xyz = sphere_grid_q16_node_xyz ( nelemx, nelemy );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The node coordinates:\n' );
  fprintf ( 1, '\n' );

  for node = 1 : node_num
    fprintf ( 1, '  %4d  %12f  %12f  %12f\n', node, node_xyz(1:3,node) );
  end
%
%  Write the elements and nodes to files.
%
  r8mat_write ( 'sphere_q16_nodes.txt', 3, node_num, node_xyz );

  i4mat_write ( 'sphere_q16_elements.txt', element_order, element_num, ...
    element_node );

  return
end


text:
function [offset, Nx, Ny] = RegdoBlockMatch(FixedImage, RegisterImage, Nx, Ny)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

global planC
global stateS

indexS = planC{end};

dim = size(FixedImage);

subFixed = cell(Nx, Ny); subMoving = cell(Nx, Ny);
Lx = ceil(dim(2)/Nx); Ly = ceil(dim(1)/Ny);

offset = cell(Nx, Ny); subOrigin = cell(Nx, Ny);
out = cell(Nx, Ny);

for i = 0 : Nx - 1
    for j = 0 : Ny - 1

        x1 = i*Lx + 1; x2 = min( dim(2), (i+1)*Lx );
        y1 = j*Ly + 1; y2 = min( dim(1), (j+1)*Ly );

        subFixed{j+1, i+1} = FixedImage(y1:y2, x1:x2);
        subMoving{j+1, i+1} = RegisterImage(y1:y2, x1:x2);

        %MImg = planC{indexS.scan}(stateS.imageRegistrationBaseDataset).scanArray;
        centerValue = RegisterImage(y1+ceil(abs(y2-y1)/2), x1+ceil(abs(x2-x1)/2));
        [n,x] = hist(double(RegisterImage(:)), 40);
        ind = find(x>=centerValue, 1, 'first');
        if ind>3

            %             switch get(handles.BMMetric,'Value')
            %                 case 1
            %
            %                     [out{j+1, i+1}, rotation, os] = MeanSquare2D(int16(subFixed{j+1, i+1}), [0 0 0], [1 1 1], ...
            %                                                                int16(subMoving{j+1, i+1}), [0 0 0], [1 1 1], ...
            %                                                                0.01, 4, 200);
            %                 case 2
            try
                [out{j+1, i+1}, rotation, os] = NormalizedCorrelation2D(int16(subFixed{j+1, i+1}), [0 0 0], [1 1 1], ...
                    int16(subMoving{j+1, i+1}), [0 0 0], [1 1 1], ...
                    0.01, 4, 200);
            catch
                [out{j+1, i+1}, rotation, os] = MeanSquare2D_64(int16(subFixed{j+1, i+1}), [0 0 0], [1 1 1], ...
                    int16(subMoving{j+1, i+1}), [0 0 0], [1 1 1], ...
                    0.01, 4, 200);
            end
            %             end

            os(2) = -os(2);
            offset{j+1, i+1} = os;
            subOrigin{j+1, i+1} = [(x2+x1)/2 (y2+y1)/2];
        else
            offset{j+1, i+1} = [0 0];
        end
    end

end

end

author: cerr


text:
function [tfg, tfg_timestamp]=extract_tfg(peak,capbuf,fc,sampling_carrier_twist, nRB)

% add 100RB support. Jiao Xianjun(putaohsu@gmail.com)
% Convert from time domain to frequency domain and create the time/frequency
% grid.

% Copyright 2012 Evrytania LLC (http://www.evrytania.com)
%
% Written by James Peroulas <james@evrytania.com>
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU Affero 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 Affero General Public License for more details.
%
% You should have received a copy of the GNU Affero General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

frame_start=peak.frame_start;
cp_type=peak.cp_type;
freq_fine=peak.freq_fine;

% nRB = 6;
if nRB == 6
    conv_idx_ratio = 1;
    decimation_ratio = 16;
elseif nRB == 100
    conv_idx_ratio = 16;
    decimation_ratio = 1;
else
    disp('nRB must be 6 or 100!');
    return;
end

fs = fs_lte/decimation_ratio;
fft_size = 2048/decimation_ratio;
len_cp_extended = 512/decimation_ratio;
len_cp_normal1 = 144/decimation_ratio;
len_cp_normal2 = 160/decimation_ratio;
nSC = nRB*12;

% fc*k_factor is the receiver's actual RX center frequency.
if sampling_carrier_twist==1
    k_factor=(fc-freq_fine)/fc;
else
    k_factor = peak.k_factor;
end

if (strcmpi(cp_type,'normal'))
  n_symb_dl=7;
  dft_location = conv_idx(frame_start+(len_cp_normal2/conv_idx_ratio), conv_idx_ratio);
elseif (strcmpi(cp_type,'extended'))
  n_symb_dl=6;
%   dft_location=frame_start+16; % wrong?
  dft_location = conv_idx(frame_start+(len_cp_extended/conv_idx_ratio), conv_idx_ratio);
else
 error('Check code...');
end

% See if we can advance the frame start
if (dft_location-k_factor*fs*.01>=0.5)
  dft_location=dft_location-k_factor*fs*.01;
end

% Perform FOC
capbuf=fshift(capbuf,-freq_fine,fs);

% Extract 6 frames + 2 slots worth of data
n_ofdm_sym=6*10*2*n_symb_dl+2*n_symb_dl;
tfg=NaN(n_ofdm_sym,fft_size);
tfg_timestamp=NaN(1,n_ofdm_sym);
sym_num=0;
for t=1:n_ofdm_sym
  indices=round(dft_location):round(dft_location)+fft_size-1;
  % Perform 2 sample coarse TOC
  %indices=[indices(3:end) indices(1:2)];
  %indices=[indices(end-2:end) indices(1:end-3)];
  tfg(t,:)=dft(capbuf(indices));
  tfg_timestamp(t)=dft_location;
  if (n_symb_dl==6)
%     dft_location=dft_location+k_factor*(128+16); % wrong?
    dft_location=dft_location+k_factor*(fft_size+len_cp_extended);
  else
    if (sym_num==6)
      dft_location=dft_location+k_factor*(fft_size+len_cp_normal2);
    else
      dft_location=dft_location+k_factor*(fft_size+len_cp_normal1);
    end
    sym_num=mod(sym_num+1,7);
  end
end

% Extract the columns of interest
tfg=[tfg(:,end-((nSC/2) - 1):end) tfg(:,2:((nSC/2) + 1))];

% Compensate for the residual time offset.
cn=[-(nSC/2):-1 1:(nSC/2)];
for t=1:n_ofdm_sym
  ideal_offset=tfg_timestamp(t);
  actual_offset=round(ideal_offset);
  % How late were we in locating the DFT
  late=actual_offset-ideal_offset;
  % Compensate for the improper location of the DFT
  tfg(t,:)=tfg(t,:).*exp(-1i*2*pi*cn*late/fft_size);
end
% At this point, tfg(t,:) contains the results of a DFT that was performed
% at time offset tfg_timestamp(t). Note that tfg_timestamp(t) is not an
% integer!

function out_idx = conv_idx(in_idx, decimation_ratio)
out_idx = (in_idx-1)*decimation_ratio + 1;

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
function [  ] = bz_piTickLabel(axis,discretization)
%bz_piTickLabel(axis) changes tick labels to 0 pi, 2pi, etc.
%
%%
if ~exist('discretization','var')
    discretization = 1;
end

switch axis
    case 'x'
        limits = xlim(gca)./pi;
    case 'y'
        limits = ylim(gca)./pi;
end
whichtick = [axis,'tick'];
whichticklabel = [axis,'ticklabel'];

ticks = round(limits(1):discretization:limits(2),1);
for ll = 1:length(ticks)
    labels{ll} = [num2str(ticks(ll)),'\pi'];
end

set(gca,whichtick,ticks.*pi)
set(gca,whichticklabel,labels)

end

author: buzsakilab


text:
%Database class (genDb)
function [db] = gen_class
   %enteries = enteries_struct();
   db= struct('name',[],'path',[],'format',[],'metaNames',[],'enteries',[]);
   db = class(db, 'gendb');
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function options = scipset(varargin)
%SCIPSET  Create or alter the options for Optimization with SCIP
%
% options = scipset('param1',value1,'param2',value2,...) creates an SCIP
% options structure with the parameters 'param' set to their corresponding
% values in 'value'. Parameters not specified will be set to the SCIP
% default.
%
% options = scipset(oldopts,'param1',value1,...) creates a copy of the old
% options 'oldopts' and then fills in (or writes over) the parameters
% specified by 'param' and 'value'.
%
% options = scipset() creates an options structure with all fields set to
% SCIPSET defaults.
%
% scipset() prints a list of all possible fields and their function.
%
% See supplied SCIP Documentation for further details of these options.

%   Copyright (C) 2013 Jonathan Currie (I2C2)

% Print out possible values of properties.
if ((nargin == 0) && (nargout == 0))
    printfields();
    return
end
%Names and Defaults
Names = {'gamsfile','testmode'}';
Defaults = {[],0}';

%Enter and check user args
try
    options = opticheckset(Names,Defaults,@checkfield,varargin{:});
catch ME
    throw(ME);
end

function checkfield(field,value)
%Check a field contains correct data type
switch lower(field)
    %Scalar 0/1
    case 'testmode'
        err = opticheckval.checkScalar01(value,field);
    %char array
    case 'gamsfile'
        err = opticheckval.checkChar(value,field);
    otherwise
        err = MException('OPTI:SetFieldError','Unrecognized parameter name ''%s''.', field);
end
if(~isempty(err)), throw(err); end

function printfields()
%Print out fields with defaults
fprintf('         gamsfile: [ Write SCIP model to GAMS file (will skip solving): {[]}, ''filename'' ] \n');
fprintf('         testmode: [ Validate the nonlinear function generation (will skip solving): {0}, 1 ] \n');
fprintf('\n');

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus