text:
%This Matlab script can be used to reproduce Figure 3.8 in the monograph:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Empty workspace and close figures
close all;
clear;

%Define number of UEs per cell
K = 10;

%Define the range of BS antennas
Mrange = 10:10:100;

%Define the pilot length
tau_p = K;

%Complexity of MMSE channel estimation, based on Table 3.1
complexity_MMSE = Mrange*tau_p + Mrange.^2;

%Complexity of EW-MMSE channel estimation, based on Table 3.1
complexity_EW_MMSE = Mrange*(tau_p + 1);

%Complexity of LS channel estimation, based on Table 3.1
complexity_LS = Mrange*tau_p;

%% Plot the simulation results
figure;
hold on; box on;

plot(Mrange,complexity_MMSE,'r--','LineWidth',1);
plot(Mrange,complexity_EW_MMSE,'k-','LineWidth',1);
plot(Mrange,complexity_LS,'b-.','LineWidth',1);

xlabel('Number of antennas (M)');
ylabel('Number of complex multiplications');
set(gca,'YScale','log');

legend('MMSE','EW-MMSE','LS','Location','NorthWest');

author: emilbjornson

repo: massivemimobook


text:
function [rectxt] = recoverbp(I,siz,b)
% RECOVERBP Recovers the text embedded in the Bth bit-plane of the image
%
% RECTXT = RECOVERBP(I,SIZ,B) returns the text string of length SIZ
% embedded in the Bth bit-plane of the image I. If B is not specified, it
% is taken as 1.
%
% See also EMBEDBP

if nargin == 2
    b = 1;
end
bsiz = 8*siz;
n = numel(I);
if bsiz > n
    error('Size of text given exceeds the maximum that can be embedded in the image')
    return
end
dim = size(I);
addl = n-bsiz;
I1 = reshape(I,1,n);
I2 = round(abs(I1(1:bsiz)));
p = 2^b;
h = 2^(b-1);
rb = zeros(1,bsiz);
for k = 1:bsiz
    I2(k) = round(I2(k));
    r = rem(I2(k),p);
    if r >= h
        rb(k) = 1;
    end
end
rbi = (dec2bin(rb,1))';
rbin = reshape(rbi,siz,8);
rectxt = (bin2dec(rbin))';
return

author: Sable

repo: mcbench-benchmarks


text:
function scanSummationMenu(command)
%"scanSummationMenu"
%   Create the GUI used to sum scans.
%
%   02/15/08  APA Based on scanSummationMenu
%
%Usage:
%   scanSummationMenu()
%
% 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;
global stateS;
indexS = planC{end};

if nargin == 0
    command = 'init';
end

hFig = findobj('tag', 'CERRScanSummationFigure');

switch upper(command)
    case 'INIT'
        if ~isempty(hFig)
            delete(hFig);
        end
        screenSize = get(0,'ScreenSize');
        units = 'pixels';
        y = 480;
        x = 480;

        dx = 20;
        dy = 20;

        hFig = figure('Name', 'Scan Summation', 'units', 'pixels', 'position',[(screenSize(3)-x)/2 (screenSize(4)-y)/2 x y], 'MenuBar', 'none', 'NumberTitle', 'off', 'resize', 'off', 'Tag', 'CERRScanSummationFigure');
        stateS.handle.scanSummationMenuFig = hFig;

        bartop = 80;
        loadbot = y-65-10;

        ud.handles.instruction = uicontrol('units',units,'HorizontalAlignment','left','Position',[dx y-dy-60 x-2*dx 20],'String', 'Select the scans to be included in the sum scan.', 'Style','text','Tag','text1');
        ud.handles.labelname = uicontrol('units',units,'HorizontalAlignment','left','Position',[dx y-dy-10 100 15],'String', 'New Scan Name:', 'Style','text','Tag','text2');
        ud.handles.newScanName = uicontrol('units',units,'Position',[dx y-dy-35 120 20],'BackgroundColor',[1 1 1],'HorizontalAlignment','left','String', '', 'Style','edit','Tag','scanName');
        ud.handles.assocname = uicontrol('units',units,'HorizontalAlignment','left','Position',[x-dx-250 y-dy-10 140 15],'String', 'Prefer Images from:', 'Style','text','Tag','scanAssoc');
        for i=1:length(planC{indexS.scan})
            scanStr{i} = [num2str(i), '. ', planC{indexS.scan}(i).scanType];
        end
        ud.handles.scanSelect = uicontrol('units',units,'Position',[x-dx-250 y-dy-35 120 20],'BackgroundColor',[1 1 1],'HorizontalAlignment','left','String', scanStr, 'Style','popup','value',1,'Tag','scanSelect');

        %Make CANCEL and SUM buttons
        ud.handles.cancelButton = uicontrol(hFig, 'callback', 'scanSummationMenu(''CANCEL'');', 'units',units,'Position',[x-dx-70 y-dy-35 70 20],'String','Cancel', 'Style','pushbutton','Tag','cancelButton');
        ud.handles.sumButton = uicontrol(hFig, 'callback', 'scanSummationMenu(''SUM'');', 'units',units, 'Position',[x-dx-70 y-dy-10 70 20],'String','Sum', 'Style','pushbutton','Tag','sumButton');

        %Plans frame and label
        planFrame  = uicontrol(hFig, 'units',units,'Position',[dx dy x-2*dx y-dy-bartop], 'style', 'frame');
        frameColor = get(planFrame, 'BackgroundColor');

        %Create scrollbar on right side. Inactive to start.
        ud.df.handles.scroll = uicontrol(hFig, 'units',units,'Position',[x-dx-20 dy 20 y-dy-bartop], 'style', 'slider', 'enable', 'off','tag','planSlider', 'callback', 'scanSummationMenu(''SLIDER'')');

        %Create plan UICONTROLS.
        uicontrol(hFig, 'units',units,'Position',[30 loadbot-46 15 15], 'style', 'text', 'String', 'No.', 'fontweight', 'bold', 'backgroundcolor', frameColor);
        %uicontrol(hFig, 'units',units,'Position',[40 loadbot-46 100 15], 'style', 'text', 'String', 'Plan Name', 'fontweight', 'bold', 'backgroundcolor', frameColor);
        uicontrol(hFig, 'units',units,'Position',[dx+100 loadbot-46 70 15], 'style', 'text', 'String', 'Assoc Scan', 'fontweight', 'bold', 'backgroundcolor', frameColor);
        uicontrol(hFig, 'units',units,'Position',[dx+210 loadbot-46 60 15], 'style', 'text', 'String', 'Select', 'fontweight', 'bold', 'backgroundcolor', frameColor);
        uicontrol(hFig, 'units',units,'Position',[dx+270 loadbot-46 100 15], 'style', 'text', 'String', 'Weighting Factor', 'fontweight', 'bold', 'backgroundcolor', frameColor);

        numRows = 50;
        for i=1:numRows
            ud.handles.scanNum(i)   = uicontrol(hFig, 'units',units,'Position',[30 359-i*20 10 15], 'style', 'text', 'String', num2str(i), 'backgroundcolor', frameColor, 'visible', 'off');
            ud.handles.scanName(i)  = uicontrol(hFig, 'units',units,'Position',[dx+100 359-i*20 40 15], 'style', 'text', 'String', 'Scan', 'backgroundcolor', frameColor, 'visible', 'off');
            ud.handles.scanCheck(i) = uicontrol(hFig, 'units',units,'Position',[dx+210 359-i*20 25 15], 'style', 'checkbox', 'backgroundcolor', frameColor, 'horizontalAlignment', 'center', 'userdata', i, 'Tag', 'ScanCheck', 'callback', 'scanSummationMenu(''Check'')', 'visible', 'off');
            ud.handles.wtfactor(i)  = uicontrol(hFig, 'units',units,'Position',[dx+270 359-i*20 45 15], 'Style','edit', 'String', '1.0', 'BackgroundColor',[1 1 1],'horizontalAlignment', 'center','Tag','wtfactor', 'callback', 'scanSummationMenu(''Check'')', 'visible', 'off');
        end
        ud.numScanRows = numRows;

        ud.firstVisScan     = 1;

        set(hFig, 'userdata', ud);

        nScans = length(planC{indexS.scan});
        ud.checkedScans = zeros(nScans, 1);

        %Set range of currently displayed scans.
        ud.df.range = 1:min(16, nScans);

        set(hFig, 'userdata', ud);
        scanSummationMenu('refresh');

    case 'REFRESH'
        ud = get(hFig, 'userdata');
        indexS = planC{end};

        %Number of visible scans.
        nvScans = length(ud.df.range);

        %Refresh scan UIelements.
        set(ud.handles.scanNum,  'visible', 'off');
        set(ud.handles.scanName,  'visible', 'off');
        set(ud.handles.scanCheck, 'visible', 'off');
        set(ud.handles.wtfactor, 'visible', 'off');
        nScans = length(ud.checkedScans);
        ud.firstVisScan = min(ud.df.range);
        for i=1:min(length(ud.handles.scanName), nvScans)
            scanNum  = ud.firstVisScan+i-1;
            scanName = [num2str(scanNum),'. ',planC{indexS.scan}(scanNum).scanType];
            set(ud.handles.scanNum(scanNum),'visible', 'on')
            set(ud.handles.scanName(scanNum), 'string', scanName, 'visible', 'on', 'Position',[20+100 359-i*20 80 15]);
            set(ud.handles.scanCheck(scanNum), 'visible', 'on', 'Position',[20+240 359-i*20 25 15]);
            set(ud.handles.wtfactor(scanNum), 'visible', 'on', 'Position',[20+300 359-i*20 45 15]);
        end

        if nScans > 16
            set(ud.df.handles.scroll, 'min', 0, 'max', nScans-nvScans, 'value', nScans-nvScans+1-min(ud.df.range), 'enable', 'on', 'sliderstep', [1/(nScans-nvScans), nvScans/(nScans-nvScans)]);
        else
            set(ud.df.handles.scroll, 'enable', 'off');
        end

        set(hFig, 'userdata', ud);

    case 'CHECK'
        checkNum = get(gcbo, 'userdata');
        type     = get(gcbo, 'Tag');
        value    = get(gcbo, 'value');

        ud = get(hFig, 'userdata');

        switch upper(type)
            case 'SCANCHECK'
                ud.checkedScans(checkNum+ud.firstVisScan-1) = value;
                if value == 1
                    set(ud.handles.newScanName,'string',planC{indexS.scan}(checkNum).scanType)
                    set(ud.handles.scanSelect,'value',checkNum)
                end
