text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% To run:  Edit user options in mtfcalc.m and execute the code.
%
% This code is intended to accompany the paper:
%
% S. N. Friedman, G. S. K. Fung, J. H. Siewerdsen, and B. M. W. Tsui.  "A
% simple approach to measure computed tomography (CT) modulation transfer
% function (MTF) and noise-power spectrum (NPS) using the American College
% of Radiology (ACR) accreditation phantom," Med. Phys. 40, 051907-1 -
% 051907-9 (2013).
% http://dx.doi.org/10.1118/1.4800795
%
% This code is free to distribute (see below).
%
% This program requires the following 2 files:
%   mtfcalc.m
%   license.txt
%
% Purpose:  To calculate the 1D radial (axial) MTF using CT data of the
%           American College of Radiology (ACR) accreditation phantom.
%
% Input:    Two consecutive scans of the phantom are needed.  The program
%           requires a data directory to be selected in which there are
%           only two subdirectories containing only CT slices corresponding
%           to the third module of the phantom.  Be careful of partial
%           volume effects with surrounding modules.
%
%           i.e.,  datadir
%                      |-> scan 1 dir
%                      |           | -> only module 3 slices
%                      |
%                      |-> scan 2 dir
%                                  |-> only module 3 slices
%
% Copyright 2012 Saul N. Friedman
%   Distributed under the terms of the "New BSD License."  Please see
%   license.txt.
%
% N.B.:  This code assumed the data are already organized such that only
%        relevant axial slices are contained within the data directory.
%        ESF = edge spread function
%        LSF = line spread function
%        MTF = modulation transfer function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all;
clear all;
clc;

