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

global planC stateS;
indexS = planC{end};

    [originF, spacingF, centerF] = getScanOriginSpacing(planC{indexS.scan}(stateS.imageRegistrationBaseDataset));

    [originM, spacingM, centerM] = getScanOriginSpacing(planC{indexS.scan}(stateS.imageRegistrationMovDataset));

    Bins = str2double(get(handles.mattes_bins, 'string'));
    Samples = str2double(get(handles.mattes_samples, 'string'));
    RelaxFactor = str2double(get(handles.mattes_rf, 'string'));
    DefaultPixelValue = str2double(get(handles.mattes_pv, 'string'));
    Iter_MinStep = str2double(get(handles.mattes_minstep, 'string'));
    Iter_MaxStep = str2double(get(handles.mattes_maxstep, 'string'));
    Iter_Num = str2double(get(handles.mattes_inum, 'string'));

    tscale = str2double(get(handles.mattes_tscale, 'string'));
    rscale = str2double(get(handles.mattes_rscale, 'string'));
    scaleFactor = str2double(get(handles.mattes_scalefactor, 'string'));

    output = cell(1, 16);

    FImg = planC{indexS.scan}(stateS.imageRegistrationBaseDataset).scanArray;
    MImg = planC{indexS.scan}(stateS.imageRegistrationMovDataset).scanArray;

    dimF = size(FImg);
    dimM = size(MImg);
    %generate cliped base dataset
    clipBox = uint16(getappdata(stateS.handle.CERRSliceViewer, 'clipBox_baseTrans'));
    if ~isempty(clipBox)
        FImg = FImg(clipBox(2):clipBox(4), clipBox(1):clipBox(3), :);
        originF(1) = originF(1) + spacingF(1)*double(clipBox(1)-1);
        originF(2) = originF(2) + spacingF(2)*double(uint16(dimF(1)-clipBox(4)));
    end
    clipBox = uint16(getappdata(stateS.handle.CERRSliceViewer, 'clipBox_baseSag'));
    if ~isempty(clipBox)
        FImg = FImg(:, :, clipBox(2):clipBox(4));
        originF(3) = originF(3) + spacingF(3)*double(clipBox(2)-1);
    end

    %generate cliped move dataset
    clipBox = uint16(getappdata(stateS.handle.CERRSliceViewer, 'clipBox_movTrans'));
    if ~isempty(clipBox)
        MImg = MImg(clipBox(2):clipBox(4), clipBox(1):clipBox(3), :);
        originM(1) = originM(1) + spacingM(1)*double(clipBox(1)-1);
        originM(2) = originM(2) + spacingM(2)*double(uint16(dimM(1)-clipBox(4)));
    end
    clipBox = uint16(getappdata(stateS.handle.CERRSliceViewer, 'clipBox_movSag'));
    if ~isempty(clipBox)
        MImg = MImg(:, :, clipBox(2):clipBox(4));
        originM(3) = originM(3) + spacingM(3)*double(clipBox(2)-1);
    end

%downsample the input datasets
    dsampled = 0;
    if (get(handles.dsampleCheck,'Value') == get(handles.dsampleCheck,'Max'))

        tic;
        xyScale =  2; zScale = 2;
        spacingF = [spacingF(1)*xyScale spacingF(2)*xyScale spacingF(3)*zScale];
        spacingM = [spacingM(1)*xyScale spacingM(2)*xyScale spacingM(3)*zScale];
        disp('downSampling ...');
%         FImg = imdesample3d(FImg,xyScale,zScale);
%         MImg = imdesample3d(MImg,xyScale,zScale);
%         FImg=GPReduce2D(FImg,1);
%         MImg=GPReduce2(MImg,1);
        FImg = FImg(1:xyScale:end, 1:xyScale:end, 1:zScale:end);
        MImg = MImg(1:xyScale:end, 1:xyScale:end, 1:zScale:end);

        toc;
        dsampled = 1;
    end

    tic;

%flip the source datasets on X for itk coordinate system
    fdim = 1;
    FImg = flipdim(FImg, fdim);
    MImg = flipdim(MImg, fdim);

% transform initializing modes 1:MomentON 0:GeometryOn 2:initail transform On
    initMode = get(handles.InitTrans, 'value');

%prepare the initial transform matrix
    transMB = planC{indexS.scan}(stateS.imageRegistrationBaseDataset).transM;
    transMM = planC{indexS.scan}(stateS.imageRegistrationMovDataset).transM;
    if isempty(transMM), transMM = eye(4); end;
    if isempty(transMB), transMB = eye(4); end;
    transM = inv(transMB) * transMM;

%flip the moving dataset to match the initial transM
    Tf = eye(4);
    flipM = 0;
    if (get(handles.flipX, 'value'))
        MImg = flipdim(MImg, 2);
        Tf = [-1 0 0 2*centerM(1); 0 1 0 0; 0 0 1 0; 0 0 0 1];
        flipM = 2;
    end
    if (get(handles.flipY, 'value'))
        MImg = flipdim(MImg, 1);
        Tf = [1 0 0 0; 0 -1 0 2*centerM(2); 0 0 1 0; 0 0 0 1];
        flipM = 1;
    end
    if (get(handles.flipZ, 'value'))
        MImg = flipdim(MImg, 3);
        Tf = [1 0 0 0; 0 1 0 0; 0 0 -1 2*centerM(3); 0 0 0 1];
        flipM = 3;
    end
    transM = transM * inv(Tf);
    rotM = transM(1:3, 1:3); transV = transM(1:3, 4);
    rotM = rotM';
    transV = -transV;

%call registration method
    switch metric
        case 'mutual information'
            metricIndex = 0;

            try
                try
                [im, Rotation, Offset] = Mattes_MI3D32_2006a(int16(FImg), originF, spacingF, ...
                                         int16(MImg), originM, spacingM, ...
                                         Bins, Samples, RelaxFactor, DefaultPixelValue, ...
                                         Iter_MinStep, Iter_MaxStep, Iter_Num, tscale, rscale, scaleFactor, ...
                                         rotM, transV, initMode);
                catch
                [im, Rotation, Offset] = Mattes_MI3D32_2006b(int16(FImg), originF, spacingF, ...
                                         int16(MImg), originM, spacingM, ...
                                         Bins, Samples, RelaxFactor, DefaultPixelValue, ...
                                         Iter_MinStep, Iter_MaxStep, Iter_Num, tscale, rscale, scaleFactor, ...
                                         rotM, transV, initMode);
                end
            catch
                try
                    [im, Rotation, Offset] = Mattes_MI3D6406b(uint16(FImg), originF, spacingF, ...
                                         int16(MImg), originM, spacingM, ...
                                         Bins, Samples, RelaxFactor, DefaultPixelValue, ...
                                         Iter_MinStep, Iter_MaxStep, Iter_Num, tscale, rscale, scaleFactor, ...
                                         rotM, transV, initMode);
                catch
                    try
                        [im, Rotation, Offset] = Mattes_MI3D6406a(uint16(FImg), originF, spacingF, ...
                                         int16(MImg), originM, spacingM, ...
                                         Bins, Samples, RelaxFactor, DefaultPixelValue, ...
                                         Iter_MinStep, Iter_MaxStep, Iter_Num, tscale, rscale, scaleFactor, ...
                                         rotM, transV, initMode);
                    catch
                        [im, Rotation, Offset] = Mattes_MI3D_debug(uint16(FImg), originF, spacingF, ...
                                         int16(MImg), originM, spacingM, ...
                                         Bins, Samples, RelaxFactor, DefaultPixelValue, ...
                                         Iter_MinStep, Iter_MaxStep, Iter_Num, tscale, rscale, scaleFactor, ...
                                         rotM, transV, initMode);
                    end
                end
            end
        case 'viola well'
            metricIndex = 1;
    end
    im = flipdim(im, fdim);
%??     if (flipM), im = flipdim(im, flipM); end;

    toc;

    output{1} = ['Rotation Angle Z = ' num2str(asin(Rotation(1, 2))*45.0/atan(1.0))];
    output{2} = ['Translation X = ' num2str(Offset(1))];
    output{3} = ['Translation Y = ' num2str(Offset(2))];
    output{4} = ['Translation Z = ' num2str(Offset(3))];

    output{5} = ['Iterations    = ' num2str(Offset(4))];
    output{6} = ['Metric Value (MI) = ' num2str(Offset(5))];

    output{7} = ['Rotation Center X = ' num2str(Offset(6))];
    output{8} = ['Rotation Center Y = ' num2str(Offset(7))];
    output{9} = ['Rotation Center Z = ' num2str(Offset(8))];

    output{10} = ['Scale = ' num2str(Offset(9))];

    output{11} = ['Offset X = ' num2str(Offset(10))];
    output{12} = ['Offset Y = ' num2str(Offset(11))];
    output{13} = ['Offset Z = ' num2str(Offset(12))];

    set(handles.OutputList, 'string', output);

%update the transM;
    scale = Offset(9);
    Tv = [Offset(1) Offset(2) Offset(3)]; %input image is a flip of source along X direction.
    Cv = Offset(6:8);
    Cv(3) = -Cv(3); %CERR z axis is opposite to the ITK.

    rot = reshape(Rotation, 3,3);
    offset = Offset(10:12); %offset related to origin rotation;

%     rot(2,1) = -rot(2,1); rot(1,2) = -rot(1,2);
%     rot(1,3) = -rot(1,3); rot(3,1) = -rot(3,1);
%     rot(2,3) = -rot(2,3); rot(3,2) = -rot(3,2);

%     offset(1) = -offset(1);
%     offset(2) = -offset(2);
%     offset(3) = -offset(3);

    TM = eye(4);
    TM(:,4) = [offset 1];

    RM = eye(4);
    RM(1:3, 1:3) = rot;

    newTransform = inv(TM*RM);
    newTransform =  transMB * newTransform * Tf;

    scanSetM = stateS.imageRegistrationMovDataset;
    scanSetF = stateS.imageRegistrationBaseDataset;
    planC{indexS.(stateS.imageRegistrationMovDatasetType)}(scanSetM).transM = newTransform;

% save the resampled dataset
    if (get(handles.saveCheck, 'value'))
        if (~dsampled)
            planC{indexS.scan}(end+1) = planC{indexS.scan}(scanSetF);
            planC{indexS.scan}(end).scanArray = im;
            planC{indexS.scan}(end).scanUID = createUID('scan');
        else
            %imReg = resample(scanSetF, scanSetM); % need to up-sample im ????
            planC{indexS.scan}(end+1) = planC{indexS.scan}(scanSetF);
            planC{indexS.scan}(end).scanArray = imReg;
            planC{indexS.scan}(end).scanUID = createUID('scan');
        end
        controlFrame('refresh');
    end

    sliceCallBack('refresh');

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% amplitude spectra
% imagery and movement

clc; close all; clear all;

dd='C:\Users\Doyeunlee\Desktop\Journal_dylee\1_ConvertedData\';
fsldx=2;
fs={'100','250','1000'};
ref='FCz';
channelMatrix={'F3','F1','Fz','F2','F4';
    'FC3','FC1','FCz','FC2','FC4';
    'C3','C1', 'Cz', 'C2', 'C4';
    'CP3','CP1','CPz','CP2','CP4'
    'P3','P1','Pz','P2','P4'};

% Motor Imagery
% filelist={'sub01','sub02','sub03','sub04','sub05','sub06','sub07','sub08','sub09','sub10'};
filelist={'sub01'};

for sub = 1:length(filelist)
    ival = [0 3000];
    % Data load - MI
    [cntReach,mrkReach,mntReach]=eegfile_loadMatlab([dd 'MI' '\' fs{fsldx} '\' filelist{sub} '_reaching_' 'MI']);
    % Data load - realMove
    [cntReach_ME,mrkReach_ME,mntReach_ME]=eegfile_loadMatlab([dd 'realMove' '\' fs{fsldx} '\' filelist{sub} '_reaching_' 'realMove']);

    % butter
    cntReach=proc_filtButter(cntReach,5,[4 40]);
    cntReach_ME=proc_filtButter(cntReach_ME,5,[4 40]);

    epoReach=cntToEpo(cntReach,mrkReach,ival);
    epoReach_ME=cntToEpo(cntReach_ME,mrkReach_ME,ival);

    epoReach=proc_selectChannels(epoReach,{'F3','F1','Fz','F2','F4',...
        'FC3','FC1','FCz','FC2','FC4',...
        'C3','C1', 'Cz', 'C2', 'C4', ...
        'CP3','CP1','CPz','CP2','CP4',...
        'P3','P1','Pz','P2','P4'});

    epoReach_ME=proc_selectChannels(epoReach_ME,{'F3','F1','Fz','F2','F4',...
        'FC3','FC1','FCz','FC2','FC4',...
        'C3','C1', 'Cz', 'C2', 'C4', ...
        'CP3','CP1','CPz','CP2','CP4',...
        'P3','P1','Pz','P2','P4'});

    fs = 250;

    %
    plot(psd(spectrum.periodogram, cntReach_ME.x, 'Fs', fs));
    psdest = psd(spectrum.periodogram, cntReach_ME.x, 'Fs', fs);
    plot(psdest.Frequencies,psdest.Data);
    xlabel('Hz'); grid on;
    hold on;
    psdest = psd(spectrum.periodogram, cntReach.x, 'Fs', fs);
    plot(psdest.Frequencies,psdest.Data);
    xlabel('Hz'); grid on;

    %
    xdft = fft(cntReach_ME.x);
    xdft = xdft(1:length(cntReach_ME.x)/2+1);
    freq = 0:fs/length(cntReach_ME.x):fs/2;
    figure();
    plot(freq,abs(xdft));
    xlabel('Hz');

    %
    y= fft(cntReach.x);
    n = length(cntReach.x);
    f = (0:n-1)*(fs/n);
    power = abs(y).^2/n;

    plot(f,power);
    xlabel('Frequency');
    ylabel('Power');

    y= fft(cntReach_ME.x);
    n = length(cntReach_ME.x);
    f = (0:n-1)*(fs/n);
    power = abs(y).^2/n;

    plot(f,power);
    xlabel('Frequency');
    ylabel('Power');

end

author: PatternRecognition

repo: OpenBMI


text:
function [DataOut, Metric] = ReceiverSD(SDchips, G, Scrambler);
%
% ReceiverSD    Soft Decision Receiver. This function performs transmitted packet data recovering,
%                          based on Viterbi Decoder.
%
% 						Block Diagram:
%						SDchips -> [DeScrambler] -> [De-Interleaver] -> [Viterbi Decoder]
%
%
% 						Inputs: SDchips - Soft Decision RAKE receiver symbols
%                         G - Viterbi Encoder generation polynom
%                         Scrambler - Scrambler sequence
%
% 						Outputs: DataOut - Received data (binary form)
%                                       Metric - The  best metric of Viterbi Decoder
%

%====================== COMMON DEFINITIONS =================
%------- Viterbi Polynom ---------
if (nargin == 1)
   G = [1 1 1 1 0 1 0 1 1; 1 0 1 1 1 0 0 0 1];
end

%====================== R E C E I V E R =================

%----------- DeScrambler---------
% Rate = 19.2 KBps
SDchips = SDchips.*sign(1/2-Scrambler);

%-------- DeInterleaver ---------
INTERL = reshape(SDchips, 16, 24);		% IN-> rows, OUT-> columns
SDchips = reshape(INTERL', length(SDchips), 1);  % Rate = 19.2 KBps

%-------- SD Input Viterbi Decoder ----------
[DataOut Metric] = SoftVitDec(G, SDchips, 1);

author: Sable

repo: mcbench-benchmarks


text:
function plot_cdfZ(vResults)

figure;i=1;cdfplot(reshape(vResults(i).mValueGrid,size(vResults(i).mValueGrid,1)*size(vResults(i).mValueGrid,2),1))
hold on;i=2;cdfplot(reshape(vResults(i).mValueGrid,size(vResults(i).mValueGrid,1)*size(vResults(i).mValueGrid,2),1))
hold on;i=3;cdfplot(reshape(vResults(i).mValueGrid,size(vResults(i).mValueGrid,1)*size(vResults(i).mValueGrid,2),1))

legend('a','b','c')
h=get(gca,'Children');

set(h(3),'Color',[0.8 0.8 0.8],'LineWidth',4)
set(h(2),'Color',[0 0 1],'LineWidth',2)
set(h(1),'Color',[1 0 0],'LineWidth',1)
set(gca,'YLim',[-0.05 1.05])

author: CelsoReyes


text:
function plot_boundary_orange(V,F)
%PLOT_BOUNDARY_ORANGE  plots the boundary vertices of a triangle mesh in
%orange, and its interior in blue.
%
% plot_boundary_orange(V,F)(V,F);
%
% Input:
%  V,F  mesh whose boundary is to be plotted

%Compute a function u that is 0 on all interior vertices and 1 on all
%boundary vertices.
u = ...

%Plot the function u with the right colors.
t = tsurf(F,V, 'CData',u);
shading interp;
axis equal;
axis off;
cm = flipud(cbrewer('RdYlBu', 500));
colormap(cm(100:450,:));
light('Position',[-1.5 1 1],'Style','local');
lights = camlight;
set(t, 'FaceLighting','gouraud', 'FaceColor','interp');
set(t, 'DiffuseStrength',0.5, 'SpecularStrength',0.2, 'AmbientStrength',0.3);
camproj('perspective');
add_shadow([t],lights);

end

author: odedstein

repo: sgi-introduction-course


text:
function [superPixSeg,superPixInd,superPixBounds]=mvg_FelzenswalbSuperpixelWrapper(img,sigma,min_area,k)
% The funtion [superPixSeg,superPixInd,superPixBounds]=FelzenswalbSuperpixels(img,sigma,min_area,k)
% is a wrapper funtion to superpixel segmentation algorithm from
% Pedro F. Felzenszwalb and Daniel P. Huttenlocher, described in
% "Efficient Graph-Based Image Segmentation"
% Pedro F. Felzenszwalb and Daniel P. Huttenlocher
% International Journal of Computer Vision, 59(2) September 2004.
% The executables are available online in
% See http://www.cs.brown.edu/~pff/segment/
%

%% Path to executables of Felzenszwalb superpixel algorithm (see above)
global configjson
soft_dir=fullfile(configjson.rahtu.rahtuPath,'/segment/'); % Change here the path to the executables

%% Default parameters
% Default parameters proposed by Felzenszwalb et al.  sigma = 0.5, k = 500, min_area = 20. (Larger values for k result in larger components in the result.)
if nargin<2
    sigma=0.5;
end
if nargin<3
    min_area=20;
end
if nargin<4
    k=500;
end

%% Initialize path
tmp_dir=soft_dir;
%tmp_dir=[soft_dir,'tmp',filesep];

%% Initialize files
tempFile=tempname(tmp_dir);
imgFile=[tempFile,'.ppm'];
segFile=[tempFile,'_superpix.ppm'];

% Make image file (in ppm format) if it does not exist already
if ~ischar(img)
    % Image matrix given, make ppm file
    imwrite(img,imgFile,'ppm');
    tempImage=1;
else
    % Image file name given (image has to be in ppm format, check)
    if ~strcmp(img(end-2:end),'ppm')
        error('If image file is given, it must be in ppm format');
    end
    imgFile=img;
    tempImage=0;
end

%% Run segmentation algorithm by Pedro F. Felzenszwalb and Daniel P. Huttenlocher
if isunix
    cmd = ['sh -c "unset LD_LIBRARY_PATH;' soft_dir 'segment ' num2str(sigma) ' ' num2str(k) ' ' num2str(min_area) ' "' imgFile '" "' segFile '""' ];
else
    cmd = [soft_dir 'segment ' num2str(sigma) ' ' num2str(k) ' ' num2str(min_area) ' ' imgFile ' ' segFile ];
end
system(cmd);

%% Read in segmentation result
superPixSeg=imread(segFile);

%% Return also index image if needed
if nargout>1
    superPixInd=mvg_numerizeLabels(superPixSeg);
end

%% Return also superpixel boundaries if needed
if nargout>2
    superPixBounds=superpixels2boundaries(superPixInd);
end

%% Delete temporary all created temporary files
if exist(imgFile,'file') && tempImage==1
    delete(imgFile);
end
if exist(segFile,'file')
    delete(segFile);
end

author: Cloud-CV

repo: object-proposals


text:
function [D]=distND(V1,V2)

D=zeros(size(V1,1),size(V2,1));

for q=1:1:size(V1,2) %For all dimensions
    A=V1(:,q); %Coordinates of first set
    B=V2(:,q); %Coordinates of second set

    AB=A(:,ones(1,size(B,1)));
    BA=B(:,ones(1,size(A,1)))';

    D=D+(AB-BA).^2;
end
D=sqrt(D);

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2020 Kevin Mattheus Moerman
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
%TRANIMATE Animate a coordinate frame
%
% TRANIMATE(P1, P2, OPTIONS) animates a 3D coordinate frame moving from pose X1
% to pose X2.  Poses X1 and X2 can be represented by:
%   - homogeneous transformation matrices (4x4)
%   - orthonormal rotation matrices (3x3)
%   - Quaternion
%
% TRANIMATE(X, OPTIONS) animates a coordinate frame moving from the identity pose
% to the pose X represented by any of the types listed above.
%
% TRANIMATE(XSEQ, OPTIONS) animates a trajectory, where XSEQ is any of
%   - homogeneous transformation matrix sequence (4x4xN)
%   - orthonormal rotation matrix sequence (3x3xN)
%   - Quaternion vector (Nx1)
%
% Options::
%  'fps', fps    Number of frames per second to display (default 10)
%  'nsteps', n   The number of steps along the path (default 50)
%  'axis',A      Axis bounds [xmin, xmax, ymin, ymax, zmin, zmax]
%  'movie',M     Save frames as files in the folder M
%
%  Additional options are passed through to TRPLOT.
%
% Notes::
% - Uses the Animate helper class to record the frames.
% - Poses X1 and X2 must both be of the same type
% - The 'movie' options saves frames as files NNNN.png.
% - To convert frames to a movie use a command like:
%        ffmpeg -r 10 -i %04d.png out.avi
%
% See also TRPLOT, Animate.

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

% TODO
%  auto detect the axis scaling
function tranimate(P2, varargin)

    opt.fps = 10;
    opt.nsteps = 50;
    opt.axis = [];
    opt.movie = [];

    [opt, args] = tb_optparse(opt, varargin);

    if ~isempty(opt.movie)
        anim = Animate(opt.movie);
    end
    P1 = [];

    % convert quaternion and rotation matrix to hom transform
    if isa(P2, 'Quaternion')
        T2 = P2.T;   % convert quaternion to transform
        if ~isempty(args) && isa(args{1},'Quaternion')
            T1 = T2;
            Q2 = args{1};
            T2 = Q2.T;
            args = args(2:end);
        else
            T1 = eye(4,4);
        end
    elseif isrot(P2)
        T2 = r2t(P2);
        if ~isempty(args) && isrot(args{1})
            T1 = T2;
            T2 = r2t(args{1});
            args = args(2:end);
        else
            T1 = eye(4,4);
        end
    elseif ishomog(P2)
        T2 = P2;
        if ~isempty(args) && ishomog(args{1})
            T1 = T2;
            T2 = args{1};
            args = args(2:end);
        else
            T1 = eye(4,4);
        end
    end

    % at this point
    %   T1 is the initial pose
    %   T2 is the final pose
    %
    %  T2 may be a sequence

    if size(T2,3) > 1
        % tranimate(Ts)
        % we were passed a homog sequence
        if ~isempty(P1)
            error('only 1 input argument if sequence specified');
        end
        Ttraj = T2;
    else
        % tranimate(P1, P2)
        % create a path between them
        Ttraj = ctraj(T1, T2, opt.nsteps);
    end

    if isempty(opt.axis)
        % create axis limits automatically based on motion of frame origin
        t = transl(Ttraj);
        mn = min(t) - 1.5;  % min value + length of axis + some
        mx = max(t) + 1.5;  % max value + length of axis + some
        axlim = [mn; mx];
        axlim = axlim(:)';
        args = [args 'axis' axlim];
    end

    hg = trplot(eye(4,4), args{:});  % create a frame at the origin

    % animate it for all poses in the sequence
    for i=1:size(Ttraj,3)
        T = Ttraj(:,:,i);
        trplot(hg, T);

        if ~isempty(opt.movie)
            anim.add();
        end

        pause(1/opt.fps);
    end

author: Allopart

repo: rbpf-gmapping


text:
% This file runs univariate linear regression to predict profits of food trucks based on previous

% actual values of profits in $10,000s in various cities with populations in 10,000s respectively.

clear ; close all; clc ;

fprintf('Plotting data\n');

data = load('data_.txt');
x = data(:, 1); Y = data(:, 2);
n = length(Y); % Number of training examples.

plotdata(x, Y);

fprintf('Program paused, press enter to continue\n');

pause;

x = [ones(n, 1), data(:,1)];
Theta = zeros(2, 1);

noi = 1500;   % Number of iterations in gradient descent.
Alpha = 0.01; % Learning rate.

fprintf('Testing the cost function\n')

j = computecost(x, Y, Theta);
fprintf('With Theta = [0 ; 0]\nCost computed = %f\n', j);
fprintf('Expected cost value (approx) 32.07\n');

j = computecost(x, Y, [-1 ; 2]);
fprintf('With theta = [-1 ; 2]\nCost computed = %f\n', j);
fprintf('Expected cost value (approx) 54.24\n');

fprintf('Program paused, press enter to continue\n');

pause;

fprintf('Running gradient descent\n');

Theta = gradientdescent(x, Y, Theta, Alpha, noi);

fprintf('Theta found by gradient descent\n');
fprintf('%f\n', Theta);
fprintf('Expected Theta vector (approx)\n');
fprintf(' -3.6303\n  1.1664\n\n');

hold on;                                                           % To plot hypothesis on data.

plot(x(:, 2), x * Theta, '-');
legend('Training data', 'Linear regression');

predict1 = [1, 3.5] * Theta;
fprintf('For population = 35,000, we predict a profit of %f\n',...
    predict1*10000);

predict2 = [1, 7] * Theta;
fprintf('For population = 70,000, we predict a profit of %f\n',...
    predict2*10000);

author: TheAlgorithms

repo: MATLAB-Octave


text:
function [packet, state]= bbci_control_ERPSpeller(cfy_out, state, event, opt)
%BBCI_CONTROL_ERP_SPELLER - Generate control signal for ERP-based Hex-o-Spell
%
%Synopsis:
%  PACKET= bbci_control_ERPSpeller(CFY_OUT, EVENT, SETTINGS, STATE)
%
%Arguments:
%  CFY_OUT - Output of the classifier
%  STATE - Internal state variable
%  EVENT - Structure that specifies the event (fields 'time' and 'desc')
%      that triggered the evaluation of this control.
%  SETTINGS - Structure specifying the following feedback paramters:
%    .nClasses - number of classes of stimuli from which to select
%    .nSequences - number of sequences (repetitions) that should be used
%        in order to determine one choice
%
%Output:
% PACKET: Variable/value list in a CELL defining the control signal that
%     is to be sent via UDP to the Speller application.
% STATE: Updated internal state variable

% 02-2011 Benjamin Blankertz

if isempty(state),
  state.counter= zeros(1, opt.nClasses);
  state.output= zeros(1, opt.nClasses);
end

this_cue= opt.mrk2feedback_fcn(event.desc);
state.counter(this_cue)= state.counter(this_cue) + 1;
state.output(this_cue)= state.output(this_cue) + cfy_out;

if sum(state.counter) >= opt.nClasses*opt.nSequences,
  idx= find(state.counter>0);  % avoid divide by zero
  state.output(idx)= state.output(idx) ./ state.counter(idx);
  [max_score, selected_class]= min(state.output);
  packet= {'i:cl_output', selected_class};
  state.counter(:)= 0;
  state.output(:)= 0;
else
  packet= [];
end

author: bbci


text:
function [cfg] = ft_layoutplot(cfg, data)

% FT_LAYOUTPLOT makes a figure with the 2-D layout of the channel positions
% for topoplotting and the individual channel axes (i.e. width and height
% of the subfigures) for multiplotting. A correct 2-D layout is a
% prerequisite  for plotting the topographical distribution of the
% potential or field distribution, or for plotting timecourses in a
% topographical arrangement.
%
% This function uses the same configuration options as prepare_layout and
% as the topoplotting and multiplotting functions. The difference is that
% this function plots the layout without any data, which facilitates
% the validation of your 2-D layout.
%
% Use as
%   ft_layoutplot(cfg, data)
%
% There are several ways in which a 2-D layout can be made: it can be read
% directly from a *.lay file, it can be created based on 3-D electrode or
% gradiometer positions in the configuration or in the data, or it can be
% created based on the specification of an electrode of gradiometer file.
%
% You can specify either one of the following configuration options
%   cfg.layout      = filename containg the layout
%   cfg.rotate      = number, rotation around the z-axis in degrees (default = [], which means automatic)
%   cfg.projection  = string, 2D projection method can be 'stereographic', 'ortographic', 'polar', 'gnomic' or 'inverse' (default = 'orthographic')
%   cfg.elec        = structure with electrode positions or filename, see FT_READ_SENS
%   cfg.grad        = structure with gradiometer definition or filename, see FT_READ_SENS
%   cfg.opto        = structure with optode definition or filename, see FT_READ_SENS
%   cfg.output      = filename to which the layout will be written (default = [])
%   cfg.montage     = 'no' or a montage structure (default = 'no')
%   cfg.image       = filename, use an image to construct a layout (e.g. usefull for ECoG grids)
%   cfg.box         = string, 'yes' or 'no' whether box should be plotted around electrode (default = 'yes')
%   cfg.mask        = string, 'yes' or 'no' whether the mask should be plotted (default = 'yes')
%   cfg.visible     = string, 'on' or 'off' whether figure will be visible (default = 'on')
%   cfg.figure      = 'yes' or 'no', whether to open a new figure. You can also specify a figure handle from FIGURE, GCF or SUBPLOT. (default = 'yes')
%   cfg.position    = location and size of the figure, specified as [left bottom width height] (default is automatic)
%   cfg.renderer    = string, 'opengl', 'zbuffer', 'painters', see RENDERERINFO (default is automatic, try 'painters' when it crashes)
%
% Alternatively the layout can be constructed from either
%   data.elec     structure with electrode positions
%   data.grad     structure with gradiometer definition
%
% Alternatively, you can specify
%   cfg.layout = 'ordered'
% which will give you a 2-D ordered layout. Note that this is only suited
% for multiplotting and not for topoplotting.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
% If you specify this option the input data will be read from a *.mat
% file on disk. This mat files should contain only a single variable named 'data',
% corresponding to the input structure.
%
% See also FT_PREPARE_LAYOUT, FT_TOPOPLOTER, FT_TOPOPLOTTFR, FT_MULTIPLOTER, FT_MULTIPLOTTFR

% Undocumented options
%   cfg.montage

% Copyright (C) 2006-2008, 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$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar data
ft_preamble provenance data

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% the data can be passed as input argument or can be read from disk
hasdata = exist('data', 'var');

if hasdata
  % check if the input data is valid for this function
  data = ft_checkdata(data);
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'renamed', {'elecfile', 'elec'});
cfg = ft_checkconfig(cfg, 'renamed', {'gradfile', 'grad'});
cfg = ft_checkconfig(cfg, 'renamed', {'optofile', 'opto'});
cfg = ft_checkconfig(cfg, 'renamed', {'newfigure', 'figure'});

% set the defaults
cfg.box       = ft_getopt(cfg, 'box', 'yes');
cfg.mask      = ft_getopt(cfg, 'mask', 'yes');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% extract or generate the layout information
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tmpcfg = keepfields(cfg, {'layout', 'channel', 'rows', 'columns', 'commentpos', 'skipcomnt', 'scalepos', 'skipscale', 'projection', 'viewpoint', 'rotate', 'width', 'height', 'elec', 'grad', 'opto', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
if hasdata
  layout = ft_prepare_layout(tmpcfg, data);
else
  layout = ft_prepare_layout(tmpcfg);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot all details pertaining to the layout in one figure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% open a new figure with the specified settings
h = open_figure(keepfields(cfg, {'figure', 'position', 'visible', 'renderer'}));

if isfield(cfg, 'image') && ~isempty(cfg.image)
  % start with the background image
  fprintf('reading background image from %s\n', cfg.image);
  img = imread(cfg.image);
  img = flipdim(img, 1); % in combination with "axis xy"

  bw = 1;

  if bw
    % convert to greyscale image
    img = mean(img, 3);
    imagesc(img);
    colormap gray
  else
    % plot as RGB image
    image(img);
  end
  axis equal
  axis off
  axis xy
end

ft_plot_layout(layout, 'point', true, 'box', cfg.box, 'label', true, 'mask', cfg.mask, 'outline', true);

% the following code can be used to verify a bipolar montage, given the
% layout of the monopolar channels
if isfield(cfg, 'montage') && ~isempty(cfg.montage)
  fprintf('plotting an arrow for each of the bipolar electrode pairs\n');
  % the arrow begins at the +1 electrode
  % the arrow ends   at the -1 electrode
  for i=1:length(cfg.montage.labelnew)
    begindx = find(cfg.montage.tra(i,:)==+1);
    endindx = find(cfg.montage.tra(i,:)==-1);
    if ~numel(begindx)==1 || ~numel(endindx)==1
      % the re-referenced channel does not seem to be a bipolar pair
      continue
    end
    % find the position of the begin and end of the arrow
    beglab = cfg.montage.labelold{begindx};
    endlab = cfg.montage.labelold{endindx};
    begindx = find(strcmp(layout.label, beglab)); % the index in the layout
    endindx = find(strcmp(layout.label, endlab)); % the index in the layout
    if ~numel(begindx)==1 || ~numel(endindx)==1
      % one of the channels in the bipolar pair does not seem to be in the layout
      continue
    end

    begpos = layout.pos(begindx,:);
    endpos = layout.pos(endindx,:);
    arrow(begpos, endpos, 'Length', 5)

  end % for all re-referenced channels
end % if montage

% this is needed for the figure title
if isfield(cfg, 'dataname') && ~isempty(cfg.dataname)
  dataname = cfg.dataname;
elseif isfield(cfg, 'inputfile') && ~isempty(cfg.inputfile)
  dataname = cfg.inputfile;
elseif nargin>1
  dataname = arrayfun(@inputname, 2:nargin, 'UniformOutput', false);
else
  dataname = {};
end

% set the figure window title
if ~isempty(dataname)
  set(h, 'Name', sprintf('%d: %s: %s', double(h), mfilename, join_str(', ', dataname)));
else
  set(h, 'Name', sprintf('%d: %s', double(h), mfilename));
end
set(h, 'NumberTitle', 'off');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% deal with the output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous data
ft_postamble provenance
ft_postamble savefig

% add a menu to the figure, but only if the current figure does not have subplots
menu_fieldtrip(gcf, cfg, false);

if ~ft_nargout
  % don't return anything
  clear cfg
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function getMetaboliteMsDistr(mets, molfiledir, msdistrdir, phs)
% Calculate microspecies distributions at all pH values in phs using
% ChemAxon's cxcalc. Metabolite mol files in molfiledir are input to cxcalc. Mol file names
% should correspond to the metabolite ID in mets (it is assumed that
% compartment assignments are appended to the end of metabolite ID in the
% format `metID[c]`). Microspecies distributions are returned as .sdf files
% in `msdistrdir`
%
% USAGE:
%
%    getMetaboliteMsDistr(mets, molfiledir, msdistrdir, phs)
%
% INPUTS:
%    mets:          metabolites
%    molfiledir:    directory with mol files
%    msdistrdir:    directory with microspecies distributions
%    phs:           pH values

d = dir(molfiledir);
molfilelist = cat(2,d.name);
molfilelist = regexp(molfilelist,'.mol','split');
molfilelist = regexprep(molfilelist,'\.','');
molfilelist = molfilelist(~cellfun('isempty',molfilelist))';

nMet = length(mets);
done = {};

model.mets = mets;

for n = 1:nMet
    disp(n)

    if ~any(strcmp(model.mets{n}(1:(end-3)),done))
        done = [done; {model.mets{n}(1:(end-3))}];

        if any(strcmp(molfilelist,mets{n}(1:(end-3))))

            for m = 1:length(phs)
                ph = phs(m);

                if ispc
                    system(['cxcalc -o ', msdistrdir, filesep, model.mets{n}(1:(end-3)), '_msdistr_ph', num2str(ph), '.sdf msdistr -H ', num2str(ph), ' -M true ', molfiledir, filesep, model.mets{n}(1:(end-3)), '.mol']);
                end

            end

        end

    end

end

author: opencobra

repo: cobratoolbox


text:
function [fMeanMag, fBValue, fStdDev, fAValue] = ...
        calc_bmemag(mag, fBinning)
    %CALC_BMEMAG Calculate Maximum Likelihood b-value
    % [fMeanMag, fBValue, fStdDev, fAValue] = ...
    %                   calc_bmemag(mag, fBinning)
    % Calculates the mean magnitute, the b-value based
    % on the maximum likelihood estimation, the a-value and the
    % standard deviation of the b-value
    %
    % Input parameters:
    %   mag             Vector of magnitudes
    %   fBinning        Bin size in magnitude (default 0.1)
    %
    % Output parameters:
    %   fMeanMag        Mean magnitude
    %   fBValue         b-value
    %   fStdDev         Standard deviation of b-value
    %   fAValue         a-value

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

    % Set the default value if not passed to the function
    if ~exist('fBinning')
        fBinning = 0.1;
    end

    % Calculate the minimum and mean magnitude, length of catalog
    nLen = length(mag);
    fMinMag = min(mag);
    fMeanMag = mean(mag);
    % Calculate the b-value (maximum likelihood)
    fBValue = (1/(fMeanMag-(fMinMag-(fBinning/2))))*log10(exp(1));
    % Calculate the standard deviation
    fStdDev = (sum((mag-fMeanMag).^2))/(nLen*(nLen-1));
    fStdDev = 2.30 * sqrt(fStdDev) * fBValue^2;
    % Calculate the a-value
    fAValue = log10(nLen) + fBValue * fMinMag;
end

author: geoscience-community-codes

repo: GISMO


text:
function successFlag=unfoldMeshFromGUI(meshFileName,grayFileName,flatFileName,startCoords,scaleFactor,perimDist,statusHandle,busyHandle,showFigures,saveExtra,truePerimDist);
% function successFlag=unfoldMeshFromGUI(meshFileName,grayFileName,flatFileName,startCoords,scaleFactor,perimDist,statusHandle,busyHandle,showFigures,saveExtra,truePerimDist);
% To be called from unfoldMeshAlpha.m
% This is pre-release code 0.1
%
% ARW 021401
% Unfolds a legally-triangulated mesh using floater/tuele's sparse linear matrix method. Then maps L1 grey nodes to the flattened
% mesh and squeezes remaining grey layers to L1 to generate a flat.mat file for use in mrLoadRet

NPERIMS=1;
SAVE_INTERMEDIATE=0;

if (length(startCoords)~=3)
	error ('Error: you must enter 3 start coords');
end

startTime=now;

% READ IN THE ORIGNIAL MRGRAY MESH
statusStringAdd(statusHandle,['Loading: ',meshFileName]);

mesh=mrReadMrM(meshFileName,0); %
statusStringAdd(statusHandle,'Finding faces');

mesh.faceIndexList=findFaces(mesh,busyHandle);
% The entries in mesh.vertices aren't unique. That is, the same
% vertex can be represented by two different numbers.  This happens
% because strips always contain shared vertices.
%
% This next series of functions finds the unique vertices and
% and builds hash tables for vert->uniqueVerts and uniqueVerts->vert.
% Using these tables, we can adjust the list of colors and faces
% to represent each vertex by a unique index.

[mesh.uniqueVertices,mesh.vertsToUnique,mesh.UniqueToVerts]= unique(mesh.vertices,'rows'); % Get rid of identical points
mesh.uniqueVertices=mesh.uniqueVertices./(repmat(scaleFactor,length(mesh.uniqueVertices),1));
mesh.uniqueNormal=mesh.normal(mesh.vertsToUnique);

mesh.uniqueCols=mesh.rgba(mesh.vertsToUnique,:); % try to assign the colors that came in from mrGray

mesh.uniqueFaceIndexList=findUniqueFaceIndexList(mesh); % this gets rid of faces with multiple duplicate vertices and different permutations of the same face

% Now we find the connection matrix:
% a sparse matrix of nxn points where M(i,j) is 1 if i and j are connected

statusStringAdd(statusHandle,'Finding connection matrix.');

[mesh.connectionMatrix]=findConnectionMatrix(mesh);
messageString=sprintf('%d connections found',length(mesh.connectionMatrix));
statusStringAdd(statusHandle,messageString);

statusStringAdd(statusHandle,'Checking group perimeter.');

% Check to make sure that this is a clean surface: no edge points yet.
edgeList=findGroupPerimeter(mesh,1:length(mesh.uniqueVertices));
if (~isempty(edgeList))
	error('Error - initial mesh is not closed!');
	break;
else
	fprintf('\nInitial mesh is closed');
end

statusStringAdd(statusHandle,'Finding closest mesh point.');

% **********************************
% THIS IS ****** CHECKPOINT 1 ******
% **********************************

[startNode,snDist]=assignToNearest(mesh.uniqueVertices,startCoords);

messageString=sprintf('Start node %d selected at %d voxel units from input coords.',startNode,sqrt(snDist));
statusStringAdd(statusHandle,messageString);

if (sqrt(snDist)>15)
	beep;
	messageString=sprintf('** Warning: mesh node far from start coord. Expect trouble.');
	statusStringAdd(statusHandle,messageString)
end

% Have replaced mrManDist with dll dijkstra
statusStringAdd(statusHandle,'Finding distances from start node');

D=find3DNeighbourDists(mesh,scaleFactor);

% find distances to all the nodes
% When we go to the newer mrGray, dimdist can
% be recovered from the mesh.parameters field, we think.
% Probably mesh.parameters.voxsize.  See mrReadMrM.m

mesh.dist=dijkstra(D,startNode);

% Now we want to generate a perimeter by thresholding these distances
% First find a good distance to take for the perimeter

mesh.perimDist=perimDist;
messageString=sprintf('Perimeter minimum distance=%d',perimDist);
statusStringAdd(statusHandle,messageString);

statusStringAdd(statusHandle,'Using threshold to find perimeter(s).');

% Find perims with simple thold
insideNodes=find(mesh.dist<=perimDist);
insideNodes=insideNodes(:);

% Eliminate 'hanging nodes' : nodes that are not part of a face:
%insideNodes=removeHangingNodes(mesh,insideNodes);

statusStringAdd(statusHandle,'Internal nodes found. Defining perimeter.');

% Find the perimeter points of this group (there can be more than perimeter at this point)
numBadNodes=9999999;
 [perimeterEdges,eulerCondition]=findGroupPerimeter(mesh,insideNodes);

while (numBadNodes>0)
    [perimeterEdges,eulerCondition]=findGroupPerimeter(mesh,insideNodes);
	length(perimeterEdges)
	length(unique(perimeterEdges,'rows'))
    fprintf('\nEuler number=%d',eulerCondition);

    badPerimNodes=findBadPerimNodes(mesh,perimeterEdges);
    numBadNodes=length(badPerimNodes)
    disp(badPerimNodes);

    if(numBadNodes)
        [insideNodes]=setdiff(insideNodes,badPerimNodes);
		%insideNodes=removeHangingNodes(mesh,insideNodes);
    end
end

% The euler number will tell you whether you've got a mesh with no holes in it.

messageString=sprintf('Euler number for this set=%d',eulerCondition);
statusStringAdd(statusHandle,messageString);
uniquePerimPoints=unique(perimeterEdges(:));
messageString=sprintf('%d unique perimeter points found.',length(uniquePerimPoints));
statusStringAdd(statusHandle,messageString);
[orderedUniquePerimeterPoints,biggest]=orderMeshPerimeterPointsAll(mesh,perimeterEdges);
nPerims=size(orderedUniquePerimeterPoints);
fprintf('\n%d different perimeters found',nPerims);
orderedUniquePerimeterPoints=orderedUniquePerimeterPoints{biggest}.points;

%insideNodes=[insideNodes(:),setdiff(perimeterEdges(:),perimeterEdges2(:))];
%perimeterEdges=perimeterEdges2;

% internal points (the ones we want)
unfoldMesh.connectionMatrix=mesh.connectionMatrix(insideNodes,insideNodes);
unfoldMesh.uniqueVertices=mesh.uniqueVertices(insideNodes,:);
unfoldMesh.dist=mesh.dist(insideNodes);

% Convert the edges to feed into orderMeshPerimeterPoints
fullEdgePointList=perimeterEdges(:);

[numEdges,x]=size(perimeterEdges);

newEdges=zeros((numEdges*2),1);
statusStringAdd(statusHandle,'Finding sub-mesh edges.');

for t=1:(numEdges*2)
	if ((~mod(t,100)) & busyHandle)
		updateBusybar(busyHandle,t);
	end
	newEdges(t)=find(insideNodes==fullEdgePointList(t));
end

newEdges=reshape(newEdges,numEdges,2);
statusStringAdd(statusHandle,'Finding sub-mesh perim.');

% Find the perimeter points.
unfoldMesh.orderedUniquePerimeterPoints=zeros(length(orderedUniquePerimeterPoints),1);

for t=1:length(orderedUniquePerimeterPoints)
	f1=find(insideNodes==orderedUniquePerimeterPoints(t));%
	unfoldMesh.orderedUniquePerimeterPoints(t)=f1;%orderMeshPerimeterPoints(newEdges);

end

% Unfolding bit...
% Now we'd like to unfold this.
% Need the following things...
% Connection matrix N (nxn) - almost the same as the connection matrix except that rows=rows/sum(rows)
% and all points on the perimeter have been removed.
% X0 - a px2 matrix containing the 2D locations of the perimeter points.
% P - The perimeter connection matrix: (nxp) 'whose (i,j)th entry is 1/mi when perimeter node j is connected
% to sample node i.

% Find the N and P connection matrices
statusStringAdd(statusHandle,'Finding sub-mesh con. mat.');
[N,P,unfoldMesh.internalNodes]=findNPConnection(unfoldMesh);
[unfoldMesh.distSQ]=find3DNeighbourDists(unfoldMesh);   % Here we find the 3D distance from each point to its neighbours.

% Assign the initial perimeter points - they're going to go in a circle for now...
numPerimPoints=length(unfoldMesh.orderedUniquePerimeterPoints);

perimeterDists=mesh.dist(orderedUniquePerimeterPoints); % Distance of each perimeter point from the start node

statusStringAdd(statusHandle,'Assigning perimeter points');

% We'd like to place the perimeter points down in an intelligent manner. We can place them at the correct distance from the start node
% and at the correct distance from each other. I think.
% We already know their distance from the start node, now we'd like to get their distances
% from each other.

% Should be able to extract this from unfoldMesh.distSQ

%  start at the first perimeter point in unfoldMesh.orderedUniquePerimeterPoints
% Find the distance between this and the next point, etc etc...
nPerimPoints=length(unfoldMesh.orderedUniquePerimeterPoints);
interPerimDists=zeros(nPerimPoints,1);

for thisPerimPoint=1:nPerimPoints
	nextIndex=mod((thisPerimPoint+1),nPerimPoints);
	if(nextIndex==0)
		nextIndex=1;
	end

	interPerimDists(thisPerimPoint)=unfoldMesh.distSQ(unfoldMesh.orderedUniquePerimeterPoints(thisPerimPoint),unfoldMesh.orderedUniquePerimeterPoints(nextIndex));
end

interPerimDists=sqrt(interPerimDists);

unfoldMesh.X_zero=assignPerimeterPositions(perimeterDists); % Can set distances around circle to match actual distances from the center.
% Angular positions will come next.

statusStringAdd(statusHandle,'Solving position equation (slow)');
X=(speye(size(N)) - N) \ (sparse(P * unfoldMesh.X_zero));

unfoldMesh.N=N;
unfoldMesh.P=P;
unfoldMesh.X=X;

% Find out the differences between
dist2DSQ=find2DNeighbourDists(unfoldMesh);
d1=sparse(unfoldMesh.distSQ)-sparse(dist2DSQ);
d1=abs(d1);

goodness=sum(sum(d1.^2));

messageString=sprintf('Current error per node: %d',full(sqrt(goodness))/length(insideNodes));
statusStringAdd(statusHandle,messageString);

% Show the mesh
if (showFigures)

	statusStringAdd(statusHandle,'Displaying unfold');
	figure(50);

	hold off;
	gplot(unfoldMesh.N,unfoldMesh.X);

	axis equal;
	axis off;
	zoom on

end

% Finally - the mapping of grey to mesh points takes place using the entire mesh.
% Therefore, we need to generate X for the mesh as well as the unfold mesh'

unfoldToOrigPerimeter=insideNodes(unfoldMesh.orderedUniquePerimeterPoints);
unfoldToOrigInside=insideNodes(unfoldMesh.internalNodes);

mesh.X=zeros(length(mesh.uniqueVertices),2);
mesh.X(unfoldToOrigPerimeter,:)=unfoldMesh.X_zero;
mesh.X(unfoldToOrigInside,:)=unfoldMesh.X;

hasCoords=[unfoldToOrigPerimeter(:);unfoldToOrigInside(:)];
coords=mesh.X(hasCoords,:);%mesh.X(unfoldToOrigInside,:);
dists=mesh.dist(hasCoords);

% use griddata to image the distance map

	mesh.distMap=makeDistanceImage(coords,dists,128);
	ZI=mesh.distMap;

	if (showFigures)
		figure(51);

		imagesc(mesh.distMap);
		axis image;

		colormap hot;
		title('Manifold distance map');
		colorbar;
	end

mesh.strain=zeros(length(mesh.uniqueVertices),1);
strain=d1.^2;
mesh.strain(unfoldToOrigPerimeter)=sum(strain(unfoldMesh.orderedUniquePerimeterPoints,:)');
mesh.strain(unfoldToOrigInside)=sum(strain(unfoldMesh.internalNodes,:)');
% Record which nodes in the big mesh are in the unfold
mesh.insideNodes=insideNodes;
if (SAVE_INTERMEDIATE)
	statusStringAdd(statusHandle,'Saving intermediate data.');
	save ('meshOutTemp.mat','mesh');
end

% (At this point, switch to testGray....)

% *********************************************************************************
% ************************************************************************************
% **********************************************************************************

statusStringAdd(statusHandle,'Reading grey graph...');

[gNodes, gEdges, gvSize] = readGrayGraph_progress(grayFileName,0);

% Get the indices for all the gnodes (all 3 layers)
l1NodeIndices=find(gNodes(6,:)==1);
l2NodeIndices=find(gNodes(6,:)==2);
l3NodeIndices=find(gNodes(6,:)==3);
l4NodeIndices=find(gNodes(6,:)==4);

% Extract the layer 1 nodes
l1gNodes=gNodes(:,l1NodeIndices);
l1mesh.vertices=l1gNodes(1:3,:);
l1mesh.indices=l1NodeIndices;

% How many gNodes are there?
nGnodes=length(gNodes);
% How many gEdges are there?
nGedges=length(gEdges);

% We want to make a grey node connection matrix - which grey nodes are connected to which other gnodes?
sp=sparse(nGnodes,nGnodes);
i=zeros(nGnodes*30,1); % no more that 30 conenctions per gNode on average!
j=i;

offset=1;

statusStringAdd(statusHandle,'Finding grey connection matrix (slow)');

for t=1:nGnodes % for each gNode...
	if ((~mod(t,1000)) & busyHandle)
		updateBusybar(busyHandle,t);
	end
	% Find its edges (the nodes of the things that it's connected to...)
	thisOffset=gNodes(5,t);
	thisNumEdges=gNodes(4,t);
	theseEdges=gEdges(thisOffset:(thisOffset-1+thisNumEdges));

	% add these to i,j - eventually we'll call sp=sparse(i,j,s,nGnodes,nGnodes)
	% i contains the y coords, j contains the x coords
	endPoint=offset+thisNumEdges-1;

	i(offset:endPoint)=ones(1,thisNumEdges)*t;
	j(offset:endPoint)=theseEdges;

	offset=endPoint+1;

	% % This takes about 18 secs on gwyrdd (nGnodes=32000)
	% and 2.1 seconds on Gwyn :)
end

i=i(1:offset-1);
j=j(1:offset-1);
s=ones(size(i));
sp=sparse(i,j,s,nGnodes,nGnodes);

clear i;
clear j;
clear s;

% We can assign layer 1 grey nodes to the white matter mesh using assignToNearest.dll (see assignToNearest.c)
% Can't do this for higher levels of grey matter 'cos they might get mis-assigned. (Also, potential problem near
% very crinkly edges. - Could we accidentally assign a l1 grey matter node to the wrong WM point?)

% So for higher grey matter points, we have to restrict the possible sub node search space by assigning them >only< to
% points they are connected to. Note that a single layer2 grey node may be connected to several l1 nodes

% The gray is defined over the entire mesh but we only want to deal with gray points over the
% unfolded part. The strategy should be....
% 1) do assignToNearest for each mesh point to find the nearest connected l1 node
% 2) Use the set of l1 nodes found in 1) to build up a list of other connected gray nodes
% 3) Proceed as before...

statusStringAdd(statusHandle,'Mapping L1 to mesh.');
% Find 3D coords of all the l1 gnodes
l1GNodeCoords=l1gNodes(1:3,:)';

% Find 3D coords of all the mesh points (not just the unfolded ones) We have to do this in order
% to shift the two sets to a common mean
meshCoords=mesh.uniqueVertices;

% These coordinate sets seem to have very similar shapes (as we'd expect) but offset to each other. Remove the means to
% center them (but there must be a better way...);
% Mesh coords are in mm, gnode coords are in voxels:
%scaleFactor=[240/256 240/256 1.2]; % This has to be sent in eventually.... IMPORTANT!

l1GNodeCoordsN=l1GNodeCoords;

% There are roughly half as many l1gnodes as there are unique vertices in the mesh.

% Now restrict the mesh coords to just those points in the unfold
meshCoordsN=meshCoords(mesh.insideNodes,:);

% And now restrict the set of l1 gray nodes so that only those that are relatively near the
% mesh are included in the search - this is done first as a simple bounds check
boundedL1NodeIndices=boundsCheck3D(min(meshCoordsN)-3,max(meshCoordsN)+3,l1GNodeCoordsN);
boundedL1GNodes=l1GNodeCoordsN(boundedL1NodeIndices,:);
boundedL1NodeIndices=l1NodeIndices(boundedL1NodeIndices); % This is now a list of indices into the full gNode array

statusStringAdd(statusHandle,'Finding nearest Ll gray points to mesh (very slow)');

% Now we >could< restrict further at this point by running
% [meshToBoundedL1Indices,sqrDist]=assignToNearest(boundedL1GNodes,meshCoordsN);
% to define a subset of l1gNodes closest to the mesh. However, this seems to bugger things up - probably becasue the
% l1gNodes are sampled more densely than the mesh in places.

% What we do instead is run
[boundedL1toMeshIndices,sqrDist]=assignToNearest(meshCoordsN,boundedL1GNodes);
% This returns a list of indices into the meshCoordsN array that links a single 3D mesh point to each l1Gnode)

% and then eliminate any l1gNodes that are more than a set distance away from the mesh - here 3.2mm

% *************************
closeEnough=find(sqrDist<2);
% *************************

boundedL1toMeshNodes=boundedL1GNodes(closeEnough,:); % remember, assignToNearest returns the squared distance
boundedL1toMeshIndices=boundedL1toMeshIndices(closeEnough);

% For each member of the bounded l1Gnodes, this tells us the index of the full mesh point that it maps to.
fullBoundedL1toMeshIndices=insideNodes(boundedL1toMeshIndices);

restL1NodeIndices=boundedL1NodeIndices(closeEnough);
statusStringAdd(statusHandle,'Setting L1 glocs');
% We can start setting gLocs
l1Glocs3d=boundedL1toMeshNodes;
l1Glocs2d=mesh.X(fullBoundedL1toMeshIndices,:);

%l1Curvature=mesh.uniqueNormal(fullBoundedL1toMeshIndices);

l1ConMat=sp(restL1NodeIndices,restL1NodeIndices);

% Prune the l1gray nodes here ? Eliminate connections that are too big (roughly >2 or so.).
% PRUNE
% PRUNE

if (showFigures);
	statusStringAdd(statusHandle,'Displaying L1 gray mesh.');
	figure(8);
	hold off;
	gplot(l1ConMat,l1Glocs2d);
	title('L1 gray mesh');
	zoom on
end

% Now we (only!) have to find l2tol1Indices, l3tol2Indices and l4tol3Indices. This is faster since for each point, we restrict its
% potential nearest neighbours to points that it is connected to in the previous layer.
% We also restrict the l2 nodes to just those that are connected to the restricted l1 nodes and the l3 nodes to those connected to the
% l2 nodes.

% Use the full connection matrix to find which l2 Gnodes are connected to the restricted l1Gnodes
statusStringAdd(statusHandle,'Mapping higher levels (2,3,4)');
L1L2sp=sp(:, restL1NodeIndices);

% We assume here that, almost by definition, the only things that l1 nodes connect to are l2 nodes and other l1 nodes
[restL2NodeIndices,dummy]=find(L1L2sp);
restL2NodeIndices=unique(restL2NodeIndices);
restL2NodeIndices=intersect(restL2NodeIndices,l2NodeIndices); % Only take potential l2 node indices

% repeat this for the l3 nodes
L2L3sp=sp(:,restL2NodeIndices);
[restL3NodeIndices,dummy]=find(L2L3sp);
restL3NodeIndices=unique(restL3NodeIndices);
restL3NodeIndices=intersect(restL3NodeIndices,l3NodeIndices);

% repeat again for the l4 nodes
L3L4sp=sp(:,restL3NodeIndices);
[restL4NodeIndices,dummy]=find(L3L4sp);
restL4NodeIndices=unique(restL4NodeIndices);
restL4NodeIndices=intersect(restL4NodeIndices,l4NodeIndices);

% For each l2 node, find the l1 nodes it's connected to, then find the
% 3D coords of the l2 node and all the connected l1 nodes and send them
% in to assignToNearest.

l2ToL1Indices=findNearestConnected(gNodes', restL2NodeIndices,restL1NodeIndices,sp);
l3ToL2Indices=findNearestConnected(gNodes', restL3NodeIndices,restL2NodeIndices,sp);
l4ToL3Indices=findNearestConnected(gNodes', restL4NodeIndices,restL3NodeIndices,sp);

statusStringAdd(statusHandle,'Setting upper levels glocs');
% Set the l2, l3 and l4 glocs3d
l2Glocs3d=gNodes(1:3,restL2NodeIndices)';
l3Glocs3d=gNodes(1:3,restL3NodeIndices)';
l4Glocs3d=gNodes(1:3,restL4NodeIndices)';

% And the gLocs2
l2Glocs2d=l1Glocs2d(l2ToL1Indices,:);
l3Glocs2d=l2Glocs2d(l3ToL2Indices,:);
l4Glocs2d=l3Glocs2d(l4ToL3Indices,:);

% and the curvature
l2Curvature=l1Curvature(l2ToL1Indices);
l3Curvature=l2Curvature(l3ToL2Indices);
l4Curvature=l3Curvature(l4ToL3Indices);

if (showFigures)
	statusStringAdd(statusHandle,'Displaying 2D glocs');
	figure(12);
	hold off;
	subplot(4,1,1);
	plot(l1Glocs2d(:,1),l1Glocs2d(:,2),'.');
	title('L1 2d glocs');
	subplot(4,1,2);
	plot(l2Glocs2d(:,1),l2Glocs2d(:,2),'g.');
	title('L2 2d glocs');
	subplot(4,1,3);
	plot(l3Glocs2d(:,1),l3Glocs2d(:,2),'r.');
	title('L3 2d glocs');
	subplot(4,1,4);
    plot(l4Glocs2d(:,1),l4Glocs2d(:,2),'k.');
	title('L4 2d glocs');

end
% old skool flat.mat structure looks like
%   curvature      47263x1         378104  double array
%   gLocs2d        47263x2         756208  double array
%   gLocs3d        47263x3        1134312  double array
%   gLocs3dfloat       0x0              0  double array
%   startPoint         1x1              8  double array
%   unfList            1x47263     378104  double array
%   xSampGray          1x701         5608  double array

% But all we really need are.....
statusStringAdd(statusHandle,'Creating flat.mat structure');
gLocs2d=[l1Glocs2d;l2Glocs2d;l3Glocs2d;l4Glocs2d];
gLocs3d=[l1Glocs3d;l2Glocs3d;l3Glocs3d;l4Glocs3d];
curvature=[l1Curvature;l2Curvature;l3Curvature;l4Curvature];

% Curvature goes from 1 to 64
curvature=normalize(curvature)*63+1;
endTime=now;

% TODO , Get this file name from a GUI
statusStringAdd(statusHandle,'Saving:');
statusStringAdd(statusHandle,flatFileName);

messageString=sprintf('Unfold started at %s\nFinished at %s',datestr(startTime),dateStr(endTime));
statusStringAdd(statusHandle,messageString);
if (saveExtra)
	statusString=char(get(statusHandle,'UserData'));

	info.perimDist=perimDist;

	infoStr.startTime=datestr(startTime);
	infoStr.endTime=datestr(endTime);

	infoStr.perimType=truePerimDist;
	infoStr.meshFile=meshFileName;
	infoStr.grayFile=grayFileName;

	save (flatFileName,'gLocs2d','gLocs3d','curvature','statusString','infoStr','ZI');
else
	save (flatFileName,'gLocs2d','gLocs3d','curvature');
end
statusStringAdd(statusHandle,'Done.');

successFlag=1;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef factorization_sparse_qrt < factorization_generic
%FACTORIZATION_SPARSE_QRT economy sparse QR of A': (p*A)*(p*A)'=R'*R

% Copyright 2009, Timothy A. Davis, University of Florida

    methods

        function F = factorization_sparse_qrt (A)
            %FACTORIZATION_SPARSE_QRT economy sparse QR: (p*A)*(p*A)'=R'*R
            m = size (A,1) ;
            C = A' ;
            p = sparse (1:m, colamd (C), 1) ;
            R = qr (C*p', 0) ;
            assert (nnz (diag (R)) == m, 'Matrix is rank deficient.') ;
            F.A = A ;
            F.R = R ;
            F.p = p ;
        end

        function disp (F)
            %DISP displays a Q-less economy sparse QR
            fprintf ('  Q-less economy sparse QR factorization of A'': ') ;
            fprintf ('(p*A)*(p*A)'' = R''*R\n') ;
            fprintf ('  A:\n') ; disp (F.A) ;
            fprintf ('  R:\n') ; disp (F.R) ;
            fprintf ('  p:\n') ; disp (F.p) ;
            fprintf ('  is_inverse: %d\n', F.is_inverse) ;
        end

        function x = mldivide_subclass (F,b)
            %MLDIVIDE_SUBCLASS x = A\b using economy sparse QR of A'
            % minimum 2-norm solution of an underdetermined system
            A = F.A ;
            R = F.R ;
            p = F.p ;
            x = A' * (p' * (R \ (R' \ (p * b)))) ;
            e = A' * (p' * (R \ (R' \ (p * (b - A * x))))) ;
            x = x + e ;
        end

        function x = mrdivide_subclass (b,F)
            %MRDIVIDE_SUBCLASS x = b/A using economy sparse QR of A'
            % least-squares solution of an overdetermined problem
            bT = b' ;
            A = F.A ;
            R = F.R ;
            p = F.p ;
            x = p' * (R \ (R' \ (p * (A * bT)))) ;
            e = p' * (R \ (R' \ (p * (A * (bT - A' * x))))) ;
            x = (x + e)' ;
        end
    end
end

repo: REST


text:
function result = saveImageFile(arg1, arg2, arg3);
% saveImageFile(IMGDIR, fname, res);
% saveImageFile(IMGFULLFILEPATH, res);
% res=200 for spectrograms
global paths PARAMS; % we need to knwo the value of PARAMS.mode

debug.printfunctionstack('>');

result = 0;
switch nargin
	case 2,
		[IMGDIR, fname] = fileparts(arg1);
		outpath = arg1;
		res = arg2;
	case 3,
		[IMGDIR, fname] = deal(arg1, arg2);
		outpath = fullfile(arg1, sprintf('%s.png',arg2));
		res = arg3;
    otherwise,
        debug.print_debug(0, 'could not save image file - wrong number of arguments')
        return;
end

if ~exist(IMGDIR,'dir')
	mkdir(IMGDIR);
end

try
	% Save the image file
	print(gcf, '-dpng', sprintf('-r%d',res), outpath );

	% Lock the output directory if in archive mode - it means we do not want to delete these files ever
	if strcmp('PARAMS.mode','archive')
		system(sprintf('touch %s/lock',IMGDIR));
	end

	% Did our image file actually get saved?
	if exist(outpath, 'file')
		debug.print_debug(1, sprintf('%s: Saved image file %s',datestr(utnow),outpath) );
		result = 1;
	else
		debug.print_debug(1, sprintf('%s: Image file %s was not created',datestr(utnow),outpath));
	end
catch
	debug.print_debug(0, sprintf('%s: Could not save the image file %s',datestr(utnow),outpath));
end

%print_debug(sprintf('< %s',mfilename),2);
debug.printfunctionstack('<');


text:
%LIBSVC Trainable classifier: LIBSVM
%
% 	[W,J] = LIBSVC(A,KERNEL,C)
% 	[W,J] = A*LIBSVC([],KERNEL,C)
% 	[W,J] = A*LIBSVC(KERNEL,C)
%
% INPUT
%   A	    Dataset
%   KERNEL  Mapping to compute kernel by A*MAP(A,KERNEL)
%           or string to compute kernel by FEVAL(KERNEL,A,A)
%           or cell array with strings and parameters to compute kernel by
%           FEVAL(KERNEL{1},A,A,KERNEL{2:END})
%           Default: linear kernel (PROXM([],'P',1))
%   C       Trade_off parameter in the support vector classifier.
%           Default C = 1;
%
% OUTPUT
%   W       Mapping: Support Vector Classifier
%   J       Object idences of support objects. Can be also obtained as W{4}
%
% DESCRIPTION
% Optimizes a support vector classifier for the dataset A by the libsvm
% package, see http://www.csie.ntu.edu.tw/~cjlin/libsvm/. LIBSVC calls the
% svmtrain routine of libsvm for training. Classifier execution for a
% test dataset B may be done by D = B*W; In D posterior probabilities are
% given as computed by svmpredict using the '-b 1' option.
%
% The kernel may be supplied in KERNEL by
% - an untrained mapping, e.g. a call to PROXM like W = LIBSVC(A,PROXM([],'R',1))
% - a string with the name of the routine to compute the kernel from A
% - a cell-array with this name and additional parameters.
% This will be used for the evaluation of a dataset B by B*W or PRMAP(B,W) as
% well.
%
% If KERNEL = 0 (or not given) it is assumed that A is already the
% kernelmatrix (square). In this also a kernel matrix should be supplied at
% evaluation by B*W or PRMAP(B,W). However, the kernel has to be computed with
% respect to support objects listed in J (the order of objects in J does matter).
%
% EXAMPLE
% a = gendatb;                     % generate banana classes
% [w,J] = a*libsvc(proxm('p',3));  % compute svm with 3rd order polynomial
% a*w*testc                        % show error on train set
% scatterd(a)                      % show scatterplot
% plotc(w)                         % plot classifier
% hold on;
% scatterd(a(J,:),'o')             % show support objcts
%
% REFERENCES
% R.-E. Fan, P.-H. Chen, and C.-J. Lin. Working set selection using the second order
% information for training SVM. Journal of Machine Learning Research 6, 1889-1918, 2005
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, SVC, PROXM

% 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 [W,J,u] = libsvc(varargin)

	checktoolbox('libsvm');

  mapname = 'LIBSVM';
  argin = shiftargin(varargin,{'prmapping','char','cell'});
  argin = setdefaults(argin,[],proxm([],'p',1),1,1);

  if mapping_task(argin,'definition')

    W = define_mapping(argin,'untrained',mapname);

	elseif mapping_task(argin,'training')			% Train a mapping.

    [a,kernel,C] = check_for_old_call(argin);
    opt = ['-s 0 -t 4 -b 1 -e 1e-3 -c ',num2str(C), ' -q'];

    islabtype(a,'crisp');
    isvaldfile(a,1,2); % at least 1 object per class, 2 classes
    a = testdatasize(a,'objects');
    [m,k,c] = getsize(a);
    nlab = getnlab(a);

    K = compute_kernel(a,a,kernel);
    K = min(K,K');   % make sure kernel is symmetric
    K = [[1:m]' K];  % as libsvm wants it
                     % call libsvm
    u = svmtrain(nlab,K,opt);
		if isempty(u)
			prwarning(1,'libsvc: no solution for SVM, pseudo-inverse will be used')
			W = lkc(prdataset(K(:,2:end),getlabels(a)),0);
			J = [1:m]';
			return
		end
                     % Store the results:
    J = full(u.SVs);
    if isempty(J) | J == 0
      % LIBSVM failed, use fisher
      W = fisherc(a);
      prwarning(1,'LIBSVC failed, Fisher used instead')
      return
    end
    if isequal(kernel,0)
      s = [];
      in_size = 0; % to allow old and new style calls
    else
      s = a(J,:);
      in_size = k;
    end

    lablist = getlablist(a);
    W = prmapping(mfilename,'trained',{u,s,kernel,J,opt},lablist(u.Label,:),in_size,c);

    W = setname(W,'LIBSVM Classifier');
    W = setcost(W,a);

  else % Evaluation

    [a,W] = deal(argin{1:2});
    [u,s,kernel,J,opt] = getdata(W);
    m = size(a,1);

    K = compute_kernel(a,s,kernel);
    k = size(K,2);
    if k ~= length(J)
      if isequal(kernel,0)
        if (k > length(J)) &  (k >= max(J))
          % precomputed kernel; old style call
          prwarning(2,'Old style execution call: The precomputed kernel was calculated on a test set and the whole training set!')
        else
          error(['Inappropriate precomputed kernel!' newline ...
              'For the execution the kernel matrix should be computed on a test set' ...
              newline 'and the set of support objects']);
        end
      else
        error('Kernel matrix has the wrong number of columns');
      end
    else
      % kernel was computed with respect to the support objects
      % we make an approprite correction in the libsvm structure
      u.SVs = sparse((1:length(J))');
    end
    K = [[1:m]' K];  % as libsvm wants it
    %[lab,acc,d] = svmpredict(getnlab(a),K,u,' -b 1');
    [lab,acc,d] = svmpredict(ones(m,1),K,u,' -b 1');
    W = setdat(a,d,W);
  end

return;

function K = compute_kernel(a,s,kernel)

	% compute a kernel matrix for the objects a w.r.t. the support objects s
	% given a kernel description

	if  isstr(kernel) % routine supplied to compute kernel
		K = feval(kernel,a,s);
	elseif iscell(kernel)
		K = feval(kernel{1},a,s,kernel{2:end});
	elseif ismapping(kernel)
		K = a*prmap(s,kernel);
	elseif kernel == 0 % we have already a kernel
		K = a;
	else
		error('Do not know how to compute kernel matrix')
	end

	K = +K;

return

function [a,kernel,C] = check_for_old_call(argin)

[a,kernel,C,par] = deal(argin{:});
if ischar(kernel) && exist(kernel,'file') ~= 2
  kernel = proxm(kernel,C);
  C = par;
end

author: marianux

repo: ecg-kit

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


text:
function m = meltpdd(dem,T,minalt,params)
% PURPOSE: calculate melt for a DEM based on a simple temperature index model
%          using positive degree days (PDD).
% -------------------------------------------------------------------
% USAGE: m = meltpdd(dem,T,minalt,params)
% where: [dem] is the input topography grid
%        [T] is the mean annual air temperature MAAT at sealevel in degC
%        [minalt] (optional) is the minimum altitude to calculate melt for
%                 (default min(dem(:)) )
%        [params] (optional) is a vector containing the 4 parameters:
%               ddf_ice: degree day factor for ice [m/d degC]
%               lrate:   temperature lapse rate in deg/100m
%               a_seas:  annual temp fluctuation in degrees.
%               tau_a:   length of the year in days to calculate annual fluctuation
%                        for, should be 365
%        [params] default to [8.0,0.6,5.0,365];
%        Values and model taken from
%        SJ Marshall and GKC Clarke 1999: Ice sheet inception: subgrid
%        hypsometric parameterisation of mass balance in an ice sheet
%        model. Climate Dynamics 15:533-550
%
%% -------------------------------------------------------------------------
% OUTPUTS:
%        m is an array with melt values for evey grid cell in dem above
%        minalt
% -------------------------------------------------------------------
% NOTE:  Seasonal variation is calculated using cos(2*pi*d/tau_a)
%        for d 1:tau_a.
%        Melt for areas lower than [minalt] is set to 0.
%
% EXAMPLE: m=meltpdd(randn(100,100)*1000,13,500,[8.0,0.6,5.0,365])
%
%% Felix Hebeler, Geography Dept., University Zurich, Jan 2007.

if nargin < 2
    error('This function needs 2 input parameters (DEM,T)');
end

% default parameter values
ddf_ice  = 8.0; % DDF f2or ice in m/d degC
%ddf_snow = 0.003;
lrate = 0.6;  % temperatur lapse rate in deg/100m
a_seas   = 5.0;   % annual amplitude of temperature fluctuations in degrees
tau_a    = 365;   %length of the year in days used in PDD

if ~exist('minalt','var')
    minalt=min(dem(:));
end

if ~exist('params','var')
    params=[ddf_ice,lrate,a_seas,tau_a];
end
params(1)=params(1)/1000;
% initialise
T = T - (dem*params(2)/100); % adjust T to be surface temperature for the DEM
%tempm=zeros(size(dem,1),size(dem,2)); % allocate matrix to hold daily melt
m = zeros(size(dem,1),size(dem,2)); % matrix for meltsum
% calc melt
for d=1:params(4)  % seasonal variation in melt
        % tempm = (T - (a_seas * cos(2*pi*d/tau_a)))*ddf_ice;
        tempm = (T - (params(3)*cos(2*pi*d/params(4))))*params(1);
        tempm(tempm<0)=0;
        m = m + tempm;
end
clear tempm;
m(dem<minalt)=0;
clear dem;


text:
function [stack, norm_value] = ReadLDRStack(dir_name, format, bNormalization, bToSingle)
%
%       [stack, norm_value] = ReadLDRStack(dir_name, format, bNormalization, bToSingle)
%
%       This function reads an LDR stack from a directory, dir_name, given
%       an image format.
%
%        Input:
%           -dir_name: the path where the stack is
%           -format: the LDR format of the images that we want to load in
%           the folder dir_name. For example, it can be 'jpg', 'jpeg',
%           'png', 'tiff', 'bmp', etc.
%           -bNormalization: is a flag for normalizing or not the stack in
%           [0, 1].
%           -bToSingle:
%
%        Output:
%           -stack: a stack of LDR images, in floating point (single)
%           format. No normalization is applied.
%           -norm_value:
%
%     This function reads a stack of images from the disk
%
%     Copyright (C) 2011-16  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

if(~exist('bNormalization', 'var'))
    bNormalization = 0;
end

if(~exist('bToSingle', 'var'))
    bToSingle = 1;
end

if(bNormalization)
    bToSingle = 1;
end

norm_value = 1.0;

list = dir([dir_name, '/*.', format]);
n = length(list);

if(n > 0)
    img_info = [];
    name = [dir_name, '/', list(1).name];

    try
        if(exist('imfinfo') == 2)
            img_info = imfinfo(name);
        end
    catch err
        disp(err);

        try
            if(exist('exifread') == 2)
                img_info = exifread(name);

                if(isfield(img_info, 'SamplesPerPixel'))
                    if(img_info.SamplesPerPixel == 3)
                        img_info.ColorType = 'truecolor';
                    end

                    if(img_info.SamplesPerPixel == 1)
                        img_info.ColorType = 'grayscale';
                    end
                else
                    img_info.ColorType = 'truecolor';
                end

                if(isfield(img_info, 'BitsPerSample'))
                    img_info.BitDepth = round(mean(img_info.BitsPerSample));
                else
                    img_info.BitDepth = 8;
                end

                if(isfield(img_info, 'ImageWidth'))
                    img_info.Width = img_info.ImageWidth;
                end

                if(isfield(img_info, 'ImageLength'))
                    img_info.Height = img_info.ImageLength;
                end

                if(isfield(img_info, 'PixelXDimension'))
                    img_info.Width = PixelXDimension;
                end

                if(isfield(img_info, 'PixelYDimension'))
                    img_info.Height = PixelYDimension;
                end
            end
        catch
            disp(err);
        end
    end

    colorChannels = 0;

    norm_value = 255.0;

    if(~isempty(img_info))
        if(isfield(img_info, 'NumberOfSamples'))
            colorChannels = img_info.NumberOfSamples;
        else
            switch img_info.ColorType
                case 'grayscale'
                    colorChannels = 1;

                    switch img_info.BitDepth
                        case 8
                            norm_value = 255.0;
                        case 16
                            norm_value = 65535.0;
                    end

                case 'truecolor'
                    colorChannels = 3;

                    switch img_info.BitDepth
                        case 24
                            norm_value = 255.0;
                        case 48
                            norm_value = 65535.0;
                    end
            end
        end

        if(img_info.Height == 0 || img_info.Width == 0)
            tmp = imread(name);
            [img_info.Height, img_info.Width, n] = size(tmp);
            clear('tmp');
        end

        stack = zeros(img_info.Height, img_info.Width, colorChannels, n, 'single');

        for i=1:n
            disp(list(i).name);
            %read an image, and convert it into floating-point
            img_tmp = imread([dir_name, '/', list(i).name]);

            %store in the stack
            if(bToSingle)
                stack(:,:,:,i) = single(img_tmp);
            else
                stack(:,:,:,i) = img_tmp;
            end
        end

        if(bNormalization)
            stack = stack / norm_value;
        end
    else
        warning('The stack is empty!');
        stack = [];
    end
else
    warning('The stack is empty!');
    stack = [];
end

end

author: banterle


text:
function [surfDice,addedPathLen,addedPathLenNorm] = ...
    calc_surfDice_addPathLen(structNum1,structNum2,margin,planC)
% function [surfDice,addedPathLen,addedPathLenNorm] = ...
%     calc_surfDice_addPathLen(structNum1,structNum2,margin,planC)

mask1M = getSurfaceRing(structNum1,margin,planC);
mask2M = getSurfaceRing(structNum2,margin,planC);
avgMask = (sum(mask1M(:))+sum(mask2M(:)))/2;
intrsctM = mask1M & mask2M;
surfDice = sum(intrsctM(:))/avgMask;

% addedPathLenNorm = sum((mask1M(:) | mask2M(:)) & ~intrsctM(:))/(avgMask*2);
% addedPathLen = sum((mask1M(:) | mask2M(:)) & ~intrsctM(:));

scanNum = getStructureAssociatedScan(structNum1,planC);
% [~,~,numSlcs] = size(intrsctM);

indexS = planC{end};

[xScanV,yScanV,zScanV] = getScanXYZVals(planC{indexS.scan}(scanNum));
[~,~,zUnifScanV] = getUniformScanXYZVals(planC{indexS.scan}(scanNum));
numSlcs = length(zScanV);

structLen = 0;
intersectStructLen = 0;
for slc = 1:numSlcs
    if length(planC{indexS.structures}(structNum1).contour(slc).segments) ~= 1
        continue
    end
    xV = planC{indexS.structures}(structNum1).contour(slc).segments(1).points(:,1);
    yV = planC{indexS.structures}(structNum1).contour(slc).segments(1).points(:,2);
    %[rowV, colV] = xytom(xV, yV, slc, planC,scanNum);
    siz = size(intrsctM(:,:,slc));
    %indV = sub2ind(siz, round(rowV), round(colV));

    structXYv = [xV';yV'];

    d = hypot(diff(xV), diff(yV)); % Distance Of Each Segment
    contourLen = sum(d);

    structLen = structLen + contourLen;

    unifSlc = findnearest(zUnifScanV,zScanV(slc));
    cc = bwconncomp(intrsctM(:,:,unifSlc));
    intersectLen = 0;
    for comp = 1:cc.NumObjects
        compIndV = cc.PixelIdxList{comp};
        [compiV,compjV] = ind2sub(siz,compIndV);
        compXYv = [xScanV(compjV);yScanV(compiV)];

        distM = sepsq(structXYv,compXYv);
        [~,indMinV] = min(distM,[],1);
        indMinV = unique(indMinV);
        segX = xV(indMinV);
        segY = yV(indMinV);

        d = hypot(diff(segX), diff(segY)); % Distance Of Segment between points
        segDist = sum(d);

        intersectLen = intersectLen + segDist;

    end

    intersectStructLen = intersectStructLen + intersectLen;

end

addedPathLen = structLen - intersectStructLen;
addedPathLenNorm = addedPathLen / structLen;


text:
%  Mex code here based on GrabBitmaps.cpp in DirectShow samples
%

MS_platform_sdk = 'C:\Program Files\Microsoft Platform SDK for Windows Server 2003 R2\';

include_dirs = [' -DMSWIN -I''', MS_platform_sdk, '\Samples\Multimedia\DirectShow\BaseClasses'' '];
include_dirs = [include_dirs, ' -I''', MS_platform_sdk, '\Include'''];
include_dirs = [include_dirs, ' -I''C:\Program Files (x86)\Microsoft DirectX SDK (November 2007)\Include'''];

libraries = [' ''', MS_platform_sdk, '\Samples\Multimedia\DirectShow\BaseClasses\Release\strmbase.lib'''];
libraries = [libraries, ' zlib.lib winmm.lib'];
libraries = [libraries, ' Filters.cpp VidHeader.cpp showErrMsgBox.cpp read-frame.cpp'];

mex_cpps = {'dxAviOpenMex.cpp', 'dxConvertToAviMex.cpp', 'dxAviReadMex.cpp', 'dxAviCloseMex.cpp'};
for mexi = 1:length(mex_cpps);
    eval(['mex -O ', mex_cpps{mexi}, ' ', include_dirs, libraries]);
end


text:
function sparse_test02 ( )

%*****************************************************************************80
%
%% SPARSE_TEST02 demonstrates a simple use of the SPARSE matrix facility.
%
%  Discussion:
%
%    This is a nice but very simple example.
%
%    Note in this example that we define three sparse matrices,
%    SUP, DIAG, and SUB, and then define A as the sum of those matrices.
%
%    However, again, each time we define a sparse matrix, we are assuming
%    we have all the information available at one time.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 August 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    George Lindfield, John Penny,
%    Numerical Methods Using MATLAB,
%    Prentice Hall, 1999
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_TEST02:\n' );
  fprintf ( 1, '  Demonstrate the use of MATLAB''s SPARSE facility\n' );
  fprintf ( 1, '  to define a sparse matrix, and solve an associated\n' );
  fprintf ( 1, '  linear system.\n' );

  n = 100;
%
%  We set up the three diagonals of the -1, 2, -1 matrix.
%
  sup  = sparse ( 1:n-1, 2:n,   -1.0, n, n );
  diag = sparse ( 1:n,   1:n,    2.0, n, n );
  sub  = sparse ( 2:n,   1:n-1, -1.0, n, n );
%
%  A is the matrix whose superdiagonal, diagonal, and subdiagonal
%  have been set above.  Because SUP, DIAG and SUB are sparse,
%  A will "inherit" this property.
%
  A = sup + diag + sub;
%
%  Set up a right hand side b that defines a linear system whose
%  solution is [ 1, 2, 3, ..., n ].
%
  b(1:n-1,1) = 0.0;
  b(n,1) = n + 1;
%
%  Use MATLAB's backslash command to solve the linear system.
%  Because MATLAB "knows" that A is a sparse matrix, it will
%  automatically do the efficient thing.
%
  x = A \ b;
%
%  Print the solution.
%
  x(1:n)

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [mbar] = hPa2mbar(hPa)
% Convert pressure from hectopascals to millibars.
% They're the same!
% Chad Greene 2012
mbar = hPa;