text:
function demprior(action);
%DEMPRIOR Demonstrate sampling from a multi-parameter Gaussian prior.
%
%	Description
%	This function plots the functions represented by a multi-layer
%	perceptron network when the weights are set to values drawn from a
%	Gaussian prior distribution. The parameters AW1, AB1 AW2 and AB2
%	control the inverse variances of the first-layer weights, the hidden
%	unit  biases, the second-layer weights and the output unit biases
%	respectively.  Their values can be adjusted on a logarithmic scale
%	using the sliders, or  by typing values into the text boxes and
%	pressing the return key.
%
%	See also
%	MLP
%

%	Copyright (c) Ian T Nabney (1996-2001)

if nargin<1,
    action='initialize';
end;

if strcmp(action,'initialize')

  aw1 = 0.01;
  ab1 = 0.1;
  aw2 = 1.0;
  ab2 = 1.0;

  % Create FIGURE
  fig=figure( ...
    'Name','Sampling from a Gaussian prior', ...
    'Position', [50 50 480 380], ...
    'NumberTitle','off', ...
    'Color', [0.8 0.8 0.8], ...
    'Visible','on');

  % The TITLE BAR frame
  uicontrol(fig,  ...
    'Style','frame', ...
    'Units','normalized', ...
    'HorizontalAlignment', 'center', ...
    'Position', [0.5 0.82 0.45 0.1], ...
    'BackgroundColor',[0.60 0.60 0.60]);

  % The TITLE BAR text
  uicontrol(fig, ...
    'Style', 'text', ...
    'Units', 'normalized', ...
    'BackgroundColor', [0.6 0.6 0.6], ...
    'Position', [0.54 0.85 0.40 0.05], ...
    'HorizontalAlignment', 'left', ...
    'String', 'Sampling from a Gaussian prior');

  % Frames to enclose sliders
  uicontrol(fig, ...
    'Style', 'frame', ...
    'Units', 'normalized', ...
    'BackgroundColor', [0.6 0.6 0.6], ...
    'Position', [0.05 0.08 0.35 0.18]);

  uicontrol(fig, ...
    'Style', 'frame', ...
    'Units', 'normalized', ...
    'BackgroundColor', [0.6 0.6 0.6], ...
    'Position', [0.05 0.3 0.35 0.18]);

  uicontrol(fig, ...
    'Style', 'frame', ...
    'Units', 'normalized', ...
    'BackgroundColor', [0.6 0.6 0.6], ...
    'Position', [0.05 0.52 0.35 0.18]);

  uicontrol(fig, ...
    'Style', 'frame', ...
    'Units', 'normalized', ...
    'BackgroundColor', [0.6 0.6 0.6], ...
    'Position', [0.05 0.74 0.35 0.18]);

  % Frame text
  uicontrol(fig, ...
    'Style', 'text', ...
    'Units', 'normalized', ...
    'HorizontalAlignment', 'left', ...
    'BackgroundColor', [0.6 0.6 0.6], ...
    'Position', [0.07 0.17 0.06 0.07], ...
    'String', 'aw1');

  % Frame text
  uicontrol(fig, ...
    'Style', 'text', ...
    'Units', 'normalized', ...
    'HorizontalAlignment', 'left', ...
    'BackgroundColor', [0.6 0.6 0.6], ...
    'Position', [0.07 0.39 0.06 0.07], ...
    'String', 'ab1');

  % Frame text
  uicontrol(fig, ...
    'Style', 'text', ...
    'Units', 'normalized', ...
    'HorizontalAlignment', 'left', ...
    'BackgroundColor', [0.6 0.6 0.6], ...
    'Position', [0.07 0.61 0.06 0.07], ...
    'String', 'aw2');

  % Frame text
  uicontrol(fig, ...
    'Style', 'text', ...
    'Units', 'normalized', ...
    'HorizontalAlignment', 'left', ...
    'BackgroundColor', [0.6 0.6 0.6], ...
    'Position', [0.07 0.83 0.06 0.07], ...
    'String', 'ab2');

  % Slider
  minval = -5; maxval = 5;
  aw1slide = uicontrol(fig, ...
    'Style', 'slider', ...
    'Units', 'normalized', ...
    'Value', log10(aw1), ...
    'BackgroundColor', [0.8 0.8 0.8], ...
    'Position', [0.07 0.1 0.31 0.05], ...
    'Min', minval, 'Max', maxval, ...
    'Callback', 'demprior update');

  % Slider
  ab1slide = uicontrol(fig, ...
    'Style', 'slider', ...
    'Units', 'normalized', ...
    'Value', log10(ab1), ...
    'BackgroundColor', [0.8 0.8 0.8], ...
    'Position', [0.07 0.32 0.31 0.05], ...
    'Min', minval, 'Max', maxval, ...
    'Callback', 'demprior update');

  % Slider
  aw2slide = uicontrol(fig, ...
    'Style', 'slider', ...
    'Units', 'normalized', ...
    'Value', log10(aw2), ...
    'BackgroundColor', [0.8 0.8 0.8], ...
    'Position', [0.07 0.54 0.31 0.05], ...
    'Min', minval, 'Max', maxval, ...
    'Callback', 'demprior update');

  % Slider
  ab2slide = uicontrol(fig, ...
    'Style', 'slider', ...
    'Units', 'normalized', ...
    'Value', log10(ab2), ...
    'BackgroundColor', [0.8 0.8 0.8], ...
    'Position', [0.07 0.76 0.31 0.05], ...
    'Min', minval, 'Max', maxval, ...
    'Callback', 'demprior update');

  % The graph box
  haxes = axes('Position', [0.5 0.28 0.45 0.45], ...
    'Units', 'normalized', ...
    'Visible', 'on');

  % Text display of hyper-parameter values

  format = '%8f';

  aw1val = uicontrol(fig, ...
    'Style', 'edit', ...
    'Units', 'normalized', ...
    'Position', [0.15 0.17 0.23 0.07], ...
    'String', sprintf(format, aw1), ...
    'Callback', 'demprior newval');

  ab1val = uicontrol(fig, ...
    'Style', 'edit', ...
    'Units', 'normalized', ...
    'Position', [0.15 0.39 0.23 0.07], ...
    'String', sprintf(format, ab1), ...
    'Callback', 'demprior newval');

  aw2val = uicontrol(fig, ...
    'Style', 'edit', ...
    'Units', 'normalized', ...
    'Position', [0.15 0.61 0.23 0.07], ...
    'String', sprintf(format, aw2), ...
    'Callback', 'demprior newval');

  ab2val = uicontrol(fig, ...
    'Style', 'edit', ...
    'Units', 'normalized', ...
    'Position', [0.15 0.83 0.23 0.07], ...
    'String', sprintf(format, ab2), ...
    'Callback', 'demprior newval');

  % The SAMPLE button
  uicontrol(fig, ...
    'Style','push', ...
    'Units','normalized', ...
    'BackgroundColor', [0.6 0.6 0.6], ...
    'Position',[0.5 0.08 0.13 0.1], ...
    'String','Sample', ...
    'Callback','demprior replot');

  % The CLOSE button
  uicontrol(fig, ...
    'Style','push', ...
    'Units','normalized', ...
    'BackgroundColor', [0.6 0.6 0.6], ...
    'Position',[0.82 0.08 0.13 0.1], ...
    'String','Close', ...
    'Callback','close(gcf)');

  % The HELP button
  uicontrol(fig, ...
    'Style','push', ...
    'Units','normalized', ...
    'BackgroundColor', [0.6 0.6 0.6], ...
    'Position',[0.66 0.08 0.13 0.1], ...
    'String','Help', ...
    'Callback','demprior help');

   % Save handles to objects

  hndlList=[fig aw1slide ab1slide aw2slide ab2slide aw1val ab1val aw2val ...
      ab2val haxes];
  set(fig, 'UserData', hndlList);

  demprior('replot')

elseif strcmp(action, 'update'),

  % Update when a slider is moved.

  hndlList   = get(gcf, 'UserData');
  aw1slide   = hndlList(2);
  ab1slide = hndlList(3);
  aw2slide  = hndlList(4);
  ab2slide = hndlList(5);
  aw1val = hndlList(6);
  ab1val = hndlList(7);
  aw2val = hndlList(8);
  ab2val = hndlList(9);
  haxes = hndlList(10);

  aw1 = 10^get(aw1slide, 'Value');
  ab1 = 10^get(ab1slide, 'Value');
  aw2 = 10^get(aw2slide, 'Value');
  ab2 = 10^get(ab2slide, 'Value');

  format = '%8f';
  set(aw1val, 'String', sprintf(format, aw1));
  set(ab1val, 'String', sprintf(format, ab1));
  set(aw2val, 'String', sprintf(format, aw2));
  set(ab2val, 'String', sprintf(format, ab2));

  demprior('replot');

elseif strcmp(action, 'newval'),

  % Update when text is changed.

  hndlList   = get(gcf, 'UserData');
  aw1slide   = hndlList(2);
  ab1slide = hndlList(3);
  aw2slide  = hndlList(4);
  ab2slide = hndlList(5);
  aw1val = hndlList(6);
  ab1val = hndlList(7);
  aw2val = hndlList(8);
  ab2val = hndlList(9);
  haxes = hndlList(10);

  aw1 = sscanf(get(aw1val, 'String'), '%f');
  ab1 = sscanf(get(ab1val, 'String'), '%f');
  aw2 = sscanf(get(aw2val, 'String'), '%f');
  ab2 = sscanf(get(ab2val, 'String'), '%f');

  set(aw1slide, 'Value', log10(aw1));
  set(ab1slide, 'Value', log10(ab1));
  set(aw2slide, 'Value', log10(aw2));
  set(ab2slide, 'Value', log10(ab2));

  demprior('replot');

elseif strcmp(action, 'replot'),

  % Re-sample from the prior and plot graphs.

  oldFigNumber=watchon;

  hndlList   = get(gcf, 'UserData');
  aw1slide   = hndlList(2);
  ab1slide = hndlList(3);
  aw2slide  = hndlList(4);
  ab2slide = hndlList(5);
  haxes = hndlList(10);

  aw1 = 10^get(aw1slide, 'Value');
  ab1 = 10^get(ab1slide, 'Value');
  aw2 = 10^get(aw2slide, 'Value');
  ab2 = 10^get(ab2slide, 'Value');

  axes(haxes);
  cla
  set(gca, ...
    'Box', 'on', ...
    'Color', [0 0 0], ...
    'XColor', [0 0 0], ...
    'YColor', [0 0 0], ...
    'FontSize', 14);
  axis([-1 1 -10 10]);
  set(gca,'DefaultLineLineWidth', 2);

  nhidden = 12;
  prior = mlpprior(1, nhidden, 1, aw1, ab1, aw2, ab2);
  xvals = -1:0.005:1;
  nsample = 10;    % Number of samples from prior.
  hold on
  plot([-1 0; 1 0], [0 -10; 0 10], 'b--');
  net = mlp(1, nhidden, 1, 'linear', prior);
  for i = 1:nsample
    net = mlpinit(net, prior);
    yvals = mlpfwd(net, xvals');
    plot(xvals', yvals, 'y');
  end

  watchoff(oldFigNumber);

elseif strcmp(action, 'help'),

  % Provide help to user.

  oldFigNumber=watchon;

  helpfig = figure('Position', [100 100 480 400], ...
    'Name', 'Help', ...
    'NumberTitle', 'off', ...
    'Color', [0.8 0.8 0.8], ...
    'Visible','on');

    % The HELP TITLE BAR frame
  uicontrol(helpfig,  ...
    'Style','frame', ...
    'Units','normalized', ...
    'HorizontalAlignment', 'center', ...
    'Position', [0.05 0.82 0.9 0.1], ...
    'BackgroundColor',[0.60 0.60 0.60]);

  % The HELP TITLE BAR text
  uicontrol(helpfig, ...
    'Style', 'text', ...
    'Units', 'normalized', ...
    'BackgroundColor', [0.6 0.6 0.6], ...
    'Position', [0.26 0.85 0.6 0.05], ...
    'HorizontalAlignment', 'left', ...
    'String', 'Help: Sampling from a Gaussian Prior');

  helpstr1 = strcat( ...
    'This demonstration shows the effects of sampling from a Gaussian', ...
     ' prior over weights for a two-layer feed-forward network. The', ...
     ' parameters aw1, ab1, aw2 and ab2 control the inverse variances of', ...
     ' the first-layer weights, the hidden unit biases, the second-layer', ...
     ' weights and the output unit biases respectively. Their values can', ...
     ' be adjusted on a logarithmic scale using the sliders, or by', ...
     ' typing values into the text boxes and pressing the return key.', ...
     '  After setting these values, press the ''Sample'' button to see a', ...
     ' new sample from the prior. ');
   helpstr2 = strcat( ...
     'Observe how aw1 controls the horizontal length-scale of the', ...
     ' variation in the functions, ab1 controls the input range over', ...
     ' such variations occur, aw2 sets the vertical scale of the output', ...
     ' and ab2 sets the vertical off-set of the output. The network has', ...
     ' 12 hidden units. ');
   hstr(1) = {helpstr1};
   hstr(2) = {''};
   hstr(3) = {helpstr2};

  % The HELP text
  helpui = uicontrol(helpfig, ...
    'Style', 'edit', ...
    'Units', 'normalized', ...
    'ForegroundColor', [0 0 0], ...
    'HorizontalAlignment', 'left', ...
    'BackgroundColor', [1 1 1], ...
    'Min', 0, ...
    'Max', 2, ...
    'Position', [0.05 0.2 0.9 0.8]);

   [hstrw , newpos] = textwrap(helpui, hstr, 70);
   set(helpui, 'String', hstrw, 'Position', [0.05, 0.2, 0.9, newpos(4)]);

  % The CLOSE button
  uicontrol(helpfig, ...
    'Style','push', ...
    'Units','normalized', ...
    'BackgroundColor', [0.6 0.6 0.6], ...
    'Position',[0.4 0.05 0.2 0.1], ...
    'String','Close', ...
    'Callback','close(gcf)');

   watchoff(oldFigNumber);

end;


text:
% initFigures
% Position 6 figures on the edges of the screen.
% [xmin ymin w h] where (0,0) = bottom left
% Numbers assume screen resolution is 1024 x 1280

global FIGNUM NUMFIGS
FIGNUM = 1; NUMFIGS = 6;

screenMain = true; % set to false if initializing figures for second screen
%screenMain = false;

if screenMain
  xoff = 0;
else
  %xoff = 1280;
  xoff = -1280;
end

% 2 x 3 design
w = 400; h = 300;
xs = [10 450 875] + xoff;
ys = [650 40];

if 0
% 3x3 design
w = 350; h = 250;
xs = [10 380 750]+xoff;
ys = [700 350 10];
end

Nfigs = length(xs)*length(ys);
if screenMain
  fig = 1;
else
  fig = Nfigs + 1;
end

for yi=1:length(ys)
  for xi=1:length(xs)
    figure(fig);
    set(gcf, 'position', [xs(xi) ys(yi) w h]);
    fig = fig + 1;
  end
end

% To plot something on the next available figure (with wrap around), use
% sfigure(FIGNUM); clf; FIGNUM = wrap(FIGNUM+1, NUMFIGS);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
classdef (Abstract) mm_shared_pfcpf_acps < mp.mm_shared_pfcpf_acp

%   MATPOWER
%   Copyright (c) 2021-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end

    methods
        function ad = build_aux_data(obj, nm, dm, mpopt)
            %% call parent
            ad = build_aux_data@mp.mm_shared_pfcpf_acp(obj, nm, dm, mpopt);

            switch mpopt.pf.alg
                case 'GS'
                    ad.Y = nm.C * nm.get_params([], 'Y') * nm.C';
                case 'ZG'
                    pvq = [ad.pv; ad.pq];
                    Y = nm.C * nm.get_params([], 'Y') * nm.C';
                    Y21 = Y(pvq, ad.ref);
                    Y22 = Y(pvq, pvq);
                    [L, U, p, q] = lu(Y22, 'vector');
                    [junk, iq] = sort(q);
                    [ad.Y, ad.Y21, ad.L, ad.U, ad.p, ad.iq] = ...
                        deal(Y, Y21, L, U, p, iq);
            end
        end

        function obj = add_system_vars_pf(obj, nm, dm, mpopt)
            %% get model variables
            vvars = nm.model_vvars();

            %% voltage angles
            obj.add_system_varset_pf(nm, vvars{1}, 'pv');
            obj.add_system_varset_pf(nm, vvars{1}, 'pq');

            %% voltage magnitudes
            obj.add_system_varset_pf(nm, vvars{2}, 'pq');
        end

        function [f, J] = node_balance_equations(obj, x, nm, fdpf)
            %% index vector
            ad = obj.aux_data;
            pvq = [ad.pv; ad.pq];

            %% update network model state ([v_; z_]) from math model state (x)
            [v_, z_] = obj.convert_x_m2n(x, nm, 1);

            %% incidence matrix
            C = nm.C;

            %% Jacobian
            if nargout > 1
                %% get port power injections with derivatives
                var_names = cellfun(@(x)x{1}, ad.var_map, 'UniformOutput', false);
                dz = any(strcmp(var_names, 'zr')) || ...
                     any(strcmp(var_names, 'zi'));
                if dz
                    [S, dS.va, dS.vm, dS.zr, dS.zi] = nm.port_inj_power([v_; z_], 1);
                else
                    [S, dS.va, dS.vm] = nm.port_inj_power([v_; z_], 1);
                end
                dS.va = C * dS.va;
                dS.vm = C * dS.vm;
                if dz
                    dS.zr = C * dS.zr;
                    dS.zi = C * dS.zi;
                end
                JJ = cell(2, length(ad.var_map));

                for k = 1:length(ad.var_map)
                    m = ad.var_map{k};
                    name = m{1};
                    if ~isempty(m{2})       %% i1:iN
                        i1 = m{2};
                        iN = m{3};
                        JJ{1, k} = real(dS.(name)(pvq,   i1:iN));
                        JJ{2, k} = imag(dS.(name)(ad.pq, i1:iN));
                    elseif isempty(m{4})    %% :
                        JJ{1, k} = real(dS.(name)(pvq,   :));
                        JJ{2, k} = imag(dS.(name)(ad.pq, :));
                    else                    %% idx
                        idx = m{4};
                        JJ{1, k} = real(dS.(name)(pvq,   idx));
                        JJ{2, k} = imag(dS.(name)(ad.pq, idx));
                    end
                end
                J = vertcat( horzcat(JJ{1, :}), ...
                             horzcat(JJ{2, :})  );
            else
                %% get port power injections (w/o derivatives)
                S = nm.port_inj_power([v_; z_], 1);
            end

            %% nodal power balance
            if nargin > 3 && fdpf
                SS = C * S ./ abs(v_);  %% for fast-decoupled formulation
            else
                SS = C * S;
            end
            f = [real(SS(pvq)); imag(SS(ad.pq))];
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [ y, m, d ] = nyt_to_ymd ( volume, issue )

%*****************************************************************************80
%
%% NYT_TO_YMD converts an NYT date to a YMD date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer VOLUME, ISSUE, the New York Times
%    volume and issue.
%
%    Output, integer Y, M, D, the year, month and day.
%
  jed = nyt_to_jed ( volume, issue );

  [ y, m, d, f ] = jed_to_ymdf_common ( jed );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [A,b,C,P] = update_spatial_components_old(Y,C,f,A_,P,options)

% update spatial footprints and background through Basis Pursuit Denoising
% for each pixel i solve the problem
%   [A(i,:),b(i)] = argmin sum(A(i,:))
%       subject to || Y(i,:) - A(i,:)*C + b(i)*f || <= sn(i)*sqrt(T);
% for each pixel the search is limited to a few spatial components

% INPUTS:
% Y:    raw data
% C:    temporal components
% f:    temporal background
% A_:   current estimate of spatial footprints (used for determining search locations only)
% P:    dataset parameters (used for noise values and interpolated entries)

% options    parameter struct (for noise values and other parameters)

% OUTPUTS:
% A:    new estimate of spatial footprints
% b:    new estimate of spatial background
% C:    temporal components (updated only when spatial components are completely removed)

% Written by:
% Eftychios A. Pnevmatikakis, Simons Foundation, 2015

warning('off', 'MATLAB:maxNumCompThreads:Deprecated');
memmaped = isobject(Y);
if memmaped
    sizY = size(Y,'Y');
    d = prod(sizY(1:end-1));
    T = sizY(end);
else
    [d,T] = size(Y);
end
if nargin < 6 || isempty(options); options = []; end
if ~isfield(options,'d1') || isempty(options.d1); d1 = input('What is the total number of rows? \n'); options.d1 = d1; else d1 = options.d1; end          % # of rows
if ~isfield(options,'d2') || isempty(options.d2); d2 = input('What is the total number of columns? \n'); options.d2 = d2; else d2 = options.d2; end          % # of columns
if ~isfield(options,'d3') || isempty(options.d3); d3 = input('What is the total number of z-planes? \n'); options.d3 = d3; else d3 = options.d3; end          % # of columns
if ~isfield(options,'show_sum'); show_sum = 0; else show_sum = options.show_sum; end            % do some plotting while calculating footprints
if ~isfield(options,'interp'); Y_interp = sparse(d,T); else Y_interp = options.interp; end      % identify missing data
if ~isfield(options,'spatial_parallel'); spatial_parallel = ~isempty(which('parpool')); else spatial_parallel = options.spatial_parallel; end % use parallel toolbox if present
if ~isfield(options,'search_method'); method = []; else method = options.search_method; end     % search method for determining footprint of spatial components
if ~isfield(options,'tsub') || isempty(options.tsub); tsub = 1; else tsub = options.tsub; end  % downsample temporally to estimate A and b

if nargin < 2 || (isempty(A_) && isempty(C))  % at least either spatial or temporal components should be provided
    error('Not enough input arguments')
else
    if ~isempty(C); K = size(C,1); elseif islogical(A_); K = size(A_,2); else K = size(A_2,2) - options.nb; end
end

if nargin < 5 || isempty(P); P = preprocess_data(Y,1); end  % etsimate noise values if not present
if nargin < 4 || isempty(A_);
    IND = ones(d,size(C,1));
else
    if islogical(A_)     % check if search locations have been provided, otherwise estimate them
        IND = A_;
        if isempty(C)
            INDav = double(IND)/diag(sum(double(IND)));
            px = (sum(IND,2)>0);
            f = mean(Y(~px,:));
            b = max(Y*f',0)/norm(f)^2;
            C = max(INDav'*Y - (INDav'*b)*f,0);
        end
    else
        IND = determine_search_location(A_(:,1:K),method,options);
    end
end

K = size(C,1);
if strcmpi(options.spatial_method,'constrained'); A_ = A_(:,1:K); end

Cf = [C;f];
if size(Cf,1) > size(A_,2) && strcmpi(options.spatial_method,'regularized');
    error('When using options.spatial_method = regularized pass [A,b] as an input and not just A');
end

if tsub ~= 1 %&& strcmpi(options.spatial_method,'constrained');
    P.sn_ds = zeros(d,1);
    Ts = floor(T/tsub);
    Cf = squeeze(mean(reshape(Cf(:,1:Ts*tsub),[],tsub,Ts),2));
    f = Cf(end-size(f,1)+1:end,:);
    T = Ts;
    if ~memmaped
        Y = squeeze(mean(reshape(Y(:,1:Ts*tsub),[],tsub,Ts),2));
        if ~isfield(P,'sn_ds');
            [P_ds,Y] = preprocess_data(Y);
            P.sn_ds = P_ds.sn;
        end
        options.sn = P.sn_ds;
    else
        Y_ds = zeros(d,Ts);
        step_size = 2e4;
        for i = 1:step_size:d
            Ytemp = double(Y.Yr(i:min(i+step_size-1,d),:));
            Y_ds(i:min(i+step_size-1,d),:) = squeeze(mean(reshape(Ytemp(:,1:Ts*tsub),[],tsub,Ts),2));
        end
        Y = Y_ds;
        options.sn = get_noise_fft(Y_ds);
    end
else
    options.sn = P.sn;
end

if strcmpi(options.spatial_method,'constrained');
    if spatial_parallel         % solve BPDN problem for each pixel
        Nthr = max(20*maxNumCompThreads,round(d*T/2^24));
        Nthr = min(Nthr,round(d/1e3));
        siz_row = [floor(d/Nthr)*ones(Nthr-mod(d,Nthr),1);(floor(d/Nthr)+1)*ones(mod(d,Nthr),1)];
        indeces = [0;cumsum(siz_row)];
        Yf = cell(Nthr,1);
        A = spalloc(d,size(Cf,1),nnz(IND)+size(f,1)*d);
        for nthr = 1:Nthr
            if memmaped
                Ytemp = double(Y.Yr(indeces(nthr)+1:indeces(nthr+1),:));
                if tsub ~= 1;
                    Ytemp = squeeze(mean(reshape(Ytemp(:,1:Ts*tsub),[],tsub,Ts),2));
                    if ~isfield(P,'sn_ds')
                        [P_ds,Ytemp] = preprocess_data(Ytemp);
                        P.sn_ds(indeces(nthr)+1:indeces(nthr+1)) = P_ds.sn;
                    end
                    sn_temp = P.sn_ds(indeces(nthr)+1:indeces(nthr+1));
                else
                    sn_temp = P.sn(indeces(nthr)+1:indeces(nthr+1));
                end
            else
                Ytemp = Y(indeces(nthr)+1:indeces(nthr+1),:);
                sn_temp = options.sn(indeces(nthr)+1:indeces(nthr+1));
            end
            IND_temp = IND(indeces(nthr)+1:indeces(nthr+1),:);
            Atemp = spalloc(siz_row(nthr),size(Cf,1),nnz(IND_temp));
            Yf{nthr} = Ytemp*f';

            parfor px = 1:siz_row(nthr)
                fn = ~isnan(Ytemp(px,:));       % identify missing data
                ind = find(IND_temp(px,:));
                if ~isempty(ind);
                    ind2 = [ind,K+(1:size(f,1))];
                    [~, ~, a, ~] = lars_regression_noise(Ytemp(px,fn)', Cf(ind2,fn)', 1, sn_temp(px)^2*T);
                    a_sparse = sparse(1,ind2,double(a'));
                    Atemp(px,:) = a_sparse';
                end
            end
            if mod(nthr,50) == 0
                fprintf('%2.1f%% of pixels completed \n', indeces(nthr+1)*100/d);
            end
            %Acell{nthr} = Atemp;
            A(indeces(nthr)+1:indeces(nthr+1),:) = Atemp;
        end
        %A = cell2mat(Acell);
        Yf = cell2mat(Yf);
    else
        A = [zeros(d,K),zeros(d,size(f,1))];
        sA = zeros(d1,d2,d3);
        Yf = Y*f';
        for px = 1:d   % estimate spatial components
            fn = ~isnan(Y(px,:));       % identify missing data
            ind = find(IND(px,:));
            if ~isempty(ind);
                ind2 = [ind,K+(1:size(f,1))];
                [~, ~, a, ~] = lars_regression_noise(Y(px,fn)', Cf(ind2,fn)', 1, options.sn(px)^2*T);
                A(px,ind2) = a';
                sA(px) = sum(a);
            end
            if show_sum
                if mod(px,d1) == 0;
                   figure(20); imagesc(max(sA,[],3)); axis square;
                   title(sprintf('Sum of spatial components (%i out of %i columns done)',round(px/d1),d2)); drawnow;
                end
            end
        end
    end
elseif strcmpi(options.spatial_method,'regularized')

    A = update_spatial_lasso(Y, A_, Cf, IND, options.sn, [], [], options);
    K = size(A,2)-options.nb;
    b = full(A(:,K+1:end));
    A = A(:,1:K);
    C = C(1:K,:);
end

A(isnan(A))=0;
A = sparse(A);
A = threshold_components(A,options);  % post-processing of components

fprintf('Updated spatial components \n');

ff = find(sum(A(:,1:K))==0);           % remove empty components
if ~isempty(ff)
    K = K - length(ff);
    A(:,ff) = [];
    C(ff,:) = [];
    Cf(ff,:) = [];
end

if strcmpi(options.spatial_method,'constrained');
    b = double(max((double(Yf) - A(:,1:K)*double(Cf(1:K,:)*f'))/(f*f'),0));
    A = A(:,1:K);
end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function test_bug1775

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_sourceparcellate ft_checkdata ft_datatype_source ft_datatype_volume ft_datatype_parcellation ft_datatype_segmentation

%% create a set of sensors
[pnt, tri] = mesh_sphere(162);
pnt = pnt .* 10; % convert to cm
sel = find(pnt(:,3)>0);

grad.pnt = pnt(sel,:) .* 1.2;
grad.ori = pnt(sel,:);
grad.tra = eye(length(sel));
for i=1:length(sel)
  grad.ori(i,:) = grad.ori(i,:) ./ norm(grad.ori(i,:));
  grad.label{i} = sprintf('magnetometer%d', i);
end
grad.unit = 'cm';
grad.type = 'meg';

grad = ft_datatype_sens(grad);

%% create a volume conductor

vol = [];
vol.r     = 10;
vol.o     = [0 0 0];
vol.unit  = 'cm';

vol = ft_datatype_headmodel(vol);

%% create some precomputed leadfields

cfg = [];
cfg.grad            = grad;
cfg.headmodel       = vol;
cfg.resolution      = 2; % cm
cfg.channel         = 'all';
grid = ft_prepare_leadfield(cfg);

%% create an anatomical parcellation
parcellation = [];
parcellation.pos          = grid.pos;
parcellation.unit         = grid.unit;
parcellation.tissue       = zeros(size(grid.pos,1),1);
parcellation.tissuelabel  = {};
height = [3 4 5 6 7 8 9];
for i=1:length(height)
  sel = parcellation.pos(:,3)==height(i);
  parcellation.tissue(sel) = i;
  parcellation.tissuelabel{i} = sprintf('%d%s', height(i), parcellation.unit);
end
parcellation.cfg = 'manual'; % to check whether the provenance is correct

%% create simulated data
cfg = [];
cfg.grad = grad;
cfg.headmodel = vol;
cfg.dip.pos = [0 0 4];
data = ft_dipolesimulation(cfg);

cfg = [];
cfg.covariance = 'yes';
timelock = ft_timelockanalysis(cfg, data);

cfg = [];
cfg.method  = 'mtmfft';
cfg.taper   = 'hanning';
freq1 = ft_freqanalysis(cfg, data);

cfg = [];
cfg.method  = 'wavelet';
cfg.toi     = data.time{1};
freq2 = ft_freqanalysis(cfg, data);

cfg = [];
cfg.grad        = grad;
cfg.headmodel   = vol;
cfg.sourcemodel = grid;
cfg.method      = 'lcmv';
source1 = ft_sourceanalysis(cfg, timelock);

cfg = [];
cfg.grad        = grad;
cfg.headmodel   = vol;
cfg.sourcemodel = grid;
cfg.method      = 'mne';
cfg.mne.lambda  = 0;
source2 = ft_sourceanalysis(cfg, timelock);

%% make some parcellations
cfg = [];
gridp    = ft_sourceparcellate(cfg, grid, parcellation);
source1p = ft_sourceparcellate(cfg, source1, parcellation);
source2p = ft_sourceparcellate(cfg, source2, parcellation);

%% construct a more complex source structure
% note that this increases memory requirements
source3 = [];
source3.pos       = source2.pos;
source3.freq      = 1:5;
source3.coh       = randn(size(source2.pos,1), size(source2.pos,1), 5);
source3.cohdimord = 'pos_pos_freq';

cfg = [];
source3p = ft_sourceparcellate(cfg, source3, parcellation);

%%
% this increases memory requirements even more
source4 = [];
source4.pos       = source2.pos;
source4.freq      = 1:5;
source4.time      = 1:3;
source4.coh       = randn(size(source2.pos,1), size(source2.pos,1), 5, 3);
source4.cohdimord = 'pos_pos_freq_time';

cfg = [];
source4p = ft_sourceparcellate(cfg, source4, parcellation);

%%
source5 = [];
source5.pos       = source2.pos;
source5.inside    = source2.inside;
source5.freq      = 1:2;
source5.coh       = cell(size(source2.pos,1), size(source2.pos,1));
for i=find(source2.inside)'
  for j=find(source2.inside)'
    source5.coh{i,j} = randn(3, 2);
  end
end
source5.cohdimord = '{pos_pos}_ori_freq';

cfg = [];
cfg.method = 'mean';
source5p = ft_sourceparcellate(cfg, source5, parcellation);
cfg.method = 'min';
source5p = ft_sourceparcellate(cfg, source5, parcellation);
cfg.method = 'max';
source5p = ft_sourceparcellate(cfg, source5, parcellation);

%%
source6 = [];
source6.pos       = source2.pos;
source6.inside    = source2.inside;
source6.time      = 1:20;
source6.mom       = randn(size(source2.pos,1),20);
source6.momdimord = 'pos_time';

cfg = [];
cfg.method = 'mean';
source6p = ft_sourceparcellate(cfg, source6, parcellation);
cfg.method = 'min';
source6p = ft_sourceparcellate(cfg, source6, parcellation);
cfg.method = 'max';
source6p = ft_sourceparcellate(cfg, source6, parcellation);
cfg.method = 'eig';
source6p = ft_sourceparcellate(cfg, source6, parcellation);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function points2d
%POINTS2D  Description of functions operating on points
%
%   A point is defined by its two cartesian coordinate, put into a row
%   vector of 2 elements:
%   P = [x y];
%
%   Several points are stores in a matrix with two columns, one for the
%   x-coordinate, one for the y-coordinate.
%   PTS = [x1 y1 ; x2 y2 ; x3 y3];
%
%   Example
%   P = [5 6];
%
%   See also:
%   centroid, midPoint, boundingBox, polarPoint
%   distancePoints, minDistancePoints, nndist, circumCenter
%   isCounterClockwise, angle2Points, angle3Points, angleSort
%   transformPoint, clipPoints, drawPoint
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2008-10-13,    using Matlab 7.4.0.287 (R2007a)
% Copyright 2008 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

help('points2d');

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function r = exp(a)
%EXP          Taylor exponential  exp(a)
%

% written  05/21/09     S.M. Rump
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                   % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  K1 = getappdata(0,'INTLAB_TAYLOR_ORDER') + 1;

  r = a;
  N = size(a.t,2);
  r.t(1,:) = exp(a.t(1,:));
  for j=2:K1
    r.t(j,:) = sum( repmat((1:j-1)',1,N).*r.t(j-1:-1:1,:).*a.t(2:j,:) , 1 ) ./ (j-1);
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function varargout = imlook3d(varargin)
% IMLOOK3D M-file for imlook3d.fig
%      IMLOOK3D, by itself, creates a new IMLOOK3D or raises the existing
%      singleton*.
%
%      H = IMLOOK3D returns the handle to a new IMLOOK3D or the handle to
%      the existing singleton*.
%
%      IMLOOK3D('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in IMLOOK3D.M with the given input arguments.
%
%      IMLOOK3D('Property','Value',...) creates a new IMLOOK3D or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before imlook3d_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to imlook3d_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% EXAMPLE:
%  Let us first load the mri image in Matlab using
%  >load mri;
%  and then run imlook3d
%  >imlook3d(squeeze(D)); % squeeze removes singleton dimensions
%
% Descriptions of Menus
%
% ImFile:
%
% -READ FROM WORKSPACE: Reads or displays images existing in the
%  workspace(base).
%
% -READ A STACK (DISK): Reads (displays)slices of images into a 3D array creating
%  a 3D image. It reads them from the disk. A file browser pups up and
%  allows the user to pick up the first file in the series. It reads all
%  the files in the directory that has similar names except a suffix
%  changing for each file.
%
%  For example: slice01.tif, slice01.tif,...,slice0N.tif
%
% -READ DICOM SERIES:
%  Reads/displays series of DICOM images into a 3D array. A file browser
%  pups up and allows users to select the first file in the series.
%  It reads all the files in the directory that has similar names
%  except a suffix changing for each file.

% -SAVE SLICE AS TIFF:
%  Saves the current slice as a tiff image.

% -SAVE INTO WORKSPACE: This saves displayed 3D image into workspace as
%  a 3D variable. This may be useful to save the 3D image that is
%  read by "read a stack (disk)".
%
% ImTools:
%
% -PIXVAL: Allows looking at pixel intensities.
% -HISTOGRAM: Displays the histogram of the current image (slice).
% -PROFILE: Allows plotting intensity profiles across the image in an
%  interactive manner.
% -ADJUST INTENSITY RANGE: Allows adjusting the intensity range for
%  the current slice.
% -GLOBAL INTENSITY WINDOWING:
%  If not checked,not global, the intensity limits are set to [min(currentslice)
%  max(currentslice)]
%
% VIEW MENU:
%
%  This pup-up menu allows users to change the view of the image
%  from transverse to sagittal or coronal views. It assumes that
%  original image is in transverse mode.
%

% Copyright(c) Omer Demirkaya
% Last Modified by GUIDE v2.5 22-May-2007 21:27:00
% Begin initialization code - DO NOT EDIT

% This entire code or part of it could be modified and used if the due
% credit to the original author is given.
%------------------------------------------------------------------

gui_Singleton = 0;
gui_State = struct('gui_Name',       mfilename, ...
    'gui_Singleton',  gui_Singleton, ...
    'gui_OpeningFcn', @imlook3d_OpeningFcn, ...
    'gui_OutputFcn',  @imlook3d_OutputFcn, ...
    'gui_LayoutFcn',  [] , ...
    'gui_Callback',   []);
if nargin & isstr(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end
if isempty(varargin)
    disp('There was no image to display');
    disp('USAGE: >>imlook3d(img) or use ImFile menu to load an image');
    varargin{1,1} = zeros(32,32);
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before imlook3d is made visible.
function imlook3d_OpeningFcn(hObject, eventdata, handles, varargin)
% Choose default command line output for imlook3d
if isempty(varargin{:})
    return;
end
handles.output = hObject;
Std.Interruptible = 'off';
Std.BusyAction = 'queue';
Ax       = Std;
Ax.Units = 'Pixels';
Ax.ydir  = 'reverse';
Ax.XLim  = [.5 128.5];
Ax.YLim  = [.5 128.5];
Ax.XTick = [];
Ax.YTick = [];
Ax.CLimMode = 'auto';

set(handles.axes1,Ax);

Img = Std;Img.CData = [];Img.Xdata = [];Img.Ydata = [];
Img.CDataMapping = 'Scaled';
Img.Erasemode = 'none';
handles.image = Img;
inpargs = varargin{:};
[r,c,z] = size(inpargs);
handles.image.CData     = inpargs;
handles.imSize          = [r,c,z];
handles.viewtype        = 1;
handles.imageModality   = [];
handles.hText           = [];
handles.global_windowing= 0;
handles.SliceThickness  = [];
handles.imageInfo       = [];
handles.currentSliceNumber = 1;
handle.image.SliceSensitivityFactor =[];

handles.oimage = handles.image;
cimg = handles.image.CData(:,:,1);
% Create image object and set the properties
handles.ImgObject = image(Img,'Parent',handles.axes1);
set(handles.SliceNumEdit,'String',1);
set(handles.figure1,'Colormap',gray,'Units','Pixels');
% Set the properties of the axes
cclim =[min(cimg(:)) max(cimg(:))];
if (cclim(2)>0 & cclim(2)-cclim(1) ~= 0)
    set(handles.axes1,'CLim',cclim,...
        'PlotBoxAspectRatio',[c r 1]);
end
set(handles.axes1,'XLim',[1 c],'Ylim',[1 r],...
    'PlotBoxAspectRatio',[c r 1]);
set(handles.ImgObject,'Xdata',[1 c],'Ydata',[1 r],'CData',cimg);
drawnow
% Set the slider's value and step size
if z > 1
    turnDisplayButtons(handles,'on')
    set(handles.SliceNumSlider,'Min',1,'Max',z,'Value',...
        1,'SliderStep',[1.0/double(z-1) 1.0/double(z-1)]);
end
% Define the line and associate it with the context menu
set(handles.ImgObject,'UIContextMenu',handles.ImMaps);
% Update handles structure
%--------------------------------------------------------------------------
%% Set appdata for the inntensity adjustment GUI
setappdata(0  ,  'hMainGui'    , gcf);
hMainGui = getappdata(0, 'hMainGui');
maxINT = max(max(max(handles.image.CData)));
minINT = min(min(min(handles.image.CData)));
intensityParam = struct('lowerLimit',cclim(1),'upperLimit',cclim(2),...
    'maxIntensity', maxINT, 'minIntensity',minINT);
handles.intensityParam = intensityParam;
setappdata(hMainGui,  'intensityParam'  ,intensityParam);
%% Read and display logo image
% ----------------------------------------------------------------
fdir = which('imlook3d');
txt = [fdir(1:end-10) '\imlook3d.png'];
imlogo = imread(txt);
Sz= size(imlogo);
flogo = figure('Position',[10000 10000 Sz(2) + 4 Sz(1) + 4],'name','Imlook3D',...
    'numbertitle','off','menubar','none');
movegui(flogo,'center');
set(flogo,'Units', 'pixels');
image(imlogo(:,:,1:3))
set(gca, ...
    'Visible', 'off', ...
    'Units', 'pixels', ...
    'Position', [2 2 Sz(2) Sz(1)]);
text(220,16,'3D Image Display Tool','Fontsize',12);
text(130,355,'Version 1.2: Copyright 2006, Omer Demirkaya','Fontsize',8);
pause(3);delete(flogo);clear imlogo
%% ------------------------------------
%--------------------------------------------------------------------------
guidata(hObject, handles);
%--------------------------------------------------------------------------
% --- Outputs from this function are returned to the command line.
function varargout = imlook3d_OutputFcn(hObject, eventdata, handles)
 varargout{1} = handles.output;
% --- Executes during object creation, after setting all properties.
function SliceNumSlider_CreateFcn(hObject, eventdata, handles)
usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor',[.9 .9 .9]);
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on slider movement.
%------------------------------------------------------------
function SliceNumSlider_Callback(hObject, eventdata, handles)
NewVal= round(get(hObject,'Value'));
cimg = handles.image.CData(:,:,NewVal);
if ~(handles.global_windowing)
    cclim = [min(cimg(:)) max(cimg(:))];
    if (cclim(2)>0 & cclim(2)-cclim(1) ~= 0)
        set(handles.axes1,'CLim',[min(cimg(:)) max(cimg(:))])
    end
else
    set(handles.axes1,'CLim',[handles.intensityParam.lowerLimit  handles.intensityParam.upperLimit])
end
set(handles.ImgObject,'Cdata',cimg);
set(handles.SliceNumEdit,'String',num2str(NewVal));
drawnow
guidata(hObject, handles);
% --- Executes during object creation, after setting all properties.
function SliceNumEdit_CreateFcn(hObject, eventdata, handles)
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
%--------------------------------------------------------------------------
function SliceNumEdit_Callback(hObject, eventdata, handles)
z= size(handles.image.CData,3);
strg = get(hObject,'String');
if str2num(strg)>=1 & str2num(strg)<=z
    set(handles.SliceNumSlider,'Value',str2num(strg));
    cimg = handles.image.CData(:,:,str2num(strg));
    set(handles.ImgObject,'Cdata',cimg);
    handles.currentSliceNumber = str2num(strg);
    drawnow
end
if ~(handles.global_windowing)
    cclim = [min(cimg(:)) max(cimg(:))];
    if (cclim(2)>0 & cclim(2)-cclim(1) ~= 0)
        set(handles.axes1,'CLim',[min(cimg(:)) max(cimg(:))])
    end
else
    set(handles.axes1,'CLim',[handles.intensityParam.lowerLimit  handles.intensityParam.upperLimit])

end
guidata(hObject, handles);
% --------------------------------------------------------------------
% --------------------------------------------------------------------
function changeColorMaps_Callback(hObject, eventdata, handles)
% --------------------------------------------------------------------
function petColormaps_Callback(hObject, eventdata, handles)
% --------------------------------------------------------------------
function CTcolormaps_Callback(hObject, eventdata, handles)
% --------------------------------------------------------------------
% --------------------------------------------------------------------
function Hot_Callback(hObject, eventdata, handles)
htable = feval('hot');
set(handles.figure1,'Colormap',htable);
if handles.imSize(3)> 1
    set(handles.SliceNumEdit,'BackgroundColor',[.1 .1 .1],'ForegroundColor','r');
end
guidata(hObject, handles);
% --------------------------------------------------------------------
function Gray_Callback(hObject, eventdata, handles)
htable = feval('gray');
set(handles.figure1,'Colormap',htable);
if handles.imSize(3)> 1
    set(handles.SliceNumEdit,'BackgroundColor','k','ForegroundColor','w');
end
guidata(hObject, handles);
% --------------------------------------------------------------------
function InvGray_Callback(hObject, eventdata, handles)
htable = 1-feval('gray');
set(handles.figure1,'Colormap',htable);
if handles.imSize(3)> 1
    set(handles.SliceNumEdit,'BackgroundColor','w','ForegroundColor','k');
end
guidata(hObject, handles);
% --------------------------------------------------------------------
function Jet_Callback(hObject, eventdata, handles)
htable = feval('jet');
set(handles.figure1,'Colormap',htable);
if handles.imSize(3)> 1
    set(handles.SliceNumEdit,'BackgroundColor','b','ForegroundColor','y');
end
guidata(hObject, handles);
% --------------------------------------------------------------------
function boneWindow_Callback(hObject, eventdata, handles)
handles.intensityParam.lowerLimit = -250;
handles.intensityParam.upperLimit = 1250;
set(handles.axes1,'CLim',[-300 1250])
guidata(hObject,handles);
% --------------------------------------------------------------------
function lungWindow_Callback(hObject, eventdata, handles)
handles.intensityParam.lowerLimit = -2500;
handles.intensityParam.upperLimit = 1100;
set(handles.axes1,'CLim',[-2500 1100])
guidata(hObject,handles);
% --------------------------------------------------------------------
function mediasitnalWindow_Callback(hObject, eventdata, handles)
handles.intensityParam.lowerLimit = -450;
handles.intensityParam.upperLimit = 350;
set(handles.axes1,'CLim',[-450 350])
guidata(hObject,handles);
% --------------------------------------------------------------------
%---------------------------------------------------------
% Callbacks for imtoolsmenu
%----------------------------------------------------------
function imToolsMenu_Callback(hObject, eventdata, handles)
% --------------------------------------------------------------------
function imPixval_Callback(hObject, eventdata, handles)
if isempty(handles.hText)
handles.hText = impixelinfoval(handles.figure1,handles.ImgObject);
set(handles.hText,'FontWeight','bold');
set(handles.hText,'FontSize',10);
else
    delete(handles.hText);
    handles.hText =[];
end
guidata(hObject, handles);
% --------------------------------------------------------------------
function imHistogram_Callback(hObject, eventdata, handles)
myhistogram(1)
%-------------------------------------------------------------
function imProfile_Callback(hObject, eventdata, handles)
improfile;
title('Intensity profile');
ylabel('Intensity');
xlabel('Pixel position');
%------------------------------------------------
% Histogram function
% -----------------------------------------------
function myhistogram(htype)
img = getimage(gcbf);
if (htype==1)
    figure;hist(double(img(:)),64);
    title('Slice histogram');
    ylabel('Frequency');
    xlabel('Intensity');
elseif (htype == 2)
    figure;hist(log(double(img(:))+1),64);
else
    htype = 0;
end
return;

% --- Executes during object creation, after setting all properties.
function popup_view_CreateFcn(hObject, eventdata, handles)
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
%-------------------------------------------------------------
% --- Executes on selection change in popup_view.
%-------------------------------------------------------------
function popup_view_Callback(hObject, eventdata, handles)
v = get(hObject,'value');
if v ~=1
    nimg = convert_view(handles.oimage.CData,v);
else
    nimg = handles.oimage.CData;
end
[r,c,z]=size(nimg);
handles.image.CData = nimg;
handles.image.Xdata = [1 c];
handles.image.Ydata = [1 r];
handles.viewtype = v;
handles.imSize = [r,c,z];
reset(handles.axes1); axis off;
set(handles.SliceNumSlider,'Min',1,'Max',z,'Value',...
    round(z/2),'SliderStep',[1.0/double(z-1) 1.0/double(z-1)]);
set(handles.ImgObject,'Xdata',[1 c],'Ydata',[1 r]);
% Set the properties of the axes

set(handles.axes1,'XLim',[1 c],'YLim',[1 r],'PlotBoxAspectRatio',[c r 1]);

if ~isempty(handles.imageInfo)
    Sthick = handles.imageInfo.SliceThickness;
    if v == 2
        psize = handles.imageInfo.PixelSpacing(1);
        set(handles.axes1,'DataAspectRatio',[1 psize/Sthick  1]);
    elseif v == 3
        psize = handles.imageInfo.PixelSpacing(2);
        set(handles.axes1,'DataAspectRatio',[1 psize/Sthick  1]);
    else
        psize = handles.imageInfo.PixelSpacing;
        set(handles.axes1,'DataAspectRatio',[psize(1) psize(2)  1]);
    end
end
set(handles.SliceNumEdit,'String',num2str(round(z/2)));
cimg = handles.image.CData(:,:,round(z/2));
handles.currentSliceNumber = round(z/2);
set(handles.ImgObject,'Cdata',cimg);
guidata(hObject, handles);
drawnow
%-------------------------------------------------------------
% function s= convert_view(img,desired_view)
% img: input image
% It assumes that the original image is in transverse view(1)
% 1: transverse, 2: coronal, 3: sagittal views
% s: output image
function s= convert_view(img,desired_view)
[r,c,z]=size(img);
switch desired_view
    case 2
        s = permute(img,[3 2 1]);
        s=flipdim(s,1);
    case 3
        s = permute(img,[3 1 2]);
        s=flipdim(s,1);
    otherwise;
end
% --------------------------------------------------------------------
function readOrdinaryImageSeries_Callback(hObject, eventdata, handles)
[filename, pathname] = uigetfile({'*.*'},'Image');
if isequal(filename,0)|isequal(pathname,0)
    setstatus('File not found');
else
    cd(pathname)
    k = findstr(filename,'.');
    fnames = dir(['*' filename(k+1:end)]);
    timg = imread(filename);
    [r,c,z] = size(timg);
    numSlices = size(fnames,1);
    nimg = zeros([r,c,numSlices]);
    if z == 3
        hw=warndlg('Color images will be converted to gray');
        uiwait;
        nimg(:,:,1) = rgb2gray(timg);
    else
        nimg(:,:,1) = timg;
    end
    hwb = waitbar(0,'Reading Images, please wait...');
    for i=2:numSlices
        waitbar(i/numSlices, hwb);
        timg = imread(fnames(i).name);
        [r,c,z] = size(timg);
        if z == 3
            nimg(:,:,i) = rgb2gray(timg);
        else
            nimg(:,:,i) = timg;
        end
    end
    close(hwb);
    clear timg;
    [r,c,z]=size(nimg);
    handles.image.CData = nimg;
    handles.image.Xdata = [1 c];
    handles.image.Ydata = [1 r];
    handles.imSize = [r,c,z];
    handles.SliceThickness = 1;
    if z==1
        sliderStep = 1;
    else
        sliderStep = 1.0/double(z-1);
    end
    set(handles.SliceNumSlider,'Min',1,'Max',z,'Value',...
        round(z/2),'SliderStep',[sliderStep sliderStep]);
    set(handles.ImgObject,'Xdata',[1 c],'Ydata',[1 r]);
    % Set the properties of the axes
    set(handles.axes1,'XLim',[1 c],'YLim',[1 r],'PlotBoxAspectRatio',[c r 1]);
    set(handles.SliceNumEdit,'String',num2str(round(z/2)));
    cimg = handles.image.CData(:,:,round(z/2));
    set(handles.ImgObject,'Cdata',cimg);
    handles.oimage = handles.image;
    guidata(hObject, handles);
    clear nimg;
    if z > 1
        turnDisplayButtons(handles,'on');
        set(handles.popupmenu_imageModality,'visible','off');
    end
    drawnow
end
% --------------------------------------------------------------------
function saveImageIntoWorkspace_Callback(hObject, eventdata, handles)
assignin('base','IIMG',handles.image.CData);
% --------------------------------------------------------------------
function imfileMenu_Callback(hObject, eventdata, handles)
% --------------------------------------------------------------------
function ImMaps_Callback(hObject, eventdata, handles)
% --------------------------------------------------------------------
% --------------------------------------------------------------------
function readWorkspaceImage_Callback(hObject, eventdata, handles)
Vars = evalin('base','who;');
if ~isempty (Vars)
    [s,v] = listdlg('PromptString','Select an image:',...
        'SelectionMode','single',...
        'Listsize',[150,200],...
        'ListString',Vars);
    if (v)
        imlook3d(evalin('base',Vars{s}));
    end
else
    disp('No workspace variable');
end

% --------------------------------------------------------------------
% Reads series of dicom images into a 3D array
% --------------------------------------------------------------------
function readDicomPETSeries_Callback(hObject, eventdata, handles)
[filename, pathname] = uigetfile({'*.dcm';'*.*'},'Image');
if isequal(filename,0)|isequal(pathname,0)
    setstatus('File not found');
else
    cd(pathname)
    k = findstr(filename,'.');
    fnames = dir([filename(1:2) '*']);
    info = dicominfo(filename);
    r = info.Width;
    c = info.Height;
    numSlices = length(fnames);
    nimg = zeros([r,c,numSlices]);
    hw = waitbar(0,'Reading volume image, please wait...');
    for i=1:numSlices
        info = dicominfo(fnames(i).name);
        nimg(:,:,info.ImageIndex) = info.RescaleSlope*dicomread(info)+info.RescaleIntercept;
        waitbar(i/numSlices,hw);
    end
    close(hw);

    fields = {'Filename','FileModDate','FileSize','ImageType','Width','Height','BitDepth',...
        'NumberOfFrames','Rows','Columns','BitsAllocated','SmallestImagePixelValue',...
        'LargestImagePixelValue','PixelDataGroupLength','StudyDescription','SeriesDescription',...
        'StudyID','ImageID','PatientID','SliceThickness','StudyDate','AcquisitionDate',...
        'AcquisitionTime','PixelSpacing'};

    myinfo = keepfield(info,fields);

    handles.imageInfo = myinfo;

    [r,c,z]=size(nimg);
    handles.image.CData = nimg;
    handles.image.Xdata = [1 c];
    handles.image.Ydata = [1 r];
    handles.imSize = [r,c,z];
    set(handles.SliceNumSlider,'Min',1,'Max',z,'Value',...
        round(z/2),'SliderStep',[1.0/double(z-1) 1.0/double(z-1)]);
    set(handles.ImgObject,'Xdata',[1 c],'Ydata',[1 r]);
    % Set the properties of the axes
    set(handles.axes1,'XLim',[1 c],'YLim',[1 r],'PlotBoxAspectRatio',[c r 1]);
    set(handles.SliceNumEdit,'String',num2str(round(z/2)));
    cimg = handles.image.CData(:,:,round(z/2));
    set(handles.ImgObject,'Cdata',cimg);
    handles.oimage = handles.image;
    guidata(hObject, handles);
    clear nimg;
    if z>1
        turnDisplayButtons(handles,'on')
    end
    handles.imageModality = 1;
    set(handles.popupmenu_imageModality,'visible','on');
    set(handles.popupmenu_imageModality,'value',1);
    drawnow
end
% --------------------------------------------------------------------
function saveCurrentSlice_Callback(hObject, eventdata, handles)
[filename, pathname,filterindex] = uiputfile({'*.png';'*.tif';'*.jpg'},'Image');
fformat = {'png','tif','jpg'};
if isequal(filename,0) | isequal(pathname,0)
    warndlg('No file name entered!');
else
    I= getframe(handles.axes1);
    I=frame2im(I);
    buffer=pwd;
    cd(pathname);
    imwrite(I,filename,fformat{filterindex});
    cd(buffer);
end
% --------------------------------------------------------------------
function adjust_intensity_window_Callback(hObject, eventdata, handles)
setappdata(0  ,  'hMainGui'    , gcf);
hMainGui = getappdata(0, 'hMainGui');
maxINT = max(max(max(handles.image.CData)));
minINT = min(min(min(handles.image.CData)));
cclim = get(handles.axes1,'CLim');
intensityParam = struct('lowerLimit',cclim(1),'upperLimit',cclim(2),'maxIntensity',maxINT, 'minIntensity',minINT);
setappdata(hMainGui,  'intensityParam', intensityParam);
adjustIntensityRangeGui;
uiwait;
hMainGui = getappdata(0,'hMainGui');
handles.intensityParam = getappdata(hMainGui,  'intensityParam');
guidata(hObject, handles);
% -----------------------------------------------------------------
function global_intensity_windowing_Callback(hObject, eventdata, handles)
str = get(gcbo,'Checked');
if strcmp(str,'off')
    handles.global_windowing =1;
    set(gcbo,'Checked','on');
else
    handles.global_windowing =0;
    set(gcbo,'Checked','off');
end
guidata(hObject,handles);
% --------------------------------------------------------------------
function renderVolume_Callback(hObject, eventdata, handles)
img= handles.image.CData;
if size(img,3) > 1
    prompt = {'Enter Lower Threshold:','Enter Upper Threshold:'};
    dlg_title = 'Enter Thresholds'; num_lines = 1; def = {'1','2'};
    answer = inputdlg(prompt,dlg_title,num_lines,def);
    T(1) = str2num(answer{1,1});
    T(2) = str2num(answer{2,1});
    nimg = img >= T(1) & img <= T(2);
    s = size(img);
    vlim=[1 s(2) 1 s(1) 1 s(3)];
    figure;
    % set(get(hndl,'CurrentAxes'),'YDir','rev');
    [x y z D] = subvolume(nimg,vlim);
    p = patch(isosurface(x,y,z,D,0),'FaceColor','r','EdgeColor','none');
    p2 = patch(isocaps(x,y,z,D,0),'FaceColor','interp','EdgeColor','none');
    set(gca,'XTick',[],'YTick',[],'ZTick',[],'Box','on','Color','c');
    isonormals(x,y,z,D,p);
    axis tight;
    daspect([1 1 1])
    lighting gouraud
else
    msgbox('This is not a volume image','Error','modal');
end

% --------------------------------------------------------------------
% READS CT DICOM IMAGES
% --------------------------------------------------------------------
function readDicomCTSeries_Callback(hObject, eventdata, handles)
[filename, pathname] = uigetfile({'*.dcm';'*.*'},'Image');
if isequal(filename,0)|isequal(pathname,0)
    setstatus('File not found');
else
    cd(pathname)
    k = findstr(filename,'.');
    fnames = dir(['*' filename(k(end)+1:end)]);
    info = dicominfo(filename);
    r = info.Width;
    c = info.Height;
    numSlices = length(fnames);
    nimg = zeros([r,c,numSlices]);
    hwb = waitbar(0,'Reading Images, please wait...');
    for i=1:numSlices
        waitbar(i/numSlices, hwb);
        info = dicominfo(fnames(i).name);
        nimg(:,:,i) = info.RescaleSlope*dicomread(info)+info.RescaleIntercept;
    end
    close(hwb);

    fields = {'Filename','FileModDate','FileSize','ImageType','Width','Height','BitDepth',...
        'NumberOfFrames','Rows','Columns','BitsAllocated','SmallestImagePixelValue',...
        'LargestImagePixelValue','PixelDataGroupLength','StudyDescription','SeriesDescription',...
        'StudyID','ImageID','PatientID','SliceThickness','StudyDate','AcquisitionDate','AcquisitionTime','PixelSpacing'};

    myinfo = keepfield(info,fields);
    handles.imageInfo = myinfo;

    [r,c,z]=size(nimg);
    handles.image.CData = nimg;
    handles.image.Xdata = [1 c];
    handles.image.Ydata = [1 r];
    handles.imSize = [r,c,z];
    set(handles.SliceNumSlider,'Min',1,'Max',z,'Value',...
        round(z/2),'SliderStep',[1.0/double(z-1) 1.0/double(z-1)]);
    set(handles.ImgObject,'Xdata',[1 c],'Ydata',[1 r]);
    % Set the properties of the axes
    set(handles.axes1,'XLim',[1 c],'YLim',[1 r],'PlotBoxAspectRatio',[c r 1]);
    set(handles.SliceNumEdit,'String',num2str(round(z/2)));
    cimg = handles.image.CData(:,:,round(z/2));
    set(handles.ImgObject,'Cdata',cimg);
    handles.oimage = handles.image;
    guidata(hObject, handles);
    clear nimg;
    handles.imageModality = 4;
    if z>1
        turnDisplayButtons(handles,'on')
    end
    handles.imageModality = 2;
    set(handles.popupmenu_imageModality,'visible','on');
    set(handles.popupmenu_imageModality,'value',2);
    drawnow
end
% --------------------------------------------------------------------
function imageProcessMenu_Callback(hObject, eventdata, handles)
% --------------------------------------------------------------------
function imageThresholdManual_Callback(hObject, eventdata, handles)
img = getimage(gcbf);
th = cell2mat(inputdlg('Enter Threshold','Threshold',1,{'100'}));
if ~isempty(th)
    th = str2num(th);
    figure;imagesc(img>th);axis image;axis off
    title('Thresholded Image');
end
% --------------------------------------------------------------------
function automatedThresholding_Callback(hObject, eventdata, handles)
img = getimage(gcbf);
th=otsu_th(single(img(:)),64);
figure;imagesc(img>th);axis image;axis off
title('Thresholded Image');
% --------------------------------------------------------------------
% READS DICOM IMAGE SINGLE OR VOLUME
% --------------------------------------------------------------------
function readSingleDicomImage_Callback(hObject, eventdata, handles)
[filename, pathname] = uigetfile({'*.dcm';'*.*'},'Image');
if isequal(filename,0)|isequal(pathname,0)
    setstatus('File not found');
else
    cd(pathname)
    info = dicominfo(filename);
    if isfield(info,'NumberOfSlices')
        z = info.NumberOfSlices;
    elseif isfield(info,'NumberOfFrames')
        z = info.NumberOfFrames;
    end
    if isfield(info,'SliceThickness')
        handles.SliceThickness= info.SliceThickness;
    end
    % Reading dicom image
    nimg = squeeze(dicomread(info));
    [r,c,z]=size(nimg);
    handles.image.CData = nimg;
    handles.image.Xdata = [1 c];
    handles.image.Ydata = [1 r];
    handles.imSize = [r,c,z];

    if z >1
        turnDisplayButtons(handles,'on');
        set(handles.SliceNumSlider,'Min',1,'Max',z,'Value',...
            round(z/2),'SliderStep',[1.0/double(z-1) 1.0/double(z-1)]);
        set(handles.SliceNumEdit,'String',num2str(round(z/2)));
        turnDisplayButtons(handles,'on');
        cimg = handles.image.CData(:,:,round(z/2));
    else
        turnDisplayButtons(handles,'off');
        cimg = handles.image.CData;
    end
    fields = {'Filename','FileModDate','FileSize','ImageType','Width','Height','BitDepth',...
        'NumberOfFrames','Rows','Columns','BitsAllocated','SmallestImagePixelValue',...
        'LargestImagePixelValue','PixelDataGroupLength','StudyDescription','SeriesDescription',...
        'StudyID','ImageID','PatientID','SliceThickness','StudyDate','AcquisitionDate','AcquisitionTime','PixelSpacing'};
    myinfo = keepfield(info,fields);
    handles.imageInfo = myinfo;

    % Set the properties of the axes and image object
    set(handles.ImgObject,'XData',[1 c],'YData',[1 r],'CData',cimg);
    set(handles.axes1,'XLim',[1 c],'YLim',[1 r],'PlotBoxAspectRatio',[c r 1],'DataAspectRatio',...
        [1 info.PixelSpacing(2)/info.PixelSpacing(1)  1]);
    handles.oimage = handles.image;
    guidata(hObject, handles);
    clear nimg;
    drawnow
end
% --------------------------------------------------------------------
% --- Executes on selection change in popup_flip.
function popup_flip_Callback(hObject, eventdata, handles)
v = get(hObject,'value');
switch v
    case 1
       handles.image.CData = flipdim(handles.image.CData,2);
    case 2
        handles.image.CData = flipdim(handles.image.CData,1);
    case 3
        handles.image.CData = flipdim(handles.image.CData,3);
    case 4
    handles.image.CData = handles.oimage.CData;
    otherwise
end
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function popup_flip_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function turnDisplayButtons(handles,bmode)
set(handles.SliceNumSlider,'visible',bmode);
set(handles.SliceNumEdit,'visible',bmode);
set(handles.popup_view,'visible',bmode);
set(handles.popup_flip,'visible',bmode);
set(handles.text1,'visible',bmode);
set(handles.text2,'visible',bmode);
set(handles.text3,'visible',bmode);
set(handles.popup_flip,'visible',bmode);
return

% --------------------------------------------------------------------
function helpCallback_Callback(hObject, eventdata, handles)
% --------------------------------------------------------------------
% --------------------------------------------------------------------
function helpMenu_Callback(hObject, eventdata, handles)
% --------------------------------------------------------------------
function helpAboutImlook3d_Callback(hObject, eventdata, handles)
s = sprintf(['IMLOOK3D: 3D Image Display Tool\n',...
    'Version 1.2, Copyright (2006 - Inf and Beyond)\n',...
    'Omer Demirkaya, KFSH&RC']);
msgbox(s,'About IMLOOK3D Tool','modal');
% --------------------------------------------------------------------
function helpDisclaimer_Callback(hObject, eventdata, handles)
s = {'This SOFTWARE DOES NOT (A) MAKE ANY WARRANTY, EXPRESS OR IMPLIED',
    'WITH RESPECT TO THE USE OF THE INFORMATION PROVIDED HEREBY; NOR (B)',
    'GUARANTEE THE ACCURACY,COMPLETENESS,USEFULNESS,OR ADEQUACY OF ANY',
    'RESOURCES, INFORMATION,SYSTEM,PRODUCT,OR PROCESS AVAILABLE - AT OR',
    'THROUGH THIS SOFTWARE.'};
msgbox(s,'Disclaimer','modal');

% --------------------------------------------------------------------
function imageInformation_Callback(hObject, eventdata, handles)
if ~isempty(handles.imageInfo)
    imageinfo(handles.imageInfo);
else
    msgbox('There is no image or no information about it','Image Info..','modal');
end

% --------------------------------------------------------------------
function readSingleOrdinaryImage_Callback(hObject, eventdata, handles)
[filename, pathname] = uigetfile({'*.*'},'Image');
if isequal(filename,0)|isequal(pathname,0)
    setstatus('File not found');
else
    cd(pathname)
    info = imfinfo(filename);
    [img, map]=imread(filename);
    [r,c,z] = size(img);
    handles.image.CData = img;
    handles.image.Xdata = [1 c];
    handles.image.Ydata = [1 r];
    handles.imSize = [r,c,z];

    if z > 1   % single slice
    set(handles.SliceNumSlider,'Min',1,'Max',z,'Value',...
        round(z/2),'SliderStep',[1.0/double(z-1) 1.0/double(z-1)]);
    set(handles.SliceNumEdit,'String',num2str(round(z/2)));
    turnDisplayButtons(handles,'on');
    cimg = handles.image.CData(:,:,round(z/2));
    else
        turnDisplayButtons(handles,'off');
        cimg = handles.image.CData;
    end
    handles.imageInfo = info;
   % Set the properties of the axes and image object
    set(handles.ImgObject,'XData',[1 c],'YData',[1 r],'CData',cimg);
    set(handles.axes1,'XLim',[1 c],'YLim',[1 r],'PlotBoxAspectRatio',[c r 1],'DataAspectRatio',[1 1 1]);
    handles.oimage = handles.image;
    guidata(hObject, handles);
    clear nimg;
    drawnow
end
% --------------------------------------------------------------------
% SMOOTH IMAGE
% --------------------------------------------------------------------
function smoothVolume_Callback(hObject, eventdata, handles)
if size(handles.image.CData,3) > 1
h = waitbar(0,'Smoothing, please wait...');
handles.image.CData= smooth3(handles.image.CData,'gaussian');
close(h);
set(handles.ImgObject,'Cdata',handles.image.CData(:,:,handles.currentSliceNumber));
guidata(hObject,handles);
end
% --------------------------------------------------------------------
% Voxel dimensions
% --------------------------------------------------------------------
function voxelSize_Callback(hObject, eventdata, handles)
prompt = {'X Size:','Y Size:', 'Z Size'};
dlg_title = 'Enter Voxel Sizes'; num_lines = 1; def = {'1','1','1'};
answer = inputdlg(prompt,dlg_title,num_lines,def);
if ~isempty(answer)
    voxdim(1) = str2num(cell2mat(answer(1,1)));
    voxdim(2) = str2num(cell2mat(answer(2,1)));
    voxdim(3) = str2num(cell2mat(answer(3,1)));
    bw = voxdim > 0; % checking for zero or negative
    if prod(double(bw))
        handles.imageInfo.PixelSpacing(1) = voxdim(1);
        handles.imageInfo.PixelSpacing(2) = voxdim(2);
        handles.imageInfo.SliceThickness  = voxdim(3);
        if handles.viewtype == 2
            set(handles.axes1,'DataAspectRatio',[1 voxdim(1)/voxdim(3)  1]);
        elseif handles.viewtype == 3
            set(handles.axes1,'DataAspectRatio',[1 voxdim(2)/voxdim(3)  1]);
        else
            set(handles.axes1,'DataAspectRatio',[1 voxdim(2)/voxdim(1)  1]);
        end

    end
end
guidata(hObject,handles);
% --------------------------------------------------------------------
% Horizontal Profile
% --------------------------------------------------------------------
function horizontalProfile_Callback(hObject, eventdata, handles)
h = warndlg('Double-click on the row','Pixel Selection');
uiwait;
[c,r,P] = impixel;
f = get(findobj(gcf,'Type','image'));
txt = ['Profile for row: ' num2str(r)];
figure;
plot(f.CData(r,:));title(txt);xlabel('Pixel location');ylabel('Intensity');
% --------------------------------------------------------------------
% Vertical Profile
% --------------------------------------------------------------------
function verticalProfile_Callback(hObject, eventdata, handles)
h = warndlg('Double-click on the column','Pixel Selection');
uiwait;
[c,r,P] = impixel;
f = get(findobj(gcf,'Type','image'));
txt = ['Profile for column: ' num2str(c)];
figure;
plot(f.CData(:,c));title(txt);xlabel('Pixel location');ylabel('Intensity');

% --------------------------------------------------------------------
% --- Executes during object creation, after setting all properties.
% --------------------------------------------------------------------
function popupmenu_imageModality_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
% --------------------------------------------------------------------
% --- Image modality selection
% --------------------------------------------------------------------
function popupmenu_imageModality_Callback(hObject, eventdata, handles)
imMod = get(hObject,'value');
handles.imageModality = imMod;

% --------------------------------------------------------------------
%  Executes MATLAB IMTOOL
% --------------------------------------------------------------------
function imtoolMATLAB_Callback(hObject, eventdata, handles)
f = get(findobj(gcf,'Type','image'));
imtool(f.CData,[min(f.CData(:)) max(f.CData(:))]);
% --------------------------------------------------------------------
%  helpMenu on Imlook3D
% --------------------------------------------------------------------
function helpIMLOOK3D_Callback(hObject, eventdata, handles)
warndlg('Under Construction and needs help :)');
% --------------------------------------------------------------------
%  Reading MRI images
% --------------------------------------------------------------------
function readDICOMMRISeries_Callback(hObject, eventdata, handles)
[filename, pathname] = uigetfile({'*.dcm';'*.*'},'Image');
if isequal(filename,0)|isequal(pathname,0)
    setstatus('File not found');
else
    cd(pathname)
    k = findstr(filename,'.');
    fnames = dir(['*' filename(k+1:end)]);
    info = dicominfo(filename);
    r = info.Width;
    c = info.Height;
    numSlices = length(fnames);
    nimg = zeros([r,c,numSlices]);
    hwb = waitbar(0,'Reading Images, please wait...');
    for i=1:numSlices
        waitbar(i/numSlices, hwb);
        info = dicominfo(fnames(i).name);
        nimg(:,:,i) = dicomread(info);
    end
    close(hwb);

    fields = {'Filename','FileModDate','FileSize','ImageType','Width','Height','BitDepth',...
        'NumberOfFrames','Rows','Columns','BitsAllocated','SmallestImagePixelValue',...
        'LargestImagePixelValue','PixelDataGroupLength','StudyDescription','SeriesDescription',...
        'StudyID','ImageID','PatientID','SliceThickness','StudyDate','AcquisitionDate','AcquisitionTime','PixelSpacing'};

    myinfo = keepfield(info,fields);
    handles.imageInfo = myinfo;

    [r,c,z]=size(nimg);
    handles.image.CData = nimg;
    handles.image.Xdata = [1 c];
    handles.image.Ydata = [1 r];
    handles.imSize = [r,c,z];
    set(handles.SliceNumSlider,'Min',1,'Max',z,'Value',...
        round(z/2),'SliderStep',[1.0/double(z-1) 1.0/double(z-1)]);
    set(handles.ImgObject,'Xdata',[1 c],'Ydata',[1 r]);
    % Set the properties of the axes
    set(handles.axes1,'XLim',[1 c],'YLim',[1 r],'PlotBoxAspectRatio',[c r 1]);
    set(handles.SliceNumEdit,'String',num2str(round(z/2)));
    cimg = handles.image.CData(:,:,round(z/2));
    set(handles.ImgObject,'Cdata',cimg);
    handles.oimage = handles.image;
    guidata(hObject, handles);
    clear nimg;
    if z>1
        turnDisplayButtons(handles,'on')
    end
    handles.imageModality = 4;
    set(handles.popupmenu_imageModality,'visible','on');
    set(handles.popupmenu_imageModality,'value',4);
    drawnow
end

% --------------------------------------------------------------------
function cropVolume_Callback(hObject, eventdata, handles)
% hObject    handle to cropVolume (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

[I,rect] = imcrop;
s = handles.image.CData;
switch handles.viewtype
    case 1
        [r,c,z]  = size(s);
        oimg = crop3d(s,rect,[1 z]);
    case 2
%         s = permute(img,[3 2 1]);
%         s = flipdim(s,1);
        [r,c,z]  = size(s);
        oimg = crop3d(s,rect,[1 z]);
    case 3
%         s = permute(img,[3 1 2]);
%         s = flipdim(s,1);
        [r,c,z]  = size(s);
        oimg = crop3d(s,rect,[1 z]);
    otherwise;
end
imlook3d(oimg);

author: Sable

repo: mcbench-benchmarks


text:
% We consider a switching Kalman filter of the kind studied
% by Zoubin Ghahramani, i.e., where the switch node determines
% which of the hidden chains we get to observe (data association).
% e.g., for n=2 chains
%
% X1 -> X1
% | X2 -> X2
% \ |
%  v
%  Y
%  ^
%  |
%  S
%
% Y is a gmux (multiplexer) node, where S switches in one of the parents.
% We differ from Zoubin by not connecting the S nodes over time (which
% doesn't make sense for data association).
% Indeed, we assume the S nodes are always observed.
%
%
% We will track 2 objects (points) moving in the plane, as in BNT/Kalman/tracking_demo.
% We will alternate between observing them.

nobj = 2;
N = nobj+2;
Xs = 1:nobj;
S = nobj+1;
Y = nobj+2;

intra = zeros(N,N);
inter = zeros(N,N);
intra([Xs S], Y) =1;
for i=1:nobj
  inter(Xs(i), Xs(i))=1;
end

Xsz = 4; % state space = (x y xdot ydot)
Ysz = 2;
ns = zeros(1,N);
ns(Xs) = Xsz;
ns(Y) = Ysz;
ns(S) = n;

bnet = mk_dbn(intra, inter, ns, 'discrete', S, 'observed', [S Y]);

% For each object, we have
% X(t+1) = F X(t) + noise(Q)
% Y(t) = H X(t) + noise(R)
F = [1 0 1 0; 0 1 0 1; 0 0 1 0; 0 0 0 1];
H = [1 0 0 0; 0 1 0 0];
Q = 1e-3*eye(Xsz);
%R = 1e-3*eye(Ysz);
R = eye(Ysz);

% We initialise object 1 moving to the right, and object 2 moving to the left
% (Here, we assume nobj=2)
init_state{1} = [10 10 1 0]';
init_state{2} = [10 -10 -1 0]';

for i=1:nobj
  bnet.CPD{Xs(i)} = gaussian_CPD(bnet, Xs(i), 'mean', init_state{i}, 'cov', 1e-4*eye(Xsz));
end
bnet.CPD{S} = root_CPD(bnet, S); % always observed
bnet.CPD{Y} = gmux_CPD(bnet, Y, 'cov', repmat(R, [1 1 nobj]), 'weights', repmat(H, [1 1 nobj]));
% slice 2
eclass = bnet.equiv_class;
for i=1:nobj
  bnet.CPD{eclass(Xs(i), 2)} = gaussian_CPD(bnet, Xs(i)+N, 'mean', zeros(Xsz,1), 'cov', Q, 'weights', F);
end

% Observe objects at random
T = 10;
evidence = cell(N, T);
data_assoc = sample_discrete(normalise(ones(1,nobj)), 1, T);
evidence(S,:) = num2cell(data_assoc);
evidence = sample_dbn(bnet, 'evidence', evidence);

% plot the data
true_state = cell(1,nobj);
for i=1:nobj
  true_state{i} = cell2num(evidence(Xs(i), :)); % true_state{i}(:,t) = [x y xdot ydot]'
end
obs_pos = cell2num(evidence(Y,:));
figure(1)
clf
hold on
styles = {'rx', 'go', 'b+', 'k*'};
for i=1:nobj
  plot(true_state{i}(1,:), true_state{i}(2,:), styles{i});
end
for t=1:T
  text(obs_pos(1,t), obs_pos(2,t), sprintf('%d', t));
end
hold off
relax_axes(0.1)

% Inference
ev = cell(N,T);
ev(bnet.observed,:) = evidence(bnet.observed, :);

engines = {};
engines{end+1} = jtree_dbn_inf_engine(bnet);
%engines{end+1} = scg_unrolled_dbn_inf_engine(bnet, T);
engines{end+1} = pearl_unrolled_dbn_inf_engine(bnet);
E = length(engines);

inferred_state = cell(nobj,E); % inferred_state{i,e}(:,t)
for e=1:E
  engines{e} = enter_evidence(engines{e}, ev);
  for i=1:nobj
    inferred_state{i,e} = zeros(4, T);
    for t=1:T
      m = marginal_nodes(engines{e}, Xs(i), t);
      inferred_state{i,e}(:,t) = m.mu;
    end
  end
end
inferred_state{1,1}
inferred_state{1,2}

% Plot results
figure(2)
clf
hold on
styles = {'rx', 'go', 'b+', 'k*'};
nstyles = length(styles);
c = 1;
for e=1:E
  for i=1:nobj
    plot(inferred_state{i,e}(1,:), inferred_state{i,e}(2,:), styles{mod(c-1,nstyles)+1});
    c = c + 1;
  end
end
for t=1:T
  text(obs_pos(1,t), obs_pos(2,t), sprintf('%d', t));
end
hold off
relax_axes(0.1)


text:
function  data = datasaver(varargin)

persistent saved

if nargin > 3
    saved.savedfunctiondata = varargin{1};
    saved.savedfunctionGradientdata = varargin{2};
    saved.savedfunctionHessiandata = varargin{3};
    saved.savedconstraintdata = varargin{4};
    saved.savedconstraintGradientdata = varargin{5};
    saved.savedconstraintHessiandata = varargin{6};
    saved.m = size(saved.savedconstraintdata.F_struc,1);
    saved.n = length(saved.savedfunctiondata.linearindicies);
    saved.monomials = [];
    saved.x = randn(saved.n,1);
    return
end

x = varargin{2};
if ~isequal(x,saved.x)
    % Recompute the monomial terms for new x
    saved.monomials=ones(1,size(saved.savedfunctiondata.monomtable,2));
    saved.monomials(saved.savedfunctiondata.linearindicies) = x(:)';
   % changed = find(x-saved.x);
   % changed = saved.savedfunctiondata.linearindicies(changed);
   % changedmonomials = find(any(saved.savedfunctiondata.monomtable(saved.savedfunctiondata.nonlinearindicies,changed),2));

    for i = 1:length(saved.savedfunctiondata.nonlinearindicies)
        saved.monomials(saved.savedfunctiondata.nonlinearindicies(i)) = prod(saved.monomials.^saved.savedfunctiondata.monomtable(saved.savedfunctiondata.nonlinearindicies(i),:));
    end
    saved.x = x;
end

switch varargin{1}
    case 1
        F = saved.savedfunctiondata.F_struc;
        x = varargin{2};
    case 2
        F = saved.savedfunctionGradientdata.F_struc;
    case 3
        F = saved.savedfunctionHessiandata.F_struc;
    case 4
        if nargin == 2
            % ipopt call
            F = saved.savedconstraintdata.F_struc;
        else
            % pennlp call
            F = saved.savedconstraintdata.F_struc(1+varargin{3},:);
        end
    case 5
        F = saved.savedconstraintGradientdata.F_struc;
        start = 1+(varargin{3}-1)*saved.n;
        F = F(varargin{3}:saved.m:end,:);
    case 6
        F = saved.savedconstraintHessiandata.F_struc;
        F = full(F((varargin{3}-1)+(1:saved.m:saved.m*saved.n^2),:));
    otherwise
end

data = (F*[1;saved.monomials(:)]);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function Y = minfilt2(X,varargin)
%  MINFILT2    Two-dimensional min filter
%
%     Y = MINFILT2(X,[M N]) performs two-dimensional minimum
%     filtering on the image X using an M-by-N window. The result
%     Y contains the minimun value in the M-by-N neighborhood around
%     each pixel in the original image.
%     This function uses the van Herk algorithm for min filters.
%
%     Y = MINFILT2(X,M) is the same as Y = MINFILT2(X,[M M])
%
%     Y = MINFILT2(X) uses a 3-by-3 neighborhood.
%
%     Y = MINFILT2(..., 'shape') returns a subsection of the 2D
%     filtering specified by 'shape' :
%        'full'  - Returns the full filtering result,
%        'same'  - (default) Returns the central filter area that is the
%                   same size as X,
%        'valid' - Returns only the area where no filter elements are outside
%                  the image.
%
%     See also : MAXFILT2, VANHERK
%

% Initialization
[S, shape] = parse_inputs(varargin{:});

% filtering
Y = vanherk(X,S(1),'min',shape);
Y = vanherk(Y,S(2),'min','col',shape);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [S, shape] = parse_inputs(varargin)
shape = 'same';
flag = [0 0]; % size shape

for i = 1 : nargin
   t = varargin{i};
   if strcmp(t,'full') & flag(2) == 0
      shape = 'full';
      flag(2) = 1;
   elseif strcmp(t,'same') & flag(2) == 0
      shape = 'same';
      flag(2) = 1;
   elseif strcmp(t,'valid') & flag(2) == 0
      shape = 'valid';
      flag(2) = 1;
   elseif flag(1) == 0
      S = t;
      flag(1) = 1;
   else
      error(['Too many / Unkown parameter : ' t ])
   end
end

if flag(1) == 0
   S = [3 3];
end
if length(S) == 1;
   S(2) = S(1);
end
if length(S) ~= 2
   error('Wrong window size parameter.')
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function u = unicycle_index_to_sequence ( n, u_index )

%*****************************************************************************80
%
%% UNICYCLE_INDEX_TO_SEQUENCE converts a unicycle from index to sequence form.
%
%  Example:
%
%    N = 4
%
%    U_INDEX = 3 1 4 2
%    U       = 1 3 4 2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the unicycles.
%
%    Output, integer U_INDEX(N), the unicycle index vector.
%
%    Input, integer U(N), the unicycle sequence vector.
%
  u(1) = 1;
  i = 1;

  for j = 2 : n

    i = u_index(i);
    u(j) = i;

    if ( i == 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'UNICYCLE_INDEX_TO_SEQUENCE - Fatal error!\n' );
      fprintf ( 1, '  The index vector does not represent a unicycle.\n' );
      fprintf ( 1, '  On step %d, u_index(%d) = 1.\n', j, i );
      error ( 'UNICYCLE_INDEX_TO_SEQUENCE - Fatal error!' );
    end

  end

  return
end


text:
function renamefolder(HOMEANNOTATIONS, annotationfolder)
%
% Writes in the folder field of the xml file the nane of the folder in
% which the annotation files are.
%
% This script is useful when moving annotation files into other folders
% manually. This function guaratees consistency.

files = dir(fullfile(HOMEANNOTATIONS, annotationfolder, '*.xml'));

Nfiles = length(files);
for i = 1:Nfiles
    v = loadXML(fullfile(HOMEANNOTATIONS, annotationfolder, files(i).name));
    v.annotation.folder = annotationfolder;
    writeXML(fullfile(HOMEANNOTATIONS, annotationfolder, files(i).name), v);
end

author: CSAILVision

repo: LabelMeToolbox


text:
function geometry_test0023 ( )

%*****************************************************************************80
%
%% TEST0023 tests ANGLE_HALF_2D;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0023\n' );
  fprintf ( 1, '  ANGLE_HALF_2D computes the half angle between two rays;\n' );
  fprintf ( 1, '  The angle is defined by the points (P1,P2,P3)\n' );
  fprintf ( 1, '  or by the rays P2-->P3, P2-->P1.\n' );

  p2(1:2,1) = [ 5.0; 3.0 ];

  angle_deg = 75.0;
  r = 3.0;
  p1(1,1) = p2(1,1) + r * cos_deg ( angle_deg );
  p1(2,1) = p2(2,1) + r * sin_deg ( angle_deg );

  angle_deg = 15.0;
  r = 2.0;
  p3(1,1) = p2(1,1) + r * cos_deg ( angle_deg );
  p3(2,1) = p2(2,1) + r * sin_deg ( angle_deg );

  r8vec_print ( 2, p1, '  Point P1:' );
  r8vec_print ( 2, p2, '  Point P2:' );
  r8vec_print ( 2, p3, '  Point P3:' );

  p4 = angle_half_2d ( p1, p2, p3 );

  r8vec_print ( 2, p4, ...
    '  End point of unit ray from P2, defining half angle, P4:' );

  angle_deg = 45.0;
  r = 1.0;
  p4(1,1) = p2(1,1) + r * cos_deg ( angle_deg );
  p4(2,1) = p2(2,1) + r * sin_deg ( angle_deg );

  r8vec_print ( 2, p4, '  Expected value of P4:' );

  return
end


text:
% This file performs the dual channel speech enhancement based on GSC beamformer and postfilter.
% Refer to "Analysis of Two-Channel Generalized Sidelobe Canceller (GSC) With Post-Filtering".

clear
clc

%% Load audio file

audioname = 'carmix';
[x, fs0 ] = audioread([audioname,'.wav']);
fs      = 16e3;
x       = resample(x,fs,fs0);
Ntime   = size(x,1);
Nmic    = size(x,2);
if Nmic ~= 2
    error('The number of channels does not equal 2!')
end

%% Set parameters

NFFT        = 256;
f           = 0:fs/NFFT:fs/2;
full        = 0:fs/NFFT:fs-fs/NFFT;
len         = 17e-2;
fragsize    = 256;
overlap     = 0.5;
c           = 340;
lamdaB0     = 1.54;
omegalow    = 1;
omegahigh   = 3;
gamma0      = 4.6;
alpha       = 0.92;
alphas      = 0.80;
alphadconst = 0.85;
KL          = 9;
KR          = 50;
phi0        = 0.25;
% b           = [1 3 4 3 1] / 12;
b           = hanning(5).';
b           = b/sum(b);
Gmin        = 10^(-20/20); % -20dB
vmax        = 5;
EPIS        = 1E-10;
angle1      = 0;
angle2      = 0.1;
xia         = 0.98;
ximin       = 10^(-15/10);

%% set window

win_1       = hanning(fragsize).';
win_2       = hanning(NFFT).';
win_pre     = sqrt(win_1);
win_post    = sqrt(win_2);
win_pre     = sqrt(win_pre);
win_post    = sqrt(win_post);
scale_fac   = sqrt(NFFT/sum(win_1));
scale_postfac = 1.0 / sqrt(NFFT/sum(win_2));
win_pre     = win_pre * scale_fac;
win_post    = win_post * scale_postfac;

%% Cutframes

% ptout = zeros(Nmic, Nframe, fragsize);
for m = 1:Nmic
    [ptout(m,:,:), Nframe] = cutframe(x(:,m), fragsize, overlap);
end

%% Process

Z       = zeros(Nmic,NFFT/2+1);
pt      = zeros(Nmic,fragsize);
pxx     = zeros(NFFT,Nmic);
q       = ones(Nframe,NFFT/2+1);
p       = zeros(Nframe,NFFT/2+1);
yxx     = zeros(Nframe,NFFT/2+1);
yt      = zeros(Nframe,NFFT);
alphad  = zeros(Nframe,NFFT/2+1);
v       = zeros(Nframe,NFFT/2+1);
gamma   = ones(Nframe,NFFT/2+1);
gammas  = ones(Nframe,NFFT/2+1);
xi      = zeros(Nframe,NFFT/2+1);
phi     = zeros(Nframe,NFFT/2+1);
SY      = zeros(Nframe,NFFT);
MY      = zeros(Nframe,NFFT);
SU      = zeros(Nframe,NFFT);
MU      = zeros(Nframe,NFFT);
lamdad  = zeros(Nframe,NFFT/2+1);
GH1     = zeros(Nframe,NFFT/2+1);
G       = ones(Nframe,NFFT/2+1);
lamdaBY = zeros(Nframe,NFFT/2+1);
lamdaBU = zeros(Nframe,NFFT/2+1);
omega   = zeros(Nframe,NFFT/2+1);
phi_tilt = zeros(1,Nframe);
YMAT    = zeros(Nframe,NFFT);
UMAT    = zeros(Nframe,NFFT);

deltak  = 2*pi*full * len * sin(angle1*pi/180) / c + angle2*pi/180;
W       = 0.5 * conj([exp(1i*deltak/2); exp(-1i*deltak/2)]);
B       = 0.5 * conj([exp(1i*deltak/2); -exp(-1i*deltak/2)]);
coh     = sinc(2*pi*full*len/c);       % diffused noise field
H       = 1i*imag(exp(1i*deltak).*coh) ./ (1-real(exp(1i*deltak).*coh));

% compute OMEGA_HIGH
% sincwl = min(sinc(2*pi*f*len/c), 0.998);
% OMEGA_HIGH = 0.57 * cot(deltak(1:NFFT/2+1)/2).^2 .* (1-sincwl).^2 ./ max((1-sincwl.*cos(deltak(1:NFFT/2+1))).^2, EPIS);
% OMEGA_HIGH = 0.57 * cot(deltak(1:NFFT/2+1)/2).^2;
% OMEGA_HIGH = 0.57 * sin(2*pi*f*len*sin(0)/c).^2 ./ (sin(deltak(1:NFFT/2+1)/2).^2 .* sin(2*pi*f*len*sin(0)/c-deltak(1:NFFT/2+1)/2).^2);

hbar = waitbar(0,'Please wait...');
for n = 1:Nframe
    waitbar(n/Nframe,hbar);

%     % ***************   GSC part   **************************************
%     pt1 = [zeros(1,(NFFT-fragsize)/2), reshape(ptout(1,n,:),1,fragsize), zeros(1,(NFFT-fragsize)/2)];
%     pt2 = [zeros(1,(NFFT-fragsize)/2), reshape(ptout(2,n,:),1,fragsize), zeros(1,(NFFT-fragsize)/2)];
%     pt  = [ pt1 .* win_pre.' ; pt2 .* win_pre.'];
    pt1 = [reshape(ptout(1,n,:),1,fragsize) .* win_pre, zeros(1,NFFT-fragsize)];
    pt2 = [reshape(ptout(2,n,:),1,fragsize) .* win_pre, zeros(1,NFFT-fragsize)];
    pt  = [ pt1 .' pt2.' ];
    pxx = fft(pt,NFFT);
    Z   = pxx.';
    Y   = sum(conj(W) .* Z) - conj(H) .* sum(conj(B) .* Z);
    U   = sum(conj(B) .* Z);

    % **********   Postfilter part   *************************************
    if n == 1
        SY(n,:)     = abs(Y).^2;
        MY(n,:)     = abs(Y).^2;
        lamdad(n,:) = abs(Y(1:NFFT/2+1)).^2;
        SU(n,:)     = abs(U).^2;
        MU(n,:)     = abs(U).^2;
        GH1(n,:)    = ones(1,NFFT/2+1);
        gamma(n,:)  = ones(1,NFFT/2+1);

        param_Y = initialise_parameters(MY(n,:).',fs,'imcra');
        param_U = initialise_parameters(MU(n,:).',fs,'imcra');

        alphad(n,:)     = alphadconst + (1-alphadconst) * p(n,:);
        lamdad(n+1,:)   = alphad(n,:) .* lamdad(n,:) + (1-alphad(n,:)) .* abs(Y(1:NFFT/2+1)).^2;

        xi(n,:) = xia + (1-xia) * max(gamma(n,:)-1,0);     % initialization of a priori SNR
    else
        Ymat    = buildmat(Y,length(b));
        Umat    = buildmat(U,length(b));
        SY(n,:) = alphas * SY(n-1,:) + (1-alphas) * b * abs(Ymat).^2;
        SU(n,:) = alphas * SU(n-1,:) + (1-alphas) * b * abs(Umat).^2;
        param_Y = noise_estimation((abs(Y).^2).','imcra',param_Y);
        param_U = noise_estimation((abs(U).^2).','imcra',param_U);
        MY(n,:) = param_Y.noise_ps.';
        MU(n,:) = param_U.noise_ps.';

        lamdaBY(n,:)    = SY(n,1:NFFT/2+1) ./ max(MY(n,1:NFFT/2+1), EPIS);
        lamdaBU(n,:)    = SU(n,1:NFFT/2+1) ./ max(MU(n,1:NFFT/2+1), EPIS);

%         omega(n,:)      = max(SY(n,1:NFFT/2+1) - MY(n,1:NFFT/2+1), 0) ./ max(SU(n,1:NFFT/2+1) - MU(n,1:NFFT/2+1), 0.001*MY(n,1:NFFT/2+1));
        omega(n,:)      = max(SY(n,1:NFFT/2+1) - MY(n,1:NFFT/2+1), 0) ./ max(SU(n,1:NFFT/2+1) - MU(n,1:NFFT/2+1), EPIS);

        % decide q according to fig.4
        for k = 1:NFFT/2+1
            if lamdaBY(n,k) > lamdaB0
                if lamdaBU(n,k) > lamdaB0
                    % ***** compute phi, using eq.33 ***********
                    if omega(n,k) <= omegalow
                        phi(n,k) = 0;
                    else
                        if omegalow < omega(n,k) && omega(n,k) <= omegahigh      % OMEGA_HIGH(k)       % replace omegahigh with OMEGA_HIGH(k)
                            phi(n,k) = (omega(n,k) - omegalow) / (omegahigh - omegalow);
%                         if omegalow < omega(n,k) <= OMEGA_HIGH(k)     % OMEGA_HIGH(k)       % replace omegahigh with OMEGA_HIGH(k)
%                             phi(n,k) = (omega(n,k) - omegalow) / (OMEGA_HIGH(k) - omegalow);
                        else
                            phi(n,k) = 1;
                        end
                    end
                    % ******************************************
                else
                    phi(n,k) = 1;
                end
            else
                phi(n,k) = 0;
            end
        end

        % compute eq.34
        phi_tilt(n) = sum(phi(n,KL:KR)) / (KR-KL+1);

        gammas(n,:) = abs(Y(1:NFFT/2+1)).^2 ./ MY(n,1:NFFT/2+1);

        % compute q, using eq.35
        if phi_tilt(n) > phi0
            for k = 1:NFFT/2+1
                if (gammas(n,k) <= 1)  || (phi_tilt(n) <= phi0)
                    q(n,k) = 1;
                else
                    q(n,k) = max( (gamma0-gammas(n,k))/(gamma0-1) , 1-phi(n,k) );
                end
            end
        else
            q(n,:) = ones(1,NFFT/2+1);
        end

        q(n,:) = min(q(n,:),1);
        q(n,:) = max(q(n,:),0);

        % compute xi, using eq.37
        gamma(n,:) = abs(Y(1:NFFT/2+1)).^2 ./ max(lamdad(n,1:NFFT/2+1), EPIS);

        xi(n,:) = alpha * GH1(n-1,:).^2 .* gamma(n-1,:) + (1-alpha) * max(gamma(n,:)-1,0);
        xi(n,:) = max(xi(n,:), ximin);
        v(n,:)  = gamma(n,:) .* xi(n,:) ./ (1+xi(n,:));

        % compute GH1, using eq.38
        GH1(n,:)        = ones(1,NFFT/2+1);
        index           = find(v(n,:) > vmax);
        GH1(n,index)    = xi(n,index) ./ (1+xi(n,index));
        index           = find(0<v(n,:) & v(n,:)<=vmax);
        GH1(n,index)    = xi(n,index) .* exp(0.5*expint(v(n,index)))./ (1+xi(n,index));

        % compute p, using eq.36
        p(n,:) = 1 ./ (1 + q(n,:) .* (1+xi(n,:)).*exp(-v(n,:)) ./ (1-q(n,:)+EPIS));
        p(n,:) = min(p(n,:), 1);
        p(n,:) = max(p(n,:), 0);

        % compute eq.39 eq.40
        alphad(n,:) = alphadconst + (1-alphadconst) * p(n,:);
        if n < Nframe
            lamdad(n+1,:) = alphad(n,:) .* lamdad(n,:) + (1-alphad(n,:)) .* abs(Y(1:NFFT/2+1)).^2;
        end

    end

    % apply filter and ifft
    G(n,:)      = (GH1(n,:).^p(n,:)) .* (Gmin.^(1-p(n,:)));
    yxx(n,:)    = G(n,:) .* Y(1:NFFT/2+1);
    yt(n,:)     = real(ifft([yxx(n,:),conj(fliplr(yxx(n,2:end-1)))],NFFT));
    yt(n,:)     = yt(n,:) .* win_post;

end
close(hbar)

%% OLA

yout = ola(yt,fragsize,overlap);

%% .wav output

audiowrite([audioname,'_tc_out_',num2str(fs),'_',num2str(NFFT),'_',num2str(fragsize),'.wav'], yout, fs);

author: XiaoxiangGao


text:
function compareOutputs_subjs(group_h2, group_h3, type)

prec = -10;
if isempty(group_h2) || isempty(group_h3)
    return;
end
for qq=1:length(group_h3.group.subjs)
    yavg_h2 = eval( sprintf('group_h2.group.subjs(qq).procResult.%s;', type) );
    yavg_h3 = eval( sprintf('group_h3.group.subjs(qq).procStream.output.%s;', type) );
    for ll=1:size(yavg_h3,2)
        for mm=1:size(yavg_h3,3)
            for nn=1:size(yavg_h3,4)
                if ~isequal_prec(yavg_h2(:,ll,mm,nn), yavg_h3(:,ll,mm,nn), prec)
                    fprintf('Subj %d: %s data for Hb type %d, Ch %d, Cond %d does not match\n', qq, type, ll, mm, nn);
                end
            end
        end
    end
end

author: BUNPC


text:
function dt = niftiClass2DataType(matlabType)
% Transform the string matlabType into a numerical codes used by NIFTI-1
%
%  dt = niftiClass2DataType(matlabType)
%
% Here is the table showing the assocation between Matlab type and NIFTI-1
% data type.
%
%  String (Matlab) ----> NIFTI-1
% -------------------------------
%            uint8 ----> 2
%            int16 ----> 4
%            int32 ----> 8
%  single, float32 ----> 16
%        complex64 ----> 32
%  double, float64 ----> 64
%            RGB24 ----> 128
%             int8 ----> 256
%           unit16 ----> 512
%           uint32 ----> 768
%       complex128 ----> 1792
%
% Franco (c) Stanford Vista team, 2012

matlabType = mrvParamFormat(matlabType);
switch matlabType
    case {'uint8'},      dt = 2;
    case {'int16'},      dt = 4;
    case {'int32'},      dt = 8;
    case {'float32','single'}, dt = 16;
    case {'complex64'},  dt = 32';
    case {'float64','double'}, dt = 64;
    case {'rgb64'},      dt = 128;
    case {'int8'},       dt = 256;
    case {'uint16'},     dt = 512;
    case {'uint32'},     dt = 768;
    case {'complex128'}, dt = 1792;
    otherwise,
        error('Unknown string %s\n',matlabType);

end

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%%train network config
param.caffe_path=fullfile(fileparts(pwd),'matcaffe');
param.solver_netfile=fullfile(fileparts(pwd),'model','solver_PQAN_image.prototxt');
param.fintune_model=fullfile(fileparts(pwd),'pretrain_model','googlenet_bn.caffemodel');
param.test_net_file=fullfile(fileparts(pwd),'model','test_PQAN_image.prototxt');
param.test_batch_size=32;
param.result_save_file='result.txt';
param.save_model_file='ilIDS_VID_model';
param.save_model_name='PQAN_322_iter';
param.test_interval=500;
param.train_maxiter=15000;
param.data_path=fullfile(fileparts(pwd),'generate_data');
param.traindata_filename='train_data_iLID';
param.testdata_filename='test_data_iLID';
param.use_data_split_index=1;
param.data_split_num=1;
param.train_person_num=150;
param.use_gpu=1;
param.gpu_id=0;
%%
for split_index=param.use_data_split_index:param.data_split_num
    if ~exist(strcat(param.save_model_file,num2str(split_index)),'file')
        mkdir(strcat(param.save_model_file,num2str(split_index)));
    end
    load(strcat(param.data_path,param.traindata_filename,num2str(split_index),'/train_data.mat'));
    %% find caffe
    cur_path=pwd;
    cd(param.caffe_path);
    addpath(genpath(param.caffe_path))
    caffe.reset_all;
    caffe.init_log(fullfile(cur_path,'log'));
    addpath(param.caffe_path);
    if param.use_gpu
        caffe.set_mode_gpu;
    else
        caffe.set_mode_cpu;
    end
    cd(cur_path);
    caffe_solver=caffe.get_solver(param.solver_netfile,param.gpu_id);
    net=caffe.get_net(param.test_net_file,'test');
    caffe_solver.use_caffemodel(param.fintune_model);
    input_data_shape=caffe_solver.nets{1}.blobs('data').shape;
    batch_size=input_data_shape(4);
    batch_data=zeros(input_data_shape,'single');
    batch_label=zeros(1,1,1,batch_size,'single');
    batch_sim=ones(1,1,1,1,'single');
    %
    batch_score=ones(1,1,3,batch_size,'single');
    %
    train_x_axis=[];
    train_y_axis=[];
    iter=0;
    while iter<param.train_maxiter
        [batch_data,batch_label]=get_train_minibatch(param,batch_data,batch_label,batch_size,train_data_cam1,train_data_cam2,label_train_cam1,label_train_cam2);
        caffe_solver.nets{1}.blobs('data').set_data(batch_data);
        caffe_solver.nets{1}.blobs('label').set_data(batch_label);
        caffe_solver.nets{1}.blobs('sim').set_data(batch_sim);
        caffe_solver.step(1);
        iter = caffe_solver.iter;
        if mod(iter,10)==0
            soft_loss=caffe_solver.nets{1}.blobs('loss3/loss').get_data;
            acc=caffe_solver.nets{1}.blobs('loss3/top-1').get_data;
            loss=caffe_solver.nets{1}.blobs('triplet_loss').get_data;
            loss1=caffe_solver.nets{1}.blobs('loss1/loss').get_data;
            loss2=caffe_solver.nets{1}.blobs('loss2/loss').get_data;
            contrastive_loss=caffe_solver.nets{1}.blobs('contrastive_loss').get_data;
            train_x_axis=[train_x_axis,iter];
            train_y_axis=[train_y_axis,soft_loss];
            plot(train_x_axis,train_y_axis);
            drawnow;
            fprintf('iter= %d ,softmaxloss= %f ,tripletloss= %f,contrastive_loss=%f.loss1=%f  loss2=%f\n',iter,soft_loss,loss,contrastive_loss,loss1,loss2);
            if soft_loss+loss+contrastive_loss+loss1+loss2 > 80
                fprintf('stop');
            end
        end
        % test in the training stage
        if iter>10000 && mod(iter,param.test_interval)==0
            caffe_solver.nets{1}.save(strcat(param.save_model_file,num2str(split_index),'/',param.save_model_name,'_',num2str(iter),'.caffemodel'));
            net.copy_from(strcat(param.save_model_file,num2str(split_index),'/',param.save_model_name,'_',num2str(iter),'.caffemodel'));
            net.blobs('data').reshape([224 224 3 param.test_batch_size]);
            load(strcat(param.data_path,param.testdata_filename,num2str(split_index),'/test_data.mat'));
            test_script;
        end
        if iter>9900 && mod(iter,100)==0
            caffe_solver.nets{1}.save(strcat(param.save_model_file,num2str(split_index),'/',param.save_model_name,'_',num2str(iter),'.caffemodel'));
        end
    end
%     test_network;
end

author: liuyuisanai

repo: Quality-Aware-Network