fprintf('****************************************************************************\n');
fprintf('CT MTF Calculation\n\n')
fprintf('This code is intended to accompany the paper:\n');
fprintf('S. N. Friedman, G. S. K. Fung, J. H. Siewerdsen, and B. M. W. Tsui.  "A\n')
fprintf('simple approach to measure computed tomography (CT) modulation transfer \n')
fprintf('function (MTF) and noise-power spectrum (NPS) using the American College \n')
fprintf('of Radiology (ACR) accreditation phantom," Med. Phys. 40, 051907-1 - \n')
fprintf('051907-9 (2013).\n')
fprintf('http://dx.doi.org/10.1118/1.4800795\n\n')
fprintf('This code can be freely distributed.  Please see license.txt for details.\n');
fprintf('****************************************************************************\n\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% User selected options
%

selectdir = 1;   % 1 = prompt for directly locations via dialogue box, 0 = hardcode locations below
    datadir = fullfile('data');
    resultsdir = fullfile('results',datadir);

saveresults = 1;  % Choose whether to save the calculated MTF to disk in text files

showfig = 3;  % decide how important a figure should be in order to be plotted, 1 = no figures, 4 = all figures

rstep = 0.1; %  Choose bin size in mm (i.e., effective pixel size) for oversampled edge spread function (ESF)

window = 1;  % Decide whether to apply a Hann window to the LSF before calculating the MTF (typically applied)

% Background and foreground values:  1 = prompted to draw rectangle to select ROI, 0 = use values below
manualROI = 1;
   % box border for background ROI in pixels
   BGxROI1 = 1;
   BGxROI2 = 100;
   BGyROI1 = 1;
   BGyROI2 = 100;

   % box border for foreground ROI in pixels
   FGxROI1 = 205;
   FGxROI2 = 305;
   FGyROI1 = 208;
   FGyROI2 = 308;

% 2D map of good pixels (1) and bad pixels (0) to use in calculation
% 1 = prompted to draw rectangle(s) to select ROI(s), 0 = use values below
manualmask = 1;
selectmask = 1; %prompt for image to use when creating data mask (only applicable if manualmask = 1)
    mask = ones(512,512);
    mask(450:512,1:125) = 0;
    mask(450:512,400:512) = 0;
    mask(400:425,84:105) = 0;
    mask(410:430,400:420) = 0;
    mask(220:230,279:288) = 0;
    mask(371:381,138:147) = 0;

% Determine which angles of data to use in radians (-pi to pi to use all)
thlim1 = -pi;
thlim2 = pi;

%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start of analysis code
%

homedr = cd;

if selectdir ==1
    fprintf('Select the data directory.\n\n');
    datadir = uigetdir(homedr,'Pick data directory');
    if saveresults == 1
        fprintf('Select the results directory.\n\n');
        resultsdir = uigetdir(homedr,'Pick results directory');
    end
else
    datadir = fullfile(homedr,datadir);
    resultsdir = fullfile(homedr,resultsdir);
end

if saveresults == 1
    mkdir(resultsdir)
end

cd(datadir)

scandir = dir;

cd(homedr)

while or(strcmp(scandir(1).name,'.'),strcmp(scandir(1).name,'..'))
    scandir = scandir(2:end);
end

Nscan = size(scandir,1);
Nscan = 1;

nfilescheck = 0;

for i = 1:Nscan
    cd(fullfile(datadir,scandir(i).name))

    % get list of filenames

    ffiles = dir;

    cd(homedr)

    while or(strcmp(ffiles(1).name,'.'),strcmp(ffiles(1).name,'..'))
        ffiles = ffiles(2:end);
    end

    Nz = size(ffiles,1);

     if nfilescheck ~= Nz && nfilescheck ~=0
        fprintf('Error.  Scan directories do not contain the same number of slices.\n')
        return
    end

    nfilescheck = Nz;

    if i == 1

        % Read DICOM header info and determine correct mapping function to get HU
        % values as well as voxel sizes.

        im = single(dicomread(fullfile(datadir,scandir(i).name,ffiles(1).name)));

        imheader = dicominfo(fullfile(datadir,scandir(i).name,ffiles(1).name));

        b = imheader.RescaleIntercept;
        m = imheader.RescaleSlope;

        im = im.* m + b;

        pixelx = imheader.PixelSpacing(1);
        pixely = imheader.PixelSpacing(2);
        pixelz = imheader.SliceThickness;

        if showfig > 2
            figure
            imagesc(im)
            colormap gray
            hold on
            axis image;
            impixelinfo;
        end

        Nx = size(im,2);
        Ny = size(im,1);

        imaxisx = [1:Nx] .* pixelx;
        imaxisy = [1:Ny] .* pixely;

        imavg = zeros(Ny,Nx);
        im3D = zeros(Ny,Nx,Nz);
    end

    % Load data into memory and calculate an average 2D image to later
    % determine phantom location

    for z = 1:Nz
        im = single(dicomread(fullfile(datadir,scandir(i).name,ffiles(z).name)));
        im = im.* m + b;
        im3D(:,:,z,i) = im;
        imavg = imavg + im;
    end

end

imavg = imavg ./(Nz*Nscan);

if selectmask == 1 && manualmask == 1
    fprintf('Select the DICOM image to use for mask creation.\n\n');
    [imfile impath] = uigetfile('*','Pick a DICOM image to use for mask creation',datadir);
    if isequal(imfile,0) || isequal(impath,0)
        im = imavg;
    else
        im = single(dicomread(fullfile(impath,imfile)));
        im = im.* m + b;
    end
else
    im = imavg;
end

% Plot sample slice.  Use to create mask image if set to manual selection.

imaxisx = [1:size(im,2)] .* pixelx;
imaxisy = [1:size(im,1)] .* pixely;

if showfig > 3 || manualmask == 1
    figure
    imagesc(im)
    xlabel('pixel')
    ylabel('pixel')
    colormap gray
    hold on
    axis image;
    impixelinfo;
end

if manualmask ==1

    mask = ones(size(im));

    selectROI = 1;

    fprintf('Select the ROIs corresponding to bad data.\n\n');

    while selectROI == 1
        choice = menu('Select additional ROI of bad data?','Yes','No');

        if choice == 2
            selectROI = 0;
            hold off
        else
            h = imrect;
            ptmask = wait(h);

            maskx1 = round(ptmask(1));
            maskx2 = maskx1 + round(ptmask(3));
            masky1 = round(ptmask(2));
            masky2 = masky1 + round(ptmask(4));

            % force selected ROI to be within the bounds of the image
            maskx1 = min(maskx1,size(im,2));
            maskx2 = min(maskx2,size(im,2));
            masky1 = min(masky1,size(im,1));
            masky2 = min(masky2,size(im,1));
            maskx1 = max(maskx1,1);
            maskx2 = max(maskx2,1);
            masky1 = max(masky1,1);
            masky2 = max(masky2,1);

            mask(masky1:masky2,maskx1:maskx2) = 0;

            plot([maskx1,maskx2,maskx2,maskx1,maskx1],[masky1,masky1,masky2,masky2,masky1],'r-')
        end
    end
end

% Plot mask image overlying sample image to verify mask

if showfig > 3
    figure
    imagesc([imaxisy(1) imaxisy(end)],[imaxisx(1) imaxisx(end)],(im+100).*mask)
    xlabel('mm')
    ylabel('mm')
    colormap gray
    hold on
    axis image;
    impixelinfo;
end

% Show averaged image if desired or needed

if showfig > 2 || manualROI == 1
    figure
    imagesc(imavg)
    xlabel('pixel')
    ylabel('pixel')
    colormap gray
    hold on
    axis image;
    impixelinfo;
end

% Determine background and foreground values to properly normalized the
% data to range [0,1]

if manualROI == 1
    fprintf('Click and drag to create a rectangular ROI representing the background.  Double click the ROI when finished.\n\n');

    h = imrect;
    ptROI = wait(h);

    BGxROI1 = round(ptROI(1));
    BGxROI2 = BGxROI1 + round(ptROI(3));
    BGyROI1 = round(ptROI(2));
    BGyROI2 = BGyROI1 + round(ptROI(4));

    % force selected ROI to be within the bounds of the image
    BGxROI1 = min(BGxROI1,size(im,2));
    BGxROI2 = min(BGxROI2,size(im,2));
    BGyROI1 = min(BGyROI1,size(im,1));
    BGyROI2 = min(BGyROI2,size(im,1));
    BGxROI1 = max(BGxROI1,1);
    BGxROI2 = max(BGxROI2,1);
    BGyROI1 = max(BGyROI1,1);
    BGyROI2 = max(BGyROI2,1);
end

BGvalue = mean(mean(imavg(BGyROI1:BGyROI2,BGxROI1:BGxROI2)));

if manualROI == 1
    fprintf('Click and drag to create a rectangular ROI representing the foreground.  Double click the ROI when finished.\n\n');

    h = imrect;
    ptROI = wait(h);

    FGxROI1 = round(ptROI(1));
    FGxROI2 = FGxROI1 + round(ptROI(3));
    FGyROI1 = round(ptROI(2));
    FGyROI2 = FGyROI1 + round(ptROI(4));

    % force selected ROI to be within the bounds of the image
    FGxROI1 = min(FGxROI1,size(im,2));
    FGxROI2 = min(FGxROI2,size(im,2));
    FGyROI1 = min(FGyROI1,size(im,1));
    FGyROI2 = min(FGyROI2,size(im,1));
    FGxROI1 = max(FGxROI1,1);
    FGxROI2 = max(FGxROI2,1);
    FGyROI1 = max(FGyROI1,1);
    FGyROI2 = max(FGyROI2,1);
end

FGvalue = mean(mean(imavg(FGyROI1:FGyROI2,FGxROI1:FGxROI2)));

if or(showfig > 2,manualROI ==1)
    plot([BGxROI1,BGxROI2,BGxROI2,BGxROI1,BGxROI1],[BGyROI1,BGyROI1,BGyROI2,BGyROI2,BGyROI1],'r-')
    plot([FGxROI1,FGxROI2,FGxROI2,FGxROI1,FGxROI1],[FGyROI1,FGyROI1,FGyROI2,FGyROI2,FGyROI1],'r-')
    hold off
end

% Normalize averaged image such that values are [0,1]

imavg = (imavg - BGvalue) ./ (FGvalue - BGvalue);

if showfig > 2
    figure
    imagesc([imaxisy(1) imaxisy(end)],[imaxisx(1) imaxisx(end)],imavg)
    xlabel('mm')
    ylabel('mm')
    title('averaged and normalized image')
    colormap gray
    hold on
    axis image;
    impixelinfo;
end

% Determine which pixels correspond to the phantom cylinder

imlog = zeros(size(imavg));

level = graythresh(imavg);
imlog = im2bw(imavg,level);

% Calculate the center of the cylinder

fprintf('Calculating the center of the phantom and determining phantom-pixel locations.\n\n');

if exist('bwconncomp')==2
    CC = bwconncomp(imlog,26);
    for i = 1:CC.NumObjects;
        Lobj(i) =  length(CC.PixelIdxList{1,i}) ;
    end
else
    [CC CCn] = bwlabel(imlog,8);
    for i = 1:CCn
        Lobj(i) = length(find(CC == i));
    end
end

L = find(Lobj == max(Lobj));

S = regionprops(CC,'Centroid');

c.y = S(L).Centroid(2);
c.x = S(L).Centroid(1);

if showfig>2
    figure
    imagesc([imaxisy(1) imaxisy(end)],[imaxisx(1) imaxisx(end)],imlog)
    xlabel('mm')
    ylabel('mm')
    title('Cylinder with center marked')
    colormap gray
    hold on
    axis image;
    impixelinfo;
    if exist('bwconncomp')==2
        [I J] = ind2sub(size(imavg),CC.PixelIdxList{1,L});
    else
        [I,J] = find(CC == L);
    end
    plot(J*pixely,I.*pixelx,'.r')
    plot(c.x*pixelx,c.y*pixely,'*y')
    hold off
end

xaxis = ([1:size(imavg,2)] - c.x) .* pixelx;
yaxis = ([1:size(imavg,1)] - c.y) .* pixely;

% Convert all image coords to polar coords with origin at center of the
% cylinder

r = zeros(size(imavg));
th = zeros(size(imavg));

for i = 1 :size(imavg,2)
    for j = 1 : size(imavg,1)
        %         r(j,i) = sqrt(xaxis(i).^2 + yaxis(j).^2);
        %         th(j,i) = atan2(yaxis(j),xaxis(i));
        [th(j,i) r(j,i)] = cart2pol(xaxis(i),yaxis(j));
    end
end

% Create oversampled ESF excluding data specified in mask

rmax = max(max(r));

esf = zeros(ceil(rmax/rstep)+1,1);
nsamp = zeros(length(esf),1);

rbin = 0:rstep:rmax+rstep;

if showfig > 3
    figure
    imagesc((imavg+1).*mask)
    colormap gray
    hold all
    axis image;
    impixelinfo;
end

% data with radius falling within a given bin are averaged together for a
% low noise approximation of the ESF at the given radius

fprintf('Calculating the MTF.\n');
fprintf('The may take several minutes. \n');
fprintf('Please be patient.\n\n');

thlim1 = min(thlim1,thlim2);
thlim2 = max(thlim1,thlim2);

for i = 1:length(rbin)
    R1 = find(r >= rbin(i));
    R2 = find(r < rbin(i) + rstep);
    R3 = intersect(R1,R2);
    R4 = find(th >= thlim1);
    R5 = find(th < thlim2);
    R6 = intersect(R5,R4);
    R7 = intersect(R3,R6);
    R8 = R7.*mask(R7);
    R=R8(R8~=0);
    if showfig > 3
        [X Y] = ind2sub(size(imavg),R);
        plot(Y,X,'.')
        hold all
    end
    esf(i) = sum(imavg(R));
    nsamp(i) = length(R);
end

i1 = find(nsamp,1,'first');
i2 = find(nsamp,1,'last');
nsamp = nsamp(i1:i2);
esf = esf(i1:i2);
rbin = rbin(i1:i2);

I = find(nsamp > 0);
esf(I) = esf(I)./nsamp(I);

% Plot oversampled ESF
if showfig > 2
    figure
    plot(rbin,esf)
    xlabel('radius (mm)')
    ylabel('normalized + oversampled ESF')
    grid on
    axis([rbin(1) rbin(end) -1.2*min(esf) 1.1*max(esf)])
end

% Calculated the LSF from the ESF
lsf = diff(esf);

% artifacts caused by empty bins likely have a value of -/+ 1, so try to
% remove them.
I = find(abs(lsf)> 0.9);
lsf(I) = 0;

if max(lsf) < max(abs(lsf))
    lsf = -lsf;
end

maxlsf = max(lsf);

rcent = find(lsf == maxlsf);
lsfaxis = rbin(1:end-1) - rbin(rcent);

% Center edge location in LSF and keep surrounding data range specified by
% pad.  N.B. data from regions greater than the FOV is removed in
% this manner.
pad = round(length(lsfaxis)/5);
lsfaxis = lsfaxis(rcent-pad:rcent+pad);
lsf = lsf(rcent-pad:rcent+pad);

nlsf = length(lsf);

%calculate Hann window and apply to data if specified.
if window == 1
    w = 0.5 .* (1 - cos(2*pi*[1:length(lsfaxis)]./length(lsfaxis)));
    w = w';
else
    w = ones(length(lsfaxis),1);
end
lsfw = lsf .* w;

% Plot windowed LSF and scaled window together
if showfig > 1
    figure
    plot(lsfaxis,lsfw)
    hold on
    plot(lsfaxis,w.*max(lsfw),'r--')
    hold off
    legend('oversampled + normalized LSF','scaled Hann window')
    xlabel('radius (mm)')
    ylabel('signal')
    axis([lsfaxis(1) lsfaxis(end) 1.1*min(lsfw) 1.1*max(lsfw)])
    grid on
end

% Calculate the MTF from the LSF

T = fftshift(fft(lsfw));
faxis = ([1:nlsf]./nlsf - 0.5) ./ rstep;

MTF = abs(T);

%Plot the MTF
if showfig > 1
    figure
    plot(faxis,MTF,'.-')
    xlabel('spatial frequency (mm^{-1})');
    ylabel('MTF_r')
    axis([0 max(faxis) 0 1.1*max(MTF)])
    grid on
end

fprintf('MTF calculated.\n\n');

% Save MTF results if indicated
if saveresults == 1
    save(fullfile(resultsdir,'axisvalues.txt'),'faxis','-ASCII')
    save(fullfile(resultsdir,'MTFvalues.txt'),'MTF','-ASCII')
    fprintf('Results saved.\n\n');
end

return;
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


text:
function preserve_history ()
% Maintains command history indefinitely in a separate history_saved.m file
% since Matlab is truncating history.m to a fixed size. The code is fast even
% when history_saved.m gets very big.
%
% Also works with the new format introduced in R2014a which is XML. Its
% contents are parsed as plain commands and merged into history_saved.m
%
% Bogdan Roman, 2014, abr28@cam.ac.uk, University of Cambridge

savpath = fullfile(prefdir,'history_saved.m');

% If History.xml exists (R2014a or later) read that, else revert to history.m
hpath = fullfile(prefdir,'History.xml');
if exist(hpath,'file')
    % Read the XML file into a string to look as the older history.m (i.e.
    % remove all XML tags)
    mathist = fileread(hpath);
    mathist = regexprep(mathist, '(<[^>]+>\s*)+', '\n', 'lineanchors');
    % translate html entities and remove leading newline
    mathist = strrep(mathist(3:end), '>', '>');
    mathist = strrep(mathist, '<', '<');
else
    mathist = fileread(fullfile(prefdir,'history.m'));
end

% replace \r and \r\n with \n (safe to copy between OSes)
mathist = regexprep(mathist, '\r(\n)?', '\n');

fprintf('Merging current command history into %s ...\n', savpath);

% if empty saved file, dump it all
if ~exist(savpath,'file')
    fid = fopen(savpath, 'w');
    fwrite(fid, mathist);
    fclose(fid);
    return;
end

% Read saved file fully and use regexp(). Faster than parsing line by line.
savhist = fileread(savpath);
%savhist = regexprep(savhist, '\r(\n)?', '\n');

% Look for last saved date in the preserved file. Must be a valid date and
% not 'Unknown date', which may match in both files and overwrite the entire
% saved history. Note the \d to that effect.
[savstart, savend] = regexp(savhist,'^%-- \d.*?$','lineanchors');
if isempty(savstart)
    % If none found, then append current history entirely
    fid = fopen(savpath,'a');
    fwrite(fid, mathist);
    fclose(fid);
else
    % Extract last save date as string and search for it in current history
    savlastdate = savhist(savstart(end):savend(end)-1);
    matstart = regexp(mathist, ['^' savlastdate], 'lineanchors');
    % If not found, append the entire current history
    if isempty(matstart)
        fid = fopen(savpath,'a');
        fwrite(fid, mathist);
        fclose(fid);
    % Otherwise merge the chunks
    else
        fid = fopen(savpath, 'r+');
        fseek(fid, savstart(end)-1, 'bof');
        fwrite(fid, mathist(matstart(1):end));
        fclose(fid);
    end
end
disp('Done');

author: alecjacobson

repo: gptoolbox


text:
function loops=extractloops(edges)
%
% loops=extractloops(edges)
%
% extract individual loop or polyline segment from a collection of edges
%
% author: Qianqian Fang, <q.fang at neu.edu>
% date: 2007/11/21
%
% input:
%    edges:  two column matrix recording the starting/ending
%             points of all edge segments
%
% output:
%    loops:  output, a single vector separated by NaN, each segment
%             is a 3D polyline or loop consisted of node IDs
%
% example:
%    edges=[1 2;2 3;1 4;3 4;7 3;1 9;5 6;6 7;10 9; 8 10;1 8;9 3;11 11;11 12];
%    loops=extractloops(edges)
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

loops=[];
edges(edges(:,1)==edges(:,2),:)=[]; % remove degenerated edges
loops=[loops,edges(1,:)];
loophead=edges(1,1);
loopend=edges(1,end);
edges(1,:)=[];

while(~isempty(edges))
    idx=[find(edges(:,1)==loopend)',find(edges(:,2)==loopend)'];
    if(length(idx)>1) % when a node with multiple connection found
        idx=idx(1);   % take the first connection and continue
    end
    if(isempty(idx)) % when an open-line segment gets to one end
        % when both open ends are found
        if(isempty([find(edges(:,1)==loophead)',find(edges(:,2)==loophead)']))
            loops=[loops,nan];
            loops=[loops,edges(1,:)];
            loophead=edges(1,1);
            loopend=edges(1,end);
            edges(1,:)=[];
        else % only the first open end is found, flip and trace the other
            [loophead, loopend]=deal(loopend, loophead);
            lp=fliplr(loops);
            seg=find(isnan(lp),1);
            if(isempty(seg))
                loops=lp;
            else
                loops=[loops(1:end-seg(1)+1) lp(1:seg(1)-1)];
            end
        end
        continue;
    end
    if(length(idx)==1) % tracing along a single line thread
        idx=idx(1);
        ed=edges(idx,:);
        ed(ed==loopend)=[];
        newend=ed(1);
        if(newend==loophead)  % when a loop is found
            loops=[loops loophead nan];
            edges(idx,:)=[];
            if(size(edges,1)==0) break; end
            loops=[loops,edges(1,:)];
            loophead=edges(1,1);
            loopend=edges(1,end);
            edges(1,:)=[];
            continue;
        else
            loops=[loops,newend];
        end
        loopend=newend;
        edges(idx,:)=[];
    end
end

author: fangq


text:
function [s,S1,S2] = makeSegment(p1,p2)

% MAKESEGMENT  Make a segment out of two endpoints
%   MAKESEGMENT(E1,E2) makes a segment out of the two endpoints E1 and E2
%   by stacking both endpoints [E1;E2];
%
%   [s,S1,S2] = MAKESEGMENT(E1,E2) returns the Jacobians of the segments
%   wrt the endpoints.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

s = [p1(:);p2(:)];

if nargout>1

    n  = numel(p1);
    S1 = [eye(n);zeros(n)];
    S2 = [zeros(n);eye(n)];

end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [bg, cga, cgb, tg] = computeDepthCues(I, colorParam)
	[mPb_nmax, mPb_nmax_rsz, bg1, bg2, bg3, cga1, cga2, cga3, cgb1, cgb2, cgb3, tg1, tg2, tg3, textons] = multiscalePb(im2double(I));

	bg = cat(4, bg1, bg2, bg3);
	tg = cat(4, tg1, tg2, tg3);
	cgb = cat(4, cgb1, cgb2, cgb3);
	cga = cat(4, cga1, cga2, cga3);
end

author: s-gupta

repo: rgbd

save_path: github-repos/MATLAB/s-gupta-rgbd


text:
%   AUTHORSHIP
%   Math Lead & Primary Developer:  Connor Meehan <connor.gw.meehan@gmail.com>
%   Secondary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Funded by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause
%
function y=probability_bin_labels(ptrs, labels, N)
assert(max(ptrs)<=N, 'Expect bin pointer can be higher than bin');
if size(ptrs,2)==1
    assert(size(ptrs,1)>N, 'Expect more bin pointers than bins');
elseif size(ptrs,1)==1
    assert(size(ptrs,2)>N, 'Expect more bin pointers than bins');
else
    assert(false, 'Expect 1 column for bin pointers');
end
y=zeros(1,N);
y(ptrs)=labels;
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function figbolden(pos,pv,m)
% embolden the current figure
% Inputs: pos = [xmin ymin width height] gives the lower left corner position and the window size in pixels
%               [width height] leaves the lower left corner alone
%               [width] has a standard aspect ratio of 4:3
%               [-width/height] leaves the area unchanged but fixes the aspect ratio
%         pv is a cell name containing attribute-value pairs.
%            default = {'FontName' 'Ariel'; 'FontSize' 16; 'LineWidth' 2; 'MarkerSize' 8}
%         m is a mode string:
%                'l' lists the changes made
%                'd' uses default pv settings
%                'x' suppresses changes

%      Copyright (C) Mike Brookes 2003
%      Version: $Id: figbolden.m,v 1.10 2011/02/28 12:52:10 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ps={'Title' 'XLabel' 'YLabel' 'Children'};
if nargin<3
    m='';
end
if nargin<2 || any(m=='d')
    pv={'FontName' 'Arial'; 'FontSize' 16; 'LineWidth' 2; 'MarkerSize' 8};
end
pp={'Symbol';'Wingdings'};      % protected fonts
if nargin<1
    pos=[];
end
mlist=any(m=='l');  % list changes
mnotx=~any(m=='x'); % do changes
scsz=get(0,'screensize');
if length(pos)
    po=get(gcf,'position');
    if length(pos)>2            % position is specified
        po(1:2)=pos(1:2);
        pos(1:2)=[];      % remove xmin,ymin
    end
    if length(pos)>1
        po(3:4)=pos(1:2);
    else
        if pos(1)>0
            po(3:4)=[1 0.75]*pos(1);
        else
            po(3:4)=[-pos(1) 1]*sqrt(-po(3)*po(4)/pos(1)); % preserve area
        end
    end
    set(gcf,'position',po);
end
hlist=get(gcf,'children');
while length(hlist)
    pl=get(hlist(1));
    %fprintf('list length = %d, handle = %f\n',length(hlist),hlist(1));
    for i=1:size(pv,1)
        if isfield(pl,pv{i,1})
            if i>1 || all(~strcmp(get(hlist(1),pv{i,1}),pp))
                pval=get(hlist(1),pv{i,1});
                if ~all(size(pval)==size(pv{i,2})) || ~all(pval(:) == pv{i,2}(:))
                    if mnotx
                        set(hlist(1),pv{i,1},pv{i,2})
                    end
                    if mlist
                        if ischar(pval)
                            fch='s';
                        else
                            fch='g';
                        end
                        fprintf(['change %f %s: %' fch ' -> %' fch '\n'],hlist(1),pv{i,1},pval,pv{i,2});
                    end
                end
            end
        end
    end
    for i=1:length(ps)
        if isfield(pl,ps{i})
            hlist=[hlist; get(hlist(1),ps{i})];
            %fprintf('add %f:%s\n',hlist(1),ps{i});
        end
    end
    hlist(1)=[];
end

author: decouples


function [ ] = loop_compare_diff_domain( cutCostsFile, max_cut_num, search_domain_1, search_domain_2, search_domain_3, max_local_dis )
%LOOP_COMPARE_DIFF_SEARCH Summary of this function goes here
%   Detailed explanation goes here

    DP_cost_list_1 = []; DP_time_list_1 = [];
    DP_cost_list_2 = []; DP_time_list_2 = [];
    DP_cost_list_3 = []; DP_time_list_3 = [];

    av_cost_list = []; av_time_list = [];

    for i = 2:max_cut_num
        [ DP_cost_1, DP_time_1, av_cost, av_time ] = DP_cut_domain( cutCostsFile, i, search_domain_1);
        [ DP_cost_2, DP_time_2, av_cost, av_time ] = DP_cut_domain( cutCostsFile, i, search_domain_2);
        [ DP_cost_3, DP_time_3, av_cost, av_time ] = DP_cut_domain( cutCostsFile, i, search_domain_3);

        DP_cost_list_1(i, :) = DP_cost_1;
        DP_time_list_1(i, :) = DP_time_1 * 1000;

        DP_cost_list_2(i, :) = DP_cost_2;
        DP_time_list_2(i, :) = DP_time_2 * 1000;

        DP_cost_list_3(i, :) = DP_cost_3;
        DP_time_list_3(i, :) = DP_time_3 * 1000;

        av_cost_list(i, :) = av_cost;
        av_time_list(i, :) = av_time * 1000;

    end

%     % for local min cut
%     local_cost_list = []; local_time_list = []; cut_num_list = [];
%     i = 2;
%     exp_num = 2;
%     for l = max_local_dis:-1:1
%         [ local_cost, local_time, cut_num ] = localMinTest( cutCostsFile, l);
%
%         if cut_num >= exp_num
%             local_cost_list(i, :) = local_cost;
%             local_time_list(i, :) = local_time * 1000;
%             cut_num_list(i, :) = cut_num;
%             exp_num = exp_num + 1;
%             i = i + 1;
%         end
%
%         if cut_num >= max_cut_num
%            break;
%         end
%
%     end

    % optimal cut costs, but not unniform obviously
    opti_cost_list = [];
    for i = 2:max_cut_num
       opti_cost_list(i, :) = get_the_optimal(cutCostsFile, i);
    end

    figure;
    plot(DP_cost_list_1, 'g-', 'LineWidth', 2);
    hold on;
    plot(DP_cost_list_2, 'b-', 'LineWidth', 2);
    hold on;
    plot(DP_cost_list_3, 'r-', 'LineWidth', 2);
    hold on;
    plot(av_cost_list, 'k-', 'LineWidth', 2);
    hold on;
    plot(opti_cost_list, 'm-', 'LineWidth', 2)
%     plot(cut_num_list, local_cost_list, 'm-', 'LineWidth', 2);
    grid on;
    axis([2 max_cut_num 0 max(av_cost_list)]);
    xlabel('Section Num', 'FontWeight', 'bold', 'FontSize', 12);
    ylabel('Cut Cost', 'FontWeight', 'bold', 'FontSize', 12);
%     legend('DP 5', 'DP 10', 'DP 30', 'Uniform Cutting', 'Local Minimum');
    legend('DP 5', 'DP 10', 'DP 30', 'Uniform Cutting', 'Minest');

%     figure;
%     plot(DP_time_list_1, 'g-', 'LineWidth', 2);
%     hold on;
%     plot(DP_time_list_2, 'b-', 'LineWidth', 2);
%     hold on;
%     plot(DP_time_list_3, 'r-', 'LineWidth', 2);
%     hold on;
%     plot(av_time_list, 'k-', 'LineWidth', 2);
%     hold on;
% %     plot(cut_num_list, local_time_list, 'c-', 'LineWidth', 2);
%     grid on;
%     axis([2 100 0 2]);
%     xlabel('Section Num', 'FontWeight', 'bold', 'FontSize', 12);
%     ylabel('Time Cost (ms)', 'FontWeight', 'bold', 'FontSize', 12);
% %     legend('Dynamic Programming', 'Uniform Cutting', 'Local Minimum');
%     legend('DP 5', 'DP 10', 'DP 20', 'Uniform Cutting');

end


text:
function [outputs] = mhe(inputs)

persistent first s b

if isempty(first)
    % initialization
    s = 'http://byu.apmonitor.com';
    b = 'mhe';
    msg = mhe_init(s,b)
    addpath('apm')
end

T_meas = inputs(1);
Q1 = inputs(2);

% input measurements
apm_meas(s,b,'TC',T_meas);
apm_meas(s,b,'Q1',Q1);

% solve MPC
output = apm(s,b,'solve');

% test for successful solution
if (apm_tag(s,b,'nlc.appstatus')==1)
    % retrieve the parameter values

    % option %2 retrieval (apm_tag)
    Kp = apm_tag(s,b,'Kp.Newval');
    tau = apm_tag(s,b,'tau.Newval');
    zeta = apm_tag(s,b,'zeta.Newval');
    TC_ss = apm_tag(s,b,'TC_ss.Newval');
    TC = apm_tag(s,b,'TC.Model');
else
    % display output for debugging
    disp(output)
    % not successful, set default parameters
    Kp = 0.359806899452;
    tau = 47.7311112863;
    zeta = 1.56206738412;
    TC_ss = 23.0;
    TC = 23.0;
end
params(1) = TC;
params(2) = Kp;
params(3) = tau;
params(4) = zeta;
params(5) = TC_ss;

outputs = params;

if isempty(first)
    apm_web(s,b);
    first = false;
end

return

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function pass = test_size(pref)
%TEST_SIZE     Test the SIZE method of the CHEBMATRIX class
%%
% Create some CHEBMATRIX objects and check their sizes:
x = chebfun(@(x) x);
D = operatorBlock.diff();
S = functionalBlock.sum();

ff = chebmatrix([x, x, x]);
A = [[D, D, D]; [S, S, S]];
B = [[D x]; [S 1]];

pass(1) = all(size(ff) == [1 3]);
pass(2) = all(size(ff') == [3 1]);
pass(3) = all(size(A) == [2 3]);
pass(4) = all(size(A') == [3 2]);
pass(5) = all(size(B) == [2 2]);
pass(6) = all(size(B) == [2 2]);
%%
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%EXPORT_FIG  Exports figures suitable for publication
%
% Examples:
%   im = export_fig
%   [im alpha] = export_fig
%   export_fig filename
%   export_fig filename -format1 -format2
%   export_fig ... -nocrop
%   export_fig ... -transparent
%   export_fig ... -native
%   export_fig ... -m<val>
%   export_fig ... -r<val>
%   export_fig ... -a<val>
%   export_fig ... -q<val>
%   export_fig ... -<renderer>
%   export_fig ... -<colorspace>
%   export_fig ... -append
%   export_fig ... -bookmark
%   export_fig(..., handle)
%
% This function saves a figure or single axes to one or more vector and/or
% bitmap file formats, and/or outputs a rasterized version to the
% workspace, with the following properties:
%   - Figure/axes reproduced as it appears on screen
%   - Cropped borders (optional)
%   - Embedded fonts (vector formats)
%   - Improved line and grid line styles
%   - Anti-aliased graphics (bitmap formats)
%   - Render images at native resolution (optional for bitmap formats)
%   - Transparent background supported (pdf, eps, png)
%   - Semi-transparent patch objects supported (png only)
%   - RGB, CMYK or grayscale output (CMYK only with pdf, eps, tiff)
%   - Variable image compression, including lossless (pdf, eps, jpg)
%   - Optionally append to file (pdf, tiff)
%   - Vector formats: pdf, eps
%   - Bitmap formats: png, tiff, jpg, bmp, export to workspace
%
% This function is especially suited to exporting figures for use in
% publications and presentations, because of the high quality and
% portability of media produced.
%
% Note that the background color and figure dimensions are reproduced
% (the latter approximately, and ignoring cropping & magnification) in the
% output file. For transparent background (and semi-transparent patch
% objects), use the -transparent option or set the figure 'Color' property
% to 'none'. To make axes transparent set the axes 'Color' property to
% 'none'. Pdf, eps and png are the only file formats to support a
% transparent background, whilst the png format alone supports transparency
% of patch objects.
%
% The choice of renderer (opengl, zbuffer or painters) has a large impact
% on the quality of output. Whilst the default value (opengl for bitmaps,
% painters for vector formats) generally gives good results, if you aren't
% satisfied then try another renderer.  Notes: 1) For vector formats (eps,
% pdf), only painters generates vector graphics. 2) For bitmaps, only
% opengl can render transparent patch objects correctly. 3) For bitmaps,
% only painters will correctly scale line dash and dot lengths when
% magnifying or anti-aliasing. 4) Fonts may be substitued with Courier when
% using painters.
%
% When exporting to vector format (pdf & eps) and bitmap format using the
% painters renderer, this function requires that ghostscript is installed
% on your system. You can download this from:
%   http://www.ghostscript.com
% When exporting to eps it additionally requires pdftops, from the Xpdf
% suite of functions. You can download this from:
%   http://www.foolabs.com/xpdf
%
%IN:
%   filename - string containing the name (optionally including full or
%              relative path) of the file the figure is to be saved as. If
%              a path is not specified, the figure is saved in the current
%              directory. If no name and no output arguments are specified,
%              the default name, 'export_fig_out', is used. If neither a
%              file extension nor a format are specified, a ".png" is added
%              and the figure saved in that format.
%   -format1, -format2, etc. - strings containing the extensions of the
%                              file formats the figure is to be saved as.
%                              Valid options are: '-pdf', '-eps', '-png',
%                              '-tif', '-jpg' and '-bmp'. All combinations
%                              of formats are valid.
%   -nocrop - option indicating that the borders of the output are not to
%             be cropped.
%   -transparent - option indicating that the figure background is to be
%                  made transparent (png, pdf and eps output only).
%   -m<val> - option where val indicates the factor to magnify the
%             on-screen figure pixel dimensions by when generating bitmap
%             outputs. Default: '-m1'.
%   -r<val> - option val indicates the resolution (in pixels per inch) to
%             export bitmap and vector outputs at, keeping the dimensions
%             of the on-screen figure. Default: '-r864' (for vector output
%             only). Note that the -m option overides the -r option for
%             bitmap outputs only.
%   -native - option indicating that the output resolution (when outputting
%             a bitmap format) should be such that the vertical resolution
%             of the first suitable image found in the figure is at the
%             native resolution of that image. To specify a particular
%             image to use, give it the tag 'export_fig_native'. Notes:
%             This overrides any value set with the -m and -r options. It
%             also assumes that the image is displayed front-to-parallel
%             with the screen. The output resolution is approximate and
%             should not be relied upon. Anti-aliasing can have adverse
%             effects on image quality (disable with the -a1 option).
%   -a1, -a2, -a3, -a4 - option indicating the amount of anti-aliasing to
%                        use for bitmap outputs. '-a1' means no anti-
%                        aliasing; '-a4' is the maximum amount (default).
%   -<renderer> - option to force a particular renderer (painters, opengl
%                 or zbuffer) to be used over the default: opengl for
%                 bitmaps; painters for vector formats.
%   -<colorspace> - option indicating which colorspace color figures should
%                   be saved in: RGB (default), CMYK or gray. CMYK is only
%                   supported in pdf, eps and tiff output.
%   -q<val> - option to vary bitmap image quality (in pdf, eps and jpg
%             files only).  Larger val, in the range 0-100, gives higher
%             quality/lower compression. val > 100 gives lossless
%             compression. Default: '-q95' for jpg, ghostscript prepress
%             default for pdf & eps. Note: lossless compression can
%             sometimes give a smaller file size than the default lossy
%             compression, depending on the type of images.
%   -append - option indicating that if the file (pdfs only) already
%             exists, the figure is to be appended as a new page, instead
%             of being overwritten (default).
%   -bookmark - option to indicate that a bookmark with the name of the
%               figure is to be created in the output file (pdf only).
%   handle - The handle of the figure, axes or uipanels (can be an array of
%            handles, but the objects must be in the same figure) to be
%            saved. Default: gcf.
%
%OUT:
%   im - MxNxC uint8 image array of the figure.
%   alpha - MxN single array of alphamatte values in range [0,1], for the
%           case when the background is transparent.
%
%   Some helpful examples and tips can be found at:
%      http://sites.google.com/site/oliverwoodford/software/export_fig
%
%   See also PRINT, SAVEAS.

% Copyright (C) Oliver Woodford 2008-2014

% The idea of using ghostscript is inspired by Peder Axensten's SAVEFIG
% (fex id: 10889) which is itself inspired by EPS2PDF (fex id: 5782).
% The idea for using pdftops came from the MATLAB newsgroup (id: 168171).
% The idea of editing the EPS file to change line styles comes from Jiro
% Doke's FIXPSLINESTYLE (fex id: 17928).
% The idea of changing dash length with line width came from comments on
% fex id: 5743, but the implementation is mine :)
% The idea of anti-aliasing bitmaps came from Anders Brun's MYAA (fex id:
% 20979).
% The idea of appending figures in pdfs came from Matt C in comments on the
% FEX (id: 23629)

% Thanks to Roland Martin for pointing out the colour MATLAB
% bug/feature with colorbar axes and transparent backgrounds.
% Thanks also to Andrew Matthews for describing a bug to do with the figure
% size changing in -nodisplay mode. I couldn't reproduce it, but included a
% fix anyway.
% Thanks to Tammy Threadgill for reporting a bug where an axes is not
% isolated from gui objects.

% 23/02/12: Ensure that axes limits don't change during printing
% 14/03/12: Fix bug in fixing the axes limits (thanks to Tobias Lamour for
%           reporting it).
% 02/05/12: Incorporate patch of Petr Nechaev (many thanks), enabling
%           bookmarking of figures in pdf files.
% 09/05/12: Incorporate patch of Arcelia Arrieta (many thanks), to keep
%           tick marks fixed.
% 12/12/12: Add support for isolating uipanels. Thanks to michael for
%           suggesting it.
% 25/09/13: Add support for changing resolution in vector formats. Thanks
%           to Jan Jaap Meijer for suggesting it.
% 07/05/14: Add support for '~' at start of path. Thanks to Sally Warner
%           for suggesting it.

function [im, alpha] = export_fig(varargin)
% Make sure the figure is rendered correctly _now_ so that properties like
% axes limits are up-to-date.
drawnow;
% Parse the input arguments
[fig, options] = parse_args(nargout, varargin{:});
% Isolate the subplot, if it is one
cls = all(ismember(get(fig, 'Type'), {'axes', 'uipanel'}));
if cls
    % Given handles of one or more axes, so isolate them from the rest
    fig = isolate_axes(fig);
else
    % Check we have a figure
    if ~isequal(get(fig, 'Type'), 'figure');
        error('Handle must be that of a figure, axes or uipanel');
    end
    % Get the old InvertHardcopy mode
    old_mode = get(fig, 'InvertHardcopy');
end
% Hack the font units where necessary (due to a font rendering bug in
% print?). This may not work perfectly in all cases. Also it can change the
% figure layout if reverted, so use a copy.
magnify = options.magnify * options.aa_factor;
if isbitmap(options) && magnify ~= 1
    fontu = findobj(fig, 'FontUnits', 'normalized');
    if ~isempty(fontu)
        % Some normalized font units found
        if ~cls
            fig = copyfig(fig);
            set(fig, 'Visible', 'off');
            fontu = findobj(fig, 'FontUnits', 'normalized');
            cls = true;
        end
        set(fontu, 'FontUnits', 'points');
    end
end
% MATLAB "feature": axes limits and tick marks can change when printing
Hlims = findall(fig, 'Type', 'axes');
if ~cls
    % Record the old axes limit and tick modes
    Xlims = make_cell(get(Hlims, 'XLimMode'));
    Ylims = make_cell(get(Hlims, 'YLimMode'));
    Zlims = make_cell(get(Hlims, 'ZLimMode'));
    Xtick = make_cell(get(Hlims, 'XTickMode'));
    Ytick = make_cell(get(Hlims, 'YTickMode'));
    Ztick = make_cell(get(Hlims, 'ZTickMode'));
end
% Set all axes limit and tick modes to manual, so the limits and ticks can't change
set(Hlims, 'XLimMode', 'manual', 'YLimMode', 'manual', 'ZLimMode', 'manual', 'XTickMode', 'manual', 'YTickMode', 'manual', 'ZTickMode', 'manual');
% Set to print exactly what is there
set(fig, 'InvertHardcopy', 'off');
% Set the renderer
switch options.renderer
    case 1
        renderer = '-opengl';
    case 2
        renderer = '-zbuffer';
    case 3
        renderer = '-painters';
    otherwise
        renderer = '-opengl'; % Default for bitmaps
end
% Do the bitmap formats first
if isbitmap(options)
    % Get the background colour
    if options.transparent && (options.png || options.alpha)
        % Get out an alpha channel
        % MATLAB "feature": black colorbar axes can change to white and vice versa!
        hCB = findobj(fig, 'Type', 'axes', 'Tag', 'Colorbar');
        if isempty(hCB)
            yCol = [];
            xCol = [];
        else
            yCol = get(hCB, 'YColor');
            xCol = get(hCB, 'XColor');
            if iscell(yCol)
                yCol = cell2mat(yCol);
                xCol = cell2mat(xCol);
            end
            yCol = sum(yCol, 2);
            xCol = sum(xCol, 2);
        end
        % MATLAB "feature": apparently figure size can change when changing
        % colour in -nodisplay mode
        pos = get(fig, 'Position');
        % Set the background colour to black, and set size in case it was
        % changed internally
        tcol = get(fig, 'Color');
        set(fig, 'Color', 'k', 'Position', pos);
        % Correct the colorbar axes colours
        set(hCB(yCol==0), 'YColor', [0 0 0]);
        set(hCB(xCol==0), 'XColor', [0 0 0]);
        % Print large version to array
        B = print2array(fig, magnify, renderer);
        % Downscale the image
        B = downsize(single(B), options.aa_factor);
        % Set background to white (and set size)
        set(fig, 'Color', 'w', 'Position', pos);
        % Correct the colorbar axes colours
        set(hCB(yCol==3), 'YColor', [1 1 1]);
        set(hCB(xCol==3), 'XColor', [1 1 1]);
        % Print large version to array
        A = print2array(fig, magnify, renderer);
        % Downscale the image
        A = downsize(single(A), options.aa_factor);
        % Set the background colour (and size) back to normal
        set(fig, 'Color', tcol, 'Position', pos);
        % Compute the alpha map
        alpha = round(sum(B - A, 3)) / (255 * 3) + 1;
        A = alpha;
        A(A==0) = 1;
        A = B ./ A(:,:,[1 1 1]);
        clear B
        % Convert to greyscale
        if options.colourspace == 2
            A = rgb2grey(A);
        end
        A = uint8(A);
        % Crop the background
        if options.crop
            [alpha, v] = crop_background(alpha, 0);
            A = A(v(1):v(2),v(3):v(4),:);
        end
        if options.png
            % Compute the resolution
            res = options.magnify * get(0, 'ScreenPixelsPerInch') / 25.4e-3;
            % Save the png
            imwrite(A, [options.name '.png'], 'Alpha', double(alpha), 'ResolutionUnit', 'meter', 'XResolution', res, 'YResolution', res);
            % Clear the png bit
            options.png = false;
        end
        % Return only one channel for greyscale
        if isbitmap(options)
            A = check_greyscale(A);
        end
        if options.alpha
            % Store the image
            im = A;
            % Clear the alpha bit
            options.alpha = false;
        end
        % Get the non-alpha image
        if isbitmap(options)
            alph = alpha(:,:,ones(1, size(A, 3)));
            A = uint8(single(A) .* alph + 255 * (1 - alph));
            clear alph
        end
        if options.im
            % Store the new image
            im = A;
        end
    else
        % Print large version to array
        if options.transparent
            % MATLAB "feature": apparently figure size can change when changing
            % colour in -nodisplay mode
            pos = get(fig, 'Position');
            tcol = get(fig, 'Color');
            set(fig, 'Color', 'w', 'Position', pos);
            A = print2array(fig, magnify, renderer);
            set(fig, 'Color', tcol, 'Position', pos);
            tcol = 255;
        else
            [A, tcol] = print2array(fig, magnify, renderer);
        end
        % Crop the background
        if options.crop
            A = crop_background(A, tcol);
        end
        % Downscale the image
        A = downsize(A, options.aa_factor);
        if options.colourspace == 2
            % Convert to greyscale
            A = rgb2grey(A);
        else
            % Return only one channel for greyscale
            A = check_greyscale(A);
        end
        % Outputs
        if options.im
            im = A;
        end
        if options.alpha
            im = A;
            alpha = zeros(size(A, 1), size(A, 2), 'single');
        end
    end
    % Save the images
    if options.png
        res = options.magnify * get(0, 'ScreenPixelsPerInch') / 25.4e-3;
        imwrite(A, [options.name '.png'], 'ResolutionUnit', 'meter', 'XResolution', res, 'YResolution', res);
    end
    if options.bmp
        imwrite(A, [options.name '.bmp']);
    end
    % Save jpeg with given quality
    if options.jpg
        quality = options.quality;
        if isempty(quality)
            quality = 95;
        end
        if quality > 100
            imwrite(A, [options.name '.jpg'], 'Mode', 'lossless');
        else
            imwrite(A, [options.name '.jpg'], 'Quality', quality);
        end
    end
    % Save tif images in cmyk if wanted (and possible)
    if options.tif
        if options.colourspace == 1 && size(A, 3) == 3
            A = double(255 - A);
            K = min(A, [], 3);
            K_ = 255 ./ max(255 - K, 1);
            C = (A(:,:,1) - K) .* K_;
            M = (A(:,:,2) - K) .* K_;
            Y = (A(:,:,3) - K) .* K_;
            A = uint8(cat(3, C, M, Y, K));
            clear C M Y K K_
        end
        append_mode = {'overwrite', 'append'};
        imwrite(A, [options.name '.tif'], 'Resolution', options.magnify*get(0, 'ScreenPixelsPerInch'), 'WriteMode', append_mode{options.append+1});
    end
