text:
function generalscript(path, lambdas, sigmas, tasks)

addpath(path);

kernelFID=fopen('all.kernel');
A=fread(kernelFID, 'float');
K=reshape(A,[sqrt(length(A)),sqrt(length(A))]);

perclassattributes = dlmread('all-perclass.attributes');
S = dlmread('all-perimage.attributes');
load experimentIndices

A=unique(perclassattributes, 'rows');
if length(dir('attrClasses*'))>0
    load attrClasses
    A=attrClasses;
end
% keyboard
attributesInput=A(trainClassesIndices,:);
attributesOutput=A(testClassesIndices,:);
nAttrs=size(A,2);
% allClass = dlmread('all.classid')+1;

% keyboard

nTrainingClasses=length(trainClassesIndices);
nTrainingInstances=length(trainInstancesIndices);

hist=[];
S=S(trainInstancesIndices,:);
Stest=attributesOutput;
% Stest(Stest==0)=-1;
% for i=1:size(Stest,1)
%     Stest(i,:)=Stest(i,:)-mean(Stest(i,:));
%     Stest(i,:)=Stest(i,:)/norm(Stest(i,:));
% end

for t=1:length(tasks)
    nNewTasks=tasks(t);

    disp('Getting labels...');

%     Y=zeros(nTrainingInstances,nTrainingClasses+nNewTasks);%-1;%-ones(nTrainingInstances,nTrainingClasses);%/(nTrainingClasses-1);
%     for i=1:nTrainingInstances
%         Y(i,trainInstancesLabels(i))=1;
%     end

    disp('Precalculating statistics...');
    KTrain=K(trainInstancesIndices, trainInstancesIndices);
    KK=KTrain'*KTrain;