%             case 'WTFACTOR'
%                 factor = str2num(get(ud.handles.wtfactor(i), 'String'));
        end

        set(hFig, 'userdata', ud);

    case 'SLIDER'
        %Slider was clicked, move ud.df.range.
        ud = get(hFig, 'userdata');
        val = round(get(gcbo, 'value'));

        nRows = 16;
        nScans  = length(planC{indexS.scan});

        lastScan = nScans - val;
        ud.df.range = max(1, lastScan-nRows+1):lastScan;
        set(hFig, 'userdata', ud);
        scanSummationMenu('REFRESH');

    case 'SUM'
        hFig = findobj('tag', 'CERRScanSummationFigure');
        ud = get(hFig, 'userdata');
        checkedScans = find(ud.checkedScans);
        if length(checkedScans) < 2
			warndlg('You must select atleast two scans','Incorrect Selection','modal')
            return
        end

        % Check for duplicate or blank name.
        newScanName = get(ud.handles.newScanName, 'String');
        if isempty(newScanName)
            warndlg('Please enter a name for the new scan.', 'Warning','modal');
            return;
        end

        %Check that no checked plan has a weighting factor that is blank.
        for i = 1:length(planC{indexS.scan})
            factor = str2num(get(ud.handles.wtfactor(i), 'String'));
            if isempty(factor) && ~isempty(find(checkedScans == i))
                warndlg('Please enter a weighting factor for all selected scans.', 'Warning','modal');
                return
            elseif isempty(factor)
                factor = 0;
            end
            wtfactor(i) = factor;
        end

        %% APA code begins
        assocScan = get(ud.handles.scanSelect,'value');
        [xAssocV, yAssocV, zAssocV] = getScanXYZVals(planC{indexS.scan}(assocScan));

        if assocScan > 0
            assocScanUID = planC{indexS.scan}(assocScan).scanUID;
            %Get associated transM
            assocTransM = planC{indexS.scan}(assocScan).transM;
            if isempty(assocTransM)
                assocTransM = eye(4);
            end
        else %No Association
            assocScanUID = '';
            assocTransM = eye(4);
        end

        scanNums = checkedScans;

        %Get the x,y,z grid for new scan
        for i = 1:length(scanNums)
            scanNum = scanNums(i);
            %Get x,y,z values for scanNum
            [xV, yV, zV] = getScanXYZVals(planC{indexS.scan}(scanNum));

            %Get the corners of the original dataset.
            [xCorn, yCorn, zCorn] = meshgrid([min(xV) max(xV)], [min(yV) max(yV)], [min(zV) max(zV)]);

            %Add ones to the corners so we can apply a transformation matrix.
            corners = [xCorn(:) yCorn(:) zCorn(:) ones(prod(size(xCorn)), 1)];

            %Apply transform to corners, so we know boundary of the slice.
            transM = getTransM(planC{indexS.scan}(scanNum),planC);
            if isempty(transM) || isequal(transM,eye(4))
                [xV,yV,zV] = getScanXYZVals(planC{indexS.scan}(scanNum));
                xGrid{scanNum} = xV(:)';
                yGrid{scanNum} = yV(:)';
                zGrid{scanNum} = zV(:)';
            else
                newCorners = inv(assocTransM) * transM * corners';
                xGrid{scanNum} = linspace(min(newCorners(1,:)), max(newCorners(1,:)), length(xV));
                yGrid{scanNum} = linspace(max(newCorners(2,:)), min(newCorners(2,:)), length(yV));
                zGrid{scanNum} = linspace(min(newCorners(3,:)), max(newCorners(3,:)), length(zV));
            end
            xRes{scanNum} = length(xV);
            yRes{scanNum} = length(yV);
            zRes{scanNum} = length(zV);

            %Get associated scan
            assocScanV{scanNum} = scanNum;
        end
        newXgrid = linspace(min(cell2mat(xGrid)),max(cell2mat(xGrid)),max(cell2mat(xRes)));
        newYgrid = linspace(max(cell2mat(yGrid)),min(cell2mat(yGrid)),max(cell2mat(yRes)));
        newZgrid = linspace(min(cell2mat(zGrid)),max(cell2mat(zGrid)),round(max(cell2mat(zRes))*1.5));

%         %Obtain scans with same grid
%         if isempty(assocTransM)
%             assocTransM = eye(4);
%         end
%         scanIndC = {};
%         scanNumsTmp = scanNums;
%         for iSortAll = 1:length(scanNums)
%             iSort = scanNums(iSortAll);
%             indRemaining = scanNums;
%             indRemaining(iSortAll) = [];
%             scanSortM = [iSort];
%             if ~ismember(iSort,[scanIndC{:}])
%                 for jSortAll = 1:length(indRemaining)
%                     jSort = indRemaining(jSortAll);
%                     if ~isempty(getTransM('scan',iSort,planC))
%                         scanITM = inv(assocTransM) * getTransM('scan',iSort,planC);
%                     else
%                         scanITM = inv(assocTransM);
%                     end
%                     if ~isempty(getTransM('scan',jSort,planC))
%                         scanJTM = inv(assocTransM) * getTransM('scan',jSort,planC);
%                     else
%                         scanJTM = inv(assocTransM);
%                     end
%                     if isequal([xRes{iSort},yRes{iSort},zRes{iSort}],[xRes{jSort},yRes{jSort},zRes{jSort}]) && isequal(scanITM,scanJTM)
%                         scanSortM(end+1) = jSort;
%                         indJsort = find(scanNumsTmp == jSort);
%                         scanNumsTmp(indJsort) = [];
%                     end
%                 end
%             end
%             scanIndC{iSort} = scanSortM;
%         end
%
%         scanNums = scanNumsTmp;

        %Loop over scans and add over new grid
        hWait = waitbar(0,['Appending one scans to Scan # ',num2str(assocScan)]);
        scanSumM = zeros([length(newYgrid),length(newXgrid),length(newZgrid)],'single');
        scanEmptyM = zeros([length(newYgrid),length(newXgrid)],'single');

        for i = 1:length(scanNums)

            scanNum = scanNums(i);

            %Check for transM
            if ~isempty(assocTransM) && ~isequal(assocTransM,eye(4))
                scanTtransM = getTransM('scan',scanNum,planC);
                if isempty(scanTtransM)
                    scanTtransM = eye(4);
                end
                inputTM = inv(assocTransM) * scanTtransM;
            else
                inputTM = getTransM('scan',scanNum,planC);
                if isempty(inputTM)
                    inputTM = eye(4);
                end
            end

%             %Get the summation for this grid
%             scanCombinedM = [];
%             for iScanAll = 1:length(scanIndC{scanNum})
%                 iScan = scanIndC{scanNum}(iScanAll);
%                 if ~isempty(scanCombinedM)
%                     scanCombinedM = scanCombinedM + wtfactor(iScan) * single(getScanArray(planC{indexS.scan}(iScan)));
%                 else
%                     scanCombinedM = wtfactor(iScan) * single(getScanArray(planC{indexS.scan}(iScan)));
%                 end
%             end

            %Check if this scan is on same grid as the new one
            %if isequal(xGrid{scanNum},newXgrid) && isequal(yGrid{scanNum},newYgrid) && isequal(zGrid{scanNum},newZgrid)
            chkLength = length(xGrid{scanNum})==length(newXgrid) && length(yGrid{scanNum})==length(newYgrid) && length(zGrid{scanNum})==length(newZgrid);
            if 0 && sum(sum((inputTM-eye(4)).^2)) < 1e-5 && chkLength && max(abs(xGrid{scanNum}-newXgrid)) < 1e-6 && max(abs(yGrid{scanNum}-newYgrid)) < 1e-6 && max(abs(zGrid{scanNum}-newZgrid)) < 1e-6

                scanSumM = scanSumM + scanCombinedM;
                waitbar((i-1)/length(scanNums),hWait,['Calculating contribution from Scan ', num2str(scanNum)])

            else %interpolation required

                %Transform this scan
                scanTmpM = [];

                for slcNum=1:length(newZgrid)
                    [xV, yV, zV] = getScanXYZVals(planC{indexS.scan}(scanNum));
                    if min(zAssocV)<newZgrid(slcNum) && newZgrid(slcNum)<max(zAssocV) && scanNum ~= assocScan
                        scanTmp = [];
                    else
                        scanTmp = slice3DVol(planC{indexS.scan}(scanNum).scanArray, xV, yV, zV, newZgrid(slcNum), 3, 'linear', inputTM, [], newXgrid, newYgrid);
                    end
                    if isempty(scanTmp)
                        scanTmpM(:,:,slcNum) = scanEmptyM;
                    else
                        scanTmpM(:,:,slcNum) = scanTmp;
                    end
                    waitbar((i-1)/length(scanNums) + (slcNum-1)/length(newZgrid)/length(scanNums) ,hWait,['Calculating contribution from Scan ', num2str(scanNum)])
                end

                scanSumM = scanSumM + wtfactor(scanNum) * scanTmpM;

            end

        end

        delete(hWait)

        %Create new scan distribution
        newScanNum = length(planC{indexS.scan}) + 1;
        newScanS = initializeCERR('scan');
        newScanS(1).scanArray = scanSumM;
        newScanS(1).scanType = newScanName;
        newScanS(1).scanUID = createUID('scan');
        newScanS(1).transM = planC{indexS.scan}(assocScan).transM;
        scanInfoS = planC{indexS.scan}(assocScan).scanInfo(1);
        scanInfoS.grid1Units = abs(newYgrid(1)-newYgrid(2));
        scanInfoS.grid2Units = abs(newXgrid(1)-newXgrid(2));
        scanInfoS.scanFileName = '';
        scanInfoS.DICOMHeaders = '';
        scanInfoS.sliceThickness =  abs(newZgrid(1)-newZgrid(2));
        scanInfoS.voxelThickness =  abs(newZgrid(1)-newZgrid(2));
        scanInfoS.sizeOfDimension1 = size(scanSumM,1);
        scanInfoS.sizeOfDimension2 = size(scanSumM,2);
        clear scanSumM
        for sInfoNum = 1:length(newZgrid)
            scanInfoNewS = scanInfoS;
            scanInfoNewS.zValue = newZgrid(sInfoNum);
            newScanS(1).scanInfo(sInfoNum) = scanInfoNewS;
        end
        planC{indexS.scan} = dissimilarInsert(planC{indexS.scan},newScanS,newScanNum);
        planC = setUniformizedData(planC);

        % Save scan statistics for fast image rendering
        for scanNum = newScanNum
            scanUID = ['c',repSpaceHyp(planC{indexS.scan}(scanNum).scanUID(max(1,end-61):end))];
            %stateS.scanStats.minScanVal.(scanUID) = single(min(planC{indexS.scan}(scanNum).scanArray(:)));
            %stateS.scanStats.maxScanVal.(scanUID) = single(max(planC{indexS.scan}(scanNum).scanArray(:)));
            minScan = single(min(planC{indexS.scan}(scanNum).scanArray(:)));
            maxScan = single(max(planC{indexS.scan}(scanNum).scanArray(:)));
            scanDiff = maxScan - minScan;
            scanCenter = (minScan + maxScan - 2*CTOffset) / 2;
            stateS.scanStats.CTLevel.(scanUID) = scanCenter;
            stateS.scanStats.CTWidth.(scanUID) = scanDiff;
            stateS.scanStats.Colormap.(scanUID) = 'gray256';
            stateS.scanStats.windowPresets.(scanUID) = 1;
        end

        %switch to new scan, with a short pause to let the dialogue clear.
        pause(.1);
        sliceCallBack('selectScan', num2str(newScanNum));

        %Refresh this GUI to include new scan distribution
        ud.df.range = [ud.df.range ud.df.range(end)+1];
        if length(ud.df.range)>16
            ud.df.range = ud.df.range(2:17);
        end
        set(hFig, 'userdata', ud);
        scanSummationMenu('refresh')

        %% APA code ends

    case 'CANCEL'
        hFig = findobj('tag', 'CERRScanSummationFigure');
        delete(hFig);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%% export_off
% Below is a demonstration of the features of the |export_off| function

%%
clear; close all; clc;

%% Syntax
% |export_off(fileName,F,V);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*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-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function Offspring = GLMO_SMPSOOperator(Problem, Particle, Pbest, Gbest, numberOfGroups, typeOfGroups)
% -----------------------------------------------------------------------
%  Copyright (C) 2020 Heiner Zille
%
%  This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0
%  International License. (CC BY-NC-SA 4.0). To view a copy of this license,
%  visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or see the
%  pdf-file "License-CC-BY-NC-SA-4.0.pdf" that came with this code.
%
%  You are free to:
%  * Share ? copy and redistribute the material in any medium or format
%  * Adapt ? remix, transform, and build upon the material
%  Under the following terms:
%  * Attribution ? You must give appropriate credit, provide a link to the
%     license, and indicate if changes were made. You may do so in any reasonable
%     manner, but not in any way that suggests the licensor endorses you or your use.
%  * NonCommercial ? You may not use the material for commercial purposes.
%  * ShareAlike ? If you remix, transform, or build upon the material, you must
%    distribute your contributions under the same license as the original.
%  * No additional restrictions ? You may not apply legal terms or technological
%    measures that legally restrict others from doing anything the license permits.
%
%  Author of this Code:
%   Heiner Zille <heiner.zille@ovgu.de> or <heiner.zille@gmail.com>
%
%  This code is based on the following publications:
%
%  1) Heiner Zille
%     "Large-scale Multi-objective Optimisation: New Approaches and a Classification of the State-of-the-Art"
%     PhD Thesis, Otto von Guericke University Magdeburg, 2019
%     http://dx.doi.org/10.25673/32063
%
%  2) Heiner Zille, Hisao Ishibuchi, Sanaz Mostaghim and Yusuke Nojima
%     "Mutation Operators Based on Variable Grouping for Multi-objective Large-scale Optimization"
%     IEEE Symposium Series on Computational Intelligence (SSCI), IEEE, Athens, Greece, December 2016
%     https://ieeexplore.ieee.org/document/7850214
%
%  This file is intended to work with the PlatEMO framework version 2.5.
%  Date of publication of this code: 06.04.2020
%  Last Update of this code: 06.04.2020
%  A newer version of this algorithm may be available. Please contact the author
%  or see http://www.ci.ovgu.de/Research/Codes.html.
%
% The files may have been modified in Feb 2021 by the authors of the Platemo framework to work with the Platemo 3.0 release.
% -----------------------------------------------------------------------
% This file is derived from its original version containied in the PlatEMO
% framework. The original copyright disclaimer can be found below.
% -----------------------------------------------------------------------

