text:
function [h,g,a,info]=wfilt_remez(L,K,B)
%WFILT_REMEZ Filters designed using Remez exchange algorithm
%   Usage: [h,g,a]=wfilt_remez(L,K,B)
%
%   Input parameters:
%         L     : Length of the filters.
%         K     : Degree of flatness (regularity) at $z=-1$.
%         B     : Normalized transition bandwidth.
%
%   `[h,g,a]=wfilt_remez(L,K,B)` calculates a set of wavelet filters.
%   Regularity, frequency selectivity, and length of the filters can be
%   controlled by *K*, *B* and *L* parameters respectivelly.
%
%   The filter desigh algorithm is based on a Remez algorithm and a
%   factorization of the complex cepstrum of the polynomial.
%
%   Examples:
%   ---------
%   :::
%
%     wfiltinfo('remez50:2:0.1');
%
%   References: rioul94remez

% Original copyright goes to:
% Copyright (C) 1994, 1995, 1996, by Universidad de Vigo
% Author: Jose Martin Garcia
% e-mail: Uvi_Wave@tsc.uvigo.es

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

complainif_notposint(L,'L',mfilename);
complainif_notposint(L,'K',mfilename);

if B>0.2
    error(['%s: Bandwidth of the transition band should not be',...
           ' bigger than 0.2.'],upper(mfilename));
end

poly=remezwav(L,K,B);
rh=fc_cceps(poly);

g{1} = flipud(rh(:));
g{2} = -(-1).^(1:length(rh)).'.*flipud(g{1});

