text:
function out = horzcat(varargin)
%HORZCAT   Horizontal concatenation.
%   [A B] horizontally concatenates the CHEBTECH objects A and B to form an
%   array-valued CHEBTECH. [A,B] does the same. Any number of CHEBTECH objects
%   can be concatenated within one pair of brackets. Vertical concatenation is
%   not supported.

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

% Remove empties:
empties = cellfun(@isempty, varargin);
if ( all(empties) )
    out = varargin{1};
    return
else
    varargin(empties) = [];
end

% Check that all objects are CHEBTECHs. If not we cannot concatenate.
if ( ~all( cellfun(@(f) isa(f, 'chebtech'), varargin) ) )
    error('CHEBFUN:CHEBTECH:horzcat:typeMismatch', ...
    'Incompatible concatenation. Ensure discretizations are of the same type.');
end

% Prolong each Chebtech to the same length:
n = max(cellfun(@length, varargin));
F = cellfun(@(f) prolong(f, n), varargin, 'UniformOutput', false);

% Extract the data and collate the an array-valued CHEBTECH:
out = varargin{1};

% Coeffs:
out.coeffs = cell2mat(cellfun(@(f) f.coeffs, F, 'UniformOutput', false));

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [kHz] = GHz2kHz(GHz)
% Convert frequency from gigahertz to kilohertz.
% Chad A. Greene 2012
kHz = GHz*1e+6;

author: Sable

repo: mcbench-benchmarks


text:
%kdtrans 'three dimensional distance transform '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros dtrans.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'First Input data object'
% OutputFile: od 'Distance Output', required: '3D Map of distances'
% OutputFile: ov 'Distance Vector ', required: '3D Map of distances (as vectors'
%
% Example: [od, ov] = kdtrans(i, {'i','';'od','';'ov',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% dtrans - three dimensional distance transform
%
%  DESCRIPTION
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1996-2003, Rainer Heintzmann,  All rights reserved.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [Population,fitness] = EnvironmentalSelection(Population,PopObj,OffObj,N)
% The environmental selection of TiGE_1

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    [FrontNo,~] = NDSort([PopObj;OffObj],N);
    fcv = Calculate_fcv(Population); % CV of hybrid population
    fitness = FrontNo' + fcv./(fcv+1);
    [~,index] = sort(fitness);
    Population = Population(index(1:N));
    fitness = fitness(index(1:N));
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function plot_ideal_deconv5(rmsd,msdstd,msdb,biasmean,meanest,min95est,max95est,INFO,hrf,snr,TR,varargin)
% :Usage:
% ::
%
%     function plot_ideal_deconv5(rmsd,msdstd,msdb,biasmean,meanest,min95est,max95est,INFO,hrf,snr,TR,[truer])
%
% delta should be matrix of column vectors
%
% optional: truer, vector of "true" responses for each trial type
%
% ..
%    Tor Wager
% ..

msdblim = max(max(abs(msdb)));
bmlim = max(max(abs(biasmean))); bmlim = [-bmlim bmlim];

if isfield(INFO,'ideal_data') & isfield(INFO,'delta')
    figure('Color','w'), myleg{1} = 'Ideal response';
    h(1) = plot(INFO.ideal_data,'LineWidth',2);, hold on;
    mycols = {'r' 'g' 'y' 'c' 'm'}; mycols = repmat(mycols,1,10);
    delta = INFO.delta;
    for i = 1:size(delta,2)
        tmp1 = find(delta(:,i));
        hh = plot([tmp1 tmp1]',[-1*ones(size(tmp1)) zeros(size(tmp1))]',mycols{i},'LineWidth',2);
        h(i+1) = hh(1);
        myleg{i+1} = ['Trial type ' num2str(i) ' onsets'];
    end
    set(gca,'FontSize',18)
    legend(myleg)
    set(gcf,'Position',[19         606        1561         360])

    %tmp = conv(hrf,delta(:,1));figure;plot(tmp),title('ttype1')
    %tmp = conv(hrf,sum(delta,2));figure;plot(tmp),title('alltypes')
end
%figure;plot([find(tmp);find(tmp)],[zeros(size(find(tmp)),1); ones(size(find(tmp)),1)])

if length(varargin) > 0
    truer = varargin{1};
else
    truer = ones(1,size(INFO.delta,2));
end

figure; subplot 221; hold on; set(gcf,'Color','w');
if ~isempty(INFO.autocorrelation)
    figure('Color','w'); set(gca,'FontSize',18)
    plot((1:length(INFO.autocorrelation))./(length(INFO.autocorrelation) / length(hrf)),INFO.autocorrelation,'k-','LineWidth',2)
end
title('Noise autocorrelation function')

% old plot of RMSD
% -----------------
%plot(snr,rmsd,'ko-','MarkerFaceColor','k','MarkerSize',3,'LineWidth',2)
%plot(snr,msdstd(:,1),'k--')
%plot(snr,msdstd(:,2),'k--')
%legend({'Mean population value' '95% CI for one regression'})
%ylabel('Root Mean Squared Deviation')
%xlabel('Signal to Noise Ratio')

%title(['Est. - True response: ' INFO.condition_tested],'FontSize',14)

subplot 222; hold on;
x = (1:length(hrf)) * TR;
plot(x, hrf.*truer(INFO.ttype),'k','LineWidth',2)
myleg = {'HRF'};
legend(myleg)
title(['Est. vs. True response: ' INFO.condition_tested],'FontSize',14)
xlabel('Time from trial onset')

mycolors = {'r' 'b' 'g' 'c' 'm'};
for i = 1:min(5,length(snr))
    plot(x(1:size(meanest,2)),meanest(i,:),mycolors{i})
    myleg{end+1} = ['Estimated, snr = ' num2str(snr(i))];
end

if ~isempty(min95est) & ~isempty(max95est)
    plot(x(1:size(min95est,2)),min95est(1,:),'r--');
    plot(x(1:size(min95est,2)),max95est(1,:),'r--');
    myleg{end+1} = ['95% CI for ind. regressions, snr = ' num2str(snr(1))];
end
legend(myleg)

subplot 223; hold on;
imagesc(msdb,[0 msdblim]), colormap copper
xlabel('Time from trial onset')
ylabel('Signal to Noise Ratio')
%set(gca,'YDir','Reverse')
set(gca,'YTick',1:length(snr))
set(gca,'YTickLabel',snr)
set(gca,'XTick',(1:2:INFO.estlength))
set(gca,'XTickLabel',get(gca,'XTick') .* TR)
title('Mean abs. dev. by time from event onset','FontSize',14)
colorbar('horiz')

subplot 224; hold on;
imagesc(biasmean,bmlim), colormap jet
xlabel('Time from trial onset')
ylabel('Signal to Noise Ratio')
set(gca,'YTick',1:length(snr))
set(gca,'YTickLabel',snr)
set(gca,'XTick',(1:2:INFO.estlength))
set(gca,'XTickLabel',get(gca,'XTick') .* TR)
title('Bias in estimation by time from event onset','FontSize',14)
colorbar('horiz')

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [ y, j ] = j_carry_islamic ( y, j )

%*****************************************************************************80
%
%% J_CARRY_ISLAMIC carries year-days to years in an Islamic date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, J, a YJ date.
%
  while ( 1 )

    days = year_length_islamic ( y );

    if ( j < days )
      break
    end

    j = j - days;
    y = y + 1;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function add_io(hRD,design)

% add AXI4 and AXI4-Lite slave interfaces
hRD.addAXI4SlaveInterface( ...
    'InterfaceConnection', 'axi_cpu_interconnect/M05_AXI', ...
    'BaseAddress',         '0x43C00000', ...
    'MasterAddressSpace',  'sys_ps7/Data');

if ~strcmp(design,'Tx')
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Rx Reference design interfaces
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    hRD.addInternalIOInterface( ...
        'InterfaceID',    'IP Data Valid OUT', ...
        'InterfaceType',  'OUT', ...
        'PortName',       'dut_data_valid', ...
        'PortWidth',      1, ...
        'InterfaceConnection', 'util_ad9361_adc_pack/adc_valid_0', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'IP Data 0 OUT', ...
        'InterfaceType',  'OUT', ...
        'PortName',       'dut_data_0', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_adc_pack/adc_data_0', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'IP Data 1 OUT', ...
        'InterfaceType',  'OUT', ...
        'PortName',       'dut_data_1', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_adc_pack/adc_data_1', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'IP Data 2 OUT', ...
        'InterfaceType',  'OUT', ...
        'PortName',       'dut_data_2', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_adc_pack/adc_data_2', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'IP Data 3 OUT', ...
        'InterfaceType',  'OUT', ...
        'PortName',       'dut_data_3', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_adc_pack/adc_data_3', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'AD9361 ADC Data I0', ...
        'InterfaceType',  'IN', ...
        'PortName',       'sys_wfifo_0_dma_wdata', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_adc_fifo/dout_data_0', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'AD9361 ADC Data Q0', ...
        'InterfaceType',  'IN', ...
        'PortName',       'sys_wfifo_1_dma_wdata', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_adc_fifo/dout_data_1', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'AD9361 ADC Data I1', ...
        'InterfaceType',  'IN', ...
        'PortName',       'sys_wfifo_2_dma_wdata', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_adc_fifo/dout_data_2', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'AD9361 ADC Data Q1', ...
        'InterfaceType',  'IN', ...
        'PortName',       'sys_wfifo_3_dma_wdata', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_adc_fifo/dout_data_3', ...
        'IsRequired',     false);
end

if ~strcmp(design,'Rx')
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Tx Reference design interfaces
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    hRD.addInternalIOInterface( ...
        'InterfaceID',    'AD9361 DAC Data I0', ...
        'InterfaceType',  'OUT', ...
        'PortName',       'axi_ad9361_dac_data_i0', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'axi_ad9361_dac_fifo/din_data_0', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'AD9361 DAC Data Q0', ...
        'InterfaceType',  'OUT', ...
        'PortName',       'axi_ad9361_dac_data_q0', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'axi_ad9361_dac_fifo/din_data_1', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'AD9361 DAC Data I1', ...
        'InterfaceType',  'OUT', ...
        'PortName',       'axi_ad9361_dac_data_i1', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'axi_ad9361_dac_fifo/din_data_2', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'AD9361 DAC Data Q1', ...
        'InterfaceType',  'OUT', ...
        'PortName',       'axi_ad9361_dac_data_q1', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'axi_ad9361_dac_fifo/din_data_3', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'IP Load Tx Data OUT', ...
        'InterfaceType',  'OUT', ...
        'PortName',       'util_dac_unpack_dac_valid_00', ...
        'PortWidth',      1, ...
        'InterfaceConnection', 'axi_ad9361_dac_fifo/din_valid_in_0', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'IP Data 0 IN', ...
        'InterfaceType',  'IN', ...
        'PortName',       'util_dac_unpack_dac_data_00', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_dac_upack/dac_data_0', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'IP Data 1 IN', ...
        'InterfaceType',  'IN', ...
        'PortName',       'util_dac_unpack_dac_data_01', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_dac_upack/dac_data_1', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'IP Data 2 IN', ...
        'InterfaceType',  'IN', ...
        'PortName',       'util_dac_unpack_dac_data_02', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_dac_upack/dac_data_2', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'IP Data 3 IN', ...
        'InterfaceType',  'IN', ...
        'PortName',       'util_dac_unpack_dac_data_03', ...
        'PortWidth',      16, ...
        'InterfaceConnection', 'util_ad9361_dac_upack/dac_data_3', ...
        'IsRequired',     false);

    hRD.addInternalIOInterface( ...
        'InterfaceID',    'IP Valid Tx Data IN', ...
        'InterfaceType',  'IN', ...
        'PortName',       'util_dac_unpack_upack_valid_00', ...
        'PortWidth',      1, ...
        'InterfaceConnection', 'util_ad9361_dac_upack/dac_valid_out_0', ...
        'IsRequired',     false);

end

author: analogdevicesinc


text:
function tbxStruct  = demos

% DEMOS Demo list for the ANN Library.

% Version 1.1
% Giampiero Campa, West Virginia University
% 1-June-2007

if nargout == 0, demo blockset; return; end

tbxStruct.Name='Adaptive Neural Networks';
tbxStruct.Type='Blockset';

tbxStruct.Help={

   'The Adaptive Neural Networks Library (BANN) consists in a collection'
   'of blocks that implement several adaptive neural networks in simulink.'
};

 tbxStruct.DemoList = { ' Adaptive Neural Network Library',       'ann';
                        ' Main Demo File ',          'anndemo' };


text:
function value = sweet4_condition ( )

%*****************************************************************************80
%
%% SWEET4_CONDITION returns the L1 condition of the SWEET4 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real VALUE, the L1 condition.
%
  a_norm = 100.3190000000000;
  b_norm = 0.510081684645161;
  value = a_norm * b_norm;

  return
end


text:
function vl_demo_print(varargin)
% VL_DEMO_PRINT
%   VL_DEMO_PRINT(NAME) prints the current figure to the documentation
%   directory with the specified filename, assuming that the global
%   variable VL_DEMO_PRINT is defined and non-empty when MATLAB is
%   started (our using SETENV() from MATLAB). Otherwise the function
%   flushes the displays and returns.
%
%   VL_DEMO_PRINT(NAME, R) specifies a magnification factor R, setting
%   the figure width relatively to the page width. If not specified, R
%   is assumed to be 1/2.
%
%   Remarks:: The figure paper type is set to letter, that has size 8.5 x
%     11 inches. When converted for web viewing, images are rasterized
%     at either 75 or 95 DPI, The documentation system converts images
%     to bitmap with a resolution of 75 DPI, which makes a letter size
%     page 637 or 808 pixels large, repsectively.
%
%     In MATLAB font sizes are usually expressed in points, where a
%     point is a 1/72 inch. Thus a 12pt font sampled at 75 DPI is
%     about 12.5 pixels high.

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% 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).

if isempty(getenv('VL_DEMO_PRINT'))
  drawnow ;
  return ;
end

if isa(varargin{1}, 'double')
  fig = varargin{1} ;
  varargin(1) = [] ;
else
  fig = gcf ;
end

name = varargin{1} ;

if length(varargin) < 2
  figurePaperSize = 0.5 ;
else
  figurePaperSize = varargin{2} ;
end

vl_printsize(fig, figurePaperSize) ;

figDir = fullfile(vl_root,'doc','demo') ;
if ~ exist(figDir, 'dir')
  mkdir(figDir) ;
end

if 0
  filePath = fullfile(figDir, [name '.eps']) ;
  print(fig, '-depsc2', filePath) ;
else
  filePath = fullfile(figDir, [name '.jpg']) ;
  print(fig, '-djpeg95', filePath, '-r95') ;
end
fprintf('%s: wrote file ''%s''\n', mfilename,  filePath) ;

author: jianxiongxiao

repo: ProfXkit


text:
function g = year_to_golden_number ( y )

%*****************************************************************************80
%
%% YEAR_TO_GOLDEN_NUMBER returns the golden number of a Common year.
%
%  Discussion:
%
%    Nineteen solar years are very close to 235 lunations.  Calendars
%    that try to keep track of both the sun and moon often make use of
%    this fact, ascribed to the Greek astronomer Meton.
%
%    While trying to determine a formula for Easter, Dionysus Exiguus
%    symbolized the place of each year in its Metonic cycle by a
%    "golden number" between 1 and 19.  The numbering began with the
%    year 1 BC, assigned the golden number of 1.  The following year,
%    1 AD, got the golden number of 2, and after that it gets easier.
%
%    The same golden year calculation is done for years in the Julian
%    or Gregorian calendar.
%
%  Example:
%
%    Year  Golden Number
%
%      1 BC     1
%      1 AD     2
%      2 AD     3
%     18 AD    19
%     19 AD     1
%     20 AD     2
%   1066 AD     3
%   1900 AD     1
%   1919 AD     1
%   1938 AD     1
%   1957 AD     1
%   1976 AD     1
%   1995 AD     1
%   2014 AD     1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year.
%
%    Output, integer G, the golden number, between 1 and 19.  This
%    records the position of the year in the 19 year Metonic cycle.
%
  if ( y == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'YEAR_TO_GOLDEN_NUMBER - Fatal error!\n' );
    fprintf ( 1, '  Illegal input Y = 0.\n' );
    error ( 'YEAR_TO_GOLDEN_NUMBER - Fatal error!' );
  end
%
%  We assume that BC years come in as negative numbers, and that
%  the year before 1 AD is 1 BC.  So add 1 to any negative value
%  so that the arithmetic works.
%
  y2 = y_common_to_astronomical ( y );

  g = i4_wrap ( y2 + 1, 1, 19 );

  return
end


text:
function [day] = s2day(s)
% Convert time from seconds to days.
% Chad Greene 2012
day = s*0.00001157407407407;


text:
%% Weighted optimization for CP tensor decomposition with incomplete data
% We explain how to use |cp_wopt| with the POBLANO toolbox.

%% Create an example problem with missing data.
% Here we have 25% missing data and 10% noise.
R = 2;
info = create_problem('Size', [15 10 5], 'Num_Factors', R, ...
    'M', 0.25, 'Noise', 0.10);
X = info.Data;
P = info.Pattern;
M_true= info.Soln;

%% Create initial guess using 'nvecs'
M_init = create_guess('Data', X, 'Num_Factors', R, ...
    'Factor_Generator', 'nvecs');

%% Set up the optimization parameters
% It's genearlly a good idea to consider the parameters of the optimization
% method. The default options may be either too stringent or not stringent
% enough. The most important options to consider are detailed here.

% Get the defaults
ncg_opts = ncg('defaults');
% Tighten the stop tolerance (norm of gradient). This is often too large.
ncg_opts.StopTol = 1.0e-6;
% Tighten relative change in function value tolearnce. This is often too large.
ncg_opts.RelFuncTol = 1.0e-20;
% Increase the number of iterations.
ncg_opts.MaxIters = 10^4;
% Only display every 10th iteration
ncg_opts.DisplayIters = 10;
% Display the final set of options
ncg_opts

%% Call the |cp_wopt| method
% Here is an example call to the cp_opt method. By default, each iteration
% prints the least squares fit function value (being minimized) and the
% norm of the gradient. The meaning of any line search warnings
% can be checked via <matlab:doc('cvsrch') doc cvsrch>.
[M,~,output] = cp_wopt(X, P, R, 'init', M_init, ...
    'alg', 'ncg', 'alg_options', ncg_opts);

%% Check the output
% It's important to check the output of the optimization method. In
% particular, it's worthwhile to check the exit flag.
% A zero (0) indicates successful termination with the gradient smaller
% than the specified StopTol, and a three (3) indicates a successful
% termination where the change in function value is less than RelFuncTol.
% The meaning of any other flags can be checked via
% <matlab:doc('poblano_params') doc poblano_params>.
exitflag = output.ExitFlag

%% Evaluate the output
% We can "score" the similarity of the model computed by CP and compare
% that with the truth. The |score| function on ktensor's gives a score in
% [0,1]  with 1 indicating a perfect match. Because we have noise, we do
% not expect the fit to be perfect. See <matlab:doc('ktensor/score') doc
% score> for more details.
scr = score(M,M_true)

%% Create a SPARSE example problem with missing data.
% Here we have 95% missing data and 10% noise.
R = 2;
info = create_problem('Size', [150 100 50], 'Num_Factors', R, ...
    'M', 0.95, 'Sparse_M', true, 'Noise', 0.10);
X = info.Data;
P = info.Pattern;
M_true= info.Soln;

%% Create initial guess using 'nvecs'
M_init = create_guess('Data', X, 'Num_Factors', R, ...
    'Factor_Generator', 'nvecs');

%% Set up the optimization parameters
% It's genearlly a good idea to consider the parameters of the optimization
% method. The default options may be either too stringent or not stringent
% enough. The most important options to consider are detailed here.

% Get the defaults
ncg_opts = ncg('defaults');
% Tighten the stop tolerance (norm of gradient). This is often too large.
ncg_opts.StopTol = 1.0e-6;
% Tighten relative change in function value tolearnce. This is often too large.
ncg_opts.RelFuncTol = 1.0e-20;
% Increase the number of iterations.
ncg_opts.MaxIters = 10^4;
% Only display every 10th iteration
ncg_opts.DisplayIters = 10;
% Display the final set of options
ncg_opts

%% Call the |cp_wopt| method
% Here is an example call to the cp_opt method. By default, each iteration
% prints the least squares fit function value (being minimized) and the
% norm of the gradient. The meaning of any line search warnings
% can be checked via <matlab:doc('cvsrch') doc cvsrch>.
[M,~,output] = cp_wopt(X, P, R, 'init', M_init, ...
    'alg', 'ncg', 'alg_options', ncg_opts);

%% Check the output
% It's important to check the output of the optimization method. In
% particular, it's worthwhile to check the exit flag.
% A zero (0) indicates successful termination with the gradient smaller
% than the specified StopTol, and a three (3) indicates a successful
% termination where the change in function value is less than RelFuncTol.
% The meaning of any other flags can be checked via
% <matlab:doc('poblano_params') doc poblano_params>.
exitflag = output.ExitFlag

%% Evaluate the output
% We can "score" the similarity of the model computed by CP and compare
% that with the truth. The |score| function on ktensor's gives a score in
% [0,1]  with 1 indicating a perfect match. Because we have noise, we do
% not expect the fit to be perfect. See <matlab:doc('ktensor/score') doc
% score> for more details.
scr = score(M,M_true)

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
%% Ex. 7 More complicated use of loop and index

b = [2 5 7 4 9 8 3];
c = [2 3 5 7];
sum1 = 0;
for k = 1:4
 sum1 = sum1+b(c(k));
end
sum1
%Output:
%       24
% Remark: This program performs the summation of
% sum1 = b(c(1))+b(c(2))+b(c(3))+b(c(4))
%      = b(2)+b(3)+b(5)+b(7)
%       = 5+7+9+3
%       = 24

author: TheAlgorithms

repo: MATLAB-Octave


text:
function [brainMask,checkSlices] = mrAnatExtractBrain(img, mmPerVox, betLevel, outFile)
%Uses the FSL BET tool to compute a brain mask for the given image volume.
%
% [brainMask,checkSlices] = mrAnatExtractBrain([img], [mmPerVox=[1,1,1]], [betLevel=0.5], [outFile])
%
% The betLevel is the BET fractional intensity threshold parameter. It
% should be [0-1]. The default of 0.5 usually works well. Smaller values
% will yield a larger brain estimate. If you get brain stuff chopped off,
% try a lower value. If you get non-brain stuff included, try a larger
% value.
%
% If you provide an array of betLevel values, you'll get a cell array of
% length(betLevel) brain masks back.
%
% checkSlices is a montage of every 3rd slice with the mask shown in blue
% and the underlying anatomy on yellow. Again, for length(betLevels)>1,
% you'll get a cell array.
%
% img can be a string, in which case it is assumed to be a NIFTI filename.
% or, it can be a NIFTI struct (as from niftiRead). In both cases,
% mmPerVox is ignored and instead gleaned from the NIFTI header.
%
% If no output arguments are captured, the brain mask is saved in the same
% place as the input file, but with '_mask' appended to the name.
%
% HISTORY:
% 2006.02.02 RFD: wrote it.
%
% Bob (c) Stanford VISTASOFT, 2006

if(~exist('betLevel','var') || isempty(betLevel)), betLevel = 0.5; end
if(~exist('mmPerVox','var') || isempty(mmPerVox)), mmPerVox = [1 1 1]; end

if(~exist('img','var') || isempty(img))
    [f,p] = uigetfile({'*.nii.gz';'*.*'},'Select a t1-weighted NIFTI file...');
    if(isnumeric(f)), disp('User canceled.'); return; end
    img = fullfile(p,f);
end

if(isstruct(img))
    % It's a nifti image
    ni = img;
    mmPerVox = ni.pixdim;
    img = ni.data;
elseif(ischar(img))
    % It's a nifti filename
    ni = niftiRead(img);
    mmPerVox = ni.pixdim;
elseif(nargout==0)
    error('If you pass in raw image data, you must capture at least one output!');
end

% TODO: more robust way to find cygwin on PC
if(ispc), bash = 'c:/cygwin/bin/bash.exe';
else      bash = 'bash';
end

% First try the shell for the preferrred FSL location, if it is not found
% then we'll test to see if the user is running OSX, if so then we'll point
% to an OSX version of bet2 (bet2_osx) else we'll point to the linux
% version.
[stat,res] = system('which bet2');
if(stat==0 && ~isempty(res))
    bet = res(1:end-1);
else
    if(ismac)
        bet = fullfile(fileparts(which(mfilename)), 'bet2_osx');
    else
        bet = fullfile(fileparts(which(mfilename)), 'bet2');
    end
end

bet = strtrim(bet); % The name of the bet funciton on the system

% Set temporary file names for the out file and bet script. These will be
% removed at the end of the function run.
out = [tempname, '_bet_temp'];
betScript = [tempname, '_bet_script.sh'];

if(ischar(img))
    out = img;
else
    img = double(img);
    mx = max(abs(img(:)));
    img = int16(round(img./mx.*32767));
    analyzeWrite(img, out, mmPerVox);   % Why not write a nifti?
end

% We could specify a better starting position for the BET surface
% sphere estimate, e.g., by using the talairach landmarks.
if(nargout==0)
    [p,f,e] = fileparts(ni.fname);
    [x,f,e2] = fileparts(f);
    imBaseName = fullfile(p,f);
end

%% Run brain extraction tool for each level
for ii=1:numel(betLevel)
    betOut = tempname;
    betCmd = [bet ' ' out ' ' betOut ' -mnf ' num2str(betLevel(ii))];
    fid = fopen(betScript,'wt');
    fprintf(fid,'#!/bin/bash\nexport FSLOUTPUTTYPE=NIFTI_GZ\n%s\n',betCmd);
    fclose(fid);
    unix([bash ' ' betScript]);
    betOut = [betOut '_mask.nii.gz'];
    if(nargout==0)
        if(numel(betLevel)>1)
            movefile(betOut, sprintf('%s_mask_f%02d',imBaseName,betLevel(ii)*100));
        else
            movefile(betOut, [imBaseName '_mask.nii.gz']);
        end
    else
        brainMaskFile{ii} = betOut;
    end
end

%% Manage outputs, but not sure about the logic here.
if(nargout>0)
    tmp = niftiRead(betOut);
    brainMask{ii} = logical(tmp.data);
    if(nargout>1)
        slices = [1:3:size(brainMask{ii},3)-6];
        if(ischar(img))
            im = double(makeMontage(ni.data, slices));
        else
            im = double(makeMontage(img, slices));
        end
        im = mrAnatHistogramClip(im, 0.4, 0.99);
        im = uint8(round(im./max(im(:)).*255));
        im(:,:,2) = im(:,:,1);
        im(:,:,3) = uint8(double(makeMontage(brainMask{ii},slices)).*255);
        checkSlices{ii} = im;
    end
end

if(nargout>0)
    if(length(betLevel)==1)
        brainMask = brainMask{1};
        if(exist('checkSlices','var')), checkSlices = checkSlices{1}; end
    end
else
    clear all;
end

%% Clean up

% Remove the betScript from the tempdir
delete(betScript)

% If img and out are the same then it is a temporary file and should be removed.
% If not, then it was passed in from a location on disk and should be preserved
if ~strcmp(img, out) && strfind(out, tempdir) && exist(out,'file')
    delete(out)
end

%% Return

return

%% Example script

bd = pwd;
d = dir('*0*');
subDirs = {d.name};
betLevel = 0.45;

for(ii=1:length(subDirs))
    t1 = fullfile(bd,subDirs{ii},'t1','t1');
    if(exist([t1 '.nii.gz'],'file'))
        t1Ni = niftiRead([t1 '.nii.gz']);
        [brainMask,checkSlices] = mrAnatExtractBrain(t1Ni.data, t1Ni.pixdim, betLevel);
        figure(99); image(checkSlices); axis image tight off;
        dtiWriteNiftiWrapper(uint8(brainMask), t1Ni.qto_xyz, [t1 '_mask.nii.gz']);
        ccPerPixel = prod(t1Ni.pixdim(1:3))/1000;
        brainVolumeCc(ii) = numel(find(brainMask)) * ccPerPixel;
    end
end

outFile = 'brainVolumes.txt';
fid = fopen(outFile,'wt');
for(ii=1:length(subDirs))
    %t1Mask = niftiRead(fullfile(bd,subDirs{ii},'t1','t1_mask.nii.gz'));
    %ccPerPixel = prod(t1Mask.pixdim(1:3))/1000;
    %brainVolumeCc(ii) = numel(find(t1Mask.data)) * ccPerPixel;
    fprintf(fid,'%s\t%0.1f\n',subDirs{ii},brainVolumeCc(ii));
end
fclose(fid);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function backtrack_binary_rc_test01 ( )

%*****************************************************************************80
%
%% BACKTRACK_BINARY_RC_TEST01 seeks binary powers that have a given sum.
%
%  Discussion:
%
%    We consider the binary powers 1, 2, 4, ... 2^(n-1).
%
%    We wish to select some of these powers, so that the sum is equal
%    to a given target value.  We are actually simply seeking the binary
%    representation of an integer.
%
%    A partial solution is acceptable if it is less than the target value.
%
%    We list the powers in descending order, so that the bactracking
%    procedure makes the most significant choices first, thus quickly
%    eliminating many unsuitable choices.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 December 2013
%
%  Author:
%
%    John Burkardt
%
  n = 8;
  test_num = 3;
  targets = [ 73, 299, -3 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BACKTRACK_BINARY_RC_TEST01\n' );
  fprintf ( 1, '  Use BACKBIN_RC to find the binary expansion of\n' );
  fprintf ( 1, '  an integer between 0 and 255.\n' );
  fprintf ( 1, '  The choices are 0/1 for the 8 digits.\n' );

  for test = 1 : test_num

    target = targets(test);
    fprintf ( 1, '\n' );
    fprintf ( 1, '  TARGET = %d\n', target );
    call_num = 0;
    reject = 0;
    n2 = -1;
    choice = [];

    while ( 1 )

      [ n2, choice ] = backbin_rc ( n, reject, n2, choice );
      call_num = call_num + 1;

      if ( n2 == -1 )
        fprintf ( 1, '  Termination without solution.\n' );
        break
      end
%
%  Evaluate the integer determined by the choices.
%
      factor = 1;
      for i = n : -1 : n2 + 1
        factor = factor * 2;
      end

      result = 0;
      for i = 1 : n2
        result = result * 2 + choice(i);
      end

      result = result * factor;
%
%  If the integer is too big, then we reject it, and
%  all the related integers formed by making additional choices.
%
      reject = ( target < result );
%
%  If we hit the target, then in this case, we can exit because
%  the solution is unique.
%
      if ( result == target )
        break
      end

    end

    fprintf ( 1, '  Number of calls = %d\n', call_num );
    fprintf ( 1, '  Binary search space = %d\n', 2 ^ n );
    fprintf ( 1, '  ' );
    for i = 1 : n
      fprintf ( 1, '%2d', choice(i) );
    end
    fprintf ( 1, '\n' );

  end

  return
end


text:
function [jump_list, lid_ko, data,  running_mean] = getJumps(data, rate_scale)
% Return outliers in a relatively constant time series
%
% INPUT
%   data          data (column, eventually it can have a second column containing std)
%   rate_scale    24/rate_scale => used for scaling window size (1 => daily, 24 => hourly ...)
%
% OUTPUT
%
% SYNTAX
%   [jump_list, id_ko, data, running_mean] = getJumps(data, rate_scale)
%

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:        Andrea Gatti
%  Contributors:      Andrea Gatti
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------
    flag_smooth = true; % Apply a spline smoothing to rate (duration 3 days)
    spline_base = 7; % in days

    if size(data,2) == 2
        data_std = data(:,2);
        data = data(:,1);
    else
        data_std = 0;
    end
    med_data = median(data, 'omitnan');
    data = data - med_data;

    % Compute thresholds - empirical but it works (outliers should be less than 75% percent)
    thr_perc = 0.75;
    tmp = diff(data);
    tmp(abs(tmp) > 3 * perc(noNaN(abs(tmp)), thr_perc)) = nan;
    thr = std(tmp, 'omitnan');
    clear tmp;

    % Outlier detection
    if any(data_std)
        [lid_ko, data] = getOutliers([data data_std], thr);
    else
        [lid_ko, data] = getOutliers(data, thr);
    end

    flag_running = true;
    % Initialization

    % Jump parameters
    run_win = 5 * rate_scale;   % window size for running mean
    cumsum_thr = 2;             % above this threshold, cumulate residuals
    cum_thr = 6;                % above this threshold the cumulate TRIGGERS A JUMP

    % Recompute a running_threshold
    % when the noise is higher the discrimination of jumps is lower
    thr = movstd(diff(data), run_win, 'omitnan');
    if not(any(thr))
        thr(:) = nan2zero(std(data, 'omitnan'));
    end

    q_win = floor(run_win / 2);
    thr((q_win + 1) : end - q_win) = min(thr((2*q_win +1):end), thr(1:end - (2*q_win)));

    mthr = median(thr, 'omitnan');

    thr(isnan(zero2nan(thr))) = mthr;
    run_win_variable = round(max(1,sqrt((thr - mthr)/mthr)) * run_win);
    run_win_variable = [run_win_variable; run_win_variable(end)];
    %thr = mthr + 0.1 * (thr - mthr);

    mu = data(1); % initial mean
    i = 1 + find(not(isnan(data(2:end))), 1, 'first');  % running index on observations
    [s_plus, s_minus] = deal(zeros(size(data,1),1)); % sensors
    [n_plus, n_minus] = deal(uint16(zeros(size(data,1),1))); % epochs since sensor cumulation started

    running_mean = data; % running mean

    j = 1;           % running index on jumps
    jump_list = 0;   % list of jumps
    while i <= length(data)
        % if jumps are not detected
        if isnan(data(i))
            s_plus(i) = s_plus(i-1);        % keep the old sensor
            s_minus(i) = s_minus(i-1);      % keep the old sensor
            n_plus(i)  = iif(s_plus(i)  > 0, n_plus(i-1)  + 1, 0);
            n_minus(i) = iif(s_minus(i) < 0, n_minus(i-1) + 1, 0);
            running_mean(i) = mu;
            i = i + 1;
        else
            % while no jumps are detected and i is running on observations
            while i <= length(data) && s_plus(i-1) <= (cum_thr * thr(i-1)) && s_minus(i-1) >= -(cum_thr * thr(i-1))
                val = iif(abs(data(i) - mu) > 0.9*(cum_thr * thr(i-1)), mu, data(i));
                s_plus(i)  = max(0, s_plus(i-1)  + (data(i) - mu) - (cumsum_thr * thr(i-1)));
                s_minus(i) = min(0, s_minus(i-1) + (data(i) - mu) + (cumsum_thr * thr(i-1)));

                n_plus(i)  = iif(s_plus(i)  > 0 , n_plus(i-1)  + 1, 0);
                n_minus(i) = iif(s_minus(i) < 0 , n_minus(i-1) + 1, 0);

                if flag_running
                    mu = nan2zero(mean(data(max(jump_list(j)+1, i + 1 - run_win_variable(i)):i), 'omitnan'));
                end
                running_mean(i) = mu;
                i = i + 1;
            end
            if i <= length(data)
                n_d = i - 1;
                j = j + 1;
                jump_list(j) = iif(s_plus(n_d) > (cum_thr * thr(i-1)), n_d - n_plus(n_d), n_d - n_minus(i-1));
                mu = nan2zero(mean(data(max(jump_list(j)+1, n_d - run_win_variable(n_d)):(n_d)), 'omitnan'));

                jmp_magnitude = running_mean(jump_list(j)) - mu;
                Core.getLogger.addMessage(sprintf('Jump detected at epoch %d, it happened in epoch %d, length: %d, old mean %f.2, new mean %f.2, for a jump of %.2f\n', n_d, jump_list(j), n_d - jump_list(j) + 1, running_mean(jump_list(j)), mu, jmp_magnitude),100);
                if (i <= length(data))
                    [s_plus(i-1), s_minus(i-1)] = deal(0);
                    [n_plus(i-1), n_minus(i-1)] = deal(0);
                end
            end
        end
        mu = nan2zero(mean(data(max(jump_list(j)+1, i - run_win_variable(i-1)):(i-1)), 'omitnan'));
        running_mean(i-1) = mu;
    end

    % compute running mean backwords
    if flag_running
        running_smooth = running_mean;
        robustness_perc = 0.90;
        n_sigma = 6;

        % flag small jumps (1-2 epochs are too small to be considered "real jumps"
        for j = fliplr(find(diff(jump_list) < 3))
            id_ko = jump_list(j)+1:jump_list(j+1);
            data(id_ko) = nan;
            lid_ko(id_ko) = true;
            jump_list(j) = jump_list(j+1);
            jump_list(j+1) = [];
        end

        last_win_id = jump_list(j)+1 : size(data,1); % last window
        % Compute a trend "robust" using the robustness_perc of data
        %step_median(last_win_id) = median(data(last_win_id), 'omitnan');
        [tmp, trend] = strongDeTrend(data(last_win_id), robustness_perc, 1-((1-robustness_perc)/2), n_sigma);
        if flag_smooth && sum(not(isnan((data(last_win_id))))) > 5
            if any(trend); running_mean(last_win_id) = movmean(tmp, run_win, 'omitnan') + trend; end
            try
                [~, ~, ~, running_smooth(last_win_id)] = splinerMat(1:numel(last_win_id), data(last_win_id), ceil(spline_base * rate_scale), 1e-5, 1:numel(last_win_id)); % medium splines
            catch
                if any(trend); running_smooth(last_win_id) = movmean(tmp, run_win, 'omitnan') + trend; end
            end
        else
            if any(trend); running_smooth(last_win_id) = movmean(tmp, run_win, 'omitnan') + trend; end
        end
        for j = 2 : numel(jump_list)
            last_win_id = jump_list(j-1)+1 : (jump_list(j));
            %step_median(last_win_id) = median(data(last_win_id), 'omitnan');
            [tmp, trend] = strongDeTrend(data(last_win_id), robustness_perc, 1-((1-robustness_perc)/2), n_sigma);
            if flag_smooth && sum(not(isnan((data(last_win_id))))) > 5
                if any(trend); running_mean(last_win_id) = movmean(tmp, run_win, 'omitnan') + trend; end
                try
                    [~, ~, ~, running_smooth(last_win_id)] = splinerMat(1:numel(last_win_id), data(last_win_id), ceil(spline_base * rate_scale), 1e-5, 1:numel(last_win_id)); % medium splines
                catch
                    if any(trend); running_smooth(last_win_id) = movmean(tmp, run_win, 'omitnan') + trend; end
                end
            else
                if any(trend); running_smooth(last_win_id) = movmean(tmp, run_win, 'omitnan') + trend; end
            end
        end

        %  remove very small jumps - recursive
        flag_rem_jmp = true;
        while flag_rem_jmp
            % Check the magnitude of the jump, with both moving mean and splines, if it is too small is a false alarm
            jmp_magnitude = min(diff(running_mean(jump_list(2:end) + repmat([0 1]', 1, numel(jump_list)-1))), ...
                                diff(running_smooth(jump_list(2:end) + repmat([0 1]', 1, numel(jump_list)-1))));
            % Jump smaller than half the noise in that window of data will be ignored
            id_small_jmp = abs(jmp_magnitude)' < max(0.4,(thr(jump_list(2:end))));
            if any(id_small_jmp)
                jump_list(find(id_small_jmp) + 1) = [];  % remove jmp

                % recompute running mean
                last_win_id = jump_list(end)+1 : size(data,1);
                %step_median(last_win_id) = median(data(last_win_id), 'omitnan');
                % remove the trend before computing the running mean (this help in case of very steep intervals of data)
                [tmp, trend] = strongDeTrend(data(last_win_id), robustness_perc, 1-((1-robustness_perc)/2), n_sigma);
                if flag_smooth && sum(not(isnan((data(last_win_id))))) > 5
                    if any(trend); running_mean(last_win_id) = movmean(tmp, run_win, 'omitnan') + trend; end
                    try
                        [~, ~, ~, running_smooth(last_win_id)] = splinerMat(1:numel(last_win_id), data(last_win_id), ceil(spline_base * rate_scale), 1e-5, 1:numel(last_win_id)); % medium splines
                    catch
                        running_smooth(last_win_id) = movmean(tmp, run_win, 'omitnan') + trend;
                    end
                else
                    if any(trend); running_smooth(last_win_id) = movmean(tmp, run_win, 'omitnan') + trend; end
                end
                for j = 2 : numel(jump_list)
                    last_win_id = jump_list(j-1)+1 : (jump_list(j));
                    [tmp, trend] = strongDeTrend(data(last_win_id), robustness_perc, 1-((1-robustness_perc)/2), n_sigma);
                    if flag_smooth && sum(not(isnan((data(last_win_id))))) > 5
                        if any(trend); running_mean(last_win_id) = movmean(tmp, run_win, 'omitnan') + trend; end
                        try
                            [~, ~, ~, running_smooth(last_win_id)] = splinerMat(1:numel(last_win_id), data(last_win_id), ceil(spline_base * rate_scale), 1e-5, 1:numel(last_win_id)); % medium splines
                        catch
                            if any(trend); running_smooth(last_win_id) = movmean(tmp, run_win, 'omitnan') + trend; end
                        end
                    else
                        if any(trend); running_smooth(last_win_id) = movmean(tmp, run_win, 'omitnan') + trend; end
                    end
                end
            else
                flag_rem_jmp = false;
            end
        end
        running_mean = running_smooth;
    else
        last_win_id = jump_list(j)+1 : size(data,1);
        running_mean(last_win_id) = mean(data(last_win_id), 'omitnan');
    end

    running_mean = running_mean + med_data;
    data = data + med_data;
    % figure; plot(-s_minus); hold on; plot(s_plus); hold on; plot(cum_thr * thr)
end

author: goGPS-Project


text:
function YESNO = factored(X)

YESNO = ~isempty(X.midfactors);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [a,c]=HermiteInterpMultiPoly(x,y,numDims)
%%HERMITEINTERPMULTIPOLYS Given a multidimensional function y(x) evaluated
%             a certain values of x (control points) along with an
%             arbitrary number of derivatives of y(x) at those points, this
%             function returns the coefficients of the Hermite
%             interpolating polynomials fitting those points across all
%             dimensions. This function differs from
%             findHermiteInterpPolySet in that this function tries to find
%             one high-order Hermite interpolation polynomial for each
%             dimension that matches all of the points given, whereas
%             findHermiteInterpPolySet finds a set of low-order
%             polynomials, which is more practical for large regions as
%             finite precision errors can dominate this function. Compare
%             this to the function HermiteInterpPoly, which can only handle
%             scalar values.
%
%INPUTS: x An NpX1 or 1XNp vector of real, scalar values at which the
%          function y and its derivative are given. The values are assumed
%          to be given in increasing order. Np>=2.
%        y A (numMoments*numDims)XNp matrix of values of the
%          multidimensional function y(x) and its derivatives evaluated at
%          the points in x. All values for a particular derivative are
%          given before any of the next derivative. For example, to
%          interpolate in 3D position and velocity, the ordering would be
%          [x;y;z;xDot;yDot;zDot].
%  numDims The scalar number of dimensions present. This will typically be
%          2 or 3 when dealing with target states. If this parameter is
%          omitted or an empty matrix is passed, then numDims=1 is used.
%
%OUTPUTS: a The numDimsXnumCoeff matrix of the interpolating polynomial
%           coefficients for each dimension and moment. This and the
%           next output can be used in the function polyValNewton to
%           interpolate values.
%         c The numDimsX(numCoeff-1) matrix of control points for the
%           interpolating polynomials.
%
%This function calls HermiteInterpPoly for each of the dimensions present
%and stacks the results by dimension.
%
%May 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(numDims))
    numDims=1;
end

if(numDims==1)
    [a,c]=HermiteInterpPoly(x,y);
    return
end

Np=length(x);
stateSize=size(y,1);

numMatch=stateSize/numDims;
numCoeff=numMatch*Np;

a=zeros(numDims,numCoeff);
c=zeros(numDims,numCoeff-1);

for curDim=1:numDims
    yCur=y(curDim:numDims:stateSize,:);
    [a(curDim,:),c(curDim,:)]=HermiteInterpPoly(x,yCur);
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
classdef Density < DesignVariable

    properties (Access = private)
        creatorSettings
        initCase
    end

    methods (Access = public)

        function obj = Density(cParams)
            obj.nVariables = 1;
            obj.init(cParams);
            obj.initCase = cParams.initialCase;
            obj.creatorSettings  = cParams.creatorSettings;
            obj.createValue();
        end

        function v = getVariablesToPlot(obj)
            v{1} = obj.value;
        end

        function [fun, funNames] = getFunsToPlot(obj)
            aa.mesh = obj.mesh;
            aa.fValues = obj.value;
            valFun = P1Function(aa);

            fun = {valFun};
            funNames = {'Density'};
        end

        function rho = computeVolumeFraction(obj)
            s.mesh   = obj.mesh;
            s.fValues = obj.value;
            f = P1Function(s);
            q = Quadrature.set(obj.mesh.type);
            q.computeQuadrature('CONSTANT');
            xV = q.posgp;
            rho = f.evaluate(xV);
        end

    end

    methods (Access = private)

        function createValue(obj)
            s = obj.creatorSettings;
            switch s.type
                case 'FromLevelSet'
                    s.ndim  = obj.mesh.ndim;
                    s.coord = obj.mesh.coord;
                    s.type  = obj.initCase;
                    lsCreator  = LevelSetCreator.create(s);
                    phi        = lsCreator.getValue();
                    obj.value  = 1 - heaviside(phi);
                case 'Given'
                    obj.value = s.rho0.*ones(size(obj.mesh.coord,1),1);
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function Y = vl_nnreshape(X,dest_size,dzdy)
% VL_NNRESHAPE CNN reshapes input to des_size
%
%    DEST_SIZE: desired size
%
%    DZDX = VL_NNRESHAPE(X, DEST_SIZE, DZDY) computes the derivative DZDX
%    of the CNN with respect to the input X given the derivative DZDY
%    with respect to the block output Y. DZDX has the same dimension
%    as X.

% Copyright (C) 2015 Tuan-Hung VU.
% All rights reserved.
%
% This file is made available under the terms of the BSD license (see the COPYING file).

if nargin <= 2
    sz = [size(X,1) size(X,2) size(X,3) size(X,4)] ;
    Y = reshape(X, [size(X,1) size(X,2) dest_size(2) size(X,3)*size(X,4)/dest_size(2)]);
else
    Y = reshape(dzdy, size(X));
end

author: aosokin


text:
classdef TestUtils
    %TestUtils

    methods (Static)
        function test_version
            info = cv.Utils.getBuildInformation();
            validateattributes(info, {'char'}, {'nonempty'});

            v = cv.Utils.version();
            validateattributes(v, {'char'}, {'row', 'nonempty'});
        end

        function test_cpu
            s = cv.Utils.checkHardwareSupport();
            validateattributes(s, {'struct'}, {'scalar'});

            name = cv.Utils.getHardwareFeatureName(int32(1));
            validateattributes(name, {'char'}, {'row', 'nonempty'});
        end

        function test_num_cpu
            n = cv.Utils.getNumberOfCPUs();
            validateattributes(n, {'numeric'}, {'scalar', 'integer', 'nonnegative'});
        end

        function test_num_threads
            n = cv.Utils.getNumThreads();
            validateattributes(n, {'numeric'}, {'scalar', 'integer'});
            cv.Utils.setNumThreads(n);
        end

        function test_optimization
            b = cv.Utils.useOptimized();
            validateattributes(b, {'logical'}, {'scalar'});
            cv.Utils.setUseOptimized(b);
        end

        function test_ipp
            str = cv.Utils.getIppVersion();
            validateattributes(str, {'char'}, {'row', 'nonempty'});

            if ~strcmpi(str, 'disabled')
                b = cv.Utils.useIPP();
                validateattributes(b, {'logical'}, {'scalar'});
                cv.Utils.setUseIPP(b);
            end
        end

        function test_ovx
            b = cv.Utils.haveOpenVX();
            validateattributes(b, {'logical'}, {'scalar'});

            if b
                b = cv.Utils.useOpenVX();
                validateattributes(b, {'logical'}, {'scalar'});
                cv.Utils.setUseOpenVX(b);
            end
        end

        function test_ocl_1
            b = cv.Utils.haveOpenCL();
            validateattributes(b, {'logical'}, {'scalar'});

            b = cv.Utils.haveAmdBlas();
            validateattributes(b, {'logical'}, {'scalar'});

            b = cv.Utils.haveAmdFft();
            validateattributes(b, {'logical'}, {'scalar'});

            b = cv.Utils.haveSVM();
            validateattributes(b, {'logical'}, {'scalar'});
        end

        function test_ocl_2
            if cv.Utils.haveOpenCL()
                b = cv.Utils.useOpenCL();
                validateattributes(b, {'logical'}, {'scalar'});
                cv.Utils.setUseOpenCL(b);
            end
        end

        function test_ocl_3
            if cv.Utils.haveOpenCL()
                p = cv.Utils.getPlatfomsInfo();
                validateattributes(p, {'struct'}, {'vector'});
            end
        end

        function test_cuda
            n = cv.Utils.getCudaEnabledDeviceCount();
            validateattributes(n, {'numeric'}, {'scalar', 'integer'});
            if n > 0
                id = cv.Utils.getDevice();
                validateattributes(id, {'numeric'}, {'scalar', 'integer'});
                cv.Utils.setDevice(id);

                s = cv.Utils.deviceSupports();
                validateattributes(s, {'struct'}, {'scalar'});

                cv.Utils.printShortCudaDeviceInfo(id);
                cv.Utils.printCudaDeviceInfo(id);

                s = cv.Utils.deviceInfo(id);
                validateattributes(s, {'struct'}, {'scalar'});

                cv.Utils.resetDevice();
            end
        end

        function test_tega
            % functions only defined for Tegra SoC
            if false
                b = cv.Utils.useTegra();
                validateattributes(b, {'logical'}, {'scalar'});
                cv.Utils.setUseTegra(b);
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
clear all; close all;
I=imread('rice.png');
J=im2double(I);
T=dctmtx(8);
K=blkproc(J, [8 8], 'P1*x*P2', T, T');
mask=[ 1  1  1  1  0  0  0  0
            1  1  1  0  0  0  0  0
            1  1  0  0  0  0  0  0
            1  0  0  0  0  0  0  0
            0  0  0  0  0  0  0  0
            0  0  0  0  0  0  0  0
            0  0  0  0  0  0  0  0
            0  0  0  0  0  0  0  0 ];
K2=blkproc(K, [8 8], 'P1.*x', mask);
L=blkproc(K2, [8 8], 'P1*x*P2', T', T);
figure;
subplot(121);
imshow(J);
subplot(122);
imshow(L);

author: HuangCongQing


text:
function FR = factorialratio(A,B)
% FACTORIALRATIO - ratio of factorial numbers
%
%   FR = FACTORIALRATIO (A,B) calculates the ratio between products of
%   factorials specified by A and B.  A and B are lists of positive
%   integers, with N and M elements, respectively. The number of values in
%   each list may differ. The factorial ratio is based on the following formula:
%
%              (A1! * A2! * ... * AN!)
%      FR  =  -------------------------
%              (B1! * B2! * ... * BM!)
%
%    However, an approach is taken that increases the accuracy when A or B
%    contain larger values, as mentioned in the help of FACTORIAL.
%
%    Examples:
%
%      f1 = factorialratio([10 8 6],[9 6 5]) ;
%      f2 = 10 * 8 * 7 * 6 ; % by hand
%      f3 = (factorial(10) * factorial(8) * factorial(6)) ./ ...
%            (factorial(9) * factorial(6) * factorial(5)) ;
%      disp('f1, f2, f3 = ') ; disp([f1 f2 f3])
%
%
%      % Also for larger values FACTORIALRATIO works correctly
%      f1 = factorialratio(150,143) ;
%      f2 = prod(144:150) ; % by hand
%      % where FACTORIAL lacks precision
%      f3 = factorial(150) ./ factorial(143) ;
%      disp('f1, f2, f3 = ') ; disp([f1 f2 f3])
%      disp('   Differences with f2 = ') ; disp([f1-f2 f2-f2 f3-f2])
%
%      % and even for extreme numbers things turn out well
%      f1 = factorialratio(30000,29998) ;
%      f2 = (30000 * 29999) ; % by hand
%      % whereas calls to FACTORIAL fail miserably ...
%      f3 = factorial(30000) ./ factorial(29998) ;
%      disp('f1, f2, f3 = ') ; disp([f1 f2 f3])
%
%    Note that, for instance, factorialratio(200,1) also fails, of course.
%    For such occasions other engines have to be used.
%
%    See also PROD, FACTORIAL, GAMMALN

% for Matlab R13 and up
% version 2.0 (oct 2012)
% (c) Jos van der Geest
% email: jos@jasen.nl

% History
% 1.0 (feb 2009) - created for (later) use in a statistical test
% 1.1 (feb 2009) - fixed minor bug in error checking
% 2.0 (oct 2012) - fixed bug "maxE > 1" which should be "maxE > 0". Thanks
%                  to 'Charles' for pointing this out.

% check the input arguments
error(nargchk(2,2,nargin)) ;
if isempty(A) || isempty(B) || ...
    ~isnumeric(A) || ~isnumeric(B) || ...
    any(A(:)<0) || any(fix(A(:)) ~= A(:)) || ...
    any(B(:)<0) || any(fix(B(:)) ~= B(:))
    error('Inputs should be numerical arrays with positive integers.')
end

% Example for terminology used in comments:
%   FR = factorialratio([6 4],[3 2 1])      call
%      = (6! * 4!) / (3! * 2!)              numerators / denominators  (ignoring ones)
%      = 6*5*4*3*2*4*3*2 / 3*2*2            expanded factorial (showing all factors > 1)
%         [6 5 4 3 2] .^ [1 1 1 2 2]
%      = ----------------------------       as powers (factors and exponents)
%         [6 5 4 3 2] .^ [0 0 0 1 2]
%
%      = [6 5 4 3 2] .^ [1 1 1 1 0]         rewritten

% only values larger than 2 have an effect as both 1! and 0! equal 1
% we can subtract 1 if we do not forget to add 1 later
A = A(A>1) - 1 ;
B = B(B>1) - 1 ;

% what is the largest element?
maxE = max([A(:) ; B(:) ; 0]) ;

if maxE > 0 && ~isequal(A,B),
    % Create a 2 by m array holding the elements of A and B:
    % - the first column refers to the numerators of the ratio (A)
    % - the second column refers to the denominators of the ratio (B)
    R = sparse(A,1,1,maxE,2) + sparse(B,2,1,maxE,2) ;

    % By flipping, cumsumming (over columns) and flipping again ...
    R = flipud(cumsum(flipud(R))) ;
    % ... we get numbers that indicate how many times a factor is present in
    % the expanded ratio. R(X,1) indicates the number of times (X+1) is present in
    % the numerator in the expanded ratio, and R(X,2) the number of times
    % (X+1) is present in the denominator. The difference will give us the
    % exponent for taking the power of the value (X+1).
    R = (R(:,1) - R(:,2)).' ;    % exponents
    X = 2:(maxE+1) ;             % base numbers (we now do need to add the 1)

    % We only used the nonzero entries, since X^0 = 1, and this may affect
    % the time used for ...
    q = find(R) ;
    % ... getting the total product
    FR = prod(X(q).^R(q)) ; %
else
    % All elements were less than 2, or A equals B
    FR = 1 ;
end

author: Sable


text:
function Y = eminus(A,B)
Y = bsxfun(@minus, A, B);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
 function g = disimKernDiagGradient(kern, x, covDiag)

% DISIMKERNDIAGGRADIENT Compute the gradient of the DISIM kernel's diagonal wrt parameters.
% FORMAT
% DESC computes the gradient of functions of the diagonal of the
% single input motif kernel matrix with respect to the parameters of the kernel. The
% parameters' gradients are returned in the order given by the
% disimKernExtractParam command.
% ARG kern : the kernel structure for which the gradients are
% computed.
% ARG x : the input data for which the gradient is being computed.
% ARG factors : partial derivatives of the function of interest with
% respect to the diagonal elements of the kernel.
% RETURN g : gradients of the relevant function with respect to each
% of the parameters. Ordering should match the ordering given in
% disimKernExtractParam.
%
% SEEALSO : disimKernParamInit, kernDiagGradient, disimKernExtractParam, disimKernGradient
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% COPYRIGHT : Antti Honkela, 2007

% KERN

error('disimKernDiagGradient not yet implemented.')

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function S_lags=bss_make_lags(S,L)

% create a matrix containing lagged versions of some signal(s).
%
% Usage: S_lags=bss_make_lags(S,L)
%
% Input:
%   - S: n x T matrix containing the input signal(s),
%   - L: number of lagged versions of the signal(s).
%
% Output:
%   - S_lags: n*L x T matrix containing lagged versions of S, S_lags(t)=
%   [s1(t) ; s1(t-1) ; ...; s1(t-L+1); ... ; sn(t) ; sn(t-1) ; ...; sn(t-L+1)]
%
% WARNINGS:
%   * S_lags is zero-padded where necessary,
%   * We use the conventions make_lags(S,0)=makes_lags(S,1)=S.
%
% Developers:  - Cedric Fevotte (fevotte@tsi.enst.fr) - Emmanuel Vincent
% (emmanuel.vincent@irisa.fr) - Remi Gribonval (remi.gribonval@irisa.fr)

[n,T]=size(S);

if L==0
    S_lags=S;
else

    N=n*L;
    S_lags=zeros(N,T);
    for i=1:N
        q=floor((i-1)/L);
        r=mod(i-1,L);
        B=zeros(1,r+1); B(end)=1;
        S_lags(i,:)=filter(B,1,S(q+1,:));
    end

end

author: posenhuang

repo: deeplearningsourceseparation


text:
addpath(genpath(pwd))
x = input('Would you like to save the path for the next session of MATLAB? (y/(n))','s');
if (isempty(x))
    x = 'n';
end

if x=='y' || x=='Y'
    savepath
    disp('Path saved.')
else
    disp('Path added, but not saved.  It will only be used for this session of MATLAB.')
end


text:
function iou=boxiou(x1,y1,w1,h1,x2,y2,w2,h2)
% compute intersection over union of two bboxes
%

    bisect=boxIntersect(x1,x1+w1,y1+h1,y1,x2,x2+w2,y2+h2,y2);
    iou=0;
    if ~bisect, return; end

    bunion=boxUnion(x1,x1+w1,y1+h1,y1,x2,x2+w2,y2+h2,y2,bisect);

    assert(bunion>0,'something wrong with union computation');
    iou=bisect/bunion;

end

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
function cvt_1d_sampling ( n, it_num, s_num )

%*****************************************************************************80
%
%% CVT_1D_SAMPLING carries out the Lloyd algorithm.
%
%  Discussion:
%
%    This program is a variation of the CVT_1D_LLOYD method.
%    Instead of using an exact technique to determine the Voronoi
%    regions, it uses sampling.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of generators.
%
%    Input, integer IT_NUM, the number of CVT iterations.
%
%    Input, integer S_NUM, the number of samples in [0,1] to use
%    when estimating the Voronoi regions.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVT_1D_SAMPLING\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Use sampling to implement an approximation of Lloyd''s algorithm\n' );
  fprintf ( 1, '  in the 1D unit interval [0,1].\n' );

  if ( nargin < 1 )
    n = input ( '  Enter number of generators: ' );
  elseif ( ischar ( n ) )
    n = str2num ( n );
  end

  if ( nargin < 2 )
    it_num = input ( '  Enter number of iterations: ' );
  elseif ( ischar ( it_num ) )
    it_num = str2num ( it_num );
  end

  if ( nargin < 3 )
    s_num = input ( '  Enter number of sample points: ' );
  elseif ( ischar ( s_num ) )
    s_num = str2num ( s_num );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of generators is %d\n', n );
  fprintf ( 1, '  Number of iterations is %d\n', it_num );
  fprintf ( 1, '  Number of samples is %d\n', s_num );
%
%  For convenience, add a 0.0 and 1.0 point.
%  Also, sort the array.
%
  g = zeros ( n + 2, 1 );
  g_new = zeros ( n + 2, 1 );

  g(1,1) = 0.0;
  if ( 0 )
    g(2:n+1,1) = rand ( n, 1 );
  else
    g(2:n+1,1) = linspace ( 0.01, 0.02, n );
  end
  g(n+2,1) = 1.0;

  g = sort ( g );
%
%  Print the initial generators.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Initial generators:\n' );
  fprintf ( 1, '\n' );
  for k = 2 : n + 1
    fprintf ( 1, '  %4d  %f\n', k - 1, g(k,1) );
  end
%
%  Initialize the plotting arrays.
%
  g_plot = zeros ( n+2, it_num + 1 );

  step = 1 : it_num;
  e = nan ( it_num, 1 );
  gm = nan ( it_num, 1 );

  for it = 1 : it_num

    g_plot ( 1:n+2, it ) = g(1:n+2,1);

    s = rand ( s_num, 1 );
%
%  We need to "assign" each entry of S to the nearest X, then
%  replace X by the average of the values assigned to it.
%
%  There's a faster way for this 1D case.
%  1) Sort the S values.
%  2) Determine the values XM that delimit the Voronoi intervals.
%  3) For each XM, find the index LEFT in S such that S(LEFT) <= XM < S(LEFT+1).
%  4) Replace X(I) by the average of the S values that you have determined
%     are within its Voronoi interval.
%
    s = sort ( s );

    xm = zeros ( n - 1, 1 );
    for j = 1 : n - 1
      xm(j) = 0.5 * ( g(j+1) + g(j+2) );
    end

    left = r8vec_bracket4 ( s_num, s, n - 1, xm );
%
%  Compute the new generators.
%
    e(it,1) = 0.0;

    j = 1;
    g_new(j) = 0.0;

    k2 = 0;
    for j = 1 : n - 1
      k1 = k2;
      k2 = left(j);
      g_new(j+1) = sum ( s(k1+1:k2) ) / ( k2 - k1 );
      for k = k1 + 1 : k2
        e(it,1) = e(it,1) + ( s(k) - g(j+1) ).^2;
      end
    end

    k1 = k2;
    k2 = s_num;
    j = n + 1;
    g_new(j) = sum ( s(k1+1:k2) ) / ( k2 - k1 );
    for k = k1 + 1 : k2
      e(it,1) = e(it,1) + ( s(k) - g(j) ).^2;
    end

    j = n + 2;
    g_new(j) = 1.0;

    e(it,1) = e(it,1) / s_num;
%
%  Display the energy.
%
    figure ( 1 )
    plot ( step, log ( e ), 'm-*' )
    title ( 'Log (Energy)' )
    xlabel ( 'Step' )
    ylabel ( 'Energy' )
    grid
%
%  Compute the generator motion.
%
    gm(it,1) = sum ( ( g_new(:) - g(:) ).^2 ) / n;
%
%  Display the generator motion.
%
    figure ( 2 )
    plot ( step, log ( gm ), 'm-*' )
    title ( 'Log (Average generator motion)' )
    xlabel ( 'Step' )
    ylabel ( 'Motion' )
    grid
%
%  Update the generators.
%
    g = g_new;

  end

  g_plot(1:n+2,it_num+1) = g(1:n+2,1);

%
%  Print the current generators.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Current generators:\n' );
  fprintf ( 1, '\n' );
  for k = 2 : n + 1
    fprintf ( 1, '  %4d  %f\n', k - 1, g(k,1) );
  end
%
%  Plot the evolution of the locations of the generators.
%
  figure ( 3 )
  y = ( 0 : it_num );
  for k = 1 : n + 2
    plot ( g_plot(k,1:it_num+1), y )
    hold on;
  end
  grid on
  hold off;

  title ( 'Generator evolution' );
  xlabel ( 'Generator positions' );
  ylabel ( 'Iterations' );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVT_1D_SAMPLING\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function left = r8vec_bracket4 ( nt, t, ns, s )

%*****************************************************************************80
%
%% R8VEC_BRACKET4 finds the interval to each of a vector of values.
%
%  Discussion:
%
%    An R8VEC is a vector of R8's.
%
%    The routine always returns the index LEFT of the sorted array
%    T with the property that either
%    *  T is contained in the interval [ T(LEFT), T(LEFT+1) ], or
%    *  T < T(LEFT) = T(1), or
%    *  T > T(LEFT+1) = T(N).
%
%    The routine is useful for interpolation problems, where
%    the abscissa must be located within an interval of data
%    abscissas for interpolation, or the "nearest" interval
%    to the (extreme) abscissa must be found so that extrapolation
%    can be carried out.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NT, length of the input array.
%
%    Input, real T(NT), an array that has been sorted
%    into ascending order.
%
%    Input, integer NS, the number of points to be bracketed.
%
%    Input, real S(NS), values to be bracketed by entries of T.
%
%    Output, integer LEFT(NS).
%    LEFT(I) is set so that the interval [ T(LEFT(I)), T(LEFT(I)+1) ]
%    is the closest to S(I); it either contains S(I), or else S(I)
%    lies outside the interval [ T(1), T(NT) ].
%

%
%  Check the input data.
%
  if ( nt < 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8VEC_BRACKET4 - Fatal error!\n' );
    fprintf ( 1, '  NT must be at least 2.\n' );
    error ( 'R8VEC_BRACKET4 - Fatal error!' );
  end

  for i = 1 : ns

    left(i) = floor ( ( nt + 1 ) / 2 );
%
%  CASE 1: S < T(LEFT):
%  Search for S in [T(I), T(I+1)] for intervals I = 1 to LEFT-1.
%
    if ( s(i) < t(left(i)) )

      if ( left(i) == 1 )
        continue
      elseif ( left(i) == 2 )
        left(i) = 1;
        continue
      elseif ( t(left(i)-1) <= s(i) )
        left(i) = left(i) - 1;
        continue
      elseif ( s(i) <= t(2) )
        left(i) = 1;
        continue
      end
%
%  ...Binary search for S in [T(I), T(I+1)] for intervals I = 2 to LEFT-2.
%
      low = 2;
      high = left(i) - 2;

      while ( 1 )

        if ( low == high )
          left(i) = low;
          break
        end

        mid = floor ( ( low + high + 1 ) / 2 );

        if ( t(mid) <= s(i) )
          low = mid;
        else
          high = mid - 1;
        end

      end
%
%  CASE2: T(LEFT+1) < S:
%  Search for S in [T(I),T(I+1)] for intervals I = LEFT+1 to N-1.
%
    elseif ( t(left(i)+1) < s(i) )

      if ( left(i) == nt - 1 )
        continue
      elseif ( left(i) == nt - 2 )
        left(i) = left(i) + 1;
        continue
      elseif ( s(i) <= t(left(i)+2) )
        left(i) = left(i) + 1;
        continue
      elseif ( t(nt-1) <= s(i) )
        left(i) = nt - 1;
        continue
      end
%
%  ...Binary search for S in [T(I), T(I+1)] for intervals I = LEFT+2 to NT-2.
%
      low = left(i) + 2;
      high = nt - 2;

      while ( 1 )

        if ( low == high )
          left(i) = low;
          break
        end

        mid = floor ( ( low + high + 1 ) / 2 );

        if ( t(mid) <= s(i) )
          low = mid;
        else
          high = mid - 1;
        end

      end
%
%  CASE3: T(LEFT) <= S <= T(LEFT+1):
%  S is in [T(LEFT), T(LEFT+1)].
%
    else

    end

  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

author: johannesgerer


text:
% Add path
addpath(genpath('../toolbox_graph/toolbox_graph/'));
addpath(genpath('../icp/'));
addpath(genpath(''));

% Load data
base_path = '';
scan_path = '';

base_name = '';
scan_name = '';
load([base_path base_name]);
base_ver = Ver_3d';
base_tri = Tri';

load([scan_path scan_name]);
scan_ver = Ver_3d';
scan_tri = Tri';

% Select face parts of scan
% base_min = min(base_ver(3, :))-0.015;
% [index, value] = find(scan_ver(3, :)>=base_min*10);
% scan_ver = scan_ver(:, value);
% test_tri = [];
% for i = 1:length(scan_tri)
% 	[judge, ind] = ismember(scan_tri(:, i), value);
% 	if judge(1) && judge(2) && judge(3) == 1
% 		test_tri = [test_tri, ind];
%     end
% end
% scan_tri = test_tri;
% a = scan_ver';
% c = double(scan_tri');
% b = base_ver';
% d = double(base_tri');
% save('output/select_part/points/scan_ver.txt', 'a', '-ascii');
% save('output/select_part/points/scan_tri.txt', 'c', '-ascii');
% save('output/select_part/points/base_ver.txt', 'b', '-ascii');
% save('output/select_part/points/base_tri.txt', 'd', '-ascii');

% Calculation normal
base_norm = NormDirection(base_ver, base_tri);
scan_norm = NormDirection(scan_ver, scan_tri);

% Save mesh struct
field1 = 'faces'; field2 = 'vertices'; field3 = 'normals';
base_s = struct(field1, base_tri', field2, base_ver', field3, base_norm');
scan_s = struct(field1, scan_tri', field2, scan_ver', field3, scan_norm');

% Specify that surface normals are available and can be used.
Options.useNormals = 1;

% Specify that the source deformations should be plotted.
Options.plot = 1;

% Specify that the landmarks term are available and can be used.
Options.landmarks = 0;

% Selected landmarks
ls_source = [];
ls_target = [];
% x = base_ver(1, :); y = base_ver(2, :); z = base_ver(3, :);
% X = scan_ver(1, :); Y = scan_ver(2, :); Z = scan_ver(3, :);
% figure;hold on;
% scatter3(x(:), y(:), z(:) ,'b.');
% scatter3(x(ls_source), y(ls_source), z(ls_source) ,'r*');
% scatter3(X(:), Y(:), Z(:) ,'y.');
% scatter3(X(ls_target), Y(ls_target), Z(ls_target) ,'ro');

% Perform non-rigid ICP
[pointsTransformed, X] = nricp_landmarks(apple_s, scan_s, Options, ls_source, ls_target);
% faces = double(base_s.faces);
% save('output/point_cloud/vertices.txt', 'pointsTransformed', '-ascii');
% save('output/point_cloud/faces.txt', 'faces', '-ascii');

author: RhythmJnh

repo: Non-rigid-ICP


text:
function p = predict(Theta1, Theta2, X)
%PREDICT Predict the label of an input given a trained neural network
%   p = PREDICT(Theta1, Theta2, X) outputs the predicted label of X given the
%   trained weights of a neural network (Theta1, Theta2)

% Useful values
m = size(X, 1);
num_labels = size(Theta2, 1);

% You need to return the following variables correctly
p = zeros(size(X, 1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned neural network. You should set p to a
%               vector containing labels between 1 to num_labels.
%
% Hint: The max function might come in useful. In particular, the max
%       function can also return the index of the max element, for more
%       information see 'help max'. If your examples are in rows, then, you
%       can use max(A, [], 2) to obtain the max for each row.
%
X = [ones(m, 1) X];
XX = sigmoid(X*Theta1');
pp = sigmoid([ones(size(XX, 1), 1) XX] * Theta2');
[a, p] = max(pp, [], 2);

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

end

author: loserChen

repo: Coursera-MachineLearning


text:
classdef MimMain < CoreBaseClass
    % MimMain. Imports and provides access to data from the TD MIM Toolkit
    %
    %     MimMain provides access to data from the TD MIM Toolkit, and allows
    %     you to import new data. Data is accessed through one or more MimDataset
    %     objects. Your code should create a single MimMain object, and then ask
    %     it to create a MimDataset object for each dataset you wish to access.
    %
    %     MimMain is essentially a class factory for MimDatasets, but shares the
    %     MimReportingInterface (error/progress reporting) objects between all
    %     datasets, so you have a single error/progress reporting pipeline for
    %     your use of the TD MIM Toolkit.
    %
    %     To import a new dataset, construct a PTKImageInfo object with the file
    %     path and file name set to the image file. For DICOM files it is only
    %     necessary to specify the path since all image files in that directory
    %     will be imported. Then call CreateDatasetFromInfo. MimMain will import
    %     the data (if it has not already been imported) and return a new
    %     MimDataset object for that dataset.
    %
    %     To access an existing dataset you can use CreateDatasetFromInfo as
    %     above, or you can use CreateDatasetFromUid to retrieve a dataset which
    %     has peviously been imported, using the UID that was associated with
    %     that dataset.
    %
    %     Example
    %     -------
    %     Replace <image path> and <filenames> with the path and filenames
    %     to your image data.
    %
    %         image_info = PTKImageInfo( <image path>, <filenames>, [], [], [], []);
    %         mim = MimMain();
    %         dataset = mim.CreateDatasetFromInfo(image_info);
    %
    %     You can then obtain results from this dataset, e.g.
    %
    %         results = dataset.GetResult('PluginName');
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (SetAccess = private)
        FrameworkAppDef
        FrameworkSingleton
        Reporting          % Object for error and progress reporting
        ReportingWithCache % For the dataset, uses the same object, but warnings and messages are cached so multiple warnings can be displayed together
    end

    methods

        function obj = MimMain(framework_app_def, reporting)
            % Creates a new main MIM object using the configuration
            % specified by the supplied framework_app_def object
            obj.FrameworkAppDef = framework_app_def;
            obj.Reporting = reporting;
            obj.ReportingWithCache = MimReportingWithCache(obj.Reporting);
            obj.FrameworkSingleton = MimFrameworkSingleton.GetFrameworkSingleton(framework_app_def, obj.Reporting);

            bin_directory = framework_app_def.GetBinDirectory();
            CoreDiskUtilities.CreateDirectoryAndAddPathIfNotExisting(bin_directory);
            files_to_compile = framework_app_def.GetFilesToCompile(reporting);
            obj.FrameworkSingleton.CompileMexFileIfRequired(files_to_compile, bin_directory, obj.Reporting);
        end

        function Recompile(obj)
            % Forces recompilation of mex files

            bin_directory = obj.FrameworkAppDef.GetBinDirectory();
            files_to_compile = obj.FrameworkAppDef.GetFilesToCompile(obj.Reporting);
            obj.FrameworkSingleton.Recompile(files_to_compile, bin_directory, obj.Reporting);
        end

        function RebuildDatabase(obj)
            % Forces rebuilding of the image database

            obj.FrameworkSingleton.RebuildDatabase(obj.Reporting);
        end

        function DeleteCacheForAllDatasets(obj)
            % Delete automatic plugin results cache files on disk and in memory. Does not remove manually genereted files (edits, markers etc.)

            % Remove memory and disk results cache files for known datasets
            obj.FrameworkSingleton.GetDatasetApiCache().RemoveAllCachedFiles(obj.Reporting);

            % Do a hard clearout of the results cache to get rid of
            % anything else that's left
            obj.FrameworkAppDef.GetFrameworkDirectories.DeleteCacheForAllDatasets();
        end

        function ReloadPlugins(obj)
            % Forces a reload of plugin classes to pick up changes in properties etc.

            obj.FrameworkSingleton.GetPluginInfoMemoryCache.Clear();
        end

        function dataset_uids = FindDataset(obj, dataset_uid_prefix)
            % Searches for the dataset with this uid. If a dataset exactly
            % matches, this will be returned.
            % Otherwise, the specified string prefix will be matched
            % against all dataset uids in the cache, and if there is a
            % single match then this will be returned.
            %
            % Note: it is possible for datasets to exist in
            % the cache but not in the database; this can be fixed by
            % rebuilding the database

            if iscell(dataset_uid_prefix)
                dataset_uid_prefix = dataset_uid_prefix{1};
            end
            dataset_exists = (7 == exist(fullfile(obj.FrameworkAppDef.GetFrameworkDirectories.GetFrameworkDatasetCacheDirectory, dataset_uid_prefix), 'dir')) || (7 == exist(fullfile(obj.FrameworkAppDef.GetFrameworkDirectories.GetCacheDirectory, dataset_uid_prefix), 'dir'));

            if dataset_exists
                dataset_uids = {dataset_uid_prefix};
            else
                all_uids = obj.FrameworkAppDef.GetFrameworkDirectories().GetUidsOfAllDatasetsInCache();
                matches = strncmp(all_uids, dataset_uid_prefix, length(dataset_uid_prefix));
                dataset_uids = all_uids(matches);
            end
        end

        function dataset = Load(obj, filename)
            % Helper method to import the data from the specified file path and create a MimDataset for the first series found

            [uids, ~] = obj.ImportData(filename);
            if isempty(uids)
                dataset = [];
            else
                dataset = obj.CreateDatasetFromUid(uids{1});
            end
        end

        function dataset = CreateDatasetFromUid(obj, dataset_uid_prefix)
            % Creates a MimDataset object for a dataset specified by the uid. The
            % dataset must already be imported.

            dataset_uid = obj.FindDataset(dataset_uid_prefix);
            if isempty(dataset_uid)
                obj.Reporting.Error(MimErrors.UidNotFoundErrorId, 'Cannot find the dataset for this UID. Try importing the image using CreateDatasetFromInfo.');
            elseif numel(dataset_uid) > 1
                obj.Reporting.Error(MimErrors.UidNotFoundErrorId, 'More than one dataset for this UID prefix. Try specifying the complete uid.');
            else
                dataset_uid = dataset_uid{1};
            end

            dataset_disk_cache = obj.FrameworkSingleton.GetDatasetApiCache.GetDatasetDiskCache(dataset_uid, obj.Reporting);

            image_info = dataset_disk_cache.LoadData(obj.FrameworkAppDef.GetFrameworkConfig.ImageInfoCacheName, obj.Reporting);
            if isempty(image_info)
                obj.Reporting.Error(MimErrors.UidNotFoundErrorId, 'Cannot find the dataset for this UID. Try importing the image using CreateDatasetFromInfo.');
            end

            dataset = MimDataset(image_info, dataset_disk_cache, obj.FrameworkSingleton.GetLinkedDatasetChooserMemoryCache, obj.FrameworkSingleton.GetClassFactory, obj.ReportingWithCache);

            obj.AutoLinkDatasets(dataset);

            obj.RunCustomPostLoadFunction(dataset_uid, dataset);
        end

        function dataset = CreateDatasetFromInfo(obj, new_image_info)
            % Creates a MimDataset object for a dataset specified by the path,
            % filenames and/or uid specified in a PTKImageInfo object. The dataset is
            % imported from the specified path if it does not already exist.
            [image_info, dataset_disk_cache] = obj.ImportDataFromInfo(new_image_info, obj.Reporting);

            obj.FrameworkSingleton.AddToDatabase(image_info.ImageUid, obj.Reporting)

            dataset = MimDataset(image_info, dataset_disk_cache, obj.FrameworkSingleton.GetLinkedDatasetChooserMemoryCache, obj.FrameworkSingleton.GetClassFactory, obj.ReportingWithCache);

            obj.AutoLinkDatasets(dataset);

            obj.RunCustomPostLoadFunction(dataset.GetImageInfo.ImageUid, dataset);
        end

        function [uids, patient_ids] = ImportDataRecursive(obj, filename)
            % Identical to ImportData()
            [uids, patient_ids] = obj.ImportData(filename);
        end

        function [uids, patient_ids] = ImportData(obj, filename)
            % Imports data into the TD MIM Toolkit so that it can be accessed
            % from the CreateDatasetFromUid() method. The input argument is a string
            % containing the path to the data file to import. If the path points to
            % a single (non-DICOM) file, then only the file will be imported. If the
            % path points to a directory, or to a DICOM file, then all image files in
            % the directory and its subdirectories will be imported.

            % Only support a string input
            if ~ischar(filename)
                obj.Reporting.Error('MimMain:FileNotAsExpected', 'The file or directory passed to MimMain.ImportData() is not of the expected type.');
            end

            % Adds the files to the image database, which groups them into
            % series
            [uids, patient_ids] = obj.FrameworkSingleton.ImportData(filename, obj.ReportingWithCache);

            % Add the necessary files to the cache
            obj.ImportSeries(uids);

            % Save changes to the database
            obj.FrameworkSingleton.SaveImageDatabase(obj.Reporting);

            % Tell the image database to fire a database changed event. We do this here
            % rather than during the import to prevent multiple events being fired
            obj.FrameworkSingleton.ReportChangesToDatabase();
        end

        function RunCustomPostLoadFunction(obj, dataset_uid, dataset)
            % This method is called to run a user-defined function for linking
            % datasets

            obj.FrameworkAppDef.NewDatasetLoaded(dataset_uid, dataset, obj.Reporting);
        end

        function image_database = GetImageDatabase(obj)
            % Returns the Framework's image database
            image_database = obj.FrameworkSingleton.GetImageDatabase;
        end

        function directories = GetDirectories(obj)
            % Returns a MimDirectories object which can be used to query and create framework directories
            directories = obj.FrameworkAppDef.GetFrameworkDirectories;
        end

        function DeleteDatasets(obj, series_uids)
            % Deletes the datasets specified by a uid or a cell array of uids
            if isempty(series_uids)
                return;
            end

            if ~iscell(series_uids)
                series_uids = {series_uids};
            end

            for series_uid_cell = series_uids
                series_uid = series_uid_cell{1};

                dataset_to_delete = [];
                try
                    dataset_to_delete = obj.CreateDatasetFromUid(series_uid);
                catch exc
                    %
                    if MimErrors.IsErrorUidNotFound(exc.identifier)
                        obj.Reporting.ShowMessage('MimMain:UidNotFound', ['Failed to delete dataset because its UID could not be found']);
                    else
                        rethrow exc
                    end
                end

                if ~isempty(dataset_to_delete)
                    dataset_to_delete.DeleteCacheForThisDataset();
                    delete(dataset_to_delete);
                end
            end

            obj.FrameworkSingleton.GetDatasetApiCache.DeleteSeries(series_uids, obj.Reporting);
            obj.FrameworkSingleton.GetLinkedDatasetChooserMemoryCache.DeleteSeries(series_uids, obj.Reporting);
            obj.GetImageDatabase.DeleteSeries(series_uids, obj.Reporting);
        end

        function output = RunScript(obj, script_name, varargin)
            if nargin < 3
                parameters = [];
            end
            try
                script_class = feval(script_name);
                output = script_class.RunScript(obj, obj.Reporting, varargin{:});
            catch ex
                obj.Reporting.ErrorFromException('MimMain:ScriptFailure', ['The script ' script_name ' failed with the following error: ' ex.message], ex);
                output = [];
            end
        end

        function AutoLinkDatasets(obj, primary_dataset)
            linked_recorder = obj.FrameworkSingleton.GetLinkedDatasetRecorder();
            link_map = linked_recorder.LinkMap;
            primary_uid = primary_dataset.GetPrimaryDatasetUid();
            if link_map.isKey(primary_uid)
               link_record = link_map(primary_uid);
               link_map = link_record.LinkMap;
               for uid = link_map.keys
                   name = link_map(uid{1});
                   dataset_to_link = obj.CreateDatasetFromUid(uid{1});
                   primary_dataset.LinkDataset(name, dataset_to_link)
               end
            end
        end
    end

    methods (Access = private)

        function ImportSeries(obj, uids)
            if ~isempty(uids)
                for series_uid = uids
                    series_info = obj.FrameworkSingleton.GetSeriesInfo(series_uid{1});
                    image_info = series_info.GetImageInfo;
                    obj.ImportDataFromInfo(image_info, obj.Reporting);
                end
            end
        end

        function [image_info, dataset_disk_cache] = ImportDataFromInfo(obj, new_image_info, reporting)
            % Imports data into the TD MIM Toolkit so that it can be accessed
            % from the CreateDatasetFromUid() method. The input argument is a
            % PTKImageInfo object containing the path, filenames and file type of
            % the data to import. If you do not know the file type, use the
            % ImportData() method instead.

            if isempty(new_image_info.ImageUid)
                [series_uid, study_uid, modality] = MimMain.GetImageUID(new_image_info, reporting);
                new_image_info.ImageUid = series_uid;
                new_image_info.StudyUid = study_uid;
                new_image_info.Modality = modality;
            end

            dataset_disk_cache = obj.FrameworkSingleton.GetDatasetApiCache.GetDatasetDiskCache(new_image_info.ImageUid, reporting);
            image_info = dataset_disk_cache.LoadData(obj.FrameworkAppDef.GetFrameworkConfig.ImageInfoCacheName, reporting);
            if isempty(image_info)
                image_info = PTKImageInfo;
            end

            [image_info, anything_changed] = image_info.CopyNonEmptyFields(image_info, new_image_info);
            if (anything_changed)
                dataset_disk_cache.SaveData(obj.FrameworkAppDef.GetFrameworkConfig.ImageInfoCacheName, image_info, reporting);
            end
        end
    end

    methods (Static, Access = private)

        function [image_uid, study_uid, modality] = GetImageUID(image_info, reporting)
            % We need a unique identifier for each dataset. For DICOM files we use
            % the series instance UID. For other files we hash the full file path.

            if ~isempty(image_info.ImageFilenames) && isa(image_info.ImageFilenames{1}, 'CoreFilename')
                first_path = image_info.ImageFilenames{1}.Path;
                first_filename = image_info.ImageFilenames{1}.Name;
            else
                filenames = CoreDiskUtilities.GetDirectoryFileList(image_info.ImagePath, '*');
                first_filename = filenames{1};
                first_path = image_info.ImagePath;
            end

            single_image_info = MimGetSingleImageInfo(first_path, first_filename, image_info.ImageFileFormat, DMDicomDictionary.GroupingDictionary, reporting);
            image_uid = single_image_info.ImageUid;
            study_uid = single_image_info.StudyUid;
            modality = single_image_info.Modality;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function fiber_r=dtiFiberResample(fiber, value, flag)
%Resamples a fiber forcing a given number of nodes or a given step length.
%
% function fiber_r=dtiFiberResample(fiber, value, [flag = 'N'])
%
% Input: flag=='N'
%        Will resample the fiber to the number of points specified by
%        variable 'value' flag=='L' Will resample the fiber to the
%        intervals of length specified by variable 'value'
% Resampling is performed using interpolation with splines.
%
% Example:
%     fiber=fg.fibers{1}; fiber_r=dtiFiberResample(fiber, 10, 'N');
%     plot3(fiber_r(1,:),fiber_r(2,:),fiber_r(3,:),'bx-', ...
%     fiber(1,:),fiber(2,:),fiber(3,:),'ro--');
%     legend('Resampled', 'Original');
%
% See Also: dtiResampleFiberGroup
%
% (c) Vistalab

% HISTORY:
%  2007 ER wrote it
%  2010 ER fixed a bug: resampling of a fiber with
%  seriously unevenly spaced out breakpoints will now return a fiber with
%  evenly spaced ones.

if ~exist('flag', 'var') || isempty(flag)
    flag='N';
end

fiberNoNan=fiber(~isnan(fiber(:, :)));
sizeI=3; sizeJ=size(fiberNoNan, 1)/3;
fiberNoNan=reshape(fiberNoNan, sizeI, sizeJ);

%This block deals with issues a) all NaNs b) oNly one point is not NaN c)
%npoints=1 TODO: revise this block
if(strcmp(flag, 'N') && value==1)||(sizeJ==0) || (sizeJ==1)
    %fiber_r=mean(fiberNoNan, 2);
    error('All NANs, only one point is not NAN, or npoints==1');
end

%Regularize the nodes: Renumerate the nodes by distance from 1st
node2nodedist = squareform(pdist(fiber'));
archcumdist = cumsum(diag(node2nodedist, 1));

F = spline([0 archcumdist'], fiberNoNan);

switch flag
    case 'N'
        % Trajectory
        stepP = archcumdist(end)/(value-1);

    case 'L'
        stepP = value;
    otherwise
        error('Flag should be either L or N');
end

t = 0:stepP:archcumdist(end);
fiber_r=ppval(F,t);
return

author: vistalab