% Particle swarm optimization in SMPSO

    %% Parameter setting
    ParticleDec = Particle.decs;
    PbestDec    = Pbest.decs;
    GbestDec    = Gbest.decs;
    [N,D]       = size(ParticleDec);
    ParticleVel = Particle.adds(zeros(N,D));

    %% Particle swarm optimization
    W  = repmat(unifrnd(0.1,0.5,N,1),1,D);
    r1 = repmat(rand(N,1),1,D);
    r2 = repmat(rand(N,1),1,D);
    C1 = repmat(unifrnd(1.5,2.5,N,1),1,D);
    C2 = repmat(unifrnd(1.5,2.5,N,1),1,D);
    OffVel = W.*ParticleVel + C1.*r1.*(PbestDec-ParticleDec) + C2.*r2.*(GbestDec-ParticleDec);
    phi    = max(4,C1+C2);
    OffVel = OffVel.*2./abs(2-phi-sqrt(phi.^2-4*phi));
    delta  = repmat((Problem.upper-Problem.lower)/2,N,1);
    OffVel = max(min(OffVel,delta),-delta);
    OffDec = ParticleDec + OffVel;

    %% Deterministic back
    Lower  = repmat(Problem.lower,N,1);
    Upper  = repmat(Problem.upper,N,1);
    repair = OffDec < Lower | OffDec > Upper;
    OffVel(repair) = 0.001*OffVel(repair);
    OffDec = max(min(OffDec,Upper),Lower);

    %% Polynomial mutation
    disM  = 20;
    Site1 = repmat(rand(N,1)<0.15,1,D);

    [outIndexList,~] = GLMO_createGroups(numberOfGroups,OffDec,D,typeOfGroups); % 3 = random groups
    chosengroups = randi(numberOfGroups,size(outIndexList,1),1);
    Site2 = outIndexList == chosengroups;
    mu    = rand(N,1);
    mu    = repmat(mu,1,D);

    temp  = Site1 & Site2 & mu<=0.5;
    OffDec(temp) = OffDec(temp)+(Upper(temp)-Lower(temp)).*((2.*mu(temp)+(1-2.*mu(temp)).*...
                   (1-(OffDec(temp)-Lower(temp))./(Upper(temp)-Lower(temp))).^(disM+1)).^(1/(disM+1))-1);
    temp  = Site1 & Site2 & mu>0.5;
    OffDec(temp) = OffDec(temp)+(Upper(temp)-Lower(temp)).*(1-(2.*(1-mu(temp))+2.*(mu(temp)-0.5).*...
                   (1-(Upper(temp)-OffDec(temp))./(Upper(temp)-Lower(temp))).^(disM+1)).^(1/(disM+1)));
    Offspring = Problem.Evaluation(OffDec,OffVel);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function prov = hcp_provenance

% HCP_PROVENANCE returns a structure with provenance information

% Copyright (C) 2011-2014 by the Human Connectome Project, WU-Minn Consortium (1U54MH091657)
%
% This file is part of megconnectome.
%
% megconnectome 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.
%
% megconnectome 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 megconnectome.  If not, see <http://www.gnu.org/licenses/>.

prov.version.matlab = ver('MATLAB');

prov.version.megconnectome = ver('megconnectome');
if isempty(prov.version.megconnectome)
  prov.version.megconnectome = 'unknown'; % otherwise struct2xml fails
end

prov.version.fieldtrip = ver('fieldtrip');
if isempty(prov.version.fieldtrip)
  prov.version.fieldtrip = 'unknown'; % otherwise struct2xml fails
end

if isdeployed
  prov.compiled   = 'true';
else
  prov.compiled   = 'false';
end

% it is not desired to keep the date and time of execution of the analysis pipeline
% see https://github.com/Washington-University/megconnectome/issues/126
% prov.time           = datestr(now);

prov.username       = getusername;
prov.hostname       = gethostname;
prov.architecture   = computer('arch');
if isdeployed
  % the buildtimestamp function is created on the fly by the compile script
  prov.buildtimestamp = buildtimestamp;
else
  prov.buildtimestamp = [];
end
prov.pwd            = pwd;
prov.matlabstack    = printstack(dbstack); % using private helper function
prov.Attributes.xmlns_colon_xsi = 'http://www.w3.org/2001XMLSchema-instance';
prov.Attributes.xsi_colon_noNamespaceSchemaLocation = 'megconnectome.xsd';