end
% Now do the vector formats
if isvector(options)
    % Set the default renderer to painters
    if ~options.renderer
        renderer = '-painters';
    end
    % Generate some filenames
    tmp_nam = [tempname '.eps'];
    if options.pdf
        pdf_nam = [options.name '.pdf'];
    else
        pdf_nam = [tempname '.pdf'];
    end
    % Generate the options for print
    p2eArgs = {renderer, sprintf('-r%d', options.resolution)};
    if options.colourspace == 1
        p2eArgs = [p2eArgs {'-cmyk'}];
    end
    if ~options.crop
        p2eArgs = [p2eArgs {'-loose'}];
    end
    try
        % Generate an eps
        print2eps(tmp_nam, fig, p2eArgs{:});
        % Remove the background, if desired
        if options.transparent && ~isequal(get(fig, 'Color'), 'none')
            eps_remove_background(tmp_nam, 1 + using_hg2(fig));
        end
        % Add a bookmark to the PDF if desired
        if options.bookmark
            fig_nam = get(fig, 'Name');
            if isempty(fig_nam)
                warning('export_fig:EmptyBookmark', 'Bookmark requested for figure with no name. Bookmark will be empty.');
            end
            add_bookmark(tmp_nam, fig_nam);
        end
        % Generate a pdf
        eps2pdf(tmp_nam, pdf_nam, 1, options.append, options.colourspace==2, options.quality);
    catch ex
        % Delete the eps
        delete(tmp_nam);
        rethrow(ex);
    end
    % Delete the eps
    delete(tmp_nam);
    if options.eps
        try
            % Generate an eps from the pdf
            pdf2eps(pdf_nam, [options.name '.eps']);
        catch ex
            if ~options.pdf
                % Delete the pdf
                delete(pdf_nam);
            end
            rethrow(ex);
        end
        if ~options.pdf
            % Delete the pdf
            delete(pdf_nam);
        end
    end
end
if cls
    % Close the created figure
    close(fig);
else
    % Reset the hardcopy mode
    set(fig, 'InvertHardcopy', old_mode);
    % Reset the axes limit and tick modes
    for a = 1:numel(Hlims)
        set(Hlims(a), 'XLimMode', Xlims{a}, 'YLimMode', Ylims{a}, 'ZLimMode', Zlims{a}, 'XTickMode', Xtick{a}, 'YTickMode', Ytick{a}, 'ZTickMode', Ztick{a});
    end
end
return

function [fig, options] = parse_args(nout, varargin)
% Parse the input arguments
% Set the defaults
fig = get(0, 'CurrentFigure');
options = struct('name', 'export_fig_out', ...
                 'crop', true, ...
                 'transparent', false, ...
                 'renderer', 0, ... % 0: default, 1: OpenGL, 2: ZBuffer, 3: Painters
                 'pdf', false, ...
                 'eps', false, ...
                 'png', false, ...
                 'tif', false, ...
                 'jpg', false, ...
                 'bmp', false, ...
                 'colourspace', 0, ... % 0: RGB/gray, 1: CMYK, 2: gray
                 'append', false, ...
                 'im', nout == 1, ...
                 'alpha', nout == 2, ...
                 'aa_factor', 0, ...
                 'magnify', [], ...
                 'resolution', [], ...
                 'bookmark', false, ...
                 'quality', []);
native = false; % Set resolution to native of an image

% Go through the other arguments
for a = 1:nargin-1
    if all(ishandle(varargin{a}))
        fig = varargin{a};
    elseif ischar(varargin{a}) && ~isempty(varargin{a})
        if varargin{a}(1) == '-'
            switch lower(varargin{a}(2:end))
                case 'nocrop'
                    options.crop = false;
                case {'trans', 'transparent'}
                    options.transparent = true;
                case 'opengl'
                    options.renderer = 1;
                case 'zbuffer'
                    options.renderer = 2;
                case 'painters'
                    options.renderer = 3;
                case 'pdf'
                    options.pdf = true;
                case 'eps'
                    options.eps = true;
                case 'png'
                    options.png = true;
                case {'tif', 'tiff'}
                    options.tif = true;
                case {'jpg', 'jpeg'}
                    options.jpg = true;
                case 'bmp'
                    options.bmp = true;
                case 'rgb'
                    options.colourspace = 0;
                case 'cmyk'
                    options.colourspace = 1;
                case {'gray', 'grey'}
                    options.colourspace = 2;
                case {'a1', 'a2', 'a3', 'a4'}
                    options.aa_factor = str2double(varargin{a}(3));
                case 'append'
                    options.append = true;
                case 'bookmark'
                    options.bookmark = true;
                case 'native'
                    native = true;
                otherwise
                    val = str2double(regexp(varargin{a}, '(?<=-(m|M|r|R|q|Q))(\d*\.)?\d+(e-?\d+)?', 'match'));
                    if ~isscalar(val)
                        error('option %s not recognised', varargin{a});
                    end
                    switch lower(varargin{a}(2))
                        case 'm'
                            options.magnify = val;
                        case 'r'
                            options.resolution = val;
                        case 'q'
                            options.quality = max(val, 0);
                    end
            end
        else
            [p, options.name, ext] = fileparts(varargin{a});
            if ~isempty(p)
                options.name = [p filesep options.name];
            end
            switch lower(ext)
                case {'.tif', '.tiff'}
                    options.tif = true;
                case {'.jpg', '.jpeg'}
                    options.jpg = true;
                case '.png'
                    options.png = true;
                case '.bmp'
                    options.bmp = true;
                case '.eps'
                    options.eps = true;
                case '.pdf'
                    options.pdf = true;
                otherwise
                    options.name = varargin{a};
            end
        end
    end
end

% Set default anti-aliasing now we know the renderer
if options.aa_factor == 0
    options.aa_factor = 1 + 2 * (~using_hg2(fig) | (options.renderer == 3));
end