% Default offset
d = [0,0];
  % Do a filter alignment according to "center of gravity"
  d(1) = -floor(sum((1:L)'.*abs(g{1}).^2)/sum(abs(g{1}).^2));
  d(2) = -floor(sum((1:L)'.*abs(g{2}).^2)/sum(abs(g{2}).^2));
  if rem(d(1)-d(2),2)==1
      % Shift d(2) just a bit
      d(2) = d(2) + 1;
  end

g = cellfun(@(gEl,dEl) struct('h',gEl,'offset',dEl),g,num2cell(d),...
            'UniformOutput',0);
h = g;

a= [2;2];
info.istight = 1;

function [p,r]=remezwav(L,K,B)

%REMEZWAV    P=REMEZWAV(L,K,B) gives impulse response of maximally
%	     frequency selective P(z), product filter of paraunitary
%	     filter bank solution H(z) of length L satisfying K flatness
%	     constraints (wavelet filter), with normalized transition
%	     bandwidth B (optional argument if K==L/2).
%
%	     [P,R]=REMEZWAV(L,K,B) also gives the roots of P(z) which can
%	     be used to determine H(z).
%
%	     See also: REMEZFLT, FC_CCEPS.
%
%	     References: O. Rioul and P. Duhamel, "A Remez Exchange Algorithm
%			 for Orthonormal Wavelets", IEEE Trans. Circuits and
%			 Systems - II: Analog and Digital Signal Processing,
%			 41(8), August 1994
%
%       Author: Olivier Rioul, Nov. 1, 1992 (taken from the
%		above reference)
%  Modified by: Jose Martin Garcia
%       e-mail: Uvi_Wave@tsc.uvigo.es
%--------------------------------------------------------

computeroots=(nargout>1);

%%%%%%%%%%%%%%%%%%%%%%%%%% STEP 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%
if rem(L,2), error('L must be even'); end
if rem(L/2-K,2), K=K+1; end
N=L/2-K;
%%%%%%%%%%%%%%%%%%%%%%%%%% STEP 2  %%%%%%%%%%%%%%%%%%%%%%%%%%
% Daubechies solution
% PK(z)=z^(-2K-1))+AK(z^2)
if K==0, AK=0;
else
   binom=pascal(2*K,1);
   AK=binom(2*K,1:K)./(2*K-1:-2:1);
   AK=[AK AK(K:-1:1)];
   AK=AK/sum(AK);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%% STEP 2' %%%%%%%%%%%%%%%%%%%%%%%%%%%
% Daubechies factor
% PK(z)=((1+z^(-1))/2)^2*K QK(z)
if computeroots && K>0
   QK=binom(2*K,1:K);
   QK=QK.*abs(QK);
   QK=cumsum(QK);
   QK=QK./abs(binom(2*K-1,1:K));
   QK=[QK QK(K-1:-1:1)];
   QK=QK/sum(QK)*2;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%% STEP 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%
% output Daubechies solution PK(z)
if K==L/2
   p=zeros(1,2*L-1);
   p(1:2:2*L-1)=AK; p(L)=1;
   if computeroots
      r=[roots(QK); -ones(L,1)];
   end
   return
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% STEP 4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Daubechies polinomial
% PK(x)=1+x*DK(x^2)
if K==0, DK=0;
else
   binom=pascal(K,1);
   binom=binom(K,:);
   DK=binom./(1:2:2*K-1);
   DK=fliplr(DK)/sum(DK);
end

wp=(1/2-B)*pi;  % cut-off frequency
gridens=16*(N+1);  % grid density
found=0;  % boolean for Remez loop

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% STEP I %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initial estimate of yk
a=min(4,K)/10;
yk=linspace(0,1-a,N+1);
yk=(yk.^2).*(3+a-(2+a)*yk);
yk=1-(1-yk)*(1-cos(wp)^2);
ykold=yk;

iter=0;
while 1  % REMEZ LOOP
iter=iter+1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% STEP II %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute delta
Wyk=sqrt(yk).*((1-yk).^K);
Dyk=(1-sqrt(yk).*polyval(DK,yk))./Wyk;
for k=1:N+1
   dy=yk-yk(k); dy(k)=[];
   dy=dy(1:N/2).*dy(N:-1:N/2+1);
   Lk(k)=prod(dy);
end
invW(1:2:N+1)=2./Wyk(1:2:N+1);
delta=sum(Dyk./Lk)/sum(invW./Lk);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% STEP III %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute R(y) on fine grid
Ryk=Dyk-delta.*invW; Ryk(N+1)=[];
Lk=(yk(1:N)-yk(N+1))./Lk(1:N);
y=linspace(cos(wp)^2,1-K*1e-7,gridens);
yy=ones(N,1)*y-yk(1:N)'*ones(1,gridens);
% yy contain y-yk on each line
ind=find(yy==0);  % avoid division by 0
if ~isempty(ind)
   yy(ind)=1e-30*ones(size(ind));
end
yy=1./yy;
Ry=((Ryk.*Lk)*yy)./(Lk*yy);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% STEP IV %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find next yk
Ey=1-delta-sqrt(y).*(polyval(DK,y)+((1-y).^K).*Ry);
k=find(abs(diff(sign(diff(Ey))))==2)+1;
% N extrema
if length(k)>N
% may happen if L and K are large
   k=k(1:N);
end
yk=[yk(1) y(k)];
% N+1 extrema including wp
if K==0, yk=[yk 1]; end
% extrema at y==1 added
if all(yk==ykold), break; end
ykold=yk;

end  % REMEZ LOOP

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  STEP A %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute impulse response
w=(0:2*N-2)*pi/(2*N-1);
y=cos(w).^2;
yy=ones(N,1)*y-yk(1:N)'*ones(1,2*N-1);
ind=find(yy==0);
if ~isempty(ind)
   yy(ind)=1e-30*ones(size(ind));
end
yy=1./yy;
Ry=((Ryk.*Lk)*yy)./(Lk*yy);
Ry(2:2:2*N-2)=-Ry(2:2:2*N-2);
r=Ry*cos(w'*(2*(0:N-1)+1));
% partial real IDFT done
r=r/(2*N-1);
r=[r r(N-1:-1:1)];
p1=[r 0]+[0 r];
pp=p1;  % save p1 for later use
for k=1:2*K
   p1=[p1 0]-[0 p1];
end
if rem(K,2), p1=-p1; end
p1=p1/2^(2*K+1);
p1(N+1:N+2*K)=p1(N+1:N+2*K)+AK;
% add Daubechies response:
p(1:2:2*L-1)=p1; p(L)=1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% STEP A' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute roots
if computeroots
   Q(1:2:2*length(pp)-1)=pp;
   for k=1:2*K
     Q=[Q 0]-[0 Q];
   end
   if rem(K,2), Q=-Q; end
   Q=Q/2;
   if K>0  % add Daubechies factor QK
      Q(2*N+1:L-1)=Q(2*N+1:L-1)+QK;
   else
      Q(L)=1;
   end
   r=[roots(Q); -ones(2*K,1)];
end

function  h=fc_cceps(poly,ro)

%FC_CCEPS    Performs a factorization using complex cepstrum.
%
%	     H = FC_CCEPS (POLY,RO) provides H that is the spectral
%	     factor of a FIR transfer function POLY(z) with non-negative
%	     frequency response. This methode let us obtain lowpass
%	     filters of a bank structure without finding the POLY zeros.
%	     The filter obtained is minimum phase (all zeros are inside
%	     unit circle).
%
%	     RO is a parameter used to move zeros out of unit circle.
%	     It is optional and the default value is RO=1.02.
%
%	     See also: INVCCEPS, MYCCEPS, REMEZWAV.
%
%	     References: P.P Vaidyanathan, "Multirate Systems and Filter
%			 Banks", pp. 849-857, Prentice-Hall, 1993

%--------------------------------------------------------
% Copyright (C) 1994, 1995, 1996, by Universidad de Vigo
%
%
% Uvi_Wave 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, or (at your option) any
% later version.
%
% Uvi_Wave 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 Uvi_Wave; see the file COPYING.  If not, write to the Free
% Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
%
%       Author: Jose Martin Garcia
%       e-mail: Uvi_Wave@tsc.uvigo.es
%--------------------------------------------------------

if nargin < 2
	ro=1.02;
end

L=4096;   % number points of fft.

N=(length(poly)-1)/2;

%% Moving zeros out of unit circle
roo=(ro).^[0:2*N];
g=poly./roo;

%% Calculate complex cepstrum of secuence g
ghat=mycceps(g,L);

%% Fold the anticausal part of ghat, add it to the causal part and divide by 2
gcausal=ghat(1 : L/2);
gaux1=ghat(L/2+1 : L);
gaux2=gaux1(L/2 :-1: 1);
gantic=[0 gaux2(1 : L/2-1)];

xhat=0.5*(gcausal+gantic);

%% Calculate cepstral inversion
h=invcceps(xhat,N+1);

%% Low-pass filter has energie sqrt(2)
h=h*sqrt(2)/sum(h);

function  x=invcceps(xhat,L)

%INVCCEPS    Complex cepstrum Inversion
%
%	     X= INVCCEPS (CX,L) recovers X from its complex cepstrum sequence
%	     CX. X has to be real, causal, and stable (X(z) has no zeros
%	     outside unit circle) and x(0)>0. L is the length of the
%	     recovered secuence.
%
%	     See also: MYCCEPS, FC_CCEPS, REMEZWAV.
%
%	     References: P.P Vaidyanathan, "Multirate Systems and Filter
%			 Banks", pp. 849-857, Prentice-Hall, 1993

%--------------------------------------------------------
% Copyright (C) 1994, 1995, 1996, by Universidad de Vigo
%
%
% Uvi_Wave 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, or (at your option) any
% later version.
%
% Uvi_Wave 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 Uvi_Wave; see the file COPYING.  If not, write to the Free
% Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
%
%       Author: Jose Martin Garcia
%       e-mail: Uvi_Wave@tsc.uvigo.es
%--------------------------------------------------------

x=zeros(1,L);

%% First point of x
x(1)=exp(xhat(1));

%% Recursion to obtain the other point of x
for muestra=1:L-1
   for k=1:muestra
	x(muestra+1)=x(muestra+1)+k/muestra*xhat(k+1)*x(muestra-k+1);
   end
end

function xhat=mycceps(x,L)

%MYCCEPS     Complex Cepstrum
%
%	     CX = MYCCEPS (X,L) calculates complex cepstrum of the
%	     real sequence X. L is the number of points of the fft
%	     used. L is optional and its default value is 1024 points.
%
%	     See also: FC_CEPS, INVCCEPS, REMEZWAV.

%--------------------------------------------------------
% Copyright (C) 1994, 1995, 1996, by Universidad de Vigo
%
%
% Uvi_Wave 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, or (at your option) any
% later version.
%
% Uvi_Wave 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 Uvi_Wave; see the file COPYING.  If not, write to the Free
% Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
%
%       Author: Jose Martin Garcia
%       e-mail: Uvi_Wave@tsc.uvigo.es
%--------------------------------------------------------

if nargin < 2
   L=1024;
end

H = fft(x,L);

%% H must not be zero
ind=find(abs(H)==0);
if length(ind) > 0
   H(ind)=H(ind)+1e-25;
end

logH = log(abs(H))+sqrt(-1)*rcunwrap(angle(H));

xhat = real(ifft(logH));

function y = rcunwrap(x)
%RCUNWRAP Phase unwrap utility used by CCEPS.
%	RCUNWRAP(X) unwraps the phase and removes phase corresponding
%	to integer lag.  See also: UNWRAP, CCEPS.

%	Author(s): L. Shure, 1988
%		   L. Shure and help from PL, 3-30-92, revised
%	Copyright (c) 1984-94 by The MathWorks, Inc.
%       $Revision: 1.4 $  $Date: 1994/01/25 17:59:42 $

n = max(size(x));
y = unwrap(x);
nh = fix((n+1)/2);
y(:) = y(:)' - pi*round(y(nh+1)/pi)*(0:(n-1))/nh;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%% Dirichlet Mixture

clear classes

ds = prtDataGenOldFaithful;

mix = prtBrvMixture('components', repmat(prtBrvMvn,5,1), 'vbVerboseText', true, 'vbVerbosePlot', true, 'vbConvergenceThreshold', 1e-11);

[mixLearned, training] = mix.vbBatch(ds);

%% DP Mixture
clear classes
ds = prtDataGenOldFaithful;

mix = prtBrvDpMixture('components',repmat(prtBrvMvn,15,1), 'vbVerboseText',true, 'vbVerbosePlot', true, 'vbConvergenceThreshold',1e-10);
[mixLearned, training] = mix.vbBatch(ds);

%% Online Dirichlet Mixture

clear classes
%ds = prtDataGenOldFaithful;
ds = prtDataGenBimodal(5000);

mix = prtBrvMixture('components',repmat(prtBrvMvn,25,1),'vbOnlineLearningRateFunctionHandle',@(t)(32 + t).^(-0.6),'vbVerbosePlot',true,'vbOnlineBatchSize',25,'vbOnlineFullDataSetSize',ds.nObservations);
mixLearned = mix.vbOnline(ds.X);

%% Online DP Mixture

clear classes
ds = prtDataGenOldFaithful;
%ds = prtDataGenBimodal(5000);

mix = prtBrvDpMixture('components',repmat(prtBrvMvn,50,1),'vbOnlineLearningRateFunctionHandle',@(t)(32 + t).^(-0.8),'vbVerbosePlot',true,'vbOnlineBatchSize',25,'vbOnlineFullDataSetSize',ds.nObservations);
mixLearned = mix.vbOnline(ds.X);

%%

%% MCMC Dirichlet
clear classes
%ds = prtDataGenOldFaithful;
ds = prtDataGenBimodal;

mix = prtBrvMixture('components', repmat(prtBrvMvn,5,1), 'mcmcVerboseText', false, 'mcmcVerbosePlot', 100, 'mcmcTotalIterations', 5000);

[mixLearned, training] = mix.mcmc(ds);

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
% This script launches all the testUnits.

d = dir(pwd);
isub = [d(:).isdir];
nameFolds = {d(isub).name}';
nameFolds(ismember(nameFolds,{'.','..'})) = [];

for i=1:length(nameFolds)
    cd([pwd filesep nameFolds{i}]);
    eval([nameFolds{i} '_test_suite']);
    cd('..');
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [origin,W] = db_get_origin_info(INorigin,dbname,varargin)

% ORIGIN = DB_GET_ORIGIN_INFO(EVENT,DATABASE_NAME) If EVENT is a vector
% of origin id numbers (orid) then the returned structure ORIGIN contains
% fields pulled from the origin table.
%
% If EVENT is a vector of string times (one time per row) then the database
% is searched for the origin occuring closest in time to each element of
% EVENT (preferred origins only). Using the ORID is a more precise way to
% pull events from the database, however, orid's have no inherent meaning
% and may change. Using (approximate) origin times avoids this.
%
% [ORIGIN,WAVEFORM] = DB_GET_ORIGIN_INFO(EVENT,DATABASE_NAME,WAVEFORM) Add
% origin information as new fields in WAVEFORM object.

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

%  FOR TESTING ONLY
%dbname = '/home/admin/databases/AU06/wf/au06wf_2006_01';
%INorigin = [
%	'1/26/2006 11:58:59'
%	'1/27/2006 04:10:13'
%	'1/11/2006 07:01:49'
%	'1/31/2006 22:02:10'
%];
%INorigin = [
%	31301880
%	31301918
%	31300878
%	31302068
%];

if isnumeric(INorigin)
	origin.orid = reshape( INorigin , numel(INorigin), 1 );
elseif ischar(INorigin)
	for n = 1:size(INorigin,1)
		tmp(n) = str2epoch(INorigin(n,:));
	end
	tmp = getorids(tmp,dbname);
	origin.orid = reshape( tmp , numel(tmp), 1 );

else
    error('First input must be either origin ids or sting-formatted times');
end

% Check for waveform
if numel(varargin) == 1
	W = varargin{1};
	if ~isa(W,'waveform')
		error('Third argument must be a WAVEFORM object');
	end
elseif numel(varargin) > 1
	error('Too many arguments');
end
if nargout==2 && nargin~=3
	error('Mismatched number of input and/or output arguments');
end

% GET ORIGIN TIME
db = dbopen(dbname,'r');
db = dblookup_table(db,'origin');
nlist = [];
for n = 1:length(origin.orid);
	db1 = dbsubset(db,['orid==' num2str(origin.orid(n))]);
	[ origin.lat(n),origin.lon(n),origin.depth(n),origin.time_epoch(n),origin.ml(n),origin.orid(n)] = dbgetv(db1,'lat','lon','depth','time','ml','orid');
end
dbclose(db);
origin.lat = reshape( origin.lat , numel(origin.lat), 1 );
origin.lon = reshape( origin.lon , numel(origin.lon), 1 );
origin.depth = reshape( origin.depth , numel(origin.depth), 1 );
origin.time_epoch = reshape( origin.time_epoch , numel(origin.time_epoch), 1 );
origin.ml = reshape( origin.ml , numel(origin.ml), 1 );
origin.time_matlab = datenum(strtime(origin.time_epoch));
origin.orid = reshape( origin.orid , numel(origin.orid), 1 );
if exist('W')
	for n = 1:numel(W)
		disp([ 'Retreiving origin info for orid: ' num2str(origin.orid(n)) ' ... '])
		W(n) = addfield(W(n),'ORIGIN_LAT',origin.lat(n));
		W(n) = addfield(W(n),'ORIGIN_LON',origin.lon(n));
		W(n) = addfield(W(n),'ORIGIN_DEPTH',origin.depth(n));
		W(n) = addfield(W(n),'ORIGIN_TIME_EPOCH',origin.time_epoch(n));
		W(n) = addfield(W(n),'ORIGIN_ML',origin.ml(n));
		W(n) = addfield(W(n),'ORIGIN_TIME_MATLAB',origin.time_matlab(n));
		W(n) = addfield(W(n),'ORIGIN_ORID',origin.orid(n));
	end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RETURN ORID FOR EVENT CLOSEST IN TIME TO INPUT
function neworigin = getorids(origin,dbname);

db = dbopen(dbname,'r');
db = dblookup_table(db,'origin');
db1 = dblookup_table(db,'event');
db = dbjoin(db,db1);
db = dbsubset(db,'orid==prefor');
nrecords = dbquery(db,'dbRECORD_COUNT');
%display(['Number of records: ' num2str(nrecords)]);
[orid,time] = dbgetv(db,'orid','time');
dbclose(db);

for n = 1:length(origin)
	[tmp,index] = min(abs(origin(n)-time));
	neworigin(n) = orid(index);
end;

author: geoscience-community-codes

repo: GISMO


text:
function [trans_head2mri] = fiff_read_coord_trans(transfile)
%
%  usage: [trans_head2mri] = fiff_read_coord_trans(transfile)
%
%  input:
%       transfile = name of transformation fif file (usually stored in
%       the subject's Freesurfer directory in /mri/T1-neuromag/sets).
%
%  output:
%       trans_head2mri = transformation structure from head to MRI coordinate systems.
%
%                    Note: the inverse transformation, from MRI to head coordinate systems
%                              can be obtained by just taking the inverse:
%                                   trans_mri2head.from=5; trans_mri2head.to=4;
%                                   trans_mri2head.trans=inv(trans_head2mri.trans);
%
%  author:  Rey Ramirez  email: rrramir@uw.edu

FIFF = fiff_define_constants;
[fid,~,dir] = fiff_open(transfile);
a = find([dir.kind] == FIFF.FIFF_COORD_TRANS);
pos = dir(a(1)).pos;
tag = fiff_read_tag(fid,pos);
trans_head2mri = tag.data;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function p = printRxnFormulaLinked(model, rxnAbbrList, printFlag, metNameFlag, flux, nCharBreak, commandUsingMet)
% print reaction formulae in which the metabolites are hyperlinked to run a
% specific command involving the metabolites when displayed in the Matlab
% command window. Called by `surfNet.m`.
%
% USAGE:
%    p = printRxnFormulaLinked(model, rxnAbbrList, printFlag, metNameFlag, flux, nCharBreak, commandUsingMet)
%
% INPUTS:
%    model:                COBRA model
%
% OPTIONAL INPUTS:
%    rxnAbbrList:          cell array of reactions to be printed (default all reactions)
%    printFlag:            true to print, otherwise return the strings only (default true)
%    metNameFlag:          true to display met names in *.metNames instead of *.mets (default false)
%    flux:                 the flux of the reaction. When non-empty, determine the arrow for the formula accordingly
%    nCharBreak:           max. number of character printed on each line
%    commandUsingMet:      string of Matlab command for the genes in the model, e.g., 'fprintf(''%s'')' where %s
%                          will be replaced by the metabolites (default [], not adding the field for hyperlinked grRules)
if nargin < 2 || isempty(rxnAbbrList)
    rxnAbbrList = model.rxns;
end
if nargin < 3 || isempty(printFlag)
    printFlag = true;
end
if nargin < 4 || isempty(metNameFlag)
    metNameFlag = false;
end
if nargin < 5
    flux = [];
end
if nargin < 6 || isempty(nCharBreak)
    nCharBreak = get(0, 'CommandWindowSize');
    nCharBreak = nCharBreak(1);
end
if nargin < 7 || isempty(commandUsingMet)
    commandUsingMet = 'fprintf(''%s\\n'')';
end

if ~ischar(rxnAbbrList) && ~(iscell(rxnAbbrList) && numel(rxnAbbrList) == 1)
    p = cell(numel(rxnAbbrList), 1);
    for j = 1:numel(rxnAbbrList)
        p(j) = printRxnFormulaLinked(model, rxnAbbrList{j}, printFlag, metNameFlag, flux, nCharBreak);
    end
    return
end

rxnID = findRxnIDs(model, rxnAbbrList);

% determine the arrow sign
arrow = '->';
if (isempty(flux) || flux == 0) && model.lb(rxnID) < 0 && model.ub(rxnID) > 0
    arrow = '<=>';
elseif (isempty(flux) || flux == 0) && model.lb(rxnID) < 0 && model.ub(rxnID) <= 0
    arrow = '<-';
elseif ~isempty(flux) && flux < 0
    arrow = '<-';
end

% reactant
metR = model.S(:, rxnID) < 0;
if ~any(metR)
    reactantSide = cell(2, 0);
    isMetR = false(2, 0);
    isMetRst = false(2, 0);
else% reactant stoichiometry
    metRst = cellfun(@num2str, num2cell(abs(model.S(metR, rxnID))), 'UniformOutput', false);
    % not printing the coefficient if it is 1
    metRst(model.S(metR, rxnID) == -1) = {''};
    reactantSide = [metRst(:)'; columnVector(model.mets(metR))'; repmat({'+'}, 1, sum(metR) - 1), {arrow}];
    isMetR = [false(1, sum(metR)); true(1, sum(metR)); false(1, sum(metR))];
    isMetRst = [true(1, sum(metR)); false(2, sum(metR))];
end

% product
metP = model.S(:, rxnID) > 0;
if ~any(metP)
    productSide = cell(2, 0);
    isMetP = false(2, 0);
    isMetPst = false(2, 0);
else
    % product stoichiometry
    metPst = cellfun(@num2str, num2cell(model.S(metP, rxnID)), 'UniformOutput', false);
    % not printing the coefficient if it is 1
    metPst(model.S(metP, rxnID) == 1) = {''};
    productSide = [metPst(:)'; columnVector(model.mets(metP))'; repmat({'+'}, 1, sum(metP) - 1), {''}];
    isMetP = [false(1, sum(metP)); true(1, sum(metP)); false(1, sum(metP))];
    isMetPst = [true(1, sum(metP)); false(2, sum(metP))];
end

formulaeParts = [reactantSide(:); productSide(:)];
isMet = [isMetR(:); isMetP(:)];
isStoich = [isMetRst(:); isMetPst(:)];
% remove empty cells
isMet(cellfun(@isempty, formulaeParts)) = [];
isStoich(cellfun(@isempty, formulaeParts)) = [];
formulaeParts(cellfun(@isempty, formulaeParts)) = [];
isArrow = strcmp(formulaeParts, arrow);

if metNameFlag
    productSide(2, :) = columnVector(model.metNames(metP))';
    reactantSide(2, :) = columnVector(model.metNames(metR))';
    formulaePartsToPrint = [reactantSide(:); productSide(:)];
    formulaePartsToPrint(cellfun(@isempty, formulaePartsToPrint)) = [];
else
    formulaePartsToPrint = formulaeParts;
end

%printFormulae(formulaeParts, dict, arrow, metNameFlag, nCharBreak, nCharBreakReal)

nChar = 0;  % number of characters printed in the current line
p = '';
[lineBreak, firstWordOnLine] = deal(false, true);

for jPart = 1:numel(formulaeParts)

    nChar = nChar + length(formulaePartsToPrint{jPart}) + 1;
    if ~isMet(jPart)  % not a metabolite name (stoich/'+'/'->'/'<=>')
        spaceAdded = 0;  % additional space for arrow
        if isArrow(jPart)
            spaceAdded = 4;
        end
        % check if a new line is needed
        if ~firstWordOnLine
            if ~isStoich(jPart) && jPart < numel(formulaeParts)
                % not a stoich ('+'/'->'/'<=>'): new line if togther with the next stoich + met exceed the limit
                if isStoich(jPart + 1) && nChar + length(formulaePartsToPrint{jPart + 1}) + length(formulaePartsToPrint{jPart + 2}) + 1 + spaceAdded >= nCharBreak
                    lineBreak = true;
                elseif isMet(jPart + 1) && nChar + length(formulaePartsToPrint{jPart + 1}) + spaceAdded >= nCharBreak
                    lineBreak = true;
                end
            end
        end
        firstWordOnLine = false;
        if lineBreak
            p = [p '\n  '];
            nChar = length(formulaePartsToPrint{jPart}) + 3;  % 2 spaces before + 1 space after
            [lineBreak, firstWordOnLine] = deal(false, true);
            if isArrow(jPart)
                spaceAdded = 2;
            end
        end
        [spaceAddLeft, spaceAddRight] = deal('');
        if isArrow(jPart)
            if jPart < numel(formulaeParts)
                spaceAddRight = '  ';
            end
            if ~firstWordOnLine
                spaceAddLeft = '  ';
            end
            nChar = nChar + length(spaceAddLeft) + length(spaceAddRight);
        end
        p = [p, spaceAddLeft, formulaeParts{jPart}, spaceAddRight, ' '];
    else
        % print hyperlink if it is a metabolite name
        p = [p, printHyperlink(strrep(commandUsingMet, '%s', formulaeParts{jPart}), formulaePartsToPrint{jPart}, 0, false), ' '];
        firstWordOnLine = false;
    end
end
if printFlag
    fprintf(p);
end
p = {p};
end

author: opencobra

repo: cobratoolbox


text:
function [lh_AR, varx] = ARMLfit(rcs,ng,xg,rc0,lag_max)

% The conditional-likelihood fit.
%
%   rcs = tan(.5*pi*rc)
%   In this way, the reflection coefficient [-1,+1] is mapped to
%   [-Inf,+Inf], allowing the use of unconstrained optimization.

%   S. de Waele, August 2001.

rc = 2/pi*atan(rcs);
if any(~isreal(rc)),lh_AR = +Inf; return; end
nseg = length(xg);
npre = length(rcs)+length(rc0);
a = rc2arset([1 rc0 rc]);
%cor = par2cor(a,lag_max+1);
cor = arma2cor(a,1,lag_max);

sum_s = 0;     %The squared error in the exponent of the normal distribution
     		   %without varx and the factor -.5;
sum_lnd = 0;   %contribution of the log of the determinant without varx and the factor .5;
nobs_used = 0;

for seg = 1:nseg,
    nc = ng{seg};
    if any(floor(nc)~=nc), error('Times must contain only integers.'), end
    xc = xg{seg};
    nobsc = length(nc);
    nobs_used = nobs_used + nobsc;

    obs_start = 1;
    if obs_start == 1,
        obs = 1;
        sum_s = sum_s+xc(1)^2;
        sum_lnd= sum_lnd+log(1);
    end
    for obs = max(2,obs_start):nobsc
        nprev_a = min(obs-1,lag_max); %nprev_a = available previous observations
        it=find(nc(obs)-nc(obs-1:-1:obs-nprev_a)<lag_max+1);
        nprev_a=length(it);
        if nprev_a
            n1 = nc(obs:-1:obs-nprev_a)*ones(1,nprev_a+1);
            n2 = ones(nprev_a+1,1)*nc(obs:-1:obs-nprev_a)';
            lags = abs(n1-n2);
            Rn = cor(1+lags);
            %Notations as in Random signals, K.S. Shanmugan and A.M. Breipohl, p.51.
            x1 = xc(obs);
            x2 = xc(obs-1:-1:obs-nprev_a);
            Rn11 = 1;
            Rn12 = Rn(1,2:end);
            Rn21 = Rn(2:end,1);
            Rn22 = Rn(2:end,2:end);
            muc = Rn12*(Rn22\x2);       % = Conditional mean of x1 given x2
            Rnc = 1-Rn12*(Rn22\Rn21);   % = scaled Conditional variance of x1.
            sum_s = sum_s + (x1-muc)^2/Rnc;
            sum_lnd = sum_lnd + log(Rnc);
        else
            sum_s = sum_s+xc(obs)^2;
            sum_lnd= sum_lnd+log(1);
        end
    end  %for obs =
end      %for seg = 1:nseg

varx = sum_s/nobs_used;

lh_AR = -nobs_used/2*log(varx)-sum_lnd/2 - sum_s/2/varx;

%Turning the Likelihood in an estimate of the K-L discrepancy:
lh_AR = -2*lh_AR;

if isnan(lh_AR), lh_AR = +Inf; end

author: Sable

repo: mcbench-benchmarks


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: SSD versus translations, splineInter, HNSP level=8
%
%==============================================================================

clear, close all, help(mfilename)

setup2DHNSPData;
level = 8; m = ML{level}.m;
imgModel('set','imgModel','splineInter');
[T,R] = imgModel('coefficients',ML{level}.T,ML{level}.R,omega);
xc    = getCellCenteredGrid(omega,m);
Rc    = imgModel(R,omega,xc);

trafo('set','trafo','translation2D');
figure(1); clf;
[w1,w2] = ndgrid(0.2*linspace(-1,1,21),0.2*linspace(-1,1,21));
dc  = zeros(size(w1));
for j=1:numel(dc),
  yc = trafo([w1(j);w2(j)],xc);
  Tc = imgModel(T,omega,yc);
  dc(j) = SSD(Tc,Rc,omega,m);
  viewImage(Tc,omega,m); FAIRpause(1/100)
end;
figure(1); clf; surf(w1,w2,dc); hold on; grid off; contour(w1,w2,dc)
title(sprintf('translation, m=[%d,%d]',m)); view(-135,33);
%==============================================================================


text:
function saveNEVSpikes(spikeStruct, newFileName)

% saveNEVSpikes
%
% Allows the user to save a modified version of the spike waveforms into a
% NEV file. This can be very useful for those who want to save MATLAB
% sorted or rethresholded NEV data back into the NEV file. It also re-saves
% the other data in the original NEV file into the new NEV. A
%
% Use saveNEVSpikes(spikeStruct, newFileName)
%
%   spikeStruct: The structure containing the spike data needed to be
%                saved. The structure format must match the one that
%                openNEV outputs.
%                spikeStruct.TimeStamp: contains all the timestamps
%                spikeStruct.Electrode: contains all the electrode #s
%                spikeStruct.Unit: contains all the sorted unit #s
%                spikeStruct.Waveform: Spike waveforms, containing 48 data
%                                      points
%
%   newFileName: The file name of the new NEV file.
%                DEFAULT: User will be prmpted for a file name.
%
%   Example 1:
%   saveNEVSpikes(spikeStruct, 'sortedNEV';
%
%   In the example above, the user will be prompted to select a NEV file.
%   The data stored in spikeStruct will be saved into sortedNEV alongside
%   the data in the user-selected NEV file. The new NEV will be saved as
%   sortedNEV.nev
%
%   Kian Torab
%   ktorab@blackrockmicro.com
%   Blackrock Microsystems
%   Version 1.2.1.0
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Version History
%
% 1.0.0.0: Initial Release.
%
% 1.2.0.0:
%   - Fixed a bug where the data saved incorrectly under the Windows OS.
%   - Sped up the processing significantly.
%
% 1.2.1.0
%   - Fixed a bug introduced in 1.2.0.0.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Validating the input argument
if ~exist('spikeStruct', 'var')
    disp('spikeStruct is a required input argument.');
    return;
end

% if isfield(spikeStruct, 'Data') && isfield(spikeStruct, 'MetaTags')
%     newFileName = fullfile(NEV.MetaTags.FilePath, [NEV.MetaTags.Filename, '.nev']);
%     spikeStruct = NEV.Data.Spikes;
% end

if ~exist('newFileName', 'var')
    newFileName = input('What file name would you like to use to save the new NEV? ', 's');
    if isempty(newFileName)
        disp('A filename is required.');
        return;
    end
end

%% Opening the file and reading header
[dataFilename dataFolder] = getFile('*.nev', 'Choose a NEV that you would like to modify.');
fileFullPath = [dataFolder dataFilename];
FID = fopen(fileFullPath, 'r', 'ieee-le');

%% Calculating the header bytes
BasicHeader = fread(FID, 20, '*uint8');
headerBytes  = double(typecast(BasicHeader(13:16), 'uint32'));
dataPacketByteLength = double(typecast(BasicHeader(17:20), 'uint32'));

%% Calculating the data file length and eeking to the beginning of the file
fseek(FID, 0, 'eof');
endOfDataByte = ftell(FID);
dataByteLength = endOfDataByte - headerBytes;
numberOfPackets = double(dataByteLength) / double(dataPacketByteLength);

%% Reading the header binaries and saving it for future
fseek(FID, 0, 'bof');
headerBinaries = fread(FID, headerBytes, '*uint8');

%% Reading the data binaries
dataBinaries = fread(FID, [dataPacketByteLength numberOfPackets], '*uint8', 0);

%% Finding what PacketIDs have the desired channels
PacketIDs = zeros(1,size(dataBinaries,2));
for IDX = 1:size(dataBinaries,2)
    PacketIDs(IDX) = typecast(dataBinaries(5:6, IDX), 'uint16');
end

%% Only capturing all the binary lines that contain non-spike events, such
% as comments, tracking, patient trigger, etc.
newDataBinaries = dataBinaries(:, PacketIDs > 256); clear dataBinaries;

%% Extracting the timestamps of all the above non-spike events
dataTimestamps = zeros(1,size(newDataBinaries,2));
for IDX = 1:size(newDataBinaries,2)
    dataTimestamps(IDX) = typecast(newDataBinaries(1:4, IDX), 'uint32');
end

%% Converting all the user supplied data into binaries for saving into the
% new structure
newSpikesBinary = zeros(size([typecast(uint32(spikeStruct.TimeStamp(1)), 'uint8'),...
                         typecast(uint16(spikeStruct.Electrode(1)), 'uint8'),...
                         typecast(uint8(spikeStruct.Unit(1)), 'uint8'),...
                         1,...
                         typecast(int16(spikeStruct.Waveform(:,1))', 'uint8')]',1), size(spikeStruct.Electrode,2));
for idx = 1:size(spikeStruct.Electrode,2)
    newSpikesBinary(:,idx) = [typecast(uint32(spikeStruct.TimeStamp(idx)), 'uint8'),...
                         typecast(uint16(spikeStruct.Electrode(idx)), 'uint8'),...
                         typecast(uint8(spikeStruct.Unit(idx)), 'uint8'),...
                         1,...
                         typecast(int16(spikeStruct.Waveform(:,idx))', 'uint8')]';
end

%% Processing data
% Concatinating the user supplied data (non-spike) and the user supplied data
newDataBinaries = [newDataBinaries, newSpikesBinary];

% Ranking the spike and re-ranking the data for saving (timestamp descending)
allTimestamps = [dataTimestamps, spikeStruct.TimeStamp];
[~, ranking] = sort(allTimestamps);
newDataBinaries = newDataBinaries(:,ranking);

%% Saving the new NEV containig the desired channels
FIDw = fopen([dataFolder newFileName '.nev'], 'w+', 'ieee-le');
fwrite(FIDw, headerBinaries, 'uint8');
fwrite(FIDw, newDataBinaries, 'uint8');
fclose(FID);
fclose(FIDw);


text:
function z = trimr(x,n1,n2)
% PURPOSE: return a matrix (or vector) x stripped of the specified rows.
% -----------------------------------------------------
% USAGE: z = trimr(x,n1,n2)
% where: x = input matrix (or vector) (n x k)
%       n1 = first n1 rows to strip
%       n2 = last  n2 rows to strip
% NOTE: modeled after Gauss trimr function
% -----------------------------------------------------
% RETURNS: z = x(n1+1:n-n2,:)
% -----------------------------------------------------

% written by:
% James P. LeSage, Dept of Economics
% University of Toledo
% 2801 W. Bancroft St,
% Toledo, OH 43606
% jpl@jpl.econ.utoledo.edu

  [n junk] = size(x);
  if (n1+n2) >= n;
     error('Attempting to trim too much in trimr');
  end;
  h1 = n1+1;
  h2 = n-n2;
  z = x(h1:h2,:);

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
% OP_UGRADU_JAC_V: assemble the second term of the Jacobian of the
% convective term (see solve_navier_stokes for usage)
%
%   mat = op_ugradu_jac_v (spu, spv, msh, coeff);
%   [rows, cols, values] = op_ugradu_jac (spu, spv, msh, uder);
%
% INPUT:
%
%  spu:   structure representing the space of trial functions (see sp_scalar/sp_evaluate_col)
%  spv:   structure representing the space of test functions  (see sp_scalar/sp_evaluate_col)
%  msh:   structure containing the domain partition and the quadrature rule (see msh_cartesian/msh_evaluate_col)
%  uder:  derivative of the velocity, evaluated in all the degrees of freedom
%
% OUTPUT:
%
%  mat:    assembled jacobian matrix
%  rows:   row indices of the nonzero entries
%  cols:   column indices of the nonzero entries
%  values: values of the nonzero entries
%
% Copyright (C) 2018 Luca Coradello, Luca Pegolotti
%
%    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/>.

function varargout = op_ugradu_jac_v (spu, spv, msh, uder)

  gradu_shpu = sum(bsxfun( @times, reshape(spu.shape_functions, ...
                   [1 size(spu.shape_functions)]),uder),2);

  shpu = reshape (gradu_shpu, spu.ncomp, msh.nqn, spu.nsh_max, msh.nel);
  shpv = reshape (spv.shape_functions, spv.ncomp, msh.nqn, spv.nsh_max, msh.nel);

  rows = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);
  cols = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);
  values = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);

  jacdet_weights = msh.jacdet .* msh.quad_weights;

  ncounter = 0;
  for iel = 1:msh.nel
    if (all (msh.jacdet(:,iel)))
      shpu_iel = reshape (shpu(:, :, 1:spu.nsh(iel), iel), spu.ncomp, msh.nqn, 1, spu.nsh(iel));
      shpv_iel = reshape (shpv(:, :, 1:spv.nsh(iel), iel), spv.ncomp, msh.nqn, spv.nsh(iel), 1);

      jacdet_iel = reshape (jacdet_weights(:,iel), [1,msh.nqn,1,1]);

      jacdet_shpu = bsxfun (@times, jacdet_iel, shpu_iel);
      tmp1 = sum (bsxfun (@times, jacdet_shpu, shpv_iel), 1);
      values(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = reshape (sum (tmp1, 2), spv.nsh(iel), spu.nsh(iel));

      [rows_loc, cols_loc] = ndgrid (spv.connectivity(:,iel), spu.connectivity(:,iel));
      rows(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = rows_loc;
      cols(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = cols_loc;
      ncounter = ncounter + spu.nsh(iel)*spv.nsh(iel);
    else
      warning ('geopdes:jacdet_zero_at_quad_node', 'op_u_v: singular map in element number %d', iel)
    end
  end

  if (nargout == 1 || nargout == 0)
    varargout{1} = sparse (rows(1:ncounter), cols(1:ncounter), ...
                           values(1:ncounter), spv.ndof, spu.ndof);
  elseif (nargout == 3)
    varargout{1} = rows(1:ncounter);
    varargout{2} = cols(1:ncounter);
    varargout{3} = values(1:ncounter);
  else
    error ('op_ugradu_jac_v: wrong number of output arguments')
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
clc;
clearvars;
close all;
rng default;

[X, K, cluster_sizes, true_labels] = prepare_data(true);

fprintf('Computing Low Rank Subspace Clustering\n');

tstart = tic;
% Apply low rank subspace clustering
[A2, C2] = spx.cluster.lrsc.noisy_relaxed(X);
elapsed_time = toc(tstart);
fprintf('SPX Version Time taken: %.2f seconds\n', elapsed_time);

fprintf('Performing clustering: \n');
tstart = tic;
% Adjacency matrix
W = abs(C2);
clustering_result = spx.cluster.spectral.simple.normalized_symmetric_sparse(W, K);
%clustering_result = spx.cluster.spectral.simple.normalized_symmetric_fast(W, K);
cluster_labels = clustering_result.labels;
elapsed_time = toc(tstart);
fprintf('Time taken: %.2f seconds \n', elapsed_time);
comparsion_result = spx.cluster.clustering_error_hungarian_mapping(cluster_labels, true_labels, K);
clustering_error_perc = comparsion_result.error_perc;
clustering_acc_perc = 100 - comparsion_result.error_perc;

fprintf('\nclustering error: %0.2f %%, clustering accuracy: %0.2f %% \n'...
    , clustering_error_perc, clustering_acc_perc);

author: indigits

repo: sparse-plex


text:
% IM = mkSine(SIZE, PERIOD, DIRECTION, AMPLITUDE, PHASE, ORIGIN)
%      or
% IM = mkSine(SIZE, FREQ, AMPLITUDE, PHASE, ORIGIN)
%
% Compute a matrix of dimension SIZE (a [Y X] 2-vector, or a scalar)
% containing samples of a 2D sinusoid, with given PERIOD (in pixels),
% DIRECTION (radians, CW from X-axis, default = 0), AMPLITUDE (default
% = 1), and PHASE (radians, relative to ORIGIN, default = 0).  ORIGIN
% defaults to the center of the image.
%
% In the second form, FREQ is a 2-vector of frequencies (radians/pixel).

% Eero Simoncelli, 6/96.

function [res] = mkSine(sz, per_freq, dir_amp, amp_phase, phase_orig, orig)

%------------------------------------------------------------
%% OPTIONAL ARGS:

if (prod(size(per_freq)) == 2)
  frequency = norm(per_freq);
  direction = atan2(per_freq(1),per_freq(2));
  if (exist('dir_amp') == 1)
    amplitude = dir_amp;
  else
    amplitude = 1;
  end
  if (exist('amp_phase') == 1)
    phase = amp_phase;
  else
    phase = 0;
  end
  if (exist('phase_orig') == 1)
    origin = phase_orig;
  end
  if (exist('orig') == 1)
    error('Too many arguments for (second form) of mkSine');
  end
else
  frequency = 2*pi/per_freq;
  if (exist('dir_amp') == 1)
    direction = dir_amp;
  else
    direction = 0;
  end
  if (exist('amp_phase') == 1)
    amplitude = amp_phase;
  else
    amplitude = 1;
  end
  if (exist('phase_orig') == 1)
    phase = phase_orig;
  else
    phase = 0;
  end
  if (exist('orig') == 1)
    origin = orig;
  end
end

%------------------------------------------------------------

if (exist('origin') == 1)
  res = amplitude*sin(mkRamp(sz, direction, frequency, phase, origin));
else
 res = amplitude*sin(mkRamp(sz, direction, frequency, phase));
end

repo: SelfExSR


text:
clear all; close all;
I=imread('cameraman.tif');
I=im2double(I);
J=fftshift(fft2(I));
[x, y]=meshgrid(-128:127, -128:127);
[M, N]=size(I);
n1=floor(M/2);
n2=floor(N/2);
z=sqrt((x-n1).^2+(y-n2).^2);
D1=10;  D2=30;
n=6;
H1=1./(1+(z/D1).^(2*n));
H2=1./(1+(z/D2).^(2*n));
K1=J.*H1;
K2=J.*H2;
L1=ifft2(ifftshift(K1));
L2=ifft2(ifftshift(K2));
figure;
subplot(131);
imshow(I);
subplot(132);
imshow(real(L1));
subplot(133);
imshow(real(L2))

author: HuangCongQing


text:
function vol = tetrahedronVolume(vertices, varargin)
%TETRAHEDRONVOLUME Signed volume of a tetrahedron.
%
%   VOL = tetrahedronVolume(TETRA)
%   Comptues the siged volume of the tetrahedron TETRA defined by a 4-by-4
%   array representing the polyhedron vertices.
%
%   Example
%     vi = [0 0 0;1 0 0;0 1 0;0 0 1];
%     tetrahedronVolume(vi)
%     ans =
%         0.1667
%
%     [V F] = createTetrahedron;
%     tetrahedronVolume(V)
%     ans =
%         -.3333
%
%   See also
%   meshes3d, createTetrahedron, meshVolume
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2012-04-05,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2012 INRA - Cepia Software Platform.

if nargin == 2
    tetras = varargin{1};
    nTetras = size(tetras, 1);
    vol = zeros(nTetras, 1);
    for i = 1:nTetras
        tetra = tetras(i,:);
        vol(i) = det(bsxfun(@minus, vertices(tetra(2:4),:), vertices(tetra(1),:))) / 6;
    end
    return;
end

% control on inputs
if nargin == 4
    vertices = [vertices ; varargin{1} ; varargin{2} ; varargin{3}];
end

if size(vertices, 1) < 4
    error('Input vertex array requires at least 4 vertices');
end

% compute volume of tetrahedron, using first vertex as origin
vol = det(vertices(2:4,:) - vertices([1 1 1],:)) / 6;

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [Msig1,Msig2,MLong]=afi_blochsim(alpha, T1, T2, TE, TR1, TR2, crushFlag, partialDephasingFlag, partialDephasing, df, Nex, inc)
%IR_BLOCHSIM Bloch simulations of the GRE-AFI pulse sequence.
% Simulates 100 spins params.Nex repetitions of the AFI pulse
% sequences.
%
% params: Struct with the following fields:
%   alpha: Excitation pulse flip angle in radians.
%   TR1: Repetition time 1 (ms).
%   TR2: Repetition time 2 (ms).
%   TE: Echo time (ms).
%   T1: Longitudinal relaxation time (ms).
%   T2: Transverse relaxation time (ms).
%   Nex: Number of excitations
%   df: Off-resonance frequency of spins relative to excitation pulse (in Hz)
%   crushFlag: Numeric flag for perfect spoiling (1) or partial spoiling (2).
%   partialDephasing: Partial dephasing fraction (between [0, 1]). 1 = no dephasing, 0 = complete dephasing (sele
%   inc: Phase spoiling increment in degrees.
%
% Outputs:
%   MLong: Longitudinal magnetization at time TI (prior to excitation pulse).
%   Msig: Complex signal produced by the transverse magnetization at time TE after excitation.
%

%% Set up spin properties
%

Nf = 100;

if partialDephasingFlag
    phi = ((1-Nf/2):Nf/2)/Nf*2*pi*partialDephasing; % Radian phase vector going from 2Pi/Nf to 2Pi in 2Pi/Nf increments.
end

%% Calculate free-precession matrices
%

%"A" is decay and phase gained due to off resonance, "B" is regrowth

% Magnetization decayed (A) and regrowth (B) between the alpha pulse and measurement.
[Ate,Bte] = free_precess(TE,T1,T2,df);

% Magnetization decayed (A) and regrowth (B) between the measurement and the next TR.
[Atr1,Btr1] = free_precess(TR1-TE,T1,T2,df);

% Magnetization decayed (A) and regrowth (B) between the measurement and the next TR.
[Atr2,Btr2] = free_precess(TR2-TE,T1,T2,df);

%% Bloch Simulation
%

M = [zeros(2,Nf);ones(1,Nf)]; % Sets initial magnetization for every spin [0;0;1]
on = ones(1,Nf); % Vector to ensure size of matrices in further calculations

Rfph = 0;       % Rf phase
Rfinc = inc;

for n=1:Nex

    %Signal 1
    MLong = mean(M(3,:)); % Longitudinal magnetization just before excitation pulse

	A = Ate * th_rot(alpha, Rfph);
	B = Bte;

	M = A*M+B*on; % M is rotated, then decayed for TE. Regrowth factor is added.

	Msig1 = sum( squeeze(M(1,:)+1i*M(2,:)) ) / Nf; % Complex signal 1 by adding up all the spins

	M=Atr1*M+Btr1*on; % Relaxation during rest of TR1 after TE

    if crushFlag
        % To make sure spoiling is ideal
        M(1:2, :) = 0;
    elseif partialDephasing
        for k=1:Nf
            M(:,k) = z_rot(phi(k))*M(:,k);  % Dephase spins.
        end
    end

    Rfph = Rfph+Rfinc; % Calculate the next RF phase
    Rfinc = Rfinc+inc; % Calculate the next RF increment

%     Msig2 = Msig1;

    %Signal 2
    MLong = mean(M(3,:)); % Longitudinal magnetization after signal 1

	A = Ate * th_rot(alpha, Rfph);
	B = Bte;

	M = A*M+B*on; % M is rotated, then decayed for TE. Regrowth factor is added.

	Msig2 = sum( squeeze(M(1,:)+1i*M(2,:)) ) / Nf; % Complex signal by adding up all the spins

	M=Atr2*M+Btr2*on; % Relaxation during rest of TR1 after TE

    if crushFlag
        % To make sure spoiling is ideal
        M(1:2, :) = 0;
    elseif partialDephasing
        for k=1:Nf
            M(:,k) = z_rot(phi(k))*M(:,k);  % Dephase spins.
        end
    end

    Rfph = Rfph+Rfinc; % Calculate the next RF phase
    Rfinc = Rfinc+inc; % Calculate the next RF increment

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function B = contiguous(A)
% CONTIGUOUS  find contiguous sequences in an input vector of integers
% B = contiguous(A)
% Take an input vector like i = [ 2 3 4 7 11 12 13 14 17 19 21] and
% return B.start = [2 11], B.end = [4 14]
% meaning that ranges 2:4 and 4:14 are contiguous in A

% using example:
% A      = [ 2  3  4  7 11 12 13 14 17 19 21]
% mask   = [ 0  1  1  0  0  1  1  1  0  0  0  0]
% starts = [ 1  0  0  0  1  0  0  0  0  0  0]
% ends   = [ 0  0  1  0  0  0  0  1  0  0  0]
% i(starts) = [2 11]
% i(ends)   = [4 14]

% AUTHOR: Glenn Thompson, University of Alaska Fairbanks
% Vectorized by Celso Reyes

% forcing i into a column ensures this won't crash for rows OR columns.
mask = [false; diff(A(:))==1; false]; %true where i(N) == i(N-1)+1
starts = mask(1:end-1) < mask(2:end);
endIdx = mask(1:end-1) > mask(2:end);
B.start = A(starts);
B.end = A(endIdx);

% tests:
% v0 = [];
% v1= ([ 2 3 4 7 11 12 13 14 17 19 20])
% v2 = [0, v1, 0];
% contiguous([ 2 3 4 7 11 12 13 14 17 19 20])
% expectedStarts = [2 11 19];
% expectedEnds = [4 14 20];
% contiguous([ 0 2 3 4 7 11 12 13 14 17 19 20 22])
% SAME RESULTS.


text:
%% Hough Lines Transform
% An example using the Hough line detector.
%
% This program demonstrates line finding with the Hough transform.
% We show how to use the OpenCV functions |cv.HoughLines| and |cv.HoughLinesP|
% to detect lines in an image.
%
% Sources:
%
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/houghlines.cpp>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/python/houghlines.py>
% * <https://docs.opencv.org/3.2.0/d9/db0/tutorial_hough_lines.html>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/tutorial_code/ImgTrans/HoughLines_Demo.cpp>
% * <https://docs.opencv.org/3.2.0/d6/d10/tutorial_py_houghlines.html>
%

%% Theory
%
% The explanation below belongs to the book *Learning OpenCV* by
% Bradski and Kaehler.
%
% The Hough Line Transform is a transform used to detect straight lines. To
% apply the Transform, first an edge detection pre-processing is desirable.
%
% As you know, a line in the image space can be expressed with two variables.
% For example:
%
% * In the *Cartesian coordinate system:* Parameters: $(m,b)$.
% * In the *Polar coordinate system:* Parameters: $(r,\theta)$
%
% <<https://docs.opencv.org/3.2.0/Hough_Lines_Tutorial_Theory_0.jpg>>
%
% For Hough Transforms, we will express lines in the _Polar system_. Hence, a
% line equation can be written as:
%
% $$y = \left ( -\frac{\cos \theta}{\sin \theta} \right ) x +
%       \left ( \frac{r}{\sin \theta} \right )$$
%
% Arranging the terms: $r = x \cos \theta + y \sin \theta$
%
% In general for each point $(x_{0}, y_{0})$, we can define the family of
% lines that goes through that point as:
%
% $$r_{\theta} = x_{0} \cdot \cos \theta + y_{0} \cdot \sin \theta$$
%
% Meaning that each pair $(r_{\theta},\theta)$ represents each line that
% passes by $(x_{0}, y_{0})$.
%
% If for a given $(x_{0}, y_{0})$ we plot the family of lines that goes
% through it, we get a sinusoid. For instance, for $x_{0} = 8$ and $y_{0} = 6$
% we get the following plot (in a plane $\theta$ - $r$):
%
% <<https://docs.opencv.org/3.2.0/Hough_Lines_Tutorial_Theory_1.jpg>>
%
% We consider only points such that $r > 0$ and $0< \theta < 2 \pi$.
%
% We can do the same operation above for all the points in an image. If the
% curves of two different points intersect in the plane $\theta$ - $r$, that
% means that both points belong to a same line. For instance, following with
% the example above and drawing the plot for two more points:
% $x_{1} = 4$, $y_{1} = 9$ and $x_{2} = 12$, $y_{2} = 3$, we get:
%
% <<https://docs.opencv.org/3.2.0/Hough_Lines_Tutorial_Theory_2.jpg>>
%
% The three plots intersect in one single point $(0.925, 9.6)$, these
% coordinates are the parameters ($\theta, r$) or the line in which
% $(x_{0}, y_{0})$, $(x_{1}, y_{1})$ and $(x_{2}, y_{2})$ lay.
%
% What does all the stuff above mean? It means that in general, a line can be
% _detected_ by finding the number of intersections between curves. The more
% curves intersecting means that the line represented by that intersection
% have more points. In general, we can define a _threshold_ of the minimum
% number of intersections needed to _detect_ a line.
%
% This is what the Hough Line Transform does. It keeps track of the
% intersection between curves of every point in the image. If the number of
% intersections is above some _threshold_, then it declares it as a line with
% the parameters $(\theta, r_{\theta})$ of the intersection point.
%
%% Standard and Probabilistic Hough Line Transform
%
% OpenCV implements two kind of Hough Line Transforms:
%
% 1) *The Standard Hough Transform*
%
% * It consists in pretty much what we just explained in the previous section.
%   It gives you as result a vector of couples $(\theta, r_{\theta})$
% * In OpenCV it is implemented with the function |cv.HoughLines|
%
% 2) *The Probabilistic Hough Line Transform*
%
% * A more efficient implementation of the Hough Line Transform. It gives as
%   output the extremes of the detected lines $(x_{0}, y_{0}, x_{1}, y_{1})$
% * In OpenCV it is implemented with the function |cv.HoughLinesP|
%

%% Code
%
% This program:
%
% * Loads an image
% * Applies either a _Standard Hough Line Transform_ or a
%   _Probabilistic Line Transform_.
% * Display the original image and the detected line in two windows.
%
% You may observe that the number of lines detected vary while you change the
% _threshold_. The explanation is sort of evident: If you establish a higher
% threshold, fewer lines will be detected (since you will need more points to
% declare a line detected).
%

%%
% Input image
if true
    fname = fullfile(mexopencv.root(), 'test', 'sudoku.jpg');
    thresh = 200;
    threshP = 100;
    minlen = 100;
else
    fname = fullfile(mexopencv.root(), 'test', 'pic1.png');
    thresh = 85;
    threshP = 50;
    minlen = 50;
end
src = cv.imread(fname, 'Color',true);

%%
% Edge Detection
gray = cv.cvtColor(src, 'RGB2GRAY');
edges = cv.Canny(gray, [50, 150], 'ApertureSize',3);
imshow(edges), title('Edges')

%%
% HoughLines: Standard Hough Line Transform
tic
lines = cv.HoughLines(edges, 'Rho',1, 'Theta',pi/180, 'Threshold',thresh);
toc

%%
% draw the lines, and display the result
lines = cat(1, lines{:});
rho = lines(:,1);
theta = lines(:,2);
a = cos(theta); b = sin(theta);
x0 = a.*rho; y0 = b.*rho;
pt1 = round([x0 + 1000*(-b), y0 + 1000*(a)]);
pt2 = round([x0 - 1000*(-b), y0 - 1000*(a)]);
out = cv.line(src, pt1, pt2, ...
    'Color',[0 255 0], 'Thickness',2, 'LineType','AA');
figure, imshow(out), title('Detected Lines')

%%
% HoughLinesP: Probabilistic Hough Line Transform
tic
linesP = cv.HoughLinesP(edges, 'Rho',1, 'Theta',pi/180, ...
    'Threshold',threshP, 'MinLineLength',minlen, 'MaxLineGap',10);
toc

%%
% draw the line segments, and display the result
linesP = cat(1, linesP{:});
outP = cv.line(src, linesP(:,1:2), linesP(:,3:4), ...
    'Color',[0 255 0], 'Thickness',2, 'LineType','AA');
figure, imshow(outP), title('Detected Line Segments')

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [sx,dsdx,dsdP] = sparsify(x,P)
% legacy code
s = warning ('on');
warning ('*** The function `sparsify` is now deprecated. Please see `VBA_sparsifyPrior` for an alternative.')
warning (s);

% fallback
[sx,dsdx,dsdP]  = VBA_sparsifyPrior(x, P);

author: MBB-team

repo: VBA-toolbox


text:
% fig_hyper1
% show hyperbola for various delta values

dlist = [1000 10 4 1];
tmax = 4;
ptype = 'hyper3';
plist = {'quad', 'huber2', 'hyper3', 'lange1', 'lange3', ...
	'cauchy', 'qgg2'};
t = tmax * linspace(-1, 1, 401)';
for ii=1:length(dlist)
%	ptype = plist{ii};
	delta = dlist(ii);
	leg{ii} = [ptype ' \delta = ' num2str(delta)];
	param = [];

	pot = potential_fun(ptype, delta, param);
	pp(:,ii) = pot.potk(t);
	pw(:,ii) = pot.wpot(t);
	pd(:,ii) = pot.dpot(t);
end

if im
	clf
	plot(t, pp), title 'potk'
	legend(leg)
	xlabel 't'
	ylabel 'p(t)'
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [rot, trans, Mf] = regEstRot4(rot, trans, scaleFac, volume, numSlices, sagSize, NCoarseIter, coarseFlag, fineFlag);
%
% regEstRotFunction -   Function to compute automatically the rotation and translation,
%                       to be used as a callback.
%                       It returns the rotation and translation in workspace variables
%                       rot and trans.
%
%  Oscar Nestares - 5/99
%  ON - Added option to use the rotation and translation in bestrotvol
%       as initial alignment (useful to recompute automatically previous
%       manual alignments).
%  AB 4/21    - added NCoarseIter = [];
%  SL 7/15/02 - Created regEstRot4 for use in mrAlign4 (based largely on regEstRot).
%               Changed it from a script into a function and eliminated the need for
%               'Usebestrotvol' parameter and the function regParamInit.
%  SL 7/29/02 - removed NCoarseIter = [];
%  SL 8/02/02 - added coarseFlag and fineFlag parameters

global INPLANE

% registering
[rot, trans, Mf]=regVolInp4(...
                 reshape(volume,[sagSize, numSlices]),...   % volume
                 INPLANE.anat,...                           % inplanes
                 scaleFac,...                               % inverse voxel size for inplanes and volume
                 rot,...                                    % initial rotation
                 trans,...                                  % initial translation
                 NCoarseIter,...                            % number of coarse iterations
                 coarseFlag,...                             % coarse iterations flag
                 fineFlag,...                               % fine iterations flag
                 'regEstFilIntGrad',...                     % function to estimate the intensity grad.
                 0);                                        % Plane by Plane flag = 0 (=>works globaly)

% done

msgbox('Automatic computation done');

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [success] = viewAcceptedRejectedCellExtraction(inputImages,inputSignals,valid,inputMovie,varargin)
	% Plots cell image, peak transients, and image from activity in movie for accepted and rejected cells, useful for publications.
	% Biafra Ahanonu
	% started: 2014.01.03 [19:13:01]
	% inputs
		%
	% outputs
		%

	% changelog
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
	% TODO
		%

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	%========================
	% Int: number of cells to show for EACH accepted and rejected
	options.nCellsShow = 3;
	% Int vector: vector of cell indexes if want to analyze a specific subset
	options.showSpecificCells = [];
	% if loading movie inside function, provide framelist
	options.frameList = [];
	% name in HDF5 file where data is stored
	options.inputDatasetName = '/1';
	% pre-compute signal peaks
	options.signalPeaks = [];
	options.signalPeaksArray = [];
	% Int: 1/2 length in pixels of crop square
	options.cropSizeLength = 20;
	% Float: Value between 0 and 1 of max to crop image to get outline
	options.thresholdOutline = 0.1;
	% ROI for peak signal plotting
	options.peakROI = [-20:20];
	% number of standard deviations above the threshold to count as spike
	options.numStdsForThresh = 3.0;
	% Binary: 1 = normalize the cell images
	options.normalizeCellImages = 1;
	% get options
	options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================

	try
		success = 0;

		% make accepted/rejected vector a logical
		valid = logical(valid);

		% check whether movie is a string or not, load
		inputMovieClass = class(inputMovie);
		if strcmp(inputMovieClass,'char')
			% inputMovie = loadMovieList(inputMovie,'inputDatasetName',options.inputDatasetName,'frameList',options.frameList);
			% [pathstr,name,ext] = fileparts(inputFilePath);
			% options.newFilename = [pathstr '\concat_' name '.h5'];
		end

		% calculate grid size based
		nPlotsPerCell = 3;
		% [xPlot yPlot] = getSubplotDimensions(nPairs);
		xPlot = options.nCellsShow;
		yPlot = nPlotsPerCell*2;

		% get list of cells to show
		if isempty(options.showSpecificCells)
			% get list of cell's SNR
			[signalSnr ~] = computeSignalSnr(inputSignals,'testpeaks',options.signalPeaks,'testpeaksArray',options.signalPeaksArray,'numStdsForThresh',options.numStdsForThresh);
			signalSnr(isnan(signalSnr)) = nanmean(signalSnr(:));
			[signalSnr newIdx] = sort(signalSnr,'descend');

			% get top SNR accepted
			topSnrAccepted = newIdx(valid==1);
			topSnrAccepted = topSnrAccepted(1:options.nCellsShow);

			% get bottom SNR rejected
			bottomSnrRejected = newIdx(valid==0);
			bottomSnrRejected = bottomSnrRejected((end-options.nCellsShow+1):end);
			% bottomSnrRejected = bottomSnrRejected(randperm(length(bottomSnrRejected),options.nCellsShow));

			cellsToShow = cat(2,topSnrAccepted(:),bottomSnrRejected(:))';
			cellsToShow = cellsToShow(:);
		else
			cellsToShow = options.showSpecificCells;
		end

		% calculate signal peaks
		if isempty(options.signalPeaks)
			[signalPeaks, signalPeaksArray] = computeSignalPeaks(inputSignals(cellsToShow,:),'makePlots', 0,'makeSummaryPlots',0,'waitbarOn',1);
		else
			signalPeaks = options.signalPeaks;
			signalPeaksArray = options.signalPeaksArray;
		end

		% get the peak statistics
		[peakOutputStat] = computePeakStatistics(inputSignals(cellsToShow,:),'waitbarOn',1,'testpeaks',signalPeaks,'testpeaksArray',signalPeaksArray,'spikeROI',options.peakROI);

		figure;

		% loop and plot each cell
		tracesHere = inputSignals(cellsToShow,:);
		minValTraces = nanmin(tracesHere(:));
		maxValTraces = nanmax(tracesHere(:));
		for cellNo = 1:length(cellsToShow)
			try
				cellIdx = cellsToShow(cellNo);

				thisTrace = inputSignals(cellIdx,:);

				% get peaks for this cell
				% peakIdxs = signalPeaksArray{cellIdx};
				peakIdxs = signalPeaksArray{cellNo};

				% get cropped view of the algorithm image
				inputImageAddObjCut = getObjCutMovie(inputImages,inputImages(:,:,cellIdx),'createMontage',0,'extendedCrosshairs',0,'crossHairsOn',0,'outlines',0,'waitbarOn',0,'cropSize',options.cropSizeLength,'addPadding',1);
				inputImageAddObjCut = inputImageAddObjCut{1};
				inputImageAddObjCut = inputImageAddObjCut(:,:,cellIdx);

				% get cropped view from the movie
				if ischar(inputMovie)
					objCutMovie = getObjCutMovie(inputMovie,inputImages(:,:,cellIdx),'createMontage',0,'extendedCrosshairs',0,'outlines',0,'waitbarOn',0,'cropSize',options.cropSizeLength,'crossHairsOn',0,'addPadding',1,'frameList',peakIdxs);
				else
					objCutMovie = getObjCutMovie(inputMovie(:,:,peakIdxs),inputImages(:,:,cellIdx),'createMontage',0,'extendedCrosshairs',0,'outlines',0,'waitbarOn',0,'cropSize',options.cropSizeLength,'crossHairsOn',0,'addPadding',1);
				end

				objCutMovie = cat(3,objCutMovie{:});
				movieFrame = nanmean(objCutMovie,3);

				% create RGB with outline of cell
				[thresholdedImages boundaryIndices] = thresholdImages(inputImageAddObjCut,'binary',1,'getBoundaryIndex',1,'threshold',options.thresholdOutline);
				E = normalizeVector(double(movieFrame),'normRange','zeroToOne');
				tmpImg = zeros([size(E)]);
				tmpImg([boundaryIndices{:}]) = 1;
				movieRGB = [];
				movieRGB(:,:,1) = E+tmpImg;
				movieRGB(:,:,2) = E;
				movieRGB(:,:,3) = E;

				subplot(xPlot,yPlot,((cellNo-1)*nPlotsPerCell+1))
					if options.normalizeCellImages==1
						inputImageAddObjCut = normalizeVector(inputImageAddObjCut,'normRange','zeroToOne');
						imagesc(inputImageAddObjCut)
						caxis([0 1]);
					else
						imagesc(inputImageAddObjCut)
					end
					box off;axis off;
					axis equal tight
					if cellNo==length(cellsToShow)
						colorbar('Location','southoutside');
					end

				subplot(xPlot,yPlot,((cellNo-1)*nPlotsPerCell+2))
					imagesc(movieRGB)
					box off;axis off;
					axis equal tight
					% colorbar('Location','southoutside');

					imgRowY = size(movieRGB,1);
					imgColX = size(movieRGB,2);
					MICRON_PER_PIXEL = 2.37;
					options.scaleBarLengthMicron = 20;
					scaleBarLengthPx = options.scaleBarLengthMicron/MICRON_PER_PIXEL;
					% [imgColX-scaleBarLengthPx-round(imgColX*0.05) imgRowY-round(imgRowY*0.05) scaleBarLengthPx 5]
					rectangle('Position',[imgColX-scaleBarLengthPx-imgColX*0.05 imgRowY-imgRowY*0.05 scaleBarLengthPx 5],'FaceColor',[1 1 1],'EdgeColor','none')

				subplot(xPlot,yPlot,((cellNo-1)*nPlotsPerCell+3))

					% spikeCenterTrace = peakOutputStat.spikeCenterTrace{cellIdx};
					% avgSpikeTrace = peakOutputStat.avgSpikeTrace(cellIdx,:);
					spikeCenterTrace = peakOutputStat.spikeCenterTrace{cellNo};
					avgSpikeTrace = peakOutputStat.avgSpikeTrace(cellNo,:);
					% spikeCenterTrace
					% avgSpikeTrace
					peakROI = options.peakROI;
					peakSignalAmplitude = thisTrace(peakIdxs(:));
					[peakSignalAmplitude peakIdx] = sort(spikeCenterTrace(:,round(end/2)+1),'descend');
					spikeCenterTrace = spikeCenterTrace(peakIdx,:);
					if size(spikeCenterTrace,1)>20
						spikeCenterTrace = spikeCenterTrace(1:20,:);
					end

					plot(repmat(peakROI, [size(spikeCenterTrace,1) 1])', spikeCenterTrace','Color',[4 4 4]/8)
					hold on;
					plot(peakROI, avgSpikeTrace,'r', 'LineWidth',3);box off;
					ylim([minValTraces maxValTraces]);
					% plot(peakROI, nanmean(spikeCenterTrace),'Color',[1 0 0 1.0], 'LineWidth',2);box off;
					% add in zero line
					% xval = 0;
					% x=[xval,xval];
					% y=[minValTraces maxValTraces];
					% plot(x,y,'r'); box off;
			catch err
				success = 0;
				display(repmat('@',1,7))
				disp(getReport(err,'extended','hyperlinks','on'));
				display(repmat('@',1,7))
			end
		end

		success = 1;
	catch err
		success = 0;
		display(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		display(repmat('@',1,7))
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [geo,angles]=readBrukerGeometry(folder_or_file,dataset_number)

% Developed by A. Biguri

if nargin==1
    dataset_number=-1;
end

if endsWith(folder_or_file,'.log')
    % is the log file itself
    fid=fopen(folder_or_file);
else
    % is the folder where it lives
    file = dir([folder_or_file,'/*.log']); %
    if isempty(file)
        error(['No .log file found in folder: ', folder_or_file]);
    end

    filename=file(1).name;
    fid=fopen([folder_or_file,'/',filename]);
    xtekctText = textscan(fid, '%s %s', 'Delimiter', '=', 'HeaderLines', 1, 'CommentStyle', '[');
    fclose(fid);

    if ~isempty(xtekctText{2}(strcmpi('Number of connected scans', xtekctText{1})))
        if dataset_number==-1
            error('This folder contains many datasets, please select which one to load with BrukerDataLoader(...,`dataset_number`, dataset_number)');
        end
        matching=[];
        for ii=1:length(file)
           if strcmpi(file(ii).name(end-5:end-4),num2str(dataset_number,'%02d'))
               matching=[matching ii];
           end
        end
        if length(matching)>1
            error('More than 1 file for the same dataset found, confused what to do, so I error')
        elseif isempty(matching)
            error('Dataset not found, check the number')
        end
        id=matching;
    else
        id=1;
    end
    filename=file(id).name;
    fid=fopen([folder_or_file,'/',filename]);
end
% check if it was oppened right
if(fid==-1)
    error('Wrong file path');
end
xtekctText = textscan(fid, '%s %s', 'Delimiter', '=', 'HeaderLines', 1, 'CommentStyle', '[');
fclose(fid);

%% Detector information
% Number of pixel in the detector
geo.nDetector=[str2double(xtekctText{2}(strcmpi('Number of Columns', xtekctText{1})));
    str2double(xtekctText{2}(strcmpi('Number of Rows', xtekctText{1})))];
% Size of pixels in the detector
geo.dDetector=[str2double(xtekctText{2}(strcmpi('Camera Pixel Size (um)', xtekctText{1})))/1000;
    str2double(xtekctText{2}(strcmpi('Camera Pixel Size (um)', xtekctText{1})))/1000];

camera_binning=cell2mat(xtekctText{2}(strcmpi('Camera binning', xtekctText{1})));
camera_binning=[str2double(camera_binning(1)), str2double(camera_binning(3))];
if ~isempty(camera_binning)
    geo.dDetector=geo.dDetector.*camera_binning';
end

xyratio= str2double(xtekctText{2}(strcmpi('CameraXYRatio', xtekctText{1})));
if ~isempty(xyratio)
    % try again
    xyratio= str2double(xtekctText{2}(strcmpi('Camera X/Y Ratio', xtekctText{1})));
    if ~isempty(xyratio)
        geo.dDetector(2)=geo.dDetector(2)/xyratio;
    end
end
% Total size of the detector
geo.sDetector=geo.nDetector.*geo.dDetector;

%% Offset of the detector:
geo.offDetector=[0;-(geo.nDetector(2)/2-str2double(xtekctText{2}(strcmpi('Optical Axis (line)', xtekctText{1}))))].*geo.dDetector;

%% Image information
% Number of pixel in the detector

% Size of each pixel
geo.dVoxel=[str2double(xtekctText{2}(strcmpi('Image Pixel Size (um)', xtekctText{1})))/1000;
    str2double(xtekctText{2}(strcmpi('Image Pixel Size (um)', xtekctText{1})))/1000;
    str2double(xtekctText{2}(strcmpi('Image Pixel Size (um)', xtekctText{1})))/1000 ];
% Size of the image in mm
geo.nVoxel=[geo.nDetector(1);geo.nDetector(1);geo.nDetector(2)];
geo.sVoxel=geo.nVoxel.*geo.dVoxel;

geo.offOrigin=[0;0;0];
%% Global geometry
geo.DSO=str2double(xtekctText{2}(strcmpi('Object to Source (mm)', xtekctText{1})));
geo.DSD=str2double(xtekctText{2}(strcmpi('Camera to Source (mm)', xtekctText{1})));

%% Detector offset
if endsWith(folder_or_file,'.log')
    folder=folder_or_file(1:end-4);
else
    folder=folder_or_file;
end
file = dir([folder,'/*.csv']); %
if ~isempty(file)
    off=csvread([file(1).folder, '/', file(1).name],5,1).'.*(geo.dDetector./camera_binning');
    geo.offDetector=geo.offDetector+off(:,1:end);
end
%% whitelevel

geo.whitelevel=2^str2double(xtekctText{2}(strcmpi('Depth (bits)', xtekctText{1})));

%%
%% angles

angle_step=str2double(xtekctText{2}(strcmpi('Rotation Step (deg)', xtekctText{1})));
initial_angle=0;
n_angles=str2double(xtekctText{2}(strcmpi('Number of Files', xtekctText{1})))-1;
angles=initial_angle:angle_step*pi/180:(initial_angle+(n_angles-1)*angle_step)*pi/180;
assert(size(angles,2)==n_angles,'Assertion failed: Inconsistent data detected. Number of projections and angle information do not match\n');

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function mColormap = gui_Colormap_Rastafari(nSize)

% If size is not specified, set it to 256
if nargin < 1
  nSize = 256;
end

mColormap = [];

nSteps = floor(nSize * 0.1);
mColormap = [mColormap; gui_Interpolate([0.036578, 0.159091, 0.015374], [0.022239, 0.234849, 0.007687], nSteps)];

nSteps = floor(nSize * 0.05);
mColormap = [mColormap; gui_Interpolate([0.022239, 0.234849, 0.007687], [0.007899, 0.310606, 0.000000], nSteps)];

nSteps = floor(nSize * 0.09);
mColormap = [mColormap; gui_Interpolate([0.007899, 0.310606, 0.000000], [0.101896, 0.443182, 0.042827], nSteps)];

nSteps = floor(nSize * 0.09);
mColormap = [mColormap; gui_Interpolate([0.101896, 0.443182, 0.042827], [0.195893, 0.575758, 0.085655], nSteps)];

nSteps = floor(nSize * 0.09);
mColormap = [mColormap; gui_Interpolate([0.195893, 0.575758, 0.085655], [0.560068, 0.663178, 0.139025], nSteps)];

nSteps = floor(nSize * 0.08);
mColormap = [mColormap; gui_Interpolate([0.560068, 0.663178, 0.139025], [0.924242, 0.750598, 0.192395], nSteps)];

nSteps = floor(nSize * 0.1);
mColormap = [mColormap; gui_Interpolate([0.924242, 0.750598, 0.192395], [0.939393, 0.495226, 0.162308], nSteps)];

nSteps = floor(nSize * 0.1);
mColormap = [mColormap; gui_Interpolate([0.939393, 0.495226, 0.162308], [0.954545, 0.239854, 0.132221], nSteps)];

nSteps = floor(nSize * 0.2);
mColormap = [mColormap; gui_Interpolate([0.954545, 0.239854, 0.132221], [0.742424, 0.279602, 0.184117], nSteps)];

nSteps = floor(nSize * 0.1);
mColormap = [mColormap; gui_Interpolate([0.742424, 0.279602, 0.184117], [0.530303, 0.319349, 0.236012], nSteps)];

mColormap(mColormap < 0) = 0;

author: CelsoReyes


text:
function gf = gradddivHexpS51(Hkt_r, varargin)
% function gf = gradddivH(Hkt_r, varargin)
% gf is 1 x k*t row vector
% Hkt_r = reshape(Hkt,1,k*t) -> row vector
% Vxt = varargin{1};  %data
% Wxk = varargin{2};  %W matrix
% Wxk_fix = varargin{3}; %fixed part of the Wxk matrix (e.g. background) ->rows
% Hkt_fix = varargin{4}; %fixed part (lines) of the H matrix (e.g. background)

alphaH=1; %for now....
Vxt = varargin{1};  %data
Wxk_tmp = varargin{2};  %W matrix
Wxk_fix = varargin{3}; %fixed part of the Wxk matrix (e.g. background) ->rows
Hkt_fix = varargin{4}; %fixed part (lines) of the H matrix (e.g. background)

t=size(Vxt,2);
k=size(Wxk_tmp,2);

Wxk = [Wxk_tmp, Wxk_fix];
Hkt=[exp(reshape(Hkt_r,k,t)); Hkt_fix];

deltasum=sum(sum(Wxk_tmp,1))-k;
if abs(deltasum)>1e-6
    error('Wxk is not correctly normalized! (sum(Wxk_tmp,1)<>1)\n sum(Wxk_tmp,1)=%f',deltasum)
end

gfkt = (1-Wxk'*(Vxt./(Wxk*Hkt)))*alphaH.*Hkt; %d/dh(d-divergence)
% one is tehre because Wxt is normalized: sum(Wxt,1)=1
gf=reshape(gfkt(1:k,:),1,k*t); %making row vector

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function varargout = stlread(file)
% STLREAD imports geometry from an STL file into MATLAB.
%    FV = STLREAD(FILENAME) imports triangular faces from the binary STL file
%    indicated by FILENAME, and returns the patch struct FV, with fields 'faces'
%    and 'vertices'.
%
%    [F,V] = STLREAD(FILENAME) returns the faces F and vertices V separately.
%
%    [F,V,N] = STLREAD(FILENAME) also returns the face normal vectors.
%
%    The faces and vertices are arranged in the format used by the PATCH plot
%    object.

%    Eric C. Johnson, 11-Dec-2008
%    Copyright 1999-2008 The MathWorks, Inc.

    if ~exist(file,'file')
        error(['File ''%s'' not found. If the file is not on MATLAB''s path' ...
               ', be sure to specify the full path to the file.'], file);
    end

    fid = fopen(file,'r');
    if ~isempty(ferror(fid))
        error(lasterror); %#ok
    end

    M = fread(fid,inf,'uint8=>uint8');
    fclose(fid);

    if( isbinary(M) )
       [f,v,n] = stlbinary(M);
    else
       [f,v,n] = stlascii(M); % Not currently supported
    end

    varargout = cell(1,nargout);
    switch nargout
        case 2
            varargout{1} = f;
            varargout{2} = v;
        case 3
            varargout{1} = f;
            varargout{2} = v;
            varargout{3} = n;
        otherwise
            varargout{1} = struct('faces',f,'vertices',v);
    end

end

function [F,V,N] = stlbinary(M)

    F = [];
    V = [];
    N = [];

    if length(M) < 84
        error('MATLAB:stlread:incorrectFormat', ...
              'Incomplete header information in binary STL file.');
    end

    % Bytes 81-84 are an unsigned 32-bit integer specifying the number of faces
    % that follow.
    numFaces = typecast(M(81:84),'uint32');
    %numFaces = double(numFaces);
    if numFaces == 0
        warning('MATLAB:stlread:nodata','No data in STL file.');
        return
    end

    T = M(85:end);
    F = NaN(numFaces,3);
    V = NaN(3*numFaces,3);
    N = NaN(numFaces,3);

    numRead = 0;
    while numRead < numFaces
        % Each facet is 50 bytes
        %  - Three single precision values specifying the face normal vector
        %  - Three single precision values specifying the first vertex (XYZ)
        %  - Three single precision values specifying the second vertex (XYZ)
        %  - Three single precision values specifying the third vertex (XYZ)
        %  - Two unused bytes
        i1    = 50 * numRead + 1;
        i2    = i1 + 50 - 1;
        facet = T(i1:i2)';

        n  = typecast(facet(1:12),'single');
        v1 = typecast(facet(13:24),'single');
        v2 = typecast(facet(25:36),'single');
        v3 = typecast(facet(37:48),'single');

        n = double(n);
        v = double([v1; v2; v3]);

        % Figure out where to fit these new vertices, and the face, in the
        % larger F and V collections.
        fInd  = numRead + 1;
        vInd1 = 3 * (fInd - 1) + 1;
        vInd2 = vInd1 + 3 - 1;

        V(vInd1:vInd2,:) = v;
        F(fInd,:)        = vInd1:vInd2;
        N(fInd,:)        = n;

        numRead = numRead + 1;
    end

end

function [F,V,N] = stlascii(A) %#ok
    warning('MATLAB:stlread:ascii','ASCII STL files currently not supported.');
    F = [];
    V = [];
    N = [];
end

function tf = isbinary(A)
% ISBINARY determines if an STL file is binary or ASCII.

    % Look for the string 'endsolid' near the end of the file
    if isempty(A) || length(A) < 16
        error('MATLAB:stlread:incorrectFormat', ...
              'File does not appear to be an ASCII or binary STL file.');
    end

    % Read final 16 characters of M
    i2  = length(A);
    i1  = i2 - 15;
    str = char( A(i1:i2)' );

    k = strfind(lower(str), 'endsolid');
    if ~isempty(k)
        tf = false; % ASCII
    else
        tf = true;  % Binary
    end
end


text:
firstPart  = fullfile( '/home','alex','Desktop');
secondPart = fullfile('TheFolder/');
outPutPath = fullfile(firstPart,secondPart);
d.fileName = 'VademecumSmoothCorner';
d.outPutPath = fullfile(outPutPath,[d.fileName,'/']);
vc = VademecumComputerAndPlotter(d);
vc.compute();

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [ r, seed ] = r8vec_uniform_01 ( n, seed )

%*****************************************************************************80
%
%% R8VEC_UNIFORM_01 returns a unit pseudorandom R8VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Second Edition,
%    Springer, 1987,
%    ISBN: 0387964673,
%    LC: QA76.9.C65.B73.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, December 1986, pages 362-376.
%
%    Pierre L'Ecuyer,
%    Random Number Generation,
%    in Handbook of Simulation,
%    edited by Jerry Banks,
%    Wiley, 1998,
%    ISBN: 0471134031,
%    LC: T57.62.H37.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, Number 2, 1969, pages 136-143.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real R(N), the vector of pseudorandom values.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  i4_huge = 2147483647;

  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8VEC_UNIFORM_01 - Fatal error!\n' );
    fprintf ( 1, '  Input SEED = 0!\n' );
    error ( 'R8VEC_UNIFORM_01 - Fatal error!' );
  end

  for i = 1 : n

    k = floor ( seed / 127773 );

    seed = 16807 * ( seed - k * 127773 ) - k * 2836;

    if ( seed < 0 )
      seed = seed + i4_huge;
    end

    r(i) = seed * 4.656612875E-10;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ nodes, weights ] = tensor_product ( n1D, w1D )

%*****************************************************************************80
%
%% TENSOR_PRODUCT generates a tensor product quadrature rule.
%
%  Discussion:
%
%    The Kronecker product of an K by L matrix A and an M by N matrix B
%    is the K*M by L*N matrix formed by
%
%      a(1,1) * B,  a(1,2) * B,  ..., a(1,l) * B
%      a(2,1) * B,  a(2,2) * B,  ..., a(2,l) * B
%      ..........   ..........   .... ..........
%      a(k,1) * B,  a(k,2) * B,  ..., a(k,l) * B
%
%    Thanks to Ivan Ukhov for pointing out a tiny but deadly typographical
%    error, 17 July 2012.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 July 2012
%
%  Author:
%
%    Original MATLAB version by Florian Heiss, Viktor Winschel.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Florian Heiss, Viktor Winschel,
%    Likelihood approximation by numerical integration on sparse grids,
%    Journal of Econometrics,
%    Volume 144, 2008, pages 62-80.
%
%  Parameters:
%
%    Input, cell array n1D{}, contains K sets of 1D nodes.
%    The I-th set has dimension N(I).
%    Each entry of the cell array should be a column vector.
%
%    Input, cell array w1D{}, contains K sets of 1D weights.
%
%    Output, real nodes(NPROD,K), the tensor product nodes.
%    NPROD = product ( N(1) * ... * N(K) ).
%
%    Output, real weights(NPROD), the tensor product weights.
%
  dimension = length ( n1D );

  nodes = n1D{1};
  nodes = nodes ( : );
  weights = w1D{1};
  weights = weights ( : );

  for j = 2 : dimension

    newnodes = n1D{j};
    newnodes = newnodes ( : );

    a = ones ( size ( newnodes, 1 ), 1 );
    b = ones ( size ( nodes, 1 ), 1 );
    c = kron ( nodes, a );
    d = kron ( b, newnodes );

    nodes = [ c, d ];

    newweights = w1D{j};
    newweights = newweights ( : );
    weights = kron ( weights, newweights );

  end

  return
end


text:
% Print CERR screen
%
% Shunde 7/14/05

% 24 Apr 2006 KU  Changed to 'print' (with no UI controls) instead of printpreview because of bug(?) in
%                 Matlab 7.2.  Printpreview hangs up if figure is in fullscreen mode.  Maybe try
%                 printpreview again in future version??
% 27 Dec 2006 KU  Switched back to zbuffer for CERR 3.0.
%
% 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/>.

% set(gcf,'renderer','painters');
% set(gcf,'PaperPositionMode','auto');
% orient landscape;
% printpreview;

set(gcf,'renderer','zbuffer');
set(gcf,'PaperPositionMode','auto');
orient landscape;

sentence1 = {'To adjust the size of the printed figure, adjust the screensize of the ',...
                'CERR window rather than changing the printer options.  (The CERR window ',...
                'will fill a standard 8.5x11 paper if it is adjusted to be fullscreen ',...
                'with a screen setting of 1024x768.)'};
            h=msgbox(sentence1, 'Printing Tip', 'help');
            waitfor(h);
if ispc
    print -v -noui;
else
    print -dsetup -noui;    %untested code for non-windows platforms
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ pr, pq ] = plane_normal_basis_3d ( pp, normal )

%*****************************************************************************80
%
%% PLANE_NORMAL_BASIS_3D finds two perpendicular vectors in a plane in 3D.
%
%  Discussion:
%
%    The normal form of a plane in 3D is:
%
%      PP is a point on the plane,
%      N is a normal vector to the plane.
%
%    The two vectors to be computed, PQ and PR, can be regarded as
%    the basis of a Cartesian coordinate system for points in the plane.
%    Any point in the plane can be described in terms of the "origin"
%    point PP plus a weighted sum of the two vectors PQ and PR:
%
%      P = PP + a * PQ + b * PR.
%
%    The vectors PQ and PR have unit length, and are perpendicular to N
%    and to each other.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real PP(3,1), a point on the plane.  (Actually,
%    we never need to know these values to do the calculation!)
%
%    Input, real NORMAL(3,1), a normal vector N to the plane.  The
%    vector must not have zero length, but it is not necessary for N
%    to have unit length.
%
%    Output, real PQ(3,1), a vector of unit length,
%    perpendicular to the vector N and the vector PR.
%
%    Output, real PR(3,1), a vector of unit length,
%    perpendicular to the vector N and the vector PQ.
%
  dim_num = 3;
%
%  Compute the length of NORMAL.
%
  normal_norm = norm ( normal );

  if ( normal_norm == 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PLANE_NORMAL_BASIS_3D - Fatal error!\n' );
    fprintf ( 1, '  The normal vector is 0.\n' );
    error ( 'PLANE_NORMAL_BASIS_3D - Fatal error!' );
  end
%
%  Find a vector PQ that is normal to NORMAL and has unit length.
%
  pq = r8vec_any_normal ( 3, normal );
%
%  Now just take the cross product NORMAL x PQ to get the PR vector.
%
  pr = r8vec_cross_product_3d ( normal, pq );

  pr_norm = norm ( pr );

  pr(1:dim_num,1) = pr(1:dim_num,1) / pr_norm;

  return
end


text:
%                                                    cpx = whichcpx(K)
% WHICHCPX  yields structure cpx.{f,q,r,x}
%
% ******************** INTERNAL FUNCTION OF SEDUMI ********************
%
% See also sedumi

function cpx = whichcpx(K) %#ok
%
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% 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.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA
%

disp('The SeDuMi binaries are not installed.')
disp('In Matlab, launch "install_sedumi" in the folder you put the SeDuMi files.')
disp('For more information see the file Install.txt.')
error(' ')

repo: REST


text:
function varargout = histogram(varargin)
% GUI for Plotting historgam and equalizing image arrays
%
% INPUT:
% images - array to plot histograms for (can be a single image as well)
% 'N' - image to display
%
% to run:
%
%   histogram(images,'1')
%
% if want equalized images as output:
%
%   im_eq=threshold(images,'1');
%
%  created by K.Artyushkova
%  January 2004

% Kateryna Artyushkova
% Postdoctoral Scientist
% Department of Chemical and Nuclear Engineering
% The University of New Mexico
% (505) 277-0750
% kartyush@unm.edu

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @histogram_OpeningFcn, ...
                   'gui_OutputFcn',  @histogram_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

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

% --- Executes just before histogram is made visible.
function histogram_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to histogram (see VARARGIN)

% Update handles structure
axes(handles.axes1);
cla
data= varargin{1};
n=varargin{2};
n=str2double(n);
iptsetpref('ImshowAxesVisible', 'on')
imagesc(data(:,:,n)), colormap(gray)
set(handles.imi,'string',n);
handles.data=data;
handles.n=n;
guidata(hObject, handles);

% UIWAIT makes histogram wait for user response (see UIRESUME)
% uiwait(handles.figure1);

uiwait(handles.figure1);

% --- Executes during object creation, after setting all properties.
function slice_select_CreateFcn(hObject, eventdata, handles)
% hObject    handle to slice_select (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background, change
%       'usewhitebg' to 0 to use default.  See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor',[.9 .9 .9]);
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on slider movement.
function slice_select_Callback(hObject, eventdata, handles)
% hObject    handle to slice_select (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
data=handles.data;
[m,p,q]=size(data);
set(handles.im1,'string',1);
set(handles.imN,'string',q);
step=1/q;
slider_step(1)=step;
slider_step(2)=step;
set(handles.slice_select, 'SliderStep', slider_step, 'Max', q, 'Min',0)
i=get(hObject,'Value');
i=round(i);
    if i==0
        i=1;
    elseif i>=q
        i=q;
    else i=i;
    end
    set(handles.imi,'string',i);
    axes(handles.axes1);
iptsetpref('ImshowAxesVisible', 'on')
imagesc(data(:,:,i)), colormap(gray)
    axes(handles.axes2);
    hist(data(:,:,i))
    handles.n=i;
    guidata(hObject,handles)

% --- Executes on button press in equalize.
function equalize_Callback(hObject, eventdata, handles)
% hObject    handle to equalize (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

data=handles.data;
[m,p,q]=size(data);
for i=1:q
    data_eq(:,:,i)=histeq(uint8(data(:,:,i)));
end
data_eq=double(data_eq);

n=handles.n;
axes(handles.axes1);
iptsetpref('ImshowAxesVisible', 'on')
imagesc(data_eq(:,:,i)), colormap(gray)
      axes(handles.axes2);
    hist(data_eq(:,:,n))
handles.data_eq=data_eq;
    guidata(hObject,handles)

% --- Executes on button press in replace.
function varargout=replace_Callback(hObject, eventdata, handles)
% hObject    handle to replace (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

data=handles.data_eq;
handles.data=data;
handles.output=data;
varargout{1} = handles.output;
guidata(hObject,handles)
uiresume(handles.figure1);

% --- Executes on button press in cancel.
function cancel_Callback(hObject, eventdata, handles)
% hObject    handle to cancel (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
uiresume(handles.figure1);

    % --- Outputs from this function are returned to the command line.
function varargout = histogram_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% Get default command line output from handles structure

data=handles.data;
handles.output=data;
varargout{1} = handles.output;
close

author: Sable

repo: mcbench-benchmarks


text:
% ir_mri_al_p2_example.m
% Example of MRI SENSE reconstruction from under-sampled k-space data
% using the AL-P2 algorithm from Mar. 2011 T-MI paper.
% Sathish Ramani, circa 2011
% 2013-07-14, modified by Jeff Fessler

%% Test image - Shepp-Logan Phantom
if ~isvar('img')
	N = 256;
	if exist('phantom') % matlab with image processing toolbox
		img = phantom('Modified Shepp-Logan', N)';
		img = 256 * (img + 0.05);
	else
		img = ellipse_im(256);
		img = 256 * (img/max(img(:)) + 0.05);
	end
	[rs, cs] = size(img); Npix = rs*cs; mn = min(img(:)); mx = max(img(:));
	params.rs = rs;
	params.cs = cs;
	params.mn = mn;
	params.mx = mx;
	params.Npix = Npix;
	im plc 2 3
	im(1, img, 'true image')
end

%% sensitivity maps
if ~isvar('smap')
	smap = mri_sensemap_sim('nx', N, 'rcoil', 175, 'dx', 1);
	ncoils = size(smap, 3);
	params.ncoils = ncoils;
	im(2, abs(smap), '|sense maps|'), cbar
end

%% load sampling pattern
if ~isvar('SP')
	load('PDiskR256x256L1R2.252.25B0.8R1D1pctg20.2972.mat', 'SP', 'sampname');
	Ncentx = 8; % sample a window of 2*N around DC along x
	Ncenty = 8; % sample a window of 2*N around DC along y

	params.samp = sampname;
	SP = coverDC_SamplingMask(SP, Ncentx, Ncenty);
	lSP = length(find(SP>0)); % # of sampling points
	printm('Sampling %0.2f%% of k-space', lSP / Npix * 100)

	% Sampling pattern for the output of each coil
	SP3 = repmat(SP, [1 1 ncoils]);
	im(3, [-N/2:N/2-1], [-N/2:N/2-1], fftshift(SP), 'sampling pattern')
end

%% Data samples
if ~isvar('Data')
	SNR = 30; % SNR of noisy data
	smap_img = smap .* repmat(img, [1 1 ncoils]); % sensitivity-weighted images
	smap_Img = fft2(smap_img); % FFT of Sensitivity-weighted images

	% Add white Gaussian noise of desired SNR
	seedr = 0;
	stdnois = zeros(1, ncoils);
	noise_vr = zeros(rs, cs, ncoils);
	for ismap = 1:ncoils
		Dtemp = smap_Img(:,:,ismap);
		sigpow = sum(abs(Dtemp(:)).^2)/(rs*cs);
		stdnois(ismap) = sqrt(sigpow*10^(-SNR/10)); % Std. deviation of noise at ith coil
		randn('state', seedr + ismap);
		noise_realz = (randn(rs,cs) + j*randn(rs,cs))/sqrt(2);
		noise_vr(:, :, ismap) = noise_realz;
	end
	SQCStd = diag(stdnois);
	noisecorr = zeros(rs, cs, ncoils);
	for irK = 1 : ncoils % Apply the above correlation to the generated i.i.d noise to create correlated noise
		dump(1, 1, :) = SQCStd(irK, :);
		noisecorr(:, :, irK) = sum(repmat(dump, [rs cs 1]) .* noise_vr, 3);
	end
	fullData = smap_Img + noisecorr; % Noisy fully-sampled Data
	Data = fullData .* SP3; % Generate undersampled k-space data by applying a mask
	fftshift2 = @(x) fftshift(fftshift(Data, 1), 2);
	im(5, max(log(abs(fftshift2(Data))), 0), 'Data')
prompt
end

%% Simple iFFT / sum-of-squares recon of raw k-space data
if ~isvar('recon_sos')
	recon_sos = ifft2(Data);
	recon_sos = sqrt(sum(abs(recon_sos).^2, 3));
	params.recon_sos = recon_sos;
	im(4, recon_sos, 'SoS recon')
end

%% Initial Estimate = SoS + slight perturbation
if ~isvar('xini')
	randn('state', seedr);
	xini = recon_sos + (randn(rs,cs) + sqrt(-1) * randn(rs,cs))*0.001;
	errini = xini - img;
	errini = sqrt(sum(abs(errini(:)).^2)/Npix)/mx;
	xlabelf('nrmse %g dB', 20*log10(nrms(xini, img)))
	im(6, abs(xini-img), '|SoS error|')
prompt
end

%% Noise decorrelation in data and redefined
% (square root of inverse of noise correlation matrix weighted) sensitivity maps
if ~isvar('smap2')
	recon_F = Npix * ifft2(SP3 .* Data);
	recon_SF = sum(conj(smap) .* recon_F, 3);

	params.smap = smap;
	smap2 = sum(abs(smap).^2, 3); % S'S
%	im(smap2)
end

%% Penalty / regularizer / prior
if ~isvar('lod')
	params.Prior.PriorType{1} = 'TV'; % Type of penalty for wavelet coef.
	params.Prior.PriorType{2} = 'TV'; % Type of penalty for wavelet coef.

	% Type of sparsifying transform
	params.Operator = 'FD'; % finite differences
	% params.Operator = 'W'; % Redundant wavelet transform
	% params.Operator = 'WFD'; % Redundant wavelet transform & finite diff.

	% Wavelet Options
	params.Wavelet.redundancy = 'undecimated'; % Undecimated wavelet transform using wavelet filters corresponding to standard orthonormal wavelets
	params.Wavelet.wname = 'haar'; % undecimated Haar wavelets
	params.Wavelet.nlev = 2; % wavelet transform with 2 levels
	params.Wavelet.includeApprox = false; % exclude approximation level in the regularizer

%	todo: make it ok in octave
%	Uw = Gwave2('mask', true(rs,cs), 'nlevel', params.Wavelet.nlev)'; % trick

	dwtmode('per', 'nodisp'); % Period boundaries for wavelet implementation
	% params.Wavelet.redundancy = 'none'; % decimatSENSERecon_MRI029_Spiral6_23_2.1_SNR20_OpFD_PriorTV_Lam13.5_variablesed Orthonomal wavelet transform

	% Wavelet filters
	[lod, hid, lor, hir] = wfilters(params.Wavelet.wname);

	% Normalize filters so as to avoid a product with 0.5 during inverse undecimated wavelet transform
	params.Wavelet.lod = lod/sqrt(2);
	params.Wavelet.hid = hid/sqrt(2);
	params.Wavelet.lor = lor/sqrt(2);
	params.Wavelet.hir = hir/sqrt(2);
end

%% Parameters of the algorithms
if ~isvar('RR')
	% Regularization parameters for the wavelet and TV prior
	params.lambda = mean(stdnois .^ 2) * [1 0.1];
	params.precon = 0;
	params.betaD = 1e-10;
	params.xinf = zeros(rs,cs);
	params.xinfnorm = 1;
	params.dcosttol = 1e18;
	params.dxtol = 0;

	% Parameters for size of vectors
	[lzALP2 sALP2 eALP2 sr er] = get_Size_AuxVar(params);

	params.AL.lzALP2 = lzALP2;
	params.AL.sALP2 = sALP2;
	params.AL.eALP2 = eALP2;

	params.AL.sr = sr;
	params.AL.er = er;

	RR = compute_RR(params); % FFT of R^T * R
	mxRR = max(RR(:));
	mnRR = min(RR(:));
	params.AL.RR = RR;

	% Other setting
	params.figno = 1;
	params.subplot_img = 5;
	params.subplot_err = 6;
	params.maxitr = 400;
	params.maxitr_in = 1;
	params.dispfig = 1;
	params.dispitr = 1;
	params.dispitrnum = 10;
end

%% prepare to run AL-P2
if ~isvar('zALP2')
	mxsmap2 = max(smap2(:));
	mnsmap2 = min(smap2(:));
	condSS =  mxsmap2 / mnsmap2; % Condition number of S'S

	kapx = 0.9 * condSS;
	kapu0 = 24; % condition number of (FF + mu*I)
	kapu2 = 12; % condition number of (RR + nu2/n1*I)

	nu2 = (mxsmap2 - mnsmap2 * kapx) / (kapx - 1);
	mu = Npix / (kapu0 - 1); % min eig val of F'F = 0, while max eig val of F'F = N where the DFT matrix F is assumed to be un-normalized
	nu1 = (kapu2 - 1) * nu2 / (mxRR - mnRR * kapu2);
	nuratio = nu2 / nu1;

	% warn if any of the mu's is negative or less than eps
	if ((mu <= 0) || (nu1 <= 0) || (nu2 <= 0))
		error('At least one of the mu, nu''s is not positive');
	end

	% Inverse of some matrices required for solving sub-problems
	iPpmu = 1 ./ (Npix * SP3 + mu); % Freq. Response of (F'F + mu * I)^-1
	iSpnu2 = 1./(smap2 + nu2); % Inverse of (S'S + nu2 * I)
	iRpnu2nu1 = 1./(RR + nuratio); % Freq. Response of (R'R + nu2/nu1 * I)^-1

	params.AL.mu = mu;
	params.AL.nu1 = nu1;
	params.AL.nu2 = nu2;
	params.AL.iPpmu = iPpmu;
	params.AL.iSpnu2= iSpnu2;
	params.AL.iRpnu2nu1 = iRpnu2nu1;

	dALP2 = zeros(rs, cs, params.AL.lzALP2);
	zALP2 = dALP2;
end

%% Run AL-P2!
if ~isvar('x_est'), printm 'running al-p2'
	[x_est CALP2 TALP2 EALP2 ERRALP2] = SENSERecon_ALP2(img, ...
		SP3, Data, recon_F, dALP2, zALP2, xini, params);
	printm('Recon time = %g secs', sum(TALP2(:)))
	printm('Recon NRMSE = %g ', ERRALP2(end))
end