text:
%-fanDTasia ToolBox------------------------------------------------------------------
% This Matlab script is part of the fanDTasia ToolBox: a Matlab library for Diffusion
% Weighted MRI (DW-MRI) Processing, Diffusion Tensor (DTI) Estimation, High-order
% Diffusion Tensor Analysis, Diffusion Kurtosis Imaging (DKI) Estimation,
% Tensor ODF estimation, Visualization and more.
%
% A Matlab Tutorial on DW-MRI can be found in:
% http://www.cise.ufl.edu/~abarmpou/lab/fanDTasia/tutorial.php
%
%-CITATION---------------------------------------------------------------------------
% If you use this software please cite the following papers on 4th-order tensors:
% 1) A. Barmpoutis et al. "Diffusion Kurtosis Imaging: Robust Estimation from DW-MRI
%    using Homogeneous Polynomials", In the Proceedings of ISBI, 2011, pp. 262-265.
% 2) A. Barmpoutis and B.C. Vemuri, "A Unified Framework for Estimating Diffusion
%    Tensors of any order with Symmetric Positive-Definite Constraints",
%    In the Proceedings of ISBI, 2010, pp. 1385-1388.
%
%-DESCRIPTION------------------------------------------------------------------------
% This demo script shows how to compute the Diffusion Kurtosis Coefficients from a given
% Diffusion-Weighted MRI dataset. The method guarantees that the estimated diffusivity as
% well as the Diffusion Tensor are positive semi-definite, using the method
% in Sec. 4.2 of the ISBI'11 paper. Here the given demo dataset consists of 5 voxels,
% 30 gradient directions x 2 b-values from the real brain dataset used in the ISBI'11 paper.
%
%-USE--------------------------------------------------------------------------------
% [D,W]=DEMO_DKI_Estimation_Method2_v2;
%
% D: is a vector with the computed Unique Coefficients of the 2nd-order Diffusion Tensor
% W: is a vector with the computed Unique Coefficients of the 4th-order Kurtosis Tensor
%
%-DISCLAIMER-------------------------------------------------------------------------
% You can use this source code for non commercial research and educational purposes
% only without licensing fees and is provided without guarantee or warrantee expressed
% or implied. You cannot repost this file without prior written permission from the
% authors. If you use this software please cite the following papers:
% 1) A. Barmpoutis et al. "Diffusion Kurtosis Imaging: Robust Estimation from DW-MRI
%    using Homogeneous Polynomials", In the Proceedings of ISBI, 2011.
% 2) A. Barmpoutis and B.C. Vemuri, "A Unified Framework for Estimating Diffusion
%    Tensors of any order with Symmetric Positive-Definite Constraints",
%    In the Proceedings of ISBI, 2010, pp. 1385-1388.
%
%-AUTHOR-----------------------------------------------------------------------------
% Angelos Barmpoutis, PhD
% Digital Worlds Institute
% University of Florida, Gainesville, FL 32611, USA
% angelbar at ufl dot edu
%------------------------------------------------------------------------------------
function [DKI_D,DKI_W]=DEMO_DKI_Estimation_Method2_v2

%%% DATA OPENING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%open data files
[S,B]=openFDT('real_data_5voxels.fdt');

%S0 signal, no diffusion weighting
S0=S(:,:,:,1);

%acquisition shell 1, 30 orientations
S_1real=S(:,:,:,[2:31]);
GradientOrientations_1=B([2:31],[1:3]);
BValue_1=B(2,4);

%acquisition shell 2, 30 orientations
S_2real=S(:,:,:,[32:61]);
GradientOrientations_2=B([32:61],[1:3]);
BValue_2=B(32,4);

%%% OPTIONAL: ADD RICIAN NOISE TO THE DATA FOR QUANTITATIVE EVALUATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%

stdv=input('Do you want to add Rician noise to the data for validation? \n If yes, then give the Std. Dev. (e.g: 0.1), otherwise type 0.\n STD.DEV.:');
S_1=S_1real;
for i=1:size(S_1real,4)
    S_1(:,:,:,i)=sqrt((S_1real(:,:,:,i)+stdv*S0.*randn(size(S_1real,1),size(S_1real,2),size(S_1real,3))).^2+(stdv*S0.*randn(size(S_1real,1),size(S_1real,2),size(S_1real,3))).^2);
end
S_2=S_2real;
for i=1:size(S_2real,4)
    S_2(:,:,:,i)=sqrt((S_2real(:,:,:,i)+stdv*S0.*randn(size(S_2real,1),size(S_2real,2),size(S_2real,3))).^2+(stdv*S0.*randn(size(S_2real,1),size(S_2real,2),size(S_2real,3))).^2);
end
%your data can have as many shells and bvalues you want