% Convert user dir '~' to full path
if numel(options.name) > 2 && options.name(1) == '~' && (options.name(2) == '/' || options.name(2) == '\')
    options.name = fullfile(char(java.lang.System.getProperty('user.home')), options.name(2:end));
end

% Compute the magnification and resolution
if isempty(options.magnify)
    if isempty(options.resolution)
        options.magnify = 1;
        options.resolution = 864;
    else
        options.magnify = options.resolution ./ get(0, 'ScreenPixelsPerInch');
    end
elseif isempty(options.resolution)
    options.resolution = 864;
end

% Check we have a figure handle
if isempty(fig)
    error('No figure found');
end

% Set the default format
if ~isvector(options) && ~isbitmap(options)
    options.png = true;
end

% Check whether transparent background is wanted (old way)
if isequal(get(ancestor(fig(1), 'figure'), 'Color'), 'none')
    options.transparent = true;
end

% If requested, set the resolution to the native vertical resolution of the
% first suitable image found
if native && isbitmap(options)
    % Find a suitable image
    list = findobj(fig, 'Type', 'image', 'Tag', 'export_fig_native');
    if isempty(list)
        list = findobj(fig, 'Type', 'image', 'Visible', 'on');
    end
    for hIm = list(:)'
        % Check height is >= 2
        height = size(get(hIm, 'CData'), 1);
        if height < 2
            continue
        end
        % Account for the image filling only part of the axes, or vice
        % versa
        yl = get(hIm, 'YData');
        if isscalar(yl)
            yl = [yl(1)-0.5 yl(1)+height+0.5];
        else
            if ~diff(yl)
                continue
            end
            yl = yl + [-0.5 0.5] * (diff(yl) / (height - 1));
        end
        hAx = get(hIm, 'Parent');
        yl2 = get(hAx, 'YLim');
        % Find the pixel height of the axes
        oldUnits = get(hAx, 'Units');
        set(hAx, 'Units', 'pixels');
        pos = get(hAx, 'Position');
        set(hAx, 'Units', oldUnits);
        if ~pos(4)
            continue
        end
        % Found a suitable image
        % Account for stretch-to-fill being disabled
        pbar = get(hAx, 'PlotBoxAspectRatio');
        pos = min(pos(4), pbar(2)*pos(3)/pbar(1));
        % Set the magnification to give native resolution
        options.magnify = (height * diff(yl2)) / (pos * diff(yl));
        break
    end
end
return

function A = downsize(A, factor)
% Downsample an image
if factor == 1
    % Nothing to do
    return
end
try
    % Faster, but requires image processing toolbox
    A = imresize(A, 1/factor, 'bilinear');
catch
    % No image processing toolbox - resize manually
    % Lowpass filter - use Gaussian as is separable, so faster
    % Compute the 1d Gaussian filter
    filt = (-factor-1:factor+1) / (factor * 0.6);
    filt = exp(-filt .* filt);
    % Normalize the filter
    filt = single(filt / sum(filt));
    % Filter the image
    padding = floor(numel(filt) / 2);
    for a = 1:size(A, 3)
        A(:,:,a) = conv2(filt, filt', single(A([ones(1, padding) 1:end repmat(end, 1, padding)],[ones(1, padding) 1:end repmat(end, 1, padding)],a)), 'valid');
    end
    % Subsample
    A = A(1+floor(mod(end-1, factor)/2):factor:end,1+floor(mod(end-1, factor)/2):factor:end,:);
end
return

function A = rgb2grey(A)
A = cast(reshape(reshape(single(A), [], 3) * single([0.299; 0.587; 0.114]), size(A, 1), size(A, 2)), class(A));
return

function A = check_greyscale(A)
% Check if the image is greyscale
if size(A, 3) == 3 && ...
        all(reshape(A(:,:,1) == A(:,:,2), [], 1)) && ...
        all(reshape(A(:,:,2) == A(:,:,3), [], 1))
    A = A(:,:,1); % Save only one channel for 8-bit output
end
return

function [A, v] = crop_background(A, bcol)
% Map the foreground pixels
[h, w, c] = size(A);
if isscalar(bcol) && c > 1
    bcol = bcol(ones(1, c));
end
bail = false;
for l = 1:w
    for a = 1:c
        if ~all(A(:,l,a) == bcol(a))
            bail = true;
            break;
        end
    end
    if bail
        break;
    end
end
bail = false;
for r = w:-1:l
    for a = 1:c
        if ~all(A(:,r,a) == bcol(a))
            bail = true;
            break;
        end
    end
    if bail
        break;
    end
end
bail = false;
for t = 1:h
    for a = 1:c
        if ~all(A(t,:,a) == bcol(a))
            bail = true;
            break;
        end
    end
    if bail
        break;
    end
end
bail = false;
for b = h:-1:t
    for a = 1:c
        if ~all(A(b,:,a) == bcol(a))
            bail = true;
            break;
        end
    end
    if bail
        break;
    end
end
% Crop the background, leaving one boundary pixel to avoid bleeding on
% resize
v = [max(t-1, 1) min(b+1, h) max(l-1, 1) min(r+1, w)];
A = A(v(1):v(2),v(3):v(4),:);
return

function eps_remove_background(fname, count)
% Remove the background of an eps file
% Open the file
fh = fopen(fname, 'r+');
if fh == -1
    error('Not able to open file %s.', fname);
end
% Read the file line by line
while count
    % Get the next line
    l = fgets(fh);
    if isequal(l, -1)
        break; % Quit, no rectangle found
    end
    % Check if the line contains the background rectangle
    if isequal(regexp(l, ' *0 +0 +\d+ +\d+ +r[fe] *[\n\r]+', 'start'), 1)
        % Set the line to whitespace and quit
        l(1:regexp(l, '[\n\r]', 'start', 'once')-1) = ' ';
        fseek(fh, -numel(l), 0);
        fprintf(fh, l);
        % Reduce the count
        count = count - 1;
    end
end
% Close the file
fclose(fh);
return

function b = isvector(options)
b = options.pdf || options.eps;
return

function b = isbitmap(options)
b = options.png || options.tif || options.jpg || options.bmp || options.im || options.alpha;
return

% Helper function
function A = make_cell(A)
if ~iscell(A)
    A = {A};
end
return

function add_bookmark(fname, bookmark_text)
% Adds a bookmark to the temporary EPS file after %%EndPageSetup
% Read in the file
fh = fopen(fname, 'r');
if fh == -1
    error('File %s not found.', fname);
end
try
    fstrm = fread(fh, '*char')';
catch ex
    fclose(fh);
    rethrow(ex);
end
fclose(fh);

% Include standard pdfmark prolog to maximize compatibility
fstrm = strrep(fstrm, '%%BeginProlog', sprintf('%%%%BeginProlog\n/pdfmark where {pop} {userdict /pdfmark /cleartomark load put} ifelse'));
% Add page bookmark
fstrm = strrep(fstrm, '%%EndPageSetup', sprintf('%%%%EndPageSetup\n[ /Title (%s) /OUT pdfmark',bookmark_text));

% Write out the updated file
fh = fopen(fname, 'w');
if fh == -1
    error('Unable to open %s for writing.', fname);
end
try
    fwrite(fh, fstrm, 'char*1');
catch ex
    fclose(fh);
    rethrow(ex);
end
fclose(fh);
return

author: Cloud-CV

repo: object-proposals


text:
%GETHANDLEFG Framegrabber device gethandle function.
%
%SYNTAX devinfo = gethandlefg('alias')
%
%Description This function will return the handle of
%the framegrabber device when you lost the handle to
%the device.
%
%INPUT The alias (alias) of the framegrabber device.
%
%REMARKS It will be verified if the input is a string.
%If it's valid string, the function will  return the
%device information structure of the device.
%
%OUTPUT If the input isn't valid, some error message will
%be shown to the user. If everything's ok, the function
%will return the device information structure of the opened
%framegrabber device with the given alias.
%
%EXAMPLE To regain the handle to the framegrabber, type
%in the Matlab Command Window :

%m=gethandlefg('DT3155');

disp('Error: gethandlefg not found.')

author: Sable

repo: mcbench-benchmarks


text:
function [f,g]=idgtreal(coef,g,a,M,varargin)
%IDGTREAL  Inverse discrete Gabor transform for real-valued signals
%   Usage:  f=idgtreal(c,g,a,M);
%           f=idgtreal(c,g,a,M,Ls);
%
%   Input parameters:
%         c     : Array of coefficients.
%         g     : Window function.
%         a     : Length of time shift.
%         M     : Number of channels.
%         Ls    : length of signal.
%   Output parameters:
%         f     : Signal.
%
%   `idgtreal(c,g,a,M)` computes the Gabor expansion of the input coefficients
%   *c* with respect to the real-valued window *g*, time shift *a* and number of
%   channels *M*. *c* is assumed to be the positive frequencies of the Gabor
%   expansion of a real-valued signal.
%
%   It must hold that `size(c,1)==floor(M/2)+1`. Note that since the
%   correct number of channels cannot be deduced from the input, `idgtreal`
%   takes an additional parameter as opposed to |idgt|.
%
%   The window *g* may be a vector of numerical values, a text string or a
%   cell array. See the help of |gabwin| for more details.
%
%   `idgtreal(c,g,a,M,Ls)` does as above but cuts or extends *f* to length *Ls*.
%
%   `[f,g]=idgtreal(...)` additionally outputs the window used in the
%   transform. This is usefull if the window was generated from a description
%   in a string or cell array.
%
%   For perfect reconstruction, the window used must be a dual window of the
%   one used to generate the coefficients.
%
%   If *g* is a row vector, then the output will also be a row vector. If *c* is
%   3-dimensional, then `idgtreal` will return a matrix consisting of one column
%   vector for each of the TF-planes in *c*.
%
%   See the help on |idgt| for the precise definition of the inverse Gabor
%   transform.
%
%   `idgtreal` takes the following flags at the end of the line of input
%   arguments:
%
%     'freqinv'  Use a frequency-invariant phase. This is the default
%                convention described in the help for |dgt|.
%
%     'timeinv'  Use a time-invariant phase. This convention is typically
%                used in filter bank algorithms.
%
%   Examples:
%   ---------
%
%   The following example demostrates the basic pricinples for getting
%   perfect reconstruction (short version):::
%
%     f=greasy;            % test signal
%     a=32;                % time shift
%     M=64;                % frequency shift
%     gs={'blackman',128}; % synthesis window
%     ga={'dual',gs};      %  analysis window
%
%     [c,Ls]=dgtreal(f,ga,a,M); % analysis
%
%     % ... do interesting stuff to c at this point ...
%
%     r=idgtreal(c,gs,a,M,Ls); % synthesis
%
%     norm(f-r)                % test
%
%   The following example does the same as the previous one, with an
%   explicit construction of the analysis and synthesis windows:::
%
%     f=greasy;     % test signal
%     a=32;         % time shift
%     M=64;         % frequency shift
%     Ls=length(f); % signal length
%
%     % Length of transform to do
%     L=dgtlength(Ls,a,M);
%
%     % Analysis and synthesis window
%     gs=firwin('blackman',128);
%     ga=gabdual(gs,a,M,L);
%
%     c=dgtreal(f,ga,a,M);  % analysis
%
%     % ... do interesting stuff to c at this point ...
%
%     r=idgtreal(c,gs,a,M,Ls); % synthesis
%
%     norm(f-r)       % test
%
%   See also:  idgt, gabwin, gabdual, dwilt

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

% Check input paramameters.

if nargin<4
  error('%s: Too few input parameters.',upper(mfilename));
end;

if ~isnumeric(g) && prod(size(g))==1
  error('g must be a vector (you probably forgot to supply the window function as input parameter.)');
end;

% Define initial value for flags and key/value pairs.
definput.keyvals.Ls=[];
definput.keyvals.lt=[0 1];
definput.flags.phase={'freqinv','timeinv'};

[flags,kv,Ls]=ltfatarghelper({'Ls'},definput,varargin);

N=size(coef,2);
W=size(coef,3);

% Make a dummy call to test the input parameters
Lsmallest=dgtlength(1,a,M,kv.lt);

M2=floor(M/2)+1;

if M2~=size(coef,1)
  error('Mismatch between the specified number of channels and the size of the input coefficients.');
end;

L=N*a;

if rem(L,Lsmallest)>0
    error('%s: Invalid size of coefficient array.',upper(mfilename));
end;

if kv.lt(2)>2
  error('Only rectangular or quinqux lattices are supported.');
end;

if kv.lt(2)~=1 && flags.do_timeinv
    error(['%s: Time-invariant phase for quinqux lattice is not ',...
           'supported.'],upper(mfilename));
end

%% ----- step 3 : Determine the window

[g,info]=gabwin(g,a,M,L,kv.lt,'callfun',upper(mfilename));

if L<info.gl
  error('%s: Window is too long.',upper(mfilename));
end;

if ~isreal(g)
  error('%s: Window must be real-valued.',upper(mfilename));
end;

% Do the actual computation.
f=comp_idgtreal(coef,g,a,M,kv.lt,flags.do_timeinv);

% Cut or extend f to the correct length, if desired.
if ~isempty(kv.Ls)
  f=postpad(f,kv.Ls);
else
  kv.Ls=L;
end;

f=comp_sigreshape_post(f,Ls,0,[0; W]);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function coef=framecoef2tf(F,coef)
%FRAMECOEF2TF  Convert coefficients to time-frequency plane
%   Usage: cout=framecoef2tf(F,cin);
%
%   `framecoef2tf(F,cin)` converts the frame coefficients *cin* into the
%   time-frequency plane layout. The frame object *F* must have been
%   created using |frame|.
%
%   The time-frequency plane layout is a matrix, where the first
%   dimension indexes frequency and the second dimension time. This is
%   similar to the output format from |dgt| and |wmdct|.
%
%   Not all types of frames support this coefficient conversion. The supported
%   types of frames are: `'dgt'`, `'dgtreal'`, `'dwilt'`, `'wmdct'`, `'ufilterbank'`,
%   `'ufwt'`,`'uwfbt'` and `'uwpfbt'`.
%
%   See also: frame, frametf2coef, framecoef2native

complainif_notenoughargs(nargin,2,'FRAMECOEF2TF');
complainif_notvalidframeobj(F,'FRAMECOEF2TF');

switch(F.type)
 case 'dgt'
  [MN,W]=size(coef);
  N=MN/F.M;
  coef=reshape(coef,[F.M,N,W]);
 case 'dgtreal'
  [MN,W]=size(coef);
  M2=floor(F.M/2)+1;
  N=MN/M2;
  coef=reshape(coef,[M2,N,W]);
 case 'dwilt'
  [MN,W]=size(coef);
  N=MN/F.M;
  coef=wil2rect(reshape(coef,[2*F.M,N/2,W]));
 case 'wmdct'
  [MN,W]=size(coef);
  N=MN/F.M;
  coef=reshape(coef,[F.M,N,W]);
 case 'ufilterbank'
  [MN,W]=size(coef);
  M=numel(F.g);
  N=MN/M;
  coef=permute(reshape(coef,[N,M,W]),[2,1,3]);
 case {'ufwt','uwfbt','uwpfbt'}
  coef = permute(F.coef2native(coef,size(coef)),[2,1,3]);
 otherwise
  error('%s: TF-plane layout not supported for this transform.',upper(mfilename));
end;


text:
clf
hold on;
x = 0:0.01:8*pi;
data = [x', sin(x)'];
plot(data(:,1), data(:,2));
axis([0, 8*pi, -3, 3]);
legend('original data');
grid on;

pause;

clf
hold on;

data_length = length(data)
for idx = 1500:2000
    data(idx,2) = cos(data(idx,1));
end;
plot(data(:,1), data(:,2));
axis([0, 8*pi, -3, 3]);
legend('human modified data');
grid on;

pause;

%in this demo, we only use best_so_far_loc while we plot the figure
n = 500;
[best_so_far_dist, best_so_far_loc] = HOTSAX(data(:,2), n, 'nseg 4')

%plot the figure
hold on;
if ~isnan(best_so_far_loc)
    plot_x = data(best_so_far_loc:best_so_far_loc+n, 1);
    plot_y = data(best_so_far_loc:best_so_far_loc+n, 2);
    plot(plot_x,plot_y, 'r');
    legend('human modified data', 'anomaly detect');
end
hold off;

author: dsmi-lab-ntust

repo: AnomalyDetectionToolbox


text:
function rts_Dyn_populateAttAxis(hFig, hAttDispAxes,attData)
%rts_Dyn_populateAttAxis Summary of this function goes here
%   Detailed explanation goes here
    persistent hX hXText hY hYText hZ hZText hCone EndPlate1 EndPlate2 hXB hXBText hYB hYBText hZB hZBText az el hRates hRatesText;
    if(ishandle(hX))
        delete(hX);
    end
    if(ishandle(hXText))
        delete(hXText);
    end
    if(ishandle(hY))
        delete(hY);
    end
    if(ishandle(hYText))
        delete(hYText);
    end
    if(ishandle(hZ))
        delete(hZ);
    end
    if(ishandle(hZText))
        delete(hZText);
    end
    if(ishandle(hCone))
        delete(hCone);
    end
    if(ishandle(EndPlate1))
        delete(EndPlate1);
    end
    if(ishandle(EndPlate2))
        delete(EndPlate2);
    end
    if(ishandle(hXB))
        delete(hXB);
    end
    if(ishandle(hYB))
        delete(hYB);
    end
    if(ishandle(hYBText))
        delete(hYBText);
    end
    if(ishandle(hXBText))
        delete(hXBText);
    end
    if(ishandle(hZB))
        delete(hZB);
    end
    if(ishandle(hZBText))
        delete(hZBText);
    end
    if(ishandle(hRates))
        delete(hRates);
    end
    if(ishandle(hRatesText))
        delete(hRatesText);
    end

    if(isempty(az) || isempty(el))
        az = -37.5;
        el= 30;
    else
        [az, el] = view(hAttDispAxes);
    end

    hold(hAttDispAxes,'on');
    axis(hAttDispAxes);
    [hX, hXText, hY, hYText, hZ, hZText, hCone, EndPlate1, EndPlate2, hXB, hXBText, hYB, hYBText, hZB, hZBText, hRates, hRatesText] = populateAttPlotWAtt(attData, hAttDispAxes);
    rts_initOrbPlot(hFig, hAttDispAxes, []);
    view(hAttDispAxes, az, el);
    hold(hAttDispAxes,'off');
end

function [hX, hXText, hY, hYText, hZ, hZText, hCone,EndPlate1,EndPlate2, hXB, hXBText, hYB, hYBText, hZB, hZBText, hRates, hRatesText] =  populateAttPlotWAtt(attData, hAttDispAxes)
    e1 = attData(16);
    e2 = attData(17);
    e3 = attData(18);
    e4 = attData(19);

    xRate = rad2deg(attData(28));
    yRate = rad2deg(attData(29));
    zRate = rad2deg(attData(30));
    rateVect = [xRate, yRate, zRate];
    if(norm(rateVect) > 0.0)
        rateVect = rateVect / norm(rateVect);
    else
        rateVect = [0,0,0];
    end

    C(1,1)=1 - 2*e2^2 - 2*e3^2;
    C(1,2)=2*(e1*e2 - e3*e4);
    C(1,3)=2*(e3*e1 + e2*e4);
    C(2,1)=2*(e1*e2 + e3*e4);
    C(2,2)=1 - 2*e3^2 - 2*e1^2;
    C(2,3)=2*(e2*e3 - e1*e4);
    C(3,1)=2*(e3*e1 - e2*e4);
    C(3,2)=2*(e2*e3 + e1*e4);
    C(3,3)=1 - 2*e1^2 - 2*e2^2;

    xRot = (C)*[1;0;0];
    xRot = xRot/norm(xRot);

    zRot = (C)*[0;0;1];
    zRot = (zRot/norm(zRot));

    R1 = [cos(pi/2) -sin(pi/2) 0;
          sin(pi/2) cos(pi/2) 0;
          0 0 1];
    R2 = [1 0 0;
          0 cos(pi/2) -sin(pi/2);
          0 sin(pi/2) cos(pi/2)];
    R3 = [cos(pi/2) 0 sin(pi/2);
          0 1 0;
          -sin(pi/2) 0 cos(pi/2)];

    R2 = R2(:,[1:2-1,3,2+1:3-1,2,3+1:end]);
    R3 = R2(:,[1:2-1,3,2+1:3-1,2,3+1:end]);
	R = R1*R2*R3;

    xRot = R*xRot;
    zRot = R*zRot;

    yRot = cross(zRot, xRot);
    yRot = -yRot/norm(yRot);

    zRotNorm = 0.7*zRot/norm(zRot);

    yRot =0.7*yRot;
    xRot =0.7*xRot;

    X1 = [0 0 0]';
    X2 = 0.5*zRot;
    R = [0.25 0.05];
    n = 30;
    cyl_color = 'r';
    closed = 1;
    lines = 0;
    [hCone,EndPlate1,EndPlate2] = Cone(hAttDispAxes, X1,X2,R,n,cyl_color,closed,lines);

    hX = plot3(hAttDispAxes, [-1 1], [0 0], [0 0],'LineWidth',1.5);
    hXB = plot3(hAttDispAxes, [0 xRot(1)], [0 xRot(2)], [0 xRot(3)],'m','LineWidth',1.5);
    hXText = text(1.05,0,0,'x','Color','k', 'FontWeight', 'bold', 'Parent', hAttDispAxes);
    hXBText = text(xRot(1),xRot(2),xRot(3),'x_b','Color','k', 'FontWeight', 'bold', 'Parent', hAttDispAxes);

    hY = plot3(hAttDispAxes, [0 0], [-1 1], [0 0],'LineWidth',1.5);
    hYB = plot3(hAttDispAxes, [0 yRot(1)], [0 yRot(2)], [0 yRot(3)],'m','LineWidth',1.5);
    hYText = text(0,1.05,0,'y','Color','k', 'FontWeight', 'bold', 'Parent', hAttDispAxes);
    hYBText = text(yRot(1),yRot(2),yRot(3),'y_b','Color','k', 'FontWeight', 'bold', 'Parent', hAttDispAxes);

    hZ = plot3(hAttDispAxes, [0 0], [0 0], [-1 1],'LineWidth',1.5);
    hZB = plot3(hAttDispAxes, [0 zRotNorm(1)], [0 zRotNorm(2)], [0 zRotNorm(3)],'m','LineWidth',1.5);
    hZText = text(0,0,1.05,'z','Color','k', 'FontWeight', 'bold', 'Parent', hAttDispAxes);
    hZBText = text(zRotNorm(1),zRotNorm(2),zRotNorm(3),'z_b','Color','k', 'FontWeight', 'bold', 'Parent', hAttDispAxes);

    hRates = plot3(hAttDispAxes, [0 rateVect(1)], [0 rateVect(2)], [0 rateVect(3)],'g','LineWidth',1.5);
    hRatesText = text(rateVect(1),rateVect(2),rateVect(3),'h_hat','Color','k', 'FontWeight', 'bold', 'Parent', hAttDispAxes);

%     sunvect = [-0.3, 0.2, 0.1];
%     sunvect = sunvect/norm(sunvect);
%     sunVectText= 1.05*sunvect;
%     sunTextColor = [218,165,32]/256;
%     plot3([0 sunvect(1)], [0 sunvect(2)], [0 sunvect(3)],'LineWidth',1.5, 'Color',sunTextColor);
%     text(sunVectText(1), sunVectText(2), sunVectText(3),'Sun','Color',sunTextColor, 'FontWeight', 'bold');
%
%     celBodyvect = -[0.75 0.75 0.75];
%     celBodyvect = celBodyvect/norm(celBodyvect);
%     celBodyvectText= 1.05*celBodyvect;
%     celBodyTextColor = [0,100,0]/256;
%     plot3([0 celBodyvect(1)], [0 celBodyvect(2)], [0 celBodyvect(3)],'LineWidth',1.5, 'Color',celBodyTextColor);
%     text(celBodyvectText(1), celBodyvectText(2), celBodyvectText(3),'Kerbin','Color',celBodyTextColor, 'FontWeight', 'bold');

end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%PRINCSE Compute principal stresses and strains.
%
%   [ SE_P ] = PRINCSE( SE_X, SE_Y, SE_Z, SE_XY, SE_YZ, SE_XZ, ICOMP ) Computes
%   principal stresses and strains (for the ICOMP component, default 1). In two
%   dimensions only inputs SE_X, SE_Y, SE_Z, SE_XY, and ICOMP are required.

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function count = hprintf (varargin)
%HPRINTF fprintf with hypertext links not highlighted in the command window.
% hprintf does this by replacing the string "href=" with "HREF=".  If the file
% descriptor is not present, the output defaults to the command window, just
% like fprintf.  Note that all browsers accept either href= or HREF=, so your
% hypertext links will not be affected except within the MATLAB Command Window.
% General usage is identical to fprintf:
%
%      hprintf (format, arg1, arg2, ...)
%      hprintf (fid, format, arg1, arg2, ...)
%      count = hprintf ( ... same as above ... )
%
% Example:
%      fprintf ('MathWorks\n') ;
%      fprintf ('%d\n', 42) ;
%      hprintf ('MathWorks\n') ;
%      hprintf ('%d\n', 42) ;
%
% For a discussion, see Kristin's blog and the comments there at
% http://blogs.mathworks.com/desktop/2007/07/09
% (Kristin's blog).
%
% NOTE: the examples above are modified by "help hprintf" so that you cannot
% see the HREF= text.  To see the examples properly (without hypertext
% highlighting) use:
%
%      edit hprintf
%
% To try the examples above, use hprintf with no inputs (note that this usage
% of hprintf also flags an error, to exactly mimic the fprintf behavior):
%
%      hprintf
%
% Here is a slightly more complex example that has the advantage of being
% printed properly by "help hprintf":
%
%      % a string template with hypertext contents:
%      str = 'MathWorks\n' ;
%      % made into an active hypertext, which will be underlined when
%      % displayed in the command window:
%      hstr = strrep (str, 'AREF', 'href') ;
%      fprintf (hstr) ;
%      %          displays: 'MathWorks'
%      hprintf (hstr) ;
%      %          displays: 'MathWorks'
%
% See also fprintf, strrep, sprintf.

% Copyright 2007, T. Davis, with thanks to 'us' (us at neurol dot unizh dot ch)
% for suggestions.

% Aug 25, 2007

if (nargin < 1)

    % try hprintf
    help hprintf
    fprintf ('\nhypertext highlighting with fprintf:\n\n') ;
    fprintf ('MathWorks\n') ;
    fprintf ('%d\n', 42) ;
    fprintf ('\nhypertext highlighting turned off with hprintf:\n\n') ;
    hprintf ('MathWorks\n') ;
    hprintf ('%d\n\n', 42) ;

    % flag an error, to mimic fprintf behavior
    error ('Not enough input arguments') ;

elseif (nargout > 1)

    % mimic fprintf
    error ('Too many output arguments') ;

else

    if (ischar (varargin {1}))
        % mimic fprintf ('hello world %d\n', 42), with no file ID
        cnt = fprintf (strrep (sprintf (varargin {:}), 'href=', 'HREF=')) ;
    else
        % mimic fprintf (fid, 'hello world %d\n', 42), with file ID given
        cnt = fprintf (varargin {1}, ...
            strrep (sprintf (varargin {2:end}), 'href=', 'HREF=')) ;
    end
    if (nargout > 0)
        % return the fprintf output
        count = cnt ;
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function triangle01_poly_integral_test ( )

%*****************************************************************************80
%
%% TRIANGLE01_POLY_INTEGRAL_TEST: polynomial integrals over the unit triangle.
%
%  Location:
%
%    http://people.sc.fsu.edu/~jburkardt/m_src/triangle_integrals/triangle01_poly_integral_test.m
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 April 2015
%
%  Author:
%
%    John Burkardt
%
  d_max = 6;
  k_max = ( ( d_max + 1 ) * ( d_max + 2 ) ) / 2;
  qm = zeros ( k_max, 1 );
  for k = 1 : k_max
    [ i, j ] = i4_to_pascal ( k );
    qm(k) = triangle01_monomial_integral ( i, j );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGLE01_POLY_INTEGRAL_TEST\n' );
  fprintf ( 1, '  TRIANGLE01_POLY_INTEGRAL returns the integral Q of\n' );
  fprintf ( 1, '  a polynomial P(X,Y) over the interior of the unit triangle.\n' );

  d = 1;
  m = ( ( d + 1 ) * ( d + 2 ) ) / 2;
  p = [ 1.0, 2.0, 3.0 ]';
  fprintf ( 1, '\n' );
  poly_print ( d, p, '  p(x,y)' );
  q = triangle01_poly_integral ( d, p );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Q =         %g\n', q );
  q2 = p(1:m)' * qm(1:m);
  fprintf ( 1, '  Q (exact) = %g\n', q2 );

  d = 2;
  m = ( ( d + 1 ) * ( d + 2 ) ) / 2;
  p = [ 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 ]';
  fprintf ( 1, '\n' );
  poly_print ( d, p, '  p(x,y)' );
  q = triangle01_poly_integral ( d, p );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Q =         %g\n', q );
  q2 = p(1:m)' * qm(1:m);
  fprintf ( 1, '  Q (exact) = %g\n', q2 );

  d = 2;
  m = ( ( d + 1 ) * ( d + 2 ) ) / 2;
  p = [ 1.0, -2.0, 3.0, -4.0, 5.0, -6.0 ]';
  fprintf ( 1, '\n' );
  poly_print ( d, p, '  p(x,y)' );
  q = triangle01_poly_integral ( d, p );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Q =         %g\n', q );
  q2 = p(1:m)' * qm(1:m);
  fprintf ( 1, '  Q (exact) = %g\n', q2 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function eta=shocksim3(cfconds,Fperiods,N,n,fmat,ortirfcell)

% prelminary tasks: initiate R and r
R=[];
r=[];

% loop over forecast periods to build the matrix R and the vector r
for ii=1:Fperiods
   % for period ii, check if there is a condition, variable after variable
   for jj=1:N*n
      % if there is a condition
      if ~isempty(cfconds{ii,jj})
      % fill r with the corresponding condition, minus forecast..
      r=[r;cfconds{ii,jj}-fmat(ii,jj)];
      % .. and R with the corresponding orthogonalised IRFs entries
      % first increment R with a row of zeros of suitable dimension (N*n*Fperiods)
      R=[R;sparse(1,N*n*Fperiods)];
         % loop over periods up to the one on which there is the constraint
         for kk=1:ii
         R(end,(kk-1)*(N*n)+1:kk*(N*n))=ortirfcell{ii-kk+1,ii}(jj,:);
         end
      % if there is no condition, don't do anything
      end
   end
end

% once the linear system is identified, draw a full vector of shocks from the Waggoner-Zha distribution
% realise the singular value decomposition of R as in (3.3.17) and obtain the corresponding matrices
% recover Q and K, the dimensions of the matrix Rtemp
[Q,K]=size(R);
% obtain the singular value decomposition
[U,S,V]=svd(full(R));
% obtain the required matrices
P=S(:,1:Q);
V1=V(:,1:Q);
V2=V(:,Q+1:end);
% draw the vector of constrained shocks from N(etabar,gammabar)
eta=V1/P*U'*r+V2*randn(K-Q,1);

author: european-central-bank

repo: BEAR-toolbox


text:
function [l] = cl2l(cl)
% Convert volume from centiliters to liters.
% Chad Greene 2012
l = cl*0.01;


text:
function buildMFun(name,mode,sobj,svar,opts)
% BUILDMFUN Build a MATLAB function for nonlinear callbacks

%Determine callback type
switch(mode)
    case 'obj'
        title = 'Objective Function Callback';
        fcall = sprintf('function j = %s(x)\n',name);
        var = 'j';
    case 'grad'
        title = 'Objective Gradient Callback';
        fcall = sprintf('function g = %s(x)\n',name);
        var = 'g';
    case 'con'
        title = 'Constraint Function Callback';
        fcall = sprintf('function c = %s(x)\n',name);
        var = 'c';
    case 'jac'
        title = 'Constraint Jacobian Callback';
        fcall = sprintf('function J = %s(x)\n',name);
        var = 'J';
    case 'hess'
        title = 'Hessian of the Lagrangian Callback';
        fcall = sprintf('function H = %s(x,sigma,lambda)\n',name);
        var = 'H';
    otherwise
        error('Unknown MFun mode: %s',mode);
end

%Open file
fp = fopen([name '.m'],'w');
%Write Header
fprintf(fp,fcall);
fprintf(fp,'%% %s\n',upper(name));
fprintf(fp,'%%\n%% %s\n\n',title);
fprintf(fp,'%% Symbolic Builder Auto-Generated Callback Function\n');
fprintf(fp,'%% Generated %s\n\n',datestr(now));

%Distribute Vars
v = SymBuilder.detVarNames(svar);
if(~strcmp(v{1,1},'x') || size(v,1) > 1)
    fprintf(fp,'%% Slice Input Array\n');
    for i = 1:size(v,1)
        if(v{i,2}(end)-v{i,2}(1) == 0)
            fprintf(fp,'%s = x(%d);\n',v{i,1},v{i,2}(1));
        else
            fprintf(fp,'%s = x(%d:%d);\n',v{i,1},v{i,2}(1),v{i,2}(end));
        end
    end
end

%Check for options
if(nargin < 5 || ~isstruct(opts))
    preallocate = true;
else
    if(isfield(opts,'preallocate'))
        preallocate = opts.preallocate;
    end
end

%Build cell array of indexed variable strings
ivar = buildVarIndex(v,length(svar));
%Ensure both columns
if(size(svar,1) > 1), svar = svar.'; end
if(size(ivar,1) > 1), ivar = ivar'; end
%Subs out individual symbolic variables into our indexed list
eq = subs(sobj,svar,ivar);

%Get equation size (matrices treated differently)
[r,c] = size(eq);

%Enter Equations (Var Type Dictates Entry Type)
switch(var)
    case 'j' %scalar
        fprintf(fp,'\n%% Equation:\n');
        fprintf(fp,'%s(1,1) = %s;\n',var,char(eq));
    case {'g','c'} %vector, dense
        %HACK
        if(~isa(eq,'sym')), eq = sym(eq); end
        %Check for row
        if(c > 1), eq = eq.'; tr = 1; else tr = 0; end
        %Convert equations to string
        str = char(eq);
        %Remove 'matrix' if is found
        ind = strfind(str,'matrix');
        if(ind)
            str = str(ind+7:end-1);
        end
        if(preallocate)
            fprintf(fp,'\n%% Preallocate\n');
            if(tr)
                fprintf(fp,'%s = zeros(1,%d);\n',var,length(eq));
            else
                fprintf(fp,'%s = zeros(%d,1);\n',var,length(eq));
            end
        end
        if(strcmp(str(1),'[')), str = str(2:end); end %remove extra square brackets
        if(strcmp(str(end-1:end),']]')), str = str(1:end-2); end
        ss = regexp(str,'], ','split');
        fprintf(fp,'\n%% Equations:\n');
        if(tr)
            for i = 1:length(ss)
                if(~strcmp(ss{i}(2:end),'0'))
                    if(strcmp(ss{i}(1),'[')), ss{i} = ss{i}(2:end); end %remove extra square brackets
                    fprintf(fp,'%s(1,%d) = %s;\n',var,i,ss{i});
                end
            end
        else
            for i = 1:length(ss)
                if(~strcmp(ss{i}(2:end),'0'))
                    if(strcmp(ss{i}(1),'[')), ss{i} = ss{i}(2:end); end %remove extra square brackets
                    fprintf(fp,'%s(%d,1) = %s;\n',var,i,ss{i});
                end
            end
        end
    case {'J','H'} %matrix, sparse
        %Get NonZero elements from equation
        nzel = logical(eq ~= 0);
        [rows,cols] = find(sparse(nzel));
        nz = nnz(nzel);
        nzsym = eq(nzel);

        %OLD WAY (Inefficient?)
        if(preallocate)
            fprintf(fp,'\n%% Preallocate\n');
            fprintf(fp,'%s = spalloc(%d,%d,%d);\n',var,r,c,nz);
        end
        fprintf(fp,'\n%% Sparse Matrix:\n');
        for i = 1:nz
            fprintf(fp,'%s(%d,%d) = %s;\n',var,rows(i),cols(i),char(nzsym(i)));
        end

        %FASTER BUT LESS CLEAR
%         %Row & Col Indices
%         fprintf(fp,'\n%% Row & Col Indices\n');
%         fprintf(fp,'rows = [');
%         for i = 1:length(rows)
%             if(i == length(rows))
%                 fprintf(fp,'%d];\n',rows(i));
%             else
%                 fprintf(fp,'%d, ',rows(i));
%             end
%         end
%         fprintf(fp,'cols = [');
%         for i = 1:length(cols)
%             if(i == length(cols))
%                 fprintf(fp,'%d];\n',cols(i));
%             else
%                 fprintf(fp,'%d, ',cols(i));
%             end
%         end
%         fprintf(fp,'\n%% Preallocate\n');
%         fprintf(fp,'vals = zeros(%d,1);\n',nz);
%
%         fprintf(fp,'\n%% Nonzero Entries\n');
%         for i = 1:nz
%             fprintf(fp,'vals(%d,1) = %s;\n',i,char(nzsym(i)));
%         end
%
%         fprintf(fp,'\n%% Create Sparse Return Matrix\n');
%         fprintf(fp,'%s = sparse(rows,cols,vals);\n',var);
end

%Close file
fclose(fp);

%Convert variables from names to indexed strings [nominally m1 to m(1)]
function ivar = buildVarIndex(v,no)

%Create cell to store strings in
ivar = cell(no,1); n = 1;

%Process v cell matrix, building strings as we go
for i = 1:size(v,1)
    name = v{i,1};
    indices = v{i,2};
    %If single variable, easy!
    if(length(indices) == 1)
        ivar{n} = name; n = n + 1;
    %Otherwise have to run through indexing list
    else
        for j = 1:length(indices)
            ivar{n} = sprintf('%s(%d)',name,j); n = n + 1;
        end
    end
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% addpath('data');
% addpath(genpath('MSTV'));
% addpath(genpath('LPP'));
% addpath(genpath('OTVCA_V3'));
% addpath('RF classifier');

%% Load data
load Houston.mat;
HSI = Houston;
load TRLabel
load TSLabel
Tr = TRLabel;
Te = TSLabel;
clear TRLabel TSLabel Houston
%% Parameter setting
dim=max(unique(Te)); % Number of Features is set to the number of classes
Trees = 200; % for RF training
%% PCA+RF
[oa, pa, K, CM] = USFE_PCA(HSI, Tr, Te, dim, Trees);
aa = mean(pa);
%% MSTV+RF
[oa, pa, K, CM] = USFE_MSTV(HSI, Tr, Te, dim, Trees);
aa = mean(pa);
%% OTVCA Rasti et al.
[oa, pa, K, CM] = USFE_OTVCA(HSI, Tr, Te, dim, Trees);
aa = mean(pa);
%% LPP
[oa, pa, K, CM] = USFE_LPP(HSI, Tr, Te, dim, Trees);
aa = mean(pa);

author: BehnoodRasti


text:
function output = callsdplr(interfacedata)

% Retrieve needed data
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
K       = interfacedata.K;
ub      = interfacedata.ub;
lb      = interfacedata.lb;
lowrankdetails =  interfacedata.lowrankdetails;

% Create the parameter structure
pars = options.sdplr;
pars.printlevel = options.verbose;

% *********************************************
% Bounded variables converted to constraints
% N.B. Only happens when caller is BNB
% *********************************************
if ~isempty(ub)
    [F_struc,K] = addStructureBounds(F_struc,K,ub,lb);
end

% rmfield is slow...
Knew.l = K.l;
Knew.s = K.s;
K = Knew;

if options.savedebug
    save sdplrdebug F_struc c K pars -V6
end

% *********************************************
%FIND LOW RANK STRUCTURES
% *********************************************
problem = 0;
lrA = [];
if ~isempty(lowrankdetails) | (options.sdplr.maxrank>0 & (K.s(1)>0))
    showprogress('Detecting low rank data',options.showprogress);
    sdploc = K.l+cumsum([1 K.s.^2]);
    k = 1;
    % FIX : Lazy code copy...
    [ix,jx,sx] = find(F_struc);
    if isempty(lowrankdetails)
        % Okay, this means we have to go through everything...
        % Get all data in F_struc for later
        for lmiid = 1:length(K.s)
            removethese = zeros(1,size(F_struc,2)-1);
            for i = 1:size(F_struc,2)-1
                Fi = reshape(F_struc(sdploc(lmiid):sdploc(lmiid+1)-1,i+1),K.s(lmiid),K.s(lmiid));
                if nnz(Fi)>0
                    [D,V] = getfactors(Fi);
                    if length(D) <= options.sdplr.maxrank
                        lrA(k).cons = i;
                        lrA(k).start = sdploc(lmiid);
                        lrA(k).D = D;
                        lrA(k).V = V;
                        k = k+1;
                        removethese(i) = 1;
                    end
                end
            end
            removethese = find(removethese);
            if ~isempty(removethese)
                these = find((sdploc(lmiid+1)-1>= ix) & (ix>=sdploc(lmiid)) & ismember(jx,1+removethese));
                sx(these) = 0;
            end
        end
    else
        % Just check those constraints declared low-rank by user
        for lrdef = 1:length(lowrankdetails)
            for lrconstraint = 1:length(lowrankdetails{lrdef}.id)
                lmiid = lowrankdetails{lrdef}.id(lrconstraint);
                removethese = zeros(1,size(F_struc,2)-1);
                checkthese = lowrankdetails{lrdef}.variables;
                if isempty(checkthese)
                    checkthese = 1:size(F_struc,2)-1;
                end
                for i = checkthese
                    Fi = reshape(F_struc(sdploc(lmiid):sdploc(lmiid+1)-1,i+1),K.s(lmiid),K.s(lmiid));
                    if nnz(Fi)>0
                        [D,V] = getfactors(Fi);
                        if (options.sdplr.maxrank == 0) | (options.sdplr.maxrank ~= 0 & (length(D) <= options.sdplr.maxrank))
                            lrA(k).cons = i;
                            lrA(k).start = sdploc(lmiid);
                            lrA(k).D = D;
                            lrA(k).V = V;
                            k = k+1;
                            removethese(i) = 1;
                        end
                    end
                end
                removethese = find(removethese);
                if ~isempty(removethese)
                    these = find((sdploc(lmiid+1)-1>= ix) & (ix>=sdploc(lmiid)) & ismember(jx,1+removethese));
                    sx(these) = 0;
                    %F_struc(sdploc(lmiid):sdploc(lmiid+1)-1,1+removethese) = 0;
                end
            end
        end
        F_struc = sparse(ix,jx,sx,size(F_struc,1),size(F_struc,2));
    end
end

% *********************************************
% CALL SDPLR
% *********************************************
if options.showprogress;showprogress(['Calling ' interfacedata.solver.tag],options.showprogress);end
solvertime = tic;
if isempty(lrA)
   [x_s,y_s,info] = sdplr(F_struc(:,2:end),c,F_struc(:,1),K);
else
   % pars.reduce = 0;
    [x_s,y_s,info] = sdplr(F_struc(:,2:end),full(c),full(F_struc(:,1)),K,pars,lrA);
end
solvertime = toc(solvertime);

% YALMIP format
D_struc = x_s;
x = -y_s;

% No error codes currently...
problem = 0;

% Save ALL data sent to solver
if options.savesolverinput
    solverinput.A = F_struc(:,2:end);
    solverinput.c = F_struc(:,1);
    solverinput.b = c;
    solverinput.K = K;
    solverinput.pars = pars;
else
    solverinput = [];
end

% Save ALL data from the solution?
if options.savesolveroutput
    solveroutput.x = x_s;
    solveroutput.y = y_s;
    solveroutput.info = info;
else
    solveroutput = [];
end

% Standard interface
output = createOutputStructure(x(:),D_struc,[],problem,interfacedata.solver.tag,solverinput,solveroutput,solvertime);

function [D,V] = getfactors(Fi)
if nnz(Fi)>0
    [v,d] = eig(full(Fi));
    d = diag(d);
    keep = find(abs(d)>1e-6);
    V = v(:,keep);
    D = d(keep);
    % lrA(k).cons = i;
    % lrA(k).start = sdploc(j);
    % lrA(k).D = D;
    % lrA(k).V = V;
    % k = k+1;
else
    D = [];
    V = [];
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP