text:
function gen_caffe_dataset_multilabel_m128(mode, batch_size)

if nargin == 0
    mode = 'train';
    batch_size = 128*1;
end

conf = config;

%% Load images
if strcmp(mode, 'train')
    images = load(fullfile(conf.root_path, ...
        'training_images_crop15_square256.mat'), 'training_images');
    images = images.training_images;
    assert(length(images)==59551);
elseif strcmp(mode, 'val')
    images = load(fullfile(conf.root_path, ...
        'validation_images_crop15_square256.mat'), 'validation_images');
    images = images.validation_images;
else
    error('unknown mode %s', mode);
end

%% Pick image pairs
assert(mod(batch_size, 2)==0);
[image_id_pairs, labels] = get_training_examples_multilabel(mode, batch_size);

%% Prep dataset for C and Caffe conventions
images = cat(4, images.img);

assert(max(image_id_pairs(:)) <= length(images));
% images must be (height x width x channels x num)
assert(ndims(images) == 4);
% assume square images
assert(size(images, 1) == size(images, 2));
% images must have 3 channels
assert(size(images, 3) == 3);

% Convert to BGR
images = images(:,:,[3 2 1],:);
% Switch width and height
images = permute(images, [2 1 3 4]);

% Caffe likes 0-1 labels
labels = int32(labels);
labels(labels == -1) = 0;

labels_cont = 1:11318;
assert(isequal(double(unique(labels)), labels_cont'));

% Subtract 1 for C indexing
image_id_pairs = int32(image_id_pairs - 1);

% check the number of paired data is divisible

%% Do a morph before the mex file

image_id_pairs_serial = int32(zeros(numel(image_id_pairs), 1));
labels_serial = int32(zeros(length(labels)*2, 1));

insert_idx = 1;
lookup_idx = 1;
for i = 1:length(labels)
    if mod(i-1, (batch_size/2)) == 0
        image_id_pairs_serial(insert_idx : insert_idx + batch_size-1, :) = ...
            reshape(image_id_pairs(lookup_idx : lookup_idx + batch_size/2-1,:), [], 1);

        labels_serial(insert_idx : insert_idx + batch_size-1) = ...
            reshape(labels(lookup_idx : lookup_idx + batch_size/2 -1, :), [], 1);

        insert_idx = insert_idx + batch_size;
        lookup_idx = lookup_idx + batch_size/2;
    end
end
assert(length(image_id_pairs_serial) == 2*length(labels));
assert(mod(length(image_id_pairs_serial), batch_size)==0);

%% Write to lmdb
fprintf('Writing level db..\n');
if strcmp(mode, 'train')
    filename = conf.training_set_path_multilabel_m128;
elseif strcmp(mode, 'val')
    filename = conf.validation_set_path_multilabel;
else
    error('unknown mode %s', mode);
end

serialized_pairs_to_leveldb(images, image_id_pairs_serial, labels_serial, filename);
fprintf('Done writing level db..\n');

author: rksltnl


text:
function tf = par_haveupdate(current_file,reference_file)
% Return true if a code update is available.
% Result = par_haveupdate(CurrentFile,ReferenceFile);
%
% This is only for expert use (when workers need to update themselves after a code change on the
% master).
%
% In:
%   CurrentFile : name of the file that is currently executing
%
%   ReferenceFile : name of a replacement file that is possibly newer than the CurrentFile
%
% Out:
%   Result : whether a newer ReferenceFile is available
%
% See also:
%   par_worker
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-08-26

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

tf = false;

try
    current_file = env_translatepath(current_file);
    reference_file = env_translatepath(reference_file);
catch
end

if ~exist(current_file,'file')
    fprintf('The currently executing code (%s) is non-existent; cannot check for updates.\n',current_file);
elseif exist(reference_file,'file')
        % reference file present: could potentially update: compare file dates
        ref_info = dir(reference_file);
        cur_info = dir(current_file);
        if isempty(ref_info)
            fprintf('No file info for ReferenceFile (%s) available. Cannot check for updates.\n',current_file);
            return;
        end
        if isempty(ref_info)
            fprintf('No file info for CurrentFile (%s) available. Cannot check for updates.\n',current_file);
            return;
        end
        % compare file dates
        tf = ref_info.datenum > cur_info.datenum;
    end
end

repo: REST


text:
function [s,t]=LFMChirp(T,fStart,fEnd,tSamp,phi)
%%LFMCHIRP Generate a complex linearly frequency modulated (LFM) signal. It
%          can be an upchirp or a downchirp. LFM chirps have a role in
%          radar and sonar signal processing. Letting fStart=-fEnd will
%          generate a chirp at baseband with bandwidth 2*fEnd.
%
%INPUTS:  T The duration of the chirp (seconds).
%    fStart The frequency at the start of the chirp (Hertz).
%      fEnd The frequency at the end of the chirp (Hertz).
%     tSamp If this is scalar, it is the sampling rate of the chirp and the
%           number of samples depends on T. Otherwise, this is a numSampX1
%           array of times at which samples are desired. Samples before
%           time 0 and after time T are set to zero. If this parameter is
%           omitted or an empty matrix is passed, the Nyquist sampling
%           rate is used.
%       phi The phase offset of the chirp. If this parameter is omitted or
%           an empty matrix is passed, then the default of pi/2 is used,
%           which makes the real part of the signal start at zero.
%
%OUTPUT: s The numSampX1 complex chirp signal.
%        t The sample times used. If tSamp was not a scalar, then this is
%          just tSamp.
%
%LFM chirps are common in radar. They are discussed in Section 8.2 of [1].
%In [2], biases relevant to tracking from LFM signals processed without
%Doppler filtering are discussed. The real component of the complex
%waveform is cos(2*pi*fStart*t+pi*alpha*t.^2+phi), where
%alpha=(fEnd-fStart)/T.
%
%EXAMPLE 1:
%Here is an example of an upchirp:
% T=1e-3;
% fStart=0;
% fEnd=10e3;
% tSamp=1/(100*2*fEnd);%100 times the Nyquist frequency.
% [s,t]=LFMChirp(T,fStart,fEnd,tSamp);
% figure(1)
% clf
% plot(t,real(s),'linewidth',2)
% h1=xlabel('t');
% h2=ylabel('s');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%EXAMPLE 2:
%Here is an example of an downchirp:
% T=1e-3;
% fStart=10e3;
% fEnd=0;
% tSamp=1/(100*2*fStart);%100 times the Nyquist frequency.
% [s,t]=LFMChirp(T,fStart,fEnd,tSamp);
% figure(2)
% clf
% plot(t,real(s),'linewidth',2)
% h1=xlabel('t');
% h2=ylabel('s');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%REFERENCES:
%[1] M. I. Skolnik, Introduction to Radar Systems, 3rd ed. Boston:
%    McGraw Hill, 2001.
%[2] R. J. Fitzgerald, "Effects of range-Doppler coupling on chirp radar
%    tracking accuracy," IEEE Transactions on Aerospace and Electronic
%    Systems, vol. 10, no. 4, pp. 528-532, Jul. 1974.
%
%November 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4||isempty(tSamp))
   %Default to the Nyquist sampling rate.
   tSamp=1/(2*max(fStart,fEnd));
end

if(nargin<5||isempty(phi))
   phi=pi/2;%Phase offset so the signal starts at zero, by default.
end

alpha=(fEnd-fStart)/T;

if(isscalar(tSamp))
    %If a scalar sampling rate is given.
    NPulse=fix(T/tSamp);
    t=(0:(NPulse-1))*tSamp;
else
    %If the samples themselves are given
    t=tSamp;
end

%The real signal corresponding to the complex waveform given below is
%cos(2*pi*fStart*t+pi*alpha*t.^2+phi);
s=exp(1j*(2*pi*fStart*t+pi*alpha*t.^2+phi));

sel=(t>T)|(t<0);
s(sel)=0;

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function tnet = vl_simplenn_tidy(net)
%VL_SIMPLENN_TIDY  Fix an incomplete or outdated SimpleNN network.
%   NET = VL_SIMPLENN_TIDY(NET) takes the NET object and upgrades
%   it to the current version of MatConvNet. This is necessary in
%   order to allow MatConvNet to evolve, while maintaining the NET
%   objects clean. This function ignores custom layers.
%
%   The function is also generally useful to fill in missing default
%   values in NET.
%
%   See also: VL_SIMPLENN().

% Copyright (C) 2014-15 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

tnet = struct('layers', {{}}, 'meta', struct()) ;

% copy meta information in net.meta subfield
if isfield(net, 'meta')
  tnet.meta = net.meta ;
end

if isfield(net, 'classes')
  tnet.meta.classes = net.classes ;
end

if isfield(net, 'normalization')
  tnet.meta.normalization = net.normalization ;
end

% copy layers
for l = 1:numel(net.layers)
  defaults = {'name', sprintf('layer%d', l), 'precious', false};
  layer = net.layers{l} ;

  % Ignore custom layers (e.g. for classes the `isfield` does not work)
  % The only interface requirement for custom layers is forward and
  % backward function.
  if strcmp(layer.type, 'custom')
    tnet.layers{l} = layer ;
    continue;
  end

  % check weights format
  switch layer.type
    case {'conv', 'convt', 'bnorm'}
      if ~isfield(layer, 'weights')
        layer.weights = {...
          layer.filters, ...
          layer.biases} ;
        layer = rmfield(layer, 'filters') ;
        layer = rmfield(layer, 'biases') ;
      end
  end
  if ~isfield(layer, 'weights')
    layer.weights = {} ;
  end

  % Check that weights include moments in batch normalization.
  if strcmp(layer.type, 'bnorm')
    if numel(layer.weights) < 3
      layer.weights{3} = ....
        zeros(numel(layer.weights{1}),2,'single') ;
    end
  end

  % Fill in missing values.
  switch layer.type
    case 'conv'
      defaults = [ defaults {...
        'pad', 0, ...
        'stride', 1, ...
        'dilate', 1, ...
        'opts', {}}] ;

    case 'pool'
      defaults = [ defaults {...
        'pad', 0, ...
        'stride', 1, ...
        'opts', {}}] ;

    case 'convt'
      defaults = [ defaults {...
        'crop', 0, ...
        'upsample', 1, ...
        'numGroups', 1, ...
        'opts', {}}] ;

    case {'pool'}
      defaults = [ defaults {...
        'method', 'max', ...
        'pad', 0, ...
        'stride', 1, ...
        'opts', {}}] ;

    case 'relu'
      defaults = [ defaults {...
        'leak', 0}] ;

    case 'dropout'
      defaults = [ defaults {...
        'rate', 0.5}] ;

    case {'normalize', 'lrn'}
      defaults = [ defaults {...
        'param', [5 1 0.0001/5 0.75]}] ;

    case {'pdist'}
      defaults = [ defaults {...
        'noRoot', false, ...
        'aggregate', false, ...
        'p', 2, ...
        'epsilon', 1e-3, ...
        'instanceWeights', []} ];

    case {'bnorm'}
      defaults = [ defaults {...
        'epsilon', 1e-5 } ] ;
  end

  for i = 1:2:numel(defaults)
    if ~isfield(layer, defaults{i})
      layer.(defaults{i}) = defaults{i+1} ;
    end
  end

  % save back
  tnet.layers{l} = layer ;
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function status = plotMultipleSingleCycleErr(vw,scan,ROIlist)
%Plot the time series data of a single cycle.
%
%  status = plotMultipleSingleCycle(vw,[scan], [ROIlist])
%
% All cycles (blocks) in a scan are collapsed into one cycle. This is meant
% to be somthing like the average over repetitions in the scan. The mean is
% plotted with standard error bars.
%
% When called from the GRAY/VOLUME or FLAT views, the ROIs are converted to
% INPLANE to retrieve the time series.
%
% Example:
%   plotMultipleSingleCycleErr(INPLANE{1},1)
%
%  08/06 KA wrote it

mrGlobals;

% Normal status is OK
status = 1;

if notDefined('vw'),        vw       = getCurView;              end
if notDefined('scan'),      scan    = viewGet(vw,'curScan');    end
if notDefined('ROIlist'),   ROIlist =[];                        end

% Get coranal parameters
nCycles        = viewGet(vw, 'num cycles', scan);
frameRate      = viewGet(vw, 'frame rate', scan);
framesToUse    = viewGet(vw, 'frames to use', scan);
nFrames        =length(framesToUse);
framesPerCycle = nFrames/nCycles;
tScan          = (0:nFrames-1)*frameRate;

% If we don't have an integer number of data points per cycle, then
%   resample the data
if isinteger(framesPerCycle)
    resample = false;
    tCycle = tScan(1:framesPerCycle+1);

else
    resample = true;

    framesPerCycle      = round(framesPerCycle);
    nFrames             = framesPerCycle * nCycles;
    resampledT          = linspace(0,max(tScan), nFrames);
    tCycle              = resampledT(1:framesPerCycle+1);

end

% Select ROIs
[selectedROIs, nROIs] = roiGetList(vw, ROIlist);

% Get coords
ROIcoords = cell(1,nROIs);
for r=1:nROIs
    ROIcoords{r}=viewGet(vw, 'ROIcoords', selectedROIs(r));
end

%% Specifics for Flat, Gray, or Inplane views - xform ROI to INPLANE view
% This whole section is duplicated in other functions, like
% plotMultipleSingleCycleErr.m Maybe it should be split off to its own
% function.

%Compute the average single cycle
for r=1:nROIs
	tSeries                     = meanTSeriesForPlotting(vw, selectedROIs(r));

    if resample
        resampledTSeries        = interp1(tScan, tSeries{1}(framesToUse), resampledT);
        multiCycle              = reshape(resampledTSeries,nFrames/nCycles,nCycles);

    else
        multiCycle              = reshape(tSeries{1}(framesToUse),nFrames/nCycles,nCycles);
    end

    singleCycle{r}              = mean(multiCycle,2);
    singleCycleStdErr{r}        = (std(multiCycle, [], 2)/sqrt(nCycles));
    singleCycle{r}(end+1)       = singleCycle{r}(1);
    singleCycleStdErr{r}(end+1) = singleCycleStdErr{r}(1);

end

framesPerCycle=framesPerCycle+1;

% Plotting section
newGraphWin;
hold on

fontSize = 14;

% ROIname = view.ROIs(view.selectedROI).name;
% headerStr = ['Mean Cycle, ROI ',ROIname,', scan ',num2str(scan)];

% set(gcf,'Name',headerStr);
for r=1:nROIs
    hh = errorbar(tCycle,singleCycle{r},singleCycleStdErr{r});

    % set the line color to be the same as the ROI color
    set(hh,'Color',vw.ROIs(selectedROIs(r)).color);

    % but if the line color and plot color are the same, the line will be
    % invisible (e.g., if the ROI color is white)
    if isequal(get(hh, 'Color'), get(gca, 'Color'))
        set(hh, 'Color', 1 - get(hh, 'Color'));
    end

    set(hh,'LineWidth',4);

    tmp{r} = viewGet(vw, 'roiName', selectedROIs(r));
end
legend(tmp)
hold off;

% nTicks = size(tSeries,1);
xtick = (0:frameRate:framesPerCycle*frameRate);

set(gca,'xtick',xtick)
set(gca,'FontSize',fontSize)
xlabel('Time (sec)','FontSize',fontSize)
ylabel('Percent modulation','FontSize',fontSize)
set(gca,'XLim',[0,framesPerCycle*frameRate]);
grid on

%Save the data in gca('UserData')
data.x = tCycle;
data.y = singleCycle;
data.e = singleCycleStdErr;
set(gca,'UserData',data);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function styles =  plotColors()

colors =  ['r' 'b' 'k' 'g' 'c' 'y' 'm' ...
	   'r' 'b' 'k' 'g' 'c' 'y' 'm'];
symbols = ['o' 'x' '+' '>' '<' '^' 'v' ...
	   '*' 'p' 'h' 's' 'd' 'o' 'x'];
for i=1:length(colors)
  styles{i} = sprintf('-%s%s', colors(i), symbols(i));
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function class_str2node  = get_class2node(hash,wnid)
	class_str2node = hash.get(wnid);
end

author: zhoubolei

repo: CAM

save_path: github-repos/MATLAB/zhoubolei-CAM


text:
function [c,info] = ufwt(f,w,J,varargin)
%UFWT  Undecimated Fast Wavelet Transform
%   Usage:  c = ufwt(f,w,J);
%           [c,info] = ufwt(...);
%
%   Input parameters:
%         f     : Input data.
%         w     : Wavelet Filterbank.
%         J     : Number of filterbank iterations.
%
%   Output parameters:
%         c     : Coefficients stored in $L \times (J+1)$ matrix.
%         info  : Transform paramaters struct.
%
%   `ufwt(f,w,J)` computes redundant time (or shift) invariant
%   wavelet representation of the input signal *f* using wavelet filters
%   defined by *w* in the "a-trous"  algorithm.
%
%   For all accepted formats of the parameter *w* see the |fwtinit| function.
%
%   `[c,info]=ufwt(f,w,J)` additionally returns the `info` struct.
%   containing the transform parameters. It can be conviniently used for
%   the inverse transform |iufwt| e.g. `fhat = iufwt(c,info)`. It is also
%   required by the |plotwavelets| function.
%
%   The coefficents *c* are so called undecimated Discrete Wavelet transform
%   of the input signal *f*, if *w* defines two-channel wavelet filterbank.
%   Other names for this version of the wavelet transform are: the
%   time-invariant wavelet transform, the stationary wavelet transform,
%   maximal overlap discrete wavelet transform or even the "continuous"
%   wavelet transform (as the time step is one sample). However, the
%   function accepts any number filters (referred to as $M$) in the basic
%   wavelet filterbank and the number of columns of *c* is then $J(M-1)+1$.
%
%   For one-dimensional input *f* of length *L*, the coefficients *c* are
%   stored as columns of a matrix. The columns are ordered with inceasing
%   central frequency of the respective subbands.
%
%   If the input *f* is $L \times W$ matrix, the transform is applied
%   to each column and the outputs are stacked along third dimension in the
%   $L \times J(M-1)+1 \times W$ data cube.
%
%   Filter scaling
%   --------------
%
%   When compared to |fwt|, |ufwt| subbands are gradually more and more
%   redundant with increasing level of the subband. If no scaling of the
%   filters is introduced, the energy of subbands tends to grow with increasing
%   level.
%   There are 3 flags defining filter scaling:
%
%      'sqrt'
%               Each filter is scaled by `1/sqrt(a)`, where *a* is the hop
%               factor associated with it. If the original filterbank is
%               orthonormal, the overall undecimated transform is a tight
%               frame.
%               This is the default.
%
%      'noscale'
%               Uses filters without scaling.
%
%      'scale'
%               Each filter is scaled by `1/a`.
%
%   If 'noscale' is used, 'scale' has to be used in |iufwt| (and vice
%   versa) in order to obtain a perfect reconstruction.
%
%   Boundary handling:
%   ------------------
%
%   `c=ufwt(f,w,J)` uses periodic boundary extension. The extensions are
%   done internally at each level of the transform, rather than doing the
%   prior explicit padding.
%
%   Examples:
%   ---------
%
%   A simple example of calling the |ufwt| function:::
%
%     [f,fs] = greasy;
%     J = 8;
%     [c,info] = ufwt(f,'db8',J);
%     plotwavelets(c,info,fs,'dynrange',90);
%
%   See also: iufwt, plotwavelets
%
%   References: holschneider1989real

% AUTHOR: Zdenek Prusa

complainif_notenoughargs(nargin,3,'UFWT');
complainif_notposint(J,'J');

definput.import = {'uwfbtcommon'};
[flags]=ltfatarghelper({},definput,varargin);

% Initialize the wavelet filters structure
w = fwtinit(w);

%% ----- step 1 : Verify f and determine its length -------
% Change f to correct shape.
[f,Ls]=comp_sigreshape_pre(f,upper(mfilename),0);
if(Ls<2)
   error('%s: Input signal seems not to be a vector of length > 1.',upper(mfilename));
end

%% ----- step 2 : Run computation
c = comp_ufwt(f,w.h,w.a,J,flags.scaling);

%% ----- Optionally : Fill info struct ----
if nargout>1
   info.fname = 'ufwt';
   info.wt = w;
   info.J = J;
   info.scaling = flags.scaling;
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function el = export_graphic_annotation_sequence(args)
%"export_graphic_annotation_sequence"
%   Subfunction to handle graphic annotation sequence.
%   Uses the same layout and principle as the parent function.
%
%   This function takes a CERR GSPS element and the index of that
%   element within the planC.GSPS array.
%
%JRA 06/23/06
%NAV 07/19/16 updated to dcm4che3
%   replaced ml2dcm_Element to data2dcmElement
%
%Usage:
%   @export_graphic_annotation_sequence(args)
%
% 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/>.

%Init output element to empty.
el = [];

%Unpack input data.
tag         = args.tag;
gspsS       = args.data{1};
scanInfoS   = args.data{2};
template    = args.template;

switch tag
    case 528704  %0008,1140  Referenced Image Sequence

        templateEl  = template.getValue(tag);
        fHandle = @export_referenced_image_sequence;

        %New null sequence
        tmp = javaObject('org.dcm4che3.data.Attributes');
        el = tmp.newSequence(tag, 0);

        i = 1; % only one layer
        dcmobj = export_sequence(fHandle, templateEl, {scanInfoS});
        el.add(i-1, dcmobj);

        %get attribute to return
        el = el.getParent();

    case 7340034  %0070,0002  Graphic layer

        graphicLayer = gspsS.presentLabel; % must be the same name used inthe graphic layer sequence
        el = data2dcmElement(graphicLayer, tag);

    case 7340040  %0070,0008  Text object sequence

        templateEl  = template.getValue(tag);
        fHandle = @export_text_object_sequence;

        %New null sequence
        tmp = javaObject('org.dcm4che3.data.Attributes');
        el = tmp.newSequence(tag, 0);

        textAnnotationS = gspsS.textAnnotationS;
        nAnnot = length(textAnnotationS);

        for i=1:nAnnot
            dcmobj = export_sequence(fHandle, templateEl, {textAnnotationS(i)});
            %dcmobj = export_sequence(fHandle, tag, {structS(i), i});
            el.add(i-1, dcmobj);
        end
        %get attribute to return
        el = el.getParent();

    case 7340041  %0070,0009  Graphic object sequence

        templateEl  = template.getValue(tag);
        fHandle = @export_graphic_object_sequence;

        %New null sequence
        tmp = javaObject('org.dcm4che3.data.Attributes');
        el = tmp.newSequence(tag, 0);

        graphicAnnotationS = gspsS.graphicAnnotationS;
        nAnnot = length(graphicAnnotationS);

        for i=1:nAnnot
            dcmobj = export_sequence(fHandle, templateEl, {graphicAnnotationS(i)});
            %dcmobj = export_sequence(fHandle, tag, {structS(i), i});
            el.add(i-1, dcmobj);
        end
        %get attribute to return
        el = el.getParent();

    otherwise
        warning(['No methods exist to populate DICOM structure_set module''s structure_set_ROI sequence field: ' dec2hex(tag,8) '.']);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% plot a 2D SLAM graph
function plot_graph(g, iteration = -1)

clf;
graphics_toolkit gnuplot
hold on;

[p, l] = get_poses_landmarks(g);

if (length(l) > 0)
  landmarkIdxX = l+1;
  landmarkIdxY = l+2;
  plot(g.x(landmarkIdxX), g.x(landmarkIdxY), '.or', 'markersize', 4);
end

if (length(p) > 0)
  pIdxX = p+1;
  pIdxY = p+2;
  plot(g.x(pIdxX), g.x(pIdxY), '.xb', 'markersize', 4);
end

% draw line segments???
if 0
  poseEdgesP1 = [];
  poseEdgesP2 = [];
  landmarkEdgesP1 = [];
  landmarkEdgesP2 = [];
  for eid = 1:length(g.edges)
    edge = g.edges(eid);
    if (strcmp(edge.type, 'P') != 0)
      poseEdgesP1 = [poseEdgesP1, g.x(edge.fromIdx:edge.fromIdx+1)];
      poseEdgesP2 = [poseEdgesP2, g.x(edge.toIdx:edge.toIdx+1)];
    elseif (strcmp(edge.type, 'L') != 0)
      landmarkEdgesP1 = [landmarkEdgesP1, g.x(edge.fromIdx:edge.fromIdx+1)];
      landmarkEdgesP2 = [landmarkEdgesP2, g.x(edge.toIdx:edge.toIdx+1)];
    end
  end

  linespointx = [poseEdgesP1(1,:); poseEdgesP2(1,:)];
  linespointy = [poseEdgesP1(2,:); poseEdgesP2(2,:)];

  plot(linespointx, linespointy, "r");
end

%plot(poseEdgesP1(1,:), poseEdgesP1(2,:), "r");

%if (columns(poseEdgesP1) > 0)
%end
%if (columns(landmarkEdges) > 0)
%end

hold off;

figure(1, "visible", "off");
drawnow;
%pause(0.1);
if (iteration >= 0)
  filename = sprintf('../plots/lsslam_%03d.png', iteration);
  print(filename, '-dpng');
end

end

author: kiran-mohan

repo: SLAM-Algorithms-Octave


text:
function outsig = auralize_ir(ir,content,usenorm,conf)
%AURALIZE_IR auralizes an impulse response with an audio file/signal
%
%   Usage: outsig = auralize_ir(ir,content,[normalize],conf)
%
%   Input parameters:
%       ir        - impulse response (IR). Also an binaural room scanning
%                   (BRS) matrix can be auralized, then only the two first
%                   channels will be used.
%       content   - content file or signal vector to be used for auralisation
%                   (mono, if it contains more than one channel, only the
%                   first will be used).
%       normalize - normalize the signal (1 or 0), default: 1
%       conf      - configuration struct (see SFS_config)
%
%   AURALIZE_IR(ir,content,normalize,conf) convolves the first two channels of
%   the given IR with the given content and returns the resulting outsig. The
%   content can be specified in the form of an explicite content file or vector.
%
%   See also: ir_wfs, ir_generic, ir_point_source

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

%% ===== Checking of input parameters and configuration =================
nargmin = 3;
nargmax = 4;
narginchk(nargmin,nargmax);
isargmatrix(ir);
if nargin<nargmax
    conf = usenorm;
    usenorm = 1;
end
isargscalar(usenorm);
isargstruct(conf);

%% ===== Configuration ==================================================
% Sampling rate
fs = conf.fs;

%% ===== Get the right content ==========================================
if isnumeric(content)
    contentfs = conf.fs;
elseif exist(content,'file')
    [content,contentfs] = audioread(content);
else
    error('%s: %s file was not found.',upper(mfilename),content);
end

%% ===== Convolution of the IR with content =============================
% Check if we have to resample
if contentfs~=fs
    content = resample(content,fs,contentfs);
end
% Check if the content is only one channel and fix it otherwise
if min(size(content))~=1
    if size(content,1)<size(content,2)
        content = content(1,:);
    else
        content = content(:,1);
    end
end
% Convolve the two
if size(ir,2)>2
    warning('Your impulse response has more than two channels.');
end
% Convolve the impulse responses with the content signal
outsig = convolution(ir,content);
% Scale output
if(usenorm)
    outsig = norm_signal(outsig);
end

author: sfstoolbox

repo: sfs-matlab


text:
function [this, hasRegressors, hasConditions]  = init_glm(this)
% This method initializes MrGlm, i.e. checks for consistency and saves
% neccessary files
%
%   Y = MrGlm()
%   Y.specify_Glm(inputs)
%
% This is a method of class MrGlm.
%
% IN
%
% OUT
%
% EXAMPLE
%   init_glm
%
%   See also MrGlm

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2014-11-07
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.
%

% save regressor file
R = struct2array(this.regressors);
if ~isempty(R)
    R = struct2array(this.regressors);
    fileNameRegressors = fullfile(this.parameters.save.path, 'Regressors');
    save(fileNameRegressors, 'R');
    hasRegressors = 1;
else
    disp('No regressors specified. Are you sure?');
    hasRegressors = 0;
end

% save conditions file
if ~isempty(this.conditions.names)
    names = this.conditions.names;
    onsets = this.conditions.onsets;
    durations = this.conditions.durations;
    fileNameConditions = fullfile(this.parameters.save.path, 'Conditions');
    save(fileNameConditions, 'names', 'onsets', 'durations');
    hasConditions = 1;
else
    disp('No conditions specified. Are you sure?');
    hasConditions = 0;

end

% make SPM directory
spmDirectory = fullfile(this.parameters.save.path, this.parameters.save.spmDirectory);
if ~exist(spmDirectory)
    mkdir(spmDirectory);
end

author: translationalneuromodeling

repo: tapas


text:
function [filename] = getfilename(ds, scnls, starttimes)
   %GETFILENAME returns a filename that is associated with a single SCNL
   %  filename= getfilename(ds, scnlobject, starttime)
   %    where:
   %    DS is the datasource object
   %    SCNLOBJECT is a scnlobject with station/channel/network/location
   %    information.
   %    STARTTIME is a matlab formatted time, or a string
   %
   % getfile is capable of traversing directory trees.
   %
   % [filename] = getfilename(ds, [],[],[]) will return a sample-string.
   % ex.
   %  >> ds = datasource('file','myfiles/%04d/%02d/thisfile%04d_%s%s.dat',...
   %      'year', 'month','day','station','channel');
   %  >> getfile(ds,[],[])
   %  ans =
   %     myfiles/[YEAR]/[MONTH]/thisfile[DAY]_[STATION][CHANNEL].dat
   %
   %  >> getfilename(ds,scnlobject('OKCF','EHZ'),'12/31/2008 04:03:21')
   %  ans =
   %     myfiles/2008/12/thisfile0031_OKCFEHZ.dat
   %
   %
   % see also SETFILENAME

   % AUTHOR: Celso Reyes, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   % Day numbers at start of each month.  Used within to calculate "jday".
   monthadd = [ 0 31 59 90 120 151 181 212 243 273 304 334 ];

   if ~ds.usefile
      warning('Datasource:noAssociatedFiles','No files have been associated with this datasource');
      filename = [];
      return
   end

   %
   % if numel(scnl) > 1
   %   warning('Datasource:multipleSCNL',...
   %     'An array of SCNL was passed as the argument.  Only search for a single SCNL at a time');
   % end
   if isnumeric(starttimes) && ~isempty(starttimes)
      for n=1:numel(starttimes)
         st(n) = {starttimes(n)};
      end
      starttimes = st;
      clear st
   end
   if ~iscell(starttimes)
      starttimes = {starttimes};
   end
   filename = cell(numel(scnls),numel(starttimes));

   if isempty(filename) %show the generic file name
      %return a THEORETICAL filename
      year = '[YEAR]';
      jday ='[JDAY]';
      month ='[MONTH]';
      day ='[DAY]';
      hour ='[HOUR]';
      minute ='[MINUTE]';
      second ='[SECOND]';
      station = '[STATION]';
      channel = '[CHANNEL]';
      location = '[LOCATION]';
      network = '[NETWORK]';
      stringToShow = regexprep(ds.file_string,'(%|%0)(\d|)d','%s');
      % build the file arg string
      if isempty(ds.file_args)
         filename = {stringToShow};
      else
         z = cellfun(@eval,ds.file_args,'UniformOutput',false);
         filename = {sprintf(stringToShow,z{:})};
      end
   else %get detailed file.for n=1:numel(scnls)

      for n=1:numel(scnls)
         scnl = scnls(n);
         station = get(scnl,'station');
         channel = get(scnl,'channel');
         location = get(scnl,'location');
         network = get(scnl,'network');

         for m = 1:numel(starttimes)
            starttime = starttimes{m};

            %pull all the data that may be required to parse out files & directories

            % If we have a datenum, use the fast built-in version of datevec
            if isnumeric(starttime)
               [year, month, day, hour, minute, second] = datevecmx(starttime);
            else
               [year, month, day, hour, minute, second] = datevec(starttime);
            end

            second = round(second);

            % The following few IF statements handle possible rounding issues
            if second == 60
               second = 0; minute = minute+1;
            end
            if minute == 60
               minute = 0; hour = hour + 1;
            end
            if hour == 24
               hour = 0; day = day + 1;
            end

            % Calculate day number
            leapadd = +((mod(year,400) == 0) | ((mod(year,4) == 0) & (mod(year,100) ~= 0)) & (month > 2));
            jday = monthadd(month) + day + leapadd;

            stringToShow = ds.file_string;
            % build the file arg string
            if isempty(ds.file_args)
               filename = {stringToShow};
            else
               z = cellfun(@eval,ds.file_args,'UniformOutput',false);
               filename(n,m) = {sprintf(stringToShow,z{:})};
            end
         end
      end

   end
   filename = filename(:);
   %basic location of all archives
end

author: geoscience-community-codes

repo: GISMO


text:
function area = domainarea( f )
%DOMAINAREA    Area of the domain of f.
%   DOMAINAREA(F) returns the area of the domain of f.

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

if ( isempty( f ) )
    area = 0;
else
    dom = f.domain;
    area = diff( dom(1:2) )*diff( dom(3:4) );
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [fe,ll] = hmmfe(data,T,hmm,Gamma,Xi,preproc,grouping)
% Computes the Free Energy of an HMM
%
% INPUT
% data          observations, either a struct with X (time series) and C (classes, optional)
%                             or just a matrix containing the time series
% T            length of series
% hmm          hmm structure
% Gamma        probability of states conditioned on data (optional)
% Xi           joint probability of past and future states conditioned on data (optional)
%
% OUTPUT
% fe         the variational free energy
% ll         log-likelihood per time point
%
% Author: Diego Vidaurre, OHBA, University of Oxford (2017)

% to fix potential compatibility issues with previous versions
hmm = versCompatibilityFix(hmm);
mixture_model = isfield(hmm.train,'id_mixture') && hmm.train.id_mixture;
p = hmm.train.lowrank; do_HMM_pca = (p > 0);

if nargin<4, Gamma = []; end
if nargin<5, Xi = []; end
if nargin<6 || isempty(preproc), preproc = 1; end
if nargin<7 , grouping = ones(length(T),1); end
if size(grouping,1)==1,  grouping = grouping'; end

if isstruct(data), data = data.X; end

options = hmm.train;
hmm.train.grouping = grouping;

stochastic_learn = isfield(options,'BIGNbatch') && ...
    (options.BIGNbatch < length(T) && options.BIGNbatch > 0);

if iscell(T)
    for i = 1:length(T)
        if size(T{i},1)==1, T{i} = T{i}'; end
    end
    if size(T,1)==1, T = T'; end
end

if ~stochastic_learn
    if iscell(T)
        T = cell2mat(T);
    end
    checkdatacell;
    %data = data2struct(data,T,hmm.train);
else
    if ~iscell(data)
        N = length(T);
        dat = cell(N,1); TT = cell(N,1);
        for i=1:N
            t = 1:T(i);
            dat{i} = data(t,:); TT{i} = T(i);
            try data(t,:) = [];
            catch, error('The dimension of data does not correspond to T');
            end
        end
        if ~isempty(data)
            error('The dimension of data does not correspond to T');
        end
        data = dat; T = TT; clear dat TT
    end
end

if preproc && ~stochastic_learn
    % Standardise data and control for ackward trials
    data = standardisedata(data,T,options.standardise);
    % Filtering
    if ~isempty(options.filter)
        data = filterdata(data,T,options.Fs,options.filter);
    end
    % Detrend data
    if options.detrend
        data = detrenddata(data,T);
    end
    % Leakage correction
    if options.leakagecorr ~= 0
        data = leakcorr(data,T,options.leakagecorr);
    end
    % Hilbert envelope
    if options.onpower
        data = rawsignal2power(data,T);
    end
    % Leading Phase Eigenvectors
    if options.leida
        data = leadingPhEigenvector(data,T);
    end
    % pre-embedded PCA transform
    if length(options.pca_spatial) > 1 || (options.pca_spatial > 0 && options.pca_spatial ~= 1)
        if isfield(options,'As')
            data.X = bsxfun(@minus,data.X,mean(data.X));
            data.X = data.X * options.As;
        else
            [options.As,data.X] = highdim_pca(data.X,T,options.pca_spatial);
            options.pca_spatial = size(options.As,2);
        end
    end
    % Embedding
    if length(options.embeddedlags)>1
        [data,T] = embeddata(data,T,options.embeddedlags);
    end
    % PCA transform
    if isfield(options,'A')
        data = bsxfun(@minus,data,mean(data)); % must center
        data = data * options.A;
        % Standardise principal components and control for ackward trials
        data = standardisedata(data,T,options.standardise_pc);
    end
    % Downsampling
    if options.downsample > 0
        [data,T] = downsampledata(data,T,options.downsample,options.Fs);
    end
end

% get residuals
if ~stochastic_learn && ~do_HMM_pca
    if isfield(hmm.state(1),'W')
        ndim = size(hmm.state(1).W.Mu_W,2);
    else
        ndim = size(hmm.state(1).Omega.Gam_rate,2);
    end
    if ~isfield(hmm.train,'Sind')
        orders = formorders(hmm.train.order,hmm.train.orderoffset,hmm.train.timelag,hmm.train.exptimelag);
        hmm.train.Sind = formindexes(orders,hmm.train.S) == 1;
        if ~hmm.train.zeromean, hmm.train.Sind = [true(1,ndim); hmm.train.Sind]; end
    end
    residuals =  getresiduals(data,T,hmm.train.S,hmm.train.maxorder,hmm.train.order,...
        hmm.train.orderoffset,hmm.train.timelag,hmm.train.exptimelag,hmm.train.zeromean);
else
    residuals = [];
end

% get state time courses
if isempty(Gamma) || isempty(Xi)
   if ~(mixture_model && ~isempty(Gamma)) % we have Gamma and Xi is not needed
       [Gamma,Xi] = hmmdecode(data,T,hmm,0,residuals,0);
   end
end

if stochastic_learn
    if hmm.train.downsample > 0
        downs_ratio = (hmm.train.downsample/hmm.train.Fs);
    else
        downs_ratio = 1;
    end
    Tmat = downs_ratio*cell2mat(T);
    if length(hmm.train.embeddedlags)>1
        maxorder = 0;
        L = -min(hmm.train.embeddedlags) + max(hmm.train.embeddedlags);
        Tmat = Tmat - L;
    else
        maxorder = hmm.train.maxorder;
    end
    %  P/Pi KL
    fe = sum(evalfreeenergy([],[],[],[],hmm,[],[],[0 0 0 1 0]));
    % Gamma entropy&LL
    fe = fe + sum(evalfreeenergy([],Tmat,Gamma,Xi,hmm,[],[],[1 0 1 0 1]));
    tacc = 0; tacc2 = 0; fell = 0; ll = [];
    for i = 1:length(T)
        [X,XX,residuals,Ti] = loadfile(data{i},T{i},options);
        if hmm.train.lowrank > 0,  XX = X; end
        t = (1:(sum(Ti)-length(Ti)*maxorder)) + tacc;
        t2 = (1:(sum(Ti)-length(Ti)*(maxorder+1))) + tacc2;
        tacc = tacc + length(t); tacc2 = tacc2 + length(t2);
        if ~isempty(Xi)
            [f,l] = evalfreeenergy(X,Ti,Gamma(t,:),Xi(t2,:,:),hmm,residuals,XX,[0 1 0 0 0]);
            fell = fell + sum(f); %  data likelihood
            ll = [l; ll];
        else
            [f,l] = evalfreeenergy(X,Ti,Gamma(t,:),[],hmm,residuals,XX,[0 1 0 0 0]);
            fell = fell + sum(f); %  data likelihood
            ll = [l; ll];
        end
    end
    fe = fe + fell;
else
    [fe,ll] = evalfreeenergy(data,T,Gamma,Xi,hmm,residuals);
    fe = sum(fe);
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function L = orderedNoiseLikelihood(noise, mu, varsigma, y)

% ORDEREDNOISELIKELIHOOD Likelihood of the data under the ORDERED noise model.
% FORMAT
% DESC returns the likelihood of a data set under the  ordered categorical noise model.
% ARG noise : the noise structure for which the likelihood is required.
% ARG mu : input mean locations for the likelihood.
% ARG varSigma : input variance locations for the likelihood.
% ARG y : target locations for the likelihood.
%
% SEEALSO : orderedNoiseParamInit, orderedNoiseLogLikelihood, noiseLikelihood
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005

% NOISE

D = size(y, 2);

L = zeros(size(mu));
L(find(isnan(y))) = 1;

c = 1./sqrt(noise.variance + varsigma);
for j = 1:D
  % Do lowest category first
  index = find(y(:, j) == 0);
  if ~isempty(index)
    mu(index, j) = mu(index, j) + noise.bias(j) ;
    L(index, j) = cumGaussian(-mu(index, j).*c(index, j));
  end

  % Intermediate categories
  index = find(y(:, j) > 0 & y(:, j) < noise.C-1);
  if ~isempty(index)
    for i = index'
      mu(i, j) = mu(i, j) + noise.bias(j) - sum(noise.widths(1:y(i, j)-1));
      L(i, j) = cumGaussian(mu(i, j).*c(i, j)) ...
          - cumGaussian((mu(i, j) ...
                         - noise.widths(y(i, j))) ...
                        .* c(i, j));
    end
  end
  % Highest category
  index = find(y(:, j) == noise.C-1);
  if ~isempty(index)
    for i = index'
      mu(i, j) = mu(i, j) + noise.bias(j) - sum(noise.widths(1:y(i, j)-1));
    end
    L(index, j) = cumGaussian(mu(index, j).*c(index, j));
  end
end
%/~ for j = 1:D
%   L(:, j) = (1-noise.C*noise.eta)*L(:, j)+noise.eta;
%~/ end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function Cube = RotateCube(hAxes,Cube,Axe,Direction)
% to rotate all cube
MyCube = Cube; % copy object
a1 = Axe; % 1st axe
a2 = Axe+1;
if a2>3,
  a2 = 1;
end % 2nd axe
a3 = 6-a1-a2; % 3rd axe
cosa = cos(Direction*pi/36); % cos(+-1 degree)
sina = sin(Direction*pi/36); % sin(+-1 degree)
T = zeros(3);
T(a1,a1) = 1;
T(a2,a2) = cosa;
T(a2,a3) = -sina;
T(a3,a2) = sina;
T(a3,a3) = cosa;
for k1=1:18, % algle loop
  for k2=1:27, % small cube loop
    MyCube.Nodes(:,:,k2) = MyCube.Nodes(:,:,k2)*T';
  end
  PlotCube(hAxes,MyCube);
  drawnow;
end
ind = 2*Axe-1+(1-Direction)/2; % index number
MyColor = zeros(size(Cube.Color));
for k1=1:27, % small cube loop
  for k2=1:6, % flat loop
    MyColor(Cube.RotCubeFlat(ind,k2),:,Cube.RotCubeCube(ind,k1)) = ...
      Cube.Color(k2,:,k1);
  end
end
Cube.Color = MyColor; % change colors
PlotCube(hAxes,Cube);
return

author: Sable

repo: mcbench-benchmarks


text:
function [criterion] = gbc_FMC(q, CS, Dl,Dr, cmaha, varargin)
%
%  cmaha      -
%  cmaha0     - lower misorientation bias (default 0.05)
%  quatmax    - quaterion variance metrix for cluster (default 5)
%  alpha      - seed selection (default 0.2)
%  beta       - probability threshold for point in cluster (default 0.3)%
%  gammaW     - edge dilution (default 25)

A_D = sparse(Dl,Dr,true,length(q),length(q));

fmc.CS = CS;
fmc.O  = q;

fmc.cmaha    = cmaha(1);
fmc.cmaha0   = get_option(varargin,{'cmaha0'},0.05,'double');
fmc.quatmax  = get_option(varargin,{'quatmax'},5,'double');
fmc.quatmax2 = cos(fmc.quatmax/2*degree);
fmc.gammaW   = get_option(varargin,{'gammaW'},25,'double');
fmc.alpha    = get_option(varargin,{'alpha'},0.2,'double');
fmc.beta     = get_option(varargin,{'beta'},0.3,'double');

fmc.A_D      = A_D;

[AllPs,AllSals,numClusters] = FMC_MTEX(fmc);

assignments = FMC_interpret(AllSals, numClusters, AllPs, A_D, fmc.beta);

criterion = assignments(Dl,1) == assignments(Dr,1) & assignments(Dl,1) > 0;

% I_DG = sparse(1:length(assignments),assignments(:,1),1);
% I_DG(:,all(I_DG==0,1)) = [];

author: mtex-toolbox

repo: mtex

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


text:
function tests = test_spm_dcm_peb_to_gcm
% Unit Tests for spm_dcm_peb_to_gcm
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% $Id: test_spm_dcm_peb_to_gcm.m 6875 2016-09-15 09:30:18Z peter $

tests = functiontests(localfunctions);

%--------------------------------------------------------------------------
function data_path = get_data_path()
data_path = fullfile( spm('Dir'), 'tests', ...
    'data', 'test_spm_dcm_peb_to_gcm');

%--------------------------------------------------------------------------
function test_2groups(testCase)
% Tests creating GCMs for two groups of subjects who differ in a single
% parameter

% Settings
n     = 100;        % Subjects
b_Hz  = 0.6;        % Strength of group difference
param = 'B(2,1,3)'; % Parameter to vary between groups
beta  = 16;         % Within:between variance ratio for simulated parameters

rng(1);

% Load template DCM
DCM_template = fullfile(get_data_path(), 'DCM_attention.mat');
DCM_template = load(DCM_template);
DCM_template = DCM_template.DCM;

% Indices of subjects in each group
group    = cell(1,2);
group{1} = 1:ceil(n/2);
group{2} = ceil(n/2)+1:n;

% PEB design matrix
X             = ones(n,2);
X(group{2},2) = -1;

% Create dummy PEB
PEB = create_peb(DCM_template, X, b_Hz, param, beta);

% Create GCM array
options = struct('nsubjects',n);
[GCM,PEB] = spm_dcm_peb_to_gcm(PEB,DCM_template,options);

% Check size
testCase.assertTrue( all(size(GCM)==[n 1]) );

% Get group arithmetic average of parameters
Eps = [];
for s = 1:n
    Eps(:,s) = spm_vec(GCM{s,1}.Ep);
end
Eps = Eps(PEB.Pind,:);

% Check subject' parameters are similar to PEB.Ep(:,1)
testCase.assertGreaterThanOrEqual(...
    one_sample_tt(mean(Eps,2)-PEB.Ep(:,1)),...
    0.05);

% Check group difference in B-parameter is similar to PEB.Ep(idx,2)
idx       = spm_fieldindices(GCM{1}.M.pE,param);
common    = PEB.Ep(PEB.Pind == idx,1);
attention = PEB.Ep(PEB.Pind == idx,2);
expected  = [common+attention common-attention]';

for g = 1:length(group)
    % Gather group's parameters
    Eps = [];
    for s = 1:length(group{g})
        Eps(:,s) = spm_vec(GCM{group{g}(s),1}.Ep);
    end

    % Check group mean is similar to PEB parameter
    testCase.assertGreaterThanOrEqual(...
        one_sample_tt(Eps(idx,:), expected(g)),...
        0.05);

    % Check that between-subjects variance of parameters is as expected
    testCase.assertGreaterThanOrEqual(...
        one_sample_tt(diag(PEB.Ce) - var(Eps(PEB.Pind,:),[],2)),...
        0.05);

end

%--------------------------------------------------------------------------
function PEB = create_peb(DCM, X, b_Hz, param, beta)
% Generate dummy PEB for simulation
%
% DCM   - Template DCM
% X     - PEB design matrix [parameters x regressors]
% b_Hz  - Group difference on a connection in Hz
% param - Parameter name to vary between groups
% beta  - Within:between variance ratio for simulated parameters

% Fields to include in the PEB
fields = {'A','B','C','transit','decay','epsilon'};

% Get within-subject priors from template DCM
[pE,pC] = spm_dcm_fmri_priors(DCM.a,DCM.b,DCM.c,DCM.d);

% Sample once from priors to get 'true' group-level parameter means
pE_struct = pE;
pE        = full(spm_vec(pE));
gP        = spm_normrnd(pE,pC,1);

% Force C-matrix (driving input) to be positive
idx = spm_fieldindices(pE_struct,'C');
gP(idx) = abs(gP(idx));

% Get indices of parameters to include in the PEB
Pind = spm_find_pC(pC,pE_struct,fields);

% Create group difference regressor
att = zeros(length(pE),1);
idx = spm_fieldindices(pE_struct,param);
att(idx) = b_Hz;

% Collate PEB
clear PEB;
PEB.Ep   = [gP(Pind) att(Pind)];
PEB.beta = beta;
PEB.Pind = Pind;
PEB.M.X  = X;

% -------------------------------------------------------------------------
function p = one_sample_tt(data,mu)
% One sample t-test

if nargin < 2, mu = 0; end
t = (mean(data)-mu) / (std(data) / sqrt(length(data)));
p = (1-spm_Tcdf(t,length(data)-1)) * 2;

author: spm


text:
function [estimate] = ft_inverse_pcc(csd, leadfield, varargin)

% FT_INVERSE_PCC
%
% Use as
%  estimate = ft_inverse_pcc(csd, leadfield, ...)

author: PatternRecognition

repo: OpenBMI


text:
function [Lc,L]=wfbtclength(Ls,wt,varargin)
%WFBTCLENGTH  WFBT subband lengths from a signal length
%   Usage: Lc=wfbtclength(Ls,wt);
%          [Lc,L]=wfbtclength(...);
%
%   `Lc=wfbtclength(Ls,wt)` returns the lengths of coefficient subbands
%   obtained from |wfbt| for a signal of length *Ls*. Please see the help
%   on |wfbt| for an explanation of the parameters *wt*.
%
%   `[Lc,L]=wfbtclength(...)` additionally returns the next legal length
%   of the input signal for the given extension type.
%
%   The function support the same boundary-handling flags as the |fwt|
%   does.
%
%   See also: wfbt, wfbtlength

% AUTHOR: Zdenek Prusa

complainif_notposint(Ls,'Ls','WFBTCLENGTH');

definput.import = {'fwt'};
[flags,kv]=ltfatarghelper({},definput,varargin);

% Initialize the wavelet filters structure
wt = wfbtinit(wt);

if(flags.do_per)
   a = treeSub(wt);
   L = filterbanklength(Ls,a);
   Lc = L./a;
else
   L = Ls;
   Lc = treeOutLen(L,0,wt);
end


text:
function printBERContribution( d, Nd, wd, N, R )
%PRINTBERCONTRIBUTION Print contribution from distance spectrum components.
% d, Nd, wd are arrays of numbers, each index corresponding to one spectral
% component. N is interleaver size. R is rate.

%close all;

upperrange = 3.5;

%SNR Range in dB
SNR_range = [0:0.01:upperrange];

ebno = 10 .^ (SNR_range ./ 10);

ber=zeros(length(d), length(ebno));

linetypes = {'r--', 'b', 'b--', 'r:', 'b:', 'r-.', 'b-.'};
leg = cell(1,length(d));
for i=1:length(d)
    ber(i+1,:) =  ber(i,:) + ((Nd(i) * wd(i)) / N) * qfunc (sqrt(d(i) * 2 * R * ebno));
    semilogy(SNR_range, ber(i+1,:), linetypes{i});
    hold on
    leg{i} = sprintf('d = %d', d(i));
end
legend(leg);


text:
function jf_rm_toolbox
% remove matlab toolboxes from path
% for testing my toolbox with vanilla matlab

rm aero
rm aeroblks
rm bioinfo
rm comm
rm commblks
rm compiler
rm control
rm curvefit
rm database
rm des
rm distcomp
rm dspblks
rm eml
rm emlcoder
rm filterdesign
rm finance
rm fixedpoint
rm fuzzy
rm geoweb
rm ident
rm idelink
rm images
rm instrument
%rm local
rm map
%rm matlab
rm nnet
rm optim
rm pde
rm physmod
rm rf
rm rfblks
rm robust
rm rtw
%rm shared
rm signal
rm simulink
rm sl3d
rm slcontrol
rm slvnv
rm splines
rm stateflow
rm stats
rm symbolic
rm vipblks
rm vision
rm wavelet

%path

function rm(box)
tmp = which('svd');
ii = strfind(tmp, 'matlab/matfun/@single/svd');
if numel(ii) ~= 1
	fail('bad "%s"', tmp)
end
tmp(ii:end) = '';
pre = strrep(tmp, 'built-in (', '');
%pre = '/Volumes/a2/pub/matlab/matlab-2010a.app/toolbox/';
tmp = genpath([pre box])
rmpath(tmp);
%keyboard

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function pdf = truncated_normal_ab_pdf ( x, mu, sigma, a, b )

%*****************************************************************************80
%
%% TRUNCATED_NORMAL_AB_PDF evaluates the truncated Normal PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 August 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%
%    Input, real MU, SIGMA, the mean and standard deviation of the
%    parent Normal distribution.
%
%    Input, real A, B, the lower and upper truncation limits.
%
%    Output, real PDF, the value of the PDF.
%
  alpha = ( a - mu ) / sigma;
  beta = ( b - mu ) / sigma;
  xi = ( x - mu ) / sigma;

  alpha_cdf = normal_01_cdf ( alpha );
  beta_cdf = normal_01_cdf ( beta );
  xi_pdf = normal_01_pdf ( xi );

  pdf = xi_pdf / ( beta_cdf - alpha_cdf ) / sigma;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [str, tag, cind, ccnt] = gencode_item(item, tag, tagctx, stoptag, tropts)

% function [str, tag, cind, ccnt] = gencode_item(item, tag, tagctx, stoptag, tropts)
% Generate code to recreate a generic item. This code does not deal with
% arrays of cfg_items, such a configuration should not exist with the
% current definition of a configuration tree.
%
% Traversal options
% struct with fields
% stopspec - match spec to stop forced setting of eflag
% dflag    - (not used here)
% clvl     - current level in tree
% mlvl     - maximum level to force settings - range 1 (top level only) to
%            Inf (all levels)
% cnt      - item count - used for unique tags
% mcnt     - (not evaluated here)
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: gencode_item.m 5678 2013-10-11 14:58:04Z volkmar $

rev = '$Rev: 5678 $'; %#ok

%% Parent object
% Generate generic object
[str, tag, cind, ccnt] = gencode_item(item.cfg_item, tag, tagctx, stoptag, tropts);
% Check whether to generate code - ccnt == 0 means that generic object did
% not return code
if (tropts.clvl > tropts.mlvl || (~isempty(tropts.stopspec) && match(item, tropts.stopspec))) || ccnt == 0
    str = {};
    cind = [];
    ccnt = 0;
    return;
end;
% Reclassify generic object
str{cind} = sprintf('%s         = %s;', tag, class(item));
%% Filter, Dir, Ufilter
% Generate string type fields
fn = {'filter','dir','ufilter'};
for k = 1:numel(fn)
    if ~isempty(item.(fn{k}))
        str1 = gencode(item.(fn{k}), sprintf('%s.%s', tag, fn{k}), ...
                       tagctx);
        str = [str(:)' str1(:)'];
    end;
end;
%% Num
% Generate num field
str{end+1} = sprintf('%s.num     = [%d %d];', tag, item.num);


text:
% std_chaninds() - look up channel indices in a STUDY
%
% Usage:
%         >> inds = std_chaninds(STUDY,  channames);
%         >> inds = std_chaninds(EEG, channames);
% Inputs:
%         STUDY - studyset structure containing a changrp substructure.
%     channames - [cell] channel names
%
% Outputs:
%       inds - [integer array] channel indices
%
% Author: Arnaud Delorme, CERCO, 2006-

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

function finalinds = std_chaninds(STUDY, channames);

    finalinds   = [];
    if isfield(STUDY, 'chanlocs')
        EEG = STUDY;
        tmpchanlocs = EEG.chanlocs;
        tmpallchans = lower({ tmpchanlocs.labels });
    else
        tmpallchans = lower({ STUDY.changrp.name });
    end;

    if isempty(channames), finalinds = [1:length(tmpallchans)]; return; end;
    for c = 1:length(channames)
        if isnumeric(channames)
            chanind = channames(c);
        else
            chanind = strmatch( lower(channames{c}), tmpallchans, 'exact');
            if isempty(chanind), error(sprintf([ 'Channel "%s" and maybe others was not' 10 'found in pre-computed data file' ], channames{c})); end;
        end;
        finalinds   = [ finalinds chanind ];
    end;


text:
function [ t_start, y_start ] = p14_start ( neqn )

%*****************************************************************************80
%
%% P14_START returns the starting point for problem p14.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_START, Y_START(NEQN,1), the initial data.
%
  y_start = zeros ( neqn, 1 );

  t_start = 0.0;
  y_start(1) = 1.0;
  y_start(2:51) = 0.0;

  return
end


text:
classdef Optimizer_PrimalDual < Optimizer

    properties (GetAccess = public, SetAccess = protected)
        unconstrainedOptimizer
    end

    properties (Access = protected)
        NSmerit
        NullSpaceSettings
        lagrangian
        lagrangianSettings
    end

    properties (Access = private)

    end

    methods (Access = public)

        function solveProblem(obj)
            obj.cost.computeFunctionAndGradient();
            obj.constraint.computeFunctionAndGradient();
            obj.lagrangian.updateBecauseOfPrimal();
            obj.unconstrainedOptimizer.startLineSearch();
            obj.printOptimizerVariable();
            obj.hasFinished = false;

            while ~obj.hasFinished
                obj.increaseIter();
                obj.update();
                obj.updateStatus();
                obj.refreshMonitoring();
                obj.printOptimizerVariable();
                %obj.printHistory();
            end
            obj.printOptimizerVariable();
            obj.printHistory();

            obj.hasConverged = 0;
            obj.printHistoryFinalValues();
        end

    end

    methods (Access = protected)

        function updateOldValues(obj)
            obj.designVariable.updateOld();
            obj.dualVariable.updateOld();
            obj.cost.updateOld();
            obj.constraint.updateOld();
            obj.updateLagrangian();
            obj.lagrangian.updateOld();
        end

        function updateConvergenceStatus(obj)
            isOptimal   = obj.unconstrainedOptimizer.isOptimal();
            isFeasible  = obj.isFeasible();
            obj.hasConverged = isOptimal && isFeasible;
        end

        function itIs = isFeasible(obj)
            active_constr    = true(size(obj.dualVariable.value));
            constraintValues = abs(obj.constraint.value(active_constr));
            constrTol        = obj.targetParameters.constr_tol(active_constr);
            isNotFeasible = any(any(constraintValues > constrTol));
            itIs = ~isNotFeasible;
        end

        function createLagrangian(obj)
            obj.createLagrangianSettings();
            cParams = obj.lagrangianSettings;
            obj.lagrangian = ObjectiveFunction.create(cParams);
        end

        function createNullSpaceMerit(obj)
            obj.createNullSpaceSettings();
            cParams = obj.NullSpaceSettings();
            obj.NSmerit = ObjectiveFunction.create(cParams);
        end

       function createOptimizerUnconstrained(obj,cParams)
            cParams.lagrangian      = obj.lagrangian;
            cParams.convergenceVars = obj.convergenceVars;
            obj.unconstrainedOptimizer = Optimizer_Unconstrained.create(cParams);
       end

       function createOptimizerUnconstrainedNS(obj,cParams)
           cParams.NSmerit         = obj.nsMerit;
           cParams.convergenceVars = obj.convergenceVars;
           obj.unconstrainedOptimizer = Optimizer_Unconstrained.create(cParams);
       end

        function updateLagrangian(obj)
            obj.lagrangian.computeFunction();
            obj.lagrangian.computeGradient();
        end

    end

    methods (Access = private)

        function increaseIter(obj)
            obj.nIter = obj.nIter+1;
        end

    end

    methods (Access = protected, Abstract)
        createLagrangianSettings(obj)
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [rphase, fh] = tapas_physio_get_respiratory_phase(pulset, rsampint,...
    verbose, thresh)

% get_respiratory_phase is a function for creating respiratory phase regressor.
% from physiological monitoring files acquired using spike, that
% can be included in SPM5 design matrices.
%
% FORMAT [rphase] = tapas_physio_get_respiratory_phase(pulset,rsampint)
%
% Inputs:
%        pulset - respiratory belt data read from spike file
%        rsampint - sampling frequency
%        verbose - 1 for graphical output (debugging mainly)
%        thresh.resp_max -  if set, all peaks above that breathing belt amplitude
%                        are ignored for respiratory phase histogram evaluation
%
% Outputs:
%        rphase - respiratory phase when each slice of each volume was
%        acquired
%
% The regressors are calculated in the following way as described in
% Glover et al, 2000, MRM, (44) 162-167:
% 1) Find the max and min amplitudes from belt response
% 2) Normalise the respiratory amplitude
% 3) Calculate the histogram from the number of occurences of specific respiratory
%    amplitudes in bins 1:100.
% 4) Calculate the running integral of the histogram for the bin
%    corresponding to each respiratory amplitude
%_______________________________________________________________________

% Author: Lars Kasper
%         Error handling for temporary breathing belt failures:
%         Eduardo Aponte, TNU Zurich
%
% Copyright (C) 2013, Institute for Biomedical Engineering, ETH/Uni Zurich.
%
% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

%% get histogram of amplitudes of breathing belt signal
if nargin < 3
    verbose = false;
end

if nargin < 4
    resp_max = inf;
elseif ~isfield(thresh, 'resp_max')
    resp_max = inf;
elseif isempty(thresh.resp_max)
    resp_max = inf;
else
    resp_max = thres.resp_max;
end

% Check input

assert(all(abs(pulset) ~= inf), ...
    'Infinity values in the respiratory regressor');
assert(~any(isnan(pulset)), 'Nan values in the respiratory regressors');

% Weird line...
overshoot = find(abs(pulset) > resp_max);
pulset(overshoot) = resp_max;

maxr = max(pulset);
minr = min(pulset);

% Compute normalized signal and the sign of the derivative

npulse = (pulset-minr)/(maxr-minr);

% Calculate derivative of normalised pulse wrt time
% over 1 sec of data as described in Glover et al.
ksize = round(0.5 * (1/rsampint));
kernel = ksize:-1:-ksize; kernel = kernel ./ sum(kernel.^2);
dpulse = tapas_physio_conv(pulset, kernel, 'symmetric');
% This uses a quadratic Savitzky-Golay filter, for which the coefficients
% have a simple linear form. See e.g.
% https://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter
% `scipy.signal.savgol_coeffs(5, polyorder=2, deriv=1, use='conv')`

% Tolerance to the derivative
dpulse(abs(dpulse) < 1e-4) = 0;
dpulse = sign(dpulse);

% number of histogram bins determined by dynamic range of detected values
% and length of input time course
nbins = min(length(unique(pulset)), floor(length(pulset)/100));

[h, rout] = hist(npulse(dpulse ~=0 & npulse < resp_max), nbins);

binnum = floor(npulse*(nbins-1)) + 1;
binnum(overshoot) = nbins;

cumsumh = cumsum(h');
sumh = cumsumh(end);

dpulse(dpulse == 0) = nan;
rphase = pi*(cumsumh(binnum)/sumh).*dpulse+pi;

if verbose
    fh = plot_traces(pulset, rsampint, rout, resp_max, cumsumh, sumh, h, ...
        npulse, dpulse, rphase);
else
    fh = [];
end

end

function fh = plot_traces(pulset, rsampint, rout, resp_max, ...
    cumsumh, sumh, h, npulse, dpulse, rphase)

nsamples = length(pulset);
t = (0:nsamples-1)*rsampint;
feqht = cumsumh/sumh*pi;

fh = tapas_physio_get_default_fig_params();
set(fh, 'Name', ...
   'Preproc: get_respiratory_phase: histogram for respiratory phase estimation');

hs(1) = subplot(2,2,1);
plot(t,pulset);
xlabel('t (s)');
ylabel('breathing amplitude (a. u.)');
title('(filtered) breathing time series');

if resp_max < inf
    hold on;
    plot(t, ones(size(t)) * resp_max, 'k--');
    hold on;
    hp = plot(t, -ones(size(t)) * resp_max, 'k--');
    legend(hp, ...
        'threshold for maximum amplitude to be considered in histogram');
    set(gcf, 'Name', ...
        [get(gcf, 'Name') ' - with amplitude overshoot-correction']);
end

hs(2) = subplot(2,2,2);
bar(rout, h);
xlabel('normalized breathing amplitude');
ylabel('counts');
title('histogram for phase mapping');
xlim([-0.1 1.1]);

hs(3) = subplot(2,2,3); plot(rout, [feqht, cos(feqht), sin(feqht)]);
xlabel('normalized breathing amplitude');
title(...
    'equalized histogram bin amplitude -> phase transfer function (f_{eqht})');
legend('f: normalized amplitude -> phase transfer function', 'cos(f)', ...
    'sin(f)', 'Location', 'NorthWest');

%figure('Name', 'Histogram: Respiration phase estimation');
hs(4) = subplot(2,2,4);
plot(t, [npulse*10, dpulse, (rphase-pi)]);
legend('10*normalized breathing belt amplitude', ...
    '-1 = exhale, 1 = inhale', 'estimated respiratory phase');
ylim([-0.2 10.2]);
title('Histogram-based respiration phase estimation');

linkaxes(hs([1 4]), 'x');
linkaxes(hs([2 3]), 'x');

end


text:
function [c,t,i,CC,LL] = interferogram(c, width, timestep, masterIdx)

   % C = INTERFEROGRAM(C)
   % This function calculates an interferogram from the waveforms in a
   % correlation object against the a single "master" trace.
   %
   % The routine interatively selects a narrow time window on all traces,
   % carries out a cross-correlation, and stores the maximum correlation value
   % and the time lag necessary to acheive this best fit for each time window.
   % All cross correlations include sub-sample interpolation of lag values
   % (see HELP XCORR for details). The results are two matrices of correlation
   % value and lag time. These can be plotted with the C-INTERFER and
   % L-INTERFER options in correlation plot.
   %
   % [C,TIME,INDEX,MAXCORR,LAG] = INTERFEROGRAM(C,...) In addition to the
   % correlation object, this use returns a vector of window center times
   % relative to trigger times (1xM), a vector of trace indices (Nx1), a
   % matrix of maximum correlation values (NxM), and a matrix of optimum lag
   % times (NxM). These values can be used to generate a custom interferogram
   % figure, for example:
   %    imagesc(TIME,INDEX,MAXCORR)
   %
   % [C ...] = INTERFEROGRAM(C,WIDTH) Specify the half width of the time window
   % explicitly. Default value is 1/20 of the mean trace length.
   %
   % [C ...] = INTERFEROGRAM(C,WIDTH,TIMESTEP) Specify the timestep interval
   % explicitly. If TIMESTEP is a scalar, it is treated as the interval
   % between adjacent time windows. If TIMESTEP is a vector, the values are
   % treated as the specific center time times of each window.
   %
   % [C ...] = INTERFEROGRAM(C,WIDTH,TIMESTEP,TRACENUM) Specify which trace
   % should be treated as the master trace against which all correlation
   % values will be based.
   %
   % ** NOTE ABOUT DATA STORAGE **
   % The interferogram routine by default places the four output terms
   % (TIME,INDEX,MAXCORR and LAG) into supplimental fields in the waveform
   % object which holds the first trace. This is embarassingly clunky. However
   % it works and allows the data to be passed directly to the plotting
   % routine. These terms can always be accessed by extracting them from the
   % first waveform using the field names INTERFEROGRAM_TIME,
   % INTERFEROGRAM_INDEX, INTERFEROGRAM_MAXCORR, INTERFEROGRAM_LAG. Because of
   % this ad hoc structure though, in many cases it may be better to simply
   % save the data straight out of the routine and handle it yourself. It is
   % automatically placed in the first waveform however to fascilitate a quick
   % plotting option. INTERFEROGRAM_TIME is given in units of seconds relative
   % to the trigger time.

   % Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   % GET INPUT PARAMETERS

   % TEST FOR ALIGNED TRACES
   if ~check(c,'OFFSET')
      disp('Data not uniformly aligned around triggers. Consider using CROP function.');
   end

   % SET WINDOW WIDTH
   srt_offset = ( c.traces.firstsampletime() - c.trig ) * 86400;
   end_offset = ( c.traces.lastsampletime() - c.trig ) * 86400;

   if ~exist('width','var') || isempty(width)
      width = ( mean(end_offset)-mean(srt_offset) ) / 20;
   end

   if ~exist('timestep','var') || isempty('length')
      timestep = ( mean(end_offset)-mean(srt_offset) ) / 50;
   end

   if length(timestep)  > 1
      % user specified center time for each window
      t = timestep;
   else
      t = mean(srt_offset) : timestep : mean(end_offset);
   end

   if ~exist('masterIdx','var')
      masterIdx = c.ntraces; % last trace
   end

   fprintf('Time step: %4.3f    Window width: %4.3f    Reference trace no.: %2d\n', timestep, width, masterIdx);

   % CALC INTERFEROGRAM VOLUME
   CC = [];
   LL = [];
   times = [];
   for i = 1:length(t);
      times = cat(1,times,t(i));
      if (mod(i,10)==1)
         disp(['Center of time window: ' num2str(t(i),'%4.3f') ' ...']);
      end
      ctmp = c.crop(t(i)-width,t(i)+width);
      ctmp = xcorr(ctmp,'row',masterIdx);
      corr = ctmp.corrmatrix;
      lag  = ctmp.lags;
      CC(:,i) = corr(:,masterIdx);
      LL(:,i) = lag(:,masterIdx);
   end;

   % PLACE OUTPUT INTO FIRST WAVEFORM FIELD
   i = [1:c.ntraces]';
   T = c.traces(1);
   T.userdata.interferogram_index = i;
   T.userdata.interferogram_time = t;
   T.userdata.interferogram_maxcorr = CC;
   T.userdata.interferogram_lag = LL;
   c.traces(1) = T;
end


text:
function plotOptiProb(prob,opts,xb,scale,dolog,npts,mode)
%PLOTOPTIPROB Plot opti problem optimization surface with constraints
%
%   plotOptiProb(prob,xb,scale) plots a contour plot of the optimization
%   surface where prob is an optiprob structure. xb is the solution vector,
%   and scale zooms the plot.

%   Copyright (C) 2011 Jonathan Currie (I2C2)

%Defaults
if(nargin < 7 || isempty(mode)), mode = 'normal'; end
if(nargin < 6), npts = []; end
if(nargin < 5 || isempty(dolog)), dolog = false; end
if(nargin < 4 || isempty(scale)), scale = 5; end
if(nargin < 3), xb = []; end
if(nargin < 2), error('This function requires two inputs (prob, opts)'); end

ptypes = {'LP','BILP','MILP','QP','QCQP','MIQP','MIQCQP','SDP','MISDP','NLS','DNLS','UNO','NLP','MINLP'};

%Check we have a 'plottable' problem
if(~any(strcmpi(prob.type,ptypes)))
    error('A Problem of type ''%s'' cannot be plotted!',prob.type);
end

%Check for SD problem in SeDuMi format
if(any(strcmpi(prob.type,{'SDP','MISDP'})))
   if(isstruct(prob.sdcone))
       prob = sedumi2opti(prob);
       prob.objective = @(x) prob.f'*x;
   end
end

%Clear existing figure
clf;

%Check if unsolved
if(isempty(xb) && ~strcmpi(mode,'multi'))
    %Check if we have valid x0
    if(~isempty(prob.x0) && ~all(isnan(prob.x0)))
        mode = 'usex0';
        xb = prob.x0;
    elseif(~isempty(prob.lb) && ~isempty(prob.ub) && ~all(isinf(prob.lb)) && ~all(isinf(prob.ub)))
        mode = 'bounded';
        %fake xb for >2D plotting purposes
        xb = (prob.ub-prob.lb)./2 + prob.lb;
        %if user has supplied scale as vector, substitute here
        if(length(scale) > 1)
            mode = 'bounded_scale';
        end
    else
        error('OPTI can only plot unsolved problems that have a specified x0, or if all variables are bounded.');
    end
end

%Only plot fit if normal plot
if(any(strcmpi(mode,{'normal','usex0'})))
    %Check for parameter estimation problem
    if(strcmpi(prob.type,'DNLS'))
        plotDNLS(prob,opts,xb);
        return;
    elseif(strcmpi(prob.type,'NLS'))
        plotDataFit(prob,xb,dolog);
        return;
    end
end

%Setup Plotting Data
data.ndec = prob.sizes.ndec;
if(dolog ~= 1 && dolog ~= 0)
    error('doLog must be 0 or 1!');
end
data.dolog = dolog;
%Error check scale input
if(length(scale) > 1)
    if(length(scale) ~= data.ndec*2)
        error('Expected a %d element vector for specifying the plot boundaries, of the form [x1min x1max x2min x2max ... xNmin xNmax]',data.ndec*2);
    end
end
data.scale = scale;
%Check number of points
if(~isempty(npts) && (npts < 0 || npts > 3000))
    error('npts must be 0 < n < 3000');
end
data.npts = npts;
data.mode = mode;
data.fixval = xb;
data.idx = 1:data.ndec;
k = 1; n = 1;
%Check for higher dim plot
if(data.ndec > 2)
    for i = 1:data.ndec
        for j = i+1:data.ndec
            subplot(data.ndec-1,data.ndec-1,k);
            data.idx = [i j];
            if(length(scale) > 1) %grab required elements
                data.scale = [scale(i*2-1:i*2) scale(j*2-1:j*2)];
            end
            plotProblem(prob,xb,data);
            title(makeTitle(data,xb));
            k = k + 1;
            n = n + 1;
        end
        k = k + i;
    end
    %Plot Minimum on Figure Title
    if(~isempty(xb))
        data.fval = prob.objective(xb);
    else
        data.fval = [];
    end
    set(gcf,'Name',plotTitle(prob,xb,data));
else
    %Standard single plot
    plotProblem(prob,xb,data);
    set(gcf,'Name','');
end

function plotProblem(prob,xb,data)
%Plot Objective Contour
data.npts = plotObj(prob,xb,data);

%If we have search space (i.e. multisolve has been run), also plot
if(strcmp(data.mode,'multi') && isfield(prob,'multi') && ~isempty(prob.multi))
    plotMultiSearch(prob,data);
end

%Plot Bounds
if(~strcmpi(prob.type,'bilp'))
    plotBounds(prob.lb,prob.ub,data);
end

%Plot Linear Constraints
if(prob.sizes.nineq + prob.sizes.neq > 0)
    if(data.ndec <= 2)
        %Standard linear patches
        if(~isempty(prob.rl))
            [A,b,Aeq,beq] = row2gen(prob.A,prob.rl,prob.ru);
            plotLinCon(A,b,Aeq,beq);
        else
            plotLinCon(prob.A,prob.b,prob.Aeq,prob.beq);
        end
    else
        %Lazy nonlinear contour method (can't quite work out above for multi-dim)
        pcopy = prob; pcopy.nlcon = []; pcopy.nlrhs = []; pcopy.nle = []; pcopy.cl = []; pcopy.cu = [];
        if(~isempty(prob.rl))
            pcopy = rowlin2nl(pcopy);
        else
            pcopy = genlin2nl(pcopy);
        end
        %Plot patches by telling algorithm constraints are linear
        plotNonlinCon(pcopy,data,true);
    end
end

%Plot Quadratic Constraints
if(prob.sizes.nqc)
    plotQuadCon(prob.Q,prob.l,prob.qrl,prob.qru,data);
end

%Plot Semidefinite Constraints
if(prob.sizes.nsdcone)
    plotSDCon(prob.sdcone,data);
end

%Plot Nonlinear Constraints
if(~isempty(prob.nle) || ~isempty(prob.cl))
    plotNonlinCon(prob,data);
end

%Check for Integer Constraints
if(any(prob.int.ind))
    plotIntCon(prob,data);
end

%Plot Optimum
if(~isempty(xb))
    idx = data.idx;
    hold on
    switch(data.mode)
        case {'normal','multi'}
            if(length(xb)==1)
                plot(xb(idx(1)),prob.objective(xb),'r.','markersize',20);
            else
                plot(xb(idx(1)),xb(idx(2)),'r.','markersize',20);
            end
        case 'usex0'
            if(length(xb)==1)
                plot(xb(idx(1)),prob.objective(xb),'b.','markersize',20);
            else
                plot(xb(idx(1)),xb(idx(2)),'b.','markersize',20);
            end
    end
    hold off
end

function str = makeTitle(data,xb)
%Construct multi-dim plot title

xs = true(data.ndec,1);
xs(data.idx) = false;
fi = find(xs);
if(length(fi) > 1)
    str = '[';
else
    str = [];
end
for i = 1:length(fi)
    str = sprintf('%sx%d, ',str,fi(i));
end
str = str(1:end-2);
if(~isempty(xb))
    if(length(fi) > 1)
        str = [str '] = ['];
    else
        str = [str ' = '];
    end
    for i = 1:length(fi)
        str = sprintf('%s%g, ',str,xb(fi(i)));
    end
    str = str(1:end-2);
    if(length(fi) > 1)
         str = [str ']'];
    end
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function Famp = vec2amp(F_cell)

chkarg(istypesizeof(F_cell, 'Scalar3dcell', [1 Axis.count]), ...
	'"E_cell should be length-%d cell array with Scalar3d as elements"', Axis.count);

Fx = F_cell{Axis.x};
grid3d = Fx.grid3d;
osc = Fx.osc;
physQcell = Fx.physQcell;
name = Fx.name(1:end-2);

li = grid3d.l(:,GT.dual);
[Xi, Yi, Zi] = ndgrid(li{:});

fi = NaN([grid3d.N, Axis.count]);
for w = Axis.elems
	Fw = F_cell{w};
	[fw, l] = Fw.data_ghost_expanded();
	[X, Y, Z] = ndgrid(l{:});
	fi(:,:,:,w) = interpn(X, Y, Z, fw, Xi, Yi, Zi);
end

famp = sqrt(sum(abs(fi).^2, Axis.count+1));
assert(all(size(famp)==grid3d.N));

% Attach extra points.
for w = Axis.elems
	famp = attach_extra_Famp(famp, w, grid3d);
end

gt_array = [GT.dual GT.dual GT.dual];

Famp = Scalar3d(famp, grid3d, gt_array, osc, physQcell, ['|', name, '|']);

function array = attach_extra_Famp(array, w, grid3d)
% This makes Painter2d.prep_data() estimates the field values at the boundary of
% the computational domain inaccurate.  For example, suppose that we have Ez on
% the xy-plane with a PEC boundary at the x-boundaries.  Then, Ez = 0 at x = 0,
% but for Eamp we have nonzero |Ez| before and after x = 0, leading to a nonzero
% estimate of |Ez| at x = 0.  We cannot handle this case correctly in the
% current visualization framework...  It might be better to estimate Eamp and
% the corners of Yee's cells rather than centers?

ind_n = {':', ':', ':'};
ind_p = {':', ':', ':'};
bc_w = grid3d.bc(w);
if bc_w == BC.p
	ind_n{w} = grid3d.N(w);
	ind_p{w} = 1;
else  % bc_d == BC.e or BC.m
	ind_n{w} = 1;
	ind_p{w} = grid3d.N(w);
end
array = cat(int(w), array(ind_n{:}), array, array(ind_p{:}));  % Bloch phases in S are ignored in amplitude

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function [res,resL2,qualMeasOut]=SIRT(proj,geo,angles,niter,varargin)
% SIRT solves Cone Beam CT image reconstruction using Oriented Subsets
%              Simultaneous Algebraic Reconstruction Technique algorithm
%
%   SIRT(PROJ,GEO,ALPHA,NITER) solves the reconstruction problem
%   using the projection data PROJ taken over ALPHA angles, corresponding
%   to the geometry described in GEO, using NITER iterations.
%
%   SIRT(PROJ,GEO,ALPHA,NITER,OPT,VAL,...) uses options and values for solving. The
%   possible options in OPT are:
%
%
%   'lambda':      Sets the value of the hyperparameter. Default is 1
%
%   'lambda_red':  Reduction of lambda. Every iteration
%                  lambda=lambdared*lambda. Default is 0.95
%
%   'Init':        Describes different initialization techniques.
%                  'none'     : Initializes the image to zeros (default)
%                  'FDK'      : Initializes image to FDK reconstruction
%                  'multigrid': Initializes image by solving the problem in
%                               small scale and increasing it when relative
%                               convergence is reached.
%                  'image'    : Initialization using a user specified
%                               image. Not recommended unless you really
%                               know what you are doing.
%   'InitImg'      an image for the 'image' initialization. Avoid.
%
%   'Verbose'      1 or 0. Default is 1. Gives information about the
%                  progress of the algorithm.
%   'QualMeas'     Asks the algorithm for a set of quality measurement
%                  parameters. Input should contain a cell array of desired
%                  quality measurement names. Example: {'CC','RMSE','MSSIM'}
%                  These will be computed in each iteration.
% 'redundancy_weighting': true or false. Default is true. Applies data
%                         redundancy weighting to projections in the update step
%                         (relevant for offset detector geometry)
%  'groundTruth'  an image as grounf truth, to be used if quality measures
%                 are requested, to plot their change w.r.t. this known
%                 data.
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Ander Biguri
%--------------------------------------------------------------------------

%% Deal with input parameters

[lambda,res,lambdared,verbose,QualMeasOpts,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,varargin);

measurequality=~isempty(QualMeasOpts) | ~any(isnan(gt(:)));
if ~any(isnan(gt(:)))
    QualMeasOpts{end+1}='error_norm';
    res_prev=gt;
    clear gt
end
if nargout<3 && measurequality
    warning("Image metrics requested but none catched as output. Call the algorithm with 3 outputs to store them")
    measurequality=false;
end
qualMeasOut=zeros(length(QualMeasOpts),niter);

resL2=zeros(1,niter);
if nargout>1
    computeL2=true;
else
    computeL2=false;
end

%% Create weighting matrices

% Projection weight, W
% Projection weight, W
W=computeW(geo,angles,gpuids);

% disp('Size of W matrix');
% disp(size(W));
clear geoaux

% Back-Projection weight, V
V=computeV(geo,angles,{angles},{1:length(angles)},'gpuids',gpuids);

if redundancy_weights
    % Data redundancy weighting, W_r implemented using Wang weighting
    % reference: https://iopscience.iop.org/article/10.1088/1361-6560/ac16bc

    num_frames = size(proj,3);
    W_r = redundancy_weighting(geo);
    W_r = repmat(W_r,[1,1,num_frames]);
    % disp('Size of redundancy weighting matrix');
    % disp(size(W_r));
    W = W.*W_r; % include redundancy weighting in W
end

%% hyperparameter stuff
nesterov=false;
if ischar(lambda)&&strcmp(lambda,'nesterov')
    nesterov=true;
    lambda=(1+sqrt(1+4))/2;
    gamma=0;
    ynesterov=zeros(size(res),'single');
    ynesterov_prev=ynesterov;
end
%% Iterate

% TODO : Add options for Stopping criteria
for ii=1:niter
    if (ii==1 && verbose==1);tic;end
    % If quality is going to be measured, then we need to save previous image
    % THIS TAKES MEMORY!
    if measurequality && ~strcmp(QualMeasOpts,'error_norm')
        res_prev = res; % only store if necesary
    end
    % --------- Memory expensive-----------
    %         proj_err=proj-Ax(res,geo,angles);                 %                                 (b-Ax)
    %         weighted_err=W.*proj_err;                         %                          W^-1 * (b-Ax)
    %         backprj=Atb(weighted_err,geo,angles);             %                     At * W^-1 * (b-Ax)
    %         weigth_backprj=bsxfun(@times,1./V,backprj);       %                 V * At * W^-1 * (b-Ax)
    %         res=res+lambda*weigth_backprj;                    % x= x + lambda * V * At * W^-1 * (b-Ax)
    % ------------------------------------
    % --------- Memory cheap(er)-----------

    if nesterov
        % The nesterov update is quite similar to the normal update, it
        % just uses this update, plus part of the last one.
        ynesterov=res + bsxfun(@times,1./V,Atb(W.*(proj-Ax(res,geo,angles,'gpuids',gpuids)),geo,angles,'gpuids',gpuids));
        res=(1-gamma)*ynesterov+gamma*ynesterov_prev;
    else
        res=res+lambda*bsxfun(@times,1./V,Atb(W.*(proj-Ax(res,geo,angles,'gpuids',gpuids)),geo,angles,'gpuids',gpuids)); % x= x + lambda * V * At * W^-1 * (b-Ax)
    end
    % ------------------------------------

    if nonneg
        res=max(res,0);
    end
    if nesterov
        gamma=(1-lambda);
        lambda=(1+sqrt(1+4*lambda^2))/2;
        gamma=gamma/lambda;
    else
        lambda=lambda*lambdared;
    end
    % If quality is being measured
    if measurequality
        % HERE GOES
        qualMeasOut(:,ii)=Measure_Quality(res_prev,res,QualMeasOpts);
    end

    if computeL2 || nesterov
        resL2(ii)=im3Dnorm(proj-Ax(res,geo,angles,'gpuids',gpuids),'L2','gpuids',gpuids); % Compute error norm2 of b-Ax
        % If the error is not minimized.
        if  ii~=1 && resL2(ii)>resL2(ii-1)
            if verbose
                disp(['Convergence criteria met, exiting on iteration number:', num2str(ii)]);
            end
            return
        end

    end

    if (ii==1 && verbose==1)
        expected_time=toc*niter;
        disp('SIRT');
        disp(['Expected duration   :    ',secs2hms(expected_time)]);
        disp(['Expected finish time:    ',datestr(datetime('now')+seconds(expected_time))]);
        disp('');
    end
end

end

function initres=init_multigrid(proj,geo,alpha)

finalsize=geo.nVoxel;
% start with 64
geo.nVoxel=[64;64;64];
geo.dVoxel=geo.sVoxel./geo.nVoxel;
if any(finalsize<geo.nVoxel)
    initres=zeros(finalsize');
    return
end
niter=100;
initres=zeros(geo.nVoxel');
while ~isequal(geo.nVoxel,finalsize)

    % solve subsampled grid
    initres=SIRT(proj,geo,alpha,niter,'Init','image','InitImg',initres,'Verbose',0,'gpuids',gpuids);

    % Get new dims.
    geo.nVoxel=geo.nVoxel*2;
    geo.nVoxel(geo.nVoxel>finalsize)=finalsize(geo.nVoxel>finalsize);
    geo.dVoxel=geo.sVoxel./geo.nVoxel;
    % Upsample!
    % (hopefully computer has enough memory............)
    [y, x, z]=ndgrid(linspace(1,size(initres,1),geo.nVoxel(1)),...
        linspace(1,size(initres,2),geo.nVoxel(2)),...
        linspace(1,size(initres,3),geo.nVoxel(3)));
    initres=interp3(initres,x,y,z);
    clear x y z
end

end

function [lambda,res,lambdared,verbose,QualMeasOpts,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,alpha,argin)
opts={'lambda','init','initimg','verbose','lambda_red','qualmeas','nonneg','gpuids','redundancy_weighting','groundtruth'};
defaults=ones(length(opts),1);
% Check inputs
nVarargs = length(argin);
if mod(nVarargs,2)
    error('TIGRE:SIRT:InvalidInput','Invalid number of inputs')
end

% check if option has been passed as input
for ii=1:2:nVarargs
    ind=find(ismember(opts,lower(argin{ii})));
    if ~isempty(ind)
        defaults(ind)=0;
    else
        error('TIGRE:SIRT:InvalidInput',['Optional parameter "' argin{ii} '" does not exist' ]);
    end
end

for ii=1:length(opts)
    opt=opts{ii};
    default=defaults(ii);
    % if one option is not default, then extract value from input
    if default==0
        ind=double.empty(0,1);jj=1;
        while isempty(ind)
            ind=find(isequal(opt,lower(argin{jj})));
            jj=jj+1;
        end
        if isempty(ind)
            error('TIGRE:SIRT:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]);
        end
        val=argin{jj};
    end

    switch opt
        % % % % % % % Verbose
        case 'verbose'
            if default
                verbose=1;
            else
                verbose=val;
            end
            if ~is2014bOrNewer
                warning('TIGRE: Verbose mode not available for older versions than MATLAB R2014b');
                verbose=false;
            end
            % % % % % % % hyperparameter, LAMBDA
        case 'lambda'
            if default
                lambda=1;
            elseif ischar(val)&&strcmpi(val,'nesterov')
                lambda='nesterov'; % just for lowercase/uppercase
            elseif length(val)>1 || ~isnumeric(val)
                error('TIGRE:SIRT:InvalidInput','Invalid lambda')
            else
                lambda=val;
            end
        case 'lambda_red'
            if default
                lambdared=1;
            else
                if length(val)>1 || ~isnumeric(val)
                    error('TIGRE:SIRT:InvalidInput','Invalid lambda')
                end
                lambdared=val;
            end
        case 'init'
            res=[];
            if default || strcmp(val,'none')
                res=zeros(geo.nVoxel','single');
                continue
            end
            if strcmp(val,'FDK')
                res=FDK(proj,geo,alpha);
                continue
            end
            if strcmp(val,'multigrid')
                res=init_multigrid(proj,geo,alpha);
                continue
            end
            if strcmp(val,'image')
                initwithimage=1;
                continue
            end
            if isempty(res)
                error('TIGRE:SIRT:InvalidInput','Invalid Init option')
            end
            % % % % % % % ERROR
        case 'initimg'
            if default
                continue
            end
            if exist('initwithimage','var')
                if isequal(size(val),geo.nVoxel')
                    res=single(val);
                else
                    error('TIGRE:SIRT:InvalidInput','Invalid image for initialization');
                end
            end
        case 'qualmeas'
            if default
                QualMeasOpts={};
            else
                if iscellstr(val)
                    QualMeasOpts=val;
                else
                    error('TIGRE:SIRT:InvalidInput','Invalid quality measurement parameters');
                end
            end
        case 'nonneg'
            if default
                nonneg=true;
            else
                nonneg=val;
            end
        case 'gpuids'
            if default
                gpuids = GpuIds();
            else
                gpuids = val;
            end
        case 'redundancy_weighting'
            if default
                redundancy_weights = true;
            else
                redundancy_weights = val;
            end
        case 'groundtruth'
            if default
                gt=nan;
            else
                gt=val;
            end
        otherwise
            error('TIGRE:SIRT:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in SIRT()']);
    end
end

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function f = extractColumns(f, colIdx)
%EXTRACTCOLUMNS   Extract columns (or rows) of an array-valued TRIGTECH.
%   G = EXTRACTCOLUMNS(F, COLIDX) extracts the columns specified by the row
%   vector COLIDX from the FUN F so that G = F(:, COLIDX). COLIDX need not be
%   increasing in order or unique, but must contain only integers in the range
%   [1, M], where F has M columns.
%
% See also MAT2CELL.

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

% Extract/re-order the columns from f.values, f.coeffs:
f.values = f.values(:,colIdx);
f.coeffs = f.coeffs(:,colIdx);
f.isReal = f.isReal(colIdx);

end