function str = printstack(s)
str = sprintf('\n');
% skip the first, it is always hcp_provenance and not interesting
for i=2:length(s)
  str = [str sprintf('In %s at %d\n', s(i).name, s(i).line)];
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef SensorTargetEnum < matlab.mixin.SetGet
    %SensorTargetEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        PointTarget('Point Target');
        LatLongRectGrid('Lat/Long Rectangular Grid')
        LatLongCircGrid('Lat/Long Circular Grid');
    end

    properties
        name(1,:) char
    end

    methods
        function obj = SensorTargetEnum(name)
            obj.name = name;
        end
    end

    methods(Static)
        function listBoxStr = getListBoxStr()
            m = enumeration('SensorTargetEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('SensorTargetEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('SensorTargetEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [xV, yV] = modifyContour(hCERR)
%function [xV, yV] = modifyContour(hCERR)
%Modify the active contour
%Button 1 Moves the nearest point to the mouse point.
%Button 2 Adds a point along the neares segment at the mouse point.
%Called by contourControls.
%JOD.
%LM:  2 Dec 02, JOD.
%
% 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/>.

figure(hCERR)

hold on
[x, y, button]  = ginput(1);

if isempty(button)
  xV = [];
  yV = [];
  return
end

hPlot = getActiveSegment;
segS = get(hPlot,'userdata');
xV = segS.GUI_xV;
yV = segS.GUI_yV;

switch button

  case 1 %replace nearest point

    rV = ((x - xV).^2 + (y - yV).^2 ).^0.5;
    indV = find(rV == min(rV));
    ind = indV(1);
    xV(ind) = x;
    yV(ind) = y;
    delete(hPlot)
    hPlot = plot([xV,xV(1)],[yV,yV(1)],'rs-');  %plot closed contour
    set(hPlot,'markersize',4)

  case 3 %add point along nearest edge

    %Cycle through all the line segments just to be sure we get the nearest perpindicular
    rV = ((x - xV).^2 + (y - yV).^2 ).^0.5;
    indV = find(rV == min(rV));
    ind0 = indV(1);
    r0V = [x, y];
    cycleLength = length(xV);
    nearest = -1;
    bestDist = 999999;
    ind_a = [];
    for i = 1 : cycleLength;
      indTest1 = ind0 + i - 1;
      [ind1] = cycle(indTest1, cycleLength);
      indTest2 = ind0 + i;
      [ind2] = cycle(indTest2, cycleLength);
      r1V = [xV(ind1),yV(ind1)];
      r2V = [xV(ind2),yV(ind2)];
      %Does the perpendicular hit on the line segment?
      [dist, flag1] = dist2seg(r0V, r1V, r2V);
      if flag1 == 0 & dist < bestDist
        ind_a = ind1;
        ind_b = ind2;
        bestDist = dist;
        nearest = i;
      end
    end

    %Compare with distances to nearest neighbors
    indTest = ind0 -1;
    cycleLength = length(xV);
    [ind1] = cycle(indTest, cycleLength);
    indTest = ind0 + 1;
    [ind2] = cycle(indTest, cycleLength);
    r1 = ((x - xV(ind1)).^2 + (y - yV(ind1)).^2 )^0.5;
    r2 = ((x - xV(ind2)).^2 + (y - yV(ind2)).^2 )^0.5;
    if  (r1 < r2 & r1 < bestDist)
        ind_a = ind0;
        ind_b = ind1;
    elseif r2 < bestDist
        ind_a = ind0;
        ind_b = ind2;
    end

    %Plug the new point in:
    if abs(ind_a - ind_b) == 1 %no end-wrapping
      insert1 = min(ind_a,ind_b);
      xV = [xV(1:insert1), x, xV(insert1+1:end)];
      yV = [yV(1:insert1), y, yV(insert1+1:end)];
    else  %case of end-wrapping
      xV = [xV, x];
      yV = [yV, y];
    end

    %redraw
    delete(hPlot)
    hPlot = plot([xV,xV(1)],[yV,yV(1)],'rs-');  %Only plot with same first and last point; store no duplicates
    set(hPlot,'markersize',4)

end

hold off

%segmentsS contains all the information about the plotted segments - it gets written into
%the userdata portion of any segment.

segS.active = 1;
stateS.activeSegment = segS.segNum;
segS.hPlot  = hPlot;
segS.GUI_xV = xV;
segS.GUI_yV = yV;

set(hPlot,'tag','newSegment')
set(hPlot,'userdata',segS)


text:
function [ projections ] = Ax(img, geo, angles, varargin )
%AX computes projections for images and geometry information

%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/blob/master/LICENSE
% Coded by:           Ander Biguri
%--------------------------------------------------------------------------
% Lets make 100% sure that data is correct
%% Optionals

ptype='Siddon';
% expected projection types. 'ray-voxel' is obsolete. Use 'Siddon'
expectedProjectionTypes = {'Siddon','ray-voxel','interpolated'};
acceptableOptionName = {'gpuids'};

if nargin > 3
    if any(strcmp(varargin{1}, expectedProjectionTypes))
        ptype = varargin{1};
        [gpuids] = parse_inputs(varargin{2:length(varargin)});
        %[ptype, gpuids] = parse_inputs1(varargin{1}, expectedProjectionTypes, varargin{2:length(varargin)});
    elseif any(strcmp(varargin{1}, acceptableOptionName))
        [gpuids] = parse_inputs(varargin{:});
        %[ptype, gpuids] = parse_inputs1(ptype, expectedProjectionTypes, varargin{:});
    else
        assert(false,'TIGRE:Ax:InvalidInput','Projection type should be either ''interpolated'' or ''Siddon''.');
    end
else
    gpuids = GpuIds();
end

%% image
assert(isa(img,'single'),'TIGRE:Ax:InvalidInput','Image should be single type');
assert(isreal(img),'TIGRE:Ax:InvalidInput','Image should be real (non-complex)');
% assert(size(img,3)>1,'TIGRE:Ax:InvalidInput', 'image should be 3D'); %TODO: needed?
%% Angles
assert(isreal(angles),'TIGRE:Ax:InvalidInput','Angles should be real (non-complex)');
assert(size(angles,1)==1 | size(angles,1)==3 ,'TIGRE:Ax:InvalidInput','Angles should be of size 1xN or 3xN');
angles=double(angles); %in case they were single.
if size(angles,1)==1
   angles=repmat(angles,[3 1]);
   angles(2,:)=0;
   angles(3,:)=0;
end
%% geometry
geo=checkGeo(geo,angles);
assert(isequal([size(img,1) size(img,2) size(img,3)],squeeze(geo.nVoxel.')),'TIGRE:Ax:BadGeometry','nVoxel does not match with provided image size');

%% Temporary (?)

s= sum(abs(angles),2);
if (s(2)~=0 || s(3)~=0) && (strcmp(ptype,'Siddon')||strcmp(ptype,'ray-voxel')) && strcmp(geo.mode,'parallel')
    warning(['''Siddon'' Not supported for parallel beam arbitrary axis rotation, changed to ''interpolated''.',char(10),'Ignore this message if you are not purposedly triying enforce its use.']);
    ptype='interpolated';
end

%% Thats it, lets call the mex fucntion
%% TODO: When Ax_mex accepts class-objects, gpuids itself will be passed.
projections=Ax_mex(img,geo,angles,ptype, gpuids.devices);

end

function [gpuids]=parse_inputs(varargin)
    %fprintf('parse_inputs0(varargin (%d))\n', length(varargin));
    if isempty(varargin)
        gpuids = GpuIds();
    else
        % create input parser
        p=inputParser;
        % add optional parameters
        addParameter(p,'gpuids', GpuIds());
        %execute
        parse(p,varargin{:});
        %extract
        gpuids=p.Results.gpuids;
    end
end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function newnode=sms(node,face,iter,alpha,method)
%
% newnode=sms(node,face,iter,useralpha,method)
%
% simplified version of surface mesh smoothing
%
% author: Qianqian Fang, <q.fang at neu.edu>
% date: 2009/10/21
%
% input:
%    node:  node coordinates of a surface mesh
%    face:  face element list of the surface mesh
%    iter:  smoothing iteration number
%    alpha: scaler, smoothing parameter, v(k+1)=alpha*v(k)+(1-alpha)*mean(neighbors)
%    method: same as in smoothsurf, default is 'laplacianhc'
%
% output:
%    newnode: output, the smoothed node coordinates
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

if(nargin<5)
   method='laplacianhc';
end
if(nargin<4)
   if(nargin<3)
      iter=10;
   end
   alpha=0.5;
end

conn=meshconn(face,size(node,1));
newnode=smoothsurf(node(:,1:3),[],conn,iter,alpha,method,alpha);


text:
function LL = dirichlet_score_family(counts, prior)
% DIRICHLET_SCORE Compute the log marginal likelihood of a single family
% LL = dirichlet_score(counts, prior)
%
% counts(a, b, ..., z) is the number of times parent 1 = a, parent 2 = b, ..., child = z
% prior is an optional multidimensional array of the same shape as counts.
% It defaults to a uniform prior.
%
% We marginalize out the parameters:
% LL = log \int \prod_m P(x(i,m) | x(Pa_i,m), theta_i) P(theta_i) d(theta_i)

% LL = log[  prod_j gamma(alpha_ij)/gamma(alpha_ij + N_ij)  *
%            prod_k gamma(alpha_ijk + N_ijk)/gamma(alpha_ijk)  ]
% Call the prod_k term U and the prod_j term  V.
% We reshape all quantities into (j,k) matrices
% This formula was first derived by Cooper and Herskovits, 1992.
% See also "Learning Bayesian Networks", Heckerman, Geiger and Chickering, MLJ 95.

ns = mysize(counts);
ns_ps = ns(1:end-1);
ns_self = ns(end);

if nargin < 2, prior = normalise(myones(ns)); end

if 1
  prior = reshape(prior(:), [prod(ns_ps) ns_self]);
  counts = reshape(counts,  [prod(ns_ps) ns_self]);
  %U = prod(gamma(prior + counts) ./ gamma(prior), 2); % mult over k
  LU = sum(gammaln(prior + counts) - gammaln(prior), 2);
  alpha_ij = sum(prior, 2); % sum over k
  N_ij = sum(counts, 2);
  %V = gamma(alpha_ij) ./ gamma(alpha_ij + N_ij);
  LV = gammaln(alpha_ij) - gammaln(alpha_ij + N_ij);
  %L = prod(U .* V);
  LL = sum(LU + LV);
else
  CPT = mk_stochastic(prior + counts);
  LL = sum(log(CPT(:) .* counts(:)));
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [filt] = bandstopfilter(dat,Fs,Fbp,N,type,dir)

% BANDSTOPFILTER filters EEG/MEG data in a specified band
%
% Use as
%   [filt] = bandstopfilter(dat, Fsample, Fbp, N, type, dir)
% where
%   dat        data matrix (Nchans X Ntime)
%   Fsample    sampling frequency in Hz
%   Fbp        frequency band, specified as [Fhp Flp]
%   N          optional filter order, default is 4 (but) or 25 (fir)
%   type       optional filter type, can be
%                'but' Butterworth IIR filter (default)
%                'fir' FIR filter using MATLAB fir1 function
%   dir        optional filter direction, can be
%                'onepass'         forward filter only
%                'onepass-reverse' reverse filter only, i.e. backward in time
%                'twopass'         zero-phase forward and reverse filter (default)
%
% Note that a one- or two-pass filter has consequences for the
% strength of the filter, i.e. a two-pass filter with the same filter
% order will attenuate the signal twice as strong.
%
% See also LOWPASSFILTER, HIGHPASSFILTER, BANDPASSFILTER

% Copyright (c) 2007, 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$

% set the default filter order later
if nargin<4
    N = [];
end

% set the default filter type
if nargin<5
  type = 'but';
end

% set the default filter direction
if nargin<6
  dir = 'twopass';
end

% Nyquist frequency
Fn = Fs/2;

% compute filter coefficients
switch type
  case 'but'
    if isempty(N)
      N = 4;
    end
    [B, A] = butter(N, [min(Fbp)/Fn max(Fbp)/Fn], 'stop');
  case 'fir'
    if isempty(N)
      N = 25;
    end
    B = fir1(N, [min(Fbp)/Fn max(Fbp)/Fn], 'stop');
    A = 1;
end

% apply filter to the data
switch dir
  case 'onepass'
    filt = filter(B, A, dat')';
  case 'onepass-reverse'
    dat  = fliplr(dat);
    filt = filter(B, A, dat')';
    filt = fliplr(filt);
  case 'twopass'
    filt = filtfilt(B, A, dat')';
end


text:
function varargout = fillSphericalPolygon(sphere, poly, germ)
%FILLSPHERICALPOLYGON  Fill a spherical polygon
%
%   fillSphericalPolygon(SPHERE, POLY, GERM)
%
%
%   Example
%   fillSphericalPolygon
%
%   See also
%   drawSphericalPolygon, fillSphericalTriangle, drawSphere
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2012-02-09,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2012 INRA - Cepia Software Platform.

nv = size(poly, 1);

h = zeros(nv, 1);
for i = 1:nv
    v1 = poly(i, :);
    v2 = poly(mod(i, nv) + 1, :);

    h(i) = fillSphericalTriangle(sphere, germ, v1, v2);
end

if nargout > 0
    varargout = {h};
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function fadingcoeff=genh(I,v,Dop,tb)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                        %
%%     Name: genh.m                                                       %
%%                                                                        %
%%     Description: We generate a unique coefficient of fading with "sum  %
%%      of sinusoides" of the Jakes model.                                %
%%                                                                        %
%%     Parameters:                                                        %
%%          I = Length of the plot                                        %
%%          v = Speed of the terminal (m/s)                               %
%%          Dop = Maximum frequency of the Doppler effect                 %
%%          tb = Symbol Duration                                          %
%%                                                                        %
%%     Authors: Bertrand Muquet, Sebastien Simoens, Shengli Zhou          %
%%      October 2000                                                      %
%%                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  fc = 2.3e9; 			% Carrier Frequency in Hertz (2.5GHz 3.2GHz)
  fdmax = Dop;

  N = 100;                      % Number of incident waves
  t = tb:tb:tb*I;               % The variable "time"

  len = length(t);
  theta = rand(1,N)*2*pi;              % Generating the uniform phases
  fd = cos(2*pi*((1:N)/N))*fdmax;      % Generate eqaul-spaced frequencies from  "-fdmax" to "+fdmax"

  E = exp(j.*(2*pi*fd(:)*t(:)'+repmat(theta(:),1,len)));
  E = E/sqrt(N);
  fadingcoeff = sum(E);
  %plot(t,abs(fadingcoeff))
  %xlabel('time (second)');ylabel('Envelope of the fading coefficient');


text:
function [ position, velocity, acceleration ] = interpolate( s, nrl, target, t )
% This function differentiates and interpolates a set of Chebyshev coefficients to give position and velocity.
%
% void interpolate (double *buf, double *t, long int ncf, long int na,
%
%                   double *position, double *velocity)
%
% ------------------------------------------------------------------------
%
%    PURPOSE:
%       This function differentiates and interpolates a set of
%       Chebyshev coefficients to give position and velocity.
%
%    REFERENCES:
%       Standish, E.M. and Newhall, X X (1988). "The JPL Export
%          Planetary Ephemeris"; JPL document dated 17 June 1988.
%
%    INPUT
%    ARGUMENTS:
%       *buf (double)
%          Array of Chebyshev coefficients of position.
%       *t (double)
%          t[0] is fractional time interval covered by coefficients at
%          which interpolation is desired (0 <= t[0] <= 1).
%          t[1] is length of whole interval in input time units.
%       ncf (long int)
%          Number of coefficients per component.
%       na (long int)
%          Number of sets of coefficients in full array
%          (i.e., number of sub-intervals in full interval).
%
%    OUTPUT
%    ARGUMENTS:
%       *position (double)
%          Position array of requested object.
%       *velocity (double)
%          Velocity array of requested object.
%
%    RETURNED
%    VALUE:
%       None.
%
%    GLOBALS
%    USED:
%       NP                eph_manager.h
%       NV                eph_manager.h
%       PC                eph_manager.h
%       VC                eph_manager.h
%       TWOT              eph_manager.h
%
%    FUNCTIONS
%    CALLED:
%       fmod              math.h
%
%    VER./DATE/
%    PROGRAMMER:
%       V1.0/03-93/WTH (USNO/AA): Convert FORTRAN to C.
%       V1.1/07-93/WTH (USNO/AA): Update to C standards.
%       V1.2/07-98/WTH (USNO/AA): Modify to make position and velocity
%                                 two distinct vector arrays.
%       V1.3/11-07/WKP (USNO/AA): Updated prolog.
%       V1.4/12-07/WKP (USNO/AA): Changed ncf and na arguments from short
%                                 int to long int.
%       V1.5/10-10/WKP (USNO/AA): Renamed function to lowercase to
%                                 comply with coding standards.
%
%    NOTES:
%       None.
%
% ------------------------------------------------------------------------

  persistent PC VC AC NP TWOT;
  if isempty(PC)
    PC = zeros(1,15); % ncf <= 15
    VC = zeros(1,15); % ncf <= 15
    AC = zeros(1,15); % ncf <= 15
    TWOT = 0.0;
    PC(1) = 1.0; % T_0(x) = 1
    PC(2) = 0.0; % T_1(x) = x
    VC(1) = 0.0; % T'_0(x) = 0
    VC(2) = 1.0; % T'_1(x) = 1
    AC(1) = 0.0; % T"_0(x) = 0
    AC(2) = 0.0; % T"_1(x) = 0
    AC(3) = 4.0; % T"_2(x) = 4
    NP = 2;
  end
  if target == Ephem.TTState
    NDIM = 1;
  elseif target == Ephem.NutationsState
    NDIM = 2;
  else
    NDIM = 3;
  end

  buf = s.IPT(1,target) + nrl*s.indexInc + s.indexOffset;
  ncf = s.IPT(2,target);
  na  = s.IPT(3,target);
  table = s.scan;

  position = zeros(1,NDIM); % row vector
  velocity = zeros(1,NDIM); % row vector
  acceleration = zeros(1,NDIM); % row vector
  if ncf <= 0
    if s.output > 0
      fprintf(s.output,'ERROR: interpolate DE%s. No coefficients available for target %d\n', ...
        s.de_number,target);
    end
    return;
  end
  %   Get correct sub-interval number for this set of coefficients and
  %   then get normalized Chebyshev time within that subinterval.
  dna = double(na);
  temp = dna * double(t);

  %   'tc' is the normalized Chebyshev time (-1 <= tc <= 1).

  tc = 2.0 * mod(temp,1.0) - 1.0;

  %   Check to see whether Chebyshev time has changed, and compute new
  %   polynomial values if it has.  (The element PC[1] is the value of
  %   t1[tc] and hence contains the value of 'tc' on the previous call.)

  if tc ~= PC(2)
    NP = 3;
    TWOT = tc + tc;
    PC(2) = tc; % T_1(x) = x
    PC(3) = TWOT * PC(2) - PC(1); % T_2(x) = 2x^2-1
    VC(3) = 2.0 * TWOT; % T'_2(x) = 4x
  end

  %   Be sure that at least 'ncf' polynomials have been evaluated and
  %   are stored in the array 'PC'.
  %   Chebyshev polynomial recurrence relationship PC[i+1] = T_i(tc)

  if NP < ncf
    for i = (NP+1):ncf
      % T_i-1(x) = 2x T_i-2(x) - T_i-3(x)
      PC(i) = TWOT * PC(i-1) - PC(i-2);
      % T'_i-1(x) = 2x T'_i-2(x) + 2 T_i-2(x) - T'_i-3(x)
      VC(i) = TWOT * VC(i-1) + PC(i-1) + PC(i-1) - VC(i-2);
      % T"_i-1(x) = 2x T"_i-2(x) + 4 T'_i-2(x) - T"_i-3(x)
      AC(i) = TWOT * AC(i-1) + 4*VC(i-1) - AC(i-2);
    end
    NP = ncf;
  end

  vfac = 2.0 * dna;
  %   Interpolate to get position for each component.

  k = buf + int32(floor(temp)) * (NDIM * ncf);
  if k + NDIM*ncf > length(table)
    if s.output > 0
      fprintf(s.output,'interpolate: past end of coefficients by %d/%d\n',...
        k + NDIM*ncf-length(table),length(table));
    end
    return;
  end
  % p = sum( a_i * T_i(x) )
  % v = sum( a_i * T'_i(x) )*2*n
  % a = sum( a_i * T"_i(x) )*4*n^2
  for i = 1:NDIM
    p = 0.0;
    v = 0.0;
    a = 0.0;
    for j = ncf:-1:1
      p = p + PC(j) * table(k + j);
      v = v + VC(j) * table(k + j);
      a = a + AC(j) * table(k + j);
    end
    position(i) = p;
    velocity(i) = v*vfac;
    acceleration(i) = a*vfac*vfac;
    k = k + ncf;
  end

end

author: goGPS-Project


text:
function [A, B] = reflexKronApprox( K, m, n, Rtol )
%
%       [A, B] = reflexKronApprox( K, m, n, Rtol );
%
% computes a kronecker sum approximation to the blurring matrix K
% that arises from the input PSF under reflexive boundary conditions.
% This approximation is done a-la Nagy-Ng-Perrone (see paper for details).
%
%  Input:
%         K - psfMatrix object
%         m - size of matrix A (assumed square)
%         n - size of matrix B (assumed square)
%         Rtol - relative tolerance ( <=1 ); determines number of kron products
%                in the approximation.  If si/s1 > Rtol,
%                where si is the ith largest
%                singular value of the weighted PSF,
%                then Ai(x)Bi will be included in the approximation.
%                IF ONLY ONE TERM IN THE SUM IS DESIRED, USE Rtol = 1.
%
%  Output:
%      Cells A and B such that K \approx \sum_i[ A{i} \otimes B{i} ].
%

%  L. Perrone, 4/28/02

%  Modifications:
%  5/25/02, J. Nagy
%           Cosmetic changes to incorporate into RestoreTools
%
%  9/??/02  L. Perrone
%           This code now conforms to the new kronMatrix class,
%           where fields K.a and K.b are cell arrays containing
%           (possibly) more than one matrix, depending on Rtol.
%  11/22/02 L. Perrone
%           The kronMatrix class should work for image processing
%           problems, where it is common to use lexicographical (row)
%           ordering.  However, the kronMatrix class was designed
%           using vec(column) ordering (unlike all other classes
%           in the RestoreTools package).  Until now, the inconsistency
%           remained hidden because the PSFs in use
%           had been symmetric or close to symmetric.
%           This discrepancy has now been fixed.

P1 = K.psf;
P2 = P1.image;
PSF = P2{1};
c1 = P1.center;
center = c1{1};

[mp, np] = size(PSF);

if ( mp ~= np )
  error('For now, we expect PSF to be square')
end

%
% Compute weighted PSF.
%
c = zeros(mp,1);
c(1) = mp;
c(2:2:end) = 1;
R = chol( toeplitz(c) );

Phat = R*PSF*R';

%
% Compute SVD of weighted PSF, which is then used to construct
% the separable approximation.
%
[U,S,V] = svd( Phat );

%
% check to make sure first column looks like
% a Gaussian, and is not inverted.
%
minU = abs(min(min(U(:,1))));
maxU = max(max(abs(U(:,1))));
if minU == maxU
  U = -U;
  V = -V;
end

%
% Construct approximation.
%
A = cell(1);
B = cell(1);
i=1;
PP = zeros(size(PSF));
while S(i,i)/S(1,1) >= Rtol
  a = R \ ( U(:,i) * sqrt(S(i,i)) );
  b = R \ ( V(:,i) * sqrt(S(i,i)) );
  PP = PP + a*b';
  % Comment out the next two statements, replace with the two that
  % follow, in order to fix the 11/22/02 problem.
  %A{i} = build_toep(a, center(1), m) + buildHank(a, center(1), m);
  %B{i} = build_toep(b, center(2), n) + buildHank(b, center(2), n);
  B{i} = build_toep(a, center(1), m) + buildHank(a, center(1), m);
  A{i} = build_toep(b, center(2), n) + buildHank(b, center(2), n);
  i=i+1;
end

%
%  The following scales the approximation so that the
%  corresponding PSF approximation has the property
%  that its sum of values is the same as the sum of the
%  original PSF values.
%
cs = sqrt(sum(PSF(:))/sum(PP(:)));
for i = 1:length(A)
  A{i} = cs*A{i};
  B{i} = cs*B{i};
  AA = A{i};
  BB = B{i};
end

repo: IRtools


text:
function show_sites(sites)
                 for ccc=1:numel(sites)
                     debug.print_debug(1,sites(ccc).channeltag.string())
                 end
end

author: geoscience-community-codes

repo: GISMO


text:
function display_legend(plot_obj)
% Displays the legend on the plot.  Call this function after all
% curves and points have been plotted.  The only function that
% should be called after this one is the one for saving the plot.
% Legend entries for curves are stored with add_legend_entry.

if length(plot_obj.legend_strings) > 0
    figure(plot_obj.fh)
    legend(plot_obj.handles_vec,plot_obj.legend_strings)
end
end

author: nesl


text:
function sigma = optimize_wrt_sigma(I, I1, options)
% set the left hand side (warp the hyperspectral image)
options.isRigid = 0;
I1 = transform(I1,eye(3),[1,1],1e-3,size(I1,2),size(I1,1),options);
options.isRigid = 1;

% configure the right hand side (transform the color image)
sigma_step = mean(options.s);
step_size = sigma_step;

pa = ParameterAnalysis();
pa.StepSizeMode = 2;
pa.NeighborhoodMode = 3;
pa.Algorithm = 'BrutalForce';
pa.BrutalForceDepth = 5;
list_params = {'sigma'};
range = [1e-6,4*mean(options.s)]';
fcn_run = @(options1) calc_val(options1,I,I1,options);

options1 = [];
options1.sigma = options.sigma;

options1 = pa.autoParamSelection(fcn_run, options1, list_params, ...
    step_size, range);
sigma = options1.sigma;

function val = calc_val(options1,I,I1,options)
sigma = options1.sigma;
T = create_T(options.degree, options.t);
val = eval_obj_fun(I, I1, T, options.s, sigma, options);

author: zhouyuanzxcv

repo: Hyperspectral


text:
function P = get_filename(dwcard,wcard,varargin)
% :Usage:
% ::
%
%     function P = get_filename(dir_wcard,file_wcard,[verbose])
%
% Start in directory above individual subject directories
%
% :Inputs:
%
%   **dwcard:**
%        Enter dwcard for wildcard of directories to look in; e.g. '02*'
%
%   **wcard:**
%        Enter wcard for image or file to get - e.g., 'beta_0001.img'
%        This can also include subdirectories (no *'s) 'anatomy/beta*img'
%
% :Output:
%
%     Returns list of all files in individual directories  in string matrix
%
% Missing files, or entries in directory that do not contain files, are
% removed from the list.
%
% NOT entering a * in dwcard seems to produce an error.
%
% :Examples:
% ::
%
%    P = get_filename('02*','beta_0001*')
%    P = get_filename('02*','beta_0001.img')
%    P = get_filename('02*','anatomy/nnhe*_seg1.img')
%    P = get_filename('020515sp*','anatomy/nnhe*_seg1.img')
%
% one * is allowed in first string, multiple *'s in second,
% as long as they are in the filename, not directory names!
%
% ..
%    Tor Wager 10/8/02
% ..

if length(varargin) > 0, verb = varargin{1};, else, verb = 0;,end

clear d, clear D
d = dir(dwcard); dind = 1;,
for i = 1:length(d),
	if d(i).isdir,
		D{dind} = d(i).name;,
		dind = dind + 1;,
	end,
end

[ld f e] = fileparts(wcard);
if isempty(ld)
	P = tor_list_files(D,wcard);

else
	for i = 1:length(D)
		D{i} = [D{i} filesep ld];
	end
	P = tor_list_files(D,[f e]);
end
P = str2mat(P);

for i = 1:size(P,1),
	if isempty(deblank(P(i,:))),
		%warning(['No file found for ' D{i}]),
		pp(i) = 1;,
	else,
		pp(i) = 0;,
	end
end

for i = 1:size(P,1)
    if isempty(deblank(P(i,:))), pp(i) = 1;, end
end

P = P(~pp,:);

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function van_der_corput_test02 ( )

%*****************************************************************************80
%
%% VAN_DER_CORPUT_TEST02 tests VAN_DER_CORPUT_SEQUENCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 February 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'VAN_DER_CORPUT_TEST02\n' );
  fprintf ( 1, '  VAN_DER_CORPUT_SEQUENCE returns N elements\n' );
  fprintf ( 1, '  of a van der Corput sequence in the current base.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I   R\n' );
  fprintf ( 1, '\n' );

  seed = 0;
  van_der_corput_seed_set ( seed );

  base = 2;
  van_der_corput_base_set ( base );

  n = 10;
  r = van_der_corput_sequence ( n );

  for i = 1: 10
    fprintf ( 1, '%d %f\n', seed+i-1, r(i) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function extra_data_info=my_net_forward(net_info, work_info_batch, data_info, net_run_config, extra_output_layer_idxes)

    gpu_mode=net_run_config.use_gpu;
    if gpu_mode
      if ~net_info.ref.net_on_gpu
           my_move_net(net_info, 'gpu') ;
      end
      data_info.ref.output_info_layers{1}=...
          move_output_info_gpu(data_info.ref.output_info_layers{1});
    else
      assert(~net_info.ref.net_on_gpu);
      data_info.ref.output_info_layers{1}=...
          move_output_info_cpu(data_info.ref.output_info_layers{1});
    end

    extra_data_info=do_forward(net_info, work_info_batch, data_info, net_run_config, extra_output_layer_idxes);

    if gpu_mode
        if ~net_info.ref.net_stay_on_gpu
            my_move_net(net_info, 'cpu') ;
        end
    end

end

function extra_data_info=do_forward(net_info, work_info_batch, data_info, net_run_config, extra_output_layer_idxes)

layer_num = numel(net_info.ref.layers) ;
gpu_mode = net_run_config.use_gpu;

if gpu_mode && net_run_config.sync
    wait(gpuDevice) ;
end

assert(check_valid_net_output(data_info.ref.output_info_layers{1}));

bp_turn_on_layer=-1;
disableDropout=true;
keep_layer_output=false;

do_bp_current_net=check_do_bp_current_net(work_info_batch, net_run_config, net_info);
if do_bp_current_net
    bp_start_layer=net_info.ref.bp_start_layer;
    bp_turn_on_layer=bp_start_layer;
end

data_info.ref.need_bp=do_bp_current_net;

extra_output_layer_flags=false(layer_num, 1);
if ~isempty(extra_output_layer_idxes)
    extra_output_layers=cell(layer_num, 1);
    extra_output_layer_flags(extra_output_layer_idxes)=true;
else
    extra_output_layers=[];
end

for layer_idx=1:layer_num

    input_info=data_info.ref.output_info_layers{layer_idx};
    l = net_info.ref.layers{layer_idx} ;
    is_simple_layer= ~strcmp(l.type, 'my_custom');

    if bp_turn_on_layer==layer_idx
        disableDropout=false;
        keep_layer_output=true;
    end

        if is_simple_layer
            output_info=data_info.ref.output_info_layers{layer_idx+1};
            assert(~input_info.is_group_data);
            assert(~output_info.is_group_data);
            switch l.type
                case 'conv'
                    input_size=size(input_info.x);
                    filter_size=size(l.filters);
                    if filter_size(3)~=input_size(3)

                        disp('filter_size:');
                        disp(filter_size);
                        disp('input_size:');
                        disp(input_size);
                        error('filter size not match input size!');
                    end
                    if any(filter_size(1:2)>input_size(1:2))
                        error('filter size larger than the input size!');
                    end
                  output_info.x = vl_nnconv(input_info.x, l.filters, l.biases, 'pad', l.pad, 'stride', l.stride) ;
                case 'pool'
                    input_size=size(input_info.x);
                    pool_size=l.pool;
                    if any(pool_size(1:2)>input_size(1:2))
                        error('pool size larger than the input size!');
                    end
                    output_info.x = vl_nnpool(input_info.x, l.pool, 'pad', l.pad, 'stride', l.stride, 'method', l.method) ;
                case 'normalize'
                  output_info.x = vl_nnnormalize(input_info.x, l.param) ;
                case 'softmax'
                  output_info.x = vl_nnsoftmax(input_info.x) ;
                case 'relu'
                  output_info.x = vl_nnrelu(input_info.x) ;
                case 'noffset'
                  output_info.x = vl_nnnoffset(input_info.x, l.param) ;
                case 'dropout'

                  if disableDropout
                    output_info.x = input_info.x ;
                  else
                    [output_info.x, output_info.aux] = vl_nndropout(input_info.x, 'rate', l.rate) ;
                  end

                otherwise
                    error('Unknown layer type %s', l.type) ;
            end
        else

            output_info= l.forward_fn(input_info, l, work_info_batch) ;

        end

        if isempty(output_info)
            break;
        end
        output_info.forward_finished=true;

    if ~keep_layer_output && layer_idx < layer_num - 1
      data_info.ref.output_info_layers{layer_idx}=[];
      input_info=[];
    end

    if gpu_mode
        if ~net_info.ref.data_stay_on_gpu
            if ~isempty(input_info)
                input_info=move_output_info_cpu(input_info);
                data_info.ref.output_info_layers{layer_idx}=input_info;
            end
        end
    end

    if gpu_mode && net_run_config.sync
        wait(gpuDevice) ;
    end

    data_info.ref.output_info_layers{layer_idx+1}=output_info;

    if extra_output_layer_flags(layer_idx)
        extra_output_layers{layer_idx}=output_info;
    end

end

if gpu_mode && ~net_info.ref.data_stay_on_gpu
    output_info=data_info.ref.output_info_layers{end};
    output_info=move_output_info_cpu(output_info);
    data_info.ref.output_info_layers{end}=output_info;
end

extra_data_info=[];
extra_data_info.output_layers=extra_output_layers;

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function [models,M,test_set] = esvm_download_models(cls)
% Download a pre-trained PASCAL VOC2007 model from my MIT CSAIL
% homepage (http://people.csail.mit.edu/tomasz/)
%
% Copyright (C) 2011-12 by Tomasz Malisiewicz
% All rights reserved.
%
% This file is part of the Exemplar-SVM library and is made
% available under the terms of the MIT license (see COPYING file).
% Project homepage: https://github.com/quantombone/exemplarsvm

f = sprintf('%s.mat',cls);
models = [];
M = [];
test_set = {};
if ~exist(f,'file')
  unix(sprintf('wget http://people.csail.mit.edu/tomasz/exemplarsvm/models/%s',f));
else
  fprintf(1,'Found %s, not downloading\n',f);
end

if nargout == 0
  return;
end

load(f);

author: quantombone

repo: exemplarsvm


text:
function [P,Uinit,output] = cp_als(X,R,varargin)
%CP_ALS Compute a CP decomposition of any type of tensor.
%
%   P = CP_ALS(X,R) computes an estimate of the best rank-R
%   CP model of a tensor X using an alternating least-squares
%   algorithm.  The input X can be a tensor, sptensor, ktensor, or
%   ttensor. The result P is a ktensor.
%
%   P = CP_ALS(X,R,'param',value,...) specifies optional parameters and
%   values. Valid parameters and their default values are:
%      'tol' - Tolerance on difference in fit {1.0e-4}
%      'maxiters' - Maximum number of iterations {50}
%      'dimorder' - Order to loop through dimensions {1:ndims(A)}
%      'init' - Initial guess [{'random'}|'nvecs'|cell array]
%      'printitn' - Print fit every n iterations; 0 for no printing {1}
%
%   [P,U0] = CP_ALS(...) also returns the initial guess.
%
%   [P,U0,out] = CP_ALS(...) also returns additional output that contains
%   the input parameters.
%
%   Note: The "fit" is defined as 1 - norm(X-full(P))/norm(X) and is
%   loosely the proportion of the data described by the CP model, i.e., a
%   fit of 1 is perfect.
%
%   Examples:
%   X = sptenrand([5 4 3], 10);
%   P = cp_als(X,2);
%   P = cp_als(X,2,'dimorder',[3 2 1]);
%   P = cp_als(X,2,'dimorder',[3 2 1],'init','nvecs');
%   U0 = {rand(5,2),rand(4,2),[]}; %<-- Initial guess for factors of P
%   [P,U0,out] = cp_als(X,2,'dimorder',[3 2 1],'init',U0);
%   P = cp_als(X,2,out.params); %<-- Same params as previous run
%
%   See also KTENSOR, TENSOR, SPTENSOR, TTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Extract number of dimensions and norm of X.
N = ndims(X);
normX = norm(X);

%% Set algorithm parameters from input or by using defaults
params = inputParser;
params.addParamValue('tol',1e-4,@isscalar);
params.addParamValue('maxiters',50,@(x) isscalar(x) & x > 0);
params.addParamValue('dimorder',1:N,@(x) isequal(sort(x),1:N));
params.addParamValue('init', 'random', @(x) (iscell(x) || ismember(x,{'random','nvecs'})));
params.addParamValue('printitn',1,@isscalar);
params.parse(varargin{:});

%% Copy from params object
fitchangetol = params.Results.tol;
maxiters = params.Results.maxiters;
dimorder = params.Results.dimorder;
init = params.Results.init;
printitn = params.Results.printitn;

%% Error checking

%% Set up and error checking on initial guess for U.
if iscell(init)
    Uinit = init;
    if numel(Uinit) ~= N
        error('OPTS.init does not have %d cells',N);
    end
    for n = dimorder(2:end);
        if ~isequal(size(Uinit{n}),[size(X,n) R])
            error('OPTS.init{%d} is the wrong size',n);
        end
    end
else
    % Observe that we don't need to calculate an initial guess for the
    % first index in dimorder because that will be solved for in the first
    % inner iteration.
    if strcmp(init,'random')
        Uinit = cell(N,1);
        for n = dimorder(2:end)
            Uinit{n} = rand(size(X,n),R);
        end
    elseif strcmp(init,'nvecs') || strcmp(init,'eigs')
        Uinit = cell(N,1);
        for n = dimorder(2:end)
            Uinit{n} = nvecs(X,n,R);
        end
    else
        error('The selected initialization method is not supported');
    end
end

%% Set up for iterations - initializing U and the fit.
U = Uinit;
fit = 0;

if printitn>0
  fprintf('\nCP_ALS:\n');
end

%% Main Loop: Iterate until convergence

if (isa(X,'sptensor') || isa(X,'tensor')) && (exist('cpals_core','file') == 3)

    %fprintf('Using C++ code\n');
    [lambda,U] = cpals_core(X, Uinit, fitchangetol, maxiters, dimorder);
    P = ktensor(lambda,U);

else

    for iter = 1:maxiters

        fitold = fit;

        % Iterate over all N modes of the tensor
        for n = dimorder(1:end)

            % Calculate Unew = X_(n) * khatrirao(all U except n, 'r').
            Unew = mttkrp(X,U,n);

            % Compute the matrix of coefficients for linear system
            Y = ones(R,R);
            for i = [1:n-1,n+1:N]
                Y = Y .* (U{i}'*U{i});
            end

            % Need to figure out which of the following lines works better.
            % Meanwhile, we'll stick witht he line from TTB 2.2 since that is
            % what seems to work best based on preliminary testing.
            %Unew = Unew * pinv(Y); %<- Line from TTB 2.3.
            Unew = (Y \ Unew')'; %<- Line from TTB 2.2.

            % Normalize each vector to prevent singularities in coefmatrix
            if iter == 1
                lambda = sqrt(sum(Unew.^2,1))'; %2-norm
            else
                lambda = max( max(Unew,[],1), 1 )'; %max-norm
            end
            Unew = Unew * spdiags(1./lambda,0,R,R);
            if issparse(Unew)
                U{n} = full(Unew);   % for the case R=1
            else
                U{n} = Unew;
            end
        end

        P = ktensor(lambda,U);
        normresidual = sqrt( normX^2 + norm(P)^2 - 2 * innerprod(X,P) );
        fit = 1 - (normresidual / normX); %fraction explained by model
        fitchange = abs(fitold - fit);

        if mod(iter,printitn)==0
            fprintf(' Iter %2d: fit = %e fitdelta = %7.1e\n', iter, fit, fitchange);
        end

        % Check for convergence
        if (iter > 1) && (fitchange < fitchangetol)
            break;
        end
    end
end

%% Clean up final result
% Arrange the final tensor so that the columns are normalized.
P = arrange(P);
% Fix the signs
P = fixsigns(P);

if printitn>0
  normresidual = sqrt( normX^2 + norm(P)^2 - 2 * innerprod(X,P) );
  fit = 1 - (normresidual / normX); %fraction explained by model
  fprintf(' Final fit = %e \n', fit);
end

output = struct;
output.params = params.Results;
output.iters = iter;

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function move(obj, device)
%MOVE Move the DagNN to either CPU or GPU
%   MOVE(obj, 'cpu') moves the DagNN obj to the CPU.
%
%   MOVE(obj, 'gpu') moves the DagNN obj to the GPU.

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

obj.reset() ;
obj.device = device ;
switch device
  case 'gpu'
    for i=1:numel(obj.params)
      obj.params(i).value = gpuArray(obj.params(i).value) ;
    end
  case 'cpu'
    for i=1:numel(obj.params)
      obj.params(i).value = gather(obj.params(i).value) ;
    end
  otherwise
    error('DEVICE must be either ''cpu'' or ''gpu''.') ;
end
for l = 1:numel(obj.layers)
  obj.layers(l).block.move(device) ;
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
clear all, close all, clc

% J = @(u,t)(25-(5-(u-t))^2);

J = @(u,t)(25-(5-(u)-sin(t)).^2);
u = 0;
y0 = J(u,0);

% Extremum Seeking Control Parameters
freq = 100; % sample frequency
dt = 1/freq;
T = 100; % total period of simulation (in seconds)
A = .2;  % amplitude
omega = 10*2*pi; % 10 Hz
phase = 0;
K = 5;   % integration gain

% high pass filter
butterorder=1;
butterfreq=2;  % in Hz for 'high'
[b,a] = butter(butterorder,butterfreq*dt*2,'high')
ys = zeros(1,butterorder+1)+y0;
HPF=zeros(1,butterorder+1);

uhat=u;
for i=1:T/dt
    t = (i-1)*dt;
    time(i) = t;
    yvals(i)=J(u,t);

    for k=1:butterorder
        ys(k) = ys(k+1);
        HPF(k) = HPF(k+1);
    end
    ys(butterorder+1) = yvals(i);
    HPFnew = 0;
    for k=1:butterorder+1
        HPFnew = HPFnew + b(k)*ys(butterorder+2-k);
    end
    for k=2:butterorder+1
        HPFnew = HPFnew - a(k)*HPF(butterorder+2-k);
    end
    HPF(butterorder+1) = HPFnew;

    xi = HPFnew*sin(omega*t + phase);
    uhat = uhat + xi*K*dt;
    u = uhat + A*sin(omega*t + phase);
    uhats(i) = uhat;
    uvals(i) = u;
end

%%
figure
subplot(2,1,1)
plot(time,uvals,time,uhats,'LineWidth',1.2)
l1=legend('$u$','$\hat{u}$')
set(l1,'interpreter','latex','Location','SouthEast')
grid on
subplot(2,1,2)
plot(time,yvals,'LineWidth',1.2)
ylim([-1 26])
grid on

set(gcf,'Position',[100 100 500 350])
set(gcf,'PaperPositionMode','auto')
print('-depsc2', '-loose', '../../../figures/ESC_ResponseVarying');

author: dynamicslab


text:
%% geoSphere
% Below is a demonstration of the features of the |geoSphere| function

%% Syntax
% |[F,V,Vs]=geoSphere(n,r,solidType);|

%% Description
% Use |geoSphere| to generate triangulated spheres with nearly geodesic
% triangle distributions. The density of the triangulation can be
% controlled through a particular choice of n (number of mesh refinement
% steps).

%% Examples

clear; close all; clc;

%%
% Plot Settings
fontSize=15;
faceAlpha=1;
edgeColor=0.2*ones(1,3);
edgeWidth=1.5;

%% Building a geodesic dome based on the icosahedron
% The function inputs are n and r which define the mesh refinement and
% radius respectively. The mesh refinement number n defines the number of
% subtriangulation (see function |subTri|) iterations performed on an
% icosahedron.

r=1; %sphere radius
n=2; %Refinements
[F,V,~]=geoSphere(n,r);
[Fi,Vi,~]=geoSphere(0,r);
%%
% Visualize sphere

cFigure; hold on;
subplot(1,2,1); hold on;
gpatch(Fi,Vi,'rw','r',0.8,2);
plotV(Vi,'r.','MarkerSize',50);
camlight headlight;
axisGeom(gca,fontSize);

subplot(1,2,2); hold on;
gpatch(F,V,'bw','k',1,2);
plotV(Vi,'r.','MarkerSize',50);
camlight headlight;
axisGeom(gca,fontSize);
drawnow;

%%
% Below is a visualisation for n=0:1:3. The function outputs the geodesic
% dome faces (F) and vertices (V) and also the spherical coordinates of the
% vertices (Vs) (this output is suppressed in the example below).

cFigure; % Open figure for plotting

%Defining triangulated geodesic domes with different densities
r=1; %sphere radius
n=0:1:3; %Refinements
pColors=gjet(numel(n));
for q=1:1:numel(n)
    [F,V,~]=geoSphere(n(q),r);
    subplot(2,2,q); hold on;
    title([num2str(n(q)),' refinement iterations'],'FontSize',fontSize);
    gpatch(F,V,pColors(q,:));
    % patchNormPlot(F,V);
    camlight headlight;
    axisGeom(gca,fontSize);
end
drawnow;

%% Using other solid types
% Other platonic solids can also be used as a starting tesselation. However
% these may not be as geodesic as the result for the icosahedron and
% dodecahedron.

%e.g. using a cube
solidTypes=1:5;

cFigure; % Open figure for plotting
titleCell={'tetrahedron','cube','octahedron','icosahedron','dodecahedron'};
pColors=gjet(numel(solidTypes));
for q=solidTypes
    [F,V,~]=geoSphere(0,r,q);
    subplot(2,3,q); hold on;
    title(['Based on: ',titleCell{q}],'FontSize',fontSize);
    gpatch(F,V,pColors(q,:));
    % patchNormPlot(F,V);
    camlight headlight;
    axisGeom(gca,fontSize);
end
drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*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-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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: gibbonCode


text:
function [V,F,V0,F0] = text_to_mesh(str,varargin)
  % TEXT_TO_MESH Create an 3D mesh the given text
  %
  % [V,F] = text_to_mesh(str)
  % [V,F] = text_to_mesh(str,'ParameterName',ParameterValue, ...)
  %
  % Inputs:
  %   str  string of text
  %   Optional:
  %     'Dilation'  followed by amount to dilate text image
  %     'FontName  followed by name of font. List available fonts using
  %       `/usr/local/bin/convert -list font | sed -n "s/^ *Font: //p"`
  %       {[]} (see `text_to_image`)
  %     'PointSize'  followed by height of font (in pixels) {200}
  % Outputs:
  %   V  #V by 3 list of mesh vertices
  %   F  #F by 3 list of face indices into V
  %
  % Example:
  %   im = text_to_image('Zapfino','FontName','Zapfino');imshow(im)
  %

  % Defaults
  fontname = [];
  pointsize = 200;
  filename = 'text_to_image-tmp.png';
  r = 1;
  triangle_flags = '';
  siters=10;
  tol = 0.2;
  levels = 1;
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Dilation','FontName','PointSize','SmoothingIters','Tol','TriangleFlags'}, ...
    {'r','fontname','pointsize','siters','tol','triangle_flags'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  im = text_to_image(str,'FontName',fontname,'PointSize',pointsize);
  [V0,F0] = bwmesh( ...
    imdilate(padarray(~im,[r,r]),strel('disk',r)), ...
    'Tol',tol, ...
    'SmoothingIters',siters, ...
    'TriangleFlags',triangle_flags);
  [V,F] = extrude(V0,F0,'Levels',levels);

end

author: alecjacobson

repo: gptoolbox


text:
function a = c4mat_test_inverse ( n )

%*****************************************************************************80
%
%% C4MAT_TEST_INVERSE computes the inverse of the test matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 April 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, complex A(N,N), the matrix.
%
  a = ( c4mat_test ( n ) )';

  return
end


text:
function params = getParamsKatie()

params = struct;
params.scale = 6;
params.npatchesx = 8;
params.npatchesy = 6;
params.psize = 20; % patch size for hog/hof
params.nbins = 8; % number of bins in hog/hof
params.wd = 0.5; % width of hog features for visualization

params.optflowwinsig = 3;
params.optflowsig = 2;
params.optreliability = 1e-4;

params.blocksize = 500; % block size for parallel features computation.

params.flow_thres = 1;
params.deepscale = 4;
params.deep_thres = 3.5;
params.deepmatching_step = 8;
params.warping_iters = 2;

params.presmoothing = true;
params.smoothing_sigma = 1;

% params.methods = {'deep-sup','hs-sup'};
% params.flownames = {'DS','hs_sup'};

author: kristinbranson

repo: JAABA


text:
function [N,E] = rentian_scaling_2d(A,XY,n,tol)
% RENTIAN_SCALING_2D    Rentian scaling for networks embedded in two dimensions.
%
% [N,E] = rentian_scaling_2d(A,XY,n,tol)
%
% Physical Rentian scaling (or more simply Rentian scaling) is a property
% of systems that are cost-efficiently embedded into physical space. It is
% what is called a "topo-physical" property because it combines information
% regarding the topological organization of the graph with information
% about the physical placement of connections. Rentian scaling is present
% in very large scale integrated circuits, the C. elegans neuronal network,
% and morphometric and diffusion-based graphs of human anatomical networks.
% Rentian scaling is determined by partitioning the system into cubes,
% counting the number of nodes inside of each cube (N), and the number of
% edges traversing the boundary of each cube (E). If the system displays
% Rentian scaling, these two variables N and E will scale with one another
% in loglog space. The Rent's exponent is given by the slope of log10(E)
% vs. log10(N), and can be reported alone or can be compared to the
% theoretical minimum Rent's exponent to determine how cost efficiently the
% network has been embedded into physical space. Note: if a system displays
% Rentian scaling, it does not automatically mean that the system is
% cost-efficiently embedded (although it does suggest that). Validation
% occurs when comparing to the theoretical minimum Rent's exponent for that
% system.
%
% INPUTS:
%
% A:              MxM adjacency matrix.
%                 Must be unweighted, binary, and symmetric.
% XY:             Matrix of node placement coordinates.
%                 Must be in the form of an Mx2 matrix [x y], where M is the
%                 number of nodes and x and y are column vectors of node
%                 coordinates.
% n:              Number of partitions to compute. Each partition is a data
%                 point. You want a large enough number to adequately
%                 estimate the Rent's exponent.
% tol:            This should be a small value (for example 1e-6).
%                 In order to mitigate the effects of boundary conditions due
%                 to the finite size of the network, we only allow partitions
%                 that are contained within the boundary of the network. This
%                 is achieved by first computing the volume of the convex
%                 hull of the node coordinates (V). We then ensure that the
%                 volume of the convex hull computed on the original node
%                 coordinates plus the coordinates of the randomly generated
%                 partition (Vnew) is within a given tolerance of the
%                 original (i.e. check abs(V - Vnew) < tol). Thus tol, should
%                 be a small value in order to make sure the partitions are
%                 contained largely within the boundary of the network, and
%                 thus the number of nodes and edges within the box are not
%                 skewed by finite size effects.
%
% OUTPUTS:
%
% N:              nx1 vector of the number of nodes in each of the n partitions.
% E:              nx1 vector of the number of edges crossing the boundary of
%                 each partition.
%
% Subsequent Analysis:
%
%     Rentian scaling plots are created by: figure; loglog(E,N,'*');
%
%     To determine the Rent's exponent, p, we need to determine
%     the slope of E vs. N in loglog space, which is the Rent's
%     exponent. There are many ways of doing this with more or less
%     statistical rigor. Robustfit in MATLAB is one such option:
%
%         [b,stats] = robustfit(log10(N),log10(E))
%
%     Then the Rent's exponent is b(1,2) and the standard error of the
%     estimation is given by stats.se(1,2).
%
% Note: n=5000 was used in Bassett et al. 2010 in PLoS CB.
%
% Reference:
% Danielle S. Bassett, Daniel L. Greenfield, Andreas Meyer-Lindenberg,
% Daniel R. Weinberger, Simon W. Moore, Edward T. Bullmore. Efficient
% physical embedding of topologically complex information processing
% networks in brains and computer circuits. PLoS Comput Biol, 2010,
% 6(4):e1000748.
%
% Modification History:
%
%     2010:     Original (Dani Bassett)
%     Dec 2016: Updated code so that both partition centers and partition
%               sizes are chosen at random. Also added in a constraint on
%               partition placement that prevents boxes from being located
%               outside the edges of the network. This helps prevent skewed
%               results due to boundary effects arising from the finite size
%               of the network. (Lia Papadopoulos)

% determine the number of nodes in the system
M = numel(XY(:,1));

% rescale coordinates so that they are all greater than unity
XYn = XY - repmat(min(XY)-1,M,1);

% compute the area of convex hull (i.e. are of the boundary) of the network
[~,V] = convhull(XYn(:,1),XYn(:,2));

% min and max network coordinates
xmin = min(XYn(:,1));
xmax = max(XYn(:,1));
ymin = min(XYn(:,2));
ymax = max(XYn(:,2));

% initialize vectors of number of nodes in box and number of edges crossing
% box

N = zeros(n,1);
E = zeros(n,1);

% create partitions, and count the number of nodes inside the partition (N)
% and the number of edges traversing the boundary of the partition (E)

nPartitions = 0;

while nPartitions<(n+1)

    % variable to check if partition center is within network boundary
    % OK if inside == 1
    inside = 0;

    while inside == 0

        % pick a random (x,y) coordinate to be the center of the box
        randx = xmin+(xmax-xmin)*rand(1);
        randy = ymin+(ymax-ymin)*rand(1);

        % make sure the point is inside the convex hull of the network
        newCoords = [XYn; [randx randy]];
        [~,Vnew] = convhull(newCoords(:,1),newCoords(:,2));

        % if the old convex hull area and new convex hull area are equal
        % then the box center must be inside the network boundary.

        if isequal(V,Vnew)==0
            inside = 0;
        else
            inside = 1;
        end

    end

    % determine the approximate maximum distance the box can extend, given
    % the center point and the  bounds of the network
    deltaY = min(abs(ymax-randy),abs(ymin-randy));
    deltaX = min(abs(xmax-randx),abs(xmin-randx));
    deltaLmin = min(deltaY,deltaX);

    % variable to check if partition is within network boundary
    % OK if inside == 1
    inside = 0;

    while inside == 0

        % pick a random (side length)/2 that is between 0 and the
        % max possible
        deltaL = deltaLmin*rand(1);

        % (x,y) coordinates for corners of box
        boxCoords = [randx - deltaL randy - deltaL; ...
            randx - deltaL randy + deltaL; ...
            randx + deltaL randy - deltaL; ...
            randx + deltaL randy + deltaL];

        % check if all corners of box are inside the convex hull of the
        % network
        newCoords = [XYn; boxCoords];
        [~,Vnew] = convhull(newCoords(:,1),newCoords(:,2));

        % make sure the new convex hull that includes the partition corners
        % is within a certain tolerance of the original convex hull area.

        if abs(V-Vnew)>tol
            inside = 0;
        else
            inside = 1;
        end
    end

    % Find nodes inside the box, edges crossing the boundary

    L = find(XYn(:,1)>(randx-deltaL) & XYn(:,1)<(randx+deltaL) ...
        & XYn(:,2)>(randy-deltaL) & XYn(:,2)<(randy+deltaL));

    if ~isempty(L) == 1
        nPartitions = nPartitions+1;
        % count edges crossing the boundary of the box
        E(nPartitions,1) = sum(sum(A(L,setdiff(1:M,L))));
        % count nodes inside of the box
        N(nPartitions,1) = numel(L);

    end

end

return;


text:
function [seed_sets] = generate_graph_seeds(graph_seed_gen_method, ...
                                            graph_seed_params, seg_obj)
% returns an NxS seed_sets logical matrix, where each column defines a
% seed. A value of 1 in a row indicates that that pixel/superpixel is
% included in the seed specified by the column.
%
% @authors:     Ahmad Humayun
% @contact:     ahumayun@cc.gatech.edu
% @affiliation: Georgia Institute of Technology
% @date:        Fall 2013 - Summer 2014

    if strcmpi(graph_seed_gen_method, 'sp_img_grid')
        [seed_sets] = gen_sp_img_grid(graph_seed_params, seg_obj);
    elseif strcmpi(graph_seed_gen_method, 'sp_all')
        [seed_sets] = gen_sp_all_seeds(graph_seed_params, seg_obj);
    elseif strcmpi(graph_seed_gen_method, 'sp_clr_seeds')
        [seed_sets] = gen_sp_clr_seeds(graph_seed_params, seg_obj);
    elseif strcmpi(graph_seed_gen_method, 'sp_seed_sampling')
        [seed_sets] = AbstractGraph.gen_sp_seed_sampling(graph_seed_params, seg_obj);
    elseif strcmpi(graph_seed_gen_method, 'gen_user_seeds')
        [seed_sets] = gen_user_seeds(graph_seed_params, seg_obj);
    else
        error('AbstractGraph:generate_graph_seeds', ...
              '''%s'' invalid seed generation methos', ...
              graph_seed_gen_method);
    end

    % this function ensures no seeds overlap
    [seed_sets] = no_overlapping_seeds(seed_sets, seg_obj);
end

function [seed_sets] = gen_sp_img_grid(graph_seed_params, seg_obj)
% generates seeds by a grid over the image. The first array in
% graph_seed_params specifies how many seeds are needed. The second array
% specifies the extent of the seed pixels around each seed grid location

    graph_num_seed_locs = graph_seed_params{1};
    graph_seed_sz = graph_seed_params{2};
    % generate the offsets from the central grid pixel
    [rect_coords] = ...
        generate_seeds_funcs('generate_rectangle_coords', ...
                             graph_seed_sz);
    % generate the seed sets over a grid (each column gives
    % pixels/superpixels selected for that seed)
    seed_sets = generate_seeds_funcs('generate_sp_img_grid', ...
                                     seg_obj.sp_data.sp_seg, ...
                                     rect_coords, ...
                                     graph_num_seed_locs, ...
                                     seg_obj.sp_data.sp_seg_szs, 0);
end

function [seed_sets] = gen_sp_all_seeds(graph_seed_params, seg_obj)
% generates seeds by choosing all sps. The single parameter in
% graph_seed_params indicate how many sps to skip when choosing a sp as a
% seed

    subsample_sps = graph_seed_params{1} + 1;

    % initially each sp is a seed
    seed_sets = logical(eye(seg_obj.sp_data.num_spx, 'int8'));

    % now subsample the seeds
    seed_sets = seed_sets(:,1:subsample_sps:end);
end

function [seed_sets] = gen_sp_clr_seeds(graph_seed_params, seg_obj)
% generates seeds by finding superpixels locations which are closest to a
% grid over the image. The first array in graph_seed_params specifies how
% many seeds are needed. The second array specifies the extent of the seed
% pixels around each seed grid location. The third argument is a string
% specifying what superpixel method to use. 'sp_seeds_caller' uses the
% superpixels already generated in seg_obj.sp_data;
% 'felzenszwalb_seeds_caller' uses FelzHutten superpixels.

    graph_num_seed_locs = graph_seed_params{1};
    graph_seed_sz = graph_seed_params{2};
    superpixel_method = graph_seed_params{3};
    % generate the offsets from the central grid pixel
    [rect_coords] = ...
        generate_seeds_funcs('generate_rectangle_coords', ...
                             graph_seed_sz);
    % generate the seed sets over a grid (each column gives
    % pixels/superpixels selected for that seed)
    seed_sets = generate_seeds_funcs('generate_sp_seeds_superpixels', ...
                                     seg_obj.sp_data, seg_obj.I, ...
                                     superpixel_method, 0, ...
                                     rect_coords, graph_num_seed_locs);
end

function [seed_sets] = gen_user_seeds(graph_seed_params, seg_obj)
    graph_seed_sz = graph_seed_params{1};
    user_centers = graph_seed_params{2};
    [rect_coords] = ...
        generate_seeds_funcs('generate_rectangle_coords', ...
                             graph_seed_sz);
    seed_sets = generate_seeds_funcs('generate_sp_user', ...
                                     seg_obj.sp_data.sp_seg, ...
                                     rect_coords, user_centers, ...
                                     seg_obj.sp_data.sp_seg_szs, 0);
end

function [seed_sets] = no_overlapping_seeds(seed_sets, seg_obj)
% this function ensures that no seeds overlap i.e. no superpixels are the
% same in two seeds

    % find number of overlapping superpixels across seeds
    ss = single(seed_sets);
    num_overlap_sps = ss' * ss;
    num_overlap_sps = triu(num_overlap_sps, 1);
    [s1,s2,os] = find(num_overlap_sps);
    num_overlap_sps = [s1, s2, os];

    num_sps_per_seed = sum(seed_sets);

    % remove overlapping superpixel for each seed pair contention (note,
    % superpixels removed from a seed in a pervious iteration, would no
    % longer be taken into consideration in later iterations - hence,
    % correctly ignoring superpixels which no longer belong to a seed)
    for idx = 1:size(num_overlap_sps,1)
        s1 = num_overlap_sps(idx,1);
        s2 = num_overlap_sps(idx,2);
        n1 = num_sps_per_seed(s1);
        n2 = num_sps_per_seed(s2);

        % find the sps which overlap
        overlap_sps = all(seed_sets(:,[s1 s2]), 2);

        % if seed1 has more superpixels, then remove overlapping sp from
        % seed1. Otherwise remove superpixel from seed2
        if n1 > n2
            seed_sets(overlap_sps,s1) = 0;
        else
            seed_sets(overlap_sps,s2) = 0;
        end
    end

    % discard any seeds which have no superpixels anymore
    num_sps_per_seed = sum(seed_sets);
    seed_sets(:,num_sps_per_seed == 0) = [];
end

author: Cloud-CV

repo: object-proposals


text:
function [img, hdr_info] = hdrimread(filename)
%
%       [img, hdr_info] = hdrimread(filename)
%
%       This function reads from a file with name filename an HDR image, if
%       the format can not be opened, it tries to open it as it was an LDR
%       image using imread from MATLAB Image Processing Toolbox.
%
%       NOTE: JPEG2000 file passed as input are meant to be compressed
%       using HDR JPEG-2000.
%
%        Input:
%           -filename: the name of the file to open
%
%        Output:
%           -img: the opened image
%           -hdr_info: RGBE format extra datum such as: exposure, gamma, etc.
%
%     Copyright (C) 2011-15  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('filename', 'var'))
    error('A filename with extension needs to be passed as input!');
end

extension = lower(fileExtension(filename));

%Radiance format can have different extensions: .hdr, .rgbe and .pic
if((strcmpi(extension,'pic')==1) || (strcmpi(extension,'rgbe')==1))
    extension = 'hdr';
end

img = [];
hdr_info = [];
bLDR = 1;

switch extension

    %PIC-HDR format by Greg Ward
    case 'hdr'
        try
	    %HDR Toolbox's HDR Reader
            [img, hdr_info] = read_rgbe(filename);
            bLDR = 0;
        catch err
            try
                %MATLAB HDR Reader
                img = double(hdrread(filename));
                bLDR = 0;
            catch err
                warning('This .hdr/.pic file cannot be read.');
            end
        end

    %OpenEXR support using TinyEXR
    case 'exr'
        try
            img = read_exr(filename);
            bLDR = 0;
        catch err
            hdr_info = struct('loaded', 0);
            warning('This .exr file cannot be read.');
        end

    %Portable float map
    case 'pfm'
        try
            img = read_pfm(filename);
            bLDR = 0;
        catch
            hdr_info = struct('loaded', 0);
            warning('This .pfm file cannot be read.');
        end

    %HDR JPEG-2000
    case 'jp2'
        try
            img = HDRJPEG2000Dec(filename);
            bLDR = 0;
        catch err
            hdr_info = struct('loaded', 0);
            warning('This .jpg2 file cannot be read as an HDR JPEG 2000 file.');
        end
end

if(isempty(img))
    if(bLDR == 1)
        img = ldrimread(filename);
        if(isempty(img))
            hdr_info = struct('loaded', 0);
            warning(['This image, ', filename,', cannot be loaded as an LDR image.']);
        else
            hdr_info = struct('loaded', 1);
            warning(['This image, ', filename,', has been loaded as an LDR image.']);
        end
    else
        hdr_info = struct('loaded', 0);
        error(['This image, ',filename,', cannot be loaded with LDR or HDR readers.']);
    end
else
    hdr_info.loaded = 1;
end

%Remove specials
img = RemoveSpecials(img);

end

author: banterle


text:
function bar_plot( res, column, tasks, dets, varargin)

opts.legend = true;
opts.legendArgs = {'Location', 'SouthOutside', 'Interpreter', 'none', 'Orientation', 'Horizontal'};
opts.sort = true;
opts.detnames = {dets.name};
opts.multiplier = 100;
opts = vl_argparse(opts, varargin);

if isfield(dets, 'normname')
  normnames = {dets.normname};
else
  normnames = [];
end
data = rproc.pick(res, tasks, {dets.name}, column, 'normname', normnames);
data = data .* opts.multiplier;
avg = mean(data, 1);
if opts.sort
  [avg, order] = sort(avg, 'ascend');
  dets = dets(order);
  data = data(:, order);
  opts.detnames = opts.detnames(order);
end

for ni = numel(dets):-1:1
  position = ni;
  barh(position, avg(ni), 'FaceColor', dets(ni).color, dets(ni).bararg{:});
  hold on;

  text(101, position, sprintf('%0.2f%%', avg(ni)),...
    'HorizontalAlignment', 'left',...
    'VerticalAlignment', 'middle');

  li = zeros(1, numel(tasks));
  for ti = 1:numel(tasks)
    li(ti) = plot(data(ti, ni), position, '.', tasks(ti).style{:});
  end
end
set(gca,'Ydir','reverse')
set(gca, 'YTick', 1:numel(dets));
set(gca, 'TickLabelInterpreter', 'none');
set(gca, 'YTickLabel', opts.detnames);

if opts.legend
  legend(li, tasks, opts.legendArgs{:});
end
axis tight;
set(gca, 'XLim', [0, 100]);
grid on;
end

author: hpatches

repo: hpatches-benchmark