text:
function outarg = meshAngle(command, name, V);
%
% outarg = meshAngle([command], [name of angle], [Volume View]);
%
% Load/Save camera angles for viewing a mesh from canonical views.
% Angles are stored in the same directory as the mesh, in the file
% 'MeshAngles.mat'. These can be saved/restored for any meshes in that
% directory, so if both right and left hemishpheres are saved, it
% might be good to name it something like 'lateral-right' or
% 'dorsolateral-left'.
%
% Volume view: view w/ the mesh attached. [defaults to selected gray]
% command: 'load', 'save', or 'list'. [defaults to 'list']
% name: name of the mesh to save.
%
% Returns a different output argument depending on the command:
%   'list': returns a struct containing the saved angles. (If meshAngle is
%   called as meshAngle(V, 'list', 'quiet'), will only return in the
%   struct and not print to the command window.)
%   'load': returns the cRot value loaded.
%   'save': returns the path of the save file.
%
% Example: meshAngle(getSelectedGray, 'save', 'ventral') saves
% the view under the name 'ventral';
% meshAngle(getSelectedGray, 'load', 'ventral') restores it.
% meshAngle(getSelectedGray, 'list') lists the saved angles for this
% mesh, if any.
%
% ras, 10/2005.
if ieNotDefined('V'), V = getSelectedGray; end
if ieNotDefined('command'), command = 'list'; end

mesh = viewGet(V, 'mesh');

anatPath = getVAnatomyPath;
meshDir = fullfile(fileparts(anatPath), 'Meshes');

switch lower(command)
    case 'load',
        % get angles (and perform file check while we're at it)
        angles = meshAngle('list', 'quiet', V);
        angleNames = {angles.name};

        % if no name specified, get from dialog
        if ieNotDefined('name')
            dlg.fieldName = 'name';
            dlg.style = 'listbox';
            dlg.string = 'Load which angle?';
            dlg.list = angleNames;
            dlg.value = 1;
            resp = generaldlg2(dlg,'Load Camera Angle...');
            name = resp.name{1};
        end

        % check that the named angle has been stored in the file
        I = cellfind(angleNames, lower(name));
        if isempty(I)
            fprintf('Angle %s not found. Saved angles include:\n',name);
            for i = 1:length(angleNames)
                fprintf('%s \n', angleNames{i});
            end
            error('Selected angle not found');
        else
            cRot = angles(I).cRot;
            mrmSet(mesh,'camerarotation',cRot);
        end

        outarg = cRot;

    case 'save',
        cRot = mrmGet(mesh,'camerarotation');

        % if no name specified, get from dialog.
        if ieNotDefined('name')
            name = inputdlg({'Name of Angle to Save:'},...
                    'Save Camera Angle...',1,{'angle1'});
            if isempty(name{1}), return; end  % user canceled
            name = lower(name{1});
        end

        % load any angles that may exist.
        pth = fullfile(meshDir, 'MeshAngles.mat');
        if exist(pth,'file'),
            load(pth,'angles');
            angles(end+1).name = lower(name);
        else
            angles.name = lower(name);
        end
        angles(end).cRot = cRot;
        angles(end).saved = datestr(clock);

        % save the angles.
        if exist(pth,'file')
            save(pth,'angles','-append');
        else
            save(pth,'angles')
        end
        fprintf('Saved camera rotation as angles(%i) in %s.\n',...
                length(angles),pth);

        outarg = pth;

    case 'list',
        % check if the Angles.mat file exists
        pth = fullfile(meshDir, 'MeshAngles.mat');
        if ~exist(pth, 'file')
            msg = sprintf(['No angles have been stored for ' ...
                           'this segmentation: %s.'], mesh.path);
            myWarnDlg(msg);
            angles = [];
            return
        else
            load(pth,'angles');
        end

        if ieNotDefined('name') | ~isequal(lower(name),'quiet')
            disp('***********************')
            disp('Angles Saved: ')
            for i = 1:length(angles)
                fprintf('%s \n', angles(i).name);
            end
            disp('***********************')
        end

        outarg = angles;

    otherwise, help(mfilename); warning('Unknown command');
end

return


text:
% pet_transmission_example.m
% Example of reconstructing a PET attenuation map via penalized-likelihood
% estimation from real 2D PET transmission data.
%
% Copyright 2002-12-19, Jeff Fessler, University of Michigan

% read raw data
if ~isvar('yi'), printm 'raw data'
    yi = ir_get_data(fullfile('pet_trans_2d_sino','phan_trans.mat'));
    bi = ir_get_data(fullfile('pet_trans_2d_sino','phan_blank.mat'));
	im plc 2 2, im(1, yi, 'yi: transmission scan')
	im(2, bi, 'bi: blank scan')

	% system model for ecat921 PET scanner geometry
	ig = image_geom('nx', 128, 'dx', 0.421875, 'dy', 0.421875); % cm
%		'center_x', 0.5, 'center_y', -0.5, ...
	sg = sino_geom('par', 'nb', size(yi,1), 'na', size(yi,2), ...
		'dr', 0.3375, 'offset_r', 0.5, 'strip_width', 'dr', ...
		'orbit_start', -15);
prompt
end

% FBP image
if ~isvar('xfbp'), printm 'do FBP'
	kernel = gaussian_kernel(3);
%	kernel = [1];
	xfbp = tr_fbp(sg, ig, max(yi,1), bi, 0*bi, 'kernel', kernel);
% fix: correct for backproject scaling problem; probably incorrect!
%	xfbp = xfbp * na/pi * f.pixel_size / f.ray_spacing;
	xfbp = max(xfbp,0);
	im(3, xfbp, 'fbp'), cbar

	ig.mask = ig.circ(26, 21, 0, 3) > 0;
	im(4, ig.mask+6*xfbp, 'mask check')
prompt
end

% strip integral system matrix with mask for iterative reconstruction.
if ~isvar('A2'), printm 'A2'
	if 0 && has_mex_jf
		A2 = Gtomo2_wtmex(sg, ig); % preferable for speed
	else
		A2 = Gtomo2_strip(sg, ig, 'strip_width', sg.dr); % slower but universal
	end
end

if ~isvar('Ab'), printm 'make Ab' % block system object for ordered-subsets
	f.nblock = 5;
	Ab = Gblock(A2, f.nblock);
end

if ~isvar('R'), printm 'make R' % regularizer object
	f.l2b = 10.5;
	f.delta = 0.03;
%	f.pot = 'huber';
	f.pot = 'hyper3';
	R = Reg1(ig.mask, 'type_denom', 'matlab', ...
		'beta', 2^f.l2b, 'pot_arg', {f.pot, f.delta});
end

% matlab iterations
if ~isvar('xmat'), printm 'matlab T-PL-OS-SPS'
	f.niter = 8+1;
	f.niter = 16+1;
	f.pixmax = 0.4;
	xinit = max(xfbp,0);

	xmat = tpl_os_sps(xinit(ig.mask), Ab, yi, bi, [], R, ...
		f.niter, f.pixmax, 'pc');
	xmat = ig.embed(xmat);

	im clf, im(xmat, 'T-PL-OSPS iterations (0th is FBP)')
prompt
end

% nice figure for book
if 1
	im clf
	im(211, yi, 'sinogram yi')
	clim = [0 0.2];
	im(223, xfbp, 'FBP', clim), cbar
	im(224, xmat(:,:,end), 'Statistical', clim), cbar
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function ierror = setpart_check ( m, nsub, s, index )

%*****************************************************************************80
%
%% SETPART_CHECK checks a set partition.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer M, the number of elements of the set.
%    M must be positive.
%
%    Input, integer NSUB, the number of nonempty subsets into
%    which the set is partitioned.  1 <= NSUB <= M.
%
%    Input, integer INDEX(NSUB), lists the location in S of the
%    last element of each subset.  Thus, the elements of subset 1
%    are S(1) through S(INDEX(1)), the elements of subset 2
%    are S(INDEX(1)+1) through S(INDEX(2)) and so on.
%
%    Input, integer S(M), contains the integers from 1 to M,
%    grouped into subsets as described by INDEX.
%
%    Output, integer IERROR, error flag.
%    0, no error.
%    -I, the I-th element of INDEX is illegal.
%    +I, the I-th element of S is illegal.
%
  ierror = 0;
%
%  Check INDEX.
%
  imin = 0;
  for i = 1 : nsub
    if ( index(i) <= imin || m < index(i) )
      ierror = -i;
      return
    end
    imin = index(i);
  end
%
%  Check the elements of S.
%
  for i = 1 : nsub

    if ( s(i) <= 0 || m < s(i) )
      ierror = i;
      return
    end

    for j = 1 : i - 1
      if ( s(j) == s(i) )
        ierror = i;
        return
      end
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function adpcm_y = adpcm_encoder(raw_y)

% This m-file is based on the app note: AN643, Adaptive differential pulse
% code modulation using PICmicro microcontrollers, Microchip Technology
% Inc. The app note is avaialbe from www.microchip.com
% Example:  Y = wavread('test.wav');
%           y = adpcm_encoder(Y);
%           YY = adpcm_decode(y);

IndexTable = [-1, -1, -1, -1, 2, 4, 6, 8, -1, -1, -1, -1, 2, 4, 6, 8];

StepSizeTable = [7, 8, 9, 10, 11, 12, 13, 14, 16, 17, 19, 21, 23, 25, 28, 31, 34, 37, 41, 45, 50, 55, 60, 66, 73, 80, 88, 97, 107, 118, 130, 143, 157, 173, 190, 209, 230, 253, 279, 307, 337, 371, 408, 449, 494, 544, 598, 658, 724, 796, 876, 963, 1060, 1166, 1282, 1411, 1552, 1707, 1878, 2066, 2272, 2499, 2749, 3024, 3327, 3660, 4026, 4428, 4871, 5358, 5894, 6484, 7132, 7845, 8630, 9493, 10442, 11487, 12635, 13899, 15289, 16818, 18500, 20350, 22385, 24623, 27086, 29794, 32767];

prevsample = 0;
previndex = 1;

Ns = length(raw_y);
n = 1;

raw_y = 32767 * raw_y;          % 16-bit operation

while (n <= Ns)
    predsample = prevsample;
    index = previndex;
    step = StepSizeTable(index);

    diff = raw_y(n) - predsample;
    if (diff >= 0)
        code = 0;
    else
        code = 8;
        diff = -diff;
    end

    tempstep = step;
    if (diff >= tempstep)
        code = bitor(code, 4);
        diff = diff - tempstep;
    end
    tempstep = bitshift(tempstep, -1);
    if (diff >= tempstep)
        code = bitor(code, 2);
        diff = diff - tempstep;
    end
    tempstep = bitshift(tempstep, -1);
    if (diff >= tempstep)
        code = bitor(code, 1);
    end

    diffq = bitshift(step, -3);
    if (bitand(code, 4))
        diffq = diffq + step;
    end
    if (bitand(code, 2))
        diffq = diffq + bitshift(step, -1);
    end
    if (bitand(code, 1))
        diffq = diffq + bitshift(step, -2);
    end

    if (bitand(code, 8))
        predsample = predsample - diffq;
    else
        predsample = predsample + diffq;
    end

    if (predsample > 32767)
        predsample = 32767;
    elseif (predsample < -32768)
        predsample = -32768;
    end

    index = index + IndexTable(code+1);

    if (index < 1)
        index = 1;
    end
    if (index > 89)
        index = 89;
    end

    prevsample = predsample;
    previndex = index;

    adpcm_y(n) = bitand(code, 15);
    %adpcm_y(n) = code;
    n = n + 1;
end

author: Sable

repo: mcbench-benchmarks


text:
function imageCha = recon_main(obj, iRep, iAvg)
% main reconstruction function
% prepare constraints and select algorithm depending on image dimensionality
%
% input:
% obj       CS reconstruction object (holding all parameters)
% iRep      current repetition
% iAvg      current average
%
% output:
% imageCha  individual reconstructed channel images
%
% (c) Thomas Kuestner
% -------------------------------------------------------------------------

nCha = obj.measPara.dim(5);
nSlices = size(obj.kSpace,1);
nTime = obj.measPara.dim(4);
imageCha = cell(nSlices,nCha);

if(strcmp(obj.measPara.dimension,'2D'))
% lMask = cellfun(@(x) abs(x) > 0, all_KSpaces, 'UniformOutput', false);

    dispProgress('Slices', 0, nSlices);
    for iSli = 1:nSlices
        % prepare k-space
        % => k_y - k_x - t - cha
        kSpaceL = cell2mat(shiftdim(obj.kSpace(iSli,:,iRep,iAvg),-2));
        if(strcmp(obj.measPara.precision,'double'))
            kSpaceL = double(kSpaceL);
        end
        % same undersampling for all coils and fully sampled in k_x direction
        % => reduce 4D problem (k_y - k_x - t - cha) to 2D problem (k_y - t)
        obj.fullMask = abs(kSpaceL(:,:,:,1)) > 0; % k_y - k_x - t

        if(obj.measPara.oversampling{2,1})
            kSpaceL = fftshift(ifft(ifftshift(kSpaceL, 2),[],2),2); % -> k_y - x - t - cha
            kSpaceL = kSpaceL(:, obj.measPara.oversampling{1,1}, :, :); % anti-aliasing
            obj.fullMask = obj.fullMask(:, obj.measPara.oversampling{1,1}, :, :);
            kSpaceL = fftshift(fft(ifftshift(kSpaceL,2),[],2),2); % -> k_y - k_x - t -cha
            kSpaceL = kSpaceL .* repmat(obj.fullMask,[1 1 1 size(kSpaceL,4)]);
            obj.measPara.dim(2) = size(kSpaceL,2);
        end

        [nPha, nFreq, nTime, nCha] = size(kSpaceL);

        %     obj.fullMask = lMask{1}(:, obj.obj.measPara.oversampling{1}, :); % k_y - k_x - t
        %     mask = abs(sum(permute(fullMask,[3 1 2]),3)) > 0;

        if(iSli == 1)
            if(isempty(obj.calibSize))
                obj.calibSize = cell(1,nTime);
                reempty = true;
            else
                % adjust calibSize due to anti-aliasing
                obj.calibSize(2) = nFreq;
                reempty = false;
            end
        else
            if(reempty)
                obj.calibSize = cell(1,nTime);
            end
        end
        obj.kCalib = cell(1,nTime);
        obj.kernel = obj.kCalib;
        obj.kernelImg = obj.kCalib;
        ESP = obj.kCalib;
%         obj.A = obj.kCalib;

        % convolution with GRAPPA kernel in k-space -> multiplication in image
        % domain
        obj.method = 'fft';

        % prepare k-space => k_y - k_x - cha - t
        kSpaceL = permute(kSpaceL, [1 2 4 3]);

        fprintf('Performing calibration\n');
        for i=1:nTime

            % perform calibration
            % get calibration data from k-space
            if(iscell(obj.calibSize))
                % automatically determine calibration region
                obj.calibSize{i} = obj.calibrationSize(obj.fullMask(:,:,i)); % k_y - k_x (cell: t)
                if(~obj.flagToolbox)
                    obj.kCalib{i} = crop(kSpaceL(:,:,:,i), [obj.calibSize{i}, nCha]); % k_y - k_x - cha (cell: t)
                end
            else
                if(~obj.flagToolbox)
                    % fixed calibration region
                    plonk = [nPha, nFreq];
                    idx = cell(1,length(plonk));
                    for iInner=1:length(plonk)
                        idx{iInner} = floor(plonk(iInner)/2)+1+ceil(-obj.calibSize(iInner)/2):floor(plonk(iInner)/2)+ceil(obj.calibSize(iInner)/2);
                    end
                    obj.kCalib{i} = kSpaceL(idx{1},idx{2},:,i);
                end
            end

            if(~obj.flagToolbox)
    %             % deprecated
    %             obj.kernel{i} = zeros([obj.kernelSize, nCha, nCha]); % k_y - k_x - cha - cha (cell: t)
    %
    %             [AtA,obj.A{i}] = obj.corrMatrix(i);
    %             for n=1:nCha
    %                 obj.kernel{i}(:,:,:,n) = obj.calibrate(AtA,nCha,n);
    %             end

                % compute eigen-value maps
                [obj.kernel{i},S] = dat2Kernel(obj.kCalib{i},obj.kernelSize);
                idx = max(find(S >= S(1)*obj.eigThresh_k));

                % crop kernels and compute eigen-value decomposition in image space to get maps
                [M,W] = kernelEig(obj.kernel{i}(:,:,:,1:idx),[nPha,nFreq]);

                % Compute Soft-SENSE ESPIRiT Maps
                % crop sensitivity maps according to eigenvalues == 1
                % weigth with n maps of eigen-values

                maps = M(:,:,:,end-(obj.n_maps-1):end);

                % Weight the eigenvectors with soft-senses eigen-values
                weights = W(:,:,end-(obj.n_maps-1):end) ;
                weights = (weights - obj.eigThresh_im)./(1-obj.eigThresh_im).* (W(:,:,end-(obj.n_maps-1):end) > obj.eigThresh_im);
                weights = -cos(pi*weights)/2 + 1/2;

                % create and ESPIRiT operator
                ESP{i} = ESPIRiT(maps,weights);
            else
                if(ispc)
                    eval(sprintf('ESP{i} = bart(''ecalib -r %d:%d -k %d:%d -m %d -1'',kSpaceL(:,:,:,i));', obj.calibSize{i}(1), obj.calibSize{i}(2), obj.kernelSize(1), obj.kernelSize(2), obj.n_maps ));

                else
                    % create temporary dir if not yet existing
                    if(~exist([obj.path,filesep,'tmpESPIRiT'],'dir'))
                        mkdir([obj.path,filesep,'tmpESPIRiT']);
                    end
                    writecfl([obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], kSpaceL(:,:,:,i)); % k_y - k_x - cha
                    eval(sprintf('! ecalib -r %d:%d -k %d:%d -m %d %s %s', obj.calibSize{i}(1), obj.calibSize{i}(2), obj.kernelSize(1), obj.kernelSize(2), obj.n_maps, [obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], [obj.path,filesep,'tmpESPIRiT',filesep,'espiritmaps_',num2str(i)]));
                    ESP{i} = [obj.path,filesep,'tmpESPIRiT',filesep,'espiritmaps_',num2str(i)];
                end
            end
        end

        % ensure that empty entries are zero
%         mask = false(nPha,nFreq,nCha,nTime);
        mask = obj.fullMask;
%         for i=1:nTime
%             mask(:,:,:,i) = repmat(obj.fullMask(:,:,i),[1,1,nCha]);
%         end
        kSpaceL(~mask) = 0;

        if(strcmp(obj.solver,'L1') && ~obj.flagToolbox)
            if(strcmp(obj.trafo.trafoType,'wavelet_lab'))
                ssx = 2^ceil(log2(nPha));
                ssy = 2^ceil(log2(nFreq));
                ss = max(ssx, ssy);
                XOP = Wavelet(obj.trafo.waveletFilter,obj.trafo.waveletFilterSize,log2(ss/2^(obj.trafo.waveletStages)));
            end
        end

        image = zeros(nPha,nFreq,obj.n_maps,nTime);
        dispProgress('Time', 0, nTime);
        fprintf('ESPIRiT reconstruction\n');
        for i=1:nTime
            if(~obj.flagToolbox)
                if(strcmp(obj.solver,'L1'))
                    FT = p2DFT(mask(:,:,:,i),[nPha,nFreq,nCha]);
                    if(strcmp(obj.trafo.trafoType,'fft'))
                        XOP = FT;
                    end
                    image(:,:,:,i) = obj.cgL1ESPIRiT(kSpaceL(:,:,:,i), zeros(size(kSpaceL(:,:,:,i))), FT, ESP{i}, obj.iNINNER, XOP, obj.trafo.wavWeight, obj.splitWeight, obj.iNIterSplit); % => y-x-cha-t
                else
                    [~,image(:,:,:,i)] = obj.cgESPIRiT(kSpaceL(:,:,:,i),ESP{i}, obj.iNINNER, obj.calibTyk, zeros(size(kSpaceL(:,:,:,i))));
                end
            else
                if(ispc)
                    if(strcmp(obj.solver,'L1'))
                        eval(sprintf('image(:,:,:,i) = squeeze(bart(''sense -l1 -r %f'', kSpaceL(:,:,:,i), ESP{i}));', obj.trafo.wavWeight));
                    else
                        image(:,:,:,i) = squeeze(bart('sense -l2', kSpaceL(:,:,:,i), ESP{i}));
                    end
                else
                    fftdims = size(kSpaceL);
                    fftdims = fftdims(1:2);
                    if(any(mod(fftdims,2) ~= 0))
                        writecfl([obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], zpad(kSpaceL(:,:,:,i),size(kSpaceL(:,:,:,i)) + [mod(fftdims,2),0])); % k_y - k_z - k_x - cha
                    else
                        writecfl([obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], kSpaceL(:,:,:,i)); % k_y - k_x - cha
                    end
                    if(strcmp(obj.solver,'L1'))
                        eval(sprintf('! sense -l1 -r %f %s %s %s', obj.trafo.wavWeight, [obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], ESP{i}, [obj.path,filesep,'tmpESPIRiT',filesep,'imgOut']));
                    else
                        eval(sprintf('! sense -l2 %s %s %s', [obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], ESP{i}, [obj.path,filesep,'tmpESPIRiT',filesep,'imgOut']));
                    end
    %                 if(obj.n_maps > 1)
    %                     % combine ESPIRiT images into one
    %                    eval(sprintf('! rss 3 %s %s', [obj.path,filesep,'tmpESPIRiT',filesep,'imgOut'], [obj.path,filesep,'tmpESPIRiT',filesep,'imgOutRSS']));
    %                    delete([obj.path,filesep,'tmpESPIRiT',filesep,'imgOut']);
    %                    movefile([obj.path,filesep,'tmpESPIRiT',filesep,'imgOutRSS'], [obj.path,filesep,'tmpESPIRiT',filesep,'imgOut']);
    %                 end
                    image(:,:,:,i) = squeeze(readcfl([obj.path,filesep,'tmpESPIRiT',filesep,'imgOut']));
                end
            end
            dispProgress('Time',i/nTime);
        end
        dispProgress('Time','Close');

        imageCha(iSli,1:obj.n_maps) = shiftdim(mat2cell(permute(image,[1 2 4 3]),nPha,nFreq,nTime,ones(1,obj.n_maps)),2);
        [imageCha{iSli,obj.n_maps+1:end}] = deal(zeros(size(imageCha{iSli,1})));

        dispProgress('Slices', iSli/nSlices);
    end
    dispProgress('Slices', 'Close');

elseif(strcmp(obj.measPara.dimension,'3D'))

    % prepare k-space
    % => k_y - k_z - k_x - cha
    kSpaceL = cell2mat(shiftdim(obj.kSpace(:,:,iRep,iAvg),-2));
    kSpaceL = permute(kSpaceL,[1 3 2 4]);
    if(strcmp(obj.measPara.precision,'double'))
        kSpaceL = double(kSpaceL);
    end
    % same undersampling for all coils and fully sampled in k_x direction
    obj.fullMask = abs(kSpaceL) > 0; % k_y - k_z - k_x
    if(obj.measPara.oversampling{2,1})
        kSpaceL = ifftnshift(kSpaceL, 3); % -> k_y - k_z - x - cha
        kSpaceL = kSpaceL(:, :, obj.measPara.oversampling{1,1}, :); % anti-aliasing
        obj.fullMask = obj.fullMask(:,:, obj.measPara.oversampling{1,1}, :);
        kSpaceL = fftnshift(kSpaceL,3); % -> k_y - k_z - k_x - cha
        kSpaceL = kSpaceL .* obj.fullMask;
        obj.measPara.dim(2) = size(kSpaceL,3);
    end
    [nPha, nZ, nFreq, nCha] = size(kSpaceL);

    if(~isempty(obj.calibSize))
        % adjust calibSize due to anti-aliasing
        obj.calibSize(3) = nFreq;
    end

    % create a 3D kernel fomr the calibration region
    fprintf('Performing calibration\n');
    if(isempty(obj.calibSize))
        % automatically determine calibration region
        obj.calibSize = obj.calibrationSize(obj.fullMask(:,:,:,1)); % k_y - k_z - k_x
        if(size(kSpaceL,3) == 1), obj.calibSize = [obj.calibSize, 1]; obj.kernelSize = obj.kernelSize([1 3 2]); end; % k_y - k_x sparse!
        if(~obj.flagToolbox)
            obj.kCalib = crop(kSpaceL, [obj.calibSize, nCha]); % k_y - k_z - k_x - cha
        end
    else
        if(~obj.flagToolbox)
            % fixed calibration region
            plonk = [nPha, nZ, nFreq];
            idx = cell(1,length(plonk));
            for iInner=1:length(plonk)
                if(mod(obj.calibSize(iInner),2) == 0)
                    idx{iInner} = floor(plonk(iInner)/2)+1+ceil(-obj.calibSize(iInner)/2):floor(plonk(iInner)/2)+ceil(obj.calibSize(iInner)/2);
                else
                    idx{iInner} = floor(plonk(iInner)/2)+ceil(-obj.calibSize(iInner)/2):floor(plonk(iInner)/2)+ceil(obj.calibSize(iInner)/2)-1;
                end

                tmp = [idx{iInner}(1) <= 0, idx{iInner}(end) > plonk(iInner)];
                if(any(tmp))
                    if(all(tmp)), error('crop(): Both index out of bounds'); end;
                    hShift = [abs(idx{iInner}(1)) + 1, idx{iInner}(end) - plonk(iInner)];
                    op = {idx{iInner}(1), idx{iInner}(end); '+', '-'; '<', '>'; plonk(iInner) + 1, 0};
                    eval(sprintf('if(op{1,~helper} %s hShift(helper) %s %d), idx{iInner} = idx{iInner} %s hShift(helper); else idx{iInner} = idx{iInner}(idx{iInner} %s %d);end;',op{2,helper}, op{3,helper}, op{4,helper}, op{2,helper}, op{3,~helper}, op{4,~helper}));
                end
            end
            obj.kCalib = kSpaceL(idx{1},idx{2},idx{3},:); % k_y - k_z - k_x - cha
        end
    end

    if(~obj.flagToolbox)
        % windowing for cut-out region
    %     [~,w1] = window2D(prop.window.type,size(obj.kCalib,1),prop.window.windowOpt{1},prop.window.windowOpt{2});
    %     [~,w2] = window2D(prop.window.type,size(obj.kCalib,2),prop.window.windowOpt{1},prop.window.windowOpt{2});
    %     window =  w1(:) * w2(:).';
        if(obj.window.windowingOn)
            window = windowND(obj.window.type,[size(obj.kCalib,1), size(obj.kCalib,2)],obj.window.windowOpt{1},obj.window.windowOpt{2});
            obj.kCalib = obj.kCalib .* repmat(window,[1 1 size(obj.kCalib,3) nCha]);
        end

        % compute eigen-value maps
        [hybridKernel,S] = dat3Kernel(obj.kCalib,obj.kernelSize);
        hybridKernel = ifftnshift(zpad(hybridKernel(:,:,end:-1:1,:,:),[obj.kernelSize(1),obj.kernelSize(2),nFreq + obj.kernelSize(3)-1,nCha,size(hybridKernel,5)]),3);
        idx = max(find(S >= S(1)*obj.eigThresh_k));

        % compute inverse fourier transform along readout direction
        kSpaceL = ifftnshift(kSpaceL,3); % k_y - k_z - x - cha

        % convolution with GRAPPA kernel in k-space -> multiplication in image
        % domain
        obj.method = 'fft';

        % prepare wavelet transformation
        if(strcmp(obj.solver,'L1'))
            if(strcmp(obj.trafo.trafoType,'wavelet_lab'))
                ssx = 2^ceil(log2(nPha));
                ssy = 2^ceil(log2(nZ));
                ss = max(ssx, ssy);
                XOP = Wavelet(obj.trafo.waveletFilter,obj.trafo.waveletFilterSize,log2(ss/2^(obj.trafo.waveletStages)));
            end
        end

        image = zeros(nPha,nFreq,nZ,obj.n_maps);
    %     AtA_sub = zeros(prod(obj.kernelSize)*nCha); % subtract this matrix from AtA_mask, due to accumulation of previous samples (A'*A multiplication)
    %     A_sub = zeros((size(obj.kCalib,1)-obj.kernelSize(1)+1) * (size(obj.kCalib,2)-obj.kernelSize(2)+1) * (size(obj.kCalib,3)-obj.kernelSize(3)+1),prod(obj.kernelSize)*nCha);

        fprintf('ESPIRiT reconstruction\n');
        dispProgress('Frequency',0,nFreq);
        for iFreq=1:nFreq

            % extract 2D kernel from complete 3D kernel
            obj.kernel = squeeze(hybridKernel(:,:,iFreq,:,:));

            % crop kernels and compute eigen-value decomposition in image space to get maps
            [M,W] = kernelEig(obj.kernel(:,:,:,1:idx),[nPha,nZ]);

            % Compute Soft-SENSE ESPIRiT Maps
            % crop sensitivity maps according to eigenvalues == 1
            % weigth with n maps of eigen-values

            maps = M(:,:,:,end-(obj.n_maps-1):end);

            % Weight the eigenvectors with soft-senses eigen-values
            weights = W(:,:,end-(obj.n_maps-1):end);
            weights = (weights - obj.eigThresh_im)./(1-obj.eigThresh_im).* (W(:,:,end-(obj.n_maps-1):end) > obj.eigThresh_im);
            weights = -cos(pi*weights)/2 + 1/2;

            % create an ESPIRiT operator
            ESP = ESPIRiT(maps,weights);

            if(strcmp(obj.solver,'L1'))
                FT = p2DFT(squeeze(obj.fullMask(:,:,iFreq,:)),[nPha,nZ,nCha]);
                if(strcmp(obj.trafo.trafoType,'fft'))
                    XOP = FT;
                end
                image(:,iFreq,:,:) = permute(obj.cgL1ESPIRiT(squeeze(kSpaceL(:,:,iFreq,:)), zeros([size(squeeze(kSpaceL(:,:,iFreq,1))),obj.n_maps]), FT, ESP, obj.iNINNER, XOP, obj.trafo.wavWeight, obj.splitWeight, obj.iNIterSplit),[1 4 2 3]); % => y-x-z-cha
            else
                [tmpK,tmp] = obj.cgESPIRiT(squeeze(kSpaceL(:,:,iFreq,:)), ESP, obj.iNINNER, obj.calibTyk, zeros(size(squeeze(kSpaceL(:,:,iFreq,:)))));
                image(:,iFreq,:,:) = permute(tmp,[1 4 2 3]);
            end
            dispProgress('Frequency', iFreq/nFreq);
        end
        dispProgress('Frequency','Close');
    %     image = ifftnshift(image,[1,3]);
        imageCha(:,1:obj.n_maps) = shiftdim(mat2cell(image,nPha,nFreq,nZ,ones(1,obj.n_maps)),2);
        [imageCha{:,obj.n_maps+1:end}] = deal(zeros(size(imageCha{1,1})));

    else
        fprintf('ESPIRiT reconstruction\n');
        if(ispc)
            eval(sprintf('ESP = bart(''ecalib -r %d:%d -k %d:%d -m %d'',kSpaceL);', obj.calibSize(1), obj.calibSize(2), obj.kernelSize(1), obj.kernelSize(2), obj.n_maps ));
            if(strcmp(obj.solver,'L1'))
                eval(sprintf('image = permute(squeeze(bart(''pics -l1 -r %f'', kSpaceL, ESP)),[1 3 2 4]);', obj.trafo.wavWeight));
            else
                image = permute(squeeze(bart('pics -l2', kSpaceL, ESP)),[1 3 2 4]);
            end
        else
            % create temporary dir if not yet existing
            if(~exist([obj.path,filesep,'tmpESPIRiT'],'dir'))
                mkdir([obj.path,filesep,'tmpESPIRiT']);
            end
            fftdims = size(kSpaceL);
            fftdims = fftdims(1:3);
            if(any(mod(fftdims,2) ~= 0))
                writecfl([obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], zpad(kSpaceL,size(kSpaceL) + [mod(fftdims,2),0])); % k_y - k_z - k_x - cha
            else
                writecfl([obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], kSpaceL); % k_y - k_z - k_x - cha
            end
            obj.calibSize = obj.calibSize - mod(obj.calibSize,2);
            eval(sprintf('! ecalib -r %d:%d:%d -k %d:%d:%d -m %d %s %s', obj.calibSize(1), obj.calibSize(2), obj.calibSize(3), obj.kernelSize(1), obj.kernelSize(2), obj.kernelSize(3), obj.n_maps, [obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], [obj.path,filesep,'tmpESPIRiT',filesep,'espiritmaps']));
            ESP = [obj.path,filesep,'tmpESPIRiT',filesep,'espiritmaps'];

            if(strcmp(obj.solver,'L1'))
                eval(sprintf('! sense -l1 -r %f %s %s %s', obj.trafo.wavWeight, [obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], ESP, [obj.path,filesep,'tmpESPIRiT',filesep,'imgOut']));
            else
                eval(sprintf('! sense -l2 %s %s %s', [obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], ESP, [obj.path,filesep,'tmpESPIRiT',filesep,'imgOut']));
            end
            image = permute(squeeze(readcfl([obj.path,filesep,'tmpESPIRiT',filesep,'imgOut'])),[1 3 2 4]);
        end
%         % OR:
%         kSpaceL = permute(ifftnshift(kSpaceL,3), [3 1 2 4]); % x - k_y - k_z - cha
%         writecfl([obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], kSpaceL); % x - k_y - k_z - cha
%         eval(sprintf('! ecalib -r %d:%d:%d -k %d:%d:%d -m %d %s %s', obj.calibSize(3), obj.calibSize(1), obj.calibSize(2), obj.kernelSize(3), obj.kernelSize(1), obj.kernelSize(2), obj.n_maps, [obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], [obj.path,filesep,'tmpESPIRiT',filesep,'espiritmaps']));
%         ESP = [obj.path,filesep,'tmpESPIRiT',filesep,'espiritmaps'];
%          if(strcmp(obj.solver,'L1'))
%             eval(sprintf('! rsense -l1 -r %f %s %s %s', obj.trafo.wavWeight, [obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], ESP, [obj.path,filesep,'tmpESPIRiT',filesep,'imgOut']));
%         else
%             eval(sprintf('! rsense -l2 %s %s %s', [obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], ESP, [obj.path,filesep,'tmpESPIRiT',filesep,'imgOut']));
%         end
%         image = permute(squeeze(readcfl([obj.path,filesep,'tmpESPIRiT',filesep,'imgOut'])), [2 1 3 4]);

        image = crop(image,[size(kSpaceL,1),size(kSpaceL,3),size(kSpaceL,2),obj.n_maps]);
        imageCha(:,1:obj.n_maps) = shiftdim(mat2cell(image,nPha,nFreq,nZ,ones(1,obj.n_maps)),2);
        [imageCha{:,obj.n_maps+1:end}] = deal(zeros(size(imageCha{1,1})));
    end

elseif(strcmp(obj.measPara.dimension,'4D'))
    % TODO

end

if(obj.flagToolbox)
    rmdir([obj.path,filesep,'tmpESPIRiT'], 's');
end

end

author: thomaskuestner


text:
function ThrowFcn(~, ~, game, FinishBtn, SelectBtn)

%% GET OFFENDER/DEFENDER DATA
offender = get(game.OffenseTxt, 'UserData');
defender = get(game.DefenseTxt, 'UserData');

offOcc = get(game.board(offender).Patch, 'UserData');
defOcc = get(game.board(defender).Patch, 'UserData');   % Defending state may already have been occupied in this turn

%% CHECK WHETHER DICE WERE THROWN AT THE APPROPRIATE MOMENT
if offOcc(1) == defOcc(1)   % same occupant
    msg = msgbox('You already invaded this state');
    uiwait(msg)
    return
end

%% NUMBER OF TROOPS
Noff = get(game.board(offender).Text, 'UserData');

%% GET NUMBER OF TROOPS THE OFFENDER WANTS TO SEND
while true
    N = str2double(inputdlg('Number of offending troops', 'How many troops?', 1, {'1'}));

    if isempty(N) || isnan(N)
        return
    end
    if N > Noff - 1
        uiwait(msgbox(sprintf('Maximum number of offending troops is %d', Noff - 1)))
        continue
    end
    if N <= 0
        uiwait(msgbox('Send at least 1 unit into battle'))
        continue
    end

    break
end

%% NO FINISH / RESELECT DURING BATTLE
set(FinishBtn, 'Enable', 'inactive')
set(SelectBtn, 'Enable', 'inactive')

%% NUMBER OF DICE TO BE USED
M(1) = N * (N <= 3) + 3 * (N > 3);

%% RESULT
res1 = round(rand(1,M(1)) * 5) + 1;
for i = 1:M(1)
    imshow(sprintf('dice/%d.jpg', res1(i)), 'Parent', game.dice(i));
end
for i = M(1) + 1 : 3
    imshow('dice/7.jpg', 'Parent', game.dice(i));
end

%% UNSET OFFENSE THROWBUTTON, SET DEFENSE
set(game.ThrowBtnOff, 'Enable', 'off')
set(game.ThrowBtnDef, 'Enable', 'on', 'Callback', {@ThrowFcn2, game, res1, N, M, FinishBtn, SelectBtn});
set(game.txtH, 'String', 'Defender may throw the dice')


text:
function msfun_khepera3(block)

% Copyright (C) 2013, Georgia Tech Research Corporation
% see the LICENSE file included with this software

%%
%% The setup method is used to set up the basic attributes of the
%% S-function such as ports, parameters, etc. Do not add any other
%% calls to the main body of the function.
%%
robot = [];

setup(block);

%% Function: setup ===================================================
%% Abstract:
%%   Set up the basic characteristics of the S-function block such as:
%%   - Input ports
%%   - Output ports
%%   - Dialog parameters
%%   - Options
%%
%%   Required         : Yes
%%   C-Mex counterpart: mdlInitializeSizes
%%
    function setup(block)

        % Register number of ports
        block.NumInputPorts  = 2;
        block.NumOutputPorts = 2;

        % Setup port properties to be inherited or dynamic
        % block.SetPreCompInpPortInfoToDynamic;
        % block.SetPreCompOutPortInfoToDynamic;
        block.SetPreCompPortInfoToDefaults;

        % Override input port properties
        block.InputPort(1).Dimensions        = 1;
        block.InputPort(1).DatatypeID  = 0;  % double
        block.InputPort(1).Complexity  = 'Real';
        block.InputPort(1).DirectFeedthrough = false;

        block.InputPort(2).Dimensions        = 1;
        block.InputPort(2).DatatypeID  = 0;  % double
        block.InputPort(2).Complexity  = 'Real';
        block.InputPort(2).DirectFeedthrough = false;

        block.InputPort(1).SamplingMode = 0;
        block.InputPort(2).SamplingMode = 0;

        % Override output port properties
        block.OutputPort(1).Dimensions       = 9;
        block.OutputPort(1).DatatypeID  = 0; % double
        block.OutputPort(1).Complexity  = 'Real';

        block.OutputPort(2).Dimensions       = 2;
        block.OutputPort(2).DatatypeID  = 0; % double
        block.OutputPort(2).Complexity  = 'Real';

        % Set up the states
        block.NumContStates = 0;
        block.NumDworks = 0;

        % Register parameters
        block.NumDialogPrms     = 3;
        block.DialogPrmsTunable = {'NonTunable', 'NonTunable', 'NonTunable'};

        % Register sample times
        %  [0 offset]            : Continuous sample time
        %  [positive_num offset] : Discrete sample time
        %
        %  [-1, 0]               : Inherited sample time
        %  [-2, 0]               : Variable sample time
        block.SampleTimes = [0.033 0];

        % Specify the block simStateCompliance. The allowed values are:
        %    'UnknownSimState', < The default setting; warn and assume DefaultSimState
        %    'DefaultSimState', < Same sim state as a built-in block
        %    'HasNoSimState',   < No sim state
        %    'CustomSimState',  < Has GetSimState and SetSimState methods
        %    'DisallowSimState' < Error out when saving or restoring the model sim state
        block.SetAccelRunOnTLC(false); % run block in interpreted mode even w/ Acceleration
        block.SimStateCompliance = 'DefaultSimState';

        %% -----------------------------------------------------------------
        %% The MATLAB S-function uses an internal registry for all
        %% block methods. You should register all relevant methods
        %% (optional and required) as illustrated below. You may choose
        %% any suitable name for the methods and implement these methods
        %% as local functions within the same file. See comments
        %% provided for each function for more information.
        %% -----------------------------------------------------------------

        block.RegBlockMethod('PostPropagationSetup',    @DoPostPropSetup);
        block.RegBlockMethod('InitializeConditions', @InitializeConditions);
        block.RegBlockMethod('Start', @Start);
        block.RegBlockMethod('Outputs', @Outputs);     % Required
        block.RegBlockMethod('Update', @Update);
        block.RegBlockMethod('Terminate', @Terminate); % Required

    end %end setup

%%
%% PostPropagationSetup:
%%   Functionality    : Setup work areas and state variables. Can
%%                      also register run-time methods here
%%   Required         : No
%%   C-Mex counterpart: mdlSetWorkWidths
%%
    function DoPostPropSetup(block)

    end

%%
%% InitializeConditions:
%%   Functionality    : Called at the start of simulation and if it is
%%                      present in an enabled subsystem configured to reset
%%                      states, it will be called when the enabled subsystem
%%                      restarts execution to reset the states.
%%   Required         : No
%%   C-MEX counterpart: mdlInitializeConditions
%%
    function InitializeConditions(block)

    end %end InitializeConditions

%%
%% Start:
%%   Functionality    : Called once at start of model execution. If you
%%                      have states that should be initialized once, this
%%                      is the place to do it.
%%   Required         : No
%%   C-MEX counterpart: mdlStart
%%
    function Start(block)

        x_0 = block.DialogPrm(1).Data;
        y_0 = block.DialogPrm(2).Data;
        theta_0 = block.DialogPrm(3).Data;

        customData = getSimulatorBlockUserData(bdroot(block.BlockHandle));
        simulator = customData('simulatorHandle');

        robot = simulator.world.add_robot('Khepera3', 'Supervisor', x_0, y_0, theta_0);

        % store info in custom data;
        customData = containers.Map('UniformValues', false);
        customData('robotHandle') = robot;
        set(block.BlockHandle, 'UserData', customData, 'UserDataPersistent', 'off');

    end %endfunction

%%
%% Outputs:
%%   Functionality    : Called to generate block outputs in
%%                      simulation step
%%   Required         : Yes
%%   C-MEX counterpart: mdlOutputs
%%
    function Outputs(block)

%         customData = get(block.BlockHandle, 'UserData');
%         robot = customData('robotHandle');

        block.OutputPort(1).Data = robot.ir_array.get_range();
        robot.encoders.ticks
        block.OutputPort(2).Data = [robot.encoders.ticks];

    end %end Outputs

%%
%% Update:
%%   Functionality    : Called to update discrete states
%%                      during simulation step
%%   Required         : No
%%   C-MEX counterpart: mdlUpdate
%%
    function Update(block)

%         customData = get(block.BlockHandle, 'UserData');
%         robot = customData('robotHandle');

        vel_r = block.InputPort(1).Data;
        vel_l = block.InputPort(2).Data;

        robot.set_wheel_speeds(vel_r,vel_l);

    end %end Update

%%
%% Terminate:
%%   Functionality    : Called at the end of simulation for cleanup
%%   Required         : Yes
%%   C-MEX counterpart: mdlTerminate
%%
    function Terminate(block)

    end %end Terminate

end

author: jdelacroix

repo: simiam

save_path: github-repos/MATLAB/jdelacroix-simiam


text:
function comp_next_tests ( )

%****************************************************************************80
%
%% COMP_NEXT_TESTS calls COMP_NEXT_TEST with various arguments.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 March 2011
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../sandia_rules' );

  test_num = 12;

  dim_num_array(1:test_num) = [ ...
    2, 2, 2, 2, 2, ...
    3, 3, 3, 3, 3, ...
    4, 4 ];
  level_max_array(1:test_num) = [ ...
    0, 1, 2, 3, 4, ...
    0, 1, 2, 3, 4, ...
    2, 3 ];

  fprintf ( 1, '\n' );
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'COMP_NEXT_TESTS\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Call COMP_NEXT_TEST with various arguments.\n' );

  for test = 1 : test_num

    dim_num = dim_num_array(test);
    level_max = level_max_array(test);
    comp_next_test ( dim_num, level_max );

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'COMP__NEXT_TESTS\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../sandia_rules' );

  return
end


text:
function test_ft_freqanalysis

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_freqanalysis

ntrl = 1000;
nchan = 32;
fs = 500;
start_time = -1; % seconds
end_time = 2.5; % seconds
nsamples = (end_time - start_time) * fs + 1;

data = [];
data.label = cellstr(num2str((1:nchan).'));
for i=1:ntrl
  data.time{i} = linspace(start_time, end_time, nsamples);
  data.trial{i} = randn(nchan,nsamples);
end

cfg = [];
cfg.method = 'mtmfft';
cfg.taper = 'hanning';
cfg.foilim = [4 30];
freq = ft_freqanalysis(cfg, data);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function rinexe(ephemerisfile, outputfile)
%RINEXE Reads a RINEX Navigation Message file version 3.03 and
%	        reformats the data into a matrix with 21 rows and a column
%           for each satellite.  The matrix is stored in outputfile

%Typical call: rinexe('pta.96n','pta.nav')

%Kai Borre 04-18-96
%Copyright (c) by Kai Borre
%$Revision: 1.0 $  $Date: 1997/09/24  $
% Revision September 11, 2015 at Samara

% Units are either seconds, meters, or radians
fide = fopen(ephemerisfile);
head_lines = 0;
while 1  % We skip header
   head_lines = head_lines+1;
   line = fgetl(fide);
   answer = findstr(line,'END OF HEADER');
   if ~isempty(answer), break;	end;
end;
head_lines;
noeph = -1;
while 1
   noeph = noeph+1;
   line = fgetl(fide);
   if line == -1, break;  end
end;
noeph = noeph/8
frewind(fide);
for i = 1:head_lines, line = fgetl(fide); end;

% Set aside memory for the input
svprn	 = zeros(1,noeph);
weekno	 = zeros(1,noeph);
t0c	 = zeros(1,noeph);
tgd	 = zeros(1,noeph);
aodc	 = zeros(1,noeph);
toe	 = zeros(1,noeph);
af2	 = zeros(1,noeph);
af1	 = zeros(1,noeph);
af0	 = zeros(1,noeph);
aode	 = zeros(1,noeph);
deltan	 = zeros(1,noeph);
M0	 = zeros(1,noeph);
ecc	 = zeros(1,noeph);
roota	 = zeros(1,noeph);
toe	 = zeros(1,noeph);
cic	 = zeros(1,noeph);
crc	 = zeros(1,noeph);
cis	 = zeros(1,noeph);
crs	 = zeros(1,noeph);
cuc	 = zeros(1,noeph);
cus	 = zeros(1,noeph);
Omega0	 = zeros(1,noeph);
omega	 = zeros(1,noeph);
i0	 = zeros(1,noeph);
Omegadot = zeros(1,noeph);
idot	 = zeros(1,noeph);
accuracy = zeros(1,noeph);
health	 = zeros(1,noeph);
fit	 = zeros(1,noeph);

for i = 1:noeph
   line = fgetl(fide);	  %
   svprn(i) = str2num(line(2:3)); % a G for GPS in front is omitted
   year = line(5:8);
   month = line(10:11);
   day = line(13:14);
   hour = line(16:17);
   minute = line(19:20);
   second = line(22:23);
   af0(i) = str2num(line(24:42));
   af1(i) = str2num(line(43:61));
   af2(i) = str2num(line(62:80));
   line = fgetl(fide);	  %
   IODE = line(5:23);
   crs(i) = str2num(line(24:42));
   deltan(i) = str2num(line(43:61));
   M0(i) = str2num(line(62:80));
   line = fgetl(fide) ;	  %
   cuc(i) = str2num(line(5:23));
   ecc(i) = str2num(line(24:42));
   cus(i) = str2num(line(43:61));
   roota(i) = str2num(line(62:80));
   line=fgetl(fide);
   toe(i) = str2num(line(5:23));
   cic(i) = str2num(line(24:42));
   Omega0(i) = str2num(line(43:61));
   cis(i) = str2num(line(62:80));
   line = fgetl(fide);	    %
   i0(i) =  str2num(line(5:23));
   crc(i) = str2num(line(24:42));
   omega(i) = str2num(line(43:61));
   Omegadot(i) = str2num(line(62:80));
   line = fgetl(fide);	    %
   idot(i) = str2num(line(5:23));
   codes = str2num(line(24:42));
   weekno = str2num(line(43:61));
   L2flag = str2num(line(62:80));
   line = fgetl(fide);	    %
   svaccur = str2num(line(5:23));
   svhealth = str2num(line(24:42));
   tgd(i) = str2num(line(43:61));
   iodc = line(62:80);
   line = fgetl(fide);	    %
   tom(i) = str2num(line(5:23));
   spare = line(24:42);
%   spare = line(43:61);
%   spare = line(62:80);
end
status = fclose(fide);

%  Description of variable eph.
eph(1,:)  = svprn;
eph(2,:)  = af2;
eph(3,:)  = M0;
eph(4,:)  = roota;
eph(5,:)  = deltan;
eph(6,:)  = ecc;
eph(7,:)  = omega;
eph(8,:)  = cuc;
eph(9,:)  = cus;
eph(10,:) = crc;
eph(11,:) = crs;
eph(12,:) = i0;
eph(13,:) = idot;
eph(14,:) = cic;
eph(15,:) = cis;
eph(16,:) = Omega0;
eph(17,:) = Omegadot;
eph(18,:) = toe;
eph(19,:) = af0;
eph(20,:) = af1;
eph(21,:) = toe;

fidu = fopen(outputfile,'w');
count = fwrite(fidu,[eph],'double');
fclose all;
%%%%%%%%% end rinexe.m %%%%%%%%%

author: yandld


text:
function [f]=comp_idwilt(coef,g)
%COMP_IDWILT  Compute Inverse discrete Wilson transform.
%
%   This is a computational routine. Do not call it
%   directly.

%   AUTHOR : Peter L. Søndergaard.
%   TESTING: OK
%   REFERENCE: OK

M=size(coef,1)/2;
N=2*size(coef,2);
W=size(coef,3);

a=M;

L=N*a;

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

% First and middle modulation are transferred unchanged.
coef2(1,1:2:N,:) = coef(1,:,:);
if mod(M,2)==0
  coef2(M+1,1:2:N,:) = coef(M+1,:,:);
else
  coef2(M+1,2:2:N,:) = coef(M+1,:,:);
end;

if M>2
  % cosine, first column.
  coef2(3:2:M,1:2:N,:)        = 1/sqrt(2)*coef(3:2:M,:,:);
  coef2(2*M-1:-2:M+2,1:2:N,:) = 1/sqrt(2)*coef(3:2:M,:,:);

  % sine, second column
  coef2(3:2:M,2:2:N,:)        = -1/sqrt(2)*i*coef(M+3:2:2*M,:,:);
  coef2(2*M-1:-2:M+2,2:2:N,:) =  1/sqrt(2)*i*coef(M+3:2:2*M,:,:);
end;

% sine, first column.
coef2(2:2:M,1:2:N,:)        = -1/sqrt(2)*i*coef(2:2:M,:,:);
coef2(2*M:-2:M+2,1:2:N,:)   =  1/sqrt(2)*i*coef(2:2:M,:,:);

% cosine, second column
coef2(2:2:M,2:2:N,:)        = 1/sqrt(2)*coef(M+2:2:2*M,:,:);
coef2(2*M:-2:M+2,2:2:N,:)   = 1/sqrt(2)*coef(M+2:2:2*M,:,:);

f = comp_isepdgt(coef2,g,L,a,2*M,0);

% Apply the final DGT
%f=comp_idgt(coef2,g,a,[0 1],0,0);

% Clean signal if it is known to be real
if (isreal(coef) && isreal(g))
  f=real(f);
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [pnt, tri, srf] = read_bv_srf(filename)

% READ_BV_SRF reads a triangulated surface from a BrainVoyager *.srf file
%
% Use as
%   [pnt, tri] = read_bv_srf(filename) or
%   [pnt, tri, srf] = read_bv_srf(filename)

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

% This documentation originates from
%   http://www.brainvoyager.com/BV2000OnlineHelp/BrainVoyagerWebHelp/mergedProjects/FileFormats/BrainVoyager_File_Formats.htm
%
% BYTES DATA TYPE DEFAULT DESCRIPTION
% 4 float 3 version number
% 4 int 0 reserved, must be '0'
% 4 int  NrOfVertices (number of vertices)
% 4 int  NrOfTriangles (number of triangles)
% 4 float 128.0 MeshCenterX
% 4 float 128.0 MeshCenterY
% 4 float 128.0 MeshCenterZ
% NrOfVertices*4 float  VertexX, sequence of X coordinates of all vertices
% NrOfVertices*4 float  VertexY, sequence of Y coordinates of all vertices
% NrOfVertices*4 float  VertexZ, sequence of Z coordinates of all vertices
% NrOfVertices*4 float  NormalX, sequence of X components of all vertex normals
% NrOfVertices*4 float  NormalY, sequence of Y components of all vertex normals
% NrOfVertices*4 float  NormalZ, sequence of Z components of all vertex normals
% 4 float 0.322 R component of convex curvature color (range: 0.0 - 1-0)
% 4 float 0.733 G component of convex curvature color (range: 0.0 - 1-0)
% 4 float 0.980 B component of convex curvature color (range: 0.0 - 1-0)
% 4 float 0.500 Alpha component of convex curvature color (range: 0.0 - 1-0)
% 4 float 0.100 R component of concave curvature color (range: 0.0 - 1-0)
% 4 float 0.240 G component of concave curvature color (range: 0.0 - 1-0)
% 4 float 0.320 B component of concave curvature color (range: 0.0 - 1-0)
% 4 float 0.500 Alpha component of concave curvature color (range: 0.0 - 1-0)
% NrOfVertices*4 int  MeshColor, sequence of color indices of all vertices (*1)
% 4 int  N, number of nearest neighbors of vertex 1
% 4 int  Nearest neighbor 1 of vertex 1
% : :  :
% 4 int  Nearest neighbor N of vertex 1
% :: ::  ::
% 4 int  N, number of nearest neighbors of vertex 'NrOfVertices'
% 4 int  Nearest neighbor 1 of vertex 'NrOfVertices'
% : :  :
% 4 int  Nearest neighbor N of vertex 'NrOfVertices'
% NrOfTriangels*3*4 int  Sequence of three indices to constituting vertices of each triangle
% 4 int  NrOfTriangleStripElements
% NrOfStripElements*4 int  Sequence of strip elements (if NrOfStripElements > 0)

fid = fopen_or_error(filename, 'rb', 'ieee-le');

srf.version_number                                        = fread(fid, 1, 'float');
srf.reserved                                              = fread(fid, 1, 'int'  );
srf.NrOfVertices                                          = fread(fid, 1, 'int'  );
srf.NrOfTriangles                                         = fread(fid, 1, 'int'  );
NrOfVertices  = srf.NrOfVertices;
NrOfTriangles = srf.NrOfTriangles;
srf.MeshCenterX                                           = fread(fid, 1, 'float');
srf.MeshCenterY                                           = fread(fid, 1, 'float');
srf.MeshCenterZ                                           = fread(fid, 1, 'float');
srf.X_coordinates                                         = fread(fid, NrOfVertices, 'float');
srf.Y_coordinates                                         = fread(fid, NrOfVertices, 'float');
srf.Z_coordinates                                         = fread(fid, NrOfVertices, 'float');
srf.X_components                                          = fread(fid, NrOfVertices, 'float');
srf.Y_components                                          = fread(fid, NrOfVertices, 'float');
srf.Z_components                                          = fread(fid, NrOfVertices, 'float');
srf.R_component_of_convex_curvature_color                 = fread(fid, 1, 'float');
srf.G_component_of_convex_curvature_color                 = fread(fid, 1, 'float');
srf.B_component_of_convex_curvature_color                 = fread(fid, 1, 'float');
srf.Alpha_component_of_convex_curvature_color             = fread(fid, 1, 'float');
srf.R_component_of_concave_curvature_color                = fread(fid, 1, 'float');
srf.G_component_of_concave_curvature_color                = fread(fid, 1, 'float');
srf.B_component_of_concave_curvature_color                = fread(fid, 1, 'float');
srf.Alpha_component_of_concave_curvature_color            = fread(fid, 1, 'float');
srf.MeshColor                                             = fread(fid, NrOfVertices, 'int'  );
for i=1:NrOfVertices
  number           = fread(fid, 1, 'int'  );
  srf.neighbour{i} = fread(fid, number, 'int'  );
end
srf.Triangles                                             = fread(fid, [3 NrOfTriangles], 'int'  );
srf.NrOfTriangleStripElements                             = fread(fid, 1, 'int'  );
srf.sequence_of_strip_elements                            = fread(fid, srf.NrOfTriangleStripElements, 'int'  );

fclose(fid);

pnt = [srf.X_coordinates(:) srf.Y_coordinates(:) srf.Z_coordinates(:)];
tri = srf.Triangles' + 1;


text:
function confidences = test_boosted_dt_mc(classifier, features)
% confidences = test_boosted_dt_mc(classifier, features)
% returns a log likelihod ratio for each class in the classifier
% confidences(ndata, nclasses)

dt = classifier.wcs(1).dt;
if size(features, 2)~=dt.npred
    error('Incorrect number of attributes')
end

wcs = classifier.wcs;
nclasses = size(wcs, 2);

ntrees = size(wcs, 1);

confidences = zeros(size(features, 1), nclasses);
for c = 1:nclasses
    for t = 1:ntrees
        if ~isempty(wcs(t,c).dt)
            dt = wcs(t,c).dt;
            [var, cut, children, catsplit] = tree_getParameters(dt);
            nodes = treevalc(int32(var), cut, int32(children(:, 1)), ...
                    int32(children(:, 2)), catsplit(:, 1), features');
            %[class_indices, nodes2, classes] = treeval(wcs(t, c).dt, features);
            %         if sum(nodes~=nodes2)>0
            %              disp('error')
            % %              disp(num2str([nodes nodes2]))
            %         end
            confidences(:, c) = confidences(:, c) + wcs(t, c).confidences(nodes);
        end
    end
    confidences(:, c) = confidences(:, c) + classifier.h0(c);
end

author: Cloud-CV

repo: object-proposals


text:
classdef MeshCreatorFromInputFile < handle

    properties (Access = public)
        backgroundMesh
        boundaryMesh
    end

    properties (Access = private)
        inputFile
    end

    methods (Access = public)

        function obj = MeshCreatorFromInputFile(cParams)
            obj.init(cParams);
            obj.createBackgroundMesh();
            obj.createBoundaryMesh();
        end

    end

    methods (Access = private)

        function init(obj,cParams)
           obj.inputFile = cParams.inputFile;
        end

        function createBackgroundMesh(obj)
            coord = [];
            connec = [];
            eval(obj.inputFile);
            s.coord  = coord(:,2:3);
            s.connec = connec(:,2:end);
            obj.backgroundMesh = Mesh(s);
        end

        function createBoundaryMesh(obj)
            eval(obj.inputFile);
            if exist('External_border_nodes','var')
                s.borderNodes    = External_border_nodes;
                s.borderElements = External_border_elements;
                s.backgroundMesh = obj.backgroundMesh;
                b = BoundaryMeshCreatorFromData(s);
                obj.boundaryMesh = b.create();
            else
            s.backgroundMesh = obj.backgroundMesh;
            s.dimension = 1:obj.backgroundMesh.ndim;
            bC = BoundaryMeshCreatorFromRectangularBox(s);
            obj.boundaryMesh = bC.create();
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [x, timing] = blendenpik(A, b, params)
% [x, timing] = blendenpik(A, b, params)
%
% Solve the equation x = arg min norm(A * x - b, 2)
% using Blendenpik.
%
% "params" - parameters governning the method.
%    params.type - type of mixing transform. Optional values: 'DCT', 'DHT', 'WHT'.
%                  Default is DHT.
%    params.gamma - gamma * min(n,m) rows/columns will be sampled (A is m-by-n).
%                   Default is 4.
%    params.preprocess_steps - number of mixing steps to do in advance.
%                               Default is 1.
%    params.maxcond - maximum condition number of the preconditioner.
%                     Default is 1 / (5 * epsilon_machine).
%    params.tol - convergence thershold for LSQR.
%                 Default is 1e-14.
%    params.maxit - maximum number of LSQR iterations.
%                   Default is 1000.
%    params.lsvec - whether to output in "timing" the LSQR residuals.
%                   Default is false.
%    params.use_full_lsqr - whether to use LSQR with full
%                           orthogonalization. Useful for
%                           preprocess_steps=0.
%                           Default is false.
%
% Output:
%   x - the solution.
%   timing - statistics on the time spent on various phases.
%
% 6-December 2009, Version 1.3
% Copyright (C) 2009, Haim Avron and Sivan Toledo.

if (nargin < 3)
    params = struct;
end

[m, n] = size(A);
if (m >= n)
    x = blendenpik_over(A, b, params);
else
    x = blendenpik_under(A, b, params);
end


text:
% Run this file to build the needed mex-files on windows

% Build merResize
mex -lopencv_core242 -lopencv_imgproc242 -L./ -I./ mexResize.cpp MxArray.cpp

% Build setnonzeros from the lightspeed matlab toolbox
mex setnonzeros.c

% Build gradientMex from Piotrs toolbox
mex gradientMex.cpp -I./

% Build mtimesx
mtimesx_build

author: flyers

repo: drone-tracking


text:
function feature3M = calcFeatureImpact(scanNum, structNum, quantizeS, ...
    patchSizeV, featureFun, featureName, planC, varargin)
% function feature3M = calcFeatureImpact(scanNum, structNum, quantizeS,...
% patchSizeV, featureFun, featureName, planC, varargin)
%
% This function computes the feature Impact at each voxel.
%
% Method:
% 1> Omit the voxel and its neighborhood.
% 2> Compute the feature with the passed featureFun and the varargin.
% Repeat 1 and 2 for each voxel in the region of interest
% 3> Calc. percentage difference between features calculated by omitting
% vs. using all voxels
%
% Example:
%
% In order to compute the impact map for 1st order kurtosis,
% scanNum = 1;
% structNum = 6;
% patchRadiusV = [2 2 0];
% funcHandle = @radiomics_first_order_stats;
% featureName = 'kurtosis';
% quantizeS.quantize = 1;
% quantizeS.nL = 32;
% quantizeS.xmin = [];
% quantizeS.xmax = [];
% quantizeS.binwidth = [];
% feature3M = calcFeatureImpact(scanNum, structNum, quantizeS, patchRadiusV,...
%   funcHandle, featureName, planC, structNum);
%
% Also, refer to call_calcFeatureImpact.m for an example to calculate
% feature impact for RLM features
%
% APA, 5/18/2017
% AI,  8/21/2018 Updates to inputs, featureFun signature

indexS = planC{end};

% Get uniformized structure Mask
struct3M = getUniformStr(structNum,planC);

% Get uniformized scan mask in HU
scan3M = getUniformizedCTScan(1, scanNum, planC);
% Convert to HU if image is of type CT
if ~isempty(planC{indexS.scan}(scanNum).scanInfo(1).CTOffset)
    scan3M = double(scan3M) - planC{indexS.scan}(scanNum).scanInfo(1).CTOffset;
end

fullScanSiz = size(scan3M);

% Crop scan within the structures bounding box
[minr, maxr, minc, maxc, mins, maxs] = compute_boundingbox(struct3M);
struct3M = struct3M(minr:maxr,minc:maxc,mins:maxs);
scan3M = scan3M(minr:maxr,minc:maxc,mins:maxs);
scan3M(struct3M==0)     = NaN;

% Initialize the feature3M matrix
feature3M = zeros(fullScanSiz);

featureCropped3M = zeros(size(scan3M));

% Get indices of non-NaN voxels
calcIndM = struct3M == 1;

% % Grid resolution
slcWindow = 2 * patchSizeV(3) + 1;
rowWindow = 2 * patchSizeV(1) + 1;
colWindow = 2 * patchSizeV(2) + 1;

% Build distance matrices
numColsPad = floor(colWindow/2);
numRowsPad = floor(rowWindow/2);
numSlcsPad = floor(slcWindow/2);

% Get number of voxels per slice
[numRows, numCols, numSlices] = size(scan3M);
numVoxels = numRows*numCols;

% Pad q, so that sliding window works also for the edge voxels
if exist('padarray.m','file')
    %scanArrayTmp3M = padarray(scanArray3M,[numRowsPad numColsPad
    %numSlcsPad],NaN,'both'); % aa commented
    q = padarray(scan3M,[numRowsPad numColsPad numSlcsPad],NaN,'both');
    structPadded3M = padarray(uint32(struct3M),[numRowsPad numColsPad numSlcsPad],NaN,'both');
    calcIndM = padarray(calcIndM,[0 0 numSlcsPad],0,'both');
else
    %scanArrayTmp3M = padarray(scanArray3M,[numRowsPad numColsPad
    %numSlcsPad],NaN,'both'); % aa commented
    q = padarray_oct(scan3M,[numRowsPad numColsPad numSlcsPad],NaN,'both');
    structPadded3M = padarray_oct(uint32(struct3M),[numRowsPad numColsPad numSlcsPad],NaN,'both');
    calcIndM = padarray_oct(calcIndM,[0 0 numSlcsPad],0,'both');
end
structPadded3M = logical(structPadded3M);

% Create indices for 2D blocks
[m,n,~] = size(q);
m = uint32(m);
n = uint32(n);
colWindow = uint32(colWindow);
rowWindow = uint32(rowWindow);
slcWindow = uint32(slcWindow);

% Index calculation adapted from
% http://stackoverflow.com/questions/25449279/efficient-implementation-of-im2col-and-col2im

%// Start indices for each block
start_ind = reshape(bsxfun(@plus,[1:m-rowWindow+1]',[0:n-colWindow]*m),[],1); %//'

%// Row indices
lin_row = permute(bsxfun(@plus,start_ind,[0:rowWindow-1])',[1 3 2]);  %//'

%// Get linear indices based on row and col indices and get desired output
indM = reshape(bsxfun(@plus,lin_row,(0:colWindow-1)*m),rowWindow*colWindow,[]);

% Feature computed using all voxels
quantized3M = q;
if ~isempty(quantizeS) && quantizeS.quantize
    quantized3M = imquantize_cerr(quantized3M,quantizeS.nL,quantizeS.xmin,quantizeS.xmax,quantizeS.binwidth);
end
featureEntireStruct = feval(featureFun,quantized3M,varargin{:});

if ~isempty(featureName) && isstruct(featureEntireStruct)
    [f1,f2] = strtok(featureName,'.');
    fieldC = {f1};
    while ~isempty(f2)
        [f1,f2] = strtok(f2,'.');
        fieldC{end+1} = f1;
    end
    for i = 1:length(fieldC)
        featureEntireStruct = featureEntireStruct.(fieldC{i});
    end
end

% Iterate over slices. Compute cooccurance for all patches per slice
for slcNum = (1+numSlcsPad):(numSlices+numSlcsPad)
    disp(['--- Feature Impact Calculation for Slice # ', num2str(slcNum), ' ----'])

    calcSlcIndV = calcIndM(:,:,slcNum);
    calcSlcIndV = calcSlcIndV(:);
    numCalcVoxs = sum(calcSlcIndV);
    indSlcM = indM(:,calcSlcIndV);
    slcV = slcNum-patchSizeV(3):slcNum+patchSizeV(3);
    featureV = zeros(numCalcVoxs,1);

    for vox = 1:size(indSlcM,2)
        scanTmp3M = quantized3M;
        for slcNumTmp = 1:length(slcV)
            slc = slcV(slcNumTmp);
            scanSlcM = scanTmp3M(:,:,slc);
            scanSlcM(indSlcM(:,vox)) = NaN;
            scanTmp3M(:,:,slc) = scanSlcM;
        end
        featureVox = feval(featureFun,scanTmp3M,varargin{:});
        for i = 1:length(fieldC)
            featureVox = featureVox.(fieldC{i});
        end
        featureV(vox) = featureVox;
    end

    indScanTmp3M = zeros(size(calcIndM),'logical');
    indScanTmp3M(:,:,slcNum) = calcIndM(:,:,slcNum);
    featuresC{slcNum} = featureV;
    calcIndC{slcNum} = indScanTmp3M;

end

for slcNum = (1+numSlcsPad):(numSlices+numSlcsPad)
    featureCropped3M(calcIndC{slcNum}) = featuresC{slcNum};
end

featureCropped3M(calcIndM) = ...
    (featureCropped3M(calcIndM) - featureEntireStruct) / featureEntireStruct * 100;
feature3M(minr:maxr,minc:maxc,mins:maxs) = featureCropped3M;

% Write feature impact as a dose in CERR

% Create Texture Scans
[xVals, yVals, zVals] = getUniformScanXYZVals(planC{indexS.scan}(scanNum));
deltaXYZv(1) = abs(xVals(2)-xVals(1));
deltaXYZv(2) = abs(yVals(2)-yVals(1));
deltaXYZv(3) = abs(zVals(2)-zVals(1));
uniqueSlicesV = mins:maxs;
zV = zVals(uniqueSlicesV);
regParamsS.horizontalGridInterval = deltaXYZv(1);
regParamsS.verticalGridInterval   = -deltaXYZv(2); %(-)ve for dose
regParamsS.coord1OFFirstPoint   = xVals(minc);
regParamsS.coord2OFFirstPoint   = yVals(minr); % for dose
regParamsS.zValues  = zV;
regParamsS.sliceThickness = [planC{indexS.scan}(scanNum).scanInfo(uniqueSlicesV).sliceThickness];
assocScanUID = planC{indexS.structures}(structNum).assocScanUID;
dose2CERR(featureCropped3M,[], featureName, [featureName,'_impact'],...
    [featureName,'_impact'],'non CT',regParamsS,'no',assocScanUID)

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [AA,bb,p,rankA] = lusolCondense(A,b,mode,printLevel)
%
%        [AA,bb,p] = lusolCondense(A,b,mode);
% extracts from A a submatrix AA = A(p,:) of full row rank, and
% the corresponding subvector bb = b(p), where A should be sparse
% and have nonzero entries not too much bigger than 1.
%
%INPUT
% A     m x n sparse matrix
% b     m x 1 vector
%
%OPTIONAL INPUT
%    mode:     If mode=1, LUSOL operates on A itself.
%              If mode=2, LUSOL operates on A'.
%
%OUTPUT
% AA        row reduced A
% bb        row reduced B
% p         AA = A(p,:),bb = b(p)
% rankA     rank of A
%

% 27 May 2007: Michael Saunders:
%              First version of lusolCondense, intended as a
%              more efficient substitute for [Q,R,P] = qr(full(A')).
%              Michael Saunders
% 10 Oct  2019: Ronan Fleming
%              64bit linux support

if ~exist('mode','var')
    mode=1; %what is the best default.
end
if ~exist('printLevel','var')
    printLevel=0;
end

if ~issparse(A)
    A=sparse(A);
end

[m,n]   = size(A);
if isempty(b)
    b=sparse(m,1);
end

archstr = computer('arch');
archstr = lower(archstr);
switch archstr
    case {'win32', 'win64'}
        error('%s\n','lusolCondense not compatible with windows OS.')
    case {'glnxa86'}

        options = lusolSet;
        options.Pivoting  = 'TRP';
        options.FactorTol = 1.5;

        if mode==1   % Factorize A itself
            [L,U,p,q,options] = lusolFactor(A,options);
            inform   = options.Inform;
            rankA    = options.Rank;

            if rankA==m
                AA = A;
                bb = b;
                p  = (1:m)';
            else
                AA = A(p(1:rankA),:);
                bb = b(p(1:rankA));
            end

        else         % Factorize A'
            [L,U,p,q,options] = lusolFactor(A',options);
            inform  = options.Inform;
            rankA    = options.Rank;

            if rankA==m
                AA = A;
                bb = b;
                p  = (1:m)';
            else
                p  = q;
                AA = A(p(1:rankA),:);
                bb = b(p(1:rankA));
            end
        end

    case {'glnxa64','maci64'}
        if ~isempty(which('lusol_obj'))
            % generate default options
            options = lusol_obj.luset();
            % |--------+----------+----------------------------------------------------|
            % | param  |  default | description                                        |
            % |--------+----------+----------------------------------------------------|
            %
            % lusol_obj options
            % |--------+----------+----------------------------------------------------|
            % | nzinit |        0 | minimum length for storage arrays                  |
            % |--------+----------+----------------------------------------------------|
            %
            % LUSOL integer parameters
            % |--------+----------+----------------------------------------------------|
            % | maxcol |        5 | max num cols searched for piv element              |
            % | pivot  |    'TPP' | pivoting method {'TPP','TRP','TCP','TSP'}          |
            % | keepLU |        1 | keep the nonzeros, if 0, permutations are computed |
            % |--------+----------+----------------------------------------------------|
            %
            % LUSOL real parameters
            % |--------+----------+----------------------------------------------------|
            % | Ltol1  |     10.0 | max Lij allowed, default depends on pivot method   |
            % | Ltol2  |     10.0 | max Lij allowed during updates                     |
            % | small  |  eps^0.8 | absolute tolerance for treating reals as zero      |
            % | Utol1  | eps^0.67 | absolute tol for flagging small diags of U         |
            % | Utol2  | eps^0.67 | rel tol for flagging small diags of U              |
            % | Uspace |      3.0 |                                                    |
            % | dens1  |      0.3 |                                                    |
            % | dens2  |      0.5 |                                                    |
            % |--------+----------+----------------------------------------------------|

            % %modification of default options
            % options.pivot  = 'TRP';
            % options.Ltol1 = 1.5;
            % options.nzinit = 1e7;
            % %factorise
            % mylu = lusol_obj(A',options);

            if mode==1   % Factorize A itself
                %factorise
                mylu = lusol_obj(A);

                %extract results
                stats = mylu.stats();
                options.Inform = stats.inform;
                options.Nsing  = stats.nsing;
                options.Growth = stats.growth;

                %matrices
                L = mylu.L0();
                U = mylu.U();
                % row permutation
                p = mylu.p();
                % column permutation
                q = mylu.q();

                L     = L(p,p);      % New L is strictly lower triangular (and square).
                % U     = U(p,q);      % New U would be upper trapezoidal the size of S'.

                %return the rank of the matrix
                rankA=mylu.rank();

                %    lu.factorize inform codes:
                switch stats.inform
                    case 0
                        if printLevel>0
                            fprintf('%s\n','The LU factors were obtained successfully.')
                        end
                    case 1
                        if printLevel>0
                            fprintf('%s\n','U appears to be singular, as judged by lu6chk.');
                        end
                    case 3
                        fprintf('%s\n','Some index pair indc(l), indr(l) lies outside the matrix dimensions 1:m , 1:n.');
                    case    4
                        fprintf('%s\n','Some index pair indc(l), indr(l) duplicates another such pair.');
                    case    7
                        fprintf('%s\n','The arrays a, indc, indr were not large enough.');
                        fprintf('%s\n','Their length "lena" should be increase to at least');
                        fprintf('%s\n','The value "minlen" given in luparm(13).');
                    case 8
                        fprintf('%s\n','There was some other fatal error.  (Shouldn''t happen!)');
                    case 9
                        fprintf('%s\n','No diagonal pivot could be found with TSP or TDP.');
                        fprintf('%s\n','The matrix must not be sufficiently definite or quasi-definite.');
                end

                if stats.inform~=0 && stats.inform~=1
                    % solve Ax=b
                    b = ones(size(A',1),1);
                    x = mylu.solve(b);
                    % multiply Ax
                    b2 = mylu.mulA(x);
                    % check the result
                    fprintf('%s\t%g\n','Check norm(b-b2) : ', norm(b-b2))
                end

                if rankA==m
                    AA = A;
                    bb = b;
                    p  = (1:m)';
                else
                    AA = A(p(1:rankA),:);
                    bb = b(p(1:rankA));
                end

            else
                %factorise A'
                mylu = lusol_obj(A');

                %extract results
                stats = mylu.stats();
                options.Inform = stats.inform;
                options.Nsing  = stats.nsing;
                options.Growth = stats.growth;

                %matrices
                L = mylu.L0();
                U = mylu.U();
                % row permutation
                p = mylu.p();
                % column permutation
                q = mylu.q();

                L     = L(p,p);      % New L is strictly lower triangular (and square).
                % U     = U(p,q);      % New U would be upper trapezoidal the size of S'.

                %return the rank of the matrix
                rankA=mylu.rank();

                %    lu.factorize inform codes:
                switch stats.inform
                    case 0
                        if printLevel>0
                            fprintf('%s\n','The LU factors were obtained successfully.')
                        end
                    case 1
                        if printLevel>0
                            fprintf('%s\n','U appears to be singular, as judged by lu6chk.');
                        end
                    case 3
                        fprintf('%s\n','Some index pair indc(l), indr(l) lies outside the matrix dimensions 1:m , 1:n.');
                    case    4
                        fprintf('%s\n','Some index pair indc(l), indr(l) duplicates another such pair.');
                    case    7
                        fprintf('%s\n','The arrays a, indc, indr were not large enough.');
                        fprintf('%s\n','Their length "lena" should be increase to at least');
                        fprintf('%s\n','The value "minlen" given in luparm(13).');
                    case 8
                        fprintf('%s\n','There was some other fatal error.  (Shouldn''t happen!)');
                    case 9
                        fprintf('%s\n','No diagonal pivot could be found with TSP or TDP.');
                        fprintf('%s\n','The matrix must not be sufficiently definite or quasi-definite.');
                end

                if stats.inform~=0 && stats.inform~=1
                    % solve Ax=b
                    b = ones(size(A',1),1);
                    x = mylu.solve(b);
                    % multiply Ax
                    b2 = mylu.mulA(x);
                    % check the result
                    fprintf('%s\t%g\n','Check norm(b-b2) : ', norm(b-b2))
                end

                if rankA==m
                    AA = A;
                    bb = b;
                    p  = (1:m)';
                else
                    p  = q;
                    AA = A(p(1:rankA),:);
                    bb = b(p(1:rankA));
                end
            end

        else
            error('%s\n','lusolCondense.m Cannot find lusol_obj.m from lusol interface. Make sure https://github.com/nwh/lusol is intalled and added to the matlab path.')
        end
end
if printLevel>1
    fprintf('%s\n',['lusolCondense: A has dimensions ' int2str(m) ' x ' int2str(n) ' and has rank ' int2str(rankA)])
end

author: opencobra

repo: cobratoolbox


text:
function loglikelihood = ComputeLogLikelihood(P, G, dataset)
% returns the (natural) log-likelihood of data given the model and graph structure
%
% Inputs:
% P: struct array parameters (explained in PA description)
% G: graph structure and parameterization (explained in PA description)
%
%    NOTICE that G could be either 10x2 (same graph shared by all classes)
%    or 10x2x2 (each class has its own graph). your code should compute
%    the log-likelihood using the right graph.
%
% dataset: N x 10 x 3, N poses represented by 10 parts in (y, x, alpha)
%
% Output:
% loglikelihood: log-likelihood of the data (scalar)
%
% Copyright (C) Daphne Koller, Stanford Univerity, 2012

N = size(dataset,1); % number of examples
K = length(P.c); % number of classes

loglikelihood = 0;
% You should compute the log likelihood of data as in eq. (12) and (13)
% in the PA description
% Hint: Use lognormpdf instead of log(normpdf) to prevent underflow.
%       You may use log(sum(exp(logProb))) to do addition in the original
%       space, sum(Prob).
% assuming data = [y x alpha]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% YOUR CODE HERE
% function val = lognormpdf(x, mu, sigma);
if size(size(G),2) == 2
	G1 = G; G2 = G;
else
	G1 = reshape(G(:,:,1),10,2);
	G2 = reshape(G(:,:,2),10,2);
end

for i = 1:N
	data = reshape(dataset(i,:,:),10,3);
	log1 = log(P.c(1));
	log2 = log(P.c(2));
	for j = 1:10
		if G1(j,1) == 1   % construct g1 and g2 in the beginning from g
			%do shit
			theta1 = P.clg(j).theta(1,:);   theta2 = P.clg(j).theta(2,:);
			parent1 = data(G1(j,2),:);   parent2 = data(G2(j,2),:);
			mu_y1 = sum(theta1(1:4).*[1, parent1]);mu_y2 = sum(theta2(1:4).*[1, parent2]);
			mu_x1 = sum(theta1(5:8).*[1, parent1]);mu_x2 = sum(theta2(5:8).*[1, parent2]);
			mu_a1 = sum(theta1(9:12).*[1, parent1]);mu_a2 = sum(theta2(9:12).*[1, parent2]);
			log1 += lognormpdf(data(j,1),mu_y1,P.clg(j).sigma_y(1));
			log1 += lognormpdf(data(j,2),mu_x1,P.clg(j).sigma_x(1));
			log1 += lognormpdf(data(j,3),mu_a1,P.clg(j).sigma_angle(1));
			log2 += lognormpdf(data(j,1),mu_y2,P.clg(j).sigma_y(2));
			log2 += lognormpdf(data(j,2),mu_x2,P.clg(j).sigma_x(2));
			log2 += lognormpdf(data(j,3),mu_a2,P.clg(j).sigma_angle(2));

		else
			log1 += lognormpdf(data(j,1),P.clg(j).mu_y(1),P.clg(j).sigma_y(1));
			log1 += lognormpdf(data(j,2),P.clg(j).mu_x(1),P.clg(j).sigma_x(1));
			log1 += lognormpdf(data(j,3),P.clg(j).mu_angle(1),P.clg(j).sigma_angle(1));
			log2 += lognormpdf(data(j,1),P.clg(j).mu_y(2),P.clg(j).sigma_y(2));
			log2 += lognormpdf(data(j,2),P.clg(j).mu_x(2),P.clg(j).sigma_x(2));
			log2 += lognormpdf(data(j,3),P.clg(j).mu_angle(2),P.clg(j).sigma_angle(2));
		end
	end
	loglikelihood += log(sum(exp(log1)+exp(log2)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: anhncs

repo: Probabilistic-Graphical-Models


text:
classdef LabelIsosurfacesOptionsDialog < handle
%LABELISOSURFACESOPTIONSDIALOG Open a dialog for 3D label isosurfaces
%
%   Class LabelIsosurfacesOptionsDialog
%
%   Example
%   LabelIsosurfacesOptionsDialog
%
%   See also
%     IsosurfaceOptionsDialog
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2013-04-19,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2013 INRA - Cepia Software Platform.

%% Properties
properties
    Parent;

    Handles;
end % end properties

%% Constructor
methods
    function obj = LabelIsosurfacesOptionsDialog(parent, varargin)
    % Constructor for LabelIsosurfacesOptionsDialog class

         % call parent constructor to initialize members
        obj = obj@handle();
        obj.Parent = parent;

        % create default figure
        fig = figure(...
            'MenuBar', 'none', ...
            'NumberTitle', 'off', ...
            'HandleVisibility', 'On', ...
            'Name', 'Label Isosurfaces Options', ...
            'CloseRequestFcn', @obj.close);
        obj.Handles.Figure = fig;

        setupLayout(fig);

        function setupLayout(hf)
            % Setup layout of figure widgets

            pos = get(hf, 'Position');
            pos(3:4) = [200 180];
            set(hf, 'Position', pos);
            set(hf, 'Units', 'normalized');

            if verLessThan('matlab', 'R2014b')
                mainPanel = uiextras.VBox('Parent', hf);
            else
                mainPanel = uix.VBox('Parent', hf);
            end

            % control panel: use vertical layout
            if verLessThan('matlab', 'R2014b')
                optionsPanel = uiextras.VButtonBox('Parent', mainPanel, ...
                    'ButtonSize', [150 30], ...
                    'VerticalAlignment', 'top', ...
                    'HorizontalAlignment', 'center', ...
                    'Spacing', 5, 'Padding', 5);
            else
                optionsPanel = uix.VButtonBox('Parent', mainPanel, ...
                    'ButtonSize', [150 30], ...
                    'VerticalAlignment', 'top', ...
                    'HorizontalAlignment', 'center', ...
                    'Spacing', 5, 'Padding', 5);
            end

            obj.Handles.ReverseZAxis = uicontrol(...
                'Style', 'checkbox', 'Parent', optionsPanel, ...
                'String', 'Reverse Z-axis', ...
                'Value', 1, ...
                'HorizontalAlignment', 'Left');

            obj.Handles.RotateOx = uicontrol(...
                'Style', 'checkbox', 'Parent', optionsPanel, ...
                'String', 'Rotate around X-axis', ...
                'Value', 1, ...
                'HorizontalAlignment', 'Left');

            obj.Handles.Smooth = uicontrol(...
                'Style', 'checkbox', 'Parent', optionsPanel, ...
                'String', 'Smooth', ...
                'Value', 1, ...
                'HorizontalAlignment', 'Left');

            obj.Handles.ShowAxisLabel = uicontrol(...
                'Style', 'checkbox', 'Parent', optionsPanel, ...
                'String', 'Show axes label', ...
                'Value', 0, ...
                'HorizontalAlignment', 'Left');

            % add control panel with "Compute" and "Close" buttons
            if verLessThan('matlab', 'R2014b')
                buttonsPanel = uiextras.HButtonBox('Parent', mainPanel, ...
                    'ButtonSize', [70 30], ...
                    'VerticalAlignment', 'top', ...
                    'Spacing', 5, 'Padding', 5, ...
                    'Units', 'normalized', ...
                    'Position', [0 0 1 1]);
            else
                buttonsPanel = uix.HButtonBox('Parent', mainPanel, ...
                    'ButtonSize', [70 30], ...
                    'VerticalAlignment', 'top', ...
                    'Spacing', 5, 'Padding', 5, ...
                    'Units', 'normalized', ...
                    'Position', [0 0 1 1]);
            end

            obj.Handles.ApplyButton = uicontrol('Style', 'PushButton', ...
                'Parent', buttonsPanel, ...
                'String', 'Compute', ...
                'Enable', 'on', ...
                'Callback', @obj.onApplyButtonClicked);
            obj.Handles.CloseButton = uicontrol('Style', 'PushButton', ...
                'Parent', buttonsPanel, ...
                'String', 'Close', ...
                'Callback', @obj.close);

            if verLessThan('matlab', 'R2014b')
                mainPanel.Sizes = [-1 40];
            else
                mainPanel.Heights = [-1 40];
            end

        end
    end

end % end constructors

%% Methods
methods
    function onApplyButtonClicked(obj, hObject, eventdata) %#ok<INUSD>

        % first, disable "compute" button to avoid multiple calls
        set(obj.Handles.ApplyButton, 'enable', 'off');
        hDlg = msgbox({'Computing isosurfaces,', 'Please wait...'}, ...
            'Isosurface computing');

        % extract options from widgets
        reverseZAxis = get(obj.Handles.ReverseZAxis, 'Value');
        rotateXAxis = get(obj.Handles.RotateOx, 'Value');
        smooth = get(obj.Handles.Smooth, 'Value');
        showAxesLabel = get(obj.Handles.ShowAxisLabel, 'Value');

        % get image data stored in parent Slicer
        imgData = obj.Parent.ImageData;
        imgSize = obj.Parent.ImageSize;

        % avoid the case of empty image
        if isempty(imgData)
            return;
        end

        % eventually rotate image around X-axis
        if rotateXAxis
            imgData = stackRotate90(imgData, 'x', 1);
            imgSize = imgSize([1 3 2]);
        end

        % compute display settings
        spacing = obj.Parent.VoxelSize;
        origin  = obj.Parent.VoxelOrigin;

        % compute grid positions
        lx = (0:imgSize(1)-1) * spacing(1) + origin(1);
        ly = (0:imgSize(2)-1) * spacing(2) + origin(2);
        lz = (0:imgSize(3)-1) * spacing(3) + origin(3);

        if rotateXAxis
            ly = ly(end:-1:1);
        end

        % number of different labels
        nLabels = double(max(imgData(:)));

        % determine the color map to use (default is empty)
        cmap = [];
        if ~isColorStack(imgData) && ~isempty(obj.Parent.ColorMap)
            cmap = obj.Parent.ColorMap;
        end
        if isempty(cmap)
            cmap = jet(256);
        end

        % if colormap has 256 entries, we need only a subset
        % otherwise, we assume colormap has as many rows as the number
        % of labels.
        if size(cmap, 1) == 256
            inds = round(linspace(2, 256, nLabels));
            cmap = cmap(inds, :);
        end

        % for binary images, colormap has two entries, so we need to keep
        % the last one
        if strcmp(obj.Parent.ImageType, 'binary') && size(cmap, 1) > 1
            cmap = cmap(2,:);
        end

        % create figure
        hf = figure(); hold on;
        set(hf, 'renderer', 'opengl');

        % compute an isosurface for each label
        for i = 1:nLabels
            im = imgData==i;
            if ~any(im)
                continue;
            end

            if smooth
                im = imGaussianFilter(double(im), [5 5 5], 2);
            end

            % display isosurface
            p = patch(isosurface(lx, ly, lz, im, .5));

            % set face color
            set(p, 'FaceColor', cmap(i,:), 'EdgeColor', 'none');
        end

        % compute display extent (add a 0.5 limit around each voxel)
        extent = stackExtent(imgSize, spacing, origin);

        % setup display
        axis equal;
        axis(extent);
        view(3);
        axis('vis3d');

        light;

        if rotateXAxis || reverseZAxis
            set(gca, 'zdir', 'reverse');
        end

        if showAxesLabel
            xlabel('X axis');
            if rotateXAxis
                ylabel('Z Axis');
                zlabel('Y Axis');
            else
                ylabel('Y Axis');
                zlabel('Z Axis');
            end
        end

        % re-enable "apply" button
        set(obj.Handles.ApplyButton, 'enable', 'on');
        if ishandle(hDlg)
            close(hDlg);
        end

    end

end % end methods

%% Figure management
methods
    function close(obj, varargin)
        delete(obj.Handles.Figure);
    end

end

end % end classdef

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
% MODSELC Model selection
%
%   V = MODSELC(A,W,N,NREP)
%   V = A*(W*MODSELC([],N,NREP))
%
% INPUT
%   A      Dataset used for training base classifiers and/or selection
%   B      Dataset used for testing (executing) the selector
%   W      Set of trained or untrained base classifiers
%   N      Number of crossvalidations, default 10
%   NREP   Number of crossvalidation repetitions
%
% OUTPUT
%   V   Selected classidfer
%
% DESCRIPTION
% This routine selects out of a set of given classifiers stored in W the
% best one on the basis of N-fold crossvalidation (see CROSSVAL), which
% might be repeated NREP times. If W contains a set of already trained
% classifiers, N and NREP are neglected and just the best classifier
% according to the evaluation set A is returned.
%
% This routine can be considered as a classifier combiner based on global
% selection. See DCSC for local, dynamic selection.
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS, STACKED, DCSC, CLASSC, TESTD, LABELD

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function OUT = modselc(par1,par2,par3,par4)

name = 'Classifier Model Selection';
DefaultN = 10;
DefaultNrep = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Empty Call MODSELC, MODSELC([],N,NREP)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin < 1 | isempty(par1)
	if nargin > 1 & (ismapping(par2) | iscell(par2))
		if nargin < 4 | isempty(par4), par4 = DefaultNrep; end
		if nargin < 3 | isempty(par3), par3 = DefaultN; end
		OUT = prmapping(mfilename,'untrained',{par2,par3,par4});
	else
		if nargin < 3 | isempty(par3), par3 = DefaultNrep; end
		if nargin < 2 | isempty(par2), par2 = DefaultN; end
		% If there are no inputs, return an untrained mapping.
		% (PRTools transfers the latter into the first)
		OUT = prmapping(mfilename,'combiner',{par2,par3});
	end
	OUT = setname(OUT,name);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Storing Base Classifiers:  W*MODSELC, MODSELC(W,N,NREP)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif ismapping(par1)
	if nargin < 3 | isempty(par3), par3 = DefaultNrep; end
	if nargin < 2 | isempty(par2), par2 = DefaultN; end
	% call like OUT = MODSELC(W,N) or OUT = W*MODSELC([],N)
	% store trained or untrained base classifiers,
	% ready for later training of the combiner
	BaseClassf = par1;
	if ~isparallel(BaseClassf) & ~isstacked(BaseClassf)
		error('Parallel or stacked set of base classifiers expected');
	end
	OUT = prmapping(mfilename,'untrained',{BaseClassf,par2,par3});
	OUT = setname(OUT,name);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Selection (and training base classifiers if needed):
%   A*(W*MODSELC), A*(W*MODSELC([],N,NREP)), MODSELC(A,W,N,NREP)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif isdataset(par1) & ismapping(par2) & ...
		(isuntrained(par2) | (istrained(par2) & (isstacked(par2) | isparallel(par2))))
	% call like OUT = MODSELC(TrainSet,W,N,NREP) or TrainSet*(W*MODSELC([],N,NREP))
	% (PRTools transfers the latter into the first)
	% W is a set of trained or untrained set classifiers.
	if nargin < 4 | isempty(par4), par4 = DefaultNrep; end
	if nargin < 3 | isempty(par3), par3 = DefaultN; end
	TrainSet = par1;
	BaseClassf = par2;
	N = par3;
	Nrep = par4;
	isvaldfile(TrainSet,1,2); % at least one object per class, 2 classes
	TrainSet = setprior(TrainSet,getprior(TrainSet,0));  % avoid many warnings
	if ~(isstacked(BaseClassf) | isparallel(BaseClassf))
		if iscombiner(BaseClassf)
			error('No base classifiers found')
		end
		Data = getdata(BaseClassf);
		BaseClassf = Data.BaseClassf; % base classifiers were already stored
  end
  if isuntrained(BaseClassf)  % base classifiers untrained, crossval needed
		randstate = randreset; % makes routine reproducing
		if isparallel(BaseClassf)
			BaseClassf = getdata(BaseClassf);
			tsize = BaseClassf{end};
			if ismapping(tsize)
				error(['Training of parallel combined untrained classifier not possible.' ...
						newline 'Feature sizes should be stored in the classifier first.'])
			end
			csize = cumsum([0 tsize(:)']);
			n = length(BaseClassf)-1;
			e = zeros(1,n);
			rstate = randreset;
			s = sprintf('Crossvalidating %i base classifiers: ',n);
			prwaitbar(n,s)
			for j=1:n
				prwaitbar(n,j,[s getname(BaseClassf{j})]);
				randreset(1);
				e(j) = crossval(TrainSet(:,csize(j)+1:csize(j+1)),BaseClassf{j},N,Nrep);
			end
			prwaitbar(0);
			[ee,k] = min(e);
			J = [csize(k)+1:csize(k+1)];
			OUT = TrainSet(:,J)*BaseClassf{k};
			OUT = featsel(csize(end),J)*OUT;
		else
			if isstacked(BaseClassf)
				BaseClassf = getdata(BaseClassf);
			end
			randreset(1);
			e = crossval(TrainSet,BaseClassf,N,Nrep);
			[ee,k] = min(e);
			OUT = TrainSet*BaseClassf{k};
		end
		randreset(randstate);
	else  % base classifiers are trained, check label lists and find best one
		BaseClassf = getdata(BaseClassf);
		n = length(BaseClassf);
		for j=1:n
    	if ~isequal(getlabels(BaseClassf{j}),getlablist(TrainSet))
     		error('Training set and base classifiers should deal with same labels')
    	end
		end
		e = testc(TrainSet,BaseClassf);
		[ee,k] = min([e{:}]);
		OUT = BaseClassf{k};
  end

else

  error('Illegal input');

end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [x,y] = snakedeform(x,y,alpha,beta,gamma,kappa,fx,fy,ITER)
% -------In this function, the initial contour of Active Coutour Model(Snake)
%  will be deformed in the given external force field.
%
%     alpha:   elasticity parameter
%     beta:    rigidity parameter
%     gamma:   viscosity parameter
%     kappa:   external force weight
%     fx,fy:   external force field

N = length(x);

alpha = alpha* ones(1,N);
beta = beta*ones(1,N);

% produce the five diagnal vectors
alpham1 = [alpha(2:N) alpha(1)];
alphap1 = [alpha(N) alpha(1:N-1)];
betam1 = [beta(2:N) beta(1)];
betap1 = [beta(N) beta(1:N-1)];

a = betam1;
b = -alpha - 2*beta - 2*betam1;
c = alpha + alphap1 +betam1 + 4*beta + betap1;
d = -alphap1 - 2*beta - 2*betap1;
e = betap1;

% generate the parameters matrix
A = diag(a(1:N-2),-2) + diag(a(N-1:N),N-2);
A = A + diag(b(1:N-1),-1) + diag(b(N), N-1);
A = A + diag(c);
A = A + diag(d(1:N-1),1) + diag(d(N),-(N-1));
A = A + diag(e(1:N-2),2) + diag(e(N-1:N),-(N-2));

invAI = inv(A + gamma * diag(ones(1,N)));

for count = 1:ITER,
   vfx = interp2(fx,x,y,'*linear');
   vfy = interp2(fy,x,y,'*linear');

   % deform snake
   x = invAI * (gamma* x + kappa*vfx);
   y = invAI * (gamma* y + kappa*vfy);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [y,newmap] = cmsort(x,map,opt)
%CMSORT   Sort color map.
%
%   NEWMAP = CMSORT(MAP) sorts the color map MAP by luminance.
%
%   NEWMAP = CMSORT(MAP,OPT) where OPT is one of the strings 'red',
%   'green', 'blue', 'hue', 'saturation', 'value' and 'luminance', sorts
%   the map according to the specified criterion.
%
%   [Y,NEWMAP] = CMSORT(X,MAP) and [Y,NEWMAP] = CMSORT(X,MAP,OPT) also
%   changes the index matrix so the image remains unchanged after the
%   sorting of the color map.
%
%   The case of the option string does not matter and the string may be
%   truncated since only the first letter is used internally.

%   Author:      Peter J. Acklam
%   Time-stamp:  1998-04-15 13:42:03
%   E-mail:      jacklam@math.uio.no (Internet)
%   URL:         http://www.math.uio.no/~jacklam

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check and identify input arguments.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ( nargin < 1 )
   error( 'Not enough input arguments.' );
elseif ( nargin == 1 )          % CMSORT(MAP)
   map = x;
   opt = 'lum';
   index_given = 0;
elseif ( nargin == 2 )
   if ischar(map)               % CMSORT(MAP,OPT)
      map = x;
      opt = map;
      index_given = 0;
      if ( nargout > 1 )
         error( 'Too many output arguments.' );
      end
   else                         % CMSORT(X,MAP)
      opt = 'lum';
      index_given = 1;
   end
elseif ( nargin == 3 )          % CMSORT(X,MAP,OPT)
   index_given = 1;
else
   error( 'Too many input arguments.' );
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now create the key by which the color map is sorted.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

swtch = lower( opt(1) );

if strcmp( swtch, 'r' )         % Red.
   key = map(:,1);
elseif strcmp( swtch, 'g' )     % Green.
   key = map(:,2);
elseif strcmp( swtch, 'b' )     % Blue.
   key = map(:,3);
elseif strcmp( swtch, 'h' )     % Hue.
   hsv = rgb2hsv( map );
   key = hsv(:,1);
elseif strcmp( swtch, 's' )     % Saturation.
   hsv = rgb2hsv( map );
   key = hsv(:,2);
elseif strcmp( swtch, 'v' )     % Value.
   hsv = rgb2hsv( map );
   key = hsv(:,3);
elseif strcmp( swtch, 'l' )     % Luminance.
   w = [ 0.298936 ; 0.587043 ; 0.114021 ];      % RGB weights.
   key = map*w;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sort color map and fix index matrix (if given).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[ dummy, idx ] = sort( key );   % Sort the keys.
newmap = map(idx,:);            % Rearrange color map.

if index_given
   n = length( idx );           % Number of colors.
   xdi = zeros( n, 1 );
   xdi(idx) = 1:n;
   y = zeros( size( x ) );
   y(:) = xdi(x);
else
   y = map;
end