%     S=[attributesInput; rand(nNewTasks, nAttrs)];
%     S=[attributesInput; sign(rand(nNewTasks, nAttrs)-0.5)];
%     S(S==0)=-1;
%     for i=1:size(S,1)
%         S(i,:)=S(i,:)-mean(S(i,:));
%         S(i,:)=S(i,:)/norm(S(i,:));
%     end

    KYS=KTrain*S;
    SS=S'*S;

    % keyboard
    for s=1:length(sigmas)
        sigma=sigmas(s);
        KYS_invSS=KYS/(SS+sigma*eye(size(S,2)));

        disp('Learning...');
        % step=10^-8;
        % for t=1:500
        % %     keyboard
        %     W = W - step * (XX*W*SS -XYS + lambda*W);
        % end
        for lambdaIndex=1:length(lambdas)
            lambda=lambdas(lambdaIndex);

            Alpha=(KK+lambda*KTrain)\KYS_invSS;

            disp('Predicting...');
            KTest=K(testInstancesIndices,trainInstancesIndices);

            pred=Stest*Alpha'*KTest';
            [scores, classPred]=max(pred',[],2);

            disp('Evaluating...');
            gt=testInstancesLabels;
            r=mean(gt==classPred')
%             confusionmat(gt, classPred')
            hist=[hist; [lambda, sigma, nNewTasks, r]];
        %     break;
            save('hist4', 'hist');
        end
    end
end
end

repo: Embarrassingly-simple-ZSL


text:
%% [core,U] = tensor_hosvd(T,t,r)
%
%%% Normal HoSVD
% [core,U] = tensor_hosvd(T)
% [core,U] = tensor_hosvd(T,0)
% [core,U] = tensor_hosvd(T,[0 0 0])
%
%%% Truncate mode-3 basis matrice (keep the first 2 eigenvalues)
% [core,U] = tensor_hosvd(T,[0 0 2])
%
%%% Perform mode-3 rank-2 partial svd
% [core,U] = tensor_hosvd(T,0,[0 0 2])
%
function [core,U] = tensor_hosvd(T,t,r)
  core = T;

  if(nargin < 3) r = zeros(1,ndims(T)); end
  if(nargin < 2) t = zeros(1,ndims(T)); end

  if(t == 0) t = zeros(1,ndims(T)); end
  if(r == 0) r = zeros(1,ndims(T)); end

  for i = 1:ndims(T)
    M{i} = double(tenmat(T,i));
    %M{i} = unfolding(double(T),i);

    %%% perform partial svd
    if(r(i) > 0)
      [U{i},S{i},V{i}] = svds(M{i},r(i),'L');
    else
      [U{i},S{i},V{i}] = svd(M{i});
    end

    %%% truncate basis matrice
    if(t(i) > 0)
      U{i}(:,t(i)+1:end) = 0;
    end

    %[m n] = size(U{i});
    %s(i) = n;
    %core = tensor(folding(U{i}'*unfolding(double(core),i), i, s));
    %core = double(ttm(tensor(core),U{i},i,'t'));
    core = ttm(core,U{i}',i);
  end
end

author: andrewssobral

repo: lrslibrary


text:
function yi=akimai(x,y,xi)
%
% Usage: yi=akimai(x,y,xi)
%
%        Given vectors x and y (of the same length)
%        and the array xi at which to interpolate,
%        fits piecewise cubic polynomials and returns
%        the interpolated values yi at xi.
%
% Ref. : Hiroshi Akima, Journal of the ACM, Vol. 17, No. 4, October 1970,
%        pages 589-602.
%
% Programmer: N. Shamsundar, University of Houston, 6/2002
%   Correction to lines 32-33, 9/2004,  motivated by Gilford Ward,
%     to make routine work correctly for linear data.
%
% Notes: Use only for precise data, as the fitted curve passes through the
%        given points exactly. This routine is useful for plotting a pleasingly
%        smooth curve through a few given points for purposes of plotting.
%
x=x(:); y=y(:); xi=xi(:); n=length(x);
if n~=length(y), error('input x and y arrays must be of same length'), end
dx=diff(x); if any(dx <= 0) error('input x-array must be in strictly ascending order'), end
if any(xi<x(1)) | any(xi > x(n))
  warning('All interpolation points xi must lie between x(1) and x(n)')
end
m=diff(y)./dx;
mm=2*m(1)-m(2);     mmm=2*mm-m(1);     % augment at left
mp=2*m(n-1)-m(n-2); mpp=2*mp-m(n-1);   % augment at right
m1=[mmm; mm; m; mp; mpp];              % slopes
dm=abs(diff(m1)); f1=dm(3:n+2); f2=dm(1:n); f12=f1+f2;
id=find(f12 > 1e-8*max(f12)); b=m1(2:n+1);
b(id)=(f1(id).*m1(id+1)+f2(id).*m1(id+2))./f12(id);
c=(3*m-2*b(1:n-1)-b(2:n))./dx;
d=(b(1:n-1)+b(2:n)-2*m)./dx.^2;

% Loop replaced by vector ops following tip from johannes.korsawe@volkswagen.de
% 1/19/2006
%
[ncnt,bin]=histc(xi,x);
bin=min(bin,n-1);
bb=bin(1:length(xi));
wj=xi-x(bb);
yi=((wj.*d(bb) +c(bb)).*wj+b(bb)).*wj+y(bb);

author: Sable

repo: mcbench-benchmarks


text:
function [] = ManyPieces()
%MANYPIECES Produce a photomosiac from an image and a directory of images
%   Given an image and a directory of images create a photomosaic
%	Paramets may be modified within m file
%	Written by Justin Dailey
%	dailejl@auburn.edu

%Mosaic Parameters
%Number of thumbs spanning width of img
thumb_ratio = 100;
%Default thumb directory (only if saving thumbs)
thumb_dir = strcat(pwd, filesep, 'Thumbs');
%Coefficient for scaling of output tiles
output_scale = 2;
%Output file name
output_file_name = 'Mosaic.jpg';
%Dithering distance
dval = 1;

%List of file types
file_types = {['*.BMP;*.GIF;*.HDF;*.JPEG;*.JPG;*.PBM;*.PCX;*.PGM;',...
    '*.PNG;*.PNM;*.PPM;*.RAS;*.TIFF;*.TIF;*.XWD'],'MATLAB Graphical Files'};
comp_file_types = {'BMP' 'GIF' 'HDF' 'JPEG' 'JPG' 'PBM' 'PCX' 'PGM' ...
    'PNG' 'PNM' 'PPM' 'RAS' 'TIFF' 'TIF' 'XWD'};

%% Get main image and directory of images
img_path = uigetfile(file_types, 'Select the main mosiac image');
img_dir = uigetdir();

%Get original img size
ref_img = imread(img_path);
img_size = size(ref_img);
img_size = img_size(1:2);

%% Calc values
%Set size for thumbnail based on original size
thumb_pixels = floor(img_size(1)/thumb_ratio);
thumb_size = [thumb_pixels thumb_pixels];
%Make sure h/w are proportional to thumb size
new_height = floor(img_size(2)/thumb_pixels);
num_tiles = [thumb_ratio new_height];
new_size = [thumb_ratio new_height].*thumb_pixels;
ref_img = imresize(ref_img, new_size);

%% Get all directory images
dir_files = dir(img_dir);
mosaic_ind = 1;
%For each image add it to array
for dir_ind = 1:length(dir_files)
    if ~dir_files(dir_ind).isdir
        file_name = dir_files(dir_ind).name;
        %Check file extension
        [~, ~, ext, ~] = fileparts(file_name);
        if max(strcmpi(ext(2:end), comp_file_types))
            mosaic_files{mosaic_ind} = file_name;
            mosaic_ind = mosaic_ind+1;
        end
    end
end

%% Resize directory images into thumbs
progress = waitbar(0, 'Creating Mosaic Thumbnails...');
num_files = length(mosaic_files);
mosaic_imgs = cell(1, num_files);
%mkdir(thumb_dir);
%Resize each image
for mosaic_ind = 1:num_files
    img = imread([img_dir, filesep, mosaic_files{mosaic_ind}]);
    %if read in grayscale img
    if size(img, 3) == 1
        img = ind2rgb(img, gray(256));
    end
    res_img = uint8(imresize(img, thumb_size*output_scale));
    thumbs{mosaic_ind} = res_img;
    %imwrite(res_img, [thumb_dir, filesep, strcat(num2str(mosaic_ind),'.jpg')], 'jpg');
    waitbar(mosaic_ind/num_files, progress);
end
close(progress);

%% Calculating thumbnail averages
progress = waitbar(0, 'Calculating Thumbnail Averages...');
for mosaic_ind = 1:num_files
    %calc average vals for thumbs
    cur_thumb = thumbs{mosaic_ind};
    RGB_vals{mosaic_ind} = mean(reshape(cur_thumb, [], 3), 1);
    waitbar(mosaic_ind/num_files, progress);
end
close(progress);

%% For each tile of image find closest matching thumb
progress = waitbar(0, 'Creating Photomosiac...');
pic_map = zeros(num_tiles);
tiles_done = 0;
for row_tile = 1:num_tiles(1)
    for col_tile = 1:num_tiles(2)
        closest = 1;
        shortest_dist = 1000;
        %get mean vals for the image tiles
        cur_tile = ref_img(thumb_pixels*(row_tile-1)+1:thumb_pixels*(row_tile), ...
        thumb_pixels*(col_tile-1)+1:thumb_pixels*(col_tile),:);
        cur_RGB = mean(reshape(cur_tile, [], 3), 1);
        %find the closest thumb to each tile
        for thumb_tile = 1:num_files
            dist = calc_distance(RGB_vals{thumb_tile}, cur_RGB);
            %if new pt is closer
            if dist < shortest_dist
                %implement dithering to limit use of tile in an area
                if isempty(find( ...
                        pic_map(max(row_tile-dval,1): ...
                        min(row_tile+dval,num_tiles(1)),...
                        max(col_tile-dval,1): ...
                        min(col_tile+dval,num_tiles(2))) ...
                        == thumb_tile, 1))
                    shortest_dist = dist;
                    pic_map(row_tile, col_tile) = thumb_tile;
                end
            end
        end
        tiles_done = tiles_done + 1;
        waitbar(tiles_done/(num_tiles(1)*num_tiles(2)), progress);
    end
end
close(progress);

%% Take mapping of thumbnails and create photomosaic
progress = waitbar(0, 'Assembling Photomosaic...');
for row_tile = 1:num_tiles(1)
    cur_row = thumbs{pic_map(row_tile, 1)};
    for col_tile = 2:num_tiles(2)
        cur_row = horzcat(cur_row, thumbs{pic_map(row_tile, col_tile)});
    end
    if row_tile == 1
        mosaic = cur_row;
    else
        mosaic = vertcat(mosaic, cur_row);
        clear cur_row;
    end
    waitbar(row_tile/num_tiles(1), progress);
end
close(progress);

imshow(mosaic)
imwrite(mosaic, output_file_name, 'jpg');
end

function distance = calc_distance(pt1, pt2)
    distance = sqrt(sum((pt1-pt2).^2));
end


text:
function rtk=udsol_sppins(rtk)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The output is centered on the INS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global glc

x=rtk.ins.x; P=rtk.ins.P;

[rr,Cen]=blh2xyz(x(7:9)); T=Dblh2Dxyz(x(7:9));
rtk.sol.pos=rr';
rtk.sol.vel=(Cen*x(4:6))';
rtk.sol.att=x(1:3)'/glc.D2R;

if rtk.sol.att(3)>=0
    rtk.sol.att(3)=360-rtk.sol.att(3);
else
    rtk.sol.att(3)=-rtk.sol.att(3);
end

posvar=P(7:9,7:9); posvar=T*posvar*T';        %blh_var to xyz_var
velvar=P(4:6,4:6); velvar=Cen*velvar*Cen';    %enu_var to xyz_var
attvar=P(1:3,1:3); attvar=attvar./glc.D2R^2;  %rad_var to deg_var

rtk.sol.posP(1)=posvar(1,1);
rtk.sol.posP(2)=posvar(2,2);
rtk.sol.posP(3)=posvar(3,3);
rtk.sol.posP(4)=posvar(1,2);
rtk.sol.posP(5)=posvar(2,3);
rtk.sol.posP(6)=posvar(1,3);

rtk.sol.velP(1)=velvar(1,1);
rtk.sol.velP(2)=velvar(2,2);
rtk.sol.velP(3)=velvar(3,3);
rtk.sol.velP(4)=velvar(1,2);
rtk.sol.velP(5)=velvar(2,3);
rtk.sol.velP(6)=velvar(1,3);

rtk.sol.attP(1)=attvar(1,1);
rtk.sol.attP(2)=attvar(2,2);
rtk.sol.attP(3)=attvar(3,3);
rtk.sol.attP(4)=attvar(1,2);
rtk.sol.attP(5)=attvar(2,3);
rtk.sol.attP(6)=attvar(1,3);

% update clk and isb
rtk.sol.dtr(1)=rtk.x(rtk.ic+1)/glc.CLIGHT;
rtk.sol.dtr(2)=rtk.x(rtk.ic+2)/glc.CLIGHT;
rtk.sol.dtr(3)=rtk.x(rtk.ic+3)/glc.CLIGHT;
rtk.sol.dtr(4)=rtk.x(rtk.ic+4)/glc.CLIGHT;
rtk.sol.dtr(5)=rtk.x(rtk.ic+5)/glc.CLIGHT;

return

repo: GINav


text:
function title = p00_title ( problem )

%*****************************************************************************80
%
%% P00_TITLE returns the title for any problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM, the index of the problem.
%
%    Output, string TITLE, the title of the problem.
%
  if ( problem == 1 )
    title = p01_title ( );
  elseif ( problem == 2 )
    title = p02_title ( );
  elseif ( problem == 3 )
    title = p03_title ( );
  elseif ( problem == 4 )
    title = p04_title ( );
  elseif ( problem == 5 )
    title = p05_title ( );
  elseif ( problem == 6 )
    title = p06_title ( );
  elseif ( problem == 7 )
    title = p07_title ( );
  elseif ( problem == 8 )
    title = p08_title ( );
  elseif ( problem == 9 )
    title = p09_title ( );
  elseif ( problem == 10 )
    title = p10_title ( );
  elseif ( problem == 11 )
    title = p11_title ( );
  elseif ( problem == 12 )
    title = p12_title ( );
  elseif ( problem == 13 )
    title = p13_title ( );
  elseif ( problem == 14 )
    title = p14_title ( );
  elseif ( problem == 15 )
    title = p15_title ( );
  elseif ( problem == 16 )
    title = p16_title ( );
  elseif ( problem == 17 )
    title = p17_title ( );
  elseif ( problem == 18 )
    title = p18_title ( );
  elseif ( problem == 19 )
    title = p19_title ( );
  elseif ( problem == 20 )
    title = p20_title ( );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_TITLE - Fatal error!\n' );
    fprintf ( 1, '  Illegal problem number = %d\n', problem );
    error ( 'P00_TITLE - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function exportModel(modelName)
% exportModel(modelName)
%
% Zips a trained model so that it can be uploaded to our homepage.
%
% Copyright by Holger Caesar, 2016

% Settings
global glBaseFolder;
if strcmp(modelName, 'frcn')
    exportName = 'FRCN_VOC2010_model';
    exportFiles = {'net-epoch-16.mat', 'net-train.pdf', 'log.txt', 'resultsEpochFinalTest.mat'};
    exportFilesPrefix = fullfile(glBaseFolder, 'Features', 'CNN-Models', 'FRCN', 'VOC2010', 'VOC2010-testRelease');
elseif strcmp(modelName, 'fcn')
    exportName = 'FCN_SiftFlow_model';
    exportFiles = {'labelings-test-epoch50', 'net-epoch-50.mat', 'net-opts.mat', 'net-train.pdf', 'log.txt', 'stats-test-epoch50.mat'};
    exportFilesPrefix = fullfile(glBaseFolder, 'Features', 'CNN-Models', 'FCN', 'SiftFlow', 'fcn16s-testRelease');
elseif strcmp(modelName, 'e2s2_fast')
    exportName = 'E2S2_SiftFlow_model_fast';
    exportFiles = {'labelings-test-epoch10', 'net-epoch-10.mat', 'net-opts.mat', 'net-train.pdf', 'log.txt', 'stats-test-epoch10.mat'};
    exportFilesPrefix = fullfile(glBaseFolder, 'Features', 'CNN-Models', 'E2S2', 'SiftFlow', 'Run1', 'SiftFlow_e2s2_run1_exp1');
elseif strcmp(modelName, 'e2s2_full')
    exportName = 'E2S2_SiftFlow_model_full';
    exportFiles = {'labelings-test-epoch25', 'net-epoch-25.mat', 'net-opts.mat', 'net-train.pdf', 'log.txt', 'stats-test-epoch25.mat'};
    exportFilesPrefix = fullfile(glBaseFolder, 'Features', 'CNN-Models', 'E2S2', 'SiftFlow', 'Run1', 'SiftFlow_e2s2_run1_exp2');
else
    error('Error: Unknown model: %s', modelName);
end

% Define paths
exportFolder = fullfile(glBaseFolder, 'Export');
if ~exist(exportFolder, 'dir')
    mkdir(exportFolder);
end

% Create zip file
exportFilePath = fullfile(exportFolder, [exportName, '.zip']);
if exist(exportFilePath, 'file')
    fprintf('Using existing zip file: %s\n', exportFilePath);
else
    fprintf('Creating zip file: %s\n', exportFilePath);
    zip(exportFilePath, strcat(exportFilesPrefix, filesep, exportFiles));
end

author: nightrome

repo: matconvnet-calvin


text:
function problem = inputValidation(problem)
%
% This function runs through the problem struct and sets any missing fields
% to the default value. If a mandatory field is missing, then it throws an
% error.
%
% INPUTS:
%   problem = a partially completed problem struct
%
% OUTPUTS:
%   problem = a complete problem struct, with validated fields
%

%%%% Check the function handles:

if ~isfield(problem,'func')
    error('Field ''func'' cannot be ommitted from ''problem''');
else
    if ~isfield(problem.func,'dynamics')
        error('Field ''dynamics'' cannot be ommitted from ''problem.func'''); end
    if ~isfield(problem.func,'pathObj'), problem.func.pathObj = []; end
    if ~isfield(problem.func,'bndObj'), problem.func.bndObj = []; end
    if ~isfield(problem.func,'pathCst'), problem.func.pathCst = []; end
    if ~isfield(problem.func,'bndCst'), problem.func.bndCst = []; end
end

%%%% Check the initial guess (also compute nState and nControl):
if ~isfield(problem, 'guess')
    error('Field ''guess'' cannot be ommitted from ''problem''');
else
    if ~isfield(problem.guess,'time')
        error('Field ''time'' cannot be ommitted from ''problem.guess'''); end
    if ~isfield(problem.guess, 'state')
        error('Field ''state'' cannot be ommitted from ''problem.guess'''); end
    if ~isfield(problem.guess, 'control')
        error('Field ''control'' cannot be ommitted from ''problem.guess'''); end

    % Compute the size of the time, state, and control based on guess
    [checkOne, nTime] = size(problem.guess.time);
    [nState, checkTimeState] = size(problem.guess.state);
    [nControl, checkTimeControl] = size(problem.guess.control);

    if nTime < 2 || checkOne ~= 1
        error('guess.time must have dimensions of [1, nTime], where nTime > 1');
    end

    if checkTimeState ~= nTime
        error('guess.state must have dimensions of [nState, nTime]');
    end
    if checkTimeControl ~= nTime
        error('guess.control must have dimensions of [nControl, nTime]');
    end

end

%%%% Check the problem bounds:
if ~isfield(problem,'bounds')
    problem.bounds.initialTime = [];
    problem.bounds.finalTime = [];
    problem.bounds.state = [];
    problem.bounds.initialState = [];
    problem.bounds.finalState = [];
    problem.bounds.control = [];
else

    if ~isfield(problem.bounds,'initialTime')
        problem.bounds.initialTime = []; end
    problem.bounds.initialTime = ...
        checkLowUpp(problem.bounds.initialTime,1,1,'initialTime');

    if ~isfield(problem.bounds,'finalTime')
        problem.bounds.finalTime = []; end
    problem.bounds.finalTime = ...
        checkLowUpp(problem.bounds.finalTime,1,1,'finalTime');

    if ~isfield(problem.bounds,'state')
        problem.bounds.state = []; end
    problem.bounds.state = ...
        checkLowUpp(problem.bounds.state,nState,1,'state');

    if ~isfield(problem.bounds,'initialState')
        problem.bounds.initialState = []; end
    problem.bounds.initialState = ...
        checkLowUpp(problem.bounds.initialState,nState,1,'initialState');

    if ~isfield(problem.bounds,'finalState')
        problem.bounds.finalState = []; end
    problem.bounds.finalState = ...
        checkLowUpp(problem.bounds.finalState,nState,1,'finalState');

    if ~isfield(problem.bounds,'control')
        problem.bounds.control = []; end
    problem.bounds.control = ...
        checkLowUpp(problem.bounds.control,nControl,1,'control');

end

end

function input = checkLowUpp(input,nRow,nCol,name)
%
% This function checks that input has the following is true:
%   size(input.low) == [nRow, nCol]
%   size(input.upp) == [nRow, nCol]

if ~isfield(input,'low')
    input.low = -inf(nRow,nCol);
end

if ~isfield(input,'upp')
    input.upp = inf(nRow,nCol);
end

[lowRow, lowCol] = size(input.low);
if lowRow ~= nRow || lowCol ~= nCol
    error(['problem.bounds.' name ...
        '.low must have size = [' num2str(nRow) ', ' num2str(nCol) ']']);
end

[uppRow, uppCol] = size(input.upp);
if uppRow ~= nRow || uppCol ~= nCol
    error(['problem.bounds.' name ...
        '.upp must have size = [' num2str(nRow) ', ' num2str(nCol) ']']);
end

if sum(sum(input.upp-input.low < 0))
    error(...
        ['problem.bounds.' name '.upp must be >= problem.bounds.' name '.low!']);
end

end

author: MatthewPeterKelly

repo: OptimTraj


text:
function structM = findStructureMatrixForOneZSlice(allSegmentsM, zSliceUniformValue, CTOriginalZValues, CTSliceThickness, CTDeltaX, CTDeltaY, reusableZerosM)
%"findStructureMatrixForOneZSlice"
% returns structM, a matrix with 1's and 0's corresponding to
% the location of the points of a structure on a particular z slice.
%
% allSegmentsM = planC{indexS.structures}(structNum).rasterSegments;
% zSliceUniformValue = desired value of the z slice of the output matrix, structM
% CTOriginalZValues = [planC{indexS.scan}.scanInfo(:).zValue];
% reusableZerosM = an all-zero matrix the size of structM output.

structM = reusableZerosM;

%Test if structure is more than 1 full slice from the requested zSlice, if
%so return the all zero matrix, as the structure is not present.
maxStructureZVal = max(allSegmentsM(:,1));
minStructureZVal = min(allSegmentsM(:,1));
CTSpacing = CTOriginalZValues(2) - CTOriginalZValues(1);
if zSliceUniformValue < (minStructureZVal-CTSpacing) | zSliceUniformValue > (maxStructureZVal+CTSpacing)
    return
end
%END

%try to find a slice with the desired z value
if maxStructureZVal == minStructureZVal
	indV = 1:size(allSegmentsM,1);
else
	indV = find(abs(allSegmentsM(:,1) - zSliceUniformValue) < 0.001);  %mask values on this slice
end
if ~isempty(indV) %no need to interpolate
  segmentsM = allSegmentsM(indV(:),7:9);     %segments
  for i = 1 : size(segmentsM,1)
    structM(segmentsM(i,1),segmentsM(i,2):segmentsM(i,3)) = 1;
  end
else %simple interpolation necessary
  %find the two slices closest to this uniform z value
  cranialSlice = max(find(CTOriginalZValues < zSliceUniformValue));
  caudalSlice  = min(find(CTOriginalZValues >= zSliceUniformValue));
  if isempty(cranialSlice)
    cranialSlice = caudalSlice;
  end
  if isempty(caudalSlice)
    caudalSlice = cranialSlice;
  end
  cranialZValue = CTOriginalZValues(cranialSlice);
  caudalZValue  = CTOriginalZValues(caudalSlice);

  sliceSpacing = caudalZValue - cranialZValue;
  spaceToCranial = zSliceUniformValue - cranialZValue;
  spaceFromCaudal = caudalZValue - zSliceUniformValue;
%   cranialIsNearest = spaceFromCaudal >= spaceToCranial;

  %If we are within the valid slice thickness of either slice, no need to
  %interpolate.
%   closeEnough = (spaceToCranial < (CTSliceThickness(cranialSlice)/2) | (spaceFromCaudal < (CTSliceThickness(caudalSlice)/2)));

  %note:  we don't really need the farther matrices except for the else of closeEnough, so speed could be optimized here.
  indCranial = find(abs(allSegmentsM(:,1) - cranialZValue) < 0.001);  %mask values on this slice
  indCaudal = find(abs(allSegmentsM(:,1) - caudalZValue ) < 0.001);  %mask values on this slice

  segmentsCranialM = allSegmentsM(indCranial(:),7:9);     %segments
  segmentsCaudalM = allSegmentsM(indCaudal(:),7:9);     %segments
%   if closeEnough & cranialIsNearest %use cranial
%     for i = 1 : size(segmentsCranialM,1)
%       structM(segmentsCranialM(i,1),segmentsCranialM(i,2):segmentsCranialM(i,3)) = 1;
%     end
%   elseif closeEnough & ~cranialIsNearest %use caudal
%     for i = 1 : size(segmentsCaudalM,1)
%       structM(segmentsCaudalM(i,1),segmentsCaudalM(i,2):segmentsCaudalM(i,3)) = 1;
%     end
%   else %do some interpolation
    cranialM = reusableZerosM;
    caudalM  = reusableZerosM;
    for i = 1 : size(segmentsCranialM,1)
      cranialM(segmentsCranialM(i,1),segmentsCranialM(i,2):segmentsCranialM(i,3)) = 1;
    end
    for i = 1:size(segmentsCaudalM,1)
      caudalM(segmentsCaudalM(i,1),segmentsCaudalM(i,2):segmentsCaudalM(i,3)) = 1;
    end

    if spaceToCranial == 0
        structM = cranialM;
        return;
    end
    if spaceFromCaudal == 0
        structM = caudalM;
        return;
    end

    %Interpolate using distance weighting.  Points interpolated to .5 or
    %greater are on.

    distanceToCranialRegion = spaceToCranial;
    distanceToCaudalRegion  = spaceFromCaudal;

    %Distance from middle point to diagonal neighbor in plane.
    diagDist = sqrt(CTDeltaX^2 + CTDeltaY^2);
    inPlaneNeighborsDistance = [diagDist CTDeltaY diagDist;CTDeltaX 0 CTDeltaX;diagDist CTDeltaY diagDist];
    distToCranialNeighbors = sqrt(inPlaneNeighborsDistance.^2 + distanceToCranialRegion^2);
    distToCaudalNeighbors  = sqrt(inPlaneNeighborsDistance.^2 + distanceToCaudalRegion^2);

    weightsCranial = 1./distToCranialNeighbors;
    weightsCaudal  = 1./distToCaudalNeighbors;
    totalWeight = sum([sum(weightsCranial(:)), sum(weightsCaudal(:))]);

    structM = conv2(single(cranialM), weightsCranial/totalWeight, 'same') + conv2(single(caudalM), weightsCaudal/totalWeight, 'same');
    structM = structM > .5;

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [solution,output] = NMPC_Solve(x0,p,options)
    tStart = Timer();

    persistent lambdaInit muInit uInit xInit zInit LAMBDAInit;

    global ParNMPCGlobalVariable
    dim = ParNMPCGlobalVariable.dim;
    N   = ParNMPCGlobalVariable.N;

    % only for the very first problem
    solutionInitGuess = ParNMPCGlobalVariable.solutionInitGuess;

    % Reshape
    sizeSeg     = N/options.DoP;
    pSplit      = reshape(p,      dim.p,        sizeSeg, options.DoP);

    if isempty(xInit)
        lambdaInit = reshape(solutionInitGuess.lambda, dim.lambda,   sizeSeg, options.DoP);
        uInit      = reshape(solutionInitGuess.u,      dim.u,        sizeSeg, options.DoP);
        xInit      = reshape(solutionInitGuess.x,      dim.x,        sizeSeg, options.DoP);
        LAMBDAInit = reshape(solutionInitGuess.LAMBDA, dim.x, dim.x, sizeSeg, options.DoP);
        % assert
        if coder.target('MATLAB') % Normal excution
            checkFeasibility(solutionInitGuess,p);
            checkOptions(options);
        end
    end
    if dim.mu ~= 0
        muInit     = reshape(solutionInitGuess.mu,     dim.mu,       sizeSeg, options.DoP);
    end
    if dim.z ~= 0
        zInit      = reshape(solutionInitGuess.z,      dim.z,        sizeSeg, options.DoP);
    end

    % Init
    lambdaSplit = lambdaInit;
    if dim.mu ~= 0
        muSplit = muInit;
    else
        muSplit = zeros(dim.mu,sizeSeg, options.DoP);
    end
    uSplit      = uInit;
    xSplit      = xInit;
    if dim.z ~= 0
        zSplit = zInit;
    else
        zSplit = zeros(dim.z,sizeSeg, options.DoP);
    end
    LAMBDASplit = LAMBDAInit;

    output   = createOutput();
    solution = createSolution(dim,N);
    costL = 0;
    error = 0;
    tLineSearch = 0;
    if (options.rhoInit == options.rhoEnd)
        options.maxIterInit = options.maxIterTotal;
        options.tolInit = options.tolEnd;
    end
    rho     = options.rhoInit;
    mode    = 1;

    % initialize a filter [LAll;xEq+C;flag]
    lineSearchFilterWidth = 10;
    lineSearchFilter = zeros(3,lineSearchFilterWidth);

    % Iteration
    for iter=1:options.maxIterTotal
        % backup
        lambdaSplit_k = lambdaSplit;
        muSplit_k     = muSplit;
        uSplit_k      = uSplit;
        xSplit_k      = xSplit;
        %% Search direction
        [lambdaSplit,muSplit,uSplit,xSplit,zSplit,LAMBDASplit,KKTError,costL,tSearchDirection] = ...
            NMPC_Solve_SearchDirection(x0,pSplit,rho,lambdaSplit,muSplit,uSplit,xSplit,zSplit,LAMBDASplit);
        output.timeElapsed.searchDirection = output.timeElapsed.searchDirection + tSearchDirection;
        %% Line search
        stepSize = 1;
        if options.isLineSearch
            switch options.lineSearchMethod
                case 'merit'
                    scaling = 1.5;
                    lambdaSplitAbs = abs(lambdaSplit);
                    lambdaSplitMax3 =  max(lambdaSplitAbs,[],2);
                    phiX = max(lambdaSplitMax3,[],3);
                    muSplitAbs = abs(muSplit);
                    muSplitMax3 =  max(muSplitAbs,[],2);
                    phiC = max(muSplitMax3,[],3);
                    [lambdaSplit,muSplit,uSplit,xSplit,stepSize,tLineSearch] = ...
                        NMPC_LineSearch_Merit(x0,pSplit,rho,lambdaSplit_k,muSplit_k,uSplit_k,xSplit_k,...
                                        lambdaSplit,muSplit,uSplit,xSplit,phiX*scaling,phiC*scaling);
                case 'filter'
                    [lambdaSplit,muSplit,uSplit,xSplit,lineSearchFilter,stepSize,tLineSearch] = ...
                        NMPC_LineSearch_Filter(x0,pSplit,rho,lambdaSplit_k,muSplit_k,uSplit_k,xSplit_k,...
                                        lambdaSplit,muSplit,uSplit,xSplit,lineSearchFilter);
                otherwise
                    if coder.target('MATLAB') % Normal excution
                            warning('Specified linesearch method is not supported!');
                    end
            end
            output.timeElapsed.lineSearch = output.timeElapsed.lineSearch + tLineSearch;
        end
        %% KKT error
        tKKTErrorCheck = 0;
        if options.checkKKTErrorAfterIteration
            [KKTError,costL,tKKTErrorCheck]   = ...
                NMPC_KKTError(x0,pSplit,rho,lambdaSplit,muSplit,uSplit,xSplit);
            KKTErrorScaling = 1;
        else
            % avoid chattering of the num of iter
            KKTErrorScaling = 5;
        end
        sMax = 100;
        lambda_L1Norm = norm(lambdaSplit(:),1);
        HuScaling =  (lambda_L1Norm + norm(muSplit(:),1))/(dim.lambda+dim.mu)/N;
        HuScaling = max(sMax,HuScaling)/sMax;
        costateEquationScaling = max(sMax,lambda_L1Norm/dim.lambda/N)/sMax;
        error = max([KKTError.stateEquation*KKTErrorScaling;...
                     KKTError.C/KKTErrorScaling;...
                     KKTError.Hu/HuScaling/KKTErrorScaling;...
                     KKTError.costateEquation/costateEquationScaling/KKTErrorScaling]);
        output.timeElapsed.KKTErrorCheck = output.timeElapsed.KKTErrorCheck + tKKTErrorCheck;
        %% print
        if coder.target('MATLAB') % Normal excution
            if options.printLevel == 1
                printMsg = ['Iter: ',num2str(iter),...
                            '   KKTError: ',num2str(error),...
                            '   Cost: ',num2str(costL),...
                            '   rho: ',num2str(rho),...
                            '   Step size: ',num2str(stepSize)];
                disp(printMsg);
            end
        end
        %% barrier parameter update
        switch mode
            case 1 % init
                if (error < options.tolInit) || (iter >= options.maxIterInit)
                    mode = 2;
                    lambdaInit = lambdaSplit;
                    if dim.mu ~= 0
                        muInit     = muSplit;
                    end
                    uInit      = uSplit;
                    xInit      = xSplit;
                    if dim.z ~= 0
                        zInit     = zSplit;
                    end
                    LAMBDAInit = LAMBDASplit;

                    output.iterInit          = iter;
                    if (options.rhoInit == options.rhoEnd) && (error < options.tolEnd)
                        output.exitflag = 1;
                        break;
                    end
                end
            case 2 % decay
                rho = rho * options.rhoDecayRate;
                rho(rho<options.rhoEnd) = options.rhoEnd;
                if (error<options.tolEnd) && (rho == options.rhoEnd)
                    output.exitflag = 1;
                    break;
                end
        end
    end
    %% Reshape
    solution.lambda = reshape(lambdaSplit, dim.lambda,   N);
    solution.mu     = reshape(muSplit,     dim.mu,       N);
    solution.u      = reshape(uSplit,      dim.u,        N);
    solution.x      = reshape(xSplit,      dim.x,        N);
    solution.z      = reshape(zSplit,      dim.z,        N);
    solution.LAMBDA = reshape(LAMBDASplit, dim.x, dim.x, N);

    output.cost            = costL;
    output.KKTError        = error;
    output.rho             = rho;
    output.iterTotal       = iter;

    tEnd = Timer();

    output.timeElapsed.total = tEnd-tStart;
end

author: deng-haoyang

repo: ParNMPC


text:
function cdf = frechet_cdf ( x, alpha )

%*****************************************************************************80
%
%% FRECHET_CDF evaluates the Frechet CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA, the parameter.
%    It is required that 0.0 < ALPHA.
%
%    Input, real X, the argument of the CDF.
%
%    Output, real CDF, the value of the CDF.
%
  if ( alpha <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FRECHET_CDF - Fatal error!\n' );
    fprintf ( 1, '  ALPHA <= 0.0.\n' );
    error ( 'FRECHET_CDF - Fatal error!' );
  end

  if ( x <= 0.0 )
    cdf = 0.0;
  else
    cdf = exp ( - 1.0 / x^alpha );
  end

  return
end


text:
function [geom_data, valid_ids] = Read_geom_files_dynamic(users, model_param_data_dir)

    geom_data = [];
    valid_ids = [];

    load('../../pca_generation/pdm_68_aligned_wild.mat');

    user_files = [];

    for i=1:numel(users)

        geom_file = [model_param_data_dir, '/' users{i} '.csv'];

        user_files = cat(1, user_files, str2num(users{i}(3:5)));

        if(i == 1)
            tab = readtable(geom_file);
            column_names = tab.Properties.VariableNames;
            valid_ind = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'confidence'));
            shape_inds = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'p_'));
        end

        res = dlmread(geom_file, ',', 1, 0);

        valid = res(:, valid_ind) > 0.7;
        res = res(:, shape_inds);

        % Do not consider global parameters
        res = res(:, 7:end);

        actual_locs = res * V';
        res = cat(2, actual_locs, res);

        valid_ids = cat(1, valid_ids, valid);

        geom_data_curr_p = res;

        geom_data = cat(1, geom_data, geom_data_curr_p);

    end

    if(numel(users) > 0)
        uq_ids = unique(user_files)';

        for u=uq_ids
            geom_data(user_files==u,:) = bsxfun(@plus, geom_data(user_files==u,:), -median( geom_data(valid_ids & user_files==u,:)));
        end
    end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function value = r8_nint ( x )

%*****************************************************************************80
%
%% R8_NINT returns the nearest integer to an R8.
%
%  Example:
%
%        X        R8_NINT
%
%      1.3         1
%      1.4         1
%      1.5         1 or 2
%      1.6         2
%      0.0         0
%     -0.7        -1
%     -1.1        -1
%     -1.6        -2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the value.
%
%    Output, integer VALUE, the nearest integer to X.
%
  if ( x < 0.0 )
    s = -1;
  else
    s = 1;
  end

  value = s * round ( abs ( x ) + 0.5 );

  return
end


text:
function [x, c, funVal, ValueL]=mtLogisticC(A, y, z, opts)
%
%%
% Function mtLogisticC:
%      Logistic Loss for Multi-task Learning
%             with the (group) L1/Lq-norm Constraint
%
%% Problem
%
%  min \sum_i - weight_i * log (p_i) + z * sum_j ||x^j||_q
%
%  p_i= 1/ (1+ exp(-y_i (x_l' * a_i + c_l) ) ) denotes the probability
%
%  a_i denotes a training sample,
%      and a_i' corresponds to the i-th row of the data matrix A
%
%  x^j denotes the j-th row of x
%  x_l denotes the l-th column of x, for the l-th task
%  c_l is the intercept for the l-th task
%
%  y_i (either 1 or -1) is the response
%
%  weight_i denotes the weight for the i-th sample
%  In this implementation, we assume weight_{il}=1/m
%
%  For the case that the multi tasks share the same data
%  matrix, please refer to the functions:
%            mcLeastR and mcLogisticR.
%
%% Input parameters:
%
%  A-         Matrix of size m x n
%                A can be a dense matrix
%                         a sparse matrix
%                         or a DCT matrix
%  y -        Response vector (of size m x 1)
%  z -        L1/Lq norm regularization parameter (z >=0)
%  opts-      Optional inputs (default value: opts=[])
%
%% Output parameters:
%  x-         The obtained weight of size n x k
%  c-         The obtained intercept if size 1 x k
%  funVal-    Function value during iterations
%
%% Copyright (C) 2009-2010 Jun Liu, and Jieping Ye
%
% You are suggested to first read the Manual.
%
% For any problem, please contact with Jun Liu via j.liu@asu.edu
%
% Last modified on February 19, 2010.
%
%% Related papers
%
% [1]  Jun Liu, Shuiwang Ji, and Jieping Ye, Multi-Task Feature Learning
%      Via Efficient L2,1-Norm Minimization, UAI, 2009
%
% [2]  Jun Liu, Lei Yuan, Songcan Chen and Jieping Ye, Multi-Task Feature Learning
%      Via Efficient L2,1-Norm Minimization, Technical Report ASU, 2009.
%
%% Related functions:
%
%  sll_opts, initFactor, pathSolutionLeast
%  mtLeastR, mtLeastC, mtLogisticR, ep21d
%
%%

%% Verify and initialize the parameters
%%
if (nargin <4)
    error('\n Inputs: A, y, z, and opts.ind should be specified!\n');
end

[m,n]=size(A);

if (length(y) ~=m)
    error('\n Check the length of y!\n');
end

if (z<=0)
    error('\n z should be positive!\n');
end

opts=sll_opts(opts); % run sll_opts to set default values (flags)

%% Detailed initialization
%%

% Initialize ind and q
if ~isfield(opts,'ind')
    error('\n In mtLeastR, .ind should be specified');
else
    ind=opts.ind;
    k=length(ind)-1;

    if ind(k+1)~=m
        error('\n Check opts.ind');
    end
end

% Initialize q
if (~isfield(opts,'q'))
    q=2; opts.q=2;
else  % currently, we only implement q=2
    q=opts.q;
    if (q~=2)
        error('\n Currently, we only implement the case q=2');
    end
end

%% Normalization

% Please refer to sll_opts for the definitions of mu, nu and nFlag
%
% If .nFlag =1, the input matrix A is normalized to
%                     A= ( A- repmat(mu, m,1) ) * diag(nu)^{-1}
%
% If .nFlag =2, the input matrix A is normalized to
%                     A= diag(nu)^{-1} * ( A- repmat(mu, m,1) )
%
% Such normalization is done implicitly
%     This implicit normalization is suggested for the sparse matrix
%                                    but not for the dense matrix
%

if (opts.nFlag~=0)
    if (isfield(opts,'mu'))
        mu=opts.mu;
        if(size(mu,2)~=n)
            error('\n Check the input .mu');
        end
    else
        mu=mean(A,1);
    end

    if (opts.nFlag==1)
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=n)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,1)/m).^(0.5); nu=nu';
        end
    else % .nFlag=2
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=m)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,2)/n).^(0.5);
        end
    end

    ind_zero=find(abs(nu)<= 1e-10);    nu(ind_zero)=1;
    % If some values in nu is typically small, it might be that,
    % the entries in a given row or column in A are all close to zero.
    % For numerical stability, we set the corresponding value to 1.
end

if (~issparse(A)) && (opts.nFlag~=0)
    fprintf('\n -----------------------------------------------------');
    fprintf('\n The data is not sparse or not stored in sparse format');
    fprintf('\n The code still works.');
    fprintf('\n But we suggest you to normalize the data directly,');
    fprintf('\n for achieving better efficiency.');
    fprintf('\n -----------------------------------------------------');
end

%% Starting point initialization

p_flag=(y==1);                  % the indices of the postive samples

for i=1:k
    ind_i=(ind(i)+1):ind(i+1);  % indices for the i-th group

    m1(1,i)=sum(p_flag(ind_i));      % the total number of the positive samples
    m2(1,i)=length(ind_i)-m1(1,i);   % the total number of the negative samples
end

% initialize a starting point
if opts.init==2
    x=zeros(n,k); c=zeros(1,k);
else
    if isfield(opts,'x0')
        x=opts.x0;
        if ( size(x,1)~=n || size(x,2)~=k )
            error('\n Check the input .x0');
        end
    else
        x=zeros(n,k);
    end

    if isfield(opts,'c0')
        c=opts.c0;

        if ( length(c)~=k )
            error('\n Check the input .c0');
        end
    else
        c=log(m1./m2);
    end
end

Ax=zeros(m,1); % m x 1
% compute Ax: Ax_i= A_i * x_i
for i=1:k
    ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
    m_i=ind(i+1)-ind(i);          % number of samples in the i-th group

    if (opts.nFlag==0)
        Ax(ind_i,1)=A(ind_i,:)* x(:,i);
    elseif (opts.nFlag==1)
        invNu=x(:,i)./nu; mu_invNu=mu * invNu;
        Ax(ind_i,1)=A(ind_i,:)*invNu -repmat(mu_invNu, m_i, 1);
    else
        Ax(ind_i,1)=A(ind_i,:)*x(:,i)-repmat(mu*x(:,i), m, 1);
        Ax(ind_i,1)=Ax./nu(ind_i,1);
    end
end

%% The main program
% The Armijo Goldstein line search schemes

if (opts.lFlag==0)

    lambda0=0;
    % a guess of the root

    bFlag=0; % this flag tests whether the gradient step only changes a little

    L=1/m; % the intial guess of the Lipschitz continuous gradient

    % assign xp with x, and Axp with Ax
    xp=x; Axp=Ax; xxp=zeros(n,k);
    cp=c;         ccp=zeros(1,k);

    alphap=0; alpha=1;

    for iterStep=1:opts.maxIter
        % --------------------------- step 1 ---------------------------
        % compute search point s based on xp and x (with beta)
        beta=(alphap-1)/alpha;    s=x + beta* xxp;   sc=c + beta* ccp;

        % --------------------------- step 2 ---------------------------
        % line search for L and compute the new approximate solution x

        % compute the gradient (g) at s
        As=Ax + beta* (Ax-Axp);

        % aa= - diag(y) * (A * s + sc)
        vec_sc=zeros(m,1);
        for i=1:k
            ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
            vec_sc(ind_i,1)=sc(i);
        end
        aa=- y.*(As+ vec_sc);

        % fun_s is the logistic loss at the search point
        bb=max(aa,0);
        fun_s= sum(sum ( log( exp(-bb) +  exp(aa-bb) ) + bb ) ) / m;

        % compute prob=[p_1;p_2;...;p_m]
        prob=1./( 1+ exp(aa) );

        % b= - diag(y) * (1 - prob)
        b= -y.*(1-prob) / m;

        gc=zeros(1,k);                     % the gradient of c
        for i=1:k
            ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
            gc(1,k)=sum(b(ind_i));
        end

        % compute g= AT b, the gradient of x
        for i=1:k
            ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group

            if (opts.nFlag==0)
                tt =A(ind_i,:)'*b(ind_i,1);
            elseif (opts.nFlag==1)
                tt= A(ind_i,:)'*b(ind_i,1) - sum(b(ind_i,1)) * mu';
                tt=tt./nu(ind_i,1);
            else
                invNu=b(ind_i,1)./nu(ind_i,1);
                tt=A(ind_i,:)'*invNu - sum(invNu)*mu';
            end

            g(:,i)= tt;
        end

        % copy x and Ax to xp and Axp
        xp=x;    Axp=Ax;
        cp=c;

        while (1)
            % let s walk in a step in the antigradient of s to get v
            % and then do the L1/Lq-norm regularized projection
            v=s-g/L; c= sc- gc/L;

            % L1/Lq-norm regularized projection
            [x, lambda, zf_step]=ep21d(v, n, k, z, lambda0);
            lambda0=lambda;

            v=x-s;  % the difference between the new approximate solution x
            % and the search point s

            % compute Ax: Ax_i= A_i * x_i
            for i=1:k
                ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
                m_i=ind(i+1)-ind(i);          % number of samples in the i-th group

                if (opts.nFlag==0)
                    Ax(ind_i,1)=A(ind_i,:)* x(:,i);
                elseif (opts.nFlag==1)
                    invNu=x(:,i)./nu; mu_invNu=mu * invNu;
                    Ax(ind_i,1)=A(ind_i,:)*invNu -repmat(mu_invNu, m_i, 1);
                else
                    Ax(ind_i,1)=A(ind_i,:)*x(:,i)-repmat(mu*x(:,i), m, 1);
                    Ax(ind_i,1)=Ax./nu(ind_i,1);
                end
            end

            % aa= - diag(y) * (A * x + c)
            vec_sc=zeros(m,1);
            for i=1:k
                ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
                vec_sc(ind_i,1)=c(i);
            end
            aa=- y.*(Ax+ vec_sc);

            % fun_s is the logistic loss at the search point
            bb=max(aa,0);
            fun_x= sum(sum ( log( exp(-bb) +  exp(aa-bb) ) + bb ) ) / m;

            r_sum=(norm(v,'fro')^2 + norm(c-sc,2)^2) / 2;
            l_sum=fun_x - fun_s - sum(sum(v.* g)) - (c-sc)* gc';

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is fun_x <= fun_s + <v, g> + <c ,gc>
            %                           + L/2 * (<v,v> + <c-sc,c-sc> )
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                % fprintf('\n L=%5.6f',L);
            end
        end

        % --------------------------- step 3 ---------------------------
        % update alpha and alphap, and check whether converge
        alphap=alpha; alpha= (1+ sqrt(4*alpha*alpha +1))/2;

        ValueL(iterStep)=L;

        xxp=x-xp;   ccp=c-cp;

        funVal(iterStep)=fun_x;

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=norm(xxp,'fro');
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=norm(xp,'fro');    norm_xxp=norm(xxp,'fro');
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end
    end
end

%% Adaptive line search scheme

if (opts.lFlag==1)

    lambda0=0;
    % a guess of the root

    bFlag=0; % this flag tests whether the gradient step only changes a little

    L=1/m; % the intial guess of the Lipschitz continuous gradient

    % assign xp with x, and Axp with Ax
    xp=x; Axp=Ax; xxp=zeros(n,k);
    cp=c;             ccp=zeros(1,k);
    % t is the upper bound of absolute value of x

    gamma=1;
    % we shall set the value of gamma = L,
    % and L is appropriate for the starting point

    for iterStep=1:opts.maxIter
        while (1)
            if (iterStep~=1)
                alpha= (-gamma+ sqrt(gamma*gamma + 4* L * gamma)) / (2*L);
                beta= (gamma - gamma* alphap) / (alphap * gamma + alphap* L * alpha);
                % beta is the coefficient for generating search point s

                s=x + beta* xxp;   sc=c + beta* ccp;
                As=Ax + beta* (Ax-Axp);
            else
                alpha= (-1+ sqrt(5)) / 2;  beta=0;
                s=x; sc= c;
                As=Ax;
            end

            % aa= - diag(y) * (A * s + sc)
            vec_sc=zeros(m,1);
            for i=1:k
                ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
                vec_sc(ind_i,1)=sc(i);
            end
            aa=- y.*(As+ vec_sc);

            % fun_s is the logistic loss at the search point
            bb=max(aa,0);
            fun_s= sum(sum ( log( exp(-bb) +  exp(aa-bb) ) + bb ) ) / m;

            % compute prob=[p_1;p_2;...;p_m]
            prob=1./( 1+ exp(aa) );

            % b= - diag(y) * (1 - prob)
            b= -y.*(1-prob) / m;

            gc=zeros(1,k);                     % the gradient of c
            for i=1:k
                ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
                gc(1,k)=sum(b(ind_i));
            end

            % compute g= AT b, the gradient of x
            for i=1:k
                ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group

                if (opts.nFlag==0)
                    tt =A(ind_i,:)'*b(ind_i,1);
                elseif (opts.nFlag==1)
                    tt= A(ind_i,:)'*b(ind_i,1) - sum(b(ind_i,1)) * mu';
                    tt=tt./nu(ind_i,1);
                else
                    invNu=b(ind_i,1)./nu(ind_i,1);
                    tt=A(ind_i,:)'*invNu - sum(invNu)*mu';
                end

                g(:,i)= tt;
            end

            % let s walk in a step in the antigradient of s to get v
            % and then do the L1/Lq-norm regularized projection
            v=s-g/L; cnew= sc- gc/L;

            % projection
            [xnew, lambda, zf_step]=ep21d(v, n, k, z, lambda0);
            lambda0=lambda;

            v=xnew-s;  % the difference between the new approximate solution x
            % and the search point s

            % compute Axnew = A *xnew: Ax_i= A_i * xnew_i
            for i=1:k
                ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
                m_i=ind(i+1)-ind(i);          % number of samples in the i-th group

                if (opts.nFlag==0)
                    Axnew(ind_i,1)=A(ind_i,:)* xnew(:,i);
                elseif (opts.nFlag==1)
                    invNu=xnew(:,i)./nu; mu_invNu=mu * invNu;
                    Axnew(ind_i,1)=A(ind_i,:)*invNu -repmat(mu_invNu, m_i, 1);
                else
                    Axnew(ind_i,1)=A(ind_i,:)*xnew(:,i)-repmat(mu*xnew(:,i), m, 1);
                    Axnew(ind_i,1)=Axnew./nu(ind_i,1);
                end
            end

            % aa= - diag(y) * (A * x + c)
            vec_sc=zeros(m,1);
            for i=1:k
                ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
                vec_sc(ind_i,1)=cnew(i);
            end
            aa=- y.*(Axnew+ vec_sc);

            % fun_x is the logistic loss at the new approximate solution
            % xnew
            bb=max(aa,0);
            fun_x= sum(sum ( log( exp(-bb) +  exp(aa-bb) ) + bb ) ) / m;

            r_sum=(norm(v,'fro')^2 + norm(cnew-sc,2)^2) / 2;
            l_sum=fun_x - fun_s - sum(sum(v.* g)) - (cnew-sc)* gc';

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is fun_x <= fun_s + <v, g> + <c ,gc>
            %                           + L/2 * (<v,v> + <c-sc,c-sc>  + <t-s_t, t-s_t>)
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                % fprintf('\n L=%5.6f',L);
            end
        end

        gamma=L* alpha* alpha;    alphap=alpha;
        % update gamma, and alphap

        ValueL(iterStep)=L;

        tao=L * r_sum / l_sum;
        if (tao >=5)
            L=L*0.8;
        end
        % decrease the value of L

        xp=x;     Axp=Ax;
        x=xnew;   Ax=Axnew;
        cp=c;      c=cnew;
        % update x, t, c, and Ax

        xxp=x-xp;  ccp=c-cp;

        funVal(iterStep)=fun_x;

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=norm(xxp,'fro');
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=norm(xp,'fro');    norm_xxp=norm(xxp,'fro');
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end
    end
end

author: thomaskuestner


text:
% Build the C/C++ files provided in the package
%
% @author: B. Schauerte
% @date:   2009
% @url:    http://cvhci.anthropomatik.kit.edu/~bschauer/

cpp_files=dir('*.cpp');
for i=1:length(cpp_files)
  fprintf('Building %d of %d: %s\n',i,length(cpp_files),cpp_files(i).name);
  mex(cpp_files(i).name);
end

author: opencobra

repo: cobratoolbox


text:
function logS = batch_uniformizer(sourceDirName, destinDirName, destinDCMDirName)
%
% function logS = batch_uniformizer(dirName)
%
% This function uniformizes CERR plans in the passed source directory and
% writes them to the destination directory
%
% APA, 08/19/2011

% sourceDirName = '/Users/aptea/Documents/MSKCC/Projects/TRISMUS/TRISMUS CERR';
% destinDCMDirName = '/Users/aptea/Documents/MSKCC/Projects/TRISMUS/TRISMUS CERR DICOM';
% destinDirName = '/Users/aptea/Documents/MSKCC/Projects/TRISMUS/TRISMUS CERR UNIFORM';

fileC = {};
if strcmpi(sourceDirName,'\') || strcmpi(sourceDirName,'/')
    filesTmp = getCERRfiles(sourceDirName(1:end-1));
else
    filesTmp = getCERRfiles(sourceDirName);
end
fileC = [fileC filesTmp];

% Initialize waitbar
hWait = waitbar(0,'Uniformizing cohort of CERR plans. Please wait...');

logS = struct('maxDiff','', 'fileName', '');

numFiles = length(fileC);

%Load CERR plan
for iFile=1:numFiles

    %Update Waitbar
    drawnow

    waitbar(iFile/numFiles,hWait)

    %Load CERR plan
    try
        planC = loadPlanC(fileC{iFile}, tempdir);
        planC = updatePlanFields(planC);
        indexS = planC{end};

        %Check for mesh representation and load meshes into memory
        currDir = cd;
        meshDir = fileparts(which('libMeshContour.dll'));
        cd(meshDir)
        for strNum = 1:length(planC{indexS.structures})
            if isfield(planC{indexS.structures}(strNum),'meshRep') && ~isempty(planC{indexS.structures}(strNum).meshRep) && planC{indexS.structures}(strNum).meshRep
                try
                    calllib('libMeshContour','loadSurface',planC{indexS.structures}(strNum).strUID,planC{indexS.structures}(strNum).meshS)
                catch
                    planC{indexS.structures}(strNum).meshRep    = 0;
                    planC{indexS.structures}(strNum).meshS      = [];
                end
            end
        end
        cd(currDir)

        stateS.optS = opts4Exe([getCERRPath,'CERROptions.json']);

        %Check color assignment for displaying structures
        [assocScanV,relStrNumV] = getStructureAssociatedScan(1:length(planC{indexS.structures}),planC);
        for scanNum = 1:length(planC{indexS.scan})
            scanIndV = find(assocScanV==scanNum);
            for i = 1:length(scanIndV)
                strNum = scanIndV(i);
                colorNum = relStrNumV(strNum);
                if isempty(planC{indexS.structures}(strNum).structureColor)
                    color = stateS.optS.colorOrder( mod(colorNum-1, size(stateS.optS.colorOrder,1))+1,:);
                    planC{indexS.structures}(strNum).structureColor = color;
                end
            end
        end

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

        %Check whether uniformized data is in cellArray format.
        if ~isempty(planC{indexS.structureArray}) && iscell(planC{indexS.structureArray}(1).indicesArray)
            planC = setUniformizedData(planC,planC{indexS.CERROptions});
            indexS = planC{end};
        end

        if length(planC{indexS.structureArrayMore}) ~= length(planC{indexS.structureArray})
            for saNum = 1:length(planC{indexS.structureArray})
                if saNum == 1
                    planC{indexS.structureArrayMore} = struct('indicesArray', {[]},...
                        'bitsArray', {[]},...
                        'assocScanUID',{planC{indexS.structureArray}(saNum).assocScanUID},...
                        'structureSetUID', {planC{indexS.structureArray}(saNum).structureSetUID});

                else
                    planC{indexS.structureArrayMore}(saNum) = struct('indicesArray', {[]},...
                        'bitsArray', {[]},...
                        'assocScanUID',{planC{indexS.structureArray}(saNum).assocScanUID},...
                        'structureSetUID', {planC{indexS.structureArray}(saNum).structureSetUID});
                end
            end
        end

        %Create uniform slices
        % planC = createUniformlySlicedPlanC(planC);

        % Check scan grid and correct for numerical noise
        zValuesV = [planC{indexS.scan}.scanInfo(:).zValue];
        zDiffV = diff(zValuesV);
        meanDiff = mean(zDiffV);
        maxDifference = max(abs((zDiffV - meanDiff)));

        % Write to log
        [jnk, fileName] = fileparts(fileC{iFile});
        logS(iFile).maxDiff = maxDifference;
        logS(iFile).fileName = fileName;

        if maxDifference < 0.05
            newZvalsV = linspace(zValuesV(1), zValuesV(end), length(zValuesV));
        else
            continue;
        end

        for sliceNum = 1:length(planC{indexS.scan}.scanInfo)

            planC{indexS.scan}.scanInfo(sliceNum).zValue = newZvalsV(sliceNum);

        end

        % Assign appropriate z-coordinates to structures
        for strNum = 1:length(planC{indexS.structures}(strNum))

            for slcNum = 1:length(planC{indexS.structures}(strNum).contour)

                for segNum = 1:length( planC{indexS.structures}(strNum).contour(slcNum).segments)

                    if ~isempty( planC{indexS.structures}(strNum).contour(slcNum).segments(segNum).points)

                         planC{indexS.structures}(strNum).contour(slcNum).segments(segNum).points(:,3) = newZvalsV(slcNum);

                    end

                end

            end

        end

        %Save plan to destination directory
        save_planC(planC,[], 'passed', fullfile(destinDirName,fileName));

        % Export DICOM
        dcmFolderName =  fullfile(destinDCMDirName,strtok(fileName,'.'));
        if ~exist(dcmFolderName, 'dir')
            mkdir(dcmFolderName)
        end

        CERRExportDICOM(fullfile(destinDirName,fileName), dcmFolderName)

    catch

        disp([fileC{iFile}, ' failed to load'])
        continue

    end

end

close(hWait)


text:
function [high_fissure_indices, ref_image] = PTKGetMaxFissurePoints(fissure_approximation, lung_mask, fissureness, image_roi, image_size)
    % PTKGetMaxFissurePoints. function for finding candidate points of high
    %     fissureness given an initial fissure approximation.
    %
    %     PTKGetMaxFissurePoints is an intermediate stage in segmenting the
    %     lobes. It is not intended to be a general-purpose algorithm.
    %
    %     For more information, see
    %     [Doel et al., Pulmonary lobe segmentation from CT images using
    %     fissureness, airways, vessels and multilevel B-splines, 2012]
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    data_points_indices = find(fissure_approximation);

    bin_size = 2;

    [high_fissure_indices, ref_image] = GetFissurePoints(data_points_indices, image_size, lung_mask, fissureness, image_roi, bin_size);
end

function [maximum_fissureness_indices, ref_image] = GetFissurePoints(indices_for_model, image_size, lung_segmentation, fissureness, image_roi, bin_size)

    % Get the coordinates of every voxel in the lung segmentation
    candidate_indices = find(lung_segmentation.RawImage);

    % An image for debugging
    ref_image = zeros(image_size, 'uint8');
    ref_image(candidate_indices) = 2;

    % Remove points with low fissureness
    max_fissureness = max(fissureness.RawImage(candidate_indices));
    fissureness_threshold = max_fissureness/3;
    intensity_threshold_hu = -900;
    intensity_threshold = image_roi.HounsfieldToGreyscale(intensity_threshold_hu);
    candidate_indices = candidate_indices((fissureness.RawImage(candidate_indices) > fissureness_threshold) & (image_roi.RawImage(candidate_indices) > intensity_threshold));
    ref_image(candidate_indices) = 4;

    if isempty(candidate_indices)
        maximum_fissureness_indices = [];
        return;
    end

    % Find a rotation matrix
    [eigv, m] = GetRotationMatrix(indices_for_model, image_size);

    [x_all, y_all, z_all] = MimImageCoordinateUtilities.FastInd2sub(image_size, candidate_indices);
    X_all = [x_all(:), y_all(:), z_all(:)]';

    % Transform to new basis
    em_all = X_all - m(:, ones(1, size(X_all, 2)));
    em_all = eigv*em_all;

    % We allocate each point to a bin according to the x-y coordinates in the
    % transformed domain
    x1_all = em_all(1, :);
    y1_all = em_all(2, :);

    % Get the coordinates of each of the model points
    [x_model, y_model, z_model] = MimImageCoordinateUtilities.FastInd2sub(image_size, indices_for_model);
    X_model = [x_model(:), y_model(:), z_model(:)]';

    % Transform to new basis
    em_model = X_model - m(:, ones(1, size(X_model, 2)));
    em_model = eigv*em_model;

    % We allocate each point to a bin according to the x-y coordinates in the
    % transformed domain
    x1_model = em_model(1, :);
    y1_model = em_model(2, :);

    % Project the candidate points onto the fissure plane, and remove those that
    % are not within the convex hull formed by the model points on the plane. In
    % effect, we are using the model coordinates (the initial guess) as a bound
    % for the x-y coordinates used to construct the fissure surface.
    xy_coords_model = [x1_model', y1_model'];
    dt = DelaunayTri(xy_coords_model);
    simplex_index = pointLocation(dt, x1_all', y1_all');
    is_valid = ~isnan(simplex_index);
    candidate_indices_ok = candidate_indices(is_valid);

    if isempty(candidate_indices_ok)
        maximum_fissureness_indices = [];
        return;
    end

    % Turn into x,y,z vectors
    [x_all, y_all, z_all] = MimImageCoordinateUtilities.FastInd2sub(image_size, candidate_indices_ok);
    X_all = [x_all(:), y_all(:), z_all(:)]';

    % Transform to new basis
    em_all = X_all - m(:, ones(1, size(X_all, 2)));
    em_all = eigv*em_all;

    % We allocate each point to a bin according to the x-y coordinates in the
    % transformed domain
    x1_all = em_all(1, :);
    y1_all = em_all(2, :);

    [maximum_fissureness_indices, all_maxima] = SortIntoBins(x1_all, y1_all, candidate_indices_ok, fissureness, bin_size, image_size);

    ref_image(all_maxima) = 6;
    ref_image(maximum_fissureness_indices) = 3;

    % Remove non-connected points (outliers)
    [x_all, y_all, z_all] = MimImageCoordinateUtilities.FastInd2sub(image_size, maximum_fissureness_indices);
    points_coords = ArraysToPoints(x_all(:), y_all(:), z_all(:));

    min_dilate_size_mm = 2.5;
    max_dilate_size_mm = 2.5;

    dilate_size_mm = min_dilate_size_mm;

    calculate_again = true;
    target_number_of_points = round(numel(x_all)/2);

    % Perform the removal of connected points with a variable dilation size
    while calculate_again
        points_coords_new = RemoveNonConnectedPoints(points_coords, lung_segmentation, image_size, dilate_size_mm);
        number_of_found_points = numel(points_coords_new);
        if (number_of_found_points >= target_number_of_points) || (dilate_size_mm > max_dilate_size_mm)
            calculate_again = false;
        else
            dilate_size_mm = dilate_size_mm + 0.5;
        end

    end

    [x_all, y_all, z_all] = PointsToArrays(points_coords_new);
    maximum_fissureness_indices = MimImageCoordinateUtilities.FastSub2ind(image_size, x_all(:), y_all(:), z_all(:));
    ref_image(maximum_fissureness_indices) = 1;

end

function [rot_matrix, m] = GetRotationMatrix(indices, image_size)
    [x, y, z] = ind2sub(image_size, indices);
    X = [x, y, z]';

    % Find a suitable basis for these points using PCA
    m = mean(X, 2);
    em = X - m(:, ones(1, size(X, 2)));
    eigv = pts_pca(em);

    rot_matrix = eigv';
end

function [maximum_fissureness_indices, all_maxima] = SortIntoBins(x1, y1, candidate_indices, fissureness, bin_size, image_size)
    binx = floor(x1/bin_size);
    binx = binx - min(binx);
    biny = floor(y1/bin_size);
    biny = biny - min(biny);

    % bin is an array containing the bin to which each element of indices has been allocated
    bin = binx + (max(binx)+1)*biny;

    % Now sort the indices and bins by fissureness
    fissureness_at_indices = fissureness.RawImage(candidate_indices);

    is_maxima = GetMaxima(candidate_indices, bin, fissureness_at_indices, fissureness, image_size);
    all_maxima = candidate_indices(is_maxima);
    bin = bin(is_maxima);
    candidate_indices = candidate_indices(is_maxima);

    fissureness_at_indices = fissureness.RawImage(candidate_indices);

    [~, sorted_indices] = sort(fissureness_at_indices, 'descend');
    indices_sorted_by_fissureness = candidate_indices(sorted_indices);
    bins_sorted_by_fissureness = bin(sorted_indices);

    % Use unique to obtain the first indices corresponding to each bin - this
    % will be the point with the largest fissureness
    maximum_fissureness_indices = [];
    for selection_index = 1 : 1
        [~, bin_indices, ~] = unique(bins_sorted_by_fissureness, 'first');
        maximum_fissureness_indices = [maximum_fissureness_indices; indices_sorted_by_fissureness(bin_indices)];

        bins_sorted_by_fissureness(bin_indices) = [];
        indices_sorted_by_fissureness(bin_indices) = [];
    end

    maximum_fissureness_indices = maximum_fissureness_indices(fissureness.RawImage(maximum_fissureness_indices) > 0);

end

function is_maxima = GetMaxima(candidate_indices, bin, fissureness_at_indices, fissureness, image_size)
    bin_allocation = zeros(image_size, 'int32');
    bin_allocation(candidate_indices) = bin;

    [linear_offsets, ~] = MimImageCoordinateUtilities.GetLinearOffsets(image_size);
    neighbours = repmat(int32(candidate_indices), 1, 6) + repmat(int32(linear_offsets), length(candidate_indices), 1);
    fissureness_neighbours = fissureness.RawImage(neighbours);
    bins_neighbours = bin_allocation(neighbours);
    bins_match = bins_neighbours == repmat(bin', 1, 6);
    fissureness_centre = repmat(fissureness_at_indices, 1, 6);

    % Force neighbouring points in different bins to have a lower fissureness
    fissureness_neighbours(~bins_match) = fissureness_centre(~bins_match) - 1;
    fissureness_less = fissureness_neighbours < fissureness_centre;
    sums = sum(fissureness_less, 2);
    is_maxima = sums == 6;
end

function points_coords = RemoveNonConnectedPoints(points_coords, template_image, image_size, dilate_size_mm)
    voxel_volume = prod(template_image.VoxelSize);
    min_component_size_mm3 = 300;
    min_component_size_voxels = round(min_component_size_mm3/voxel_volume);

    image_mask = false(image_size);
    indices = sub2ind(image_size, points_coords(:,1), points_coords(:,2), points_coords(:,3));
    image_mask(indices) = true;

    tci = template_image.BlankCopy;
    tci.ChangeRawImage(image_mask);
    tci.BinaryMorph(@imdilate, dilate_size_mm);
    connected_image = tci.RawImage;

    connected_components = bwconncomp(connected_image, 26);
    num_components = connected_components.NumObjects;
    for component = 1 : num_components
        pixels = connected_components.PixelIdxList{component};
        if length(pixels) < min_component_size_voxels
            connected_image(pixels) = false;
        end
    end
    image_mask = image_mask & connected_image;
    indices = find(image_mask(:));
    [i2_coords, j2_coords, k2_coords] = ind2sub(image_size, indices);
    points_coords = ArraysToPoints(i2_coords, j2_coords, k2_coords);
end

function points = ArraysToPoints(i, j, k)
    points = [i(:), j(:), k(:), ones(size(i(:)))];
end

function [is, js, ks] = PointsToArrays(points)
    is = points(:, 1);
    js = points(:, 2);
    ks = points(:, 3);
end

author: tomdoel

repo: pulmonarytoolkit


text:
% Function of Msg(tri)
function a = msg1(t,Ta)

t1 = -0.02:1.e-4:0;
t2 = 0:1.e-4:0.02;

m1 = 1 - abs((t1+Ta)/Ta);
m1 = [zeros([1 200]),m1,zeros([1 400])];
m2 = 1 - abs((t2-Ta)/Ta);
m2 = [zeros([1 400]),m2,zeros([1 200])];

a = m1 - m2;

end


text:
%SETUP  Compile the SPGL1 MEX interfaces

%   setup.m
%   $Id: spgsetup.m 1019 2008-06-17 19:33:48Z mpf $
%
%   ----------------------------------------------------------------------
%   This file is part of SPGL1 (Spectral Projected Gradient for L1).
%
%   Copyright (C) 2007 Ewout van den Berg and Michael P. Friedlander,
%   Department of Computer Science, University of British Columbia, Canada.
%   All rights reserved. E-mail: <{ewout78,mpf}@cs.ubc.ca>.
%
%   SPGL1 is free software; you can redistribute it and/or modify it
%   under the terms of the GNU Lesser General Public License as
%   published by the Free Software Foundation; either version 2.1 of the
%   License, or (at your option) any later version.
%
%   SPGL1 is distributed in the hope that it will be useful, but WITHOUT
%   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
%   or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
%   Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with SPGL1; if not, write to the Free Software
%   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
%   USA
%   ----------------------------------------------------------------------

root = pwd;
try
    cd('private')
    mex oneProjectorMex.c oneProjectorCore.c heap.c -output oneProjectorMex -DNDEBUG
    fprintf('Successfully compiled oneProjector.\n');
    cd(root)
catch err
    cd(root)
    fprintf('Could not compile oneProjector.');
    fprintf('You can still use the slower ".m" version.');
    rethrow(err);
end

try
    cd('private')
    mex productBMex.c -DNDEBUG
    fprintf('Successfully compiled productB.\n');
    cd(root)
catch err
    cd(root)
    fprintf('Could not compile productB.');
    rethrow(err);
end

author: mpf


text:
function mean = discrete_mean ( a, b )

%*****************************************************************************80
%
%% DISCRETE_MEAN evaluates the mean of the Discrete PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer A, the number of probabilities assigned.
%
%    Input, real B(A), the relative probabilities of
%    outcomes 1 through A.  Each entry must be nonnegative.
%
%    Output, real MEAN, the mean of the PDF.
%
  b_sum = sum ( b(1:a) );

  mean = 0.0;
  for j = 1 : a
    mean = mean + j * b(j);
  end

  mean = mean / b_sum;

  return
end


text:
function failed=demo_blockproc_header(demo_name,demo_nargin)
failed = 0;
if demo_nargin<1
   fprintf(['\n%s:\nTo run the demo, use one of the following:\n\n',...
          '%s(''gspi.wav'') to play gspi.wav (any wav file will do).\n',...
          '%s(''dialog'') to choose the wav file via file chooser dialog GUI.\n',...
          '%s(f,''fs'',fs) to play from a column vector f using sampling frequency fs.\n',...
          '%s(''playrec'') to record from a mic and play simultaneously.\n\n',...
          'Avalable input and output devices can be listed by |blockdevices|.\n',...
          'Particular device can be chosen by passing additional key-value pair ''devid'',devid.\n',...
          'Output channels of the device cen be selected by additional key-value pair ''playch'',[ch1,ch2].\n',...
          'Input channels of the device cen be selected by additional key-value pair ''recch'',[ch1].\n\n',...
          ]...
          ,upper(demo_name),demo_name,demo_name,demo_name,demo_name);
    failed=1;
end

try
   playrec('isInitialised');
catch
   error('%s: playrec or portaudio are not properly compiled. ',demo_name);
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
classdef gammaSections < ODFSections

  properties
    gamma
    sR
    maxGamma
  end

  methods

    function oS = gammaSections(varargin)
      % gamma sections for ODF and orientation plotting
      %
      % Syntax
      %
      %   oS = gammaSections(cs1,cs2,'sections',5)
      %   oS = gammaSections(cs1,cs2,'gamma',(0:15:90)*degree)
      %
      % Input
      %  cs1, cs2 - @crystalSymmetry, @specimenSymmetry
      %
      % Options
      %  sections - number of sections
      %  gamma - explicite section values
      %

      oS = oS@ODFSections(varargin{:});

      % get fundamental plotting region
      [alpha,beta,oS.maxGamma] = fundamentalRegionEuler(oS.CS1,oS.CS2,varargin{:}); %#ok<*PROP>
      oS.sR = sphericalRegion('maxTheta',beta,'maxRho',alpha);

      % get sections
      nsec = get_option(varargin,'sections',6);
      oS.gamma = linspace(0,oS.maxGamma,nsec+1);
      oS.gamma(end) = [];
      oS.gamma = get_option(varargin,'gamma',oS.gamma,'double');

      oS.updateTol(oS.gamma);
    end

    function ori = makeGrid(oS,varargin)

      oS.plotGrid = plotS2Grid(oS.sR,varargin{:});
      oS.gridSize = (0:numel(oS.gamma)) * length(oS.plotGrid);
      alpha = repmat(oS.plotGrid.rho,[1,1,numel(oS.gamma)]);
      beta = repmat(oS.plotGrid.theta,[1,1,numel(oS.gamma)]);
      gamma = repmat(reshape(oS.gamma,1,1,[]),[size(oS.plotGrid) 1]);

      ori = orientation.byEuler(alpha,beta,gamma,'ZYZ',oS.CS,oS.SS);

    end

    function ori = quiverGrid(oS,varargin)

      maxbeta = oS.sR.thetaMax;
      maxalpha = oS.sR.rhoMax;
      res = get_option(varargin,'resolution',15*degree);

      [alpha,beta,gamma] = meshgrid(res/2:res:maxalpha,res/2:res:maxbeta,oS.gamma);

      ori = orientation.byEuler(alpha,beta,gamma,'ZYZ',oS.CS,oS.SS);

    end

    function n = numSections(oS)
      n = numel(oS.gamma);
    end

    function [S2Pos,secPos] = project(oS,ori,varargin)

      % maybe this can be done more efficiently
      ori = ori.symmetrise('proper').';
      [alpha,beta,gamma] = Euler(ori,'ZYZ'); %#ok<*PROPLC>

      secPos = oS.secList(gamma,oS.gamma);
      S2Pos = vector3d.byPolar(beta,alpha);

    end

    function ori = iproject(oS,alpha,beta,igamma)
      ori = orientation.byEuler(alpha,beta,oS.gamma(igamma),'ZYZ',oS.CS,oS.SS);
    end

    function h = plotSection(oS,ax,sec,v,data,varargin)

      % plot data
      h = plot(v,data{:},oS.sR,'TR',[int2str(oS.gamma(sec)./degree),'^\circ'],...
        'parent',ax,'projection','plain','xAxisDirection','east',...
        'xlabel','$\alpha$','ylabel','$\beta$','dynamicMarkerSize',...
        'zAxisDirection','intoPlane',varargin{:},'doNotDraw');

    end

    function h = quiverSection(oS,ax,sec,v,data,varargin)

      % translate rotational data into tangential data
      if iscell(data) && isa(data{1},'quaternion')
        [v2,sec2] = project(oS,data{1},'noSymmetry','preserveOrder');
        data{1} = v2 - v;
        data{1}(sec2 ~= sec)=NaN;
      end
      if check_option(varargin,'normalize'), data{1} = normalize(data{1}); end

      % plot data
      h = quiver(v,data{:},oS.sR,'TR',[int2str(oS.gamma(sec)./degree),'^\circ'],...
      'parent',ax,'projection','plain','xAxisDirection','east',...
        'xlabel','$\alpha$','ylabel','$\beta$',...
        'zAxisDirection','intoPlane',varargin{:},'doNotDraw');

    end

  end
end

author: mtex-toolbox

repo: mtex

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


text:
function value = p28_autonomous ( )

%*****************************************************************************80
%
%% P28_AUTONOMOUS reports whether problem P28 is autonomous.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer VALUE, is 1 if the equation is autonomous.
%
  value = 1;

  return
end


text:
function nntest6(action,flag);
%NNTEST6 View  and Control Neural network.
%
%
%	See also SIMUFF.

% Koos j. den Oudsten, 1-20-99
% koos@phil.uu.nl
% Copyright (c) 1998-99 by KoosSoft vof
% $Revision: 0.4 $  $Date: 1999/01/27 22:29:28 $

global Xx w1 b1 w2 b2 hidden mlp

if nargin < 1
    action = 'start';
end

%On recursive calls get all necessary handles and data.
if ~strcmp(action,'start')
   childList = allchild(0);
   nn_fig = childList(find(childList == gcf));
   ud = gco;

%   popupvalue = get(ud.popup,'Value');

   switch action
      case {'Sw1', 'Sw2', 'Sw3', 'Sw4', 'Sw5', 'Sw6', 'Sb1', 'Sb2', 'Sb3', 'Sb4'}
         val=get(ud,'value');b=allchild(nn_fig);
	 c=findall(b,'tag',strrep(action,'S',''));set(c,'string',val);
		switch action
			case 'Sw1',w1(1)=val;
			case 'Sw2',w1(2)=val;
			case 'Sw3',w1(3)=val;
			case 'Sb1',b1(1)=val;
			case 'Sb2',b1(2)=val;
			case 'Sb3',b1(3)=val;
			case 'Sw4',w2(1)=val;
			case 'Sw5',w2(2)=val;
			case 'Sw6',w2(3)=val;
			case 'Sb4',b2(1)=val;
		end
	 action='plot';
      case 'train',
         x=0:0.05:5';plot(x,exp((-0.1.*x.^2)).*sin(4.*x));
      case {'w1', 'w2', 'w3', 'w4', 'w5', 'w6', 'b1', 'b2', 'b3', 'b4'}
         val=str2num(get(ud,'String'));
	 if val > 15  val=15; end; if val < -15  val =-15; end;
	 set(ud,'String',val);
	 b=allchild(nn_fig);
	 c=findall(b,'tag',strcat('S',action));set(c,'value',val);
		switch action
			case 'w1',w1(1)=val;
			case 'w2',w1(2)=val;
			case 'w3',w1(3)=val;
			case 'b1',b1(1)=val;
			case 'b2',b1(2)=val;
			case 'b3',b1(3)=val;
			case 'w4',w2(1)=val;
			case 'w5',w2(2)=val;
			case 'w6',w2(3)=val;
			case 'b4',b2(1)=val;
		end
	 action='plot';
      case 'reset',
         close(gcf);action='start';
      case 'set',
	 if mlp == 1,
         w1=[0.72 -0.83 0.83]';b1=[-1.8 4 -0.15]';w2=[-2 -1.7 1.7];b2=[0];
	 else w1=[0.5 2.5 4.6]';b1=[1.65 1.35 1]';w2=[3.1 2.4 1.9];b2=[-1.65];
	 end;
	 action='update';
      case 'grid',
         grid;
      case 'mlprbf',
	 if mlp == 1, mlp=0;
	 set(nn_fig,'Name','Neural Network Viewer (rbf)                  by  Koos j. den Oudsten');
	 else mlp =1;set(nn_fig,'Name','Neural Network Viewer (mlp)                  by  Koos j. den Oudsten');
	 end;
	 action='plot';
      case 'hidden',
         if hidden ==0, hidden=1; else hidden =0; end;
	 action='plot';
      case 'hold',
         hold;
   end
 disp(action);disp(ud);% test

end
% -----------------------------------------------------------------------
if strcmp(action,'update')
	b=allchild(nn_fig);

	set(findall(b,'tag','Sw1'),'value',w1(1));
	set(findall(b,'tag','Sw2'),'value',w1(2));
	set(findall(b,'tag','Sw3'),'value',w1(3));
	set(findall(b,'tag','Sb1'),'value',b1(1));
	set(findall(b,'tag','Sb2'),'value',b1(2));
	set(findall(b,'tag','Sb3'),'value',b1(3));
	set(findall(b,'tag','Sw4'),'value',w2(1));
	set(findall(b,'tag','Sw5'),'value',w2(2));
	set(findall(b,'tag','Sw6'),'value',w2(3));
	set(findall(b,'tag','Sb4'),'value',b2(1));

	set(findall(b,'tag','w1'),'string',w1(1));
	set(findall(b,'tag','w2'),'string',w1(2));
	set(findall(b,'tag','w3'),'string',w1(3));
	set(findall(b,'tag','b1'),'string',b1(1));
	set(findall(b,'tag','b2'),'string',b1(2));
	set(findall(b,'tag','b3'),'string',b1(3));
	set(findall(b,'tag','w4'),'string',w2(1));
	set(findall(b,'tag','w5'),'string',w2(2));
	set(findall(b,'tag','w6'),'string',w2(3));
	set(findall(b,'tag','b4'),'string',b2(1));

action='plot';
end

if strcmp(action,'plot')
	set(nn_fig,'DefaultAxesColorOrder',[1 0 0;0 1 0;0 0 1;1 1 0;1 0 1;0 0 0;0 1 1]);
	if mlp ==1, [hid,out]=simuff(Xx,w1,b1,'tansig',w2,b2,'purelin');
	else [hid,out]=simurb(Xx,w1,b1,w2,b2);
	end
	if hidden ==1,plot(Xx,out,'k',Xx,0,'r',Xx,hid);
		else plot(Xx,out,'k',Xx,0);
	end
	axis([0 5 -2 2]);axis manual;
end

%  Initialize all GUI objects.
if strcmp(action,'start')
% Set positions of graphic objects

Xx=0:0.05:5;w1=[0 0 0]';b1=[0 0 0]';w2=[0 0 0];b2=[0];% init weights en biases
hidden =0;mlp=1;% mlp mode

a = figure('Color',[0.8 0.8 0.8], ...
	'Position',[23 27 756 504], ...
'Name','Neural Network Viewer (mlp)                  by  Koos j. den Oudsten', ...
        'NumberTitle','off','Tag','NN-fig');
b = uimenu('Parent',a, ...
	'Label','Options', ...
	'Tag','NNmenu');
uimenu(b,'label','mlp <-> rbf', 'callback','nntest6(''mlprbf'')');
uimenu(b,'label','grid', 'callback','nntest6(''grid'')');
uimenu(b,'label','test', 'callback','nntest6(''test'')');
uimenu(b,'label','train','callback','nntest6(''train'')');
uimenu(b,'label','reset','callback','nntest6(''reset'')');
uimenu(b,'label','hold', 'callback','nntest6(''hold'')');
b = axes('Parent',a, ...
	'Units','points', ...
	'CameraUpVector',[0 1 0], ...
	'CameraUpVectorMode','manual', ...
	'Color',[1 1 1], ...
	'Position',[40.9655 125.379 410.897 180.621], ...
	'Tag','NetFigure', ...
	'XColor',[0 0 0], ...
	'YColor',[0 0 0], ...
	'ZColor',[0 0 0]);
c = text('Parent',b, ...
	'Color',[0 0 0], ...
	'HandleVisibility','callback', ...
	'HorizontalAlignment','center', ...
	'Position',[0.498489 -0.0965517 0], ...
	'Tag','Axes1Text4', ...
	'VerticalAlignment','cap');
set(get(c,'Parent'),'XLabel',c);
c = text('Parent',b, ...
	'Color',[0 0 0], ...
	'HandleVisibility','callback', ...
	'HorizontalAlignment','center', ...
	'Position',[-0.0528701 0.496552 0], ...
	'Rotation',90, ...
	'Tag','Axes1Text3', ...
	'VerticalAlignment','baseline');
set(get(c,'Parent'),'YLabel',c);
c = text('Parent',b, ...
	'Color',[0 0 0], ...
	'HandleVisibility','callback', ...
	'HorizontalAlignment','right', ...
	'Position',[-0.0996979 1.03448 0], ...
	'Tag','Axes1Text2', ...
	'Visible','off');
set(get(c,'Parent'),'ZLabel',c);
c = text('Parent',b, ...
	'Color',[0 0 0], ...
	'HandleVisibility','callback', ...
	'HorizontalAlignment','center', ...
	'Position',[0.498489 1.02759 0], ...
	'Tag','Axes1Text1', ...
	'VerticalAlignment','bottom');
set(get(c,'Parent'),'Title',c);
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','nntest6(''Sw1'')', ...
	'Position',[22.3448 31.1724 16.1379 63.931], ...
	'Style','slider','SliderStep',[0.0025 0.01], ...
	'Tag','Sw1','Min',-15.0,'Max',15.0, ...
	'Value',0.0);
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','nntest6(''Sw2'')', ...
	'Position',[63.5862 31.1724 16.1379 63.931], ...
	'Style','slider','SliderStep',[0.0025 0.01], ...
	'Tag','Sw2','Min',-15.0,'Max',15.0, ...
	'Value',0.0);
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','nntest6(''Sw3'')', ...
	'Position',[104.828 31.1724 16.1379 63.931], ...
	'Style','slider','SliderStep',[0.0025 0.01], ...
	'Tag','Sw3','Min',-15.0,'Max',15.0, ...
	'Value',0.0);
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','nntest6(''Sb1'')', ...
	'Position',[146.069 31.1724 16.1379 63.931], ...
	'Style','slider','SliderStep',[0.0025 0.01], ...
	'Tag','Sb1','Min',-15.0,'Max',15.0, ...
	'Value',0.0);
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','nntest6(''Sb2'')', ...
	'Position',[187.621 31.1724 16.1379 63.931], ...
	'Style','slider','SliderStep',[0.0025 0.01], ...
	'Tag','Sb2','Min',-15.0,'Max',15.0, ...
	'Value',0.0);
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','nntest6(''Sb3'')', ...
	'Position',[229.483 31.1724 16.1379 63.931], ...
	'Style','slider','SliderStep',[0.0025 0.01], ...
	'Tag','Sb3','Min',-15.0,'Max',15.0, ...
	'Value',0.0);
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','nntest6(''Sw4'')', ...
	'Position',[307.241 31.1724 16.1379 63.931], ...
	'Style','slider','SliderStep',[0.0025 0.01], ...
	'Tag','Sw4','Min',-15.0,'Max',15.0, ...
	'Value',0.0);
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','nntest6(''Sw5'')', ...
	'Position',[348.483 31.1724 16.1379 63.931], ...
	'Style','slider','SliderStep',[0.0025 0.01], ...
	'Tag','Sw5','Min',-15.0,'Max',15.0, ...
	'Value',0.0);
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','nntest6(''Sw6'')', ...
	'Position',[389.724 31.1724 16.1379 63.931], ...
	'Style','slider','SliderStep',[0.0025 0.01], ...
	'Tag','Sw6','Min',-15.0,'Max',15.0, ...
	'Value',0.0);
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','nntest6(''Sb4'')', ...
	'Position',[430.966 31.1724 16.1379 63.931], ...
	'Style','slider','SliderStep',[0.0025 0.01], ...
	'Tag','Sb4','Min',-15.0,'Max',15.0, ...
	'Value',0.0);
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','nntest6(''w1'')', ...
	'BackgroundColor',[1 .7 .7], ...
	'Position',[11.7931 17.3793 37.2414 12.4138], ...
	'Style','edit','String','w1', ...
	'Tag','w1');
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','nntest6(''w2'')', ...
	'BackgroundColor',[.7 1 .7], ...
	'Position',[53.0345 17.3793 37.2414 12.4138], ...
	'Style','edit','String','w2', ...
	'Tag','w2');
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','nntest6(''w3'')', ...
	'BackgroundColor',[.7 .7 1], ...
	'Position',[94.2759 17.3793 37.2414 12.4138], ...
	'Style','edit','String','w3', ...
	'Tag','w3');
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','nntest6(''b1'')', ...
	'BackgroundColor',[1 .7 .7], ...
	'Position',[135.517 17.3793 37.2414 12.4138], ...
	'Style','edit','String','b1', ...
	'Tag','b1');
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','nntest6(''b2'')', ...
	'BackgroundColor',[.7 1 .7], ...
	'Position',[176.759 16.7586 37.8621 13.0345], ...
	'Style','edit','String','b2', ...
	'Tag','b2');
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','nntest6(''b3'')', ...
	'BackgroundColor',[.7 .7 1], ...
	'Position',[218.621 16.7586 37.8621 13.0345], ...
	'Style','edit','String','b3', ...
	'Tag','b3');
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','nntest6(''w4'')', ...
	'BackgroundColor',[1 1 1], ...
	'Position',[296.69 17.3793 37.2414 12.4138], ...
	'Style','edit','String','w4', ...
	'Tag','w4');
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','nntest6(''w5'')', ...
	'BackgroundColor',[1 1 1], ...
	'Position',[337.931 17.3793 37.2414 12.4138], ...
	'Style','edit','String','w5', ...
	'Tag','w5');
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','nntest6(''w6'')', ...
	'BackgroundColor',[1 1 1], ...
	'Position',[379.172 17.3793 37.2414 12.4138], ...
	'Style','edit','String','w6', ...
	'Tag','w6');
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','nntest6(''b4'')', ...
	'BackgroundColor',[1 1 1], ...
	'Position',[420.414 16.7586 37.2414 13.0345], ...
	'Style','edit','String','b4', ...
	'Tag','b4');
b = uicontrol('Parent',a, ...
	'Units','points','string','hold', ...
	'Callback','nntest6(''hold'')', ...
	'Position',[260 80 34 16], ...
	'Tag','button1');
b = uicontrol('Parent',a, ...
	'Units','points','string','grid', ...
	'Callback','nntest6(''grid'')', ...
	'Position',[260 60 34 16], ...
	'Tag','button2');
b = uicontrol('Parent',a, ...
	'Units','points','string','set', ...
	'Callback','nntest6(''set'')', ...
	'Position',[260 40 34 16], ...
	'Tag','button3');
b = uicontrol('Parent',a, ...
	'Units','points','string','hidden', ...
	'Callback','nntest6(''hidden'')', ...
	'Position',[260 20 34 16], ...
	'Tag','button4');
end


text:
function Dnew = spm_eeg_remove_spikes(S)
% Use Will Pennys robust GLM code to remove 'spikes' from continuous data.
% Such spikes occur in EEG data recorded with the CTF MEG system at FIL
% due to some obscure electrical problem.
%
% FORMAT Dnew = spm_eeg_remove_spikes(S)
%
% S         - struct (optional)
% (optional) fields of S:
% S.D      - meeg object or filename
% S.logbf  - clean a block if log bayes factor in favour of spike model is
%            bigger than this (default - 3)
% S.hpf    - high-pass frequency above which to look for spikes (default 40 Hz)
% S.fast   - option to speed up the function by only using GLM if there is
%            threshold crossing ('yes', or check all the data with GLM - 'no')
% S.fasthresh - threshold for the fast option (in STD) - default 4
% S.trialbased - use trials in the data as they are ('yes') or break them
%              into sub-blocks ('no' - default)
% S.channels - channels to clean up (default 'gui' - brings up a GUI for
%             channel choice.
%
% Output:
% Dnew  - MEEG object with data cleaned of spikes.
%
%
% Disclaimer: this code is provided as an example and is not guaranteed to work
% with data on which it was not tested. If it does not work for you, feel
% free to improve it and contribute your improvements to the MEEGtools toolbox
% in SPM (http://www.fil.ion.ucl.ac.uk/spm)
%
% _______________________________________________________________________
% Copyright (C) 2008 Institute of Neurology, UCL

% Vladimir Litvak, Will Penny
% $Id: spm_eeg_remove_spikes.m 5640 2013-09-18 12:02:29Z vladimir $

if nargin == 0
    S = [];
end

% clean a block if log bayes factor in favour of spike model is bigger than this
if ~isfield(S, 'logbf'),             S.logbf = 3;                            end
% Have signals up to this freq in design matrix
if ~isfield(S, 'hpf'),               S.hpf = 40;                             end
% Use 1s blocks
if ~isfield(S, 'fast'),              S.fast='yes';                           end
if ~isfield(S, 'fasthresh'),         S.fasthresh=4;                          end
if ~isfield(S, 'trialbased'),        S.trialbased='no';                      end
if ~isfield(S, 'channels'),          S.channels='gui';                       end

try
    D = S.D;
catch
    [D,sts] = spm_select(1, 'mat', 'Select M/EEG mat file');
    if ~sts, D = []; return; end
    S.D = D;
end

D = spm_eeg_load(D);

if ~isfield(S, 'blocksize'),         S.blocksize = D.fsample;                end

spm('Pointer', 'Watch');drawnow;

if ~isfield(S, 'channels')
    chansel = D.indchantype({'MEEG', 'LFP'});
    S.channels = D.chanlabels(chansel);
elseif ischar(S.channels) && strcmp(S.channels, 'gui')
    chansel = listdlg('ListString', D.chanlabels, 'SelectionMode', 'multiple' ,'Name', 'Select channels' , 'ListSize', [400 300]);
    S.channels = D.chanlabels(chansel);
else
    chansel     = spm_match_str(D.chanlabels, S.channels);
end

nchan = length(chansel);    % number of channels

if strcmpi(S.trialbased, 'no') || (D.ntrials == 1)
    Nblocks=floor(D.nsamples./S.blocksize);
    block_size=S.blocksize;
    trialbased=0;
else
    Nblocks=D.ntrials;
    block_size=D.nsamples;
    trialbased=1;
end

% Make design matrix
K(1).row=[1:block_size];
K(1).RT=1/D.fsample;
K(1).HParam=1/S.hpf;
K = spm_filter(K);
X=K.X0;
X=[X ,ones(size(X,1),1)];

if ~trialbased && (D.ntrials > 1)
    rptnum = D.ntrials;
else
    rptnum = 1;
end

% generate new meeg object with new filenames
Dnew = clone(D, ['S' fnamedat(D)], [D.nchannels D.nsamples D.ntrials]);

for i = 1:D.nchannels
    Dnew(i, :, :) = D(i, :, :);
end

for r = 1:rptnum
    for ch=1:nchan
        for n=1:Nblocks,
            if trialbased
                lfp=squeeze(D(chansel(ch), :, n));
            else
                si=(n-1)*block_size+1;
                ei=si+block_size-1;
                lfp=squeeze(D(chansel(ch), si:ei, r));
            end
            m_lfp=mean(lfp);
            s_lfp=std(lfp);
            if s_lfp==0
                continue;
            end

            lfp=(lfp-m_lfp)/s_lfp;

            if ~(strcmpi(S.fast, 'yes') && max(abs(lfp))<S.fasthresh)
                glm=spm_rglm(lfp,X,1);
                [rglm,lfp_clean]=spm_rglm(lfp,X,2);
                if rglm.fm-S.logbf > glm.fm
                    if trialbased
                        Dnew(chansel(ch), :, n)=s_lfp*lfp_clean'+m_lfp;
                    else
                        Dnew(chansel(ch), si:ei, r) = s_lfp*lfp_clean'+m_lfp;
                    end
                    descr='spiky';
                else
                    descr='clean';
                end
                disp(sprintf('Channel %d out of %d, trial %d out of %d, block %d out of %d %s',ch,nchan, r, rptnum, n,Nblocks,descr));
            end
        end
    end
end

Dnew = history(Dnew, 'spm_eeg_remove_spikes', S);

save(Dnew);

spm('Pointer', 'Arrow');drawnow;

author: spm


text:
function [ a, b, c, d ] = plane_normal2imp_3d ( pp, normal )

%*****************************************************************************80
%
%% PLANE_NORMAL2IMP_3D converts a normal form plane to implicit form in 3D.
%
%  Discussion:
%
%    The normal form of a plane in 3D is
%
%      PP, a point on the plane, and
%      N, the unit normal to the plane.
%
%    The implicit form of a plane in 3D is
%
%      A * X + B * Y + C * Z + D = 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real PP(3), a point on the plane.
%
%    Input, real NORMAL(3), the unit normal vector to the plane.
%
%    Output, real A, B, C, D, the implicit plane parameters.
%
  a = normal(1);
  b = normal(2);
  c = normal(3);
  d = - a * pp(1) - b * pp(2) - c * pp(3);

  return
end


text:
function [dy,dh] = checkgrad(X_struct, f, e, varargin)

% checkgrad checks the derivatives in a function, by comparing them to finite
% differences approximations. The partial derivatives and the approximation
% are printed and the norm of the diffrence divided by the norm of the sum is
% returned as an indication of accuracy.
%
% usage: checkgrad('f', X, e, P1, P2, ...)
%
% where X is the argument and e is the small perturbation used for the finite
% differences. and the P1, P2, ... are optional additional parameters which
% get passed to f. The function f should be of the type
%
% [fX, dfX] = f(X, P1, P2, ...)
%
% where fX is the function value and dfX is a vector of partial derivatives.
%
% d - norm error
% dy - analytic
% dh - numerical
%
% Carl Edward Rasmussen, 2001-08-01.

[X_vec, X_template] = struct2vector(X_struct);

% [y, dy] = feval(f, vector2struct(X_vec + val*grad, X_template), varargin{:});
tic
[y, dy_struct] = feval(f, X_struct, varargin{:});
y = sum(y);
dy = struct2vector(dy_struct);
time1 = toc;

% keyboard

fig = figure;
dh = nan(length(X_vec),1) ;
t = cputime;
n_char = 0;
time_step = 5;
DIRECTION = 'forwards';
% DIRECTION = 'backwards';
% if length(X_vec) > 50000
%
% %   SKIP = 10;
% %   START = 1;%14400;%1;
%
%   SKIP = 2000;
%   START = 1;%14400;%1;
%
% %   SKIP = 1000;
% %   START = 30000;%14400;%1;
% else
%   SKIP = 1;%100;
%   START = 1;%4300;%1;
% end

START = 1;
n = max(40, 30/time1) % Evalulate for ~15 seconds
SKIP = ceil(length(X_vec)/n);

% START = 20210;
% SKIP = 1;

if strcmp(DIRECTION, 'backwards');
  idx = length(X_vec):-SKIP:1;
elseif strcmp(DIRECTION, 'forwards');
  idx = START:SKIP:length(X_vec);
end

% idx = length(X_vec):-1:(length(X_vec)-8);
% DIRECTION = 'backwards';
% SKIP = 1;

for j = idx

  dx = zeros(length(X_vec),1);
  dx(j) = dx(j) + e;                               % perturb a single dimension

  [y2] = feval(f, vector2struct(X_vec + dx, X_template), varargin{:});
  [y1] = feval(f, vector2struct(X_vec - dx, X_template), varargin{:});

  y1 = sum(y1(:));
  y2 = sum(y2(:));

  dh(j) = (y2 - y1)/(2*e);

  if ((cputime - t) > time_step) || (j == idx(end))
    t = cputime;
    figure(fig)
    clf
    subplot(2,1,1);
    if strcmp(DIRECTION, 'backwards');
      plot(j:SKIP:length(dy), dy(j:SKIP:end), 'bx-'); hold on;
      plot(j:SKIP:length(dy), dh(j:SKIP:end), 'r-');
    elseif strcmp(DIRECTION, 'forwards');
      plot(START:SKIP:j, dy(START:SKIP:j), 'bx-'); hold on;
      plot(START:SKIP:j, dh(START:SKIP:j), 'r-');
    end
    legend('analytical', 'numerical');
    ylabel('gradient');
    subplot(2,1,2);
    if strcmp(DIRECTION, 'backwards');
      plot(j:SKIP:length(dy), dy(j:SKIP:end) - dh(j:SKIP:end), 'gx-');
    elseif strcmp(DIRECTION, 'forwards');
      plot(START:SKIP:j, dy(START:SKIP:j) - dh(START:SKIP:j), 'gx-');
    end
    ylabel('error')
    drawnow

    s = [num2str(100*j/length(X_vec), '%0.1f'), '%, '];
    if n_char + length(s) >= 80
      fprintf('\n')
      n_char = 0;
    end
    fprintf('%s', s);
    n_char = n_char + length(s);

  end
end
fprintf('\n');

% 'numerical:'
% dh(idx)'
%
% 'analytical:'
% dy(idx)'

'delta:'
err = dy(idx)' - dh(idx)';

'mult:'
mult = dy(idx)' ./ dh(idx)';

for p = [1, 5, 20, 50, 90, 95, 99, 99.9, 99.99]
  fprintf('%g percentile err: \t%e\n', p, prctile(abs(err), p));
end

% % disp([dy dh])                                          % print the two vectors
% d = norm(dh-dy)/norm(dh+dy);       % return norm of diff divided by norm of sum

repo: CategoryShapes


text:
classdef SymmetricTensorInverter < Inverter

    properties (Access = private)
        voigtTensor
        invVoigtTensor
    end

    methods (Access = public)

        function obj = SymmetricTensorInverter(tensor)
            obj.compute(tensor);
        end
    end

    methods (Access = protected)

        function computeInverse(obj)
            obj.transformTensor2Voigt()
            obj.makeInverseOfVoigtTensor()
            obj.transformVoigt2Tensor()
        end
    end

    methods (Access = private)

        function transformTensor2Voigt(obj)
            obj.voigtTensor = Tensor2VoigtConverter.convert(obj.tensor);
        end

        function makeInverseOfVoigtTensor(obj)
            obj.invVoigtTensor = Inverter.invert(obj.voigtTensor);
        end

        function transformVoigt2Tensor(obj)
            invT = obj.invVoigtTensor;
            obj.invertedTensor = Voigt2TensorConverter.convert(invT);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function hRD = plugin_rd_tx
% Reference design definition

%   Copyright 2014-2015 The MathWorks, Inc.

% Call the common reference design definition function
hRD = AnalogDevices.fmcomms2.common.plugin_rd('ZC706', 'Tx');
AnalogDevices.fmcomms2.zc706.add_io(hRD,'Tx');

author: analogdevicesinc