%%% INITIALIZATION - COMPUTING AUXILIARY DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Construct a constant set of polynomial coefficients C
C_order2=constructSetOf81Polynomials(2)'; %computes C from section 5.1 (ISBI'10)
C_order4=constructSetOf321Polynomials(4)'; %computes C from section 5.1 (ISBI'10)
A_1=D4toDKImatrix(BValue_1);A_1=[A_1(:,[1:6])*C_order2 A_1(:,6+[1:15])];%Computes the matrix A from Table 1 (ISBI'11)
A_2=D4toDKImatrix(BValue_2);A_2=[A_2(:,[1:6])*C_order2 A_2(:,6+[1:15])];%Computes the matrix A from Table 1 (ISBI'11)

%shell 1
G_1_order2=constructMatrixOfMonomials(GradientOrientations_1, 2); %computes G from section 5.1 (ISBI'10)
G_1_order4=constructMatrixOfMonomials(GradientOrientations_1, 4); %computes G from section 5.1 (ISBI'10)
Gbig_1=[-BValue_1*G_1_order2 BValue_1*BValue_1/6*G_1_order4];  %all the monomials for orders 2 and 4.

%shell 2
G_2_order2=constructMatrixOfMonomials(GradientOrientations_2, 2); %computes G from section 5.1 (ISBI'10)
G_2_order4=constructMatrixOfMonomials(GradientOrientations_2, 4); %computes G from section 5.1 (ISBI'10)
Gbig_2=[-BValue_2*G_2_order2 BValue_2*BValue_2/6*G_2_order4];  %all the monomials for orders 2 and 4.

%your data can have as many shells and bvalues you want
Gbig=[Gbig_1; Gbig_2]; %all the monomials for orders 2 and 4 and bvalues b1 and b2.

%%%%%% MAIN LOOP - METHOD: LINEAR FITTING - NO CONSTRAINTS %%%%%%%%%%%%%%%%%%%%%%%%%
for x=1:size(S,1)
    for y=1:size(S,2)
        for z=1:size(S,3)

            logS_1=log(squeeze(S_1(x,y,z,:))/S0(x,y,z));
            logS_2=log(squeeze(S_2(x,y,z,:))/S0(x,y,z));

            %The following 2 steps implement the method in Sec. 4.2 of the ISBI'11 paper.
            %Step 1: Compute a positive-definite 4th-order tensor for each b-value
            D4_1=C_order4*lsqnonneg(-G_1_order4*C_order4, logS_1);%computes a positive-definite tensor for b1
            D4_2=C_order4*lsqnonneg(-G_2_order4*C_order4, logS_2);%computes a positive-definite tensor for b2

            %Step 2: Compute DKI from the positive definite 4th-order tensors.
            x1=zeros(size(C_order2,2),1); %the initialization is the zero vector
            x2=zeros(15,1); %the initialization is the zero vector
            for i=1:100 %we perform a kind of gradient descent for a number of iterations
                x1=lsqnonneg([A_1(:,[1:size(C_order2,2)]);A_2(:,[1:size(C_order2,2)])],[D4_1;D4_2]-[A_1(:,[size(C_order2,2)+1:size(A_1,2)]);A_2(:,[size(C_order2,2)+1:size(A_1,2)])]*x2);
                x2_new=pinv([A_1(:,[size(C_order2,2)+1:size(A_1,2)]);A_2(:,[size(C_order2,2)+1:size(A_1,2)])]) * ([D4_1;D4_2]-[A_1(:,[1:size(C_order2,2)]);A_2(:,[1:size(C_order2,2)])]*x1);
                %sum(abs(x2-x2_new)) %I just put this here for convergence check
                x2=x2_new;
            end
            dki=[C_order2*x1;x2];

            %Optional Validation if user adds noise to the data
            if stdv>0
                logS_1=log(squeeze(S_1real(x,y,z,:))/S0(x,y,z));
                logS_2=log(squeeze(S_2real(x,y,z,:))/S0(x,y,z));
                err(:,x,y,z)=abs(Gbig*dki-[logS_1; logS_2]);
            end

            %Store the data
            DKI_D(:,x,y,z)=dki([1:6]); %The 6 unique coefficients of the diffusion tensor D
            %If you want you can put the result in the form of a 3x3 matrix
            %D=[dki(6) dki(5)/2 dki(4)/2
            %   dki(5)/2 dki(3) dki(2)/2
            %   dki(4)/2 dki(2)/2 dki(1)];

            DKI_W(:,x,y,z)=dki([7:21]); %The 15 unique coefficients of the kurtosis tensor W
            %You can see which coefficient is which you can use the function: printTensor(DKI_W(:,x,y,z),4)
            % or if you want to plot a tensor or a tensor field as spherical functions
            % you have to download the plotTensors.m function developed by Angelos Barmpoutis, Ph.D.
            % and then uncomment the following lines.
            %
            % plotTensors(DKI_D(:,1,1,1),1,[321 1]); %3D ellipsoidal plot of D
            % plotTensors(DKI_W(:,1,1,1),1,[321 1]); %3D ellipsoidal plot of W

            %Optional Calculation of Dapp and Kapp
            Dapp(:,x,y,z) = G_1_order2*dki(1:6);
            Kapp(:,x,y,z) = (G_1_order4*dki(7:21))./((G_1_order2*dki(1:6)).^2);

        end
    end
end

%%%%%% OPTIONAL: PRINT RESULTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
all_err=[];
all_Dapp=[];
all_Kapp=[];
counter=0;
meanD=zeros(6,1);
meanW=zeros(15,1);
for x=1:size(S,1)
    for y=1:size(S,2)
        for z=1:size(S,3)
            if stdv>0
               all_err=[all_err; err(:,x,y,z)];
            end
            all_Dapp=[all_Dapp; Dapp(:,x,y,z)];
            all_Kapp=[all_Kapp; Kapp(:,x,y,z)];
            meanD=meanD+DKI_D(:,x,y,z);
            meanW=meanW+DKI_W(:,x,y,z);
            counter=counter+1;
        end
    end
end
meanD=meanD/counter;
meanW=meanW/counter;

fprintf(1,'\n-----RESULTS:-----\n\n');
fprintf('Number of fitted voxels: %d\n',counter);
if stdv>0
    fprintf(1,'Fitting Error: %.4f (std. dev. %.4f)\n',mean(all_err),std(all_err));
end
fprintf(1,'Mean Dapp: %.4f (std. dev. %.4f)\n',mean(all_Dapp),std(all_Dapp));
fprintf(1,'Mean Kapp: %.4f (std. dev. %.4f)\n',mean(all_Kapp),std(all_Kapp));
fprintf(1,'\nMean Diffusion Tensor D:\n');
printTensor(meanD,2);
fprintf(1,'\nMean Kurtosis Tensor W:\n');
printTensor(meanW,4);

% If you want to plot a Diffusion or Kurtosis tensor as spherical functions
% you have to download the plotTensors.m function developed by Angelos Barmpoutis, Ph.D.
% and then uncomment the following lines.
%
% subplot(1,2,1)
% plotTensors(meanD,1,[321 1]);
% title('Mean Diffusion Tensor');
% subplot(1,2,2);
% plotTensors(meanW,1,[321 1]);
% title('Mean Kurtosis Tensor');

fprintf(1,'\nIf you use this software please cite the following papers on DKI and DTI estimation:\n');
fprintf(1,'1) A. Barmpoutis et al. "Diffusion Kurtosis Imaging: Robust Estimation from DW-MRI\n');
fprintf(1,'   using Homogeneous Polynomials", In the Proceedings of ISBI, 2011, pp. 262-265.\n');
fprintf(1,'2) A. Barmpoutis and B.C. Vemuri, "A Unified Framework for Estimating Diffusion Tensors\n');
fprintf(1,'   of any order with Symmetric Positive-Definite Constraints",\n');
fprintf(1,'   In the Proceedings of ISBI, 2010, pp. 1385-1388.\n');


text:
function [Hist] = clipHistogram(Hist,NrBins,ClipLimit,NrX,NrY)
%  This function performs clipping of the histogram and redistribution of bins.
%  The histogram is clipped and the number of excess pixels is counted. Afterwards
%  the excess pixels are equally redistributed across the whole histogram (providing
%  the bin count is smaller than the cliplimit).

for i = 1:NrX
    for j = 1:NrY
        %   Calculate the total number of excess pixels.
        NrExcess = 0;
        for nr = 1:NrBins
            excess=Hist(i,j,nr) - ClipLimit;
            if excess > 0
                NrExcess = NrExcess + excess;
            end
        end

        %  Clip histogram and redistribute excess pixels in each bin
        binIncr = NrExcess / NrBins;
        upper = ClipLimit - binIncr;
        for nr = 1:NrBins
            if Hist(i,j,nr) > ClipLimit
                Hist(i,j,nr) = ClipLimit;
            else
                if Hist(i,j,nr) > upper
                    NrExcess = NrExcess + upper - Hist(i,j,nr);
                    Hist(i,j,nr) = ClipLimit;
                else
                    NrExcess = NrExcess - binIncr;
                    Hist(i,j,nr) = Hist(i,j,nr) + binIncr;
                end
            end
        end

        if NrExcess > 0
            stepSize = max(1,fix(1+NrExcess/NrBins));
            for nr = 1:NrBins
                NrExcess = NrExcess - stepSize;
                Hist(i,j,nr) = Hist(i,j,nr) + stepSize;
                if NrExcess < 1
                    break;
                end
            end
        end

    end
end

author: Sable

repo: mcbench-benchmarks


text:
%minmax - Show min and max values for any number of variables.
%
%  USAGE
%
%    minmax(x,y,...)
%
%    x,y,...        variables
%

% Copyright (C) 2012 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

function minmax_FMAT(varargin)

if nargin < 1,
  error('Incorrect number of parameters (type ''help minmax'' for details).');
end

for i = 1:length(varargin),
	str = '';
	% Variable index (min)
	str = [str '   [' int2str(i) '-] '];
	% Variable values (min)
	s = min(varargin{i});
	for j = 1:length(s),
		str = [str num2str(s(j)) ' '];
	end
	str(end) = [];
	disp(str);
	str = '';
	% Variable index (max)
	str = [str '   [' int2str(i) '+] '];
	% Variable values (max)
	s = max(varargin{i});
	for j = 1:length(s),
		str = [str num2str(s(j)) ' '];
	end
	str(end) = [];
	disp(str);
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
report_this_filefun(mfilename('fullpath'));

[lat,lon] = meshgrat(tmap,tmapleg);
%gx = rex(1,:);
%gy = rey(:,1)';

% tmap = km2deg(tmap/1);
% [X , Y]  = meshgrid(gx,gy);

ren = interp2(XB,YB,ZB,lon,lat);

l = ren < -8.8; ren(l) = -8.8;
mi = min(min(ren));

l = isnan(ren);
ren(l) = mi-2;

figure_w_normalized_uicontrolunits('pos',[150 100 600 500])

hold on; axis off
axm1 = axesm('MapProjection','eqaconic','MapParallels',[],...
    'MapLatLimit',[s4_south s3_north],'MapLonLimit',[s2_west s1_east])

meshm(ren,tmapleg,size(tmap),tmap);

daspectm('m',3);
tightmap
view([0 90])
camlight; lighting phong
set(gca,'projection','perspective');

load worldlo
%h = displaym(POline); set(h(1),'color',[0.9 0.9 0.9],'Linewidth',1.7)
% h2 = displaym(PPpoint);
%   h = displaym(PPtext); trimcart(h);

pl = plotm(ms(:,1), ms(:,2),'ow','Linewidth',1.4);
set(pl,'LineWidth',1,'MarkerSize',3,...
    'MarkerFaceColor',[0.8 0.8 0.8 ],'MarkerEdgeColor',[ 0.8 0.8 0.8  ])

plmag = 4.
l = a.Date > t0+4/365 & a.Date < t0+34/365 & a.Magnitude >= plmag & a.Magnitude < 5;

pl = plotm(a(l,2),a(l,1),'ow','Markersize',9);
set(pl,'LineWidth',1,'MarkerSize',6,...
    'MarkerFaceColor',[1 1 1 ],'MarkerEdgeColor',[ 0 0 0  ])

l = a.Date > t0+4/365 & a.Date < t0+34/365 & a.Magnitude >= 5.0 ;
pl = plotm(a(l,2),a(l,1),'^w','Markersize',9);
set(pl,'LineWidth',1,'MarkerSize',10,...
    'MarkerFaceColor',[1 1 1 ],'MarkerEdgeColor',[ 0 0 0  ])

% pl = plotm(a.Latitude,a.Longitude,'+k');
%set(pl,'LineWidth',0.5,'MarkerSize',2,...
%   'MarkerFaceColor','k','MarkerEdgeColor','k')
pl = plotm(main(:,2),main(:,1),'hw');
set(pl,'LineWidth',1,'MarkerSize',20,...
    'MarkerFaceColor','w','MarkerEdgeColor','k')

zdatam(handlem('allline'),2000) % keep line on surface
%zdatam(handlem('alltext'),10000) % keep line on surface
caxis([0.0 0.22])
j = hsv;
%j = j(64:-1:1,:);
j = [ [ 0.9 0.9 0.9] ; j];

colormap(j); brighten(0.0);

axis off; set(gcf,'color','w')

setm(gca,'ffacecolor','w')
setm(gca,'fedgecolor','k','flinewidth',3);
setm(gca,'mlabellocation',0.25)
setm(gca,'meridianlabel','on')
setm(gca,'plabellocation',0.25)
setm(gca,'parallellabel','on')
setm(gca,'Fontcolor','k','Fontweight','normal','FontSize',10,'Labelunits','dm')

h5 = colorbar;
set(h5,'position',[0.71 0.44 0.017 0.45],'TickDir','out','Ycolor','k','Xcolor','k',...
    'Fontweight','bold','FontSize',14,'Ticklength',[0.02 0.08],'Yticklabel',[]);
set(gcf,'Inverthardcopy','off');

sc =   scaleruler('RulerStyle','lines','MajorTick',0:15:30,'Linewidth',1,'color','k',...
    'MajorTickLength',2,'XLoc',-0.003,'YLoc',0.5663,'Zloc',40000,'FontSize',9,'Fontweight','normal',...
    'MinorTick',0:5:10,'color',[0 0 0 ],'TickDir','down')

return

scaleruler

setm(handlem('scaleruler'),'XLoc',-0.0,'YLoc',0.566,'Zloc',4000)
setm(handlem('scaleruler'),'MajorTick',0:10:30,...
    'MinorTick',0:5:10,'TickDir','down',...
    'MajorTickLength',(3),...
    'MinorTickLength',(3),'Linewidth',8)
%setm(handlem('scaleruler'),'RulerStyle','ruler')
setm(handlem('scaleruler1'),'RulerStyle','lines')

author: CelsoReyes


text:
classdef Display_Abstract < handle

    properties (Access = protected)
        handle
        iterationArray
        valueArray
        style
    end

    properties (Access = private)
        figTitle
    end

    methods (Access = protected, Abstract)

        setChartType(obj)

    end

    methods (Access = public)

        function obj = Display_Abstract(figTitle)
            obj.figTitle = figTitle;
        end

        function show(obj,nRows,nCols,i,margins)
            obj.setStyle(nRows,nCols,i,margins);
            obj.setChartType();
            title(obj.figTitle);
            grid on
        end

        function updateParams(obj,it,value,monitor)
            obj.iterationArray(end+1) = it;
            if ~isempty(value)
                obj.valueArray(end+1,:) = value;
            end
            monitor.increaseRefreshIterator();
        end

        function refresh(obj)
            if ~isempty(obj.valueArray) && ~isempty(obj.iterationArray)
                set(obj.handle,'XData',obj.iterationArray,'YData',obj.valueArray);
                set(obj.style,'XLim',[0 obj.iterationArray(end)])
                drawnow
            end
        end

    end

    methods (Access = private)

        function setStyle(obj,nRows,nCols,i,margins)
            obj.style = subplot_tight(nRows,nCols,i,margins);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function A = randherm(n,cnd)
%RANDHERM     Random Hermitian matrix
%
%   res = randherm(n,cnd)
%
% cnd, if specified, is the approximate condition number of the generated matrix
%

% written  11/05/12     S.M. Rump
% modified 12/06/12     S.M. Rump  lower case, typo
%

  if nargin==1
    A = randn(n) + sqrt(-1)*randn(n);
  else
    s = exp( [ 0 rand(1,n-2) 1 ] * log(cnd) ) / cnd;
    A = orth(randn(n)+sqrt(-1)*randn(n));
    A = A' * diag(s.*sign(randn(1,n))) * A;
  end

  % make sure A is Hermitian
  A = A - sqrt(-1)*diag(imag(diag(A)));   % make sure diagonal is real
  A = tril(A)+tril(A,-1)';

repo: OpenMAS


text:
clear all;close all;clc
%%
M_dir = {'E:\Janelia2014\Fish1_16states_30frames';
    'E:\Janelia2014\Fish2_20140714_2_4_16states_10frames';
    'E:\Janelia2014\Fish3_20140715_1_1_16_states_10frames';
    'E:\Janelia2014\Fish4_20140721_1_8_16states_20frames';
    'E:\Janelia2014\Fish5_20140722_1_2_16states_30frames';
    'E:\Janelia2014\Fish6_20140722_1_1_3states_30,40frames';
    'E:\Janelia2014\Fish7_20140722_2_3_3states_30,50frames';
    'E:\Janelia2014\Fish8_20141222_2_2_7d_PT_3OMR_shock_lowcut';};

for i_fish = 1:6,
    datadir = M_dir{i_fish};

    %% load anatomy
    tiffname = fullfile(datadir,'ave.tif');
    info = imfinfo(tiffname,'tiff');
    nPlanes = length(info);
    s1 = info(1).Height;
    s2 = info(1).Width;
    ave_stack = zeros(s1,s2,nPlanes);
    for i=1:nPlanes,
        ave_stack(:,:,i) = imread(fullfile(datadir,'ave.tif'),i);
    end

    % x-z view
    im = squeeze(max(ave_stack,[],1));
    im=double(im);
    temp=sort(im(:),'descend');
    th1=temp(round(length(im(:))/100));
    th2=min(im(:));
    out=(im-th2)/(th1-th2);
    out(out>1)=1;
    out = flipud(out');
    anat_zx = repmat(out,[1 1 3]);

    %%
    temp = fullfile(datadir,['Fish' num2str(i_fish) '_direct_load.mat']);
    save(temp,'anat_zx','-append');

end
% %%
% for i_fish = 1:7,
%     datadir = M_dir{i_fish};
%
%     load(['CONST_F' num2str(i_fish) '.mat'],'CONST');
%     CONST2 = CONST;
%     CONST2.stim_full = CONST.photostate;
%     CONST2.CellRespAvr = CONST.CRAZ;
%
%     CONST = [];
%     names = {'ave_stack','anat_yx','anat_yz','anat_zx','CInfo','periods','shift','dshift',...
%         'CellResp','CellRespAvr','Fictive','FictiveAvr','stim_full','stimAvr',...
%         'tlists','stimrangenames'};
%     for i = 1:length(names), % use loop to save variables into fields of CONST
%         eval(['CONST.',names{i},'=', 'CONST2.', names{i},';']);
%     end
%
%     %%
%     save(['CONST_F' num2str(i_fish) '.mat'],'CONST');
%
% end

author: xiuyechen

repo: FishExplorer


text:
function flat = vol2flatCurROILevels(volume,flat);
%
% flat = vol2flatCurROILevels(volume,flat)
%
% Calls vol2flatROI with the currently selected ROI.  Called from
% callback in xformFlatMenu.
%
% djh, 8/98
%
% djh, 2/2001, replaced globals with local variables
% ras, 10/2004, version for flat level view

if ~volume.selectedROI
  myErrorDlg('Must have a selected ROI in the Volume window before it can be transformed to the Inplane window.');
end

volROI = volume.ROIs(volume.selectedROI);
flatROI = vol2flatROILevels(volROI,volume,flat);
flat = addROI(flat,flatROI,1);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function train_id_net_vgg16(varargin)
% -------------------------------------------------------------------------
% Part 4.1: prepare the data
% -------------------------------------------------------------------------

imdb = load('./dataset/Flickr30k-prepareurl_data.mat');
imdb = imdb.imdb;
load('./dataset/Flickr30k-prepare/dense_feature_word2.1.mat');
imdb.charcnn = wordcnn;
%imdb.charmean = mean(imdb.charcnn(:,:,:,imdb.images.set==1),4);
% -------------------------------------------------------------------------
% Part 4.2: initialize a CNN architecture
% -------------------------------------------------------------------------
net = resnet52_new_hope_word2_pool_vgg19();
net.conserveMemory = true;
im_mean = imdb.rgbMean;
net.meta.normalization.averageImage = im_mean;
%net.meta.normalization.charmean = imdb.charmean;
% -------------------------------------------------------------------------
% Part 4.3: train and evaluate the CNN

% -------------------------------------------------------------------------
opts.train.averageImage = net.meta.normalization.averageImage;
opts.train.batchSize = 32;
opts.train.continue = true;
opts.train.gpus = 2;
opts.train.prefetch = false ;
opts.train.nesterovUpdate = true ;
opts.train.expDir = './data/vgg19_batch32_new_hope_word2.1_pool_shift_both_drop0.75';
opts.train.derOutputs = {'objective_img',1,'objective_txt',1} ;
%opts.train.gamma = 0.9;
opts.train.momentum = 0.9;
%opts.train.constraint = 100;
opts.train.learningRate = [0.1*ones(1,240)] ;
opts.train.weightDecay = 0.0001;
opts.train.numEpochs = numel(opts.train.learningRate) ;
[opts, ~] = vl_argparse(opts.train, varargin) ;
% Call training function in MatConvNet
[net,info] = cnn_train_dag(net, imdb, @getBatch,opts) ;

% --------------------------------------------------------------------
function inputs = getBatch(imdb,batch,opts)
% --------------------------------------------------------------------
%-- img data
im_url = imdb.images.data(batch) ;
im = vl_imreadjpeg(im_url,'Pack','Resize',[224,224],'Flip',...
    'CropLocation','random','CropSize',[0.8,1],...
    'Interpolation', 'bicubic','NumThreads',24,... %'Brightness', double(0.1*imdb.rgbCovariance),...
    'SubtractAverage',imdb.rgbMean,...
    'CropAnisotropy',[3/4,4/3]);
oim = im{1}; %bsxfun(@minus,im{1},opts.averageImage);
label_img =  imdb.images.label(batch);

%-- txt data
batchsize = numel(batch);
label_txt =  imdb.images.label(batch);
txt_branch = bsxfun(@times,(batch-1),5) + randi(5,batchsize,1);
txt = single(imdb.charcnn(:,txt_branch));
% code for faster construct one-hot vector (remove loop)   but seems not
% faster.... so I give up it.
%{
txt_v = txt(:);
txt_v(:,2) = reshape(repmat(0:(batchsize-1),32,1),1,[])';
txt_v(:,3) = reshape(repmat(0:31,1,batchsize),1,[])';
txt_v = txt_v(txt_v(:,1)>0,:);
txt_v = txt_v(:,1)+ bsxfun(@times,txt_v(:,2),642368)... %20074*32
    +bsxfun(@times,txt_v(:,3),20074);
txtinput = gpuArray(zeros(1,20074,32,batchsize,'single'));
txtinput(txt_v) = 1;
txtinput = permute(txtinput,[1,3,2,4]);
%}
txtinput = zeros(1,32,20074,batchsize,'single');
for i=1:batchsize
    len = numel(find(txt(:,i)>0));
    location = randi(33-len);
    for j=1:len
        v = txt(j,i);
        txtinput(1,location,v,i)=1;
        location = location+1;
    end
end
txtinput = gpuArray(txtinput);
%}
%--
inputs = {'x0',gpuArray(oim),'data2',txtinput,'label_img',label_img,'label_txt',label_txt};

author: layumi

repo: Image-Text-Embedding


text:
 function ob = Gdsft(om, Nd, varargin)
%function ob = Gdsft(om, Nd, varargin)
%| Construct Gdsft object, which computes (nonuniform) FT samples
%| of signals with dimensions [(Nd)] exactly, e.g., for testing Gnufft
%| For faster computation, use Gnufft instead.
%| See Gdsft_test.m for example usage.
%|
%| in
%|	om	[M D]		frequency locations (radians / sample)
%|	Nd	[1 D]		signal dimensions
%|
%| options
%|	mask	[(Nd)]		logical support array
%|	n_shift	[1 D]		see nufft_init
%|	nthread	[]		# of processor threads
%|	class	''		fatrix2 (default) or Fatrix or 'exact'
%|	use_mex	0|1		use jf_mex ? (default: has_mex_jf)
%|
%| out
%|	ob	[M np]		object
%|
%| Copyright 2005-7-22, Jeff Fessler, University of Michigan

if nargin == 1 && streq(om, 'test'), Gdsft_test, return, end
if nargin < 2, ir_usage, end

% defaults
arg.om_t = double(om)'; % [D M], as required by dtft_mex
arg.Nd = Nd;
arg.mask = [];
arg.n_shift = '';
arg.nthread = 1;
arg.class = 'fatrix2';
arg.use_mex = has_mex_jf;
arg.show = false; % used in Gdsft_gram()

% options
arg = vararg_pair(arg, varargin);
arg.ndim = numel(arg.Nd);

if isempty(arg.n_shift)
	arg.n_shift = zeros(1, arg.ndim);
end

if size(arg.om_t,1) ~= arg.ndim
	error 'dimension mismatch'
end

switch arg.class
case 'Fatrix'
	ob = Gdsft_Fatrix(arg);

case 'fatrix2'
	ob = Gdsft_fatrix2(arg);

case 'exact'
	ob = Gdsft_fatrix2(arg);
	ob = full(ob); % trick

otherwise
	fail 'class'
end

% Gdsft_fatrix2()
function ob = Gdsft_fatrix2(arg)
if ~arg.use_mex
	forw = @(arg, x) dtft(x, arg.om_t', 'n_shift', arg.n_shift);
	back = @(arg, y) fatrix2_maskit(arg.mask, ...
		dtft_adj(y, arg.om_t', arg.Nd, arg.n_shift));

elseif any(arg.n_shift ~= 0)
	arg.phasor = exp(1i * (arg.om_t' * arg.n_shift(:))); % [M 1]
	forw = @(arg, x) arg.phasor .* ...
		cast(jf_mex('dtft,forward', arg.om_t, ...
			double(x), int32(arg.nthread)), class(x));
	back = @(arg, y) fatrix2_maskit(arg.mask, ...
		cast(jf_mex('dtft,adjoint', arg.om_t, ...
			complexify(double(conj(arg.phasor) .* y)), ...
			int32(arg.Nd), int32(arg.nthread)), class(y)));
else
	forw = @(arg, x) ...
		cast(jf_mex('dtft,forward', arg.om_t, ...
			double(x), int32(arg.nthread)), class(x));
	back = @(arg, y) fatrix2_maskit(arg.mask, ...
		cast(jf_mex('dtft,adjoint', arg.om_t, ...
			complexify(double(y)), ...
			int32(arg.Nd), int32(arg.nthread)), class(y)));
end
odim = size(arg.om_t,2);
ob = fatrix2('arg', arg, 'imask', arg.mask, ...
	'idim', arg.Nd, 'odim', odim, ...
	'gram', @Gdsft_gram, 'forw', forw, 'back', back);

% Gdsft_Fatrix()
function ob = Gdsft_Fatrix(arg)

arg = Gdsft_init_Fatrix(arg); % initialize

if isempty(arg.mask)
	arg.mask = true([arg.Nd 1]); % [(Nd)]
end
arg.np = sum(arg.mask(:));
arg.dim = [size(arg.om_t,2) arg.np]; % [M np]
ob = Fatrix(arg.dim, arg, ...
	'gram', @Gdsft_gram_fail,  ...
	'forw', @Gdsft_forw_Fatrix, 'back', @Gdsft_back_Fatrix);

% Gdsft_init_Fatrix()
function arg = Gdsft_init_Fatrix(arg)

if ~isempty(arg.n_shift)
% fix: - ?
	arg.phasor = exp(1i * (arg.om_t' * arg.n_shift(:))); % [M 1]
	arg.phasor = diag_sp(arg.phasor); % trick: to handle multiples
else
	arg.phasor = 1;
end

% Gdsft_forw_Fatrix(): y = A * x
% in
%	x	[np L] or [(Nd) L]
% out
%	y	[M L]
%
function y = Gdsft_forw_Fatrix(arg, x)

if size(x,1) == arg.np
	x = embed(x, arg.mask);	% [(Nd) (L)]
end

if arg.use_mex
	y = jf_mex('dtft,forward', arg.om_t, double(x), int32(arg.nthread));
else
	y = dtft(x, arg.om_t', 'n_shift', arg.n_shift);
end
y = arg.phasor * y;

% Gdsft_back_Fatrix(): x = A' * y
% in
%	y	[M L]
% out
%	x	[np L]
%
function x = Gdsft_back_Fatrix(arg, y)

y = arg.phasor' * y;
if isreal(y)
	y = complexify(y);
end
if arg.use_mex
	x = jf_mex('dtft,adjoint', arg.om_t, y, int32(arg.Nd), int32(arg.nthread));
else
	x = dtft_adj(y, arg.om_t', arg.Nd, arg.n_shift);
end
x = reshape(x, prod(arg.Nd), []);
x = x(arg.mask(:),:);

% Gdsft_gram()
function [T, reuse] = Gdsft_gram_fail(ob, W, reuse)
% T = dsft_gram(ob, W);
error 'not done'

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [img,cropVals] = mrAnatAutoCrop(img, borderPix, thresh)
%
% img = mrAnatAutoCrop(img, [borderPix=1], [thresh=0])
%
% 2008.02.27 RFD: wrote it.

if(~exist('borderPix','var')||isempty(borderPix))
    borderPix = 1;
end
if(~exist('thresh','var')||isempty(thresh))
    thresh = 0;
end

sz = size(img);
m = img>thresh;
if(numel(sz)>2)
    tmp = sum(m,3);
    x = find(sum(tmp,1)); x = [x(1) x(end)];
    y = find(sum(tmp,2)); y = [y(1) y(end)];
    tmp = squeeze(sum(m,1));
    z = find(sum(tmp,1)); z = [z(1) z(end)];
    clear tmp;
    pad = [-borderPix borderPix];
    x = x+pad; y = y+pad; z = z+pad;
    x(1) = max(1,x(1)); y(1) = max(1,y(1)); z(1) = max(1,z(1));
    x(2) = min(sz(2),x(2)); y(2) = min(sz(1),y(2)); z(2) = min(sz(3),z(2));
    img = img(y(1):y(2),x(1):x(2),z(1):z(2));
    cropVals = [x; y; z];
else
    x = find(sum(~isnan(img)&img~=0,1)); x = [x(1) x(end)];
    y = find(sum(~isnan(img)&img~=0,2)); y = [y(1) y(end)];
    pad = [-borderPix borderPix];
    x = x+pad; y = y+pad;
    x(1) = max(1,x(1)); y(1) = max(1,y(1));
    x(2) = min(sz(2),x(2)); y(2) = min(sz(1),y(2));
    img = img(y(1):y(2),x(1):x(2));
    cropVals = [x; y; 1,1];
end
return;


text:
function result = hlp_cpuload
% Get the relative CPU load

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

osbean = java.lang.management.ManagementFactory.getOperatingSystemMXBean();
result = osbean.getSystemLoadAverage() / osbean.getAvailableProcessors();

repo: REST


text:
%% housekeeping
clear
close all
clc

%% rise the model
m=rise('lwz_model','solve_linear',true);

%% get the parameters

[p,priors]=create_parameters(true);

m=set(m,'parameters',p);

%% bring in the data

data=create_data();

%% estimate the model

ms=estimate(m,'data',data,...
    'estim_priors',priors,...
    'kf_presample',3,...
    'kf_init_variance',10);

% 'optimizer','bee_gate','optimset',struct('MaxTime',30*60)

author: jmaih


text:
function [filt, B, A] = ft_preproc_bandpassfilter(dat, Fs, Fbp, order, type, dir, instabilityfix, df, wintype, dev, plotfiltresp, usefftfilt)

% FT_PREPROC_BANDPASSFILTER applies a band-pass filter to the data and thereby
% removes the spectral components in the data except for the ones in the
% specified frequency band.
%
% Use as
%   [filt] = ft_preproc_bandpassfilter(dat, Fs, Fbp, order, type, dir, instabilityfix, df, wintype, dev, plotfiltresp, usefftfilt)
% where
%   dat             data matrix (Nchans X Ntime)
%   Fs              sampling frequency in Hz
%   Fbp             frequency band, specified as [Fhp Flp] in Hz
%   order           optional filter order, default is 4 (but) or dependent on frequency band and data length (fir/firls)
%   type            optional filter type, can be
%                     'but'       Butterworth IIR filter (default)
%                     'firws'     FIR filter with windowed sinc
%                     'fir'       FIR filter using MATLAB fir1 function
%                     'firls'     FIR filter using MATLAB firls function (requires MATLAB Signal Processing Toolbox)
%                     'brickwall' frequency-domain filter using forward and inverse FFT
%   dir             optional filter direction, can be
%                     'onepass'                   forward filter only
%                     'onepass-reverse'           reverse filter only, i.e. backward in time
%                     'onepass-zerophase'         zero-phase forward filter with delay compensation (default for firws, linear-phase symmetric FIR only)
%                     'onepass-reverse-zerophase' zero-phase reverse filter with delay compensation
%                     'onepass-minphase'          minimum-phase converted forward filter (non-linear, only for firws)
%                     'twopass'                   zero-phase forward and reverse filter (default, except for firws)
%                     'twopass-reverse'           zero-phase reverse and forward filter
%                     'twopass-average'           average of the twopass and the twopass-reverse
%   instabilityfix  optional method to deal with filter instabilities
%                     'no'       only detect and give error (default)
%                     'reduce'   reduce the filter order
%                     'split'    split the filter in two lower-order filters, apply sequentially
%   df              optional transition width (only for firws)
%   wintype         optional window type (only for firws), can be
%                     'hamming' (default)    maximum passband deviation 0.0022 [0.22%], stopband attenuation -53dB
%                     'hann'                 maximum passband deviation 0.0063 [0.63%], stopband attenuation -44dB
%                     'blackman'             maximum passband deviation 0.0002 [0.02%], stopband attenuation -74dB
%                     'kaiser'
%   dev             optional max passband deviation/stopband attenuation (only for firws with kaiser window, default = 0.001 [0.1%, -60 dB])
%   plotfiltresp    optional, 'yes' or 'no', plot filter responses (only for firws, default = 'no')
%   usefftfilt      optional, 'yes' or 'no', use fftfilt instead of filter (only for firws, default = 'no')
%
% Note that a one- or two-pass filter has consequences for the strength of the
% filter, i.e. a two-pass filter with the same filter order will attenuate the signal
% twice as strong.
%
% Further note that the filter type 'brickwall' filters in the frequency domain,
% but may have severe issues. For instance, it has the implication that the time
% domain signal is periodic. Another issue pertains to that frequencies are
% not well defined over short time intervals; particularly for low frequencies.
%
% If the data contains NaNs, these will affect the output. With an IIR
% filter, and/or with FFT-filtering, local NaNs will spread to the whole
% time series. With a FIR filter, local NaNs will spread locally, depending
% on the filter order.
%
% See also PREPROC

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

% set the default filter order later
if nargin<4 || isempty(order)
  order = [];
end

% set the default filter type
if nargin<5 || isempty(type)
  type = 'but';
end

% set the default filter direction
if nargin<6 || isempty(dir)
  if strcmp(type, 'firws')
    dir = 'onepass-zerophase';
  else
    dir = 'twopass';
  end
end

% set the default method to deal with filter instabilities
if nargin<7|| isempty(instabilityfix)
  instabilityfix = 'no';
end

% Set default transition width later
if nargin < 8 || isempty(df)
  df = [];
end

% Set default window type
if nargin < 9 || isempty(wintype)
  wintype = 'hamming';
end

% Set default passband deviation/stopband attenuation for Kaiser window
if nargin < 10 || isempty(dev)
  if strcmp(wintype, 'kaiser')
    dev = 0.001;
  else
    dev = [];
  end
end

% Set default passband deviation/stopband attenuation for Kaiser window
if nargin < 11 || isempty(plotfiltresp)
  plotfiltresp = 'no';
end

% Set default filter function
if nargin < 12 || isempty(usefftfilt)
  usefftfilt = false;
else
  % convert to boolean value
  usefftfilt = istrue(usefftfilt);
end

% Filtering does not work on integer data
if ~isa(dat, 'double') && ~isa(dat, 'single')
  dat = cast(dat, 'double');
end

% preprocessing fails on channels that contain NaN
if any(isnan(dat(:)))
  ft_warning('FieldTrip:dataContainsNaN', 'data contains NaN values');
end

% Nyquist frequency
Fn = Fs/2;

% demean the data before filtering
meandat = nanmean(dat,2);
dat = bsxfun(@minus, dat, meandat);

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

  case 'firws'
    % Input arguments
    if length(Fbp) ~= 2
      ft_error('Two cutoff frequencies required.')
    end

    % Filter order AND transition width set?
    if ~isempty(order) && ~isempty(df)
      ft_warning('firws:dfOverridesN', 'Filter order AND transition width set - transition width setting will override filter order.')
    elseif isempty(order) && isempty(df) % Default transition width heuristic
      df = fir_df(Fbp, Fs);
    end

    % Compute filter order from transition width
    [foo, maxDf] = fir_df(Fbp, Fs); %#ok<ASGLU>
    isOrderLow = false;
    if ~isempty(df)
      if df > maxDf
        ft_error('Transition band too wide. Maximum transition width is %.2f Hz.', maxDf)
      end
      [order, dev] = firwsord(wintype, Fs, df, dev);
    else % Check filter order otherwise
      [df, dev] = invfirwsord(wintype, Fs, order, dev);
      if df > maxDf
        nOpt = firwsord(wintype, Fs, maxDf, dev);
        ft_warning('firws:filterOrderLow', 'Filter order too low. For better results a minimum filter order of %d is recommended. Effective cutoff frequency might deviate from requested cutoff frequency.', nOpt)
        isOrderLow = true;
      end
    end

    % Window
    if strcmp(wintype, 'kaiser')
      beta = kaiserbeta(dev);
      win = windows('kaiser', order + 1, beta);
    else
      win = windows(wintype, order + 1);
    end

    % Impulse response
    Fbp = sort(Fbp);
    B = firws(order, sort(Fbp) / Fn, 'band', win);
    A = 1;

    % Convert to minimum phase
    if strcmp(dir, 'onepass-minphase')
      B = minphaserceps(B);
    end

    % Twopass filtering
    if strncmp(dir, 'twopass', 7)
      pbDev = (dev + 1)^2 - 1;
      sbAtt = 20 * log10(dev^2);
      order = 2 * (length(B) - 1);
      isTwopass = true;
    else
      pbDev = dev;
      sbAtt = 20 * log10(dev);
      order = length(B) - 1;
      isTwopass = false;
    end

    % Reporting
    ft_info once
    ft_info('Bandpass filtering data: %s, order %d, %s-windowed sinc FIR\n', dir, order, wintype);
    if ~isTwopass && ~isOrderLow % Do not report shifted cutoffs
      ft_info('  cutoff (-6 dB) %g Hz and %g Hz\n', Fbp(1), Fbp(2));
      tb = [max([Fbp(1) - df / 2 0]), Fbp(1) + df / 2, Fbp(2) - df / 2, min([Fbp(2) + df / 2 Fn])]; % Transition band edges
      ft_info('  transition width %.1f Hz, stopband 0-%.1f Hz, passband %.1f-%.1f Hz, stopband %.1f-%.0f Hz\n', df, tb, Fn);
    end
    if ~isOrderLow
      ft_info('  maximum passband deviation %.4f (%.2f%%), stopband attenuation %.0f dB\n', pbDev, pbDev * 100, sbAtt);
    end

  case 'fir'
    if isempty(order)
      order = 3*fix(Fs / Fbp(1));
    end
    if order > floor( (size(dat,2) - 1) / 3)
      order=floor(size(dat,2)/3) - 1;
    end
    B = fir1(order, [min(Fbp)/Fn max(Fbp)/Fn]);
    A = 1;

  case 'firls' % from NUTMEG's implementation
    % Deprecated: see bug 2453
    ft_warning('The filter type you requested is not recommended for neural signals, only proceed if you know what you are doing.')
    if isempty(order)
      order = 3*fix(Fs / Fbp(1));
    end
    if order > floor( (size(dat,2) - 1) / 3)
      order=floor(size(dat,2)/3) - 1;
    end
    f = 0:0.001:1;
    if rem(length(f),2)~=0
      f(end)=[];
    end
    z = zeros(1,length(f));
    if(isfinite(min(Fbp)))
      [val,pos1] = min(abs(Fs*f/2 - min(Fbp)));
    else
      [val,pos2] = min(abs(Fs*f/2 - max(Fbp)));
      pos1=pos2;
    end
    if(isfinite(max(Fbp)))
      [val,pos2] = min(abs(Fs*f/2 - max(Fbp)));
    else
      pos2 = length(f);
    end
    z(pos1:pos2) = 1;
    A = 1;
    B = firls(order,f,z); % requires MATLAB signal processing toolbox

  case 'brickwall'
    ix = round((0:Fs) ./ (Fs ./ size(dat,2)));
    ax = ix ./ (size(dat,2)./Fs); % frequency axis, including the other end of the spectrum

    a    = ones(1, size(dat,2));
    fbin1 = nearest(ax, [min(Fbp)    max(Fbp)]);
    fbin2 = nearest(ax, [Fs-max(Fbp) Fs-min(Fbp)]); % same band at the other end of the spectrum

    a(1:(fbin1(1)-1))            = 0;
    a((fbin1(2)+1):(fbin2(1)-1)) = 0;
    a((fbin2(2)+1):end)          = 0;

    f    = fft(dat,[],2);             % FFT
    f    = f.*a(ones(size(dat,1),1),:); % brickwall
    filt = real(ifft(f,[],2));        % iFFT

  otherwise
    ft_error('unsupported filter type "%s"', type);
end

% Plot filter responses
if strcmp(plotfiltresp, 'yes')
  plotfresp(B, A, [], Fs, dir)
end

if ~isequal(type, 'brickwall')
  try
    filt = filter_with_correction(B,A,dat,dir,usefftfilt);
  catch
    switch instabilityfix
      case 'no'
        rethrow(lasterror);
      case 'reduce'
        ft_warning('off','backtrace');
        ft_warning('instability detected - reducing the %dth order filter to an %dth order filter', order, order-1);
        ft_warning('on','backtrace');
        filt = ft_preproc_bandpassfilter(dat,Fs,Fbp,order-1,type,dir,instabilityfix,df,wintype,dev,plotfiltresp,usefftfilt);
      case 'split'
        N1 = ceil(order/2);
        N2 = floor(order/2);
        ft_warning('off','backtrace');
        ft_warning('instability detected - splitting the %dth order filter in a sequential %dth and a %dth order filter', order, N1, N2);
        ft_warning('on','backtrace');
        filt = ft_preproc_bandpassfilter(dat ,Fs,Fbp,N1,type,dir,instabilityfix,df,wintype,dev,plotfiltresp,usefftfilt);
        filt = ft_preproc_bandpassfilter(filt,Fs,Fbp,N2,type,dir,instabilityfix,df,wintype,dev,plotfiltresp,usefftfilt);
      otherwise
        ft_error('incorrect specification of instabilityfix');
    end % switch
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


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

global planC stateS;
indexS = planC{end};

    [xV1, yV1, zV1] = getScanXYZVals(planC{indexS.scan}(scanSetF));
    [xV2, yV2, zV2] = getScanXYZVals(planC{indexS.scan}(scanSetM));

    XYZRes1  = {length(xV1), length(yV1), length(zV1)};

    xLims1 = [min(xV1) max(xV1)];
    yLims1 = [max(yV1) min(yV1)];
    zLims1 = [min(zV1) max(zV1)];

    sliceXVals1 = linspace(xLims1(1), xLims1(2), XYZRes1{1});
    sliceYVals1 = linspace(yLims1(1), yLims1(2), XYZRes1{2});
    sliceZVals1 = linspace(zLims1(1), zLims1(2), XYZRes1{3});

    for i=1:length(sliceZVals1)
        [xM1, yM1, zM1] = meshgrid(sliceXVals1, sliceYVals1, sliceZVals1(i));

        %Apply transformation to the limits if necessary.
        if isfield(planC{indexS.scan}(scanSetM), 'transM')
            mat = [xM1(:) yM1(:) zM1(:) ones(prod(size(xM1)), 1)]';
            mat = inv(planC{indexS.scan}(scanSetM).transM) * mat;
            xM = mat(1,:);
            yM = mat(2,:);
            zM = mat(3,:);
        end

        %Find corners of scan data included in this slice.
        [minX, jnk] = findnearest(xV2, min(xM(:)));
        [jnk, maxX] = findnearest(xV2, max(xM(:)));
        [minY, jnk] = findnearest(yV2, min(yM(:)));
        [jnk, maxY] = findnearest(yV2, max(yM(:)));

        [minZ, jnk] = findnearest(zV2, min(zM(:)));
        [jnk, maxZ] = findnearest(zV2, max(zM(:)));

        %Prepare the x,y,z vector inputs for finterp3.
        xVec = [xV2(1) xV2(2)-xV2(1) xV2(end)];
        yVec = [yV2(1) yV2(2)-yV2(1) yV2(end)];
        zVec = zV2(minZ:maxZ);

        %Interpolate to get slice.
        im = finterp3(xM(:), yM(:), zM(:), planC{indexS.scan}(scanSetM).scanArray(maxY:minY, minX:maxX, minZ:maxZ), xVec, yVec, zVec,0);

        imgSize      = [length(sliceYVals1) length(sliceXVals1) 1];
        im = reshape(im, imgSize);
        im = squeeze(im);
        slc(:,:,i) = im;
    end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [node,elem,HB] = cubemesh(box,h)
%% CUBEMESH a uniform mesh of a cube
%
% [node,elem,HB] = cubemesh([x0,x1,y0,y1,z0,z1],h) enerates a uniform mesh
% of the cube [x0,x1]*[y0,y1]*[z0,z1] with mesh size h.
%
% Example
%
%  [node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],0.5);
%  showmesh3(node,elem);
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

x0 = box(1); x1 = box(2);
y0 = box(3); y1 = box(4);
z0 = box(5); z1 = box(6);
node = [x0,y0,z0; x1,y0,z0; x1,y1,z0; x0,y1,z0; ...
        x0,y0,z1; x1,y0,z1; x1,y1,z1; x0,y1,z1];
elem = [1 2 3 7; 1 4 3 7; 1 5 6 7; 1 5 8 7; 1 2 6 7; 1 4 8 7];
n = ceil(log2(abs(x1-x0)/h));
for k = 1:n
    [node,elem] = uniformrefine3(node,elem);
end
% Set this as an initial grid
N0 = size(node,1);
HB(1:N0,1:3) = repmat((1:N0)',1,3);
HB(1:N0,4) = 0;

repo: Aero-matlab


text:
function w = cplxdual3D(x, J, Faf, af)

% 3D Complex Dual-Tree Discrete Wavelet Transform
%
% USAGE:
%   w = cplxdual3D(x, J, Faf, af)
% INPUT:
%   x - 3D array
%   J - number of stages
%   Faf - first stage filters
%   af - filters for remaining stages
% OUPUT:
%   w{j}{m}{n}{p}{d} - wavelet coefficients
%       j = 1..J, m = 1..2, n = 1..2, p = 1..2, d = 1..7
%   w{J+1}{m}{n}{d} - lowpass coefficients
%       m = 1..2, n = 1..2, p = 1..2, d = 1..7
% EXAMPLE:
%   x = rand(64,64,64);
%   J = 3;
%   [Faf, Fsf] = FSfarras;
%   [af, sf] = dualfilt1;
%   w = cplxdual3D(x, J, Faf, af);
%   y = icplxdual3D(w, J, Fsf, sf);
%   err = x - y;
%   max(max(max(abs(err))))
%
% WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
% http://taco.poly.edu/WaveletSoftware/

% normalization
x = x/sqrt(8);

for m = 1:2
    for n = 1:2
        for p = 1:2
            [lo w{1}{m}{n}{p}] = afb3D(x, Faf{m}, Faf{n}, Faf{p});
            for j = 2:J
                [lo w{j}{m}{n}{p}] = afb3D(lo, af{m}, af{n}, af{p});
            end
            w{J+1}{m}{n}{p} = lo;
        end
    end
end

for j = 1:J
    for m = 1:7
        [w{j}{1}{1}{1}{m} w{j}{2}{2}{1}{m} w{j}{2}{1}{2}{m} w{j}{1}{2}{2}{m}] = ...
            pm4(w{j}{1}{1}{1}{m}, w{j}{2}{2}{1}{m}, w{j}{2}{1}{2}{m}, w{j}{1}{2}{2}{m});
         [w{j}{2}{2}{2}{m} w{j}{1}{1}{2}{m} w{j}{1}{2}{1}{m} w{j}{2}{1}{1}{m}] = ...
            pm4(w{j}{2}{2}{2}{m}, w{j}{1}{1}{2}{m}, w{j}{1}{2}{1}{m}, w{j}{2}{1}{1}{m});
    end
end

author: thomaskuestner


text:
hydro = struct();

hydro = readCAPYTAINE(hydro,'ellipsoid_full.nc');
hydro = radiationIRF(hydro,15,[],[],[],[]);
hydro = radiationIRFSS(hydro,[],[]);
hydro = excitationIRF(hydro,15,[],[],[],[]);
writeBEMIOH5(hydro)
plotBEMIO(hydro)

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function [ know, x ] = p13_sol ( n )

%*****************************************************************************80
%
%% P13_SOL returns the solution for problem 13.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the problem.  This value
%    is only needed for those problems with variable N.
%
%    Output, integer KNOW.
%    If KNOW is 0, then the solution is not known.
%    If KNOW is positive, then the solution is known, and is returned in X.
%
%    Output, real X(N), the solution, if known.
%
  know = 0;

  x = zeros ( n, 1 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function data=nDexample(cx,N,n,seedi)
rand('seed',seedi);
mu = rand(cx,n)*10;
x=[];
C=[];
for i=1:cx
    r=[];
    dum=rand(n,n)-0.5+diag(rand(n,1))+diag(ones(n,1))/3;
    R=dum*dum';
 %   R=eye(n)/3;
    for j=1:N
         r(j,:) = randn(1,size(R,1)) * R + mu(i,:);
    end
    C=[C;ones(N,1)*i];
    x=[x;r];
end
%plot(x(:,1),x(:,2),'.');
data=[x];


text:
function out_channel_ascii( BstFile, OutputFile, Format, isEEG, isHeadshape, isHeader, Factor, Transf)
% OUT_CHANNEL_ASCII: Exports a Brainstorm channel file in an ascii file.
%
% USAGE:  out_channel_ascii( BstFile, OutputFile, Format={X,Y,Z}, isEEG=1, isHeadshape=1, isHeader=0, Factor=.01, Transf=[]);
%
% INPUT:
%     - BstFile    : full path to Brainstorm file to export
%     - OutputFile : full path to output file
%     - Format     : Cell-array describing the columns in the ASCII file
%                    => 'X','Y','Z' : 3D position (float)
%                    => 'indice'    : indice (integer)
%                    => 'name'      : name (string)
%                    => otherwise   : do not use this column
%     - isEEG       : Writes the coordinates of the electrodes
%     - isHeadshape : Writes the coordinates of the headshape points
%     - isHeader    : Writes header (number of EEG points)
%     - Factor      : Factor to convert the positions values in meters.
%     - Transf      : 4x4 transformation matrix to apply to the 3D positions before saving
%                     or entire MRI structure for conversion to MNI space

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2012-2020

%% ===== PARSE INPUTS =====
if (nargin < 3) || isempty(Format)
    Format = {'X','Y','Z'};
end
if (nargin < 4) || isempty(isEEG)
    isEEG = 1;
end
if (nargin < 5) || isempty(isHeadshape)
    isHeadshape = 1;
end
if (nargin < 6) || isempty(isHeader)
    isHeader = 0;
end
if (nargin < 7) || isempty(Factor)
    Factor = .01;
end
if (nargin < 8) || isempty(Transf)
    Transf = [];
end

% Load brainstorm channel file
BstMat = in_bst_channel(BstFile);
% Get all the positions
Loc    = zeros(3,0);
Label  = {};
if isEEG && isfield(BstMat, 'Channel') && ~isempty(BstMat.Channel)
    for i = 1:length(BstMat.Channel)
        if ~isempty(BstMat.Channel(i).Loc) && ~all(BstMat.Channel(i).Loc(:) == 0)
            Loc(:,end+1) = BstMat.Channel(i).Loc(:,1);
            Label{end+1} = strrep(BstMat.Channel(i).Name, ' ', '_');
        end
    end
end
if isHeadshape && isfield(BstMat, 'HeadPoints') && ~isempty(BstMat.HeadPoints) && ~isempty(BstMat.HeadPoints.Loc)
    Loc   = [Loc, BstMat.HeadPoints.Loc];
    Label = cat(2, Label, BstMat.HeadPoints.Label);
end

% Apply transformation
if ~isempty(Transf)
    % MNI coordinates: the entire MRI is passed in input
    if isstruct(Transf)
        Loc = cs_convert(Transf, 'scs', 'mni', Loc')';
    % World coordinates
    else
        R = Transf(1:3,1:3);
        T = Transf(1:3,4);
        Loc = R * Loc + T * ones(1, size(Loc,2));
    end
end
% Apply factor
Loc = Loc ./ Factor;

% Open output file
fid = fopen(OutputFile, 'w');
if (fid < 0)
   error('Cannot open file');
end
% Write header
nLoc = size(Loc,2);
if isHeader
    fwrite(fid, sprintf('%d\n', nLoc), 'char');
end
% Write file: one line per location
for i = 1:nLoc
    % Write each format entry
    for iF = 1:length(Format)
        % Entry types
        switch lower(Format{iF})
            case 'x',       str = sprintf('%1.4f', Loc(1,i));
            case '-x',      str = sprintf('%1.4f', -Loc(1,i));
            case 'y',       str = sprintf('%1.4f', Loc(2,i));
            case '-y',      str = sprintf('%1.4f', -Loc(2,i));
            case 'z',       str = sprintf('%1.4f', Loc(3,i));
            case '-z',      str = sprintf('%1.4f', -Loc(3,i));
            case 'indice',  str = sprintf('%d', i);
            case 'name',    str = Label{i};
            otherwise,      str = ' ';
        end
        % Write value
        fwrite(fid, str);
        % Add separator (space)
        if (iF ~= length(Format))
            fwrite(fid, sprintf('\t'));
        % Terminate line
        else
            fwrite(fid, 10);
        end
    end
end
% Close file
fclose(fid);

author: brainstorm-tools


text:
function plot_pcolor(vResults,j)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% example: plot_pcolor(vResults,85)
%
%  vResults : result matrix
%  i : position in result matrix
%  j : no of simulation to plot
%
%
figure;
set(gca,'LineWidth',2);
set(gca,'FontSize',14);
set(get(gca,'XLabel'),'FontSize',16);
set(get(gca,'YLabel'),'FontSize',16);

string=sprintf('load mCatBkgr%03.0f.mat',j);eval(string);
mBkgr=mCatalog;
string=sprintf('load vMain%03.0f.mat',j);eval(string);
vMain=logical(vMain);
string=sprintf('load mCatETAS%03.0f.mat',j);eval(string);
mETAS=mCatalog;
string=sprintf('load vDeclus%03.0f.mat',j);eval(string);
vDeclus=logical(vDeclus);

vTmp = ones(length(vResults(1).vX) * length(vResults(1).vY), 1) * nan;

for i=1:3
vTmp(vResults(i).vUsedNodes) = vResults(i).mValueGrid(:,j);
mPlotValues = reshape(vTmp, length(vResults(i).vY), length(vResults(i).vX));

mX=repmat(vResults(i).vX,length(vResults(i).vY),1);
mY=repmat(vResults(i).vY',1,length(vResults(i).vX));
% plot pcolor
subplot(3,1,i);
pcolor(mX,mY,mPlotValues);
shading interp;
xlabel('Lon');ylabel('Lat');

set(gca,'CLim',[-4 4]);
if i==3
    colorbar('location','southoutside')
end

switch i
    case 1
        hold on;plot(mETAS(vMain,1),mETAS(vMain,2),'k.');
    case 2
        hold on;plot(mETAS(vDeclus,1),mETAS(vDeclus,2),'k.');
    case 3
        hold on;plot(mETAS(vDeclus,1),mETAS(vDeclus,2),'k.');
        hold on;plot(mETAS(~vDeclus,1),mETAS(~vDeclus,2),'rx');
end
end


text:
function varargout = lz4hcencode(varargin)
%
% output = lz4hcencode(input)
%    or
% [output, info] = lz4hcencode(input)
%
% Compress a string or a numerical array using LZ4HC-compression
%
% This function depends on the ZMat toolbox (http://github.com/fangq/zmat)
%
% authors:Qianqian Fang (q.fang <at> neu.edu)
%
% input:
%      input: the original data, can be a string, a numerical vector or array
%
% output:
%      output: the compressed byte stream stored in a uint8 vector
%      info: (optional) a struct storing the metadata of the input, see "help zmat" for details
%
% examples:
%      [bytes, info]=lz4hcencode(eye(10));
%      orig=lz4hcdecode(bytes,info);
%
% license:
%     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details
%
% -- this function is part of JSONLab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

if(nargin==0)
    error('you must provide at least 1 input');
end

if(exist('zmat','file')==2 || exist('zmat','file')==3)
    [varargout{1:nargout}]=zmat(varargin{1}, 1,'lz4hc',varargin{2:end});
    return;
else
    error('you must install ZMat toolbox to use this feature: http://github.com/fangq/zmat')
end


text:
function denoised = WNNM_WRAP(image_path, sigma)

    N_Img = double(rgb2gray(imread(image_path)));

    Par   = ParSet(sigma);
    Denoised_Image = WNNM_DeNoising(N_Img, Par);

    recombinedRGBImage = cat(3, Denoised_Image, Denoised_Image, Denoised_Image);

    denoised = uint8(recombinedRGBImage);

end

author: lbasek

repo: image-denoising-benchmark


text:
function [termMats,alreadyHomogenized]=homogenizeMultiDimPolys(termMats,makeProjTransform)
%%HOMOGENIZEMULTIDIMPOLYS Given a system of multivariate polynomials,
%              produce a homogenized system. The degree of a monomial term
%              c*x1^a1*x2^a2*...xn^an for a1...an>=0 is the sum of all
%              a1...an. Most multivariate polynomials will contain terms of
%              different powers. This function adds an extra variable so
%              that all of the terms have the same power. That is for each
%              polynomial f(x1,x2,...,xn), the homogenized polynomial is
%              x0^d*f(x1/x0,x2/x0,...,xn/x0), where d is the highest degree
%              of any monomial term in the polynomial.
%
%INPUTS: termMats Either a single term matrix of a cell array containing
%           multiple term matrices. Each term matrix is An (n+1)XnumTerms
%           matrix such that termMat(:,i)=[c,a1,a2,...,an] is a monomial
%           term where c is the value of of the monomial coefficient and
%           the monomial is x1^(a1-1)*x2^(a2-1)*x3^(a3-1)...xn^(an-1). The
%           number of variables and the maximum degree need not be the same
%           for each term matrix. Differing numbers of variables implies
%           that the missing coefficients all have zero exponents.
% makeProjTransform An optional  boolean variable indicating whether the
%           final output system should be a projective transformation of
%           the input system. In such an instance, a cell array must be
%           returned as termMats is expanded by one equation, which is a
%           linear function consisting of random complex numbers times each
%           of the variables. As noted in Section 8 of [1], for a system of
%           n simultaneous polynomial equations in n unknowns, the
%           homogenized system/ projective transform system has a number of
%           solutions exactly equal to the B�zout bound and has no
%           solutions at infinity. If omitted or an empty matrix is passed,
%           the default is false.
%
%OUTPUTS: termMats The homogenized equations. If a sinle term matrix is
%            passed, then this is a matrix. otherwise, this is a cell array
%            of term matrices like the input. All of the term matrices now
%            contain the same number of variables (with zeros for unused
%            variables) with one more variable than the equation with the
%            most variables on the input. This new variable is the
%            homogenizing variable.
%  alreadyHomogenized A boolean value indicating that the input system was
%            already homogenized. In such an instance, all exponents of
%            the added homogenization variable in the retunred system are
%            zero.
%
%Homogenization of multivariate polynomial systems is often used in
%multivariate polynomial solving algorithms to better handle "solutions at
%infinite", which are generally not solutions that one cares about. For
%example, an application of homogenizing multivariate polynomials is in
%Section 8 of [1].
%
%EXAMPLE:
%Consider four equations in four unknowns.
% termMats=cell(4,1);
% termMats{1}=string2Terms('x1*x3-2*x2*x3-6*x1-x2-4*x4+12');
% termMats{2}=string2Terms('x1*x3-3*x2*x3-7*x1+12');
% termMats{3}=string2Terms('x1*x4^3-2*x1*x2-6*x4-x2+12');
% termMats{4}=string2Terms('x1^2+x2^2+x3^2+x4^2-6');
% [termMats1,alreadyHomogenized]=homogenizeMultiDimPolys(termMats)
% %The new system will correspond to
% %'x1*x3-2*x2*x3-6*x1*x5-x2*x5-4*x4*x5+12*x5^2'
% %'x2*x3-3*x3*x4-7*x1*x5+12*x5^2'
% %'x1*x4^3-2*x1*x2*x5^2-6*x4*x5^3-1*x2*x5^3+12*x5^4'
% %'x1^2+x2^2+x3^2+x4^2-6*x5^2'
% %The projective transformed solution is
% [termMats2,alreadyHomogenized]=homogenizeMultiDimPolys(termMats,true)
% %It is the same as the above solution except an extra equation consisting
% %of
% % 'c1*x1+c2*x2+c3*x3+c4*x4+c5*x5'
% %has been added, where c1...c5 are random complex numbers (with both
% %components between 0 and 1).
%
%REFERENCES:
%[1] L. T. Watson, S. C. Billups, and A. P. Morgan, "Algorithm 652 HOMPACK:
%    A suit of codes for globally convergent homotopy algorithms," ACM
%    Transactions on Mathematical Software, vol. 13, no. 3, pp. 281-310,
%    Sep. 1987.
%
%March 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(makeProjTransform))
   makeProjTransform=false;
end

isCell=true;
%If only a single matrix is passed.
if(~isa(termMats,'cell'))
    isCell=false;
    termMats={termMats};%Put it in a cell array.
end

numPoly=length(termMats);

%Determine the maximum number of variables.
n=0;
for curPoly=1:numPoly
   n=max([n,size(termMats{curPoly},1)-1]);
end

alreadyHomogenized=false;
for curPoly=1:numPoly
    termMat=termMats{curPoly};

    %The degrees of the terms.
    dTerms=sum(termMat(2:end,:),1);

    %The total degree of the polynomial.
    d=max(dTerms);

    numTerms=size(termMat,2);

    nCur=size(termMats,1)-1;

    %Augment to make have the same number of variables as the equation with
    %the most variables, plus one extra for the homogenization.
    termMat=[termMat;zeros(n-nCur,numTerms)];

    for i=1:numTerms
        %The degree of the term.
        deg=dTerms(i);

        %The homogenization makes all terms of degree d
        diff=d-deg;
        termMat(end,i)=diff;
        alreadyHomogenized=alreadyHomogenized||(d-deg~=0);
    end

    termMats{curPoly}=termMat;
end
alreadyHomogenized=~alreadyHomogenized;

%If only a single equation was passed as a matrix, not a cell array, then
%return a matrix.
if(isCell==false&&makeProjTransform==false)
    termMats=termMats{1};
    return;
end

if(makeProjTransform)
    termMat=zeros(n+2,n+1);

    termMat(1,:)=rand(1,n+1)+1j*rand(1,n+1);
    termMat(2:end,:)=eye(n+1,n+1);

    termMats{end+1}=termMat;
end

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
classdef TestStructuredEdgeDetection
    %TestStructuredEdgeDetection

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','balloon.jpg');
    end

    methods (Static)
        function test_1
            img = imread(TestStructuredEdgeDetection.im);
            img = single(img) / 255.0;

            modelFilename = get_model_file();
            pDollar = cv.StructuredEdgeDetection(modelFilename);

            E = pDollar.detectEdges(img);
            validateattributes(E, {'single'}, ...
                {'size',[size(img,1) size(img,2)], '>=',0, '<=',1});

            O = pDollar.computeOrientation(E);
            validateattributes(O, {'single'}, ...
                {'size',[size(img,1) size(img,2)]});

            E_nms = pDollar.edgesNms(E, O);
            validateattributes(E_nms, {'single'}, ...
                {'size',[size(img,1) size(img,2)], '>=',0, '<=',1});
        end

        function test_custom_feat_extract
            % skip test if external M-file is not found on the path
            if ~exist('myRFFeatureGetter.m', 'file')
                error('mexopencv:testskip', 'undefined function');
            end

            img = imread(TestStructuredEdgeDetection.im);
            img = single(img) / 255.0;

            modelFilename = get_model_file();
            pDollar = cv.StructuredEdgeDetection(modelFilename, ...
                'myRFFeatureGetter');
            E = pDollar.detectEdges(img);
        end

        function test_get_features
            img = imread(TestStructuredEdgeDetection.im);
            img = single(img) / 255.0;

            opts = struct();
            opts.normRad = 4;
            opts.grdSmooth = 0;
            opts.shrink = 2;
            opts.nChns = 13;
            opts.nOrients = 4;

            features = cv.StructuredEdgeDetection.getFeatures(img, opts);
            validateattributes(features, {'numeric'}, {'real', 'ndims',3});
            sz = size(features);
            %assert(sz(1) * opts.shrink == size(img,1));
            %assert(sz(2) * opts.shrink == size(img,2));
            assert(sz(3) == opts.nChns);
        end
    end

end

function features = myRFFeatureGetter(src, opts)
    if false
        nsize = fix([size(src,1) size(src,2)] ./ opts.shrink);
        features = zeros([nsize opts.nChns], 'single');
        %TODO: ... compute features
    else
        % call opencv's implementation
        features = cv.StructuredEdgeDetection.getFeatures(src, opts);
    end
end

function fname = get_model_file()
    fname = fullfile(mexopencv.root(),'test','model.yml.gz');
    if exist(fname, 'file') ~= 2
        % download model from GitHub
        url = 'https://cdn.rawgit.com/opencv/opencv_extra/3.2.0/testdata/cv/ximgproc/model.yml.gz';
        urlwrite(url, fname);
    end
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function  [ fx ] = f_VBfree( x,P,u,in )
% IN:
% - x_t : two posterior moments of the mean and variance of u^(o)
% - P : volatility of the moments (2)
% - u_t : previous feedback
% - in : []

theta = exp(P);

o = u(2);
ia = 4*(1-u(1));

mu1 = x(1+ia);
s1 = x(2+ia);
mu2 = x(3+ia);
s2 = x(4+ia);

Eexpx2 = exp(-mu2);%+s2/2);
s1 = 1./(Eexpx2+1./(s1+theta(1)));
mu1 = mu1 + s1.*Eexpx2*(o-mu1);

Esquarederr = (o-mu1).^2 + s1;
s2 = 1./(exp(-mu2)*Esquarederr/2+1./(s2+theta(2)));
mu2 = mu2 - s2.*(1-exp(-mu2)*Esquarederr)/2;

fx = x;
fx(ia+[1:4]') = [mu1;s1;mu2;s2];

author: MBB-team

repo: VBA-toolbox


text:
function [ts_smooth] = circ_smoothTS(varargin)
% USAGE
% [ts_smooth] = circ_smoothTS(varargin)
%
% Given a timeseries input of circular data, this function smooths over the
% desired number of bins and returns a circularly smooth timeseries
%
% INPUTS
%
%   ts         a time series of circular data
%   nBins      integer number of bins you would like to smooth over
%   method     string argument that determines the smoothing method,
%              options are 'median', 'mean' [default: 'median']
%   exclude    vector of values in ts to exclude from smoothing (can be
%              used to exclude 0 values)
%
% OUTPUTS
%
%   ts_smooth  a time series vector of smoothed circular data
%
% HELP
%
% Written by David Tingley, 2017
% TODO error handling when nBins > 1/2 ts

p = inputParser;
addRequired(p,'ts',@isvector);
addRequired(p,'nBins',@isnumeric);
addParameter(p,'method','median',@isstr)
addParameter(p,'exclude',[],@isvector);

parse(p,varargin{:});

ts = p.Results.ts;
if size(ts,1) == 1;
    ts = ts';
end

nBins = p.Results.nBins;
method = p.Results.method;
exclude = p.Results.exclude;

if length(exclude) == length(ts)
   ts_smooth = ts;
   return
end
if nBins == 1
    ts_smooth = ts;
    return
end

if ~isempty(exclude)
    list = find(ts==exclude);
    ts(list) = nan;
end

exclude = find(isnan(ts));
keep = find(~isnan(ts));

f = find(diff(keep)<nBins);
ff = find(diff(keep)>=nBins);
if ~isempty(ff)
    ff(end+1) = length(keep);  % include last ts
end
if length(keep) == 0
    ts(isnan(ts))=0;
    ts_smooth = ts;
    return
end

ts_smooth = zeros(length(ts),1);

for i =1:length(ff)  % populate list with single spikes that occur sparsely
    if keep(ff(i))>ceil(nBins/2) & keep(ff(i))+ceil(nBins/2) < length(ts) % prevents negative indices from being added
        ts_smooth(keep(ff(i))-ceil(nBins/2):keep(ff(i))+ceil(nBins/2)) = ts(keep(ff(i)));
    elseif keep(ff(i))+ceil(nBins/2) < length(ts)
        ts_smooth(1:keep(ff(i))+ceil(nBins/2)) = ts(keep(ff(i)));
    elseif keep(ff(i))>ceil(nBins/2)
        ts_smooth(keep(ff(i))-ceil(nBins/2):end) = ts(keep(ff(i)));
    end
end

for i=1:length(f) % populate list with spikes that occur within smoothing window
    ind = (keep(f(i))-ceil(nBins/2):keep(f(i))+ceil(nBins/2));
    ind(ind<1) = [];
    ind(ind>length(ts)) = [];
    keep = [keep; ind'];
%     if keep(f(i)) > ceil(nBins/2) & keep(f(i)+1)+ceil(nBins/2) < length(ts) % prevents negative indices from being added
%         keep = [keep; [keep(f(i))-ceil(nBins/2):keep(f(i)+1)+ceil(nBins/2)]'];
%     elseif ~(keep(f(i)) > ceil(nBins/2)) &  keep(f(i)+1)+ceil(nBins/2) < length(ts)
%         keep = [keep; [ceil(nBins/2) + 1:keep(f(i)+1)+ceil(nBins/2)]'];
%     elseif keep(f(i)) > ceil(nBins/2) & ~(keep(f(i)+1)+ceil(nBins/2) < length(ts))
%         keep = [keep; [keep(f(i))-ceil(nBins/2):length(ts)- ceil(nBins/2)]'];
%     end
end

keep = sort(unique(keep));
% keep(keep>0)=[];

%% start smoothing
for ii = 1:length(keep)
    i = keep(ii);
    ind = (i-ceil(nBins/2):i+ceil(nBins/2));
    ind(ind<1) = [];
    ind(ind>length(ts)) = [];
    [loc] = ~ismember(ind,exclude);

%     if ~isempty(loc) & sum(loc) > 0

        if strcmp(method,'median')
                ts_smooth(i) = circ_median(ts(ind(loc)));
        elseif strcmp(method,'mean')
                ts_smooth(i) = circ_mean(ts(ind(loc)));
        elseif strcmp(method,'gaussian')
                g = gauss(length(ind),1)';
                ts_smooth(i) = circ_mean(ts(ind(loc)),g(loc));
%                 ts_smooth(i) = circ_mean(ts(ind(loc)));
        elseif strcmp(method,'interp')
                error('interp not implented yet...')
        else
                error('couldnt find smoothing method')
        end
%     else
%         ts_smooth(ind) = ts(i);
%     end

end

ts_smooth(isnan(ts_smooth)) = 0;
% there's a lot of indexing going on, so let's double check the returned
% time series didn't change length...
if length(ts_smooth) ~= length(ts)
   error('output TS is the wrong length!')
end

author: buzsakilab

repo: buzcode


text:
function DoUpdateBatch(Simuh)

global VSeq
global VObj
global VCtl
global VMag
global VCoi
global VMot
global VVar
global VSig

eval(['VSeq' num2str(Simuh.ScanSeriesInd) '=VSeq;']);
eval(['VObj' num2str(Simuh.ScanSeriesInd) '=VObj;']);
eval(['VCtl' num2str(Simuh.ScanSeriesInd) '=VCtl;']);
eval(['VMag' num2str(Simuh.ScanSeriesInd) '=VMag;']);
eval(['VCoi' num2str(Simuh.ScanSeriesInd) '=VCoi;']);
eval(['VMot' num2str(Simuh.ScanSeriesInd) '=VMot;']);
eval(['VVar' num2str(Simuh.ScanSeriesInd) '=VVar;']);
eval(['VSig' num2str(Simuh.ScanSeriesInd) '=VSig;']);

if ~exist([Simuh.MRiLabPath filesep 'Tmp' filesep 'BatchData.mat'],'file')
    Notes='temporary file used for batch process';
    save([Simuh.MRiLabPath filesep 'Tmp' filesep 'BatchData'], 'Notes');
end

save([Simuh.MRiLabPath filesep 'Tmp' filesep 'BatchData'], ...
    ['VSeq' num2str(Simuh.ScanSeriesInd)], ...
    ['VObj' num2str(Simuh.ScanSeriesInd)], ...
    ['VCtl' num2str(Simuh.ScanSeriesInd)], ...
    ['VMag' num2str(Simuh.ScanSeriesInd)], ...
    ['VCoi' num2str(Simuh.ScanSeriesInd)], ...
    ['VMot' num2str(Simuh.ScanSeriesInd)], ...
    ['VVar' num2str(Simuh.ScanSeriesInd)], ...
    ['VSig' num2str(Simuh.ScanSeriesInd)], ...
    '-append');

Simuh.BatchList{end+1,1}=Simuh.ScanSeriesInd;
Simuh.BatchList{end,2}=Simuh.SimName;
Simuh.BatchList{end,3}=Simuh.Engine;
Simuh.BatchList{end,4}='Dx';
Simuh.BatchListIdx(end+1)=Simuh.ScanSeriesInd;

guidata(Simuh.SimuPanel_figure,Simuh);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
% POP_FILEIODIR - import directory into EEGLAB using FileIO
%
% Usage:
%   >> OUTEEG = pop_fileiodir; % pop up window
%   >> OUTEEG = pop_fileiodir( folder );
%
% Inputs:
%   folder - [string] folder name
%
% Optional inputs:
%   'channels'   - [integer array] list of channel indices
%   'samples'    - [min max] sample point limits for importing data.
%   'trials'     - [min max] trial's limit for importing data.
%
% Outputs:
%   OUTEEG   - EEGLAB data structure
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, 2012-
%
% Note: FILEIO toolbox must be installed.

% Copyright (C) 2012 Arnaud Delorme, SCCN, INC, UCSD, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [EEG, command] = pop_fileiodir(folder, varargin);
EEG = [];
command = '';

if nargin < 1
	% ask user
	folder = uigetdir('*.*', 'Choose a directory -- pop_fileiodir()');
	if folder == 0 return; end
    drawnow;

    % open file to get infos
    % ----------------------
    disp('Reading data file header...');
    dat = ft_read_header(folder);
    uilist   = { { 'style' 'text' 'String' 'Channel list (default all):' } ...
                 { 'style' 'edit' 'string' '' } ...
                 { 'style' 'text' 'String' [ 'Data range (in sample points) (default all [1 ' int2str(dat.nSamples) '])' ] } ...
                 { 'style' 'edit' 'string' '' }  };
    geom = { [3 1] [3 1] };
    if dat.nTrials > 1
        uilist{end+1} = { 'style' 'text' 'String' [ 'Trial range (default all [1 ' int2str(dat.nTrials) '])' ] };
        uilist{end+1} = { 'style' 'edit' 'string' '' };
        geom = { geom{:} [3 1] };
    end
    result = inputgui( geom, uilist, 'pophelp(''pop_fileiodir'')', 'Load data using FILE-IO -- pop_fileiodir()');
    if length(result) == 0 return; end

    options = {};
    result = { result{:} '' };
    if ~isempty(result{1}), options = { options{:} 'channels' eval( [ '[' result{1} ']' ] ) }; end
    if ~isempty(result{2}), options = { options{:} 'samples'  eval( [ '[' result{2} ']' ] ) }; end
    if ~isempty(result{3}), options = { options{:} 'trials'   eval( [ '[' result{3} ']' ] ) }; end
else
    dat = ft_read_header(folder);
    options = varargin;
end

[EEG command] = pop_fileio(folder, options{:});

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
%*****************************************************************************************************
%NAME: esemilogy.m
%AUTHOR: Andri M. Gretarsson
%DATE: 01/21/98
%
%SYNTAX:	esemilogy(X,Y, 'colour')
%
%Note that 'colour' is NOT an optional argument.
%
%This function acts just like the built-in function 'semilogy' but plots error-bars. The error-bars
%are plotted in the colour given by 'colour'.  'X' and 'Y' are Nx2 matrixes, the first column
%representing the values of the coordinate (x or y), the second column representing the uncertainty
%('error') of those values.  'colour' is a one-letter string which must be one of the letters allowed
%in the built-in matlab function 'plot', to specify the plot colour.  Note that the function exits
%with "hold" set to "off".
%
%This function does not print points in addition to the error bars.  Where the error bars cross, is
%the coordinate point.  This means that if both error bars are exceedingly small complared to the
%coordinate values, the mark will be correspondingly small.  In such situations, it may be better to
%use 'semilogy' directly and specify that the error is smaller than the size of the mark.
%
%EXAMPLE:
%
%X=[1.0	0.2
%   2.0	0.2]
%Y=[1.0 	0.25
%   2.0	0.25]
%esemilogy(X,Y,'g')
%
%plots a green cross of width 0.2 and height 0.25 at coordinate (1.0,1.0), and a cross of width 0.2
%and height 0.25 at coordinate( (2.0,2.0), on a linear x, log y scale
%
%LAST MODIFIED:  01/21/98
%*****************************************************************************************************

function esemilogx=semilogxplot(x,y,colourstring)

    xvalue=x(:,1);												%For clarity
    xerror=x(:,2);
    yvalue=y(:,1);
    yerror=y(:,2);

    semilogy(xvalue-xerror,yvalue-yerror,'w-',xvalue+xerror,yvalue+yerror,'w-'); hold on;
    %Sets appropriate axes but otherwise invisible on a white background.

    for i=1:length(xvalue)
        semilogy([xvalue(i)-xerror(i) xvalue(i)+xerror(i)], [yvalue(i) yvalue(i)], colourstring);
        semilogy([xvalue(i) xvalue(i)], [yvalue(i)-yerror(i) yvalue(i)+yerror(i)], colourstring);
    end
    hold off;

author: CelsoReyes


text:
addpath matlab;
vl_setupnn;

author: jiangqy


text:
function [D, names, counts] = LMcreateObjectIndexField(D)
%
% [D, names] = LMcreateObjectIndexField(D);
%
% This function adds a new field to the LabelMe index that provides a
% numeric index into the object class.
%
% The new field is:
%    annotation.object.namendx
% and this is a pointer into 'names'

% Create list of objects:
[names, counts, imagendx, objectndx, objclass_ndx] = LMobjectnames(D, 'name');
% add new field
for k = 1:length(objclass_ndx)
    D(imagendx(k)).annotation.object(objectndx(k)).namendx = objclass_ndx(k);
end

author: CSAILVision

repo: LabelMeToolbox


text:
%% FUNCTION Least_Trace
%   Trace-Norm Regularized Learning with Least Squares Loss.
%
%% OBJECTIVE
%   argmin_W { sum_i^t (0.5 * norm (Y{i} - X{i}' * W(:, i))^2)
%            + rho1 \|W\|_*}
%   where \|W\|_* = sum(svd(W, 0)) is the trace norm
%
%% INPUT
%   X: {n * d} * t - input matrix
%   Y: {n * 1} * t - output matrix
%   rho1: trace norm regularization parameter
%
%% OUTPUT
%   W: model: d * t
%   funcVal: function value vector.
%
%% LICENSE
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%   Copyright (C) 2011 - 2012 Jiayu Zhou and Jieping Ye
%
%   You are suggested to first read the Manual.
%   For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
%   Last modified on June 3, 2012.
%
%% RELATED PAPERS
%
%   [1] Ji, S. and Ye, J. An Accelerated Gradient Method for Trace Norm Minimization, ICML 2009
%
%% RELATED FUNCTIONS
%   Logistic_Trace, init_opts

%% Code starts here
function [W, funcVal] = Least_Trace(X, Y, rho1, opts)

if nargin <3
    error('\n Inputs: X, Y, and rho1 should be specified!\n');
end
X = multi_transpose(X);

if nargin <4
    opts = [];
end

% initialize options.
opts=init_opts(opts);

if isfield(opts, 'rho_L2')
    rho_L2 = opts.rho_L2;
else
    rho_L2 = 0;
end

task_num  = length (X);
dimension = size(X{1}, 1);
funcVal = [];

% precomputation.
XY = cell(task_num, 1);
W0_prep = [];
for t_idx = 1: task_num
    XY{t_idx} = X{t_idx}*Y{t_idx};
    W0_prep = cat(2, W0_prep, XY{t_idx});
end

% initialize a starting point
if opts.init==2
    W0 = zeros(dimension, task_num);
elseif opts.init == 0
    W0 = W0_prep;
else
    if isfield(opts,'W0')
        W0=opts.W0;
        if (nnz(size(W0)-[dimension, task_num]))
            error('\n Check the input .W0');
        end
    else
        W0=W0_prep;
    end
end

bFlag=0; % this flag tests whether the gradient step only changes a little

Wz= W0;
Wz_old = W0;

t = 1;
t_old = 0;

iter = 0;
gamma = 1;
gamma_inc = 2;

while iter < opts.maxIter
    alpha = (t_old - 1) /t;

    Ws = (1 + alpha) * Wz - alpha * Wz_old;

    % compute function value and gradients of the search point
    gWs  = gradVal_eval(Ws);
    Fs   = funVal_eval  (Ws);

    while true
        [Wzp Wzp_tn] = trace_projection(Ws - gWs/gamma, 2 * rho1 / gamma);
        Fzp = funVal_eval  (Wzp);

        delta_Wzp = Wzp - Ws;
        r_sum = norm(delta_Wzp, 'fro')^2;
        %Fzp_gamma = Fs + trace(delta_Wzp' * gWs) + gamma/2 * norm(delta_Wzp, 'fro')^2;
        Fzp_gamma = Fs + sum(sum(delta_Wzp .* gWs)) + gamma/2 * norm(delta_Wzp, 'fro')^2;

        if (r_sum <=1e-20)
            bFlag=1; % this shows that, the gradient step makes little improvement
            break;
        end

        if (Fzp <= Fzp_gamma)
            break;
        else
            gamma = gamma * gamma_inc;
        end
    end

    Wz_old = Wz;
    Wz = Wzp;

    %funcVal = cat(1, funcVal, Fzp + rho1 * sum( svd(Wzp, 0) ));
    funcVal = cat(1, funcVal, Fzp + rho1 * Wzp_tn);

    if (bFlag)
        % fprintf('\n The program terminates as the gradient step changes the solution very small.');
        break;
    end

    % test stop condition.
    switch(opts.tFlag)
        case 0
            if iter>=2
                if (abs( funcVal(end) - funcVal(end-1) ) <= opts.tol)
                    break;
                end
            end
        case 1
            if iter>=2
                if (abs( funcVal(end) - funcVal(end-1) ) <=...
                        opts.tol* funcVal(end-1))
                    break;
                end
            end
        case 2
            if ( funcVal(end)<= opts.tol)
                break;
            end
        case 3
            if iter>=opts.maxIter
                break;
            end
    end

    iter = iter + 1;
    t_old = t;
    t = 0.5 * (1 + (1+ 4 * t^2)^0.5);

end

W = Wzp;

% private functions
    function [grad_W] = gradVal_eval(W)
        if opts.pFlag
            grad_W = zeros(size(W));
            parfor t_ii = 1:task_num
                XWi = X{t_ii}' * W(:,t_ii);
                XTXWi = X{t_ii}* XWi;
                grad_W(:, t_ii) = XTXWi - XY{t_ii};
                %grad_W = cat(2, grad_W, X{t_ii}*(X{t_ii}' * W(:,t_ii)-Y{t_ii}) );
            end
        else
            grad_W = [];
            for t_ii = 1:task_num
                XWi = X{t_ii}' * W(:,t_ii);
                XTXWi = X{t_ii}* XWi;
                grad_W = cat(2, grad_W, XTXWi - XY{t_ii});
                %grad_W = cat(2, grad_W, X{t_ii}*(X{t_ii}' * W(:,t_ii)-Y{t_ii}) );
            end
        end
        grad_W = grad_W + rho_L2 * 2 * W;
    end

    function [funcVal] = funVal_eval (W)
        funcVal = 0;
        if opts.pFlag
            parfor i = 1: task_num
                funcVal = funcVal + 0.5 * norm (Y{i} - X{i}' * W(:, i))^2;
            end
        else
            for i = 1: task_num
                funcVal = funcVal + 0.5 * norm (Y{i} - X{i}' * W(:, i))^2;
            end
        end
        funcVal = funcVal + rho_L2 * norm(W, 'fro')^2;
    end

end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
% pop_rmdat() - Remove continuous data around specific events
%
% Usage:
%   >> OUTEEG = pop_rmdat( EEG); % pop-up a data entry window
%   >> OUTEEG = pop_rmdat( EEG, typerange, timelimits, invertselection);
%
% Graphic interface:
%   "Time-locking event type(s)" - [edit box] Select 'Edit > Event values'
%                to see a list of event.type values; else use the push button.
%                To use event types containing spaces, enter in single-quotes.
%   "..."      - [push button] scroll event types.
%   "Time limits" - [edit box] epoch latency range [start, end] in seconds relative
%                to the event type latency.
%
% Inputs:
%   EEG        - Input dataset. Data may already be epoched; in this case,
%                extract (shorter) subepochs time locked to epoch events.
%   typerange  - Cell array of event types to time lock to.
%                (Note: An event field called 'type' must be defined in
%                the 'EEG.event' structure).
%   timelimits - Epoch latency limits [start end] in seconds relative to the time-locking event
%                {default: [-1 2]}%
%   invertselection - [0|1] Invert selection {default:0 is no}
%
% Outputs:
%   OUTEEG     - output dataset
%
% Authors: Arnaud Delorme, SCCN, INC, UCSD, 2009-

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

function [EEG, com] = pop_rmdat( EEG, events, timelims, invertsel );

if nargin < 1
   help pop_rmdat;
	return;EEG = pop_rmdat( EEG,  {  'animal'  'bug' 'response'     },[-1 1] ,0);

end;
com = '';

if isempty(EEG.event)
    error( [ 'No event. This function removes data' 10 'based on event latencies' ]);
end;
if isempty(EEG.trials)
    error( [ 'This function only works with continuous data' ]);
end;
if ~isfield(EEG.event, 'latency'),
    error( 'Absent latency field in event array/structure: must name one of the fields ''latency''');
end;

if nargin < 3
	% popup window parameters
	% -----------------------
   promptstr    = { strvcat('Time-locking event type(s) ([]=all):', ...
                    'Select ''Edit > Event values'' to see type values.'), ...
                    'Epoch limits [start, end] in seconds:', ...
                    'Name for the new dataset:', ...
					'Out-of-bounds EEG rejection limits ([min max], []=none):'  };

   cbevent = ['if ~isfield(EEG.event, ''type'')' ...
				   '   errordlg2(''No type field'');' ...
				   'else' ...
                   '   tmpevent = EEG.event;' ...
                   '   if isnumeric(EEG.event(1).type),' ...
				   '        [tmps,tmpstr] = pop_chansel(unique([ tmpevent.type ]));' ...
				   '   else,' ...
                   '        [tmps,tmpstr] = pop_chansel(unique({ tmpevent.type }));' ...
                   '   end;' ...
				   '   if ~isempty(tmps)' ...
				   '       set(findobj(''parent'', gcbf, ''tag'', ''events''), ''string'', tmpstr);' ...
				   '   end;' ...
				   'end;' ...
				   'clear tmps tmpv tmpstr tmpevent tmpfieldnames;' ];

   geometry = { [2 1 1.2] [2 1 1.2] };
   uilist = { { 'style' 'text'       'string' 'Event type(s) ([]=all)' } ...
              { 'style' 'edit'       'string' '' 'tag' 'events' } ...
              { 'style' 'pushbutton' 'string' '...' 'callback' cbevent } ...
              { 'style' 'text'       'string' 'Time limits [start, end] in sec.' } ...
              { 'style' 'edit'       'string' '-1 1' } ...
              { 'style' 'popupmenu'  'string' 'Keep selected|Remove selected' } };

   result = inputgui( geometry, uilist, 'pophelp(''pop_rmdat'')', 'Remove data portions around events - pop_rmdat()');
   if length(result) == 0 return; end;

   if strcmpi(result{1}, '[]'), result{1} = ''; end;
   if ~isempty(result{1})
       if strcmpi(result{1}(1),'''')   % If event type appears to be in single-quotes, use comma
                                       % and single-quote as delimiter between event types. toby 2.24.2006
                                       % fixed Arnaud May 2006
            events = eval( [ '{' result{1} '}' ] );
       else events = parsetxt( result{1});
       end;
   else events = {};
   end
   timelims = eval( [ '[' result{2} ']' ] );
   invertsel = result{3}-1;

end;

tmpevent = EEG.event;
alltypes = { tmpevent.type };

% compute event indices
% ---------------------
allinds = [];
for index = 1:length(events)
    inds = strmatch(events{index},alltypes, 'exact');
    allinds = [allinds(:); inds(:) ]';
end;
allinds = sort(allinds);
if isempty(allinds)
    disp('No event found');
    return;
end;

% compute time limits
% -------------------
array = [];
tmpevent = EEG.event;
bnd = strmatch('boundary', lower({tmpevent.type }));
bndlat = [ tmpevent(bnd).latency ];
for bind = 1:length(allinds)
    evtlat = EEG.event(allinds(bind)).latency;
    evtbeg = evtlat+EEG.srate*timelims(1);
    evtend = evtlat+EEG.srate*timelims(2);
    if any(bndlat > evtbeg & bndlat < evtend)
        % find the closer upper and lower boundaries
        bndlattmp = bndlat(bndlat > evtbeg & bndlat < evtend);
        diffbound = bndlattmp-evtlat;
        allneginds = find(diffbound < 0);
        allposinds = find(diffbound > 0);
        if ~isempty(allneginds), evtbeg = bndlattmp(allneginds(1)); end;
        if ~isempty(allposinds), evtend = bndlattmp(allposinds(1)); end;
        fprintf('Boundary found: time limits for event %d reduced from %3.2f to %3.2f\n', allinds(bind), ...
            (evtbeg-evtlat)/EEG.srate, (evtend-evtlat)/EEG.srate);
    end
    if ~isempty(array) && evtbeg < array(end)
        array(end) = evtend;
    else
        array = [ array; evtbeg  evtend];
    end;
end;
array

if ~isempty(array) && array(1) < 1, array(1) = 1; end;
if ~isempty(array) && array(end) > EEG.pnts, array(end) = EEG.pnts; end;
if isempty(array)
    disp('No event found');
    return;
end;

if invertsel
    EEG = pop_select(EEG, 'notime', (array-1)/EEG.srate);
else
    EEG = pop_select(EEG, 'time', (array-1)/EEG.srate);
end;

% generate output command
% -----------------------
com = sprintf('EEG = pop_rmdat( EEG, %s);', vararg2str( { events timelims invertsel } ));

author: PatternRecognition

repo: OpenBMI


text:
function calcmcg( config )
rmpath(genpath([pwd '/dependencies/structuredEdges/release']));
images = dir(config.imageLocation);
images = regexpi({images.name}, '.*jpg|.*jpeg|.*png|.*bmp', 'match');
images = [images{:}];

if(~exist([config.outputLocation '/mcg'], 'dir'))
        mkdir(config.outputLocation,'/mcg')
end

for i=1:length(images)
    imname = char(images(i));
    impath = fullfile(config.imageLocation, imname);
    im=imread(impath);

    fprintf('Calculating MCG proposals for %s\n', imname);
    proposals=calcmcgForIm(im,config.mcg);
    saveFile=[imname '.mat'];
    save([config.outputLocation '/mcg/' saveFile], 'proposals');
end
addpath(genpath([pwd '/dependencies/structuredEdges/release']));
end

author: Cloud-CV

repo: object-proposals


text:
function [data] = load_DATA_block( file_pre, output_path, ...
            Nrg_lines_blk, Nrg_cells, AGC_atten_dB, block, UseMATfiles )
%
% This program - reads /loads data for a block
%              - converts to floating point
%              - compansates for the receiver attenuation
% -------------------------------------------------------------------------
% file_pre      : File prefix (path + prefix + block number)
% output_path   : Path where the SAR data MAT files are stored
% Nrg_lines_blk : Number of range lines in a block
% Nrg_cells     : Number of range cells
% AGC_atten_dB  : Nominal attenuation (Rx_AGC_setting) for each line (dB)
% block         : block number
% UseMATfiles   : Data is stored in MAT files
% -------------------------------------------------------------------------
% Created:  Nov 01, 2004  by Kaan Ersahin
% Modified: Nov 29, 2004  by Kaan Ersahin
% Modified: Dec  4, 2004  by Ian Cumming
%           - added single precision option for MATLAB 7
%           - added UseMATfiles option
%           - moved AGC correction to this function
% -------------------------------------------------------------------------

vv        = version;   vers = str2num(vv(1));  % Find MATLAB version

%  Read 4-bit unsigned data as single (MATLAB 7) or double
%  The data is stored as [I Q I Q I Q...] 4-bit words

if UseMATfiles
    % fprintf('\nload_DATA_block:  from MAT file  CDdata%1.0f.mat\n', block )
    disp ' '
    mat_file = strcat( output_path, 'CDdata', num2str(block) )
    load( mat_file )   % The data is already decoded in this case

else   %  Use fread to get all the data and parameters
    % fprintf('\nload_DATA_block:  from binary files\n')

    file_in  = strcat( file_pre, '_data.dat' )
    fid2     = fopen( file_in, 'r' );
    if vers < 7
      data = fread( fid2, [2*Nrg_cells, Nrg_lines_blk], 'ubit4=>double' );
    else
      data = fread( fid2, [2*Nrg_cells, Nrg_lines_blk], 'ubit4=>single' );
    end
    data = data';  %  Arrange the range lines as rows of the matrix

    %  Compensate for packed data format --> convert to signed numbers
    data = 2*( data - 16*(data > 7)) + 1;

    %  Separate the I and Q channels and make a complex array
    data = complex( data(:, 1:2:2*Nrg_cells), data(:, 2:2:2*Nrg_cells) );
    fclose(fid2);
end

% Apply gain correction to the data.  The attenuation varies from 2 to 17
% dB for this CD, so the linear gain factor varies from 1.26 to 7.08.
% If you want to store the decoded data in one byte arrays, use an additional
%  factor of 1.5, so that the maximum abs value is less than 127.
fact = 1.5;
if vers < 7      %  For MATLAB version 6
    linear_gain_factor = double( fact * 10.^(AGC_atten_dB/20) );
    ones_array         = ones(1,Nrg_cells);
    data = (linear_gain_factor * ones_array) .* double(data);
else             %  For MATLAB version 7
    linear_gain_factor = single( fact * 10.^(AGC_atten_dB/20) );
    ones_array         = single(ones(1,Nrg_cells));
    data = (linear_gain_factor * ones_array) .* single(data);
end

author: denkywu

repo: SAR-Synthetic-Aperture-Radar


text:
% REORIENT_FACETS Reorient faces of a triangle mesh (V,F) so that the left-hand
% rule normal of each face (consistently) points outward.
%
% [FF,I] = reorient_facets(V,F)
% [FF,I] = reorient_facets(V,F,'ParameterName',ParameterValue, ...)
%
% Inputs:
%   V  #V by 3 list of vertex positions
%   F  #F by 3 list of triangle indicies into V
%   Optional:
%     'NumRays'  followed by total number of rays {#F*100}
%     'MinRays'  followed by minimum number of rays per patch/face {10}
%     'Facetwise'  followed by whether each facet should be considered
%       independently, could lead to inconsistent orientation of manifoldly
%       neighboring facets {false}
%     'UseParity'  Whether to use parity(?) {false}
% Outputs:
%   FF   #F by 3 list of reoriented facets
%   I  #F list of whether each face was flipped
%

author: alecjacobson

repo: gptoolbox


text:
function test_issue1991

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY ft_read_mri ft_hastoolbox ft_filetype

[ft_ver, ft_path] = ft_version;

file0 = fullfile(ft_path, 'template', 'anatomy', 'single_subj_T1_1mm.nii');
file1 = dccnpath('/home/common/matlab/fieldtrip/data/test/original/jnifti/samples/colin27/colin27_zlib.jnii');
file2 = dccnpath('/home/common/matlab/fieldtrip/data/test/original/jnifti/samples/colin27/colin27_zlib.bnii');

%%

mri0 = ft_read_mri(file0);
mri1 = ft_read_mri(file1);
mri2 = ft_read_mri(file2);

%%

% do not compare the anatomy in the samples to the original, as that is obviously different
% it appears that the jnifti/samples represent a segmentation rather than the original MRI
fn = {'dim', 'unit', 'transform'};

for i=1:length(fn)
  assert(isequal(mri0.(fn{i}), mri1.(fn{i})), ['difference in jnii field ' fn{i}]);
  assert(isequal(mri0.(fn{i}), mri2.(fn{i})), ['difference in bnii field ' fn{i}]);
end

%%

assert(isequal(mri1.anatomy, mri2.anatomy), 'difference in anatomy')

%%

cfg = [];
figure; ft_sourceplot(cfg, mri0);
figure; ft_sourceplot(cfg, mri1);
figure; ft_sourceplot(cfg, mri2);

author: fieldtrip

repo: fieldtrip


text:
function f = schwefelsfcn(x)
% Schwefel's function.

if strcmp(x,'init')
    f.Aineq = [] ;
    f.bineq = [] ;
    f.Aeq = [] ;
    f.beq = [] ;
    f.LB = -500*ones(1,2) ; f.UB = 500*ones(1,2) ;
    f.nonlcon = [] ;
    f.options.PopulationSize = 500 ;
    f.options.PopInitRange = [-500; 500] ;
    f.options.ConstrBoundary = 'absorb' ;
    f.options.KnownMin = 420.9687*ones(1,2) ;
else
    x = reshape(x,1,[]) ;
    f = sum(-x.*sin(sqrt(abs(x))),2